This is the mail archive of the gsl-discuss@sourceware.org mailing list for the GSL project.
Index Nav: | [Date Index] [Subject Index] [Author Index] [Thread Index] | |
---|---|---|
Message Nav: | [Date Prev] [Date Next] | [Thread Prev] [Thread Next] |
Other format: | [Raw text] |
Am Freitag, den 17.03.2006, 16:13 +0000 schrieb Brian Gough: > Ivo Alxneit-Kamber writes: > > i fear you missunderstood (maybe my comments were not really sufficient) > > or i do not really get your point here. > > > > Just that you don't need to modify the API if your algorithm only > needs an initial scalar step instead of a vector -- use the magnitude > of the vector (and ignore the direction if necessary, or put it to > some use). That fits the existing API. > ok, done so (plus your previous remarks) - supplied own (vax) random number generator. - API kept (old initializer still supported). - updated test case. both initializer are tested now. - updated documentation (brian, could you please have a look at the language). -- Dr. Ivo Alxneit Laboratory for Solar Technology phone: +41 56 310 4092 Paul Scherrer Institute fax: +41 56 310 2688 CH-5232 Villigen http://solar.web.psi.ch Switzerland gnupg key: 0x515E30C7
--- /tmp/gsl-1.7/multimin/simplex.c 2005-06-26 15:25:35.000000000 +0200 +++ gsl-1.7/multimin/simplex.c 2006-03-20 08:51:58.000000000 +0100 @@ -32,7 +32,7 @@ */ #include <config.h> -#include <stdlib.h> +//#include <stdlib.h> #include <gsl/gsl_blas.h> #include <gsl/gsl_multimin.h> @@ -45,6 +45,287 @@ } nmsimplex_state_t; + +/* The following is a copy of gsl_rng_type vax. */ + +typedef struct + { + unsigned long int x; + } +my_rng_state_t; + +static void +my_rng_set(my_rng_state_t *s, const unsigned long int seed) +{ + s->x = seed; +} + +static double +my_rng_get(my_rng_state_t *s) +{ + s->x = (69069 * s->x + 1) & 0xffffffffUL; + return (s->x / 4294967296.0); +} + + +static int +old_init(void *vstate, gsl_multimin_function * f, + const gsl_vector * x, const gsl_vector * step_size) +{ + int status; + size_t i; + double val; + gsl_vector *xtemp; + + nmsimplex_state_t *state = (nmsimplex_state_t *) vstate; + + xtemp = state->ws1; + + /* first point is the original x0 */ + + val = GSL_MULTIMIN_FN_EVAL (f, x); + gsl_matrix_set_row (state->x1, 0, x); + gsl_vector_set (state->y1, 0, val); + + /* following points are initialized to x0 + step_size */ + + for (i = 0; i < x->size; i++) + { + status = gsl_vector_memcpy (xtemp, x); + + if (status != 0) + { + GSL_ERROR ("vector memcopy failed", GSL_EFAILED); + } + + val = gsl_vector_get (xtemp, i) + gsl_vector_get (step_size, i); + gsl_vector_set (xtemp, i, val); + val = GSL_MULTIMIN_FN_EVAL (f, xtemp); + gsl_matrix_set_row (state->x1, i + 1, xtemp); + gsl_vector_set (state->y1, i + 1, val); + } + return GSL_SUCCESS; +} + + +static double +a_m(const int m, const double c) +{ + return (c * sqrt(1 - c) / sqrt((1 + (m - 1) * c) * (1 + m * c))); +} + +static double +b_m(const size_t m, const double c) +{ + return (sqrt((1 - c) * (1 + m * c) / (1 + (m - 1) * c))); +} + +static gsl_matrix * +init_npod(const size_t n) +{ + /* Given the following n+1 vectors e_i: 0<i<=n in n dimensions + + e_0(a0, 0, 0, ...) + e_1(a0, b1, 0, ..) + e_2(a0, a1, b2, .) + e_n-1(a0, ......, bn) + e_n(a0, ....., -bn) + + that form a regular normalized n+1_pod in n dimensions, i.e. + a regular triangle in 2D or a regulart tetrahedron in 3D. + + Requesting that + all (e_i, e_i) = 1 + all (e_i, e_j) = c */ + + size_t i, j; + const double c = -1.0 / n; + gsl_vector_view v; + gsl_matrix *e; + + + e = gsl_matrix_calloc(n + 1, n); + if (e == NULL) + { + return (NULL); + } + + /* Initialize e_0. */ + + gsl_matrix_set(e, 0, 0, 1.0); + + /* Initialize e_1 to e_n-1. */ + + for (i = 1; i < n; i++) + { + for (j = 0; j < i; j++) + gsl_matrix_set(e, i, j, a_m((int) j, c)); + + gsl_matrix_set(e, i, j, b_m(j, c)); + + for (j = i + 1; j < n; j++) + gsl_matrix_set(e, i, j, 0.0); + } + + /* Initialize e_n. */ + + v = gsl_matrix_row(e, n - 1); + gsl_matrix_set_row(e, n, &v.vector); + + gsl_matrix_set(e, n, n - 1, -gsl_matrix_get(e, n - 1, n - 1)); + + return (e); + +} + +static gsl_matrix * +init_basis(const size_t n) +{ + /* orthonolmal basis, randomly oriented */ + + size_t i, j; + char *val; + unsigned int seed; + my_rng_state_t *r; + gsl_matrix *a, *basis; + + + a = gsl_matrix_alloc(n, n); + if (a == NULL) + { + return (NULL); + } + + basis = gsl_matrix_alloc(n, n); + if (basis == NULL) + { + gsl_matrix_free(a); + return (NULL); + } + + r = (my_rng_state_t *)malloc(sizeof(my_rng_state_t)); + if (r == NULL) + { + gsl_matrix_free(basis); + gsl_matrix_free(a); + return (NULL); + } + + val = getenv("GSL_NM_SIMPLEX_SEED"); + + if (val) + seed = strtoul(val, 0, 0); + else + seed = 0; + + my_rng_set(r, seed); + + + /* Get n random vectors of length n. */ + + for (i = 0; i < n; i++) + for (j = 0; j < n; j++) + { + double value = my_rng_get(r); + gsl_matrix_set(a, i, j, 2.0 * value - 1.0); + } + + free(r); + + /* Orthogonalize matrix a, result in matrix basis. + FIXME: we do not check for colinear vectors. */ + + gsl_matrix_memcpy(basis, a); + + for (i = 0; i < n; i++) + { + gsl_vector_view v1 = gsl_matrix_row(a, i); + gsl_vector_view v2 = gsl_matrix_row(basis, i); + + for (j = 0; j < i; j++) + { + double t1, t2; + gsl_vector_view v3 = gsl_matrix_row(basis, j); + + gsl_blas_ddot(&v1.vector, &v3.vector, &t1); + gsl_blas_ddot(&v3.vector, &v3.vector, &t2); + + gsl_blas_daxpy(-t1 / t2, &v3.vector, &v2.vector); + + } + } + + gsl_matrix_free(a); + + /* Normalize basis vectors. */ + + for (i = 0; i < n; i++) + { + gsl_vector_view v1 = gsl_matrix_row(basis, i); + double t1 = gsl_blas_dnrm2(&v1.vector); + + gsl_vector_scale(&v1.vector, 1 / t1); + + } + + return (basis); + +} +static int +new_init(void *vstate, gsl_multimin_function * f, + const gsl_vector * x, const double step_size) +{ + size_t i; + gsl_matrix *n_pod, *basis; + + nmsimplex_state_t *state = (nmsimplex_state_t *) vstate; + + size_t n = x->size; + + n_pod = init_npod(n); + if (n_pod == NULL) + { + GSL_ERROR("failed to allocate space for n_pod", GSL_ENOMEM); + } + + /* Construct a randomly oriented, orthonormal basis in n dimensions. */ + + basis = init_basis(n); + if (basis == NULL) + { + gsl_matrix_free(n_pod); + GSL_ERROR("failed to allocate space for basis", GSL_ENOMEM); + } + + /* Scale n_pod by step_size and express it in basis. result in state->x1 */ + + gsl_blas_dgemm(CblasNoTrans, CblasTrans, step_size, n_pod, basis, 0.0, + state->x1); + + gsl_matrix_free(basis); + gsl_matrix_free(n_pod); + + /* Now state->x1 contains n+1 vectors forming a n+1 dimensional, randomly + oriented simplex of size gsl_vector_get(step_size, 0) centered at the + origin. we now offset it by the initial vector x and initialize the + corresponding y values in state->y1. */ + + for (i = 0; i <= n; i++) + { + double val; + gsl_vector_view v = gsl_matrix_row(state->x1, i); + + gsl_vector_add(&v.vector, x); + + val = GSL_MULTIMIN_FN_EVAL(f, &v.vector); + gsl_vector_set(state->y1, i, val); + } + + return GSL_SUCCESS; + +} + + static double nmsimplex_move_corner (const double coeff, const nmsimplex_state_t * state, size_t corner, gsl_vector * xc, @@ -221,35 +502,53 @@ double *size, const gsl_vector * step_size) { int status; - size_t i; - double val; + size_t i, n_non_zero; nmsimplex_state_t *state = (nmsimplex_state_t *) vstate; - gsl_vector *xtemp = state->ws1; - - /* first point is the original x0 */ - - val = GSL_MULTIMIN_FN_EVAL (f, x); - gsl_matrix_set_row (state->x1, 0, x); - gsl_vector_set (state->y1, 0, val); - - /* following points are initialized to x0 + step_size */ + size_t n = x->size; + + /* Find out which method to use for setting the initial simplex. + If only one element of step_size is non-zero we set the initial + simplex to a randomly oriented n-pod. */ + + n_non_zero = 0; + for(i=0; i<n; i++) + { + if(gsl_vector_get(step_size, i) != 0.0) + { + n_non_zero++; + } + } - for (i = 0; i < x->size; i++) + if(n_non_zero == n) { - status = gsl_vector_memcpy (xtemp, x); + /* Use old methode. */ + + status = old_init(vstate, f, x, step_size); + if (status != GSL_SUCCESS) + { + return status; + } - if (status != 0) + } + else if(n_non_zero == 1) + { + /* New initialization method. + Prepare regular, normalized n_pod: + equilateral triangle 2D + equilateral tetrahedron 3D. */ + + status = new_init(vstate, f, x, gsl_vector_get(step_size, 0)); + if (status != GSL_SUCCESS) { - GSL_ERROR ("vector memcopy failed", GSL_EFAILED); - } + return status; + } - val = gsl_vector_get (xtemp, i) + gsl_vector_get (step_size, i); - gsl_vector_set (xtemp, i, val); - val = GSL_MULTIMIN_FN_EVAL (f, xtemp); - gsl_matrix_set_row (state->x1, i + 1, xtemp); - gsl_vector_set (state->y1, i + 1, val); + } + else + { + GSL_ERROR("invalid step size", GSL_EINVAL); } /* Initialize simplex size */
--- /tmp/gsl-1.7/multimin/test.c 2005-06-26 15:25:35.000000000 +0200 +++ gsl-1.7/multimin/test.c 2006-03-18 19:20:24.000000000 +0100 @@ -33,7 +33,8 @@ initpt_function initpt, const gsl_multimin_fdfminimizer_type *T); int -test_f(const char * desc, gsl_multimin_function *f, initpt_function initpt); +test_f(const char * desc, gsl_multimin_function *f, initpt_function initpt, + const int method); int main (void) @@ -59,9 +60,15 @@ T++; } - test_f("Roth", &roth_fmin, roth_initpt); - test_f("Wood", &wood_fmin, wood_initpt); - test_f("Rosenbrock", &rosenbrock_fmin, rosenbrock_initpt); +/* test old initialization method */ + test_f("Roth", &roth_fmin, roth_initpt, 0); + test_f("Wood", &wood_fmin, wood_initpt, 0); + test_f("Rosenbrock", &rosenbrock_fmin, rosenbrock_initpt, 0); + +/* test new initialization method */ + test_f("Roth", &roth_fmin, roth_initpt, 1); + test_f("Wood", &wood_fmin, wood_initpt, 1); + test_f("Rosenbrock", &rosenbrock_fmin, rosenbrock_initpt, 1); T = fdfminimizers; @@ -133,7 +140,8 @@ } int -test_f(const char * desc, gsl_multimin_function *f, initpt_function initpt) +test_f(const char * desc, gsl_multimin_function *f, initpt_function initpt, + const int method) { /* currently this function tests only nmsimplex */ @@ -143,15 +151,23 @@ gsl_vector *x = gsl_vector_alloc (f->n); - gsl_vector *step_size = gsl_vector_alloc (f->n); + gsl_vector *step_size = gsl_vector_calloc (f->n); gsl_multimin_fminimizer *s; (*initpt) (x); - for (i = 0; i < f->n; i++) - gsl_vector_set (step_size, i, 1); - + if (method == 0) + { + for (i = 0; i < f->n; i++) + gsl_vector_set (step_size, i, 1); + } + + if (method == 1) + { + gsl_vector_set(step_size, 0, 1); + } + s = gsl_multimin_fminimizer_alloc(T, f->n); gsl_multimin_fminimizer_set (s, f, x, step_size);
--- /tmp/gsl-1.7/doc/multimin.texi 2005-05-21 15:28:05.000000000 +0200 +++ gsl-1.7/doc/multimin.texi 2006-03-20 09:19:27.000000000 +0100 @@ -457,7 +457,18 @@ @end ifinfo @noindent These vectors form the @math{n+1} vertices of a simplex in @math{n} -dimensions. On each iteration the algorithm tries to improve +dimensions. Alternatively, if only the first component of @var{step_size} +is non zero, a randomly oriented regular @var{n+1}-pod of the size given +in @var{step_size}, centered at @var{x} is used as simplex to start +with. The orientation of this initial simplex can be controlled by +setting the environment variable @code{GSL_NM_SIMPLEX_SEED} to some +integer value. The initial simplex then corresponds to +an equilateral triangle (two dimensional problem) or to +a regular tetrahedron (three dimensional problem). Note, that +the parameters to be determined should be are of similar magnitude +if this initialize is applied. + +On each iteration the algorithm tries to improve the parameter vector @math{p_i} corresponding to the highest function value by simple geometrical transformations. These are reflection, reflection followed by expansion, contraction and multiple
Attachment:
signature.asc
Description: Dies ist ein digital signierter Nachrichtenteil
Index Nav: | [Date Index] [Subject Index] [Author Index] [Thread Index] | |
---|---|---|
Message Nav: | [Date Prev] [Date Next] | [Thread Prev] [Thread Next] |