annotate opt.cpp @ 13:de3961f74f30 tip

Add Linux/gcc Makefile; build fix

/*
    Harmonic sinusoidal modelling and tools

    C++ code package for harmonic sinusoidal modelling and relevant signal processing.
    Centre for Digital Music, Queen Mary, University of London.
    This file copyright 2011 Wen Xue.

    This program is free software; you can redistribute it and/or
    modify it under the terms of the GNU General Public License as
    published by the Free Software Foundation; either version 2 of the
    License, or (at your option) any later version.
*/
//---------------------------------------------------------------------------

#include <math.h>
#include <memory.h>
#include <stdlib.h>
#include "opt.h"
#include "matrix.h"

/** @file opt.h - optimization routines */

//---------------------------------------------------------------------------
//macro nsign: judges if two double-precision floating pointer numbers have different signs
#define nsign(x1, x2) (0x80&(((char*)&x1)[7]^((char*)&x2)[7]))

//---------------------------------------------------------------------------
/**
  function GradientMax: gradient method maximizing F(x)

  In: function F(x) and its derivative dF(x)
      start[dim]: initial value of x
      ep: a stopping threshold
      maxiter: maximal number of iterates
  Out: start[dim]: the final x

  Returnx max F(final x)
*/
double GradientMax(void* params, double (*F)(int, double*, void*), void (*dF)(double*, int, double*, void*),
  int dim, double* start, double ep, int maxiter)
{
  double oldG, G=F(dim, start, params), *dG=new double[dim], *X=new double[dim];

  int iter=0;
  do
  {
    oldG=G;

    dF(dG, dim, start, params);
    Search1Dmaxfrom(params, F, dim, start, dG, 1e-2, 1e-6, 1e6, X, &G, ep);
    memcpy(start, X, sizeof(double)*dim);

    iter++;
  }
  while (iter<maxiter && G>oldG*(1+ep));
  return G;
}//GradientMax

/**
  function GradientOneMax: gradient method maximizing F(x), which searches $dim dimensions one after
  another instead of taking the gradient descent direction.

  In: function F(x) and its derivative dF(x)
      start[dim]: initial value of x
      ep: a stopping threshold
      maxiter: maximal number of iterates
  Out: start[dim]: the final x

  Returnx max F(final x)
*/
double GradientOneMax(void* params, double (*F)(int, double*, void*), void (*dF)(double*, int, double*, void*),
  int dim, double* start, double ep, int maxiter)
{
  double oldG, G=F(dim, start, params), *dG=new double[dim], *X=new double[dim];

  int iter=0;
  do
  {
    oldG=G;
    for (int d=0; d<dim; d++)
    {
      dF(dG, dim, start, params);
      for (int i=0; i<dim; i++) if (d!=i) dG[i]=0;
      if (dG[d]!=0)
      {
        Search1Dmaxfrom(params, F, dim, start, dG, 1e-2, 1e-6, 1e6, X, &G, ep);
        memcpy(start, X, sizeof(double)*dim);
      }
    }
    iter++;
  }
  while (iter<maxiter && G>oldG*(1+ep));
  return G;
}//GradientOneMax

//---------------------------------------------------------------------------
/**
  function mindenormaldouble: returns the minimal denormal double-precision floating point number.
*/
double mindenormaldouble()
{
  double result=0;
  *((unsigned char*)&result)=1;
  return result;
}//mindenormaldouble

/**
  function minnormaldouble: returns the minimal normal double-precision floating point number.
*/
double minnormaldouble()
{
  double result=0;
  ((unsigned char*)&result)[6]=16;
  return result;
}//minnormaldouble

/*
  //function nextdouble: returns the next double value after x in the direction of dir
*/
double nextdouble(double x, double dir)
{
  if (x==0)
  {
    if (dir<0) return -mindenormaldouble();
    else return mindenormaldouble();
  }
  else if ((x>0)^(dir>0)) *((__int64*)&x)-=1;
  else *((__int64*)&x)+=1;
  return x;
}//nextdouble

/**
  function Newton: Newton method for finding the root of y(x)

  In: function y(x) and its derivative y'(x)
      params: additional argument for calling y and y'
      x0: inital value of x
      epx, epy: stopping thresholds on x and y respectively
      maxiter: maximal number of iterates
  Out: x0: the root

  Returns 0 if successful.
*/
double Newton(double& x0, double (*y)(double, void*), double (*y1)(double, void*), void* params, double epx, int maxiter, double epy)
{
  double y0=y(x0, params), x2, y2, dx, dy, x3, y3;
  if (y0==0) return 0;

  int iter=0;
  while (iter<maxiter)
  {
    dy=y1(x0, params);
    if (dy==0) return -1;
    x2=x0-y0/dy;
    y2=y(x2, params);
    while (fabs(y2)>fabs(y0))
    {
      x2=(x2+x0)/2;
      y2=y(x2, params);
    }
    if (y2==0)
    {
      x0=x2;
      return 0;
    }
    dx=fabs(x0-x2);
    if (dx==0)
    {
      if (fabs(y0)<epy) return 0;
      else return -1;
    }
    else if (dx<epx)
    {
      if nsign(y2, y0)
      {
        x0=x2;
        return dx;
      }
      else
      {
        x3=2*x2-x0;
        y3=y(x3, params);
        if (y3==0)
        {
          x0=x3;
          return 0;
        }
        if nsign(y2, y3)
        {
          x0=x2;
          return dx;
        }
        else if (fabs(y3)<fabs(y2))
        {
          x2=x3, y2=y3;
        }
      }
    }
    x0=x2;
    y0=y2;
    iter++;
  }
  if (fabs(y0)<epy) return 0;
  else return -1;
}//Newton

