I am using the multifit_nlin module from pygsl for nonlinear least squares fitting. pygsl is a python binding of the c numerical library gsl. The problem that I am experiencing does not seem to be related to pygsl or gsl, but it appears in this context only.
I am fitting parameters of a function to some data. To use pygsl for parameters fitting I need to define the function and its jacobian. Then multifit_nlin's fitter lmsder calls these two function when needed in the fitting process. When, I make a call to the jacobian, it produces a matrix of numbers. I can output this matrix to screen and I see that the number are correct. Next, I define a lmsder class and initialize it with the lmsder.set command. I output the jacobian matrix with the lmsder.getJ() command to screen and I see the same numbers as before. Of course, this is not what I want to do with my code but for illustrative and debugging purposes only.
The agreement between the outputs of jacobian and lmsder.getJ() are what you would expect since lmsder.getJ() accesses the jacobian matrix in memory which was produced by the jacobian function. However, if I insert a line a code, say print 'bob" (or anything else), as in the following
system = gsl_multifit_function_fdf(...) # jacobian is passed here
solver = lmsder(...) # system is passed here
solver.set(...) # first call to jacobian is in here
print "bob"
print solver.getJ()
where ... means the appropriate arguments. Then the print solver.getJ() prints a matrix which is a transpose of the jacobian matrix with lower rows filled with random content. So again, this only happens when there are extra lines of code between the set() and getJ() calls.
If I execute my code normally, i.e. the entire fitting process that I have, the code works error free. If the jacobian matrix was indeed what the getJ() command shows then, there would pretty of places where an exception could be raised. So, I know for certain that my code works and also because the values that I get for the parameters are reasonable.
I have also tracked the chain of calls that pygsl does all the way to the gsl's c library. There is nothing that causes this problem. Also, gsl has been round for ages and something as simple as displaying a matrix would have been fixed ages ago.
Any suggestions to what might be the cause of this problem? Garbage collector, incorrect ordering of import statements, multicore? What tools can I use to check for memory leaks, garbage collection process?
Thanks, Alexander
