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

[Bruce Southey]

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



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
>