This is the mail archive of the gsl-discuss@sources.redhat.com 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]

Preliminary Patch for Genetic Algorithm in GSL

Hi everybody,

I join a patch I use against the cvs trunk of gsl 03242005.

It implements differential evolution (a kind of genetic algorithm very efficient, used in mathematica 5) to globally minimize objective function knowing only the function values not the derivatives.

It work very well for difficult function with a lot of local minimum.

Please try it and let me know. I am available to any comment.

apply patch inside gsl directory with:
patch -Np1 -i path-to-patch

should be no problem because it is all inside a new dir called gamin. configure.ac and Makefile.am are also modified to build the new lib.

in the directory gamin there is a test.c which serves its name and as an example. the header file is gsl_gamin.h .

Happy Hacking.

Guillaume MICHEL

-------------------------------------------------------------------------------------------------------- 




diff -Nur ./Makefile.am ../gsl-gamin/Makefile.am
--- ./Makefile.am    2004-12-24 14:55:10.000000000 +0100
+++ ../gsl-gamin/Makefile.am    2005-03-24 16:59:17.700383600 +0100
@@ -2,9 +2,9 @@

# AUTOMAKE_OPTIONS = readme-alpha

-SUBDIRS = gsl utils sys test err const complex cheb block vector matrix permutation combination sort ieee-utils cblas blas linalg eigen specfunc dht qrng rng randist fft poly fit multifit statistics siman sum integration interpolation histogram ode-initval roots multiroots min multimin monte ntuple diff deriv cdf wavelet doc
+SUBDIRS = gsl utils sys test err const complex cheb block vector matrix permutation combination sort ieee-utils cblas blas linalg eigen specfunc dht qrng rng randist fft poly fit multifit statistics siman sum integration interpolation histogram ode-initval roots multiroots min multimin monte ntuple diff deriv cdf wavelet gamin doc

-SUBLIBS = block/libgslblock.la blas/libgslblas.la complex/libgslcomplex.la cheb/libgslcheb.la dht/libgsldht.la diff/libgsldiff.la deriv/libgslderiv.la eigen/libgsleigen.la err/libgslerr.la fft/libgslfft.la fit/libgslfit.la histogram/libgslhistogram.la ieee-utils/libgslieeeutils.la integration/libgslintegration.la interpolation/libgslinterpolation.la linalg/libgsllinalg.la matrix/libgslmatrix.la min/libgslmin.la monte/libgslmonte.la multifit/libgslmultifit.la multimin/libgslmultimin.la multiroots/libgslmultiroots.la ntuple/libgslntuple.la ode-initval/libgslodeiv.la permutation/libgslpermutation.la combination/libgslcombination.la poly/libgslpoly.la qrng/libgslqrng.la randist/libgslrandist.la rng/libgslrng.la roots/libgslroots.la siman/libgslsiman.la sort/libgslsort.la specfunc/libgslspecfunc.la statistics/libgslstatistics.la sum/libgslsum.la sys/libgslsys.la test/libgsltest.la utils/libutils.la vector/libgslvector.la cdf/libgslcdf.la wavelet/libgslwavelet.la
+SUBLIBS = block/libgslblock.la blas/libgslblas.la complex/libgslcomplex.la cheb/libgslcheb.la dht/libgsldht.la diff/libgsldiff.la deriv/libgslderiv.la eigen/libgsleigen.la err/libgslerr.la fft/libgslfft.la fit/libgslfit.la histogram/libgslhistogram.la ieee-utils/libgslieeeutils.la integration/libgslintegration.la interpolation/libgslinterpolation.la linalg/libgsllinalg.la matrix/libgslmatrix.la min/libgslmin.la monte/libgslmonte.la multifit/libgslmultifit.la multimin/libgslmultimin.la multiroots/libgslmultiroots.la ntuple/libgslntuple.la ode-initval/libgslodeiv.la permutation/libgslpermutation.la combination/libgslcombination.la poly/libgslpoly.la qrng/libgslqrng.la randist/libgslrandist.la rng/libgslrng.la roots/libgslroots.la siman/libgslsiman.la sort/libgslsort.la specfunc/libgslspecfunc.la statistics/libgslstatistics.la sum/libgslsum.la sys/libgslsys.la test/libgsltest.la utils/libutils.la vector/libgslvector.la cdf/libgslcdf.la wavelet/libgslwavelet.la gamin/libgslgamin.la

pkginclude_HEADERS = gsl_math.h gsl_pow_int.h gsl_nan.h gsl_machine.h gsl_mode.h gsl_precision.h gsl_types.h gsl_version.h

@@ -40,3 +40,4 @@
#dummy_LDADD = $(SUBLIBS)
#main_SOURCES = version.c env.c
#main_LDADD = libgsl.la
+
diff -Nur ./configure.ac ../gsl-gamin/configure.ac
--- ./configure.ac    2005-03-05 12:52:52.000000000 +0100
+++ ../gsl-gamin/configure.ac    2005-03-24 16:59:17.284446832 +0100
@@ -356,5 +356,5 @@
fi

