text:
/* multimin/vector_bfgs.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Fabrice Rossi
 *
 * 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.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* vector_bfgs.c -- Limited memory Broyden-Fletcher-Goldfarb-Shanno
   conjugate gradient method */

/* Modified by Brian Gough to use single iteration structure */

#include <config.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_blas.h>

#include "directional_minimize.c"

typedef struct
{
  int iter;
  double step;
  double max_step;
  double tol;
  gsl_vector *x1;
  gsl_vector *dx1;
  gsl_vector *x2;
  double g0norm;
  double pnorm;
  gsl_vector *p;
  gsl_vector *x0;
  gsl_vector *g0;
  gsl_vector *dx0;
  gsl_vector *dg0;
}
vector_bfgs_state_t;

static int
vector_bfgs_alloc (void *vstate, size_t n)
{
  vector_bfgs_state_t *state = (vector_bfgs_state_t *) vstate;

  state->x1 = gsl_vector_calloc (n);

  if (state->x1 == 0)
    {
      GSL_ERROR ("failed to allocate space for x1", GSL_ENOMEM);
    }

  state->dx1 = gsl_vector_calloc (n);

  if (state->dx1 == 0)
    {
      gsl_vector_free (state->x1);
      GSL_ERROR ("failed to allocate space for dx1", GSL_ENOMEM);
    }

  state->x2 = gsl_vector_calloc (n);

  if (state->x2 == 0)
    {
      gsl_vector_free (state->dx1);
      gsl_vector_free (state->x1);
      GSL_ERROR ("failed to allocate space for x2", GSL_ENOMEM);
    }

  state->p = gsl_vector_calloc (n);

  if (state->p == 0)
    {
      gsl_vector_free (state->x2);
      gsl_vector_free (state->dx1);
      gsl_vector_free (state->x1);
      GSL_ERROR ("failed to allocate space for p", GSL_ENOMEM);
    }

  state->x0 = gsl_vector_calloc (n);

  if (state->x0 == 0)
    {
      gsl_vector_free (state->p);
      gsl_vector_free (state->x2);
      gsl_vector_free (state->dx1);
      gsl_vector_free (state->x1);
      GSL_ERROR ("failed to allocate space for g0", GSL_ENOMEM);
    }

  state->g0 = gsl_vector_calloc (n);

  if (state->g0 == 0)
    {
      gsl_vector_free (state->x0);
      gsl_vector_free (state->p);
      gsl_vector_free (state->x2);
      gsl_vector_free (state->dx1);
      gsl_vector_free (state->x1);
      GSL_ERROR ("failed to allocate space for g0", GSL_ENOMEM);
    }

  state->dx0 = gsl_vector_calloc (n);

  if (state->dx0 == 0)
    {
      gsl_vector_free (state->g0);
      gsl_vector_free (state->x0);
      gsl_vector_free (state->p);
      gsl_vector_free (state->x2);
      gsl_vector_free (state->dx1);
      gsl_vector_free (state->x1);
      GSL_ERROR ("failed to allocate space for g0", GSL_ENOMEM);
    }

  state->dg0 = gsl_vector_calloc (n);

  if (state->dg0 == 0)
    {
      gsl_vector_free (state->dx0);
      gsl_vector_free (state->g0);
      gsl_vector_free (state->x0);
      gsl_vector_free (state->p);
      gsl_vector_free (state->x2);
      gsl_vector_free (state->dx1);
      gsl_vector_free (state->x1);
      GSL_ERROR ("failed to allocate space for g0", GSL_ENOMEM);
    }

  return GSL_SUCCESS;
}

