simplex minimization

Tim F fenn@agora.dhs.org
Thu Jan 2 02:40:00 GMT 2003


On Wed, 1 Jan 2003, Brian Gough wrote:

>
> Yes please send an example.  Thanks.
> 

I hope this example isn't too terse... let me know if other information 
will help.

	Tim F

#include <stdio.h>
#include <stdlib.h>

#include <glib.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_sf_exp.h>
#include <gsl/gsl_vector.h>

// num iterations for solver (max)
#define NUM_ITER 100

struct func_data{
  gint low, high;
  gdouble *gstar;
  reflection_p data;
};

/*
 * perform the following fit of k and [M]:
 * Fo = k*Fc*e^-([T]^T[M][T])
 * (T = 2pi[hkl])
 */
void calc_scale_factor(void){
  gint i, j, status;
  // number of parameters
  const size_t p = 7;
  gdouble param_init[7];
  gsl_vector_view param;
  struct func_data d;

  // f only minimizer
  gdouble ssval;
  const gsl_multimin_fminimizer_type *T;
  gsl_multimin_fminimizer *s;
  gsl_multimin_function f;
  gsl_vector *ss;
 
  // our initial M is just 0
  // order: u11 u22 u33 u12 u13 u23
  for (i=0; i<6; i++)
    param_init[i] = 0.0;

  // initial k (<FoFc> / <Fc^2>)
  param_init[6] = initial_scale_factor(curdata->sfo_hkls);

  param = gsl_vector_view_array(param_init, p);

  // set up data for solver
  /*
   * gstar is a 3x3 symmetric tensor
   * sfo_hkls is the Fo and Fc data
   */
  d.low = curdata->low;
  d.high = curdata->high;
  d.gstar = curdata->gstar;
  d.data = curdata->sfo_hkls;

  // set up function
  f.f = &scale_f;
  f.n = p;
  f.params = (void *)&d;

  // set up solver
  T = gsl_multimin_fminimizer_nmsimplex;
  s = gsl_multimin_fminimizer_alloc(T, p);
  ss = gsl_vector_alloc(p);
  gsl_multimin_fminimizer_set(s, &f, &param.vector, ss);

  for (i=1; i<NUM_ITER; i++){
    status = gsl_multimin_fminimizer_iterate(s);
    if (status){
      g_print("GSL minimizer: %s\n", gsl_strerror(status));
      break;
    }

    status = gsl_multimin_test_size(gsl_multimin_fminimizer_size(s), 1e-2);
    ssval = gsl_multimin_fminimizer_size(s);
    if (status == GSL_SUCCESS)
      g_print("minimum found at:\n");

    g_print("iteration: %d\nparams:\n", i);
    for (j=0; j<p; j++)
      g_print("%9.5f\n", gsl_vector_get(s->x, j));
    g_print(" f(x): %f tot. size: %.12f\n\n", s->fval, ssval);

    if (status != GSL_CONTINUE)
      break;
  }

  // set values based on minimization
  ...
  // free up minimizer
  gsl_multimin_fminimizer_free(s);
  gsl_vector_free(ss);
}

// function
gdouble scale_f(const gsl_vector *x, void *params){
  gint i;
  gint low = ((struct func_data *)params)->low;
  gint high = ((struct func_data *)params)->high;
  gdouble *gstar = ((struct func_data *)params)->gstar;
  reflection_p data = ((struct func_data *)params)->data;
  gdouble m[] = {gsl_vector_get(x, 0), gsl_vector_get(x, 1), 
                 gsl_vector_get(x, 2), gsl_vector_get(x, 3), 
                 gsl_vector_get(x, 4), gsl_vector_get(x, 5)};
  gdouble scale_k = gsl_vector_get(x, 6);

  gdouble u;
  guint64 sum;

  for (i=low, sum=0.0; i<=high; i++){
    if (!data[i].mate)
      continue;

    /*
     * in some cases, the off diagonal elements of gstar (elements 3-5)
     * may be zero, forcing the m[3] to m[5] parameters to zero
     * during the minimization
     */
    // calculate [T]^T[M][T]
    u = gsl_pow_2(data[i].hkl[0]) * gstar[0] * m[0] +
      gsl_pow_2(data[i].hkl[1]) * gstar[1] * m[1] +
      gsl_pow_2(data[i].hkl[2]) * gstar[2] * m[2] +
      2.0 * data[i].hkl[0] * data[i].hkl[1] * gstar[3] * m[3] +
      2.0 * data[i].hkl[0] * data[i].hkl[2] * gstar[4] * m[4] +
      2.0 * data[i].hkl[1] * data[i].hkl[2] * gstar[5] * m[5];

    // calculate sum: (Fo-k*Fc*e^(-[T]^T[M][T]))^2 / (sigFo)^2
    sum += (guint64) ((gsl_pow_2(fabs(data[i].data.fo[0]) - 
                       fabs(gsl_sf_exp_mult(-2.0*M_PI*M_PI*u, 
                       scale_k * gsl_complex_abs(data[i].mate->data.fcphi))))) / 
                       gsl_pow_2(data[i].data.fo[1]));
  }

  return((gdouble) sum);
}




More information about the Gsl-discuss mailing list