Repository: pydata/xarray

Issue:
Handle scipy fitting errors in `xarray.curve_fit`
### Is your feature request related to a problem?

I'm trying to fit a bunch of exponential decays to a dataset of dimensions `(Time, experimental_index)`, along the `Time` dimension. This works nicely... except in some cases, `scipy.optimize` raises an error such as `RuntimeError: Optimal parameters not found: The maximum number of function evaluations is exceeded.`. That's, of course, reasonable, except in my use case I'm fine with a bunch of `experimental_index` failing to fit, and I'd still like to access the fits that made it through.

### Describe the solution you'd like

I'm proposing adding a flag for `xarray.curvefit`, for example `allow_failures: bool = False`, that would take this part:

https://github.com/pydata/xarray/blob/f9037c41e36254bfe7efa9feaedd3eae5512bd11/xarray/core/dataset.py#L7738-L7752

and, if active, wrap `popt, pcov = curve_fit(func, x, y, **kwargs)` in a `try/except` block that would output `nan` just as the `not len(y)` conditional clause does if `scipy.optimize.curve_fit` fails with a `RuntimeError`.

### Describe alternatives you've considered

At the moment, I sidestep the issue like this:

```python
fits = []
for experimental_index in timetraces.experimental_index:
    trace = traces.sel(experimental_index=experimental_index)
                                                           
    try:
        fit = trace.curvefit(
            "Time",
            fitfunc,
            p0=(...),
            bounds={...},
        )
    except (RuntimeError, TypeError) as e:
        logger.debug(...)
        continue

    fits.append(fit)

fits = xarray.concat(fits, dim="experimental_index", coords="all")
```
but that seems like a pretty ugly hack...	

### Additional context

I would of course be happy to implement this!

Discussion:
mgunyho: Oh no, the doctest failure is because the test is flaky, this was introduced in #7821, see here: https://github.com/pydata/xarray/pull/7821#issuecomment-1537142237 and here: https://github.com/pydata/xarray/pull/7821/commits/a0e6659ca01188378f29a35b418d6f9e2b889d2e. I'll submit another patch to fix it soon.
Illviljan: You could use `DataArray.round` to round to significant decimals.
headtr1ck: > You could use `DataArray.round` to round to significant decimals.

Better to use a tolerance in the assterion testing.

Patch:
diff --git a/doc/whats-new.rst b/doc/whats-new.rst
index 5c0d3c3c..420b6c55 100644
--- a/doc/whats-new.rst
+++ b/doc/whats-new.rst
@@ -25,6 +25,10 @@ New Features
 
 - Added support for multidimensional initial guess and bounds in :py:meth:`DataArray.curvefit` (:issue:`7768`, :pull:`7821`).
   By `András Gunyhó <https://github.com/mgunyho>`_.
+- Add an ``errors`` option to :py:meth:`Dataset.curve_fit` that allows
+  returning NaN for the parameters and covariances of failed fits, rather than
+  failing the whole series of fits (:issue:`6317`, :pull:`7891`).
+  By `Dominik Stańczak <https://github.com/StanczakDominik>`_ and `András Gunyhó <https://github.com/mgunyho>`_.
 
 Breaking changes
 ~~~~~~~~~~~~~~~~
@@ -46,6 +50,8 @@ Bug fixes
   By `Mattia Almansi <https://github.com/malmans2>`_.
 - Don't call ``CachingFileManager.__del__`` on interpreter shutdown (:issue:`7814`, :pull:`7880`).
   By `Justus Magin <https://github.com/keewis>`_.
+- Ensure dtype of reindex result matches dtype of the original DataArray (:issue:`7299`, :pull:`7917`)
+  By `Anderson Banihirwe <https://github.com/andersy005>`_.
 
 Documentation
 ~~~~~~~~~~~~~
diff --git a/xarray/core/dataarray.py b/xarray/core/dataarray.py
index bb6a9b13..9635d678 100644
--- a/xarray/core/dataarray.py
+++ b/xarray/core/dataarray.py
@@ -6162,6 +6162,7 @@ class DataArray(
         p0: dict[str, float | DataArray] | None = None,
         bounds: dict[str, tuple[float | DataArray, float | DataArray]] | None = None,
         param_names: Sequence[str] | None = None,
+        errors: ErrorOptions = "raise",
         kwargs: dict[str, Any] | None = None,
     ) -> Dataset:
         """
@@ -6206,6 +6207,10 @@ class DataArray(
             this will be automatically determined by arguments of `func`. `param_names`
             should be manually supplied when fitting a function that takes a variable
             number of parameters.
+        errors : {"raise", "ignore"}, default: "raise"
+            If 'raise', any errors from the `scipy.optimize_curve_fit` optimization will
+            raise an exception. If 'ignore', the coefficients and covariances for the
+            coordinates where the fitting failed will be NaN.
         **kwargs : optional
             Additional keyword arguments to passed to scipy curve_fit.
 
@@ -6312,6 +6317,7 @@ class DataArray(
             p0=p0,
             bounds=bounds,
             param_names=param_names,
+            errors=errors,
             kwargs=kwargs,
         )
 
diff --git a/xarray/core/dataset.py b/xarray/core/dataset.py
index 0f271e9d..f4ba9d4f 100644
--- a/xarray/core/dataset.py
+++ b/xarray/core/dataset.py
@@ -8631,6 +8631,7 @@ class Dataset(
         p0: dict[str, float | DataArray] | None = None,
         bounds: dict[str, tuple[float | DataArray, float | DataArray]] | None = None,
         param_names: Sequence[str] | None = None,
+        errors: ErrorOptions = "raise",
         kwargs: dict[str, Any] | None = None,
     ) -> T_Dataset:
         """
