123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186 |
- typedef struct
- {
- gsl_function_fdf fdf_linear;
- gsl_multimin_function_fdf *fdf;
- /* fixed values */
- const gsl_vector *x;
- const gsl_vector *g;
- const gsl_vector *p;
- /* cached values, for x(alpha) = x + alpha * p */
- double f_alpha;
- double df_alpha;
- gsl_vector *x_alpha;
- gsl_vector *g_alpha;
- /* cache "keys" */
- double f_cache_key;
- double df_cache_key;
- double x_cache_key;
- double g_cache_key;
- }
- wrapper_t;
- static void
- moveto (double alpha, wrapper_t * w)
- {
- if (alpha == w->x_cache_key) /* using previously cached position */
- {
- return;
- }
- /* set x_alpha = x + alpha * p */
- gsl_vector_memcpy (w->x_alpha, w->x);
- gsl_blas_daxpy (alpha, w->p, w->x_alpha);
- w->x_cache_key = alpha;
- }
- static double
- slope (wrapper_t * w) /* compute gradient . direction */
- {
- double df;
- gsl_blas_ddot (w->g_alpha, w->p, &df);
- return df;
- }
- static double
- wrap_f (double alpha, void *params)
- {
- wrapper_t *w = (wrapper_t *) params;
- if (alpha == w->f_cache_key) /* using previously cached f(alpha) */
- {
- return w->f_alpha;
- }
- moveto (alpha, w);
- w->f_alpha = GSL_MULTIMIN_FN_EVAL_F (w->fdf, w->x_alpha);
- w->f_cache_key = alpha;
- return w->f_alpha;
- }
- static double
- wrap_df (double alpha, void *params)
- {
- wrapper_t *w = (wrapper_t *) params;
- if (alpha == w->df_cache_key) /* using previously cached df(alpha) */
- {
- return w->df_alpha;
- }
- moveto (alpha, w);
- if (alpha != w->g_cache_key)
- {
- GSL_MULTIMIN_FN_EVAL_DF (w->fdf, w->x_alpha, w->g_alpha);
- w->g_cache_key = alpha;
- }
- w->df_alpha = slope (w);
- w->df_cache_key = alpha;
- return w->df_alpha;
- }
- static void
- wrap_fdf (double alpha, void *params, double *f, double *df)
- {
- wrapper_t *w = (wrapper_t *) params;
- /* Check for previously cached values */
- if (alpha == w->f_cache_key && alpha == w->df_cache_key)
- {
- *f = w->f_alpha;
- *df = w->df_alpha;
- return;
- }
- if (alpha == w->f_cache_key || alpha == w->df_cache_key)
- {
- *f = wrap_f (alpha, params);
- *df = wrap_df (alpha, params);
- return;
- }
- moveto (alpha, w);
- GSL_MULTIMIN_FN_EVAL_F_DF (w->fdf, w->x_alpha, &w->f_alpha, w->g_alpha);
- w->f_cache_key = alpha;
- w->g_cache_key = alpha;
- w->df_alpha = slope (w);
- w->df_cache_key = alpha;
- *f = w->f_alpha;
- *df = w->df_alpha;
- }
- static void
- prepare_wrapper (wrapper_t * w, gsl_multimin_function_fdf * fdf,
- const gsl_vector * x, double f, const gsl_vector *g,
- const gsl_vector * p,
- gsl_vector * x_alpha, gsl_vector *g_alpha)
- {
- w->fdf_linear.f = &wrap_f;
- w->fdf_linear.df = &wrap_df;
- w->fdf_linear.fdf = &wrap_fdf;
- w->fdf_linear.params = (void *)w; /* pointer to "self" */
- w->fdf = fdf;
- w->x = x;
- w->g = g;
- w->p = p;
- w->x_alpha = x_alpha;
- w->g_alpha = g_alpha;
- gsl_vector_memcpy(w->x_alpha, w->x);
- w->x_cache_key = 0.0;
-
- w->f_alpha = f;
- w->f_cache_key = 0.0;
-
- gsl_vector_memcpy(w->g_alpha, w->g);
- w->g_cache_key = 0.0;
- w->df_alpha = slope(w);
- w->df_cache_key = 0.0;
- }
- static void
- update_position (wrapper_t * w, double alpha, gsl_vector *x, double *f, gsl_vector *g)
- {
- /* ensure that everything is fully cached */
- { double f_alpha, df_alpha; wrap_fdf (alpha, w, &f_alpha, &df_alpha); } ;
- *f = w->f_alpha;
- gsl_vector_memcpy(x, w->x_alpha);
- gsl_vector_memcpy(g, w->g_alpha);
- }
- static void
- change_direction (wrapper_t * w)
- {
- /* Convert the cache values from the end of the current minimisation
- to those needed for the start of the next minimisation, alpha=0 */
- /* The new x_alpha for alpha=0 is the current position */
- gsl_vector_memcpy (w->x_alpha, w->x);
- w->x_cache_key = 0.0;
- /* The function value does not change */
- w->f_cache_key = 0.0;
- /* The new g_alpha for alpha=0 is the current gradient at the endpoint */
- gsl_vector_memcpy (w->g_alpha, w->g);
- w->g_cache_key = 0.0;
- /* Calculate the slope along the new direction vector, p */
- w->df_alpha = slope (w);
- w->df_cache_key = 0.0;
- }
|