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

[dmitrey]

Bruce Southey wrote:
> Hi,
> My 0.5 cents on this, as I don't use it, is that  I do think the issue
> lies elsewhere than here.
>
> In my 'old' scipy0.5.2 code, MATRIXC2F is ONLY used if
> multipack_jac_transpose ==1
>
> by the three functions:
> Lib/optimize/__minpack.h, h  jac_multipack_lm_function and
> jac_multipack_calling_function
> Lib/integrate/__odepack.h the ode_jacobian_function
>
> Also,  it may be that the expectation is not correct ie that  case 2)
> with gradient info, col_deriv=0 is switched with 3) with gradient
> info, col_deriv=1 since, my limited take on the C code, is that
> col_deriv changes multipack_jac_transpose.
>
> Bruce
>   
the func MATRIXC2F is defined twice:
in Lib/optimize/minpack.h for Lib/optimize/__minpack.h
and
in Lib/integrate/multipack.h for Lib/integrate/__odepack.h
I didn't know anything about latter, but my changes didn't affect the 
Lib/integrate package. On the other hand, now scipy has 2 separate 
definition of MATRIXC2F, one with params (jac, data, m,n) and other with 
(jac,dataq, n, m). I don't know anything is there are any mistakes in 
integrate package, but having defined two funcs MATRIXC2F with different 
args isn't very good idea.
Afaik there are no bugs for now, I checked all 3 cases mentioned:

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)

Or do you have an example of incorrect leastsq work?
If yes, please send the one.
Regards, D.


>
>
> On 7/30/07, Andy Jennings <scipy2mdjhs78c@jenningsstory.com> wrote:
>   
>> On 7/26/07, Alan G Isaac <aisaac@american.edu> wrote:
>>     
>>>  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
>>>
>>>
>>>
>>>
>>>
>>>       
>>  Thanks for looking at this bug.  The fix looks fine to me.
>>
>>  I don't really have an opinion on the jac_multipack_calling_function
>> question.  I think it's likely that it has the same issue, but on the
>> chance that it doesn't, you hate to risk breaking something that's
>> working correctly.
>> _______________________________________________
>> Scipy-dev mailing list
>> Scipy-dev@scipy.org
>> http://projects.scipy.org/mailman/listinfo/scipy-dev
>>
>>     
> _______________________________________________
> Scipy-dev mailing list
> Scipy-dev@scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-dev
>
>
>
>