@@ -8675,6 +8676,10 @@ class Dataset(
             this will be automatically determined by arguments of `func`. `param_names`
             should be manually supplied when fitting a function that takes a variable
             number of parameters.
+        errors : {"raise", "ignore"}, default: "raise"
+            If 'raise', any errors from the `scipy.optimize_curve_fit` optimization will
+            raise an exception. If 'ignore', the coefficients and covariances for the
+            coordinates where the fitting failed will be NaN.
         **kwargs : optional
             Additional keyword arguments to passed to scipy curve_fit.
 
@@ -8757,6 +8762,9 @@ class Dataset(
                             f"dimensions {preserved_dims}."
                         )
 
+        if errors not in ["raise", "ignore"]:
+            raise ValueError('errors must be either "raise" or "ignore"')
+
         # Broadcast all coords with each other
         coords_ = broadcast(*coords_)
         coords_ = [
@@ -8793,7 +8801,15 @@ class Dataset(
                     pcov = np.full([n_params, n_params], np.nan)
                     return popt, pcov
             x = np.squeeze(x)
-            popt, pcov = curve_fit(func, x, y, p0=p0_, bounds=(lb, ub), **kwargs)
+
+            try:
+                popt, pcov = curve_fit(func, x, y, p0=p0_, bounds=(lb, ub), **kwargs)
+            except RuntimeError:
+                if errors == "raise":
+                    raise
+                popt = np.full([n_params], np.nan)
+                pcov = np.full([n_params, n_params], np.nan)
+
             return popt, pcov
 
         result = type(self)()
diff --git a/xarray/core/dtypes.py b/xarray/core/dtypes.py
index f5a18428..7ac342e3 100644
--- a/xarray/core/dtypes.py
+++ b/xarray/core/dtypes.py
@@ -74,7 +74,10 @@ def maybe_promote(dtype):
     else:
         dtype = object
         fill_value = np.nan
-    return np.dtype(dtype), fill_value
+
+    dtype = np.dtype(dtype)
+    fill_value = dtype.type(fill_value)
+    return dtype, fill_value
 
 
 NAT_TYPES = {np.datetime64("NaT").dtype, np.timedelta64("NaT").dtype}


Test patch:
diff --git a/xarray/tests/test_dataarray.py b/xarray/tests/test_dataarray.py
index e3138a32..cee5afa5 100644
--- a/xarray/tests/test_dataarray.py
+++ b/xarray/tests/test_dataarray.py
@@ -1698,6 +1698,19 @@ class TestDataArray:
         assert_identical(expected, actual)
         assert actual.dtype == expected.dtype
 
+    def test_reindex_empty_array_dtype(self) -> None:
+        # Dtype of reindex result should match dtype of the original DataArray.
+        # See GH issue #7299
+        x = xr.DataArray([], dims=("x",), coords={"x": []}).astype("float32")
+        y = x.reindex(x=[1.0, 2.0])
+
+        assert (
+            x.dtype == y.dtype
+        ), "Dtype of reindexed DataArray should match dtype of the original DataArray"
+        assert (
+            y.dtype == np.float32
+        ), "Dtype of reindexed DataArray should remain float32"
+
     def test_rename(self) -> None:
         da = xr.DataArray(
             [1, 2, 3], dims="dim", name="name", coords={"coord": ("dim", [5, 6, 7])}
@@ -4571,6 +4584,48 @@ class TestDataArray:
                 bounds={"a": (0, DataArray([1], coords={"foo": [1]}))},
             )
 
+    @requires_scipy
+    @pytest.mark.parametrize("use_dask", [True, False])
+    def test_curvefit_ignore_errors(self, use_dask: bool) -> None:
+        if use_dask and not has_dask:
+            pytest.skip("requires dask")
+
+        # nonsense function to make the optimization fail
+        def line(x, a, b):
+            if a > 10:
+                return 0
+            return a * x + b
+
+        da = DataArray(
+            [[1, 3, 5], [0, 20, 40]],
+            coords={"i": [1, 2], "x": [0.0, 1.0, 2.0]},
+        )
+
+        if use_dask:
+            da = da.chunk({"i": 1})
+
+        expected = DataArray(
+            [[2, 1], [np.nan, np.nan]], coords={"i": [1, 2], "param": ["a", "b"]}
+        )
+
+        with pytest.raises(RuntimeError, match="calls to function has reached maxfev"):
+            da.curvefit(
+                coords="x",
+                func=line,
+                # limit maximum number of calls so the optimization fails
+                kwargs=dict(maxfev=5),
+            ).compute()  # have to compute to raise the error
+
+        fit = da.curvefit(
+            coords="x",
+            func=line,
+            errors="ignore",
+            # limit maximum number of calls so the optimization fails
+            kwargs=dict(maxfev=5),
+        ).compute()
+
+        assert_allclose(fit.curvefit_coefficients, expected)
+
 
 class TestReduce:
     @pytest.fixture(autouse=True)