dnl
-AC_CONFIG_FILES([gsl-config gsl.pc gsl_version.h gsl.spec gsl/Makefile test/Makefile err/Makefile sys/Makefile utils/Makefile const/Makefile min/Makefile multimin/Makefile ieee-utils/Makefile fft/Makefile specfunc/Makefile dht/Makefile fit/Makefile multifit/Makefile statistics/Makefile sum/Makefile roots/Makefile multiroots/Makefile ntuple/Makefile poly/Makefile qrng/Makefile rng/Makefile randist/Makefile siman/Makefile integration/Makefile interpolation/Makefile doc/Makefile block/Makefile vector/Makefile matrix/Makefile histogram/Makefile monte/Makefile ode-initval/Makefile cblas/Makefile blas/Makefile linalg/Makefile eigen/Makefile permutation/Makefile combination/Makefile sort/Makefile complex/Makefile diff/Makefile deriv/Makefile cheb/Makefile cdf/Makefile wavelet/Makefile Makefile])
+AC_CONFIG_FILES([gsl-config gsl.pc gsl_version.h gsl.spec gsl/Makefile test/Makefile err/Makefile sys/Makefile utils/Makefile const/Makefile min/Makefile multimin/Makefile ieee-utils/Makefile fft/Makefile specfunc/Makefile dht/Makefile fit/Makefile multifit/Makefile statistics/Makefile sum/Makefile roots/Makefile multiroots/Makefile ntuple/Makefile poly/Makefile qrng/Makefile rng/Makefile randist/Makefile siman/Makefile integration/Makefile interpolation/Makefile doc/Makefile block/Makefile vector/Makefile matrix/Makefile histogram/Makefile monte/Makefile ode-initval/Makefile cblas/Makefile blas/Makefile linalg/Makefile eigen/Makefile permutation/Makefile combination/Makefile sort/Makefile complex/Makefile diff/Makefile deriv/Makefile cheb/Makefile cdf/Makefile wavelet/Makefile gamin/Makefile Makefile])
AC_OUTPUT
diff -Nur ./gamin/Makefile.am ../gsl-gamin/gamin/Makefile.am
--- ./gamin/Makefile.am 1970-01-01 01:00:00.000000000 +0100
+++ ../gsl-gamin/gamin/Makefile.am 2005-03-24 16:59:32.363154520 +0100
@@ -0,0 +1,19 @@
+noinst_LTLIBRARIES = libgslgamin.la
+
+pkginclude_HEADERS = gsl_gamin.h
+
+INCLUDES= -I$(top_builddir)
+
+libgslgamin_la_SOURCES = minimizer.c best1exp.c best2exp.c rand1exp.c \
+rand2exp.c randtobest1exp.c terminator_convergence.c \
+terminator_generation.c terminator_pop_conv.c
+
+check_PROGRAMS = test #demo
+
+TESTS = $(check_PROGRAMS)
+
+test_SOURCES = test.c
+test_LDADD = libgslgamin.la ../vector/libgslvector.la \
+../block/libgslblock.la ../ieee-utils/libgslieeeutils.la \
+../err/libgslerr.la ../test/libgsltest.la ../sys/libgslsys.la \
+../utils/libutils.la ../rng/libgslrng.la ../randist/libgslrandist.la
diff -Nur ./gamin/best1exp.c ../gsl-gamin/gamin/best1exp.c
--- ./gamin/best1exp.c 1970-01-01 01:00:00.000000000 +0100
+++ ../gsl-gamin/gamin/best1exp.c 2005-03-24 16:59:32.363154520 +0100
@@ -0,0 +1,65 @@
+/* #include <config.h> */
+
+#include <stdlib.h>
+#include <string.h>
+
+#include <stdio.h>
+
+#include "gsl_gamin.h"
+
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_randist.h>
+
+
+static void
+best1exp_get_candidate(const size_t i,
+ const size_t n,
+ const size_t NP,
+ const gsl_vector ** parents,
+ const double F,
+ const double CR,
+ const gsl_rng * r,
+ const size_t best,
+ gsl_vector * candidate)
+{
+ size_t r1, r2, k, L;
+ + do /* Pick a random population member */
+ { /* Endless loop for NP < 2 !!! */
+ r1 = gsl_rng_uniform_int(r, NP);
+ }
+ while(r1==i); + + do /* Pick a random population member */
+ { /* Endless loop for NP < 3 !!! */
+ r2 = gsl_rng_uniform_int(r, NP);
+ }
+ while((r2==i) || (r2==r1));
+ + gsl_vector_memcpy(candidate, parents[i]);
+
+ k = gsl_rng_uniform_int(r, n);
+ L = 0;
+ do
+ {
+ gsl_vector_set(candidate, k,
+ gsl_vector_get(parents[best],k) +
+ F*(gsl_vector_get(parents[r1],k)-
+ gsl_vector_get(parents[r2],k)));
+ k = (k+1) % n;
+ L++;
+ }
+ while( (gsl_rng_uniform(r)< CR) && (L < n) );
+}
+
+
+
+static const gsl_gamin_fminimizer_type best1exp_type =
+ { "DE/best/1/exp", /* name */
+ &best1exp_get_candidate
+ };
+
+const gsl_gamin_fminimizer_type *
+gsl_gamin_fminimizer_best1exp = &best1exp_type;
+
diff -Nur ./gamin/best2exp.c ../gsl-gamin/gamin/best2exp.c
--- ./gamin/best2exp.c 1970-01-01 01:00:00.000000000 +0100
+++ ../gsl-gamin/gamin/best2exp.c 2005-03-24 16:59:32.363154520 +0100
@@ -0,0 +1,79 @@
+/* #include <config.h> */
+
+#include <stdlib.h>
+#include <string.h>
+
+#include <stdio.h>
+
+#include "gsl_gamin.h"
+
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_randist.h>
+
+
+static void
+best2exp_get_candidate(const size_t i,
+ const size_t n,
+ const size_t NP,
+ const gsl_vector ** parents,
+ const double F,
+ const double CR,
+ const gsl_rng * r,
+ const size_t best,
+ gsl_vector * candidate)
+{
+ size_t r1, r2, r3, r4, k, L;
+ + do /* Pick a random population member */
+ { /* Endless loop for NP < 2 !!! */
+ r1 = gsl_rng_uniform_int(r, NP);
+ }
+ while(r1==i); + + do /* Pick a random population member */
+ { /* Endless loop for NP < 3 !!! */
+ r2 = gsl_rng_uniform_int(r, NP);
+ }
+ while((r2==i) || (r2==r1));
+ + do /* Pick a random population member */
+ { /* Endless loop for NP < 4 !!! */
+ r3 = gsl_rng_uniform_int(r, NP);
+ }
+ while((r3==i) || (r3==r1)|| (r3==r2) );
+ + do /* Pick a random population member */
+ { /* Endless loop for NP < 5 !!! */
+ r4 = gsl_rng_uniform_int(r, NP);
+ }
+ while((r4==i) || (r4==r1) || (r4==r2) || (r4==r3));
+
+ gsl_vector_memcpy(candidate, parents[i]);
+
+ k = gsl_rng_uniform_int(r, n);
+ L = 0;
+ do
+ {
+ gsl_vector_set(candidate, k,
+ gsl_vector_get(parents[best],k) +
+ F*(gsl_vector_get(parents[r1],k)+
+ gsl_vector_get(parents[r2],k)-
+ gsl_vector_get(parents[r3],k)-
+ gsl_vector_get(parents[r4],k)));
+ k = (k+1) % n;
+ L++;
+ }
+ while( (gsl_rng_uniform(r)< CR) && (L < n) );
+}
+
+
+
+static const gsl_gamin_fminimizer_type best2exp_type =
+ { "DE/best/2/exp", /* name */
+ &best2exp_get_candidate
+ };
+
+const gsl_gamin_fminimizer_type *
+gsl_gamin_fminimizer_best2exp = &best2exp_type;
+
diff -Nur ./gamin/gsl_gamin.h ../gsl-gamin/gamin/gsl_gamin.h
--- ./gamin/gsl_gamin.h 1970-01-01 01:00:00.000000000 +0100
+++ ../gsl-gamin/gamin/gsl_gamin.h 2005-03-24 16:59:32.363154520 +0100
@@ -0,0 +1,283 @@
+#ifndef __GSL_GAMIN_H__
+#define __GSL_GAMIN_H__
+
+/**
+ \file
+ \brief main include file for GAMIN.
+*/
+
+#include <stdlib.h>
+
+#include <gsl/gsl_types.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_rng.h>
+
+#undef __BEGIN_DECLS
+#undef __END_DECLS
+#ifdef __cplusplus
+# define __BEGIN_DECLS extern "C" {
+# define __END_DECLS }
+#else
+# define __BEGIN_DECLS /* empty */
+# define __END_DECLS /* empty */
+#endif
+
+
+__BEGIN_DECLS
+/*====================================================================
+ * FUNCTION
+ *====================================================================*/
+
+/**
+ \struct gsl_gamin_function_struct
+ \brief Definition of an arbitrary real-valued function with gsl_vector
+ input and some parameters.
+ */
+struct gsl_gamin_function_struct
+{
+ /** objective function to minimize. */
+ double (* f) (const gsl_vector * x, void * params);
+
+ /** Size of the objective function. */
+ size_t n;
+
+ /** Possible parameters to the objective function. */
+ void * params;
+};
+
+/**
+ \typedef gsl_gamin_function
+ */
+typedef struct gsl_gamin_function_struct gsl_gamin_function;
+
+#define GSL_GAMIN_FN_EVAL(F,x) (*((F)->f))(x,(F)->params)
+
+
+/*====================================================================
+ * FMINIMIZER
+ *====================================================================*/
+
+/**
+ \struct gsl_gamin_fminimizer_type
+ \brief Defines general interface of the minimizer type.
+ */
+typedef struct
+{
+ /** Minimizer type name. */
+ const char *name;
+
+ void (*get_candidate)(const size_t i,
+ const size_t n,
+ const size_t NP,
+ const gsl_vector ** parents,
+ const double F,
+ const double CR,
+ const gsl_rng * r,
+ const size_t best,
+ gsl_vector * candidate);
+}
+gsl_gamin_fminimizer_type;
+
+/**
+ \struct gsl_gamin_fminimizer
+ \brief TODO.
+*/
+typedef struct
+{
+ const gsl_gamin_fminimizer_type * type;
+
+ size_t NP;
+
+ const gsl_rng_type * rngT;
+ gsl_rng * r;
+
+ gsl_vector ** parents;
+ gsl_vector ** children;
+
+ double * cost;
+
+ gsl_vector * candidate;
+
+ gsl_gamin_function * f;
+
+ size_t n;
+
+ double F;
+ double CR;
+
+ double low;
+ double high;
+
+ size_t best;
+}
+gsl_gamin_fminimizer;
+
+
+/**
+ * This function returns a pointer to a newly allocated instance of a
+ * minimizer of type T for an n-dimension function. If there is
+ * insufficient memory to create the minimizer then the function returns
+ * a null pointer and the error handler is invoked with an error code of
+ * GSL_ENOMEM.
+ */
+gsl_gamin_fminimizer *
+gsl_gamin_fminimizer_alloc(const gsl_gamin_fminimizer_type *T,
+ size_t n, size_t NP);
+
+/**
+ * This function initializes the minimizer s to minimize the function f,
+ * starting from the initial population normally distributed around
+ * point x. If x is NULL, uniform distribution is choosen.
+*/
+int
+gsl_gamin_fminimizer_initialize(gsl_gamin_fminimizer * s,
+ gsl_gamin_function * f,
+ const gsl_vector * x,
+ const double F,
+ const double CR,
+ const double low,
+ const double high);
+
+/**
+ * Free Storage.
+ */
+void
+gsl_gamin_fminimizer_free(gsl_gamin_fminimizer *s);
+
+/**
+ * Returns name of the method.
+ */
+const char *
+gsl_gamin_fminimizer_name(const gsl_gamin_fminimizer * s);
+
+/**
+ * This function performs a single generation iteration of the minimizer s. If
+ * the iteration encounters an unexpected problem then an error code
+ * will be returned.
+ */
+int
+gsl_gamin_fminimizer_iterate(gsl_gamin_fminimizer *s);
+
+/**
+ * Returns current best x found.
+ */
+gsl_vector *
+gsl_gamin_fminimizer_xmin(const gsl_gamin_fminimizer * s);
+
+/**
+ * Returns current best minimum found.
+ */
+double
+gsl_gamin_fminimizer_fmin(const gsl_gamin_fminimizer * s);
+
+
+/* STRATEGIES */
+
+/**
+ DE/Best/1/Exp (aka DE0) strategy. May misconverge.
+ */
+GSL_VAR const gsl_gamin_fminimizer_type * gsl_gamin_fminimizer_best1exp;
+
+/**
+ DE/Rand/1/Exp (aka DE1) strategy.
+ */
+GSL_VAR const gsl_gamin_fminimizer_type * gsl_gamin_fminimizer_rand1exp;
+
+/**
+ DE/Rand-to-Best/1/Exp strategy.one of the best strategies.
+ */
+GSL_VAR const gsl_gamin_fminimizer_type * gsl_gamin_fminimizer_randtobest1exp;
+
+/**
+ DE/Best/2/Exp strategy. another powerful strategy worth trying.
+ */
+GSL_VAR const gsl_gamin_fminimizer_type * gsl_gamin_fminimizer_best2exp;
+
+/**
+ DE/Rand/2/Exp strategy. robust optimizer for many functions.
+ */
+GSL_VAR const gsl_gamin_fminimizer_type * gsl_gamin_fminimizer_rand2exp;
+
+
+/*====================================================================
+ * TERMINATOR_GENERATION
+ *====================================================================*/
+
+/**
+ \struct gsl_gamin_terminator_generation
+ \brief Use it when you want a fixed number of generation.
+*/
+typedef struct
+{
+ size_t genCountMax;
+}
+gsl_gamin_terminator_generation;
+
+gsl_gamin_terminator_generation *
+gsl_gamin_terminator_generation_new(const size_t genCountMax);
+
+void
+gsl_gamin_terminator_generation_free(gsl_gamin_terminator_generation * t);
+
+int
+gsl_gamin_terminator_generation_check(gsl_gamin_terminator_generation * t,
+ const size_t genCount);
+
+/*====================================================================
+ * TERMINATOR_CONVERGENCE
+ *====================================================================*/
+
+/**
+ \struct gsl_gamin_terminator_convergence
+ \brief Use it when you want to stop when you are close enough from
+ a known value for the minimum.
+*/
+typedef struct
+{
+ double fconvergence;
+ double epsilon;
+}
+gsl_gamin_terminator_convergence;
+
+gsl_gamin_terminator_convergence *
+gsl_gamin_terminator_convergence_new(const double fconvergence,
+ const double epsilon);
+
+void
+gsl_gamin_terminator_convergence_free(gsl_gamin_terminator_convergence * t);
+
+int
+gsl_gamin_terminator_convergence_check(gsl_gamin_terminator_convergence * t,
+ const double fcurrent);
+
+/*====================================================================
+ * TERMINATOR_POPULATION_CONVERGENCE
+ *====================================================================*/
+
+/**
+ \struct gsl_gamin_terminator_population_convergence
+ \brief Use it when you want the population to converge.
+*/
+typedef struct
+{
+ size_t size;
+ double * fbuffer;
+ double precision;
+ size_t internal_counter;
+}
+gsl_gamin_terminator_pop_conv;
+
+gsl_gamin_terminator_pop_conv *
+gsl_gamin_terminator_pop_conv_new(const size_t size,
+ const double precision);
+
+void
+gsl_gamin_terminator_pop_conv_free(gsl_gamin_terminator_pop_conv * t);
+
+int
+gsl_gamin_terminator_pop_conv_check(gsl_gamin_terminator_pop_conv * t,
+ const double fcurrent);
+
+__END_DECLS
+
+#endif /* __GSL_GAMIN_H__ */
diff -Nur ./gamin/minimizer.c ../gsl-gamin/gamin/minimizer.c
--- ./gamin/minimizer.c 1970-01-01 01:00:00.000000000 +0100
+++ ../gsl-gamin/gamin/minimizer.c 2005-03-24 16:59:32.363154520 +0100
@@ -0,0 +1,294 @@
+/*#include <config.h>*/
+
+#include <time.h>
+
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_machine.h>
+#include <gsl/gsl_types.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_randist.h>
+
+#include "gsl_gamin.h"
+
+gsl_gamin_fminimizer *
+gsl_gamin_fminimizer_alloc(const gsl_gamin_fminimizer_type * T,
+ size_t n, size_t NP)
+{
+ size_t i,j;
+
+ gsl_gamin_fminimizer *s =
+ (gsl_gamin_fminimizer *) malloc(sizeof(gsl_gamin_fminimizer));
+
+ if (!s)
+ {
+ GSL_ERROR_VAL("failed to allocate space for minimizer struct",
+ GSL_ENOMEM, 0);
+ }
+
+ s->type = T;
+
+ s->NP = NP;
+ s->n = n;
+
+ /* Allocate Random number generator */
+ s->rngT = gsl_rng_ranlxd2;
+ s->r = gsl_rng_alloc(s->rngT);
+ + if(!s->r)
+ {
+ free(s);
+ GSL_ERROR_VAL("failed to allocate space for RNG r",
+ GSL_ENOMEM, 0);
+ }
+
+ /* Allocate temporary storage for parents population. */
+ s->parents = malloc((s->NP) * sizeof(gsl_vector*));
+ + if(!s->parents)
+ {
+ gsl_rng_free(s->r);
+ free(s);
+ GSL_ERROR_VAL("failed to allocate space for parents",
+ GSL_ENOMEM, 0);
+ }
+
+ for(i=0; i<(s->NP); ++i)
+ {
+ s->parents[i] = gsl_vector_alloc(s->n);
+
+ if(!s->parents[i])
+ {
+ for(j=0; j<i; ++j) gsl_vector_free(s->parents[j]);
+ free(s->parents);
+ gsl_rng_free(s->r);
+ free(s);
+ GSL_ERROR_VAL("failed to allocate space for parents[]",
+ GSL_ENOMEM, 0);
+ }
+ }
+
+
+ /* Allocate temporary storage for children population. */
+ s->children = malloc((s->NP) * sizeof(gsl_vector*));
+
+ if(!s->children)
+ {
+ for(j=0; j<(s->NP); ++j) gsl_vector_free(s->parents[j]);
+ free(s->parents);
+ gsl_rng_free(s->r);
+ free(s);
+ GSL_ERROR_VAL("failed to allocate space for children",
+ GSL_ENOMEM, 0);
+ }
+
+ for(i=0; i<s->NP; ++i)
+ {
+ s->children[i] = gsl_vector_alloc(s->n);
+
+ if(!s->children[i])
+ {
+ for(j=0; j<i; ++j) gsl_vector_free(s->children[j]);
+ free(s->children);
+ for(j=0; j<(s->NP); ++j) gsl_vector_free(s->parents[j]);
+ free(s->parents);
+ gsl_rng_free(s->r);
+ free(s);
+ GSL_ERROR_VAL("failed to allocate space for children[]",
+ GSL_ENOMEM, 0);
+ }
+ }
+
+
+ /* Allocate cost values */
+ s->cost = malloc((s->NP) * sizeof(double));
+
+ if(!s->cost)
+ {
+ for(j=0; j<(s->NP); ++j) gsl_vector_free(s->children[j]);
+ free(s->children);
+ for(j=0; j<(s->NP); ++j) gsl_vector_free(s->parents[j]);
+ free(s->parents);
+ gsl_rng_free(s->r);
+ free(s);
+ GSL_ERROR_VAL("failed to allocate space for cost",
+ GSL_ENOMEM, 0);
+ }
+
+
+ /* Allocate temporary storage for candidate */
+ s->candidate = gsl_vector_alloc(s->n);
+
+ if(!s->candidate)
+ {
+ free(s->cost);
+ for(j=0; j<(s->NP); ++j) gsl_vector_free(s->children[j]);
+ free(s->children);
+ for(j=0; j<(s->NP); ++j) gsl_vector_free(s->parents[j]);
+ free(s->parents);
+ gsl_rng_free(s->r);
+ free(s);
+ GSL_ERROR_VAL("failed to allocate space for candidate",
+ GSL_ENOMEM, 0);
+ }
+
+ return s;
+}
+
+int
+gsl_gamin_fminimizer_initialize(gsl_gamin_fminimizer * s,
+ gsl_gamin_function * f,
+ const gsl_vector * x,
+ const double F,
+ const double CR,
+ const double low,
+ const double high)
+{
+ size_t i,j;
+ double best_cost;
+
+ if(s->n != f->n)
+ {
+ GSL_ERROR ("function incompatible with solver size", GSL_EBADLEN);
+ }
+ + if(x != NULL && x->size != s->n)
+ {
+ GSL_ERROR("vector length not compatible with solver size", GSL_EBADLEN);
+ } +
+/* Initialize the random generator. */
+ gsl_rng_set(s->r, (unsigned long int)time(0));
+ + s->f = f;
+ s->F = F;
+ s->CR = CR;
+ s->low = low;
+ s->high = high;
+ + /* Create random individu in the population */
+ if(x == NULL)
+ { + /* uniform initialization */
+ for(i=0; i<(s->NP); ++i)
+ for(j=0; j<(s->n); ++j)
+ gsl_vector_set(s->parents[i], j,
+ low + gsl_rng_uniform(s->r)*(high-low));
+ }
+ else
+ {
+ /* Normal distribution around x
+ if x[j] = 0 then standard deviation is 1 else is x[j]/10.*/
+ for(i=0; i<(s->NP); ++i)
+ for(j=0; j<(s->n); ++j)
+ {
+ if(gsl_vector_get(x,j) > GSL_DBL_EPSILON)
+ gsl_vector_set(s->parents[i], j,
+ gsl_vector_get(x,j) +
+ gsl_ran_gaussian(s->r,
+ gsl_vector_get(x,j)/10.));
+ else
+ gsl_vector_set(s->parents[i], j,
+ gsl_ran_gaussian(s->r, 1.));
+ }
+ }
+ + /* Compute cost for the current population */
+ for(i=0; i<(s->NP); ++i)
+ s->cost[i] = GSL_GAMIN_FN_EVAL(s->f, s->parents[i]);
+ +
+ /* Find best individual */
+ best_cost = s->cost[0];
+ s->best = 0;
+
+ for(i=1; i<(s->NP); ++i)
+ {
+ if(s->cost[i] < best_cost)
+ {
+ best_cost = s->cost[i];
+ s->best = i;
+ }
+ }
+
+return GSL_SUCCESS;
+}
+
+void
+gsl_gamin_fminimizer_free(gsl_gamin_fminimizer * s)
+{
+ size_t j;
+ + gsl_rng_free(s->r);
+ + for(j=0; j<(s->NP); ++j) gsl_vector_free(s->parents[j]);
+ free(s->parents);
+ + for(j=0; j<(s->NP); ++j) gsl_vector_free(s->children[j]);
+ free(s->children);
+
+ free(s->cost);
+
+ gsl_vector_free(s->candidate);
+
+ free(s);
+}
+
+int
+gsl_gamin_fminimizer_iterate(gsl_gamin_fminimizer * s)
+{
+ size_t i;
+ double candidate_cost;
+ gsl_vector ** tmp_swap;
+
+ /* Population Loop */
+ for(i=0; i<(s->NP); ++i)
+ {
+ (s->type->get_candidate)(i, s->n, s->NP,
+ (const gsl_vector **)s->parents,
+ s->F, s->CR, s->r, s->best, s->candidate);
+
+ candidate_cost = GSL_GAMIN_FN_EVAL(s->f, s->candidate);
+ + if(candidate_cost < s->cost[i])
+ {
+ /* candidate is better so we replace */
+ s->cost[i] = candidate_cost;
+ + gsl_vector_swap(s->candidate, s->children[i]);
+ + if(candidate_cost < s->cost[s->best])
+ s->best = i;
+ }
+ else
+ /* individual i is better so the take him */
+ gsl_vector_memcpy(s->children[i], s->parents[i]);
+ }
+ + /* swap population*/
+ tmp_swap = s->parents;
+ s->parents = s->children;
+ s->children = tmp_swap;
+ + return GSL_CONTINUE;
+
+}
+
+const char *
+gsl_gamin_fminimizer_name(const gsl_gamin_fminimizer * s)
+{
+ return s->type->name;
+}
+
+
+gsl_vector *
+gsl_gamin_fminimizer_xmin(const gsl_gamin_fminimizer * s)
+{
+ return s->parents[s->best];
+}
+
+double
+gsl_gamin_fminimizer_fmin(const gsl_gamin_fminimizer * s)
+{
+ return s->cost[s->best];
+}
diff -Nur ./gamin/rand1exp.c ../gsl-gamin/gamin/rand1exp.c
--- ./gamin/rand1exp.c 1970-01-01 01:00:00.000000000 +0100
+++ ../gsl-gamin/gamin/rand1exp.c 2005-03-24 16:59:32.363154520 +0100
@@ -0,0 +1,77 @@
+/* #include <config.h> */
+
+#include <stdlib.h>
+#include <string.h>
+
+#include <stdio.h>
+
+#include "gsl_gamin.h"
+
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_randist.h>
+
+
+static void
+rand1exp_get_candidate(const size_t i,
+ const size_t n,
+ const size_t NP,
+ const gsl_vector ** parents,
+ const double F,
+ const double CR,
+ const gsl_rng * r,
+ const size_t best,
+ gsl_vector * candidate)
+{
+ size_t r1, r2, r3, k, L;
+
+ /* Avoid unnecessary warning for unused best parameter */
+ size_t idum_;
+ idum_= best;
+ /* End of trick */
+
+ do /* Pick a random population member */
+ { /* Endless loop for NP < 2 !!! */
+ r1 = gsl_rng_uniform_int(r, NP);
+ }
+ while(r1==i); + + do /* Pick a random population member */
+ { /* Endless loop for NP < 3 !!! */
+ r2 = gsl_rng_uniform_int(r, NP);
+ }
+ while((r2==i) || (r2==r1));
+ + do /* Pick a random population member */
+ { /* Endless loop for NP < 4 !!! */
+ r3 = gsl_rng_uniform_int(r, NP);
+ }
+ while((r3==i) || (r3==r1)|| (r3==r2) );
+ + + gsl_vector_memcpy(candidate, parents[i]);
+
+ k = gsl_rng_uniform_int(r, n);
+ L = 0;
+ do
+ {
+ gsl_vector_set(candidate, k,
+ gsl_vector_get(parents[r1],k) +
+ F*(gsl_vector_get(parents[r2],k)-
+ gsl_vector_get(parents[r3],k)));
+ k = (k+1) % n;
+ L++;
+ }
+ while( (gsl_rng_uniform(r)< CR) && (L < n) );
+}
+
+
+
+static const gsl_gamin_fminimizer_type rand1exp_type =
+ { "DE/rand/1/exp", /* name */
+ &rand1exp_get_candidate
+ };
+
+const gsl_gamin_fminimizer_type *
+gsl_gamin_fminimizer_rand1exp = &rand1exp_type;
+
diff -Nur ./gamin/rand2exp.c ../gsl-gamin/gamin/rand2exp.c
--- ./gamin/rand2exp.c 1970-01-01 01:00:00.000000000 +0100
+++ ../gsl-gamin/gamin/rand2exp.c 2005-03-24 16:59:32.363154520 +0100
@@ -0,0 +1,90 @@
+/* #include <config.h> */
+
+#include <stdlib.h>
+#include <string.h>
+
+#include <stdio.h>
+
+#include "gsl_gamin.h"
+
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_randist.h>
+
+
+static void
+rand2exp_get_candidate(const size_t i,
+ const size_t n,
+ const size_t NP,
+ const gsl_vector ** parents,
+ const double F,
+ const double CR,
+ const gsl_rng * r,
+ const size_t best,
+ gsl_vector * candidate)
+{
+ size_t r1, r2, r3, r4, r5, k, L;
+ + /* Avoid unnecessary warning for unused best parameter */
+ size_t idum_;
+ idum_= best;
+ /* End of trick */
+
+ do /* Pick a random population member */
+ { /* Endless loop for NP < 2 !!! */
+ r1 = gsl_rng_uniform_int(r, NP);
+ }
+ while(r1==i); + + do /* Pick a random population member */
+ { /* Endless loop for NP < 3 !!! */
+ r2 = gsl_rng_uniform_int(r, NP);
+ }
+ while((r2==i) || (r2==r1));
+ + do /* Pick a random population member */
+ { /* Endless loop for NP < 4 !!! */
+ r3 = gsl_rng_uniform_int(r, NP);
+ }
+ while((r3==i) || (r3==r1)|| (r3==r2) );
+ + do /* Pick a random population member */
+ { /* Endless loop for NP < 5 !!! */
+ r4 = gsl_rng_uniform_int(r, NP);
+ }
+ while((r4==i) || (r4==r1) || (r4==r2) || (r4==r3));
+
+ do /* Pick a random population member */
+ { /* Endless loop for NP < 6 !!! */
+ r5 = gsl_rng_uniform_int(r, NP);
+ }
+ while((r5==i) || (r5==r1) || (r5==r2) || (r5==r3) || (r5==r4));
+
+ gsl_vector_memcpy(candidate, parents[i]);
+
+ k = gsl_rng_uniform_int(r, n);
+ L = 0;
+ do
+ {
+ gsl_vector_set(candidate, k,
+ gsl_vector_get(parents[r1],k) +
+ F*(gsl_vector_get(parents[r2],k)+
+ gsl_vector_get(parents[r3],k)-
+ gsl_vector_get(parents[r4],k)-
+ gsl_vector_get(parents[r5],k)));
+ k = (k+1) % n;
+ L++;
+ }
+ while( (gsl_rng_uniform(r)< CR) && (L < n) );
+}
+
+
+
+static const gsl_gamin_fminimizer_type rand2exp_type =
+ { "DE/rand/2/exp", /* name */
+ &rand2exp_get_candidate
+ };
+
+const gsl_gamin_fminimizer_type *
+gsl_gamin_fminimizer_rand2exp = &rand2exp_type;
+
diff -Nur ./gamin/randtobest1exp.c ../gsl-gamin/gamin/randtobest1exp.c
--- ./gamin/randtobest1exp.c 1970-01-01 01:00:00.000000000 +0100
+++ ../gsl-gamin/gamin/randtobest1exp.c 2005-03-24 16:59:32.363154520 +0100
@@ -0,0 +1,67 @@
+/* #include <config.h> */
+
+#include <stdlib.h>
+#include <string.h>
+
+#include <stdio.h>
+
+#include "gsl_gamin.h"
+
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_randist.h>
+
+
+static void
+randtobest1exp_get_candidate(const size_t i,
+ const size_t n,
+ const size_t NP,
+ const gsl_vector ** parents,
+ const double F,
+ const double CR,
+ const gsl_rng * r,
+ const size_t best,
+ gsl_vector * candidate)
+{
+ size_t r1, r2, k, L;
+ + do /* Pick a random population member */
+ { /* Endless loop for NP < 2 !!! */
+ r1 = gsl_rng_uniform_int(r, NP);
+ }
+ while(r1==i); + + do /* Pick a random population member */
+ { /* Endless loop for NP < 3 !!! */
+ r2 = gsl_rng_uniform_int(r, NP);
+ }
+ while((r2==i) || (r2==r1));
+ + gsl_vector_memcpy(candidate, parents[i]);
+
+ k = gsl_rng_uniform_int(r, n);
+ L = 0;
+ do
+ {
+ gsl_vector_set(candidate, k,
+ gsl_vector_get(candidate,k) +
+ F*(gsl_vector_get(parents[best],k)-
+ gsl_vector_get(candidate,k))+
+ F*(gsl_vector_get(parents[r1],k)-
+ gsl_vector_get(parents[r2],k)));
+ k = (k+1) % n;
+ L++;
+ }
+ while( (gsl_rng_uniform(r)< CR) && (L < n) );
+}
+
+
+
+static const gsl_gamin_fminimizer_type randtobest1exp_type =
+ { "DE/rand-to-best/1/exp", /* name */
+ &randtobest1exp_get_candidate
+ };
+
+const gsl_gamin_fminimizer_type *
+gsl_gamin_fminimizer_randtobest1exp = &randtobest1exp_type;
+
diff -Nur ./gamin/terminator_convergence.c ../gsl-gamin/gamin/terminator_convergence.c
--- ./gamin/terminator_convergence.c 1970-01-01 01:00:00.000000000 +0100
+++ ../gsl-gamin/gamin/terminator_convergence.c 2005-03-24 16:59:32.363154520 +0100
@@ -0,0 +1,42 @@
+/*#include <config.h>*/
+#include <gsl/gsl_machine.h>
+
+#include "gsl_gamin.h"
+
+#include <math.h>
+
+gsl_gamin_terminator_convergence *
+gsl_gamin_terminator_convergence_new(const double fconvergence,
+ const double epsilon)
+{
+ gsl_gamin_terminator_convergence * t =
+ malloc(sizeof(gsl_gamin_terminator_convergence));
+
+ if(!t)
+ {
+ GSL_ERROR_VAL("failed to allocate space for terminator_convergence",
+ GSL_ENOMEM, 0);
+ }
+
+ t->fconvergence = fconvergence;
+ t->epsilon = epsilon;
+
+ return t;
+}
+
+void
+gsl_gamin_terminator_convergence_free(gsl_gamin_terminator_convergence * t)
+{
+ free(t);
+}
+
+int
+gsl_gamin_terminator_convergence_check(gsl_gamin_terminator_convergence * t,
+ const double fcurrent)
+{
+ if(fabs(fcurrent - t->fconvergence)/
+ (GSL_DBL_EPSILON + fabs(t->fconvergence)) > t->epsilon)
+ return GSL_CONTINUE;
+
+ return GSL_SUCCESS;
+}
diff -Nur ./gamin/terminator_generation.c ../gsl-gamin/gamin/terminator_generation.c
--- ./gamin/terminator_generation.c 1970-01-01 01:00:00.000000000 +0100
+++ ../gsl-gamin/gamin/terminator_generation.c 2005-03-24 16:59:32.363154520 +0100
@@ -0,0 +1,34 @@
+/*#include <config.h>*/
+
+#include "gsl_gamin.h"
+
+gsl_gamin_terminator_generation *
+gsl_gamin_terminator_generation_new(const size_t genCountMax)
+{
+ gsl_gamin_terminator_generation * t =
+ malloc(sizeof(gsl_gamin_terminator_generation));
+
+ if(!t)
+ {
+ GSL_ERROR_VAL("failed to allocate space for terminator_generation",
+ GSL_ENOMEM, 0);
+ }
+
+ t->genCountMax = genCountMax;
+
+ return t;
+}
+
+void
+gsl_gamin_terminator_generation_free(gsl_gamin_terminator_generation * t)
+{
+ free(t);
+}
+
+int
+gsl_gamin_terminator_generation_check(gsl_gamin_terminator_generation * t,
+ const size_t genCount)
+{
+ if(genCount < t->genCountMax) return GSL_CONTINUE;
+ return GSL_SUCCESS;
+}
diff -Nur ./gamin/terminator_pop_conv.c ../gsl-gamin/gamin/terminator_pop_conv.c
--- ./gamin/terminator_pop_conv.c 1970-01-01 01:00:00.000000000 +0100
+++ ../gsl-gamin/gamin/terminator_pop_conv.c 2005-03-24 16:59:32.363154520 +0100
@@ -0,0 +1,66 @@
+/*#include <config.h>*/
+#include <gsl/gsl_machine.h>
+#include <gsl/gsl_statistics.h>
+#include "gsl_gamin.h"
+
+#include <math.h>
+
+gsl_gamin_terminator_pop_conv *
+gsl_gamin_terminator_pop_conv_new(const size_t size,
+ const double precision)
+{
+ gsl_gamin_terminator_pop_conv * t =
+ malloc(sizeof(gsl_gamin_terminator_pop_conv));
+
+ if(!t)
+ {
+ GSL_ERROR_VAL("failed to allocate space for terminator_pop_conv",
+ GSL_ENOMEM, 0);
+ }
+
+ t->size = size;
+ t->precision = precision;
+
+ t->fbuffer = (double *)malloc(t->size * sizeof(double));
+
+ if(!t->fbuffer)
+ {
+ free(t);
+ GSL_ERROR_VAL("failed to allocate space for fbuffer",
+ GSL_ENOMEM, 0);
+ }
+ + t->internal_counter = 0;
+
+ return t;
+}
+
+void
+gsl_gamin_terminator_pop_conv_free(gsl_gamin_terminator_pop_conv * t)
+{
+ free(t->fbuffer);
+ free(t);
+}
+
+int
+gsl_gamin_terminator_pop_conv_check(gsl_gamin_terminator_pop_conv * t,
+ const double fcurrent)
+{
+ double mean, std;
+
+ /* circulate into fbuffer */
+ t->fbuffer[t->internal_counter % (t->size)] = fcurrent;
+
+ ++(t->internal_counter);
+
+ if(t->internal_counter >= t->size)
+ {
+ mean = gsl_stats_mean(t->fbuffer, (size_t)1, t->size);
+ std = gsl_stats_sd_m(t->fbuffer, (size_t)1, t->size, mean);
+
+ if(std/(GSL_DBL_EPSILON+fabs(mean)) < t->precision) return GSL_SUCCESS;
+ }
+
+
+ return GSL_CONTINUE;
+}
diff -Nur ./gamin/test.c ../gsl-gamin/gamin/test.c
--- ./gamin/test.c 1970-01-01 01:00:00.000000000 +0100
+++ ../gsl-gamin/gamin/test.c 2005-03-24 16:59:32.363154520 +0100
@@ -0,0 +1,79 @@
+#include <config.h>
+#include <math.h>
+#include <gsl/gsl_test.h>
+#include <gsl/gsl_gamin.h>
+
+double twod(const gsl_vector * x, void * params)
+{
+ double result=0;
+
+ double pi=3.1415926536,sigma2=0.15;
+ int nn=9;
+
+ double x0 = gsl_vector_get(x, (size_t)0);
+ double x1 = gsl_vector_get(x, (size_t)1);
+
+ if(x0>1 || x1>1) return 10000000.;
+
+ result = sqrt( (0.5-x0)*(0.5-x0)+ (0.5-x1)*(0.5-x1) );
+
+ result = -cos(result*nn*pi)*cos(result*nn*pi)*exp(-result*result/sigma2);
+
+ + return result;
+}
+
+int main()
+{
+ int iter = 0;
+ int status;
+ const size_t NG = 5000;
+ const double F = 0.8;
+ const double CR = 0.9;
+ const double low = -100;
+ const double high = 100.;
+
+ const gsl_gamin_fminimizer_type * T = gsl_gamin_fminimizer_rand1exp;
+ + gsl_gamin_fminimizer *s;
+
+ gsl_gamin_terminator_generation * t1 =
+ gsl_gamin_terminator_generation_new(NG);
+
+ gsl_gamin_function my_func;
+ my_func.f = &twod;
+ my_func.n = 2;
+ my_func.params = NULL;
+
+ s = gsl_gamin_fminimizer_alloc(T, 2, 60);
+ + gsl_gamin_fminimizer_initialize(s, &my_func, NULL,
+ F, CR,
+ low, high);
+ + do
+ {
+ iter++;
+ status = gsl_gamin_fminimizer_iterate(s);
+ }
+ while(status == GSL_CONTINUE &&
+ gsl_gamin_terminator_generation_check(t1, iter) == GSL_CONTINUE);
+
+
+ if( fabs(gsl_gamin_fminimizer_fmin(s)-(-1.) < 1.e-3) )
+ status = GSL_SUCCESS;
+ else status = GSL_FAILURE;
+
+ gsl_test(status, "%s on twod in %d iterations.",
+ gsl_gamin_fminimizer_name(s),
+ iter);
+
+/* printf("twod(%f,%f) = %f\n", */
+/* gsl_vector_get(gsl_gamin_fminimizer_xmin(s), 0), */
+/* gsl_vector_get(gsl_gamin_fminimizer_xmin(s), 1), */
+/* gsl_gamin_fminimizer_fmin(s)); */
+
+ gsl_gamin_fminimizer_free(s);
+
+ return 0;
+}





Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]