/**
  function Newton: Newton method for finding the root of y(x), for use when it's more efficient for y'
  and y to be calculated in one function call

  In: function dy(x) that computes y and y' at x
      params: additional argument for calling dy
      yshift: the byte shift in params where dy(x) put the value of y(x) as double type
      x0: inital value of x
      epx, epy: stopping thresholds on x and y respectively
      maxiter: maximal number of iterates
  Out: x0: the root

  Returns 0 if successful.
*/
int Newton(double& x0, double (*dy)(double, void*), void* params, int yshift, double epx, int maxiter, double epy, double xmin, double xmax)
{
  double y0, x2, y2, dx, dy0, dy2, dy3, x3, y3;
  dy0=dy(x0, params);
  y0=*(double*)((char*)params+yshift);
  if (y0==0) return 0;

  int iter=0;
  while (iter<maxiter)
  {
    if (dy0==0) return -1;
    dx=-y0/dy0;

    x2=x0+dx;
    dy2=dy(x2, params);
    y2=*(double*)((char*)params+yshift);
    //make sure the iteration yields a y closer to 0
    while (fabs(y2)>fabs(y0))
    {
      dx/=2;
      x2=x0+dx;
      dy2=dy(x2, params);
      y2=*(double*)((char*)params+yshift);
    }
    //the above may fail due to precision problem, i.e. x0+dx=x0
    if (x2==x0)
    {
      if (fabs(y0)<epy) return 0;
      else
      {
        x2=nextdouble(x0, dx);
        dy2=dy(x2, params);
        y2=*(double*)((char*)params+yshift);
        if (y2==0)
        {
          x0=x2;
          return 0;
        }
        else
        {
          dx=x2-x0;
          while (y2==y0)
          {
            dx*=2;
            x2=x0+dx;
            dy2=dy(x2, params);
            y2=*(double*)((char*)params+yshift);
          }
          if nsign(y2, y0)
            return fabs(dx);
          else if (fabs(y2)>fabs(y0))
            return -1;
        }
      }
    }
    else if (y2==0)
    {
      x0=x2;
      return 0;
    }
    else
    {
      double fabsdx=fabs(dx);
      if (fabsdx<epx)
      {
        if nsign(y2, y0)
        {
          x0=x2;
          return fabsdx;
        }
        else
        {
          x3=x2+dx;
          dy3=dy(x3, params);
          y3=*(double*)((char*)params+yshift);
          if (y3==0){x0=x3; return 0;}
          if nsign(y2, y3)
          {
            x0=x2;
            return dx;
          }
          else if (fabs(y3)<fabs(y2))
          {
            x2=x3, y2=y3, dy2=dy3;
          }
        }
      }
    }
    if (x2<xmin || x2>xmax) return -1;
    x0=x2;
    y0=y2;
    dy0=dy2;
    iter++;
  }
  if (fabs(y0)<epy) return 0;
  else return -1;
}//Newton

/**
  function init_Newton_1d: finds an initial x and interval (xa, xb) from which to search for an extremum

  On calling, xa<x0<xb, y(x0)>=y(xa), y(x0)>=y(xb) (or y(x0)<=y(xa), y(x0)<=y(xb)).
  On return, y(x0)>=y(xa), y(x0)>=y(xb) (or y(x0)<=y(xa), y(x0)<=y(xb)), and either y'(xa)>0, y'(xb)<0
  (or y'(xa)<0, y'(xb)>0), or xb-xa<epx.

  In: function dy(x) that computes y and y' at x
      params: additional argument for calling dy
      yshift: the byte shift in params where dy(x) put the value of y(x) as double type
      x0: inital value of x
      (xa, xb): initial search range
      epx: tolerable error of extremum
  Out: improved initial value x0 and search range (xa, xb)

  Returns xb-xa if successful, -0.5 or -0.6 if y(x0) is not bigger(smaller) than both y(xa) and y(xb)
*/
double init_Newton_1d(double& x0, double& xa, double& xb, double (*dy)(double, void*), void* params, int yshift, double epx)
{
  double dy0, dya, dyb, y0, ya, yb;
  dy0=dy(x0, params); y0=*(double*)((char*)params+yshift);
  dya=dy(xa, params); ya=*(double*)((char*)params+yshift);
  dyb=dy(xb, params); yb=*(double*)((char*)params+yshift);
  if (y0>=ya && y0>=yb)
  {
    //look for xa<x0<xb, ya<y0>yb, so that dya>0, dyb<0
    while ((dya<=0 || dyb>=0) && xb-xa>=epx)
    {
      double x1=(x0+xa)/2;
      double dy1=dy(x1, params);
      double y1=*(double*)((char*)params+yshift);
      double x2=(x0+xb)/2;
      double dy2=dy(x2, params);
      double y2=*(double*)((char*)params+yshift);
      if (y0>=y1 && y0>=y2) xa=x1, ya=y1, dya=dy1, xb=x2, yb=y2, dyb=dy2;
      else if (y1>=y0 && y1>=y2) {xb=x0, yb=y0, dyb=dy0; x0=x1, y0=y1, dy0=dy1;}
      else {xa=x0, ya=y0, dya=dy0; x0=x2, y0=y2, dy0=dy2;}
    }
    return xb-xa;
  }
  else if (y0<=ya && y0<=yb)
  {
    //look for xa<x0<xb, ya>y0<yb, so that dya<0, dyb>0
    while ((dya>=0 || dyb<=0) && xb-xa>=epx)
    {
      double x1=(x0+xa)/2;
      double dy1=dy(x1, params);
      double y1=*(double*)((char*)params+yshift);
      double x2=(x0+xb)/2;
      double dy2=dy(x2, params);
      double y2=*(double*)((char*)params+yshift);
      if (y0<=y1 && y0<=y2) xa=x1, ya=y1, dya=dy1, xb=x2, yb=y2, dyb=dy2;
      else if (y1<=y0 && y1<=y2) {xb=x0, yb=y0, dyb=dy0; x0=x1, y0=y1, dy0=dy1;}
      else {xa=x0, ya=y0, dya=dy0; x0=x2, y0=y2, dy0=dy2;}
    }
    return xb-xa;
  }
  else//meaning y(x0) is not bigger(smaller) than both y(xa) and y(xb)
  {
    if (y0<=yb) return -0.5; //ya<=y0<=yb
    if (y0<=ya) return -0.6; //ya>=y0>=yb
  }
  return -0;
}//init_Newton_1d

