Repository: scipy/scipy
Issue category: Bug Report

Issue:
BUG: `optimize.curve_fit` with `nan_policy="omit"` fails to handle sigma properly
### Describe your issue.

Hello,
I have encountered an issue with curve_fit. It is bypassable, but I'd say it is unexpected behaviour. Generally speaking, it occurs when both of the following is true:
-- I have xdata or ydata with NaNs and use nan_policy = "omit" to handle them;
-- I have specified sigma.
In this case, I get ValueError "Sigma has incorrect shape" in optimize._minpack_py.py despite sigma being of the same size as xdata and ydata.

I checked what is going in _minpack_py.py and figured out the following.
It occurs because "nan_policy = omit" leads to the removal of NaN indices from xdata and ydata (lines 960-967), but not from sigma. Then, after line 969, the shape of sigma is checked by comparing it with ydata. And even if sigma and ydata were originally of the same size, they end up to be different after ydata truncation, throwing an error.

So, to bypass this problem, one has either to remove indices corresponding NaNs in xdata and ydata from sigma, or not using nan_policy = "omit" at all and remove NaN indices manually. I suggest improving nan_policy = "omit" in a way that it would remove the corresponding indices from sigma, if needed.

Sorry if my bug report is incorrect, I am first time doing this.

### Reproducing Code Example

```python
import numpy as np
from scipy.optimize import curve_fit

def exponential(x, a, b):
	return b * np.exp(a * x)

np.random.seed(0)
N = 90
x = np.linspace(1, 10, N, endpoint = False)
sigma = x * 0.05
y = exponential(x, 0.2, 0.5) + np.random.normal(0, sigma, N)

p0 = [0.1, 1.0]

# without nans, works perfect
popt, pcov = curve_fit(exponential, x, y, p0 = p0, sigma = sigma)
perr = np.diag(np.sqrt(pcov))
print("without nans:")
print(f"popt: {popt}")
print(f"perr: {perr}")

# with nans, have to set nan_policy = "omit"...
x[11] = np.nan
y[22] = np.nan
popt, pcov = curve_fit(exponential, x, y, p0 = p0, nan_policy = "omit")
perr = np.diag(np.sqrt(pcov))
print("with nans:")
print(f"popt: {popt}")
print(f"perr: {perr}")

#...but when using sigma, it says "sigma has incorrect shape"
x[11] = np.nan
y[22] = np.nan
popt, pcov = curve_fit(exponential, x, y, p0 = p0, sigma = sigma, nan_policy = "omit")
```


### Error message

```shell
Traceback (most recent call last):
  File "/root/Science/Spiral_shape_analysis/fit_results/../spiralcutter_advanced/scipy_test.py", line 34, in <module>
    popt, pcov = curve_fit(exponential, x, y, p0 = p0, sigma = sigma, nan_policy = "omit")
  File "/usr/local/lib/python3.10/dist-packages/scipy/optimize/_minpack_py.py", line 985, in curve_fit
    raise ValueError("`sigma` has incorrect shape.")
ValueError: `sigma` has incorrect shape.
```


### SciPy/NumPy/Python version and system information

```shell
>>> print(scipy.__version__, numpy.__version__, sys.version_info)
1.14.1 2.1.2 sys.version_info(major=3, minor=10, micro=12, releaselevel='final', serial=0)
>>> scipy.show_config()
Build Dependencies:
  blas:
    detection method: pkgconfig
    found: true
    include directory: /opt/_internal/cpython-3.10.14/lib/python3.10/site-packages/scipy_openblas32/include
    lib directory: /opt/_internal/cpython-3.10.14/lib/python3.10/site-packages/scipy_openblas32/lib
    name: scipy-openblas
    openblas configuration: OpenBLAS 0.3.27.dev DYNAMIC_ARCH NO_AFFINITY Zen MAX_THREADS=64
    pc file directory: /project
    version: 0.3.27.dev
  lapack:
    detection method: pkgconfig
    found: true
    include directory: /opt/_internal/cpython-3.10.14/lib/python3.10/site-packages/scipy_openblas32/include
    lib directory: /opt/_internal/cpython-3.10.14/lib/python3.10/site-packages/scipy_openblas32/lib
    name: scipy-openblas
    openblas configuration: OpenBLAS 0.3.27.dev DYNAMIC_ARCH NO_AFFINITY Zen MAX_THREADS=64
    pc file directory: /project
    version: 0.3.27.dev
  pybind11:
    detection method: config-tool
    include directory: unknown
    name: pybind11
    version: 2.12.0
Compilers:
  c:
    commands: cc
    linker: ld.bfd
    name: gcc
    version: 10.2.1
  c++:
    commands: c++
    linker: ld.bfd
    name: gcc
    version: 10.2.1
  cython:
    commands: cython
    linker: cython
    name: cython
    version: 3.0.11
  fortran:
    commands: gfortran
    linker: ld.bfd
    name: gcc
    version: 10.2.1
  pythran:
    include directory: ../../tmp/pip-build-env-h_xz8lfs/overlay/lib/python3.10/site-packages/pythran
    version: 0.16.1
Machine Information:
  build:
    cpu: x86_64
    endian: little
    family: x86_64
    system: linux
  cross-compiled: false
  host:
    cpu: x86_64
    endian: little
    family: x86_64
    system: linux
Python Information:
  path: /opt/python/cp310-cp310/bin/python
  version: '3.10'
```



