Feature request: pass closure into levmar

[jameshume]

Hi David,

Not sure if/how you accept feature requests. I'm happy to branch off the trunk and keep local modifications if you don't like the suggestion or do branch off then pull requests etc etc if you do. Anyway if these are preferred options let me know or if you don't accept feature requests...

I'd like to be able to pass closure functions into levmar (rather than have to use global variables)...
As a first pass I have change the levmar function header to read as follows. It only deals with closures on functions... I haven't considered object functions.

{
  /* automatically detect input function prototype */
  func_type = is_func(f);

  /* Determine the function type... */
  if (allof(func_type != [1, 5])) error, "levmar needs interpreted or closure function F";

  args = [];
  args_already_bound = 0;
  if( func_type == 5 ) {

    /* If we are a closure then remove the first argument from the list of arguments 
     * as this is bound... Also if it is a closure of a closure ... of a closure 
     * figure this out! */
    f2 = f;
    do {
      f2 = f2.function;
      args_already_bound += 1;
      func_type = is_func(f2);
    } while(func_type == 5);

    if( func_type != 1 ) error, "closures given to levmar must derive from iterpreted function F";
  }
  else
    f2 = f;

  args = strtok(print(f2)(1),"(");  delim = " ,)"; /* args(2) is the parameters of the function */
  keyd = strmatch(args(2), "deriv=");                 

  for(i = 1; i <= args_already_bound; ++i) args = strtok(args(2),delim);
  print, "## Args remaining are...", args;

  for (i=1;i<=3;++i) args = strtok(args(2),delim);
  hasd = strpart(args(1),1:1) == "&";
  ... rest of function continues without modification....

Used the rather non-comprehensive test to veryify:

func line(d1, x, a) { return a(1)*x + a(2); }
f = closure(line, 10)
x = span(1,101,200)
y = line(10, x, [10, 10])
levmar(y, x, f, [9,9]); 

func line2(d1, d2, x, a) { return a(1)*x + a(2); }
f2 = closure(closure(line2, 10), 20)
levmar(y, x, f2, [9,9]); 

[dhmunro]

I am certainly willing to consider feature requests. This is an odd one, because the "prototype detection" hack I built into levmar is already ugly. After agonizing over it for quite some time (sorry for the delay), I decided to give you a much simpler option to satisfy the need you identified: commit ae06f82 adds the levmar0 (and less useful levmar1) API, and modifies levmar itself to not complain about things which are not interpreted functions. The calling sequence for levmar0 (and levmar1) is identical to that for levmar; the difference being that levmar0 assumes the input function F conforms to the F(x,a) prototype. So all you need to do is to call levmar0 instead of levmar.

The levmar algorithm internally requires the F(x,a,&dfda) prototype; it needs the derivatives dfda once each major cycle, but subcycles use only F(x,a) without the derivative. Thus, the "natural" interface for this algorithm is rather complicated -- you need to supply a function which computes f=F(x,a) as quickly as possible, but only occasionally require the more expensive dfda calculation. I clean implementation would simply require you to supply F(x,a,&dfda) and be done with it. That is exactly what the new levmar1 function does. It's also what the levmar function itself now does in the case that the "prototype detector" hack fails, as it does with a closure, a built-in function, a wrap_args function, or anything other than a simple interpreted function. If your function doesn't adhere to this API, you simply write a small wrapper (usually one line will do it) that does.

However, this "natural" API doesn't actually cover the case you are complaining about, because your closure function doesn't support the full natural levmar API. I also provide a levmar_partial function, which attempts to compute numerical partial derivatives for F(x,a). This algorithm is not always appropriate, so using it blindly is not a very good idea. There are several external variables which you potentially need to tune for each F(x,a) (see help,levmar_partial), and the algorithm must respect the amin=, amax=, and fit= keywords to levmar itself. Instead of forcing you to write a wrapper function for each F(x,a), I supply a generic wrapper _levmar_f. I should have built levmar with a keyword specifying that you wanted to use levmar_partial (or even let you supply your own version), but instead I went with the "prototype detection" hack.

Your suggested fix would make the hack even hackier, and it still doesn't handle all the possible things you might want to pass to levmar. So I'm offering you levmar0, which I claim gets you what you need. Before this commit, levmar was certainly broken, in the sense that it forced you to write a small interpreted wrapper in order to use things that arguably were written to one of the supported protocols. You aren't forced now, but I've made the less dangerous protocol the default for levmar, and given you levmar0 to remind you that you are using the dangerous and expensive levmar_partial derivative calculator.

Feel free to add another comment if you disagree, and thank you for pointing out this problem.

[jameshume]

Hi David,

Thanks for taking the time to look at this. I agree my "fix" would make the hack even hackier - it was itself a bit of hack - and the solution proposed sounds good to me. Cheers!