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

dmitrey openopt at ukr.net
Mon Jul 30 15:26:52 EDT 2007


I guess yes, but didn't try that ones, because redefining 
multipack_jac_transpose will yield even more misunderstood than having 2 
MATRIXC2F funcs from different scipy toolkits. multipack_jac_transpose = 
1 if and only if user-defined variable col_deriv is set to 1, and I'm 
sure this should be remained.
from leastsq docstring:

      col_deriv -- non-zero to specify that the Jacobian function
                   computes derivatives down the columns (faster, because
                   there is no transpose operation).

(default col_deriv value is zero)

so col_deriv=0 means "no transpose", col_deriv=1 means "do transpose".
Regards, D.


Bruce Southey wrote:
> Hi,
>
> I would agree that having multiple different MATRIXC2F is incorrect
> which is why I made the reply.  To me (and which is probably
> incorrect), MATRIXC2F is just converting a C matrix to a Fortran
> matrix.  I do think that the multipack_jac_transpose flag is being
> real issue so it would be easier to change the condition that calls
> MATRIXC2F rather than it's definition.
>
> So does changing the condition (multipack_jac_transpose==1) prior to
> calling MATRIXC2F have the same result as redefining MATRIXC2F?
>
> Bruce
>
>
> On 7/30/07, dmitrey <openopt at ukr.net> wrote:
>   
>> 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 at jenningsstory.com> wrote:
>>>
>>>       
>>>> On 7/26/07, Alan G Isaac <aisaac at 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 at scipy.org
>>>> http://projects.scipy.org/mailman/listinfo/scipy-dev
>>>>
>>>>
>>>>         
>>> _______________________________________________
>>> Scipy-dev mailing list
>>> Scipy-dev at scipy.org
>>> http://projects.scipy.org/mailman/listinfo/scipy-dev
>>>
>>>
>>>
>>>
>>>       
>> _______________________________________________
>> Scipy-dev mailing list
>> Scipy-dev at scipy.org
>> http://projects.scipy.org/mailman/listinfo/scipy-dev
>>
>>     
> _______________________________________________
> Scipy-dev mailing list
> Scipy-dev at scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-dev
>
>
>
>   




More information about the SciPy-Dev mailing list