lbfgs.cpp

pocket_crf_0.45

#include <math.h>
#include <string.h>
#ifdef _WIN32
#include <io.h>
#else
#include <unistd.h>
#endif
#include "lbfgs.h"

const double ETA=1e-7;
const double XTOL=1e-7;
const int MP = 6;
const int LP = 6;
const double GTOL = .9;
const double STPMIN = 1e-20;
const double STPMAX = 1e20;
const double FTOL=1e-4;
const int MAXFEV = 20;
const double P5=.5;
const double P66=.66;
const double XTRAPF=4;


inline double inner_product(int n,double *a, double *b){
	int i;
	double s=0;
	for(i=0;i<n;i++)
		s+=a[i]*b[i];
	return s;
}
inline void daxpy(int n,double *a, double b,double *c){//a+=b*c
	int i;
	for(i=0;i<n;i++)
		a[i]+=b*c[i];
}
LBFGS::~LBFGS()
{
	char fn[1000];
	int i;
	for(i=0;i<m;i++)
	{
		sprintf(fn,"__s_%d",i);
		if(!access(fn,0))//file exists
			unlink(fn);
		sprintf(fn,"__y_%d",i);
		if(!access(fn,0))//file exists
			unlink(fn);
	}
	if(!access("__x",0))
		unlink("__x");
	if(!access("__q",0))
		unlink("__q");
	if(!access("__w0",0))
		unlink("__w0");
	if(!access("__w1",0))
		unlink("__w1");

}
void LBFGS::save(double *buf, char *file_name, int index)
{
	char fn[1000];
	if(index>=0)
		sprintf(fn,"%s_%d",file_name,index);
	else
		strcpy(fn,file_name);
	FILE *fp=fopen(fn,"wb");
	fwrite(buf,sizeof(double),n,fp);

	fclose(fp);
}
void LBFGS::load(double *buf, char *file_name, int index)
{
	char fn[1000];
	if(index>=0)
		sprintf(fn,"%s_%d",file_name,index);
	else
		strcpy(fn,file_name);
	FILE *fp=fopen(fn,"rb");
	fread(buf,sizeof(double),n,fp);
	fclose(fp);
}
bool LBFGS::init(int parameter_num, int depth, int pri){
	int i;
	prior=pri;
	n = parameter_num;
	m = depth;
	s.resize(m);
	y.resize(m);
	if(prior==0)
	{
		w.resize(n + 2*m*n);
		if(!w.size())
			return false;
		for(i=0;i<m;i++)
		{
			s[i]=&w[n+i*n];
			y[i]=&w[n+m*n+i*n];
		}
		q=&w[0];
		diag.resize(n);
		if(!diag.size())
			return false;
	}
	alpha.resize(m);
	rho.resize(m);

	if(!s.size() || !y.size() || !alpha.size() || !rho.size())
		return false;
	
	return true;
}