/**
  function init_Newton_1d_max: finds an initial x and interval (xa, xb) from which to search for an
  maximum.

  On calling xa<x0<xb, and it is assumed that y(x0)>=y(xa), y(x0)>=y(xb). If the conditions on y are met,
  it find a new set of xa<x0<xb, so that y(x0)>y(xa), y(x0)>=y(xb), y'(xa)>0, y'(xb)<0, or xb-xa<epx,
  and returns xb-xa>0; otherwise it returns a negative value showing the initial order of y(x0), y(xa)
  and y(xb) as
    -0.5 for ya<=y0<yb
    -0.6 for yb<=y0<ya
    -0.7 for y0<ya<yb
    -0.8 for y0<yb<=ya

  In: function dy(x) that computes y and y' at x
      params: additional argument for calling dy
      yshift: the byte shift in params where dy(x) put the value of y(x) as double type
      x0: inital value of x
      (xa, xb): initial search range
      epx: tolerable error of extremum
  Out: improved initial value x0 and search range (xa, xb)

  Returns xb-xa if successful, a negative value as above if not.
*/
double init_Newton_1d_max(double& x0, double& xa, double& xb, double (*dy)(double, void*), void* params, int yshift, double epx)
{
  double dy0, dya, dyb, y0, ya, yb;
  dy0=dy(x0, params); y0=*(double*)((char*)params+yshift);
  dya=dy(xa, params); ya=*(double*)((char*)params+yshift);
  dyb=dy(xb, params); yb=*(double*)((char*)params+yshift);
  if (y0>=ya && y0>=yb)
  {
    //look for xa<x0<xb, ya<y0>yb, so that dya>0, dyb<0
    while ((dya<=0 || dyb>=0) && xb-xa>=epx)
    {
      double x1=(x0+xa)/2;
      double dy1=dy(x1, params);
      double y1=*(double*)((char*)params+yshift);
      double x2=(x0+xb)/2;
      double dy2=dy(x2, params);
      double y2=*(double*)((char*)params+yshift);
      if (y0>=y1 && y0>=y2) xa=x1, ya=y1, dya=dy1, xb=x2, yb=y2, dyb=dy2;
      else if (y1>=y0 && y1>=y2) {xb=x0, yb=y0, dyb=dy0; x0=x1, y0=y1, dy0=dy1;}
      else {xa=x0, ya=y0, dya=dy0; x0=x2, y0=y2, dy0=dy2;}
    }
    return xb-xa;
  }
  else//meaning y(x0) is not bigger than both y(xa) and y(xb)
  {
    if (y0<yb && y0>=ya) return -0.5; //ya<=y0<yb
    if (y0<ya && y0>=yb) return -0.6; //ya>y0>=yb
    if (ya<yb) return -0.7; //y0<ya<yb
    return -0.8; //y0<yb<=ya
  }
}//init_Newton_1d_max

/**
  function init_Newton_1d_min: finds an initial x and interval (xa, xb) from which to search for an
  manimum.

  On calling xa<x0<xb, and it is assumed that y(x0)<=y(xa), y(x0)<=y(xb). If the conditions on y are met,
  it find a new set of xa<x0<xb, so that y(x0)<y(xa), y(x0)<=y(xb), y'(xa)<0, y'(xb)>0, or xb-xa<epx,
  and returns xb-xa>0; otherwise it returns a negative value showing the initial order of y(x0), y(xa)
  and y(xb) as
    -0.5 for ya<=y0<yb
    -0.6 for yb<=y0<ya
    -0.7 for y0>yb>ya
    -0.8 for y0>ya>=yb

  In: function dy(x) that computes y and y' at x
      params: additional argument for calling dy
      yshift: the byte shift in params where dy(x) put the value of y(x) as double type
      x0: inital value of x
      (xa, xb): initial search range
      epx: tolerable error of extremum
  Out: improved initial value x0 and search range (xa, xb)

  Returns xb-xa if successful, a negative value as above if not.
*/
double init_Newton_1d_min(double& x0, double& xa, double& xb, double (*dy)(double, void*), void* params, int yshift, double epx)
{
  double dy0, dya, dyb, y0, ya, yb;
  dy0=dy(x0, params); y0=*(double*)((char*)params+yshift);
  dya=dy(xa, params); ya=*(double*)((char*)params+yshift);
  dyb=dy(xb, params); yb=*(double*)((char*)params+yshift);
  if (y0<=ya && y0<=yb)
  {
    //look for xa<x0<xb, ya>y0<yb, so that dya<0, dyb>0
    while ((dya>=0 || dyb<=0) && xb-xa>=epx)
    {
      double x1=(x0+xa)/2;
      double dy1=dy(x1, params);
      double y1=*(double*)((char*)params+yshift);
      double x2=(x0+xb)/2;
      double dy2=dy(x2, params);
      double y2=*(double*)((char*)params+yshift);
      if (y0<=y1 && y0<=y2) xa=x1, ya=y1, dya=dy1, xb=x2, yb=y2, dyb=dy2;
      else if (y1<=y0 && y1<=y2) {xb=x0, yb=y0, dyb=dy0; x0=x1, y0=y1, dy0=dy1;}
      else {xa=x0, ya=y0, dya=dy0; x0=x2, y0=y2, dy0=dy2;}
    }
    return xb-xa;
  }
  else//meaning y(x0) is not bigger than both y(xa) and y(xb)
  {
    if (y0<yb && y0>=ya) return -0.5; //ya<=y0<yb
    if (y0<ya && y0>=yb) return -0.6; //ya>y0>=yb
    if (ya<yb) return -0.7; //ya<yb<y0
    return -0.8; //yb<=ya<y0
  }
}//init_Newton_1d_min

