update docstrings on `analyze.py` and add some disclaimer about the usage of `5pt_stencil`

[fjaviersanchez]

I thought we already had done this, but it might have gotten mixed up in some update. We should try to update the docstrings in analyze.py. In particular the get_derivatives method needs to be updated to include a mention of derivkit as an option.

As a stretch goal for this issue we should review the rest of the docstrings and see where updates or streamlining is needed and remove spurious / stale comments.

Finally, as suggested by @nikosarcevic, we should add a note (directly on the README or adding a linked document to the README) about the limitations of using the internal 5pt_stencil routine mentioning:

1. The code is intended to be used internally in augur, and is not intended as a general numerical differentiation routine, for that we refer the user to general numerical differentiation packages such as numdifftools and derivkit.
2. The default (constant / scalar) step-size value was chosen arbitrarily (it works for a certain set of 3x2pt parameter + data vector configurations, but that's irrelevant here). The user needs to change its value accordingly for their own experiments, and verify the accuracy of the derivatives themselves. I guess something that use at your own risk -- we expect expert usage and the user needs to verify their results and cannot just trust blindly that this will work for all circumstances and set ups.

[nikosarcevic]

Hi @fjaviersanchez
thanks for bringing this up.
My suggestion is to to open up an arg (and kwarg) that chooses one of 2 or 3 derivative channels:

1. derivkit
2. numdifftools
3. native augur from diff utils

3 should be an option if and only if the issue with params for differnt shapes is fixed as i pointed out in our private conversations. Otherwise, this is giving incorrect results.

Kwargs for the two deivative tools should be an option with some defaults set in place. I think numdifftools has an option to change deriv order as well as to return some metadata. I cannot speak up, it has been years since I have used it

For derivkit, since people are very much used to stencils, I would choose a 5 pt stencil WITH riddders extrapolation turned on. That way the user does not need to mess with step size at all -- the extrapolation method will handle it.

We should aboslutely be clear about the limits of each method -- and explicitly state that it is user's responsibility to approach any derivative based inference with some basic knowledge about these problems (not blindly press play).

[nikosarcevic]

@bruno-moraes and Pedro asked about reproducing the errors with the 5pt stencil so I will paste the script here

import numpy as np


def five_pt_stencil(f, x0, h=1e-4):
    if (not isinstance(x0, np.ndarray)) and (not isinstance(x0, list)):
        der = 1.0 / (12.0 * h) * (
            f(x0 - 2 * h) - 8.0 * f(x0 - h) + 8.0 * f(x0 + h) - f(x0 + 2 * h)
        )
        return der
    elif isinstance(x0, list):
        x0 = np.array(x0)

    if isinstance(x0, np.ndarray):
        if x0.ndim == 1:
            xp2 = x0 + 2.0 * h * np.identity(len(x0), dtype=np.float64)[:, None]
            xm2 = x0 - 2.0 * h * np.identity(len(x0), dtype=np.float64)[:, None]
            xp1 = x0 + 1.0 * h * np.identity(len(x0), dtype=np.float64)[:, None]
            xm1 = x0 - 1.0 * h * np.identity(len(x0), dtype=np.float64)[:, None]
            xp2 = xp2[:, 0, :]
            xm2 = xm2[:, 0, :]
            xp1 = xp1[:, 0, :]
            xm1 = xm1[:, 0, :]
            der = 1.0 / (12.0 * h) * (f(xm2) - 8.0 * f(xm1) + 8.0 * f(xp1) - f(xp2))
            return der


def show(name, func):
    print(f"\n--- {name} ---")
    try:
        func()
    except Exception as e:
        print("raised:", repr(e))


# fixed absolute step size fails badly for mixed-scale paramters
show(
    "mixed scale fail",
    lambda: print(
        "actual:",
        five_pt_stencil(
            lambda x: np.array([row[0] ** 2 + row[1] for row in x]),
            np.array([1e-8, 1e8]),
            h=1e-4,
        ),
        "expected:",
        np.array([2e-8, 1.0]),
    ),
)


# func contract for f is unclear (this is really imprtant!!!)
# a natural 1D-only user function fails because the implemntation
# passes here a batched (2, n) array in the vector case
def f_requires_1d(x):
    x = np.asarray(x)
    if x.ndim != 1:
        raise ValueError(f"expected 1d input here, got shape {x.shape} "
                         f"isntead")
    return x[0] ** 2 + 3.0 * x[1]


show(
    "1D-only function fails",
    lambda: print(
        five_pt_stencil(f_requires_1d, np.array([1.0, 2.0]), h=1e-6)
    ),
)


# unsupported ndarray rank siletnly falls thru instead of raising
show(
    "2D input silently unsupported",
    lambda: print(
        "result:",
        five_pt_stencil(lambda x: x, np.array([[1.0, 2.0], [3.0, 4.0]])),
        "expected: ValueError",
    ),
)


# tiny param scale: deriv is numerically wrong because h is absurdly large
show(
    "tiny-scale derivative is wrong",
    lambda: print(
        "actual:",
        five_pt_stencil(
            lambda x: np.array([row[0] ** 2 + 3.0 * row[1] for row in x]),
            np.array([1e-12, 2.0]),
            h=1e-4,
        ),
        "expected:",
        np.array([2e-12, 3.0]),
    ),
)


# constant func in one coordinate should give exactly zero derivative
# there, but fixed h + floating-point effects can give garbage at bad scales
show(
    "derivtive should be zero but is numerically corrupted",
    lambda: print(
        "actual:",
        five_pt_stencil(
            lambda X: np.array([row[1] for row in X]),   # independent of x[0]
            np.array([1e-8, 1e8]),
            h=1e-4,
        ),
        "expected:",
        np.array([0.0, 1.0]),
    ),
)


# strong h-dependence on tiny-scale parameter
def f_scale(x):
    return np.array([row[0] ** 2 + 3.0 * row[1] for row in x])

x_tiny = np.array([1e-12, 2.0])

show(
    "strong h-dependence (tiny-scale parameter)",
    lambda: print(
        "h=1e-2 :", five_pt_stencil(f_scale, x_tiny, h=1e-2),
        "\nh=1e-4 :", five_pt_stencil(f_scale, x_tiny, h=1e-4),
        "\nh=1e-6 :", five_pt_stencil(f_scale, x_tiny, h=1e-6),
        "\nh=1e-8 :", five_pt_stencil(f_scale, x_tiny, h=1e-8),
        "\nexpected:", np.array([2e-12, 3.0]),
    ),
)


# this returns a deriv but it is wrong
def f_bad_scale(X):
    return np.array([row[0] ** 2 + row[1] for row in X])

x_bad = np.array([1e-8, 1e16])

show(
    "returns a derivative but it is nonsense",
    lambda: print(
        "atcual:",
        five_pt_stencil(f_bad_scale, x_bad, h=1e-4),
        "expected:",
        np.array([2e-8, 1.0]),
    ),
)

[fjaviersanchez]

Passing kwargs is now allowed for both numdifftools and derivkit. Added more docstrings / context + additional safeguards in five_pt_stencil in #138. For v2.0+ we might want to remove the method all along.