[SciPy-dev] BUG gradient in leastsq (http://projects.scipy.org/scipy/scipy/ticket/416)

Alan G Isaac aisaac@american....
Thu Jul 26 10:15:40 CDT 2007


On Thu, 26 Jul 2007, dmitrey apparently wrote:
> 1. What's the difference between these 2 funcs from __minpack.h:
> int jac_multipack_calling_function(int *n, double *x, double *fvec, 
> double *fjac, int *ldfjac, int *iflag)
> int jac_multipack_lm_function(int *m, int *n, double *x, double *fvec, 
> double *fjac, int *ldfjac, int *iflag)

> They have same description. 

>   /* This is the function called from the Fortran code it should
>         -- use call_python_function to get a multiarrayobject result
>     -- check for errors and return -1 if any
>     -- otherwise place result of calculation in *fvec or *fjac.

>      If iflag = 1 this should compute the function.
>      If iflag = 2 this should compute the jacobian (derivative matrix)
>   */

> 2. So patch assigned to the ticket proposes to rewrite the line 152 
> MATRIXC2F(fjac, result_array->data, *n, *ldfjac)
> as 
> MATRIXC2F(fjac, result_array->data, *ldfjac, *n)

> however, line 92 is the same. so maybe it needs same patch. 

> 3. The MATRIXC2F is defined in minpack.h w/o any description:
> #define MATRIXC2F(jac,data,n,m) {double *p1=(double *)(jac), *p2, 
> *p3=(double *)(data);\ 
> int i,j;\ 
> for (j=0;j<(m);p3++,j++) \ 
>   for (p2=p3,i=0;i<(n);p2+=(m),i++,p1++) \
>     *p1 = *p2; }

> I have no idea what does it do. 
> So I replaced (jac,data,n,m) by (jac,data,m,n), and user's example works 
> correctly for all cases 1-3:
> 1) w/o gradient info 
> 2) with gradient info, col_deriv=0 
> 3) with gradient info, col_deriv=1 (in the case I modified the user's 
> gradient func so that it returns transposed gradient)

> scipy.test(1) also didn't yield any bugs related to leastsq. 



Andy, and others, can you comment on this?
Also, I would hope for Travis to comment
before committing these changes, since he
contributed the code.

Here is the situation: 
http://projects.scipy.org/scipy/scipy/ticket/416
Andy identified a bug and a possible fix.
(But see details above.)
Neither Dmitrey nor I are familiar with this code.
Dmitrey is willing to apply the patch if there is support 
for his doing so, but it will be done "mechancially".
In any case, he will add a unit test exposing the problem.

Cheers,
Alan Isaac






More information about the Scipy-dev mailing list