/**
  function Newton_1d_max: Newton method for maximizing y(x). On calling y(x0)>=y(xa), y(x0)>=y(xb),
  y'(xa)>0, y'(xb)<0

  In: function ddy(x) that computes y, y' and y" at x
      params: additional argument for calling ddy
      y1shift: the byte shift in params where ddy(x) put the value of y'(x) as double type
      yshift: the byte shift in params where ddy(x) put the value of y(x) as double type
      x0: inital value of x
      (xa, xb): initial search range
      epx, epdy: stopping threshold for x and y', respectively
      maxiter: maximal number of iterates
  Out: x0: the maximum

  Returns the error bound for the maximum x0, or
       -2 if y(*) do not satisfy input conditions
       -3 if y'(*) do not satisfy input conditions
*/
double Newton_1d_max(double& x0, double xa, double xb, double (*ddy)(double, void*), void* params, int y1shift, int yshift, double epx, int maxiter, double epdy, bool checkinput)
{
  double x1, y1, x2, y2, dx, dy1, dy2, dy3, x3, y3, ddy1, ddy2, ddy3;
  double y0, ya, yb, dy0, dya, dyb, ddy0;
  if (xb-xa<epx) {return xb-xa;}
  ddy(xa, params);
  dya=*(double*)((char*)params+y1shift), ya=*(double*)((char*)params+yshift);
  ddy(xb, params);
  dyb=*(double*)((char*)params+y1shift), yb=*(double*)((char*)params+yshift);
  ddy0=ddy(x0, params);
  dy0=*(double*)((char*)params+y1shift), y0=*(double*)((char*)params+yshift);

  if (checkinput)
  {
    if (y0>=ya && y0>=yb)
    {
      if (dya<=0 || dyb>=0) return -3;
    }
    else return -2;
  }

  if (dy0==0) return 0;

  int iter=0;
  while (iter<maxiter)
  {
    if (ddy0>=0) //newton method not applicable
    {
      do
      {
        x1=(x0+xa)/2;
        ddy1=ddy(x1, params);
        dy1=*(double*)((char*)params+y1shift), y1=*(double*)((char*)params+yshift);
        x2=(x0+xb)/2;
        ddy2=ddy(x2, params);
        dy2=*(double*)((char*)params+y1shift), y2=*(double*)((char*)params+yshift);
        if (y0>=y1 && y0>=y2) xa=x1, ya=y1, dya=dy1, xb=x2, yb=y2, dyb=dy2;
        else if (y1>=y0 && y1>=y2) {xb=x0, yb=y0, dyb=dy0; x0=x1, y0=y1, dy0=dy1, ddy0=ddy1;}
        else {xa=x0, ya=y0, dya=dy0; x0=x2, y0=y2, dy0=dy2, ddy0=ddy2;}
      }
      while(ddy0>=0 && (dya<=0 || dyb>=0) && xb-xa>=epx);
      if (xb-xa<epx) return xb-xa;
    }
    if (fabs(dy0)<epdy) return 0;
    //x2: Newton point
    //dx is the same direction as dy0 if ddy0<0
    dx=-dy0/ddy0;
    x2=x0+dx;
    if (x2==x0) return 0;

    //make sure x2 lies in [xa, xb]
    if (x2<xa) {x2=(xa+x0)/2; dx=x2-x0;}
    else if (x2>xb) {x2=(xb+x0)/2; dx=x2-x0;}
    ddy2=ddy(x2, params);
    y2=*(double*)((char*)params+yshift);

    //make sure the iteration yields a larger y
    if (y2<y0)
    {
      while (y2<y0 && dx>epx)
      {
        dx/=2;
        x2=x0+dx;
        ddy2=ddy(x2, params);
        y2=*(double*)((char*)params+yshift);
      }
      if (y2<y0) return fabs(dx);
    }

    //new x2 lies between xa and xb and y2>=y0
    dy2=*(double*)((char*)params+y1shift);
    if (fabs(dy2)<epdy){x0=x2; return 0;}

    double fabsdx=fabs(dx);
    if (fabsdx<epx)
    {
      if nsign(dy0, dy2){x0=x2; return fabsdx;}
      x3=x2+dx;
      ddy3=ddy(x3, params);
      dy3=*(double*)((char*)params+y1shift), y3=*(double*)((char*)params+yshift);
      if (dy3==0) {x0=x3; return 0;}
      if (y2>=y3) {x0=x2; return fabsdx;}
      else x2=x3, y2=y3, dy2=dy3, ddy2=ddy3;
    }

    if (x0<x2) xa=x0, ya=y0, dya=dy0;
    else xb=x0, yb=y0, dyb=dy0;
    x0=x2, y0=y2, dy0=dy2, ddy0=ddy2;
    iter++;
  }
  if (fabs(dy0)<epdy) return 0;
  else return -1;
}//Newton_1d_max