static int
vector_bfgs_set (void *vstate, gsl_multimin_function_fdf * fdf,
		 const gsl_vector * x, double *f, gsl_vector * gradient,
		 double step_size, double tol)
{
  vector_bfgs_state_t *state = (vector_bfgs_state_t *) vstate;

  state->iter = 0;
  state->step = step_size;
  state->max_step = step_size;
  state->tol = tol;

  GSL_MULTIMIN_FN_EVAL_F_DF (fdf, x, f, gradient);

  /* Use the gradient as the initial direction */

  gsl_vector_memcpy (state->x0, x);
  gsl_vector_memcpy (state->p, gradient);
  gsl_vector_memcpy (state->g0, gradient);

  {
    double gnorm = gsl_blas_dnrm2 (gradient);
    state->pnorm = gnorm;
    state->g0norm = gnorm;
  }

  return GSL_SUCCESS;
}

static void
vector_bfgs_free (void *vstate)
{
  vector_bfgs_state_t *state = (vector_bfgs_state_t *) vstate;

  gsl_vector_free (state->dg0);
  gsl_vector_free (state->dx0);
  gsl_vector_free (state->g0);
  gsl_vector_free (state->x0);
  gsl_vector_free (state->p);
  gsl_vector_free (state->x2);
  gsl_vector_free (state->dx1);
  gsl_vector_free (state->x1);
}

static int
vector_bfgs_restart (void *vstate)
{
  vector_bfgs_state_t *state = (vector_bfgs_state_t *) vstate;

  state->iter = 0;
  return GSL_SUCCESS;
}

static int
vector_bfgs_iterate (void *vstate, gsl_multimin_function_fdf * fdf,
		     gsl_vector * x, double *f,
		     gsl_vector * gradient, gsl_vector * dx)
{
  vector_bfgs_state_t *state = (vector_bfgs_state_t *) vstate;

  gsl_vector *x1 = state->x1;
  gsl_vector *dx1 = state->dx1;
  gsl_vector *x2 = state->x2;
  gsl_vector *p = state->p;
  gsl_vector *g0 = state->g0;
  gsl_vector *x0 = state->x0;

  double pnorm = state->pnorm;
  double g0norm = state->g0norm;

  double fa = *f, fb, fc;
  double dir;
  double stepa = 0.0, stepb, stepc = state->step, tol = state->tol;

  double g1norm;
  double pg;

  if (pnorm == 0.0 || g0norm == 0.0)
    {
      gsl_vector_set_zero (dx);
      return GSL_ENOPROG;
    }

  /* Determine which direction is downhill, +p or -p */

  gsl_blas_ddot (p, gradient, &pg);

  dir = (pg >= 0.0) ? +1.0 : -1.0;

  /* Compute new trial point at x_c= x - step * p, where p is the
     current direction */

  take_step (x, p, stepc, dir / pnorm, x1, dx);

  /* Evaluate function and gradient at new point xc */

  fc = GSL_MULTIMIN_FN_EVAL_F (fdf, x1);

  if (fc < fa)
    {
      /* Success, reduced the function value */
      state->step = stepc * 2.0;
      *f = fc;
      gsl_vector_memcpy (x, x1);
      GSL_MULTIMIN_FN_EVAL_DF (fdf, x1, gradient);
      state->g0norm = gsl_blas_dnrm2 (gradient);
      return GSL_SUCCESS;
    }

#ifdef DEBUG
  printf ("got stepc = %g fc = %g\n", stepc, fc);
#endif

  /* Do a line minimisation in the region (xa,fa) (xc,fc) to find an
     intermediate (xb,fb) satisifying fa > fb < fc.  Choose an initial
     xb based on parabolic interpolation */

  intermediate_point (fdf, x, p, dir / pnorm, pg,
		      stepa, stepc, fa, fc, x1, dx1, gradient, &stepb, &fb);

  minimize (fdf, x, p, dir / pnorm,
	    stepa, stepb, stepc, fa, fb, fc, tol,
	    x1, dx1, x2, dx, gradient, &(state->step), f, &g1norm);

  gsl_vector_memcpy (x, x2);

  /* Choose a new conjugate direction for the next step */

  state->iter = (state->iter + 1) % x->size;

  if (state->iter == 0)
    {
      gsl_vector_memcpy (p, gradient);
      state->pnorm = g1norm;
    }
  else
    {
      /* This is the BFGS update: */
      /* p' = g1 - A dx - B dg */
      /* A = - (1+ dg.dg/dx.dg) B + dg.g/dx.dg */
      /* B = dx.g/dx.dg */

      gsl_vector *dx0 = state->dx0;
      gsl_vector *dg0 = state->dg0;

      double dxg, dgg, dxdg, dgnorm, A, B;

      /* dx0 = x - x0 */
      gsl_vector_memcpy (dx0, x);
      gsl_blas_daxpy (-1.0, x0, dx0);

      /* dg0 = g - g0 */
      gsl_vector_memcpy (dg0, gradient);
      gsl_blas_daxpy (-1.0, g0, dg0);

      gsl_blas_ddot (dx0, gradient, &dxg);
      gsl_blas_ddot (dg0, gradient, &dgg);
      gsl_blas_ddot (dx0, dg0, &dxdg);

      dgnorm = gsl_blas_dnrm2 (dg0);

      B = dxg / dxdg;
      A = -(1.0 + dgnorm * dgnorm / dxdg) * B + dgg / dxdg;

      gsl_vector_memcpy (p, gradient);
      gsl_blas_daxpy (-A, dx0, p);
      gsl_blas_daxpy (-B, dg0, p);

      state->pnorm = gsl_blas_dnrm2 (p);
    }

  gsl_vector_memcpy (g0, gradient);
  gsl_vector_memcpy (x0, x);
  state->g0norm = gsl_blas_dnrm2 (g0);

#ifdef DEBUG
  printf ("updated conjugate directions\n");
  printf ("p: ");
  gsl_vector_fprintf (stdout, p, "%g");
  printf ("g: ");
  gsl_vector_fprintf (stdout, gradient, "%g");
#endif

  return GSL_SUCCESS;
}