Hints:
Oops, thanks for reporting!

@AtsushiSakai Looks like we missed this. Do you have a moment to submit a PR?

@Vallastro Just to confirm, this would also be a problem when `sigma` is 2D. In that case, we would remove both the rows and columns where there are NaNs in either `x` or `y`?
> Just to confirm, this would also be a problem when `sigma` is 2D. In that case, we would remove both the rows and columns where there are NaNs in either `x` or `y`?

I guess, yes. In this case, if `ydata[i]` is invalidated, then for any `j` covariation between `xdata[i]` and `xdata[j]`, represented in `sigma[i, j]` and `sigma[j, i]` has no sense. And if `N` elements of `ydata` out of original `M` are remaining after removing NaNs, then removing columns and rows from `sigma` will make it `N` x `N` from `M` x `M`, which is expected.
@mdhaber Is this issue accepting PRs? If so, could you please point me to the files that need to be changed? I would love to contribute.

Patch:
diff --git a/scipy/optimize/_minpack_py.py b/scipy/optimize/_minpack_py.py
index 8c0564c06bbe..8cfce8aae21b 100644
--- a/scipy/optimize/_minpack_py.py
+++ b/scipy/optimize/_minpack_py.py
@@ -964,12 +964,21 @@ def curve_fit(f, xdata, ydata, p0=None, sigma=None, absolute_sigma=False,
             xdata = xdata[..., ~has_nan]
             ydata = ydata[~has_nan]
 
+            # Also omit the corresponding entries from sigma
+            if sigma is not None:
+                sigma = np.asarray(sigma)
+                if sigma.ndim == 1:
+                    sigma = sigma[~has_nan]
+                elif sigma.ndim == 2:
+                    sigma = sigma[~has_nan, :]
+                    sigma = sigma[:, ~has_nan]
+
     # Determine type of sigma
     if sigma is not None:
         sigma = np.asarray(sigma)
 
         # if 1-D or a scalar, sigma are errors, define transform = 1/sigma
-        if sigma.size == 1 or sigma.shape == (ydata.size, ):
+        if sigma.size == 1 or sigma.shape == (ydata.size,):
             transform = 1.0 / sigma
         # if 2-D, sigma is the covariance matrix,
         # define transform = L such that L L^T = C


Test patch:
diff --git a/scipy/optimize/tests/test_minpack.py b/scipy/optimize/tests/test_minpack.py
index dce5776d27cb..ef107c5692c2 100644
--- a/scipy/optimize/tests/test_minpack.py
+++ b/scipy/optimize/tests/test_minpack.py
@@ -811,6 +811,61 @@ def test_maxfev_and_bounds(self):
         assert_allclose(popt1, 2, atol=1e-14)
         assert_allclose(popt2, 2, atol=1e-14)
 
+    @pytest.mark.parametrize("sigma_dim", [0, 1, 2])
+    def test_curvefit_omitnan(self, sigma_dim):
+        def exponential(x, a, b):
+            return b * np.exp(a * x)
+
+        rng = np.random.default_rng(578285731148908)
+        N = 100
+        x = np.linspace(1, 10, N)
+        y = exponential(x, 0.2, 0.5)
+
+        if (sigma_dim == 0):
+            sigma = 0.05
+            y += rng.normal(0, sigma, N)
+
+        elif (sigma_dim == 1):
+            sigma = x * 0.05
+            y += rng.normal(0, sigma, N)
+
+        elif (sigma_dim == 2):
+            # The covariance matrix must be symmetric positive-semidefinite
+            a = rng.normal(0, 2, (N, N))
+            sigma = a @ a.T
+            y += rng.multivariate_normal(np.zeros_like(x), sigma)
+        else:
+            assert False, "The sigma must be a scalar, 1D array or 2D array."
+
+        p0 = [0.1, 1.0]
+
+        # Choose indices to place NaNs.
+        i_x = rng.integers(N, size=5)
+        i_y = rng.integers(N, size=5)
+
+        # Add NaNs and compute result using `curve_fit`
+        x[i_x] = np.nan
+        y[i_y] = np.nan
+        res_opt, res_cov = curve_fit(exponential, x, y, p0=p0, sigma=sigma,
+                                     nan_policy="omit")
+
+        # Manually remove elements that should be eliminated, and
+        # calculate reference using `curve_fit`
+        i_delete = np.unique(np.concatenate((i_x, i_y)))
+        x = np.delete(x, i_delete, axis=0)
+        y = np.delete(y, i_delete, axis=0)
+
+        sigma = np.asarray(sigma)
+        if sigma.ndim == 1:
+            sigma = np.delete(sigma, i_delete)
+        elif sigma.ndim == 2:
+            sigma = np.delete(sigma, i_delete, axis=0)
+            sigma = np.delete(sigma, i_delete, axis=1)
+        ref_opt, ref_cov = curve_fit(exponential, x, y, p0=p0, sigma=sigma)
+
+        assert_allclose(res_opt, ref_opt, atol=1e-14)
+        assert_allclose(res_cov, ref_cov, atol=1e-14)
+
     def test_curvefit_simplecovariance(self):
 
         def func(x, a, b):