/**
  function Newton_1d_min: Newton method for minimizing y(x). On calling y(x0)<=y(xa), y(x0)<=y(xb),
  y'(xa)<0, y'(xb)>0

  In: function ddy(x) that computes y, y' and y" at x
      params: additional argument for calling ddy
      y1shift: the byte shift in params where ddy(x) put the value of y'(x) as double type
      yshift: the byte shift in params where ddy(x) put the value of y(x) as double type
      x0: inital value of x
      (xa, xb): initial search range
      epx, epdy: stopping threshold for x and y', respectively
      maxiter: maximal number of iterates
  Out: x0: the maximum

  Returns the error bound for the maximum x0, or
       -2 if y(*) do not satisfy input conditions
       -3 if y'(*) do not satisfy input conditions
*/
double Newton_1d_min(double& x0, double xa, double xb, double (*ddy)(double, void*), void* params, int y1shift, int yshift, double epx, int maxiter, double epdy, bool checkinput)
{
  double x1, y1, x2, y2, dx, dy1, dy2, dy3, x3, y3, ddy1, ddy2, ddy3;
  double y0, ya, yb, dy0, dya, dyb, ddy0;
  if (xb-xa<epx) {return xb-xa;}
  ddy(xa, params);
  dya=*(double*)((char*)params+y1shift), ya=*(double*)((char*)params+yshift);
  ddy(xb, params);
  dyb=*(double*)((char*)params+y1shift), yb=*(double*)((char*)params+yshift);
  ddy0=ddy(x0, params);
  dy0=*(double*)((char*)params+y1shift), y0=*(double*)((char*)params+yshift);

  if (checkinput)
  {
    if (y0<=ya && y0<=yb)
    {
      if (dya>=0 || dyb<=0) return -3;
    }
    else return -2;
  }

  if (dy0==0) return 0;

  int iter=0;
  while (iter<maxiter)
  {
    if (ddy0<=0) //newton method not applicable
    {
      do
      {
        x1=(x0+xa)/2;
        ddy1=ddy(x1, params);
        dy1=*(double*)((char*)params+y1shift), y1=*(double*)((char*)params+yshift);
        x2=(x0+xb)/2;
        ddy2=ddy(x2, params);
        dy2=*(double*)((char*)params+y1shift), y2=*(double*)((char*)params+yshift);
        if (y0<=y1 && y0<=y2) xa=x1, ya=y1, dya=dy1, xb=x2, yb=y2, dyb=dy2;
        else if (y1<=y0 && y1<=y2) {xb=x0, yb=y0, dyb=dy0; x0=x1, y0=y1, dy0=dy1, ddy0=ddy1;}
        else {xa=x0, ya=y0, dya=dy0; x0=x2, y0=y2, dy0=dy2, ddy0=ddy2;}
      }
      while((dya>=0 || dyb<=0) && xb-xa>=epx);
    }
    else
    {
      //Newton point
      //dx is the same direction as dy0 if ddy0<0
      if (fabs(dy0)<epdy) return 0;
      dx=-dy0/ddy0;
      x2=x0+dx;
      if (x2==x0) return 0;
      //make sure x2 lies in [xa, xb]
      if (x2<xa) {x2=(xa+x0)/2; dx=x2-x0;}
      else if (x2>xb) {x2=(xb+x0)/2; dx=x2-x0;}
      ddy2=ddy(x2, params);
      y2=*(double*)((char*)params+yshift);
      //make sure the iteration yields a larger y
      if (y2>y0)
      {
        while (y2>y0)
        {
          dx/=2;
          x2=x0+dx;
          ddy2=ddy(x2, params);
          y2=*(double*)((char*)params+yshift);
        }
        if (x2==x0) return 0; //dy0<0 but all y(x+(ep)dx)<0
      }
      dy2=*(double*)((char*)params+y1shift);
      if (fabs(dy2)<epdy)
      {
        x0=x2;
        return 0;
      }
      else
      {
        double fabsdx=fabs(dx);
        if (fabsdx<epx)
        {
          if nsign(dy2, dy0)
          {
            x0=x2;
            return fabsdx;
          }
          else
          {
            x3=x2+dx;
            ddy3=ddy(x3, params);
            dy3=*(double*)((char*)params+y1shift), y3=*(double*)((char*)params+yshift);
            if (dy3==0) {x0=x3; return 0;}
            if nsign(dy2, dy3)
            {
              x0=x2;
              return dx;
            }
            else if (y3<y2)
            {
              x2=x3, y2=y3, dy2=dy3, ddy2=ddy3;
            }
          }
        }
      }
      if (x0<x2) xa=x0, ya=y0, dya=dy0;
      else xb=x0, yb=y0, dyb=dy0;

      x0=x2, y0=y2, dy0=dy2, ddy0=ddy2;
    }

    iter++;
  }
  if (fabs(dy0)<epdy) return 0;
  else return -1;
}//Newton_1d_min

/**
  function Newton1dmax: uses init_Newton_1d_max and Newton_1d_max to search for a local maximum of y.

  By default ddy and dy use the same parameter structure and returns the lower order derivatives at
  specific offsets within that structure. However, since the same params is used, it's possible the same
  variable may take different offsets when returned from ddy() and dy().

  On calling xa<x0<xb, and it is assumed that y0>=ya, y0>=yb. If the conditions on y are not met, it
  returns a negative value showing the order of y0, ya and yb, i.e.
    -0.5 for ya<=y0<yb
    -0.6 for yb<=y0<ya
    -0.7 for y0<ya<yb
    -0.8 for y0<yb<=ya

  In: function ddy(x) that computes y, y' and y" at x
      function dy(x) that computes y and y' at x
      params: additional argument for calling ddy and dy
      y1shift: the byte shift in params where ddy(x) put the value of y'(x) as double type
      yshift: the byte shift in params where ddy(x) put the value of y(x) as double type
      dyshift: the byte shift in params where dy(x) put the value of y(x) as double type
      x0: inital value of x
      (xa, xb): initial search range
      epx, epdy: stopping threshold for x and y', respectively
      maxiter: maximal number of iterates
  Out: x0: the maximum

  Returns the error bound for maximum x0, or a negative value if the initial point does not satisfy
  relevent assumptions.
*/
double Newton1dmax(double& x0, double xa, double xb, double (*ddy)(double, void*), void* params, int y1shift, int yshift, double (*dy)(double, void*), int dyshift, double epx, int maxiter, double epdy, bool checkinput)
{
  double ep=init_Newton_1d_max(x0, xa, xb, dy, params, dyshift, epx);
  if (ep>epx) ep=Newton_1d_max(x0, xa, xb, ddy, params, y1shift, yshift, epx, maxiter, epdy, true);
  return ep;
}//Newton1dmax

//function Newton1dmin: same as Newton1dmax, in the other direction.
double Newton1dmin(double& x0, double xa, double xb, double (*ddy)(double, void*), void* params, int y1shift, int yshift, double (*dy)(double, void*), int dyshift, double epx, int maxiter, double epdy, bool checkinput)
{
  double ep=init_Newton_1d_min(x0, xa, xb, dy, params, dyshift, epx);
  if (ep>epx) ep=Newton_1d_min(x0, xa, xb, ddy, params, y1shift, yshift, epx, maxiter, epdy, true);
  return ep;
}//Newton1dmin