static const gsl_multimin_fdfminimizer_type vector_bfgs_type = {
  "vector_bfgs",		/* name */
  sizeof (vector_bfgs_state_t),
  &vector_bfgs_alloc,
  &vector_bfgs_set,
  &vector_bfgs_iterate,
  &vector_bfgs_restart,
  &vector_bfgs_free
};

const gsl_multimin_fdfminimizer_type
  * gsl_multimin_fdfminimizer_vector_bfgs = &vector_bfgs_type;


text:
/* randist/discrete.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007, 2009 James Theiler, Brian Gough
 *
 * 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 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

/*
   Random Discrete Events

   Given K discrete events with different probabilities P[k]
   produce a value k consistent with its probability.

   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 3 of the
   License, or (at your option) any later version.

   This program is distributed in the hope that it will be useful, but
   WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   General Public License for more details.  You should have received
   a copy of the GNU General Public License along with this program;
   if not, write to the Free Foundation, Inc., 59 Temple Place, Suite
   330, Boston, MA 02111-1307 USA
*/

/*
 * Based on: Alastair J Walker, An efficient method for generating
 * discrete random variables with general distributions, ACM Trans
 * Math Soft 3, 253-256 (1977).  See also: D. E. Knuth, The Art of
 * Computer Programming, Volume 2 (Seminumerical algorithms), 3rd
 * edition, Addison-Wesley (1997), p120.

 * Walker's algorithm does some preprocessing, and provides two
 * arrays: floating point F[k] and integer A[k].  A value k is chosen
 * from 0..K-1 with equal likelihood, and then a uniform random number
 * u is compared to F[k].  If it is less than F[k], then k is
 * returned.  Otherwise, A[k] is returned.

 * Walker's original paper describes an O(K^2) algorithm for setting
 * up the F and A arrays.  I found this disturbing since I wanted to
 * use very large values of K.  I'm sure I'm not the first to realize
 * this, but in fact the preprocessing can be done in O(K) steps.

 * A figure of merit for the preprocessing is the average value for
 * the F[k]'s (that is, SUM_k F[k]/K); this corresponds to the
 * probability that k is returned, instead of A[k], thereby saving a
 * redirection.  Walker's O(K^2) preprocessing will generally improve
 * that figure of merit, compared to my cheaper O(K) method; from some
 * experiments with a perl script, I get values of around 0.6 for my
 * method and just under 0.75 for Walker's.  Knuth has pointed out
 * that finding _the_ optimum lookup tables, which maximize the
 * average F[k], is a combinatorially difficult problem.  But any
 * valid preprocessing will still provide O(1) time for the call to
 * gsl_ran_discrete().  I find that if I artificially set F[k]=1 --
 * ie, better than optimum! -- I get a speedup of maybe 20%, so that's
 * the maximum I could expect from the most expensive preprocessing.
 * Folding in the difference of 0.6 vs 0.75, I'd estimate that the
 * speedup would be less than 10%.

 * I've not implemented it here, but one compromise is to sort the
 * probabilities once, and then work from the two ends inward.  This
 * requires O(K log K), still lots cheaper than O(K^2), and from my
 * experiments with the perl script, the figure of merit is within
 * about 0.01 for K up to 1000, and no sign of diverging (in fact,
 * they seemed to be converging, but it's hard to say with just a
 * handful of runs).

 * The O(K) algorithm goes through all the p_k's and decides if they
 * are "smalls" or "bigs" according to whether they are less than or
 * greater than the mean value 1/K.  The indices to the smalls and the
 * bigs are put in separate stacks, and then we work through the
 * stacks together.  For each small, we pair it up with the next big
 * in the stack (Walker always wanted to pair up the smallest small
 * with the biggest big).  The small "borrows" from the big just
 * enough to bring the small up to mean.  This reduces the size of the
 * big, so the (smaller) big is compared again to the mean, and if it
 * is smaller, it gets "popped" from the big stack and "pushed" to the
 * small stack.  Otherwise, it stays put.  Since every time we pop a
 * small, we are able to deal with it right then and there, and we
 * never have to pop more than K smalls, then the algorithm is O(K).

 * This implementation sets up two separate stacks, and allocates K
 * elements between them.  Since neither stack ever grows, we do an
 * extra O(K) pass through the data to determine how many smalls and
 * bigs there are to begin with and allocate appropriately.  In all
 * there are 2*K*sizeof(double) transient bytes of memory that are
 * used than returned, and K*(sizeof(int)+sizeof(double)) bytes used
 * in the lookup table.

 * Walker spoke of using two random numbers (an integer 0..K-1, and a
 * floating point u in [0,1]), but Knuth points out that one can just
 * use the integer and fractional parts of K*u where u is in [0,1].
 * In fact, Knuth further notes that taking F'[k]=(k+F[k])/K, one can
 * directly compare u to F'[k] without having to explicitly set
 * u=K*u-int(K*u).

 * Usage:

 * Starting with an array of probabilities P, initialize and do
 * preprocessing with a call to:

 *    gsl_rng *r;
 *    gsl_ran_discrete_t *f;
 *    f = gsl_ran_discrete_preproc(K,P);

 * Then, whenever a random index 0..K-1 is desired, use

 *    k = gsl_ran_discrete(r,f);

 * Note that several different randevent struct's can be
 * simultaneously active.

 * Aside: A very clever alternative approach is described in
 * Abramowitz and Stegun, p 950, citing: Marsaglia, Random variables
 * and computers, Proc Third Prague Conference in Probability Theory,
 * 1962.  A more accesible reference is: G. Marsaglia, Generating
 * discrete random numbers in a computer, Comm ACM 6, 37-38 (1963).
 * If anybody is interested, I (jt) have also coded up this version as
 * part of another software package.  However, I've done some
 * comparisons, and the Walker method is both faster and more stingy
 * with memory.  So, in the end I decided not to include it with the
 * GSL package.

 * Written 26 Jan 1999, James Theiler, jt@lanl.gov
 * Adapted to GSL, 30 Jan 1999, jt

 */

