This is the mail archive of the
gsl-discuss@sources.redhat.com
mailing list for the GSL project.
Re: simplex minimization
- From: Tim F <fenn at agora dot dhs dot org>
- To: Brian Gough <brian dot gough at network-theory dot co dot uk>
- Cc: gsl-discuss at sources dot redhat dot com
- Date: Wed, 1 Jan 2003 21:40:35 -0500 (EST)
- Subject: Re: simplex minimization
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, ¶m.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);
}