/**
  function NewtonMax: looks for the maximum of y within [x0, x1] using Newton method.

  In: function ddy(x) that computes y, y' and y" at x
      params: additional argument for calling ddy
      y1shift: the byte shift in params where ddy(x) put the value of y'(x) as double type
      yshift: the byte shift in params where ddy(x) put the value of y(x) as double type
      (x0, x1): initial search range
      epx, epy: stopping threshold for x and y, respectively
      maxiter: maximal number of iterates
  Out: maximum x0

  Returns an error bound on maximum x0 if successful, a negative value if not.
*/
double NewtonMax(double& x0, double x1, double (*ddy)(double, void*), void* params, int y1shift, int yshift, double epx, int maxiter, double epy)
{
  double y0, x2, y2, dx, dy0, dy2, dy3, x3, y3;
  dy0=ddy(x0, params);
  y0=*(double*)((char*)params+yshift);
  if (y0==0) return 0;

  int iter=0;
  while (iter<maxiter)
  {
    if (dy0==0) return -1;
    dx=-y0/dy0;

    x2=x0+dx;
    dy2=ddy(x2, params);
    y2=*(double*)((char*)params+yshift);
    //make sure the iteration yields a y closer to 0
    while (fabs(y2)>fabs(y0))
    {
      dx/=2;
      x2=x0+dx;
      dy2=ddy(x2, params);
      y2=*(double*)((char*)params+yshift);
    }
    //the above may fail due to precision problem, i.e. x0+dx=x0
    if (x2==x0)
    {
      if (fabs(y0)<epy) return 0;
      else
      {
        x2=nextdouble(x0, dx);
        dy2=ddy(x2, params);
        y2=*(double*)((char*)params+yshift);
        if (y2==0)
        {
          x0=x2;
          return 0;
        }
        else
        {
          dx=x2-x0;
          while (y2==y0)
          {
            dx*=2;
            x2=x0+dx;
            dy2=ddy(x2, params);
            y2=*(double*)((char*)params+yshift);
          }
          if nsign(y2, y0)
            return fabs(dx);
          else if (fabs(y2)>fabs(y0))
            return -1;
        }
      }
    }
    else if (y2==0)
    {
      x0=x2;
      return 0;
    }
    else
    {
      double fabsdx=fabs(dx);
      if (fabsdx<epx)
      {
        if nsign(y2, y0)
        {
          x0=x2;
          return fabsdx;
        }
        else
        {
          x3=x2+dx;
          dy3=ddy(x3, params);
          y3=*(double*)((char*)params+yshift);
          if (y3==0){x0=x3; return 0;}
          if nsign(y2, y3)
          {
            x0=x2;
            return dx;
          }
          else if (fabs(y3)<fabs(y2))
          {
            x2=x3, y2=y3, dy2=dy3;
          }
        }
      }
    }
    x0=x2;
    y0=y2;
    dy0=dy2;
    iter++;
  }
  if (fabs(y0)<epy) return 0;
  else return -1;
}//NewtonMax

//---------------------------------------------------------------------------
/**
  function rootsearchhalf: searches for the root of y(x) in (x1, x2) by half-split method. On calling
  y(x1) and y(x2) must not have the same sign, or -1 is returned signaling a failure.

  In: function y(x),
      params: additional arguments for calling y
      (x1, x2): inital range
      epx: stopping threshold on x
      maxiter: maximal number of iterates
  Out: root x1

  Returns a positive error estimate if a root is found, -1 if not.
*/
double rootsearchhalf(double&x1, double x2, double (*y)(double, void*), void* params, double epx, int maxiter)
{
  if (x1>x2) {double tmp=x1; x1=x2; x2=tmp;}
  double y1=y(x1, params), y3, x3, ep2=epx*2;
  if (y1==0) return 0;
  else
  {
    double y2=y(x2, params);
    if (y2==0)
    {
      x1=x2;
      return 0;
    }
    else if nsign(y1, y2) //y1 and y2 have different signs
    {
      int iter=0;
      while (x2-x1>ep2 && iter<maxiter)
      {
        x3=(x1+x2)/2;
        y3=y(x3, params);
        if (y3==0)
        {
          x1=x3;
          return 0;
        }
        else
        {
          if nsign(y1, y3) x2=x3, y2=y3;
          else x1=x3, y1=y3;
        }
        iter++;
      }
      x1=(x1+x2)/2;
      return x2-x1;
    }
    else
      return -1;
  }
}//rootsearchhalf

//---------------------------------------------------------------------------
/**
  function rootsearchsecant: searches for the root of y(x) in (x1, x2) by secant method.

  In: function y(x),
      params: additional arguments for calling y
      (x1, x2): inital range
      epx: stopping threshold on x
      maxiter: maximal number of iterates
  Out: root x1

  Returns a positive error estimate if a root is found, -1 if not.
*/
double rootsearchsecant(double& x1, double x2, double (*y)(double, void*), void* params, double epx, int maxiter)
{
  double y1=y(x1, params), y2=y(x2, params), y3, x3, dx, x0, y0;
  if (y1==0) return 0;
  if (y2==0)
  {
    x1=x2;
    return 0;
  }
  if (fabs(y1)>fabs(y2))
  {
    x3=x1, y3=y1;
  }
  else
  {
    x3=x2, y3=y2;
    x2=x1, y2=y1;
  }

  int iter=0;
  while (iter<maxiter)
  {
    x1=(y3*x2-y2*x3)/(y3-y2);
    y1=y(x1, params);

    dx=fabs(x1-x2);
    if (dx<epx)
    {
      if nsign(y1, y2)
      {
        return dx;
      }
      else
      {
        x0=2*x1-x2;
        y0=y(x0, params);
        if nsign(y0, y1)
        {
          return dx;
        }
      }
    }
    x3=x2, y3=y2;
    x2=x1, y2=y1;
    iter++;
  }
  return -1;
}//rootsearchsecant