#include <config.h>
#include <stdio.h>              /* used for NULL, also fprintf(stderr,...) */
#include <stdlib.h>             /* used for malloc's */
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#define DEBUG 0
#define KNUTH_CONVENTION 1      /* Saves a few steps of arithmetic
                                 * in the call to gsl_ran_discrete()
                                 */

/*** Begin Stack (this code is used just in this file) ***/

/* Stack code converted to use unsigned indices (i.e. s->i == 0 now
   means an empty stack, instead of -1), for consistency and to give a
   bigger allowable range. BJG */

typedef struct {
    size_t N;                      /* max number of elts on stack */
    size_t *v;                     /* array of values on the stack */
    size_t i;                      /* index of top of stack */
} gsl_stack_t;

static gsl_stack_t *
new_stack(size_t N) {
    gsl_stack_t *s;
    s = (gsl_stack_t *)malloc(sizeof(gsl_stack_t));
    s->N = N;
    s->i = 0;                  /* indicates stack is empty */
    s->v = (size_t *)malloc(sizeof(size_t)*N);
    return s;
}

static int
push_stack(gsl_stack_t *s, size_t v)
{
    if ((s->i) >= (s->N)) {
      return -1; /* stack overflow (shouldn't happen) */
    }
    (s->v)[s->i] = v;
    s->i += 1;
    return 0;
}

