BUG: run method of scipy.odr.ODR class fails when delta0 parameter is set

[jsdodge]

Describe your issue.

It's possible that I've misunderstood how to use scipy.odr, and if that is the case, then perhaps this bug report can be repurposed to improve the documentation.

The traceback indicates that ODR._gen_work expects data.we to be 2D. But in this example, data.we is 1D with shape (256,). (N.B.: In the traceback, I've replaced the path to my local virtual environment with ***.)

As a workaround, I tried changing the line

data = odr.RealData(x, y=y, sx=np.full_like(x, sigma_x), sy=np.full_like(x, sigma_y))

to

data = odr.RealData(x, y=y, covx=sigma_x ** 2 * np.eye(n), covy=sigma_y ** 2 * np.eye(n))

but this also fails. The traceback for this error is below (also with the virtual environment path replaced by ***), and indicates that the we array does not have the shape that odr expects.

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[3], line 25
     23 fit = odr.ODR(data, model, beta0=p_true, delta0=-x_noise)
     24 fit.set_job(del_init=1)
---> 25 result = fit.run()

File ~***/lib/python3.11/site-packages/scipy/odr/_odrpack.py:1115, in ODR.run(self)
   1112     if obj is not None:
   1113         kwds[attr] = obj
-> 1115 self.output = Output(odr(*args, **kwds))
   1117 return self.output

ValueError: could not convert we to a suitable array

Reproducing Code Example

import numpy as np
from numpy.random import default_rng
from scipy import odr


def fun(_p, _x):
    return _p[0] + _p[1] * _x

rng = default_rng(0)
p_true = (0.0, 1.0)
n = 256
x_true = np.linspace(0, 1, n)
y_true = fun(p_true, x_true)
sigma_x = 0.01
sigma_y = 0.02
x_noise = sigma_x * rng.standard_normal(n)
y_noise = sigma_y * rng.standard_normal(n)
x = x_true + x_noise
y = y_true + y_noise

model = odr.Model(fun)
data = odr.RealData(x, y=y, sx=np.full_like(x, sigma_x), sy=np.full_like(x, sigma_y))
fit = odr.ODR(data, model, beta0=p_true, delta0=-x_noise)
fit.set_job(del_init=1)
result = fit.run()

Error message

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Cell In[1], line 25
     23 fit = odr.ODR(data, model, beta0=p_true, delta0=-x_noise)
     24 fit.set_job(del_init=1)
---> 25 result = fit.run()

File ~***/lib/python3.11/site-packages/scipy/odr/_odrpack.py:1091, in ODR.run(self)
   1085 kwd_l = ['ifixx', 'ifixb', 'job', 'iprint', 'errfile', 'rptfile',
   1086          'ndigit', 'taufac', 'sstol', 'partol', 'maxit', 'stpb',
   1087          'stpd', 'sclb', 'scld', 'work', 'iwork']
   1089 if self.delta0 is not None and (self.job // 10000) % 10 == 0:
   1090     # delta0 provided and fit is not a restart
-> 1091     self._gen_work()
   1093     d0 = numpy.ravel(self.delta0)
   1095     self.work[:len(d0)] = d0

File ~***/lib/python3.11/site-packages/scipy/odr/_odrpack.py:901, in ODR._gen_work(self)
    897 else:
    898     # Okay, this isn't precisely right, but for this calculation,
    899     # it's fine
    900     ldwe = 1
--> 901     ld2we = self.data.we.shape[1]
    903 if self.job % 10 < 2:
    904     # ODR not OLS
    905     lwork = (18 + 11*p + p*p + m + m*m + 4*n*q + 6*n*m + 2*n*q*p +
    906              2*n*q*m + q*q + 5*q + q*(p+m) + ldwe*ld2we*q)

IndexError: tuple index out of range

SciPy/NumPy/Python version and system information

1.10.1 1.24.3 sys.version_info(major=3, minor=11, micro=3, releaselevel='final', serial=0)
Build Dependencies:
  blas:
    detection method: cmake
    found: true
    include directory: unknown
    lib directory: unknown
    name: OpenBLAS
    openblas configuration: unknown
    pc file directory: unknown
    version: 0.3.18
  lapack:
    detection method: cmake
    found: true
    include directory: unknown
    lib directory: unknown
    name: OpenBLAS
    openblas configuration: unknown
    pc file directory: unknown
    version: 0.3.18
Compilers:
  c:
    commands: cc
    linker: ld64
    name: clang
    version: 13.1.6
  c++:
    commands: c++
    linker: ld64
    name: clang
    version: 13.1.6
  cython:
    commands: cython
    linker: cython
    name: cython
    version: 0.29.33
  fortran:
    commands: gfortran
    linker: ld64
    name: gcc
    version: 12.1.0
  pythran:
    include directory: /private/var/folders/_f/lyvxf0v13gs7984d7sf7j83c0000gn/T/pip-build-env-c1_tgrdc/overlay/lib/python3.11/site-packages/pythran
    version: 0.12.1
Machine Information:
  build:
    cpu: aarch64
    endian: little
    family: aarch64
    system: darwin
  cross-compiled: false
  host:
    cpu: aarch64
    endian: little
    family: aarch64
    system: darwin
Python Information:
  path: /private/var/folders/_f/lyvxf0v13gs7984d7sf7j83c0000gn/T/cibw-run-s3_k_ke5/cp311-macosx_arm64/build/venv/bin/python
  version: '3.11'

[rkern]

You're right, I was probably thinking of it being fully-expanded at the time. The C code has the full logic handling all of the cases and should be replicated over here.

[jsdodge]

Thanks for the confirmation. I'm working on a project that needs to use that parameter, so I'd be interested in the guesstimated timescale for fixing it. If it is likely to be fixed soon, I'll work on other things. Otherwise, can you suggest any workarounds, other than calling scipy.odr.odr directly?

[rkern]

You can manually apply the fix in gh-18741 to your copy. It might be some time before the next release cycle.

[jsdodge]

Thanks for addressing this so quickly! I don't have experience with manual changes to a large external library like Scipy, so I f you could elaborate on how to do that or point me to any resources, I'd appreciate it.

[rkern]

It's just a few lines in scipy/odr/_odrpack.py. It's all pure Python, so there is no need to rebuild; you can edit the installed file in site-packages directly.

[jsdodge]

Thanks again. I can do that for now, but I'm aiming to distribute my code as a package, so it would be better to have a way to store the bug fix in my source directory until it appears in an official release.

[rkern]

You can always subclass ODR with the fixed copy of the _get_work() method.

[jsdodge]

Perfect. Thanks again!