/**
  function rootsearchbisecant: searches for the root of y(x) in (x1, x2) by bi-secant method. On calling
  given that y(x1) and y(x2) have different signs. The bound interval (x1, x2) converges to 0, as
  compared to standared secant method.

  In: function y(x),
      params: additional arguments for calling y
      (x1, x2): inital range
      epx: stopping threshold on x
      maxiter: maximal number of iterates
  Out: root x1

  Returns an estimated error bound if a root is found, a negative value if not.
*/
double rootsearchbisecant(double&x1, double x2, double (*y)(double, void*), void* params, double epx, int maxiter)
{
  if (x1>x2) {double tmp=x1; x1=x2; x2=tmp;}
  double y1=y(x1, params), y3, y4, dy2, dy3, x3, x4, ep2=epx*2;
  if (y1==0) return 0;
  else
  {
    double y2=y(x2, params);
    if (y2==0)
    {
      x1=x2;
      return 0;
    }
    else if nsign(y1, y2) //y1 and y2 have different signs
    {
      int iter=0;
      while (x2-x1>ep2 && iter<maxiter)
      {
        dy2=y1-y2;
        x3=(y1*x2-y2*x1)/dy2;
        y3=y(x3, params);
        if (y3==0)
        {
          x1=x3;
          return 0;
        }
        else
        {
          if nsign(y1, y3)
          {
            dy3=y3-y2;
            if (dy3==0 || nsign(dy3, dy2)) x2=x3, y2=y3;
            else
            {
              x4=(y3*x2-y2*x3)/dy3; //x4<x3
              if (x4<=x1) x2=x3, y2=y3;
              else
              {
                y4=y(x4, params);
                if (y4==0)
                {
                  x1=x4;
                  return 0;
                }
                else
                {
                  if nsign(y4, y3) x1=x4, y1=y4, x2=x3, y2=y3;
                  else x2=x4, y2=y4;
                }
              }
            }
          }
          else
          {
            dy3=y1-y3;
            if (dy3==0 || nsign(dy3, dy2)) x1=x3, y1=y3;
            else
            {
              x4=(y1*x3-y3*x1)/dy3; //x4>x3
              if (x4>=x2) x1=x3, y1=y3;
              else
              {
                y4=y(x4, params);
                if (y4==0)
                {
                  x1=x4;
                  return 0;
                }
                else
                {
                  if nsign(y4, y3) x1=x3, y1=y3, x2=x4, y2=y4;
                  else x1=x4, y1=y4;
                }
              }
            }
          }
        }
        iter++;
      }
      x1=(x1+x2)/2;
      return x2-x1;
    }
    else
      return -1;
  }
}//rootsearchbisecant

//---------------------------------------------------------------------------
/**
  function Search1Dmax: 1-dimensional maximum search within interval [inf, sup] using 0.618 method

  In: function F(x),
      params: additional arguments for calling F
      (inf, sup): inital range
      epx: stopping threshold on x
      maxiter: maximal number of iterates
      convcheck: specifies if convexity is checked at each interate, false recommended.
  Out: x: the maximum
       maxresult: the maximal value, if specified

  Returns an estimated error bound of x if successful, a negative value if not.
*/
double Search1Dmax(double& x, void* params, double (*F)(double, void*), double inf, double sup, double* maxresult, double ep, int maxiter, bool convcheck)
{
  double v0618=(sqrt(5.0)-1)/2, iv0618=1-v0618;
  double l=sup-inf;
  double startsup=sup, startinf=inf, startl=l;
  double x1=inf+l*(1-v0618);
  double x2=inf+l*v0618;
  double result;

  double v0=F(inf, params);
  double v1=F(x1, params);
  double v2=F(x2, params);
  double v3=F(sup, params);
  int iter=0;

  if (v1<v0 && v1<v3 && v2<v0 && v2<v3) goto return1;
  if (convcheck && (v1<v0*v0618+v3*iv0618 || v2<v0*iv0618+v3*v0618 || v1<v0*iv0618+v2*v0618 || v2<v1*v0618+v3*iv0618)) goto return1;

  while (iter<maxiter)
  {
    if (l<startl*ep) break;
    if (v1>v2)
    {
      sup=x2, v3=v2;
      l=sup-inf;
      x2=x1, v2=v1;
      x1=inf+l*(1-v0618);
      v1=F(x1, params);
      if (convcheck && (v1<v0*v0618+v3*iv0618 || v1<v0*iv0618+v2*v0618 || v2<v1*v0618+v3*iv0618)) goto return1;
    }
    else
    {
      inf=x1, v0=v1;
      l=sup-inf;
      x1=x2, v1=v2;
      x2=inf+l*v0618;
      v2=F(x2, params);
      if (convcheck && (v2<v0*iv0618+v3*v0618 || v1<v0*iv0618+v2*v0618 || v2<v1*v0618+v3*iv0618)) goto return1;
    }
    iter++;
  }
  if (v1>v2) {x=x1; if (maxresult) maxresult[0]=v1; result=x1-inf;}
  else {x=x2; if (maxresult) maxresult[0]=v2; result=sup-x2;}
  return result;

return1: //where at some stage v(.) fails to be convex
  if (v0>=v1 && v0>=v2 && v0>=v3){x=inf; if (maxresult) maxresult[0]=v0; result=(x==startinf)?(-1):(sup-x1);}
  else if (v1>=v0 && v1>=v2 && v1>=v3){x=x1; if (maxresult) maxresult[0]=v1; result=x1-inf;}
  else if (v2>=v0 && v2>=v1 && v2>=v3){x=x2; if (maxresult) maxresult[0]=v2; result=sup-x2;}
  else {x=sup; if (maxresult) maxresult[0]=v3; result=(x==startsup)?(-1):(x2-inf);}
  return result;
}//Search1Dmax

