htswanalysis/MACS/lib/gsl/gsl-1.11/multimin/linear_wrapper.c

   1 typedef struct
   2 {
   3   gsl_function_fdf fdf_linear;
   4   gsl_multimin_function_fdf *fdf;
   5   /* fixed values */
   6   const gsl_vector *x;
   7   const gsl_vector *g;
   8   const gsl_vector *p;
   9
  10   /* cached values, for x(alpha) = x + alpha * p */
  11   double f_alpha;
  12   double df_alpha;
  13   gsl_vector *x_alpha;
  14   gsl_vector *g_alpha;
  15
  16   /* cache "keys" */
  17   double f_cache_key;
  18   double df_cache_key;
  19   double x_cache_key;
  20   double g_cache_key;
  21 }
  22 wrapper_t;
  23
  24 static void
  25 moveto (double alpha, wrapper_t * w)
  26 {
  27   if (alpha == w->x_cache_key)  /* using previously cached position */
  28     {
  29       return;
  30     }
  31
  32   /* set x_alpha = x + alpha * p */
  33
  34   gsl_vector_memcpy (w->x_alpha, w->x);
  35   gsl_blas_daxpy (alpha, w->p, w->x_alpha);
  36
  37   w->x_cache_key = alpha;
  38 }
  39
  40 static double
  41 slope (wrapper_t * w)           /* compute gradient . direction */
  42 {
  43   double df;
  44   gsl_blas_ddot (w->g_alpha, w->p, &df);
  45   return df;
  46 }
  47
  48 static double
  49 wrap_f (double alpha, void *params)
  50 {
  51   wrapper_t *w = (wrapper_t *) params;
  52   if (alpha == w->f_cache_key)  /* using previously cached f(alpha) */
  53     {
  54       return w->f_alpha;
  55     }
  56
  57   moveto (alpha, w);
  58
  59   w->f_alpha = GSL_MULTIMIN_FN_EVAL_F (w->fdf, w->x_alpha);
  60   w->f_cache_key = alpha;
  61
  62   return w->f_alpha;
  63 }
  64
  65 static double
  66 wrap_df (double alpha, void *params)
  67 {
  68   wrapper_t *w = (wrapper_t *) params;
  69   if (alpha == w->df_cache_key) /* using previously cached df(alpha) */
  70     {
  71       return w->df_alpha;
  72     }
  73
  74   moveto (alpha, w);
  75
  76   if (alpha != w->g_cache_key)
  77     {
  78       GSL_MULTIMIN_FN_EVAL_DF (w->fdf, w->x_alpha, w->g_alpha);
  79       w->g_cache_key = alpha;
  80     }
  81
  82   w->df_alpha = slope (w);
  83   w->df_cache_key = alpha;
  84
  85   return w->df_alpha;
  86 }
  87
  88 static void
  89 wrap_fdf (double alpha, void *params, double *f, double *df)
  90 {
  91   wrapper_t *w = (wrapper_t *) params;
  92
  93   /* Check for previously cached values */
  94
  95   if (alpha == w->f_cache_key && alpha == w->df_cache_key)
  96     {
  97       *f = w->f_alpha;
  98       *df = w->df_alpha;
  99       return;
 100     }
 101
 102   if (alpha == w->f_cache_key || alpha == w->df_cache_key)
 103     {
 104       *f = wrap_f (alpha, params);
 105       *df = wrap_df (alpha, params);
 106       return;
 107     }
 108
 109   moveto (alpha, w);
 110   GSL_MULTIMIN_FN_EVAL_F_DF (w->fdf, w->x_alpha, &w->f_alpha, w->g_alpha);
 111   w->f_cache_key = alpha;
 112   w->g_cache_key = alpha;
 113
 114   w->df_alpha = slope (w);
 115   w->df_cache_key = alpha;
 116
 117   *f = w->f_alpha;
 118   *df = w->df_alpha;
 119 }
 120
 121 static void
 122 prepare_wrapper (wrapper_t * w, gsl_multimin_function_fdf * fdf,
 123                  const gsl_vector * x, double f, const gsl_vector *g,
 124                  const gsl_vector * p,
 125                  gsl_vector * x_alpha, gsl_vector *g_alpha)
 126 {
 127   w->fdf_linear.f = &wrap_f;
 128   w->fdf_linear.df = &wrap_df;
 129   w->fdf_linear.fdf = &wrap_fdf;
 130   w->fdf_linear.params = (void *)w;  /* pointer to "self" */
 131
 132   w->fdf = fdf;
 133
 134   w->x = x;
 135   w->g = g;
 136   w->p = p;
 137
 138   w->x_alpha = x_alpha;
 139   w->g_alpha = g_alpha;
 140
 141   gsl_vector_memcpy(w->x_alpha, w->x);
 142   w->x_cache_key = 0.0;
 143
 144   w->f_alpha = f;
 145   w->f_cache_key = 0.0;
 146
 147   gsl_vector_memcpy(w->g_alpha, w->g);
 148   w->g_cache_key = 0.0;
 149
 150   w->df_alpha = slope(w);
 151   w->df_cache_key = 0.0;
 152 }
 153
 154 static void
 155 update_position (wrapper_t * w, double alpha, gsl_vector *x, double *f, gsl_vector *g)
 156 {
 157   /* ensure that everything is fully cached */
 158   { double f_alpha, df_alpha; wrap_fdf (alpha, w, &f_alpha, &df_alpha); } ;
 159
 160   *f = w->f_alpha;
 161   gsl_vector_memcpy(x, w->x_alpha);
 162   gsl_vector_memcpy(g, w->g_alpha);
 163 }
 164
 165 static void
 166 change_direction (wrapper_t * w)
 167 {
 168   /* Convert the cache values from the end of the current minimisation
 169      to those needed for the start of the next minimisation, alpha=0 */
 170
 171   /* The new x_alpha for alpha=0 is the current position */
 172   gsl_vector_memcpy (w->x_alpha, w->x);
 173   w->x_cache_key = 0.0;
 174
 175   /* The function value does not change */
 176   w->f_cache_key = 0.0;
 177
 178   /* The new g_alpha for alpha=0 is the current gradient at the endpoint */
 179   gsl_vector_memcpy (w->g_alpha, w->g);
 180   w->g_cache_key = 0.0;
 181
 182   /* Calculate the slope along the new direction vector, p */
 183   w->df_alpha = slope (w);
 184   w->df_cache_key = 0.0;
 185 }