2

I'm working from the template program given here:

https://www.gnu.org/software/gsl/manual/html_node/Trivial-example.html

The program as they give it compiles and runs perfectly, which is nice. What I would like to do is generalise this method to find the minimum of a function with an arbitrary number of parameters.

Some cursory reading suggests that the metric function (M1) is only used in certain diagnostic and printing situations and so can more or less be ignored. All that remains is then to define E1 and S1 appropriately. Unfortunately my knowledge of using pointers and void is incomplete, so I'm stuck trying to upgrade the configuration 'xp' to be an array of parameters, rather than a single double.

In my naivete tried moving from 

double x =  *((double *) xp);

to 

double x =  (*((double *) xp))[0];

where appropriate, but obviously that didn't work. I'm sure I'm missing something stupid, so any hints would be nice! I will obviously be defining my own E1 output function which will take these N parameters and return a number.

zylatis
  • 448
  • 3
  • 14

1 Answers1

0

The underlying algorithm, gsl_siman_solve() from the link provided, is generalized to work with any data type. This is why the ubiquitous xp parameter is always being cast to a double pointer before use. It should be straightforward to use any struct or array or array of arrays instead of simply doubles provided all the callbacks are coded properly.

The problem is that gsl_siman_solve() only seems to support a scalar double step size, initial guess, and 'uniform' value (from gsl_rng_uniform()), so you would need to map scalar double values into what are naturally multidimensional quantities. This can be done, but it is messy and not very flexible. In your case, the mapping would be done in S1().

This is akin to mapping the digits of a decimal number into a multidimensional space: the ones digit represents the X axis, the tens digit represents the Y axis, and the hundreds digit represents the Z axis, for example. By incrementing an integer, one can walk the entire 3D space from (0, 0, 0) to (9, 9, 9). You don't have to use integers and powers of 10, and the components don't even have to have the same range, but there is an inherent limit in the range of each component of the packed value. You would actually do this in reverse: taking a scalar double and unpacking it into multiple quantities.

Lastly, your code double x = (*((double *) xp))[0]; won't work because you are attempting to dereference a double as an array, not a pointer to a double, which would be OK. In other words, it's that first * that is the problem.

Randall Cook
  • 6,728
  • 6
  • 33
  • 68
  • Thanks Randall. Given that it is, as you say, messy and not very flexible, could you recommend any libraries that make it a bit easier to search a higher dimensional parameter space? – zylatis Apr 15 '15 at 11:41
  • Ack, can't edit comments. I found the quick tutorial on ADA given here very useful. Loads of options for fiddling, but also the ability to just bung in your own function and bounds and run with it, very cool! http://www.quantcode.com/modules/mydownloads/singlefile.php?lid=504 – zylatis Apr 15 '15 at 13:19
  • @zylatis, glad you found something helpful. I think you are right to go with a more flexible library. After all, it is the Simulated Annealing algorithm you want, not the library it is packaged in. And for reference, you can edit comments on Stack Exchange, but only for a short time, like 5 minutes. – Randall Cook Apr 15 '15 at 18:45
  • Ah right, got it thanks. Interestingly, and also annoyingly, this new library sometimes gets a lower solution than my mathematica code, sometimes not, depending on the model. Time to tinker I guess. – zylatis Apr 15 '15 at 18:55