/* Subroutine */ static int mcstep_(double *stx, double *fx, double *dx, double *sty, double *fy, double *dy, double *stp, 
									double *fp, double *dp, int *brackt,
     double *stpmin, double *stpmax, int *info)
{
  /* System generated locals */
  double d__1, d__2, d__3;



  /* Local variables */
  static double sgnd, stpc, stpf, stpq, p, q, gamma, r__, s, theta;
  static int bound;

  *info = 0;

  /*     CHECK THE INPUT PARAMETERS FOR ERRORS. */

  if (*brackt && ((*stp <= min(*stx,*sty) || *stp >= max(*stx,*sty)) || *dx *
                  (*stp - *stx) >=(float)0. || *stpmax < *stpmin)) {
    return 0;
  }

  /*     DETERMINE IF THE DERIVATIVES HAVE OPPOSITE SIGN. */

  sgnd = *dp *(*dx / fabs(*dx));

  /*     FIRST CASE. A HIGHER FUNCTION VALUE. */
  /*     THE MINIMUM IS BRACKETED. IF THE CUBIC STEP IS CLOSER */
  /*     TO STX THAN THE QUADRATIC STEP, THE CUBIC STEP IS TAKEN, */
  /*     ELSE THE AVERAGE OF THE CUBIC AND QUADRATIC STEPS IS TAKEN. */

  if (*fp > *fx) {
    *info = 1;
    bound = 1;
    theta =(*fx - *fp) * 3 / (*stp - *stx) + *dx + *dp;
    /* Computing MAX */
    d__1 = fabs(theta), d__2 = fabs(*dx), d__1 = max(d__1,d__2), d__2 = fabs(
                                                                          *dp);
    s = max(d__1,d__2);
    /* Computing 2nd power */
    d__1 = theta / s;
    gamma = s * sqrt(d__1 * d__1 - *dx / s *(*dp / s));
    if (*stp < *stx) {
      gamma = -gamma;
    }
    p = gamma - *dx + theta;
    q = gamma - *dx + gamma + *dp;
    r__ = p / q;
    stpc = *stx + r__ *(*stp - *stx);
    stpq = *stx + *dx /((*fx - *fp) / (*stp - *stx) + *dx) / 2 * (*stp -
                                                                  *stx);
    if ((d__1 = stpc - *stx, fabs(d__1)) < (d__2 = stpq - *stx, fabs(d__2)))
      {
        stpf = stpc;
      } else {
        stpf = stpc + (stpq - stpc) / 2;
      }
    *brackt = 1;

    /*     SECOND CASE. A LOWER FUNCTION VALUE AND DERIVATIVES OF */
    /*     OPPOSITE SIGN. THE MINIMUM IS BRACKETED. IF THE CUBIC */
    /*     STEP IS CLOSER TO STX THAN THE QUADRATIC(SECANT) STEP, */
    /*     THE CUBIC STEP IS TAKEN, ELSE THE QUADRATIC STEP IS TAKEN. */

  } else if (sgnd < (float)0.) {
    *info = 2;
    bound = 0;
    theta =(*fx - *fp) * 3 / (*stp - *stx) + *dx + *dp;
    /* Computing MAX */
    d__1 = fabs(theta), d__2 = fabs(*dx), d__1 = max(d__1,d__2), d__2 = fabs(
                                                                          *dp);
    s = max(d__1,d__2);
    /* Computing 2nd power */
    d__1 = theta / s;
    gamma = s * sqrt(d__1 * d__1 - *dx / s *(*dp / s));
    if (*stp > *stx) {
      gamma = -gamma;
    }
    p = gamma - *dp + theta;
    q = gamma - *dp + gamma + *dx;
    r__ = p / q;
    stpc = *stp + r__ *(*stx - *stp);
    stpq = *stp + *dp /(*dp - *dx) * (*stx - *stp);
    if ((d__1 = stpc - *stp, fabs(d__1)) > (d__2 = stpq - *stp, fabs(d__2)))
      {
        stpf = stpc;
      } else {
        stpf = stpq;
      }
    *brackt = 1;

    /*     THIRD CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE */
    /*     SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DECREASES. */
    /*     THE CUBIC STEP IS ONLY USED IF THE CUBIC TENDS TO INFINITY */
    /*     IN THE DIRECTION OF THE STEP OR IF THE MINIMUM OF THE CUBIC */
    /*     IS BEYOND STP. OTHERWISE THE CUBIC STEP IS DEFINED TO BE */
    /*     EITHER STPMIN OR STPMAX. THE QUADRATIC(SECANT) STEP IS ALSO */
    /*     COMPUTED AND IF THE MINIMUM IS BRACKETED THEN THE THE STEP */
    /*     CLOSEST TO STX IS TAKEN, ELSE THE STEP FARTHEST AWAY IS TAKEN. */

  } else if (fabs(*dp) < fabs(*dx)) {
    *info = 3;
    bound = 1;
    theta =(*fx - *fp) * 3 / (*stp - *stx) + *dx + *dp;
    /* Computing MAX */
    d__1 = fabs(theta), d__2 = fabs(*dx), d__1 = max(d__1,d__2), d__2 = fabs(
                                                                          *dp);
    s = max(d__1,d__2);

    /*        THE CASE GAMMA = 0 ONLY ARISES IF THE CUBIC DOES NOT TEND */
    /*        TO INFINITY IN THE DIRECTION OF THE STEP. */

    /* Computing MAX */
    /* Computing 2nd power */
    d__3 = theta / s;
    d__1 = 0., d__2 = d__3 * d__3 - *dx / s *(*dp / s);
    gamma = s * sqrt((max(d__1,d__2)));
    if (*stp > *stx) {
      gamma = -gamma;
    }
    p = gamma - *dp + theta;
    q = gamma + (*dx - *dp) + gamma;
    r__ = p / q;
    if (r__ < (float)0. && gamma != (float)0.) {
      stpc = *stp + r__ *(*stx - *stp);
    } else if (*stp > *stx) {
      stpc = *stpmax;
    } else {
      stpc = *stpmin;
    }
    stpq = *stp + *dp /(*dp - *dx) * (*stx - *stp);
    if (*brackt) {
      if ((d__1 = *stp - stpc, fabs(d__1)) < (d__2 = *stp - stpq, fabs(
                                                                     d__2))) {
        stpf = stpc;
      } else {
        stpf = stpq;
      }
    } else {
      if ((d__1 = *stp - stpc, fabs(d__1)) > (d__2 = *stp - stpq, fabs(
                                                                     d__2))) {
        stpf = stpc;
      } else {
        stpf = stpq;
      }
    }

    /*     FOURTH CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE */
    /*     SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DOES */
    /*     NOT DECREASE. IF THE MINIMUM IS NOT BRACKETED, THE STEP */
    /*     IS EITHER STPMIN OR STPMAX, ELSE THE CUBIC STEP IS TAKEN. */

  } else {
    *info = 4;
    bound = 0;
    if (*brackt) {
      theta =(*fp - *fy) * 3 / (*sty - *stp) + *dy + *dp;
      /* Computing MAX */
      d__1 = fabs(theta), d__2 = fabs(*dy), d__1 = max(d__1,d__2), d__2 =
        fabs(*dp);
      s = max(d__1,d__2);
      /* Computing 2nd power */
      d__1 = theta / s;
      gamma = s * sqrt(d__1 * d__1 - *dy / s *(*dp / s));
      if (*stp > *sty) {
        gamma = -gamma;
      }
      p = gamma - *dp + theta;
      q = gamma - *dp + gamma + *dy;
      r__ = p / q;
      stpc = *stp + r__ *(*sty - *stp);
      stpf = stpc;
    } else if (*stp > *stx) {
      stpf = *stpmax;
    } else {
      stpf = *stpmin;
    }
  }

  /*     UPDATE THE INTERVAL OF UNCERTAINTY. THIS UPDATE DOES NOT */
  /*     DEPEND ON THE NEW STEP OR THE CASE ANALYSIS ABOVE. */

  if (*fp > *fx) {
    *sty = *stp;
    *fy = *fp;
    *dy = *dp;
  } else {
    if (sgnd < (float)0.) {
      *sty = *stx;
      *fy = *fx;
      *dy = *dx;
    }
    *stx = *stp;
    *fx = *fp;
    *dx = *dp;
  }

  /*     COMPUTE THE NEW STEP AND SAFEGUARD IT. */

  stpf = min(*stpmax,stpf);
  stpf = max(*stpmin,stpf);
  *stp = stpf;
  if (*brackt && bound) {
    if (*sty > *stx) {
      /* Computing MIN */
      d__1 = *stx + (*sty - *stx) * (float).66;
      *stp = min(d__1,*stp);
    } else {
      /* Computing MAX */
      d__1 = *stx + (*sty - *stx) * (float).66;
      *stp = max(d__1,*stp);
    }
  }
  return 0;

  /*     LAST LINE OF SUBROUTINE MCSTEP. */

} /* mcstep_ */



static int mcsrch_(int *n, double *x, double *f, double *g, double *s, double* stp, int *info, int *nfev, double *wa, int orthant)
{

  /* System generated locals */
  int i__1;
  double d__1;

  /* Local variables */
  static double dgxm, dgym;
  static int j, infoc;
  static double finit, width, stmin, stmax;
  static int stage1;
  static double width1, ftest1, dg, fm, fx, fy;
  static int brackt;
  static double dginit, dgtest;
  static double dgm, dgx, dgy, fxm, fym, stx, sty;

  /* Parameter adjustments */
  --wa;
  --s;
  --g;
  --x;


  /* Function Body */
  if (*info == -1) {
    goto L45;
  }