/**
  function Search1DmaxEx: 1-dimensional maximum search using Search1Dmax, but allows setting a $len
  argument so that a pre-location of the maximum is performed by sampling (inf, sup) at interval $len,
  so that Search1Dmax will work on an interval no longer than 2*len.

  In: function F(x),
      params: additional arguments for calling F
      (inf, sup): inital range
      len: rough sampling interval for pre-locating maximum
      epx: stopping threshold on x
      maxiter: maximal number of iterates
      convcheck: specifies if convexity is checked at each interate, false recommended.
  Out: x: the maximum
       maxresult: the maximal value, if specified

  Returns an estimated error bound of x if successful, a negative value if not.
*/
double Search1DmaxEx(double& x, void* params, double (*F)(double, void*), double inf, double sup, double* maxresult, double ep, int maxiter, double len)
{
  int c=ceil((sup-inf)/len+1);
  if (c<8) c=8;
  double* vs=new double[c+1], step=(sup-inf)/c;
  for (int i=0; i<=c; i++) vs[i]=F(inf+step*i, params);
  int max=0; for (int i=1; i<=c; i++) if (vs[max]<vs[i]) max=i;
  *maxresult=vs[max]; delete[] vs;
  if (max>0 && max<c) return Search1Dmax(x, params, F, inf+step*(max-1), inf+step*(max+1), maxresult, ep, maxiter);
  else if (max==0)
  {
    if (step<=ep) {x=inf; return ep;}
    else return Search1DmaxEx(x, params, F, inf, inf+step, maxresult, ep, maxiter, len);
  }
  else
  {
    if (step<=ep) {x=sup; return ep;}
    else return Search1DmaxEx(x, params, F, sup-step, sup, maxresult, ep, maxiter, len);
  }
}//Search1DmaxEx

/**
  function Search1Dmaxfrom: 1-dimensional maximum search from start[] towards start+direct[]. The
  function tries to locate an interval ("step") for which the 4-point convexity rule confirms the
  existence of a convex maximum.

  In: function F(x)
      params: additional arguments for calling F
      start[dim]: initial x[]
      direct[dim]: search direction
      step: initial step size
      minstep: minimal step size in 4-point scheme between x0 and x3
      maxstep: maximal step size in 4-point scheme between x0 and x3
  Out: opt[dim]: the maximum, if specified
       max: maximal value, if specified

  Returns x>0 that locally maximizes F(start+x*direct)
    or x|-x>0 that maximizes F(start-x*direct) at calculated points if convex check fails.
*/
double Search1Dmaxfrom(void* params, double (*F)(int, double*, void*), int dim, double* start, double* direct,
  double step, double minstep, double maxstep, double* opt, double* max, double ep, int maxiter, bool convcheck)
{
  bool record=true;
  if (!opt)
  {
    record=false;
    opt=new double[dim];
  }
  
  double v0618=(sqrt(5.0)-1)/2, iv0618=1-v0618;
  double x3=step, x0=0, x1=-1, x2=-1, l=x3-x0;

  //prepares x0~x3, v0~v3, x0=0, v2>v3
  memcpy(opt, start, sizeof(double)*dim);
  double v0=F(dim, opt, params);
  MultiAdd(dim, opt, start, direct, x3);
  double v1, v2, v3=F(dim, opt, params);
  if (v0>=v3)
  {
    x1=x0+l*iv0618;
    MultiAdd(dim, opt, start, direct, x1);
    v1=F(dim, opt, params);
    x2=x0+l*v0618;
    MultiAdd(dim, opt, start, direct, x2);
    v2=F(dim, opt, params);
    while (l>minstep && v0>v1)
    {
      x3=x2, v3=v2;
      x2=x1, v2=v1;
      l=x3-x0;
      x1=x0+l*iv0618;
      MultiAdd(dim, opt, start, direct, x1);
      v1=F(dim, opt, params);
    }
    if (l<=minstep)
    {
      if (max) max[0]=v0;
      if (!record) delete[] opt;
      return x0;
    }
  }
  else//v0<=v3
  {
    do
    {
      x1=x2, v1=v2;
      x2=x3, v2=v3;
      x3+=l*v0618;
      l=x3-x0;
      MultiAdd(dim, opt, start, direct, x3);
      v3=F(dim, opt, params);
    } while (v2<v3 && x3-x0<maxstep);
    if (l>=maxstep) //increasing all the way to maxstep
    {
      if (max) max[0]=v3;
      if (!record) delete opt;
      return x3;
    }
    else if (x1<0)
    {
      x1=x0+l*iv0618;
      MultiAdd(dim, opt, start, direct, x1);
      v1=F(dim, opt, params);
    }
  }

  double oldl=l;

  int iter=0;

  if (convcheck && (v1<v0*v0618+v3*iv0618 || v2<v0*iv0618+v3*v0618
    || v1<v0*iv0618+v2*v0618 || v2<v1*v0618+v3*iv0618)) goto return1;

  while (iter<maxiter)
  {
    if (l<oldl*ep) break;
    if (v1>v2)
    {
      x3=x2, v3=v2;
      l=x3-x0;
      x2=x1, v2=v1;
      x1=x0+l*iv0618;
      MultiAdd(dim, opt, start, direct, x1);
      v1=F(dim, opt, params);
      if (convcheck && (v1<v0*v0618+v3*iv0618 || v1<v0*iv0618+v2*v0618 || v2<v1*v0618+v3*iv0618)) goto return1;
    }
    else
    {
      x0=x1, v0=v1;
      l=x3-x0;
      x1=x2, v1=v2;
      x2=x0+l*v0618;
      MultiAdd(dim, opt, start, direct, x2);
      v2=F(dim, opt, params);
      if (convcheck && (v2<v0*iv0618+v3*v0618 || v1<v0*iv0618+v2*v0618 || v2<v1*v0618+v3*iv0618)) goto return1;
    }
    iter++;
  }

  if (!record) delete[] opt;

  if (max) max[0]=v1;
  return x1;

return1:
  if (!record) delete[] opt;
  if (v1>v0) v0=v1,x0=x1; if (v2>v0) v0=v2, x0=x2; if (v3>v0) v0=v3, x0=x3;
  if (max) max[0]=v0;
  return -x0;
}//Search1Dmaxfrom