Merge pull request #2 from HugoMVale/odr_fit

flip arg order in fcn signature

Author: HugoMVale

README.md

@@ -35,27 +35,26 @@ pip install odrpack
 The following example demonstrates a simple use of the package. For more comprehensive examples and explanations, please refer to the [documentation](https://hugomvale.github.io/odrpack-python/) pages.
 
 ```py
-from odrpack import odr
+from odrpack import odr_fit
 import numpy as np
 
-x = np.array([0.982, 1.998, 4.978, 6.01])
-y = np.array([2.7, 7.4, 148.0, 403.0])
+xdata = np.array([0.982, 1.998, 4.978, 6.01])
+ydata = np.array([2.7, 7.4, 148.0, 403.0])
 
 beta0 = np.array([2.0, 0.5])
-lower = np.array([0.0, 0.0])
-upper = np.array([10.0, 0.9])
+bounds = (np.array([0.0, 0.0]), np.array([10.0, 0.9]))
 
 def f(beta: np.ndarray, x: np.ndarray) -> np.ndarray:
     "Model function."
     return beta[0] * np.exp(beta[1]*x)
 
-sol = odr(f, beta0, y, x, lower=lower, upper=upper, iprint=1001)
+sol = odr_fit(f, xdata, ydata, beta0, bounds=bounds)
 
 print("beta:", sol.beta)
 print("delta:", sol.delta)
 ```
 
 ```sh
-beta: [1.63337057 0.9       ]
-delta: [-0.36885787 -0.31272733  0.02928942  0.11031791]
+beta: [1.63336897 0.9       ]
+delta: [-0.36885696 -0.31272648  0.02929022  0.11031872]
 ```

docs/examples/explicit-model.ipynb

@@ -88,7 +88,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "def f(beta: np.ndarray, X: np.ndarray) -> np.ndarray:\n",
+    "def f(X: np.ndarray, beta: np.ndarray) -> np.ndarray:\n",
     "    x0, y0, a, b, theta = beta\n",
     "    x, y = X\n",
     "    return ((x - x0)*np.cos(theta) + (y - y0)*np.sin(theta))**2 / a**2 \\\n",
@@ -140,7 +140,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 6,
    "metadata": {},
    "outputs": [],
    "source": [

docs/examples/implicit-model.ipynb

@@ -147,7 +147,7 @@ int odr_wrapper(int n,
         try {
             // Evaluate model function
             if (*ideval % 10 > 0) {
-                auto f_object = fcn_f_holder(beta_ndarray, xplusd_ndarray);
+                auto f_object = fcn_f_holder(xplusd_ndarray, beta_ndarray);
                 auto f_ndarray = nb::cast<nb::ndarray<const double, nb::c_contig>>(f_object);
                 auto f_ndarray_ptr = f_ndarray.data();
                 for (auto i = 0; i < (*q) * (*n); i++) {
@@ -157,7 +157,7 @@ int odr_wrapper(int n,
 
             // Model partial derivatives wrt `beta`
             if ((*ideval / 10) % 10 != 0) {
-                auto fjacb_object = fcn_fjacb_holder(beta_ndarray, xplusd_ndarray);
+                auto fjacb_object = fcn_fjacb_holder(xplusd_ndarray, beta_ndarray);
                 auto fjacb_ndarray = nb::cast<nb::ndarray<const double, nb::c_contig>>(fjacb_object);
                 auto fjacb_ndarray_ptr = fjacb_ndarray.data();
                 for (auto i = 0; i < (*q) * (*npar) * (*n); i++) {
@@ -167,7 +167,7 @@ int odr_wrapper(int n,
 
             // Model partial derivatives wrt `delta`
             if ((*ideval / 100) % 10 != 0) {
-                auto fjacd_object = fcn_fjacd_holder(beta_ndarray, xplusd_ndarray);
+                auto fjacd_object = fcn_fjacd_holder(xplusd_ndarray, beta_ndarray);
                 auto fjacd_ndarray = nb::cast<nb::ndarray<const double, nb::c_contig>>(fjacd_object);
                 auto fjacd_ndarray_ptr = fjacd_ndarray.data();
                 for (auto i = 0; i < (*q) * (*npar) * (*n); i++) {
@@ -261,13 +261,13 @@ ldstpd : int
 ldscld : int
     Leading dimension of the `scld` array, must be in `{1, n}`.
 f : Callable
-    User-supplied function for evaluating the model, `f(beta, x)`.
+    User-supplied function for evaluating the model, `f(x, beta)`.
 fjacb : Callable
     User-supplied function for evaluating the Jacobian w.r.t. `beta`,
-    `fjacb(beta, x)`.
+    `fjacb(x, beta)`.
 fjacd : Callable
     User-supplied function for evaluating the Jacobian w.r.t. `delta`,
-    `fjacd(beta, x)`.
+    `fjacd(x, beta)`.
 beta : np.ndarray[float64]
     Array of function parameters with shape `(npar)`.
 y : np.ndarray[float64]

src/odrpack/__odrpack.py.cpp

@@ -79,13 +79,13 @@ def odr(n: int, m: int, q: int, npar: int, ldwe: int, ld2we: int, ldwd: int, ld2
     ldscld : int
         Leading dimension of the `scld` array, must be in `{1, n}`.
     f : Callable
-        User-supplied function for evaluating the model, `f(beta, x)`.
+        User-supplied function for evaluating the model, `f(x, beta)`.
     fjacb : Callable
         User-supplied function for evaluating the Jacobian w.r.t. `beta`,
-        `fjacb(beta, x)`.
+        `fjacb(x, beta)`.
     fjacd : Callable
         User-supplied function for evaluating the Jacobian w.r.t. `delta`,
-        `fjacd(beta, x)`.
+        `fjacd(x, beta)`.
     beta : np.ndarray[float64]
         Array of function parameters with shape `(npar)`.
     y : np.ndarray[float64]

src/odrpack/__odrpack.pyi

@@ -253,7 +253,7 @@ def odr(f: Callable[[NDArray[np.float64], NDArray[np.float64]], NDArray[np.float
     ...     return beta[0] * np.exp(beta[1]*x)
     >>> sol = odr(f, beta0, y, x, lower=lower, upper=upper, iprint=0)
     >>> sol.beta
-    array([1.63337057, 0.9       ])
+    array([1.63336897, 0.9       ])
     """
 
     # Future deprecation warning
@@ -434,29 +434,33 @@ def odr(f: Callable[[NDArray[np.float64], NDArray[np.float64]], NDArray[np.float
 
     # Check model function and jacobians
     f0 = f(beta0, x)
+    def f_(x, beta): return f(beta, x)
+
     if f0.shape != y.shape:
         raise ValueError(
             "Function `f` must return an array with the same shape as `y`.")
 
-    def fdummy(beta, x): return np.array([np.nan])  # will never be called
+    def fdummy(x, beta): return np.array([np.nan])  # will never be called
 
     if has_jac and fjacb is not None:
         fjacb0 = fjacb(beta0, x)
+        def fjacb_(x, beta): return fjacb(beta, x)
         if fjacb0.shape[-1] != n or fjacb0.size != n*npar*q:
             raise ValueError(
                 "Function `fjacb` must return an array with shape `(n, npar, q)` or compatible.")
     elif not has_jac and fjacb is None:
-        fjacb = fdummy
+        fjacb_ = fdummy
     else:
         raise ValueError("Inconsistent arguments for `job` and `fjacb`.")
 
     if has_jac and fjacd is not None:
         fjacd0 = fjacd(beta0, x)
+        def fjacd_(x, beta): return fjacd(beta, x)
         if fjacd0.shape[-1] != n or fjacd0.size != n*m*q:
             raise ValueError(
                 "Function `fjacd` must return an array with shape `(n, m, q)` or compatible.")
     elif not has_jac and fjacd is None:
-        fjacd = fdummy
+        fjacd_ = fdummy
     else:
         raise ValueError("Inconsistent arguments for `job` and `fjacd`.")
 
@@ -483,7 +487,7 @@ def fdummy(beta, x): return np.array([np.nan])  # will never be called
                 ldwd=ldwd, ld2wd=ld2wd,
                 ldifx=ldifx,
                 ldstpd=ldstpd, ldscld=ldscld,
-                f=f, fjacb=fjacb, fjacd=fjacd,
+                f=f_, fjacb=fjacb_, fjacd=fjacd_,
                 beta=beta, y=y, x=x,
                 delta=delta,
                 we=we, wd=wd, ifixb=ifixb, ifixx=ifixx,

src/odrpack/odr_fortran.py

@@ -48,7 +48,7 @@ def odr_fit(f: Callable[[NDArray[np.float64], NDArray[np.float64]], NDArray[np.f
     Parameters
     ----------
     f : Callable[[NDArray[np.float64], NDArray[np.float64]], NDArray[np.float64]]
-        Function to be fitted, with the signature `f(beta, x)`. It must return
+        Function to be fitted, with the signature `f(x, beta)`. It must return
         an array with the same shape as `y`.
     xdata : NDArray[np.float64]
         Array of shape `(n,)` or `(m, n)` containing the observed values of the
@@ -117,13 +117,13 @@ def odr_fit(f: Callable[[NDArray[np.float64], NDArray[np.float64]], NDArray[np.f
         all `xdata` values are automatically treated as fixed.
     jac_beta : Callable[[NDArray[np.float64], NDArray[np.float64]], NDArray[np.float64]] | None
         Jacobian of the function to be fitted with respect to `beta`, with the
-        signature `jac_beta(beta, x)`. It must return an array with shape 
+        signature `jac_beta(x, beta)`. It must return an array with shape 
         `(n, npar, q)` or a compatible shape. By default, the Jacobian is
         approximated numerically using the finite difference scheme specified
         by `diff_scheme`.
     jac_x : Callable[[NDArray[np.float64], NDArray[np.float64]], NDArray[np.float64]] | None
         Jacobian of the function to be fitted with respect to `x`, with the
-        signature `jac_x(beta, x)`. It must return an array with shape 
+        signature `jac_x(x, beta)`. It must return an array with shape 
         `(n, m, q)` or a compatible shape. By default, the Jacobian is approximated
         numerically using the finite difference scheme specified by `diff_scheme`.
     delta0 : NDArray[np.float64] | None
@@ -229,7 +229,7 @@ def odr_fit(f: Callable[[NDArray[np.float64], NDArray[np.float64]], NDArray[np.f
     >>> ydata = np.array([2.7, 7.4, 148.0, 403.0])
     >>> beta0 = np.array([2., 0.5])
     >>> bounds = (np.array([0., 0.]), np.array([10., 0.9]))
-    >>> def f(beta: np.ndarray, x: np.ndarray) -> np.ndarray:
+    >>> def f(x: np.ndarray, beta: np.ndarray) -> np.ndarray:
     ...     return beta[0] * np.exp(beta[1]*x)
     >>> sol = odr_fit(f, xdata, ydata, beta0, bounds=bounds)
     >>> sol.beta
@@ -409,25 +409,25 @@ def odr_fit(f: Callable[[NDArray[np.float64], NDArray[np.float64]], NDArray[np.f
         ld2we = 1
 
     # Check model function
-    f0 = f(beta0, xdata)
+    f0 = f(xdata, beta0)
     if f0.shape != ydata.shape:
         raise ValueError(
             "Function `f` must return an array with the same shape as `ydata`.")
 
     # Check model jacobians
-    def fdummy(beta, x): return np.array([np.nan])  # will never be called
+    def fdummy(x, beta): return np.array([np.nan])  # will never be called
 
     if jac_beta is None and jac_x is None:
         has_jac = False
         jac_beta = fdummy
         jac_x = fdummy
     elif jac_beta is not None and jac_x is not None:
         has_jac = True
-        jac0_beta = jac_beta(beta0, xdata)
+        jac0_beta = jac_beta(xdata,  beta0)
         if jac0_beta.shape[-1] != n or jac0_beta.size != n*npar*q:
             raise ValueError(
                 "Function `jac_beta` must return an array with shape `(n, npar, q)` or compatible.")
-        jac0_x = jac_x(beta0, xdata)
+        jac0_x = jac_x(xdata, beta0)
         if jac0_x.shape[-1] != n or jac0_x.size != n*m*q:
             raise ValueError(
                 "Function `jac_x` must return an array with shape `(n, m, q)` or compatible.")

src/odrpack/odr_scipy.py

@@ -42,19 +42,19 @@ def test_dimension_consistency():
 def test_odr():
     "example5 from odrpack"
 
-    def f(beta: np.ndarray, x: np.ndarray) -> np.ndarray:
+    def f(x: np.ndarray, beta: np.ndarray) -> np.ndarray:
         return beta[0] * np.exp(beta[1]*x)
 
-    def fjacb(beta: np.ndarray, x: np.ndarray) -> np.ndarray:
+    def fjacb(x: np.ndarray, beta: np.ndarray) -> np.ndarray:
         jac = np.zeros((beta.size, x.size))
         jac[0, :] = np.exp(beta[1]*x)
         jac[1, :] = beta[0]*x*np.exp(beta[1]*x)
         return jac
 
-    def fjacd(beta: np.ndarray, x: np.ndarray) -> np.ndarray:
+    def fjacd(x: np.ndarray, beta: np.ndarray) -> np.ndarray:
         return beta[0] * beta[1] * np.exp(beta[1]*x)
 
-    def fdummy(beta: np.ndarray, x: np.ndarray) -> np.ndarray:
+    def fdummy(x: np.ndarray, beta: np.ndarray) -> np.ndarray:
         return np.array([42.0])
 
     beta0 = np.array([2., 0.5])
@@ -77,7 +77,7 @@ def fdummy(beta: np.ndarray, x: np.ndarray) -> np.ndarray:
                f, fdummy, fdummy, beta, y, x, delta,
                lower=lower, upper=upper, job=0)
     assert info == 1
-    np.allclose(beta, beta_ref)
+    assert np.allclose(beta, beta_ref)
 
     # solution with jacobian
     beta = beta0.copy()

tests/test_bindings.py

@@ -10,27 +10,27 @@
 # %% Functions need to be defined outside the test function
 
 
-def f1(beta: np.ndarray, x: np.ndarray) -> np.ndarray:
+def f1(x: np.ndarray, beta: np.ndarray) -> np.ndarray:
     return beta[0] + beta[1] * x + beta[2] * x**2 + beta[3] * x**3
 
 
 beta1 = np.array([1, -2., 0.1, -0.1])
 x1 = np.linspace(-10., 10., 21, dtype=np.float64)
-y1 = f1(beta1, x1)
+y1 = f1(x1, beta1)
 
 
-def f2(beta: np.ndarray, x: np.ndarray) -> np.ndarray:
+def f2(x: np.ndarray, beta: np.ndarray) -> np.ndarray:
     time.sleep(np.random.uniform(0, DELAY))
     return (beta[0] * x[0, :])**3 + x[1, :]**beta[1]
 
 
 beta2 = np.array([2., 2.])
 x2 = np.linspace(-10., 10., 41, dtype=np.float64)
 x2 = np.vstack((x2, 10+x2/2))
-y2 = f2(beta2, x2)
+y2 = f2(x2, beta2)
 
 
-def f3(beta: np.ndarray, x: np.ndarray) -> np.ndarray:
+def f3(x: np.ndarray, beta: np.ndarray) -> np.ndarray:
     time.sleep(np.random.uniform(0, DELAY))
     y = np.zeros((2, x.shape[-1]))
     y[0, :] = (beta[0] * x[0, :])**3 + x[1, :]**beta[1] + np.exp(x[2, :]/2)
@@ -41,7 +41,7 @@ def f3(beta: np.ndarray, x: np.ndarray) -> np.ndarray:
 beta3 = np.array([1., 2., 3.])
 x3 = np.linspace(-1., 1., 31, dtype=np.float64)
 x3 = np.vstack((x3, np.exp(x3), x3**2))
-y3 = f3(beta3, x3)
+y3 = f3(x3, beta3)
 
 case1 = (f1, x1, y1, np.ones_like(beta1))
 case2 = (f2, x2, y2, np.ones_like(beta2))

tests/test_multiprocessing.py

@@ -19,12 +19,12 @@ def add_noise(array, noise, seed):
 @pytest.fixture
 def case1():
     # m=1, q=1
-    def f(beta: np.ndarray, x: np.ndarray) -> np.ndarray:
+    def f(x: np.ndarray, beta: np.ndarray) -> np.ndarray:
         return beta[0] + beta[1] * x + beta[2] * x**2 + beta[3] * x**3
 
     beta_star = np.array([1, -2., 0.1, -0.1])
     x = np.linspace(-10., 10., 21)
-    y = f(beta_star, x)
+    y = f(x, beta_star)
 
     x = add_noise(x, 5e-2, SEED)
     y = add_noise(y, 10e-2, SEED)
@@ -35,13 +35,13 @@ def f(beta: np.ndarray, x: np.ndarray) -> np.ndarray:
 @pytest.fixture
 def case2():
     # m=2, q=1
-    def f(beta: np.ndarray, x: np.ndarray) -> np.ndarray:
+    def f(x: np.ndarray, beta: np.ndarray) -> np.ndarray:
         return (beta[0] * x[0, :])**3 + x[1, :]**beta[1]
 
     beta_star = np.array([2., 2.])
     x1 = np.linspace(-10., 10., 41)
     x = np.vstack((x1, 10+x1/2))
-    y = f(beta_star, x)
+    y = f(x, beta_star)
 
     x = add_noise(x, 5e-2, SEED)
     y = add_noise(y, 10e-2, SEED)
@@ -52,7 +52,7 @@ def f(beta: np.ndarray, x: np.ndarray) -> np.ndarray:
 @pytest.fixture
 def case3():
     # m=3, q=2
-    def f(beta: np.ndarray, x: np.ndarray) -> np.ndarray:
+    def f(x: np.ndarray, beta: np.ndarray) -> np.ndarray:
         y = np.zeros((2, x.shape[-1]))
         y[0, :] = (beta[0] * x[0, :])**3 + x[1, :]**beta[1] + np.exp(x[2, :]/2)
         y[1, :] = (beta[2] * x[0, :])**2 + x[1, :]**beta[1]
@@ -61,7 +61,7 @@ def f(beta: np.ndarray, x: np.ndarray) -> np.ndarray:
     beta_star = np.array([1., 2., 3.])
     x1 = np.linspace(-1., 1., 31)
     x = np.vstack((x1, np.exp(x1), x1**2))
-    y = f(beta_star, x)
+    y = f(x, beta_star)
 
     x = add_noise(x, 5e-2, SEED)
     y = add_noise(y, 10e-2, SEED)
@@ -460,16 +460,16 @@ def test_jacobians():
     beta0 = np.array([2., 0.5])
     bounds = (np.array([0., 0.]), np.array([10., 0.9]))
 
-    def f(beta: np.ndarray, x: np.ndarray) -> np.ndarray:
+    def f(x: np.ndarray, beta: np.ndarray) -> np.ndarray:
         return beta[0] * np.exp(beta[1]*x)
 
-    def jac_beta(beta: np.ndarray, x: np.ndarray) -> np.ndarray:
+    def jac_beta(x: np.ndarray, beta: np.ndarray) -> np.ndarray:
         jac = np.zeros((beta.size, x.size))
         jac[0, :] = np.exp(beta[1]*x)
         jac[1, :] = beta[0]*x*np.exp(beta[1]*x)
         return jac
 
-    def jac_x(beta: np.ndarray, x: np.ndarray) -> np.ndarray:
+    def jac_x(x: np.ndarray, beta: np.ndarray) -> np.ndarray:
         return beta[0] * beta[1] * np.exp(beta[1]*x)
 
     beta_ref = np.array([1.63337602, 0.9])
@@ -539,7 +539,7 @@ def test_implicit_model():
     xdata = np.array(x).T
     ydata = np.full(xdata.shape[-1], 0.0)
 
-    def f(beta: np.ndarray, x: np.ndarray) -> np.ndarray:
+    def f(x: np.ndarray, beta: np.ndarray) -> np.ndarray:
         v, h = x
         return beta[2]*(v-beta[0])**2 + 2*beta[3]*(v-beta[0])*(h-beta[1]) \
             + beta[4]*(h-beta[1])**2 - 1