static size_t pop_stack(gsl_stack_t *s)
{
    if ((s->i) == 0) {
      GSL_ERROR ("internal error - stack exhausted", GSL_ESANITY);
    }
    s->i -= 1;
    return ((s->v)[s->i]);
}

static inline size_t size_stack(const gsl_stack_t *s)
{
    return s->i;
}

static void free_stack(gsl_stack_t *s)
{
    free((char *)(s->v));
    free((char *)s);
}

/*** End Stack ***/

/*** Begin Walker's Algorithm ***/

gsl_ran_discrete_t *
gsl_ran_discrete_preproc(size_t Kevents, const double *ProbArray)
{
    size_t k,b,s;
    gsl_ran_discrete_t *g;
    size_t nBigs, nSmalls;
    gsl_stack_t *Bigs;
    gsl_stack_t *Smalls;
    double *E;
    double pTotal = 0.0, mean, d;

    if (Kevents < 1) {
      /* Could probably treat Kevents=1 as a special case */

      GSL_ERROR_VAL ("number of events must be a positive integer",
                        GSL_EINVAL, 0);
    }

    /* Make sure elements of ProbArray[] are positive.
     * Won't enforce that sum is unity; instead will just normalize
     */

    for (k=0; k<Kevents; ++k) {
        if (ProbArray[k] < 0) {
          GSL_ERROR_VAL ("probabilities must be non-negative",
                            GSL_EINVAL, 0) ;
        }
        pTotal += ProbArray[k];
    }

    /* Begin setting up the main "object" (just a struct, no steroids) */
    g = (gsl_ran_discrete_t *)malloc(sizeof(gsl_ran_discrete_t));
    g->K = Kevents;
    g->F = (double *)malloc(sizeof(double)*Kevents);
    g->A = (size_t *)malloc(sizeof(size_t)*Kevents);

    E = (double *)malloc(sizeof(double)*Kevents);

    if (E==NULL) {
      GSL_ERROR_VAL ("Cannot allocate memory for randevent", GSL_ENOMEM, 0);
    }

    for (k=0; k<Kevents; ++k) {
        E[k] = ProbArray[k]/pTotal;
    }

    /* Now create the Bigs and the Smalls */
    mean = 1.0/Kevents;
    nSmalls=nBigs=0;
    {
      /* Temporarily use which[k] = g->A[k] to indicate small or large */
      size_t * const which = g->A;

      for (k=0; k<Kevents; ++k) {
        if (E[k] < mean) {
          ++nSmalls; which[k] = 0;
        } else {
          ++nBigs; which[k] = 1;
        }
      }

      Bigs   = new_stack(nBigs);
      Smalls = new_stack(nSmalls);
      for (k=0; k<Kevents; ++k) {
        gsl_stack_t * Dest = which[k] ? Bigs : Smalls;
        int status = push_stack(Dest,k);
        if (status)
          GSL_ERROR_VAL ("failed to build stacks", GSL_EFAILED, 0);
      }
    }

    /* Now work through the smalls */
    while (size_stack(Smalls) > 0) {
        s = pop_stack(Smalls);
        if (size_stack(Bigs) == 0) {
            (g->A)[s]=s;
            (g->F)[s]=1.0;
            continue;
        }
        b = pop_stack(Bigs);
        (g->A)[s]=b;
        (g->F)[s]=Kevents*E[s];
#if DEBUG
        fprintf(stderr,"s=%2d, A=%2d, F=%.4f\n",s,(g->A)[s],(g->F)[s]);
#endif
        d = mean - E[s];
        E[s] += d;              /* now E[s] == mean */
        E[b] -= d;
        if (E[b] < mean) {
            push_stack(Smalls,b); /* no longer big, join ranks of the small */
        }
        else if (E[b] > mean) {
            push_stack(Bigs,b); /* still big, put it back where you found it */
        }
        else {
            /* E[b]==mean implies it is finished too */
            (g->A)[b]=b;
            (g->F)[b]=1.0;
        }
    }
    while (size_stack(Bigs) > 0) {
        b = pop_stack(Bigs);
        (g->A)[b]=b;
        (g->F)[b]=1.0;
    }
    /* Stacks have been emptied, and A and F have been filled */

    if ( size_stack(Smalls) != 0) {
      GSL_ERROR_VAL ("Smalls stack has not been emptied",
                     GSL_ESANITY, 0 );
    }

#if 0
    /* if 1, then artificially set all F[k]'s to unity.  This will
     * give wrong answers, but you'll get them faster.  But, not
     * that much faster (I get maybe 20%); that's an upper bound
     * on what the optimal preprocessing would give.
     */
    for (k=0; k<Kevents; ++k) {
        (g->F)[k] = 1.0;
    }
#endif

#if KNUTH_CONVENTION
    /* For convenience, set F'[k]=(k+F[k])/K */
    /* This saves some arithmetic in gsl_ran_discrete(); I find that
     * it doesn't actually make much difference.
     */
    for (k=0; k<Kevents; ++k) {
        (g->F)[k] += k;
        (g->F)[k] /= Kevents;
    }
#endif

    free_stack(Bigs);
    free_stack(Smalls);
    free((char *)E);

    return g;
}

