Broadcast curve_fit across other dimensions

[TheFermiSea]

Hello,

I have some incredibly large datasets over which I need to apply curve fitting. I am able to use scipp.scipy.optimize.curve_fit to perform curve fitting on a single 1D DataArray, but need to iterate this over 2048*2048 pixels (eventually other dimensions as well). Sure, I can naively dump this into a for loop, but this kinda defeats the purpose of using Scipp in the first place. My current solution is to use scipp.compat.to_xarray, perform the curve fitting in xarray, then convert back. In doing so, I lose units, strong typing, and variances, which is managable but unfortunate. I am wondering if there are plans to implement broadcasting of curve_fit? This is not a deal breaker for me, I still love everything about Scipp and will continue to use it for most things, but it would be incredible to be able to do this all in Scipp.

[SimonHeybrock]

Can you clarify which of these you are trying to achieve?

1. A single 1-D fit (returning a single set of fit parameters) with extra dimensions that do not participate in the fit (I think scipy.curve_fit can do this, see the note on (k,M)-shaped inputs, I suppose k could label the pixel).
2. A single N-D fit, i.e., what xarray.curvefit seems to be able to do by passing multiple coords.
3. A different 1-D fit (returning different fit parameters) for every pixel, i.e., you would like to map a 1-D fit over all pixels?

[TheFermiSea]

In this specific case, I am trying to do 3), a map of individual, uncoupled 1D fits, that returns a DataArray with post and pcov for every element in the other dimensions (pixels in this case). This is trivial in xarray, but I don't know how that translates to the Scipp backend.

[SimonHeybrock]

I think there are two options: 1. enhance our SciPy wrapper, 2. provide a wrapper for xarray's curvefit. We will look into providing a solution for this, as this seems generally useful!

Could you provide an example of how you to this with xarray? Is it as simple as providing a single coord to xarray.DataArray.curvefit?

[TheFermiSea]

Yes, in xarray the curvefit method uses apply_ufunc to wrap scipy.curve_fit. For a single independent variable, the coordinate is passed as the first argument. This seems to work for more than one independent variable (coordinate), but I could see how that would be much more complicated to implement.

[nvaytet]

Hi there, @SimonHeybrock has just gone on parental leave, so I will pick up this thread. As he suggested, could you provide a MWE in xarray, and I can then look into how we could make an equivalent using Scipp? Thanks!

[SimonHeybrock]

We have just released Scipp 23.12.0, which adds scipp.curve_fit. This works similar to xarray.DataArray.curvefit. @TheFermiSea please give it a try if you can, any feedback (including potential UX or performance issues) is welcome!