size_t
gsl_ran_discrete(const gsl_rng *r, const gsl_ran_discrete_t *g)
{
    size_t c=0;
    double u,f;
    u = gsl_rng_uniform(r);
#if KNUTH_CONVENTION
    c = (u*(g->K));
#else
    u *= g->K;
    c = u;
    u -= c;
#endif
    f = (g->F)[c];
    /* fprintf(stderr,"c,f,u: %d %.4f %f\n",c,f,u); */
    if (f == 1.0) return c;

    if (u < f) {
        return c;
    }
    else {
        return (g->A)[c];
    }
}

void gsl_ran_discrete_free(gsl_ran_discrete_t *g)
{
    RETURN_IF_NULL (g);
    free((char *)(g->A));
    free((char *)(g->F));
    free((char *)g);
}

double
gsl_ran_discrete_pdf(size_t k, const gsl_ran_discrete_t *g)
{
    size_t i,K;
    double f,p=0;
    K= g->K;
    if (k>K) return 0;
    for (i=0; i<K; ++i) {
        f = (g->F)[i];
#if KNUTH_CONVENTION
        f = K*f-i;
#endif
        if (i==k) {
            p += f;
        } else if (k == (g->A)[i]) {
            p += 1.0 - f;
        }
    }
    return p/K;
}

ext: c