From b0a9c8245d1732bf95918ffb1507abaa140b219d Mon Sep 17 00:00:00 2001
From: andrewgsavage <andrewgsavage@gmail.com>
Date: Mon, 8 Jul 2024 20:46:53 +0100
Subject: [PATCH 1/7] add ufloatnumpy

---
 tests/helpers.py             |  12 ++
 tests/test_ufloatnumpy.py    | 239 ++++++++++++++++++++++++++++
 uncertainties/core.py        |   8 +-
 uncertainties/ufloatnumpy.py | 294 +++++++++++++++++++++++++++++++++++
 4 files changed, 552 insertions(+), 1 deletion(-)
 create mode 100644 tests/test_ufloatnumpy.py
 create mode 100644 uncertainties/ufloatnumpy.py

diff --git a/tests/helpers.py b/tests/helpers.py
index eec53dff..1f432cba 100644
--- a/tests/helpers.py
+++ b/tests/helpers.py
@@ -187,6 +187,18 @@ def numbers_close(x, y, tolerance=1e-6):
     else:  # Either x or y is zero
         return abs(x or y) < tolerance
 
+def nominal_and_std_dev_close(x, y, tolerance=1e-6):
+    '''
+    Tests if two numbers with uncertainties are close, NOT as random
+    variables. Checks whether the magnitude of the nominal
+    values and standard deviations are close.
+
+    The tolerance is applied to both the nominal value and the
+    standard deviation of the difference between the numbers.
+    '''
+    return (numbers_close(x.n, y.n, tolerance)
+            and numbers_close(x.s, y.s, tolerance))
+
 def ufloats_close(x, y, tolerance=1e-6):
     '''
     Tests if two numbers with uncertainties are close, as random
diff --git a/tests/test_ufloatnumpy.py b/tests/test_ufloatnumpy.py
new file mode 100644
index 00000000..64b86f2e
--- /dev/null
+++ b/tests/test_ufloatnumpy.py
@@ -0,0 +1,239 @@
+from uncertainties import unumpy, umath, ufloat, UFloat
+from helpers import (power_special_cases, power_all_cases, power_wrt_ref, numbers_close,
+    ufloats_close, compare_derivatives, uarrays_close, nominal_and_std_dev_close)
+import numpy as np
+import pytest 
+
+a = ufloat(1, 0.1)
+b = ufloat(2, 0.2)
+
+
+class TestArithmetic():
+    @pytest.mark.parametrize(
+        "first, second, expected",
+        [
+            (a, b, ufloat(3.0, 0.223606797749979)),
+            (a, a, ufloat(2.0, 0.2)),
+        ],
+    )
+    def test_add(self, first, second, expected):
+        result = first + second
+        assert nominal_and_std_dev_close(result, expected)
+
+        result = np.add(first, second)
+        assert nominal_and_std_dev_close(result, expected)
+
+    @pytest.mark.parametrize(
+        "first, second, expected",
+        [
+            (a, b, ufloat(-1.00, 0.223606797749979)),
+            (a, a, ufloat(0.0, 0.0)),
+        ],
+    )
+    def test_subtact(self, first, second, expected):
+        result = first - second
+        assert nominal_and_std_dev_close(result, expected)
+
+        result = np.subtract(first, second)
+        assert nominal_and_std_dev_close(result, expected)
+
+    @pytest.mark.parametrize(
+        "first, second, expected",
+        [
+            (a, b, ufloat(2.0, 0.28284271247461906)),
+            (a, a, ufloat(1.0, 0.2)),
+        ],
+    )
+    def test_multiply(self, first, second, expected):
+        result = first * second
+        assert nominal_and_std_dev_close(result, expected)
+
+        result = np.multiply(first, second)
+        assert nominal_and_std_dev_close(result, expected)
+
+    @pytest.mark.parametrize(
+        "first, second, expected",
+        [
+            (a, b, ufloat(0.5, 0.07071067811865477)),
+            (a, a, ufloat(1.0, 0.0)),
+        ],
+    )
+    def test_divide(self, first, second, expected):
+        result = first / second
+        assert nominal_and_std_dev_close(result, expected)
+
+        result = np.divide(first, second)
+        assert nominal_and_std_dev_close(result, expected)
+
+        result = np.true_divide(first, second)
+        assert nominal_and_std_dev_close(result, expected)
+
+        
+    @pytest.mark.parametrize(
+        "first, second, expected",
+        [
+            (a, b, ufloat(0.0, 0.0)),
+            (a, a, ufloat(1.0, 0.0)),
+        ],
+    )
+    def test_floor_divide(self, first, second, expected):
+        result = first // second
+        assert nominal_and_std_dev_close(result, expected)
+
+        result = np.floor_divide(first, second)
+        assert nominal_and_std_dev_close(result, expected)
+    
+
+class TestComparative():
+    @pytest.mark.parametrize(
+        "first, second, expected",
+        [
+            (a, b, False),
+            (a, a, True),
+        ],
+    )
+    def test_equal(self, first, second, expected):
+        result = first == second
+        assert result == expected
+
+        result = np.equal(first, second)
+        assert result == expected
+
+    @pytest.mark.parametrize(
+        "first, second, expected",
+        [
+            (a, b, True),
+            (a, a, False),
+        ],
+    )
+    def test_not_equal(self, first, second, expected):
+        result = first != second
+        assert result == expected
+
+        result = np.not_equal(first, second)
+        assert result == expected
+
+    @pytest.mark.parametrize(
+        "first, second, expected",
+        [
+            (a, b, True),
+            (a, a, False),
+        ],
+    )
+    def test_less(self, first, second, expected):
+        result = first < second
+        assert result == expected
+
+        result = np.less(first, second)
+        assert result == expected
+
+    @pytest.mark.parametrize(
+        "first, second, expected",
+        [
+            (a, b, True),
+            (a, a, True),
+        ],
+    )
+    def test_less_equal(self, first, second, expected):
+        result = first <= second
+        assert result == expected
+
+        result = np.less_equal(first, second)
+        assert result == expected
+
+    @pytest.mark.parametrize(
+        "first, second, expected",
+        [
+            (a, b, False),
+            (a, a, False),
+        ],
+    )
+    def test_greater(self, first, second, expected):
+        result = first > second
+        assert result == expected
+
+        result = np.greater(first, second)
+        assert result == expected
+
+    @pytest.mark.parametrize(
+        "first, second, expected",
+        [
+            (a, b, False),
+            (a, a, True),
+        ],
+    )
+    def test_greater_equal(self, first, second, expected):
+        result = first >= second
+        assert result == expected
+
+        result = np.greater_equal(first, second)
+        assert result == expected
+
+
+class TestUfuncs():
+    zero = ufloat(0.0, 0.1)
+    one = ufloat(1.0, 0.1)
+    pi_4 = ufloat(0.7853981633974483, 0.1)  # pi/4
+    pi_2 = ufloat(1.5707963267948966, 0.1)  # pi/2
+    @pytest.mark.parametrize(
+        "numpy_func, umath_func, arg, expected",
+        [
+            ('cos', 'cos',  zero, ufloat(1.0, 0.0)),
+            ('cos', 'cos',  pi_4, ufloat(0.7071067811865476, 0.07071067811865477)),
+            ('cos', 'cos',  pi_2, ufloat(6.123233995736766e-17, 0.1)),
+            ('cosh', 'cosh',  zero, ufloat(1.0, 0.0)),
+            ('cosh', 'cosh',  pi_4, ufloat(1.324609089252006, 0.08686709614860096)),
+            ('cosh', 'cosh',  pi_2, ufloat(2.5091784786580567, 0.2301298902307295)),
+            ('sin', 'sin',  zero, ufloat(0.0, 0.1)),
+            ('sin', 'sin',  pi_4, ufloat(0.7071067811865476, 0.07071067811865477)),
+            ('sin', 'sin',  pi_2, ufloat(1.0, 6.123233995736766e-18)),
+            ('sinh', 'sinh',  zero, ufloat(0.0, 0.1)),
+            ('sinh', 'sinh',  pi_4, ufloat(0.8686709614860095, 0.1324609089252006)),
+            ('sinh', 'sinh',  pi_2, ufloat(2.3012989023072947, 0.2509178478658057)),
+            ('tan', 'tan',  zero, ufloat(0.0, 0.1)),
+            ('tan', 'tan',  pi_4, ufloat(0.9999999999999999, 0.19999999999999998)),
+            ('tan', 'tan',  pi_2, ufloat(1.633123935319537e+16, 2.6670937881135717e+31)),
+            ('tanh', 'tanh',  zero, ufloat(0.0, 0.1)),
+            ('tanh', 'tanh',  pi_4, ufloat(0.6557942026326724, 0.05699339637933774)),
+            ('tanh', 'tanh',  pi_2, ufloat(0.9171523356672744, 0.015883159318006324)),
+            ('arccos', 'acos',  zero, ufloat(1.5707963267948966, 0.1)),
+            ('arccos', 'acos',  one, ufloat(0.0, float("nan"))),
+            ('arccosh', 'acosh',  one, ufloat(0.0, float("nan"))),
+            ('arcsin', 'asin',  zero, ufloat(0.0, 0.1)),
+            ('arcsin', 'asin',  one, ufloat(1.5707963267948966, float("nan"))),
+            ('arcsinh', 'asinh',  zero, ufloat(0.0, 0.1)),
+            ('arcsinh', 'asinh',  one, ufloat(0.8813735870195429, 0.07071067811865475)),
+            ('arctan', 'atan',  zero, ufloat(0.0, 0.1)),
+            ('arctan', 'atan',  one, ufloat(0.7853981633974483, 0.05)),
+            ('arctanh', 'atanh',  zero, ufloat(0.0, 0.1)),
+            ('exp', 'exp',  zero, ufloat(1.0, 0.1)),
+            ('exp', 'exp',  one, ufloat(2.718281828459045, 0.27182818284590454)),
+            ('exp2', None,  zero, ufloat(1.0, 0.06931471805599453)),
+            ('exp2', None,  one, ufloat(2.0, 0.13862943611198905)),
+            ('expm1', 'expm1',  zero, ufloat(0.0, 0.1)),
+            ('expm1', 'expm1',  one, ufloat(1.718281828459045, 0.27182818284590454)),
+            ('log10', 'log10',  one, ufloat(0.0, 0.04342944819032518)),
+            ('log1p', 'log1p',  zero, ufloat(0.0, 0.1)),
+            ('log1p', 'log1p',  one, ufloat(0.6931471805599453, 0.05)),
+            ('degrees', 'degrees',  zero, ufloat(0.0, 5.729577951308233)),
+            ('degrees', 'degrees',  one, ufloat(57.29577951308232, 5.729577951308233)),
+            ('radians', 'radians',  zero, ufloat(0.0, 0.0017453292519943296)),
+            ('radians', 'radians',  one, ufloat(0.017453292519943295, 0.0017453292519943296)),
+            ('rad2deg', 'degrees',  zero, ufloat(0.0, 5.729577951308233)),
+            ('rad2deg', 'degrees',  one, ufloat(57.29577951308232, 5.729577951308233)),
+            ('deg2rad', 'radians',  zero, ufloat(0.0, 0.0017453292519943296)),
+            ('deg2rad', 'radians',  one, ufloat(0.017453292519943295, 0.0017453292519943296)),
+            ('sqrt', 'sqrt',  zero, ufloat(0.0, float("nan"))),
+            ('sqrt', 'sqrt',  one, ufloat(1.0, 0.05)),
+        ],
+    )
+    def test_single_arg(self, numpy_func, umath_func, arg, expected):
+        func = getattr(np, numpy_func)
+        result = func(arg)
+        assert nominal_and_std_dev_close(result, expected)
+
+        if umath_func:
+            func = getattr(umath, umath_func)
+            result = func(arg)
+            assert nominal_and_std_dev_close(result, expected)
+        
\ No newline at end of file
diff --git a/uncertainties/core.py b/uncertainties/core.py
index c4815462..796b508e 100644
--- a/uncertainties/core.py
+++ b/uncertainties/core.py
@@ -41,6 +41,8 @@ def isinfinite(x):
 import numbers
 import collections
 
+from .ufloatnumpy import UFloatNumpy
+
 # The following restricts the local function getargspec() to the common
 # features of inspect.getargspec() and inspect.getfullargspec():
 if sys.version_info < (3,):  # !! Could be removed when moving to Python 3 only
@@ -1492,7 +1494,7 @@ def __getstate__(self):
     def __setstate__(self, state):
         (self.linear_combo,) = state
 
-class AffineScalarFunc(object):
+class AffineScalarFunc(UFloatNumpy):
     """
     Affine functions that support basic mathematical operations
     (addition, etc.).  Such functions can for instance be used for
@@ -2616,6 +2618,10 @@ def raise_error(self):
         setattr(AffineScalarFunc, '__%s__' % coercion_type, raise_error)
 
 add_operators_to_AffineScalarFunc()  # Actual addition of class attributes
+AffineScalarFunc._to_affine_scalar = to_affine_scalar
+AffineScalarFunc._add_numpy_arithmetic_ufuncs()
+AffineScalarFunc._add_numpy_comparative_ufuncs()
+AffineScalarFunc.wrap = wrap
 
 class NegativeStdDev(Exception):
     '''Raise for a negative standard deviation'''
diff --git a/uncertainties/ufloatnumpy.py b/uncertainties/ufloatnumpy.py
new file mode 100644
index 00000000..d12f198f
--- /dev/null
+++ b/uncertainties/ufloatnumpy.py
@@ -0,0 +1,294 @@
+import numpy as np
+import math
+# ufuncs are listed at https://numpy.org/doc/stable/reference/ufuncs.html
+from . import core as ops
+# from .umath_core import log_der0,_deriv_copysign, _deriv_fabs, _deriv_pow_0, _deriv_pow_1
+
+# from .core import nan_if_exception
+def nan_if_exception(f):
+    '''
+    Wrapper around f(x, y) that let f return NaN when f raises one of
+    a few numerical exceptions.
+    '''
+
+    def wrapped_f(*args, **kwargs):
+        try:
+            return f(*args, **kwargs)
+        except (ValueError, ZeroDivisionError, OverflowError):
+            return float('nan')
+
+    return wrapped_f
+
+
+def log_der0(*args):
+    """
+    Derivative of math.log() with respect to its first argument.
+
+    Works whether 1 or 2 arguments are given.
+    """
+    if len(args) == 1:
+        return 1/args[0]
+    else:
+        return 1/args[0]/math.log(args[1])  # 2-argument form
+
+    # The following version goes about as fast:
+
+    ## A 'try' is used for the most common case because it is fast when no
+    ## exception is raised:
+    #try:
+    #    return log_1arg_der(*args)  # Argument number check
+    #except TypeError:
+    #    return 1/args[0]/math.log(args[1])  # 2-argument form
+
+def _deriv_copysign(x,y):
+    if x >= 0:
+        return math.copysign(1, y)
+    else:
+        return -math.copysign(1, y)
+
+def _deriv_fabs(x):
+    if x >= 0:
+        return 1
+    else:
+        return -1
+
+def _deriv_pow_0(x, y):
+    if y == 0:
+        return  0.
+    elif x != 0 or y % 1 == 0:
+        return y*math.pow(x, y-1)
+    else:
+        return float('nan')
+
+def _deriv_pow_1(x, y):
+    if x == 0 and y > 0:
+        return 0.
+    else:
+        return math.log(x) * math.pow(x, y)
+
+def is_upcast_type(t):
+    # This can be used to allow downstream modules to overide operations; see pint
+    # TODO add upcast_type list or dict to a public interface
+    return False
+
+def implements(numpy_func_string, func_type):
+    """Register an __array_function__/__array_ufunc__ implementation for UArray
+    objects.
+    """
+    print(numpy_func_string, func_type)
+
+    def decorator(func):
+        if func_type == "function":
+            HANDLED_FUNCTIONS[numpy_func_string] = func
+        elif func_type == "ufunc":
+            HANDLED_UFUNCS[numpy_func_string] = func
+        else:
+            raise ValueError(f"Invalid func_type {func_type}")
+        return func
+
+    return decorator
+
+HANDLED_FUNCTIONS = {}
+HANDLED_UFUNCS = {}
+
+
+def apply_func_elementwise(func, inputs, kwargs, result_dtype="object"):
+    if len(inputs) == 1:
+        result = func(*inputs, **kwargs)
+    elif isinstance(inputs[0], np.ndarray):
+        result = np.empty_like(inputs[0], dtype=result_dtype)
+        for index, x in np.ndenumerate(inputs[0]):
+            inputs_ = [x if i == 0 else inputs[i] for i in range(len(inputs))]
+            result[index] = func(*inputs_, **kwargs)
+    elif isinstance(inputs[1], np.ndarray):
+        result = np.empty_like(inputs[1], dtype=result_dtype)
+        for index, x in np.ndenumerate(inputs[1]):
+            inputs_ = [x if i == 1 else inputs[i] for i in range(len(inputs))]
+            result[index] = func(*inputs_, **kwargs)
+    else:
+        result = func(*inputs, **kwargs)
+    return result
+
+def numpy_wrap(func_type, func, args, kwargs, types):
+    """Return the result from a NumPy function/ufunc as wrapped by uncertainties."""
+
+    if func_type == "function":
+        handled = HANDLED_FUNCTIONS
+        # Need to handle functions in submodules
+        name = ".".join(func.__module__.split(".")[1:] + [func.__name__])
+    elif func_type == "ufunc":
+        handled = HANDLED_UFUNCS
+        # ufuncs do not have func.__module__
+        name = func.__name__
+    else:
+        raise ValueError(f"Invalid func_type {func_type}")
+
+    if name not in handled or any(is_upcast_type(t) for t in types):
+        print("NotImplemented L54")
+        raise TypeError
+        return NotImplemented
+    return handled[name](*args, **kwargs)
+
+class UFloatNumpy(object):
+    # NumPy function/ufunc support
+    __array_priority__ = 17
+
+    def __array_ufunc__(self, ufunc, method, *inputs, **kwargs):
+        if method != "__call__":
+            # Only handle ufuncs as callables
+            return NotImplemented
+
+        # Replicate types from __array_function__
+        types = {
+            type(arg)
+            for arg in list(inputs) + list(kwargs.values())
+            if hasattr(arg, "__array_ufunc__")
+        }
+
+        return numpy_wrap("ufunc", ufunc, inputs, kwargs, types)
+
+    def __array_function__(self, func, types, args, kwargs):
+        return numpy_wrap("function", func, args, kwargs, types)
+
+    # original code for _add_numpy_ufuncs. may be helpful for writing a generic wraps
+    # this can be deleted:
+    # @classmethod
+    # def _add_numpy_ufuncs(cls):
+    #     def implement_ufunc(func_str, derivatives):
+    #         func = getattr(np, func_str)
+    #         @implements(func_str, "ufunc")
+    #         def implementation(*inputs, **kwargs):
+    #             if isinstance(inputs[0], np.ndarray):
+    #                 result = np.empty_like(inputs[0], dtype=object)
+    #                 for index, x in np.ndenumerate(inputs[0]):
+    #                     inputs_ = (x if i == 0 else inputs[i] for i in range(len(inputs)))
+    #                     result[index] = cls.wrap(func, derivatives)(*inputs_, **kwargs)
+    #             elif isinstance(inputs[1], np.ndarray):
+    #                 result = np.empty_like(inputs[1], dtype=object)
+    #                 for index, x in np.ndenumerate(inputs[1]):
+    #                     inputs_ = (x if i == 1 else inputs[i] for i in range(len(inputs)))
+    #                     result[index] = cls.wrap(func, derivatives)(*inputs_, **kwargs)
+    #             else:
+    #                 result = cls.wrap(func, derivatives)(*inputs, **kwargs)
+    #             return result
+
+    #         return implementation
+
+    #     for func_str, derivatives in ufunc_derivatives.items():
+    #         implement_ufunc(func_str, derivatives)
+
+
+    @classmethod
+    def _add_numpy_arithmetic_ufuncs(cls):
+        def implement_ufunc(func_str, derivatives):
+            func = getattr(np, func_str)
+            @implements(func_str, "ufunc")
+            def implementation(*inputs, **kwargs):
+                return apply_func_elementwise(
+                    cls.wrap(func, derivatives), inputs, kwargs)
+            return implementation
+
+        ufunc_derivatives = {
+            'add': [lambda x, y: 1., lambda x, y: 1.],
+            'subtract': [lambda x, y: 1., lambda x, y: -1.],
+            'multiply': [lambda x, y: y, lambda x, y: x],
+            'divide': [lambda x, y: 1./y, lambda x, y: -x/y**2],
+            'true_divide': [lambda x, y: 1./y, lambda x, y: -x/y**2],
+            'floor_divide': [lambda x, y: 0., lambda x, y: 0.],
+            
+            'arccos': [nan_if_exception(lambda x: -1/math.sqrt(1-x**2))],
+            'arccosh': [nan_if_exception(lambda x: 1/math.sqrt(x**2-1))],
+            'arcsin': [nan_if_exception(lambda x: 1/math.sqrt(1-x**2))],
+            'arcsinh': [nan_if_exception(lambda x: 1/math.sqrt(1+x**2))],
+            'arctan': [nan_if_exception(lambda x: 1/(1+x**2))],
+            'arctan2': [nan_if_exception(lambda y, x: x/(x**2+y**2)),  # Correct for x == 0
+                        nan_if_exception(lambda y, x: -y/(x**2+y**2))],  # Correct for x == 0
+            'arctanh': [nan_if_exception(lambda x: 1/(1-x**2))],
+            'cos': [lambda x: -math.sin(x)],
+            'cosh': [math.sinh],
+            'sin': [math.cos],
+            'sinh': [math.cosh],
+            'tan': [nan_if_exception(lambda x: 1+math.tan(x)**2)],
+            'tanh': [nan_if_exception(lambda x: 1-math.tanh(x)**2)],
+            'exp': [math.exp],
+            "exp2": [lambda y: _deriv_pow_1(2, y)],
+            'expm1': [math.exp],
+            'log10': [nan_if_exception(lambda x: 1/x/math.log(10))],
+            'log1p': [nan_if_exception(lambda x: 1/(1+x))],
+            'degrees': [lambda x: math.degrees(1)],
+            'rad2deg': [lambda x: math.degrees(1)],
+            'radians': [lambda x: math.radians(1)],
+            'deg2rad': [lambda x: math.radians(1)],
+            'power': [_deriv_pow_0, _deriv_pow_1],
+            'sqrt': [nan_if_exception(lambda x: 0.5/math.sqrt(x))],
+            'hypot': [lambda x, y: x/math.hypot(x, y),
+                      lambda x, y: y/math.hypot(x, y)],
+            }
+        # TODO: test arctan2, power, hypot
+        for func_str, derivatives in ufunc_derivatives.items():
+            implement_ufunc(func_str, derivatives)
+
+    @classmethod
+    def _add_numpy_comparative_ufuncs(cls):
+        def implement_ufunc(func_str, func):
+            @implements(func_str, "ufunc")
+            def implementation(*inputs, **kwargs):
+                inputs = [cls._to_affine_scalar(i) for i in inputs]
+                return apply_func_elementwise(func, inputs, kwargs, result_dtype=bool)
+
+            return implementation
+
+        ufunc_comparatives = {
+            'equal': ops.eq_on_aff_funcs,
+            'not_equal': ops.ne_on_aff_funcs,
+            'less': ops.lt_on_aff_funcs,
+            'less_equal': ops.le_on_aff_funcs,
+            'greater': ops.gt_on_aff_funcs,
+            'greater_equal': ops.ge_on_aff_funcs,
+        }
+        for func_str, func in ufunc_comparatives.items():
+            implement_ufunc(func_str, func)
+
+    # 'copysign': [_deriv_copysign,
+    #              lambda x, y: 0],
+    # 'degrees': [lambda x: math.degrees(1)],
+    # 'erf': [lambda x: math.exp(-x**2)*erf_coef],
+    # 'erfc': [lambda x: -math.exp(-x**2)*erf_coef],
+    # 'fabs': [_deriv_fabs],
+    # 'hypot': [lambda x, y: x/math.hypot(x, y),
+    #           lambda x, y: y/math.hypot(x, y)],
+    # 'log': [log_der0,
+    #         lambda x, y: -math.log(x, y)/y/math.log(y)],
+    # 'pow': [_deriv_pow_0, _deriv_pow_1],
+    # 'radians': [lambda x: math.radians(1)],
+    # 'sqrt': [lambda x: 0.5/math.sqrt(x)],
+    # def _add_trig_ufuncs(cls):
+
+    #         "cumprod": ("", ""),
+    #         "arccos": ("", "radian"),
+    #         "arcsin": ("", "radian"),
+    #         "arctan": ("", "radian"),
+    #         "arccosh": ("", "radian"),
+    #         "arcsinh": ("", "radian"),
+    #         "arctanh": ("", "radian"),
+    #         "exp": ("", ""),
+    #         "expm1": ("", ""),
+    #         "exp2": ("", ""),
+    #         "log": ("", ""),
+    #         "log10": ("", ""),
+    #         "log1p": ("", ""),
+    #         "log2": ("", ""),
+    #         "sin": ("radian", ""),
+    #         "cos": ("radian", ""),
+    #         "tan": ("radian", ""),
+    #         "sinh": ("radian", ""),
+    #         "cosh": ("radian", ""),
+    #         "tanh": ("radian", ""),
+    #         "radians": ("degree", "radian"),
+    #         "degrees": ("radian", "degree"),
+    #         "deg2rad": ("degree", "radian"),
+    #         "rad2deg": ("radian", "degree"),
+    #         "logaddexp": ("", ""),
+    #         "logaddexp2": ("", ""),
+
+            
\ No newline at end of file

From a0bbea5f1d7a3846a49f162c2426fe9b6451b5d0 Mon Sep 17 00:00:00 2001
From: andrewgsavage <andrewgsavage@gmail.com>
Date: Mon, 8 Jul 2024 20:46:53 +0100
Subject: [PATCH 2/7] add ufloatnumpy

---
 .coveragerc                  |   2 +
 tests/helpers.py             |  11 ++
 tests/test_ufloatnumpy.py    | 239 ++++++++++++++++++++++++++++
 uncertainties/core.py        |   5 +-
 uncertainties/ufloatnumpy.py | 294 +++++++++++++++++++++++++++++++++++
 5 files changed, 549 insertions(+), 2 deletions(-)
 create mode 100644 .coveragerc
 create mode 100644 tests/test_ufloatnumpy.py
 create mode 100644 uncertainties/ufloatnumpy.py

diff --git a/.coveragerc b/.coveragerc
new file mode 100644
index 00000000..5852b6c5
--- /dev/null
+++ b/.coveragerc
@@ -0,0 +1,2 @@
+[run]
+omit = uncertainties_pandas/testsuite/*
diff --git a/tests/helpers.py b/tests/helpers.py
index 7dc7fcea..79d1b6fc 100644
--- a/tests/helpers.py
+++ b/tests/helpers.py
@@ -190,6 +190,17 @@ def numbers_close(x, y, tolerance=1e-6):
     else:  # Either x or y is zero
         return abs(x or y) < tolerance
 
+def nominal_and_std_dev_close(x, y, tolerance=1e-6):
+    '''
+    Tests if two numbers with uncertainties are close, NOT as random
+    variables. Checks whether the magnitude of the nominal
+    values and standard deviations are close.
+
+    The tolerance is applied to both the nominal value and the
+    standard deviation of the difference between the numbers.
+    '''
+    return (numbers_close(x.n, y.n, tolerance)
+            and numbers_close(x.s, y.s, tolerance))
 
 def ufloats_close(x, y, tolerance=1e-6):
     """
diff --git a/tests/test_ufloatnumpy.py b/tests/test_ufloatnumpy.py
new file mode 100644
index 00000000..64b86f2e
--- /dev/null
+++ b/tests/test_ufloatnumpy.py
@@ -0,0 +1,239 @@
+from uncertainties import unumpy, umath, ufloat, UFloat
+from helpers import (power_special_cases, power_all_cases, power_wrt_ref, numbers_close,
+    ufloats_close, compare_derivatives, uarrays_close, nominal_and_std_dev_close)
+import numpy as np
+import pytest 
+
+a = ufloat(1, 0.1)
+b = ufloat(2, 0.2)
+
+
+class TestArithmetic():
+    @pytest.mark.parametrize(
+        "first, second, expected",
+        [
+            (a, b, ufloat(3.0, 0.223606797749979)),
+            (a, a, ufloat(2.0, 0.2)),
+        ],
+    )
+    def test_add(self, first, second, expected):
+        result = first + second
+        assert nominal_and_std_dev_close(result, expected)
+
+        result = np.add(first, second)
+        assert nominal_and_std_dev_close(result, expected)
+
+    @pytest.mark.parametrize(
+        "first, second, expected",
+        [
+            (a, b, ufloat(-1.00, 0.223606797749979)),
+            (a, a, ufloat(0.0, 0.0)),
+        ],
+    )
+    def test_subtact(self, first, second, expected):
+        result = first - second
+        assert nominal_and_std_dev_close(result, expected)
+
+        result = np.subtract(first, second)
+        assert nominal_and_std_dev_close(result, expected)
+
+    @pytest.mark.parametrize(
+        "first, second, expected",
+        [
+            (a, b, ufloat(2.0, 0.28284271247461906)),
+            (a, a, ufloat(1.0, 0.2)),
+        ],
+    )
+    def test_multiply(self, first, second, expected):
+        result = first * second
+        assert nominal_and_std_dev_close(result, expected)
+
+        result = np.multiply(first, second)
+        assert nominal_and_std_dev_close(result, expected)
+
+    @pytest.mark.parametrize(
+        "first, second, expected",
+        [
+            (a, b, ufloat(0.5, 0.07071067811865477)),
+            (a, a, ufloat(1.0, 0.0)),
+        ],
+    )
+    def test_divide(self, first, second, expected):
+        result = first / second
+        assert nominal_and_std_dev_close(result, expected)
+
+        result = np.divide(first, second)
+        assert nominal_and_std_dev_close(result, expected)
+
+        result = np.true_divide(first, second)
+        assert nominal_and_std_dev_close(result, expected)
+
+        
+    @pytest.mark.parametrize(
+        "first, second, expected",
+        [
+            (a, b, ufloat(0.0, 0.0)),
+            (a, a, ufloat(1.0, 0.0)),
+        ],
+    )
+    def test_floor_divide(self, first, second, expected):
+        result = first // second
+        assert nominal_and_std_dev_close(result, expected)
+
+        result = np.floor_divide(first, second)
+        assert nominal_and_std_dev_close(result, expected)
+    
+
+class TestComparative():
+    @pytest.mark.parametrize(
+        "first, second, expected",
+        [
+            (a, b, False),
+            (a, a, True),
+        ],
+    )
+    def test_equal(self, first, second, expected):
+        result = first == second
+        assert result == expected
+
+        result = np.equal(first, second)
+        assert result == expected
+
+    @pytest.mark.parametrize(
+        "first, second, expected",
+        [
+            (a, b, True),
+            (a, a, False),
+        ],
+    )
+    def test_not_equal(self, first, second, expected):
+        result = first != second
+        assert result == expected
+
+        result = np.not_equal(first, second)
+        assert result == expected
+
+    @pytest.mark.parametrize(
+        "first, second, expected",
+        [
+            (a, b, True),
+            (a, a, False),
+        ],
+    )
+    def test_less(self, first, second, expected):
+        result = first < second
+        assert result == expected
+
+        result = np.less(first, second)
+        assert result == expected
+
+    @pytest.mark.parametrize(
+        "first, second, expected",
+        [
+            (a, b, True),
+            (a, a, True),
+        ],
+    )
+    def test_less_equal(self, first, second, expected):
+        result = first <= second
+        assert result == expected
+
+        result = np.less_equal(first, second)
+        assert result == expected
+
+    @pytest.mark.parametrize(
+        "first, second, expected",
+        [
+            (a, b, False),
+            (a, a, False),
+        ],
+    )
+    def test_greater(self, first, second, expected):
+        result = first > second
+        assert result == expected
+
+        result = np.greater(first, second)
+        assert result == expected
+
+    @pytest.mark.parametrize(
+        "first, second, expected",
+        [
+            (a, b, False),
+            (a, a, True),
+        ],
+    )
+    def test_greater_equal(self, first, second, expected):
+        result = first >= second
+        assert result == expected
+
+        result = np.greater_equal(first, second)
+        assert result == expected
+
+
+class TestUfuncs():
+    zero = ufloat(0.0, 0.1)
+    one = ufloat(1.0, 0.1)
+    pi_4 = ufloat(0.7853981633974483, 0.1)  # pi/4
+    pi_2 = ufloat(1.5707963267948966, 0.1)  # pi/2
+    @pytest.mark.parametrize(
+        "numpy_func, umath_func, arg, expected",
+        [
+            ('cos', 'cos',  zero, ufloat(1.0, 0.0)),
+            ('cos', 'cos',  pi_4, ufloat(0.7071067811865476, 0.07071067811865477)),
+            ('cos', 'cos',  pi_2, ufloat(6.123233995736766e-17, 0.1)),
+            ('cosh', 'cosh',  zero, ufloat(1.0, 0.0)),
+            ('cosh', 'cosh',  pi_4, ufloat(1.324609089252006, 0.08686709614860096)),
+            ('cosh', 'cosh',  pi_2, ufloat(2.5091784786580567, 0.2301298902307295)),
+            ('sin', 'sin',  zero, ufloat(0.0, 0.1)),
+            ('sin', 'sin',  pi_4, ufloat(0.7071067811865476, 0.07071067811865477)),
+            ('sin', 'sin',  pi_2, ufloat(1.0, 6.123233995736766e-18)),
+            ('sinh', 'sinh',  zero, ufloat(0.0, 0.1)),
+            ('sinh', 'sinh',  pi_4, ufloat(0.8686709614860095, 0.1324609089252006)),
+            ('sinh', 'sinh',  pi_2, ufloat(2.3012989023072947, 0.2509178478658057)),
+            ('tan', 'tan',  zero, ufloat(0.0, 0.1)),
+            ('tan', 'tan',  pi_4, ufloat(0.9999999999999999, 0.19999999999999998)),
+            ('tan', 'tan',  pi_2, ufloat(1.633123935319537e+16, 2.6670937881135717e+31)),
+            ('tanh', 'tanh',  zero, ufloat(0.0, 0.1)),
+            ('tanh', 'tanh',  pi_4, ufloat(0.6557942026326724, 0.05699339637933774)),
+            ('tanh', 'tanh',  pi_2, ufloat(0.9171523356672744, 0.015883159318006324)),
+            ('arccos', 'acos',  zero, ufloat(1.5707963267948966, 0.1)),
+            ('arccos', 'acos',  one, ufloat(0.0, float("nan"))),
+            ('arccosh', 'acosh',  one, ufloat(0.0, float("nan"))),
+            ('arcsin', 'asin',  zero, ufloat(0.0, 0.1)),
+            ('arcsin', 'asin',  one, ufloat(1.5707963267948966, float("nan"))),
+            ('arcsinh', 'asinh',  zero, ufloat(0.0, 0.1)),
+            ('arcsinh', 'asinh',  one, ufloat(0.8813735870195429, 0.07071067811865475)),
+            ('arctan', 'atan',  zero, ufloat(0.0, 0.1)),
+            ('arctan', 'atan',  one, ufloat(0.7853981633974483, 0.05)),
+            ('arctanh', 'atanh',  zero, ufloat(0.0, 0.1)),
+            ('exp', 'exp',  zero, ufloat(1.0, 0.1)),
+            ('exp', 'exp',  one, ufloat(2.718281828459045, 0.27182818284590454)),
+            ('exp2', None,  zero, ufloat(1.0, 0.06931471805599453)),
+            ('exp2', None,  one, ufloat(2.0, 0.13862943611198905)),
+            ('expm1', 'expm1',  zero, ufloat(0.0, 0.1)),
+            ('expm1', 'expm1',  one, ufloat(1.718281828459045, 0.27182818284590454)),
+            ('log10', 'log10',  one, ufloat(0.0, 0.04342944819032518)),
+            ('log1p', 'log1p',  zero, ufloat(0.0, 0.1)),
+            ('log1p', 'log1p',  one, ufloat(0.6931471805599453, 0.05)),
+            ('degrees', 'degrees',  zero, ufloat(0.0, 5.729577951308233)),
+            ('degrees', 'degrees',  one, ufloat(57.29577951308232, 5.729577951308233)),
+            ('radians', 'radians',  zero, ufloat(0.0, 0.0017453292519943296)),
+            ('radians', 'radians',  one, ufloat(0.017453292519943295, 0.0017453292519943296)),
+            ('rad2deg', 'degrees',  zero, ufloat(0.0, 5.729577951308233)),
+            ('rad2deg', 'degrees',  one, ufloat(57.29577951308232, 5.729577951308233)),
+            ('deg2rad', 'radians',  zero, ufloat(0.0, 0.0017453292519943296)),
+            ('deg2rad', 'radians',  one, ufloat(0.017453292519943295, 0.0017453292519943296)),
+            ('sqrt', 'sqrt',  zero, ufloat(0.0, float("nan"))),
+            ('sqrt', 'sqrt',  one, ufloat(1.0, 0.05)),
+        ],
+    )
+    def test_single_arg(self, numpy_func, umath_func, arg, expected):
+        func = getattr(np, numpy_func)
+        result = func(arg)
+        assert nominal_and_std_dev_close(result, expected)
+
+        if umath_func:
+            func = getattr(umath, umath_func)
+            result = func(arg)
+            assert nominal_and_std_dev_close(result, expected)
+        
\ No newline at end of file
diff --git a/uncertainties/core.py b/uncertainties/core.py
index fdd64db7..06739926 100644
--- a/uncertainties/core.py
+++ b/uncertainties/core.py
@@ -39,6 +39,8 @@ def isinfinite(x):
     modified_operators,
     modified_ops_with_reflection,
 )
+from .ufloatnumpy import UFloatNumpy
+
 
 # Attributes that are always exported (some other attributes are
 # exported only if the NumPy module is available...):
@@ -331,8 +333,7 @@ def __getstate__(self):
     def __setstate__(self, state):
         (self.linear_combo,) = state
 
-
-class AffineScalarFunc(object):
+class AffineScalarFunc(UFloatNumpy):
     """
     Affine functions that support basic mathematical operations
     (addition, etc.).  Such functions can for instance be used for
diff --git a/uncertainties/ufloatnumpy.py b/uncertainties/ufloatnumpy.py
new file mode 100644
index 00000000..d12f198f
--- /dev/null
+++ b/uncertainties/ufloatnumpy.py
@@ -0,0 +1,294 @@
+import numpy as np
+import math
+# ufuncs are listed at https://numpy.org/doc/stable/reference/ufuncs.html
+from . import core as ops
+# from .umath_core import log_der0,_deriv_copysign, _deriv_fabs, _deriv_pow_0, _deriv_pow_1
+
+# from .core import nan_if_exception
+def nan_if_exception(f):
+    '''
+    Wrapper around f(x, y) that let f return NaN when f raises one of
+    a few numerical exceptions.
+    '''
+
+    def wrapped_f(*args, **kwargs):
+        try:
+            return f(*args, **kwargs)
+        except (ValueError, ZeroDivisionError, OverflowError):
+            return float('nan')
+
+    return wrapped_f
+
+
+def log_der0(*args):
+    """
+    Derivative of math.log() with respect to its first argument.
+
+    Works whether 1 or 2 arguments are given.
+    """
+    if len(args) == 1:
+        return 1/args[0]
+    else:
+        return 1/args[0]/math.log(args[1])  # 2-argument form
+
+    # The following version goes about as fast:
+
+    ## A 'try' is used for the most common case because it is fast when no
+    ## exception is raised:
+    #try:
+    #    return log_1arg_der(*args)  # Argument number check
+    #except TypeError:
+    #    return 1/args[0]/math.log(args[1])  # 2-argument form
+
+def _deriv_copysign(x,y):
+    if x >= 0:
+        return math.copysign(1, y)
+    else:
+        return -math.copysign(1, y)
+
+def _deriv_fabs(x):
+    if x >= 0:
+        return 1
+    else:
+        return -1
+
+def _deriv_pow_0(x, y):
+    if y == 0:
+        return  0.
+    elif x != 0 or y % 1 == 0:
+        return y*math.pow(x, y-1)
+    else:
+        return float('nan')
+
+def _deriv_pow_1(x, y):
+    if x == 0 and y > 0:
+        return 0.
+    else:
+        return math.log(x) * math.pow(x, y)
+
+def is_upcast_type(t):
+    # This can be used to allow downstream modules to overide operations; see pint
+    # TODO add upcast_type list or dict to a public interface
+    return False
+
+def implements(numpy_func_string, func_type):
+    """Register an __array_function__/__array_ufunc__ implementation for UArray
+    objects.
+    """
+    print(numpy_func_string, func_type)
+
+    def decorator(func):
+        if func_type == "function":
+            HANDLED_FUNCTIONS[numpy_func_string] = func
+        elif func_type == "ufunc":
+            HANDLED_UFUNCS[numpy_func_string] = func
+        else:
+            raise ValueError(f"Invalid func_type {func_type}")
+        return func
+
+    return decorator
+
+HANDLED_FUNCTIONS = {}
+HANDLED_UFUNCS = {}
+
+
+def apply_func_elementwise(func, inputs, kwargs, result_dtype="object"):
+    if len(inputs) == 1:
+        result = func(*inputs, **kwargs)
+    elif isinstance(inputs[0], np.ndarray):
+        result = np.empty_like(inputs[0], dtype=result_dtype)
+        for index, x in np.ndenumerate(inputs[0]):
+            inputs_ = [x if i == 0 else inputs[i] for i in range(len(inputs))]
+            result[index] = func(*inputs_, **kwargs)
+    elif isinstance(inputs[1], np.ndarray):
+        result = np.empty_like(inputs[1], dtype=result_dtype)
+        for index, x in np.ndenumerate(inputs[1]):
+            inputs_ = [x if i == 1 else inputs[i] for i in range(len(inputs))]
+            result[index] = func(*inputs_, **kwargs)
+    else:
+        result = func(*inputs, **kwargs)
+    return result
+
+def numpy_wrap(func_type, func, args, kwargs, types):
+    """Return the result from a NumPy function/ufunc as wrapped by uncertainties."""
+
+    if func_type == "function":
+        handled = HANDLED_FUNCTIONS
+        # Need to handle functions in submodules
+        name = ".".join(func.__module__.split(".")[1:] + [func.__name__])
+    elif func_type == "ufunc":
+        handled = HANDLED_UFUNCS
+        # ufuncs do not have func.__module__
+        name = func.__name__
+    else:
+        raise ValueError(f"Invalid func_type {func_type}")
+
+    if name not in handled or any(is_upcast_type(t) for t in types):
+        print("NotImplemented L54")
+        raise TypeError
+        return NotImplemented
+    return handled[name](*args, **kwargs)
+
+class UFloatNumpy(object):
+    # NumPy function/ufunc support
+    __array_priority__ = 17
+
+    def __array_ufunc__(self, ufunc, method, *inputs, **kwargs):
+        if method != "__call__":
+            # Only handle ufuncs as callables
+            return NotImplemented
+
+        # Replicate types from __array_function__
+        types = {
+            type(arg)
+            for arg in list(inputs) + list(kwargs.values())
+            if hasattr(arg, "__array_ufunc__")
+        }
+
+        return numpy_wrap("ufunc", ufunc, inputs, kwargs, types)
+
+    def __array_function__(self, func, types, args, kwargs):
+        return numpy_wrap("function", func, args, kwargs, types)
+
+    # original code for _add_numpy_ufuncs. may be helpful for writing a generic wraps
+    # this can be deleted:
+    # @classmethod
+    # def _add_numpy_ufuncs(cls):
+    #     def implement_ufunc(func_str, derivatives):
+    #         func = getattr(np, func_str)
+    #         @implements(func_str, "ufunc")
+    #         def implementation(*inputs, **kwargs):
+    #             if isinstance(inputs[0], np.ndarray):
+    #                 result = np.empty_like(inputs[0], dtype=object)
+    #                 for index, x in np.ndenumerate(inputs[0]):
+    #                     inputs_ = (x if i == 0 else inputs[i] for i in range(len(inputs)))
+    #                     result[index] = cls.wrap(func, derivatives)(*inputs_, **kwargs)
+    #             elif isinstance(inputs[1], np.ndarray):
+    #                 result = np.empty_like(inputs[1], dtype=object)
+    #                 for index, x in np.ndenumerate(inputs[1]):
+    #                     inputs_ = (x if i == 1 else inputs[i] for i in range(len(inputs)))
+    #                     result[index] = cls.wrap(func, derivatives)(*inputs_, **kwargs)
+    #             else:
+    #                 result = cls.wrap(func, derivatives)(*inputs, **kwargs)
+    #             return result
+
+    #         return implementation
+
+    #     for func_str, derivatives in ufunc_derivatives.items():
+    #         implement_ufunc(func_str, derivatives)
+
+
+    @classmethod
+    def _add_numpy_arithmetic_ufuncs(cls):
+        def implement_ufunc(func_str, derivatives):
+            func = getattr(np, func_str)
+            @implements(func_str, "ufunc")
+            def implementation(*inputs, **kwargs):
+                return apply_func_elementwise(
+                    cls.wrap(func, derivatives), inputs, kwargs)
+            return implementation
+
+        ufunc_derivatives = {
+            'add': [lambda x, y: 1., lambda x, y: 1.],
+            'subtract': [lambda x, y: 1., lambda x, y: -1.],
+            'multiply': [lambda x, y: y, lambda x, y: x],
+            'divide': [lambda x, y: 1./y, lambda x, y: -x/y**2],
+            'true_divide': [lambda x, y: 1./y, lambda x, y: -x/y**2],
+            'floor_divide': [lambda x, y: 0., lambda x, y: 0.],
+            
+            'arccos': [nan_if_exception(lambda x: -1/math.sqrt(1-x**2))],
+            'arccosh': [nan_if_exception(lambda x: 1/math.sqrt(x**2-1))],
+            'arcsin': [nan_if_exception(lambda x: 1/math.sqrt(1-x**2))],
+            'arcsinh': [nan_if_exception(lambda x: 1/math.sqrt(1+x**2))],
+            'arctan': [nan_if_exception(lambda x: 1/(1+x**2))],
+            'arctan2': [nan_if_exception(lambda y, x: x/(x**2+y**2)),  # Correct for x == 0
+                        nan_if_exception(lambda y, x: -y/(x**2+y**2))],  # Correct for x == 0
+            'arctanh': [nan_if_exception(lambda x: 1/(1-x**2))],
+            'cos': [lambda x: -math.sin(x)],
+            'cosh': [math.sinh],
+            'sin': [math.cos],
+            'sinh': [math.cosh],
+            'tan': [nan_if_exception(lambda x: 1+math.tan(x)**2)],
+            'tanh': [nan_if_exception(lambda x: 1-math.tanh(x)**2)],
+            'exp': [math.exp],
+            "exp2": [lambda y: _deriv_pow_1(2, y)],
+            'expm1': [math.exp],
+            'log10': [nan_if_exception(lambda x: 1/x/math.log(10))],
+            'log1p': [nan_if_exception(lambda x: 1/(1+x))],
+            'degrees': [lambda x: math.degrees(1)],
+            'rad2deg': [lambda x: math.degrees(1)],
+            'radians': [lambda x: math.radians(1)],
+            'deg2rad': [lambda x: math.radians(1)],
+            'power': [_deriv_pow_0, _deriv_pow_1],
+            'sqrt': [nan_if_exception(lambda x: 0.5/math.sqrt(x))],
+            'hypot': [lambda x, y: x/math.hypot(x, y),
+                      lambda x, y: y/math.hypot(x, y)],
+            }
+        # TODO: test arctan2, power, hypot
+        for func_str, derivatives in ufunc_derivatives.items():
+            implement_ufunc(func_str, derivatives)
+
+    @classmethod
+    def _add_numpy_comparative_ufuncs(cls):
+        def implement_ufunc(func_str, func):
+            @implements(func_str, "ufunc")
+            def implementation(*inputs, **kwargs):
+                inputs = [cls._to_affine_scalar(i) for i in inputs]
+                return apply_func_elementwise(func, inputs, kwargs, result_dtype=bool)
+
+            return implementation
+
+        ufunc_comparatives = {
+            'equal': ops.eq_on_aff_funcs,
+            'not_equal': ops.ne_on_aff_funcs,
+            'less': ops.lt_on_aff_funcs,
+            'less_equal': ops.le_on_aff_funcs,
+            'greater': ops.gt_on_aff_funcs,
+            'greater_equal': ops.ge_on_aff_funcs,
+        }
+        for func_str, func in ufunc_comparatives.items():
+            implement_ufunc(func_str, func)
+
+    # 'copysign': [_deriv_copysign,
+    #              lambda x, y: 0],
+    # 'degrees': [lambda x: math.degrees(1)],
+    # 'erf': [lambda x: math.exp(-x**2)*erf_coef],
+    # 'erfc': [lambda x: -math.exp(-x**2)*erf_coef],
+    # 'fabs': [_deriv_fabs],
+    # 'hypot': [lambda x, y: x/math.hypot(x, y),
+    #           lambda x, y: y/math.hypot(x, y)],
+    # 'log': [log_der0,
+    #         lambda x, y: -math.log(x, y)/y/math.log(y)],
+    # 'pow': [_deriv_pow_0, _deriv_pow_1],
+    # 'radians': [lambda x: math.radians(1)],
+    # 'sqrt': [lambda x: 0.5/math.sqrt(x)],
+    # def _add_trig_ufuncs(cls):
+
+    #         "cumprod": ("", ""),
+    #         "arccos": ("", "radian"),
+    #         "arcsin": ("", "radian"),
+    #         "arctan": ("", "radian"),
+    #         "arccosh": ("", "radian"),
+    #         "arcsinh": ("", "radian"),
+    #         "arctanh": ("", "radian"),
+    #         "exp": ("", ""),
+    #         "expm1": ("", ""),
+    #         "exp2": ("", ""),
+    #         "log": ("", ""),
+    #         "log10": ("", ""),
+    #         "log1p": ("", ""),
+    #         "log2": ("", ""),
+    #         "sin": ("radian", ""),
+    #         "cos": ("radian", ""),
+    #         "tan": ("radian", ""),
+    #         "sinh": ("radian", ""),
+    #         "cosh": ("radian", ""),
+    #         "tanh": ("radian", ""),
+    #         "radians": ("degree", "radian"),
+    #         "degrees": ("radian", "degree"),
+    #         "deg2rad": ("degree", "radian"),
+    #         "rad2deg": ("radian", "degree"),
+    #         "logaddexp": ("", ""),
+    #         "logaddexp2": ("", ""),
+
+            
\ No newline at end of file

From 637e3d4c8d64f6800ef0bca19b6e9dcecdf80da8 Mon Sep 17 00:00:00 2001
From: Andrew <andrewgsavage@gmail.com>
Date: Fri, 12 Jul 2024 09:24:46 +0100
Subject: [PATCH 3/7] tests

---
 uncertainties/core.py        |  2 ++
 uncertainties/ufloatnumpy.py | 25 ++++++++-----------------
 2 files changed, 10 insertions(+), 17 deletions(-)

diff --git a/uncertainties/core.py b/uncertainties/core.py
index 06739926..37436aec 100644
--- a/uncertainties/core.py
+++ b/uncertainties/core.py
@@ -658,6 +658,8 @@ def __setstate__(self, data_dict):
 ops.add_arithmetic_ops(AffineScalarFunc)
 ops.add_comparative_ops(AffineScalarFunc)
 to_affine_scalar = AffineScalarFunc._to_affine_scalar
+AffineScalarFunc._add_numpy_arithmetic_ufuncs()
+AffineScalarFunc._add_numpy_comparative_ufuncs()
 
 # Nicer name, for users: isinstance(ufloat(...), UFloat) is
 # True. Also: isinstance(..., UFloat) is the test for "is this a
diff --git a/uncertainties/ufloatnumpy.py b/uncertainties/ufloatnumpy.py
index d12f198f..3c6be355 100644
--- a/uncertainties/ufloatnumpy.py
+++ b/uncertainties/ufloatnumpy.py
@@ -1,23 +1,10 @@
 import numpy as np
 import math
 # ufuncs are listed at https://numpy.org/doc/stable/reference/ufuncs.html
-from . import core as ops
+from . import ops
 # from .umath_core import log_der0,_deriv_copysign, _deriv_fabs, _deriv_pow_0, _deriv_pow_1
 
-# from .core import nan_if_exception
-def nan_if_exception(f):
-    '''
-    Wrapper around f(x, y) that let f return NaN when f raises one of
-    a few numerical exceptions.
-    '''
-
-    def wrapped_f(*args, **kwargs):
-        try:
-            return f(*args, **kwargs)
-        except (ValueError, ZeroDivisionError, OverflowError):
-            return float('nan')
-
-    return wrapped_f
+from .ops import nan_if_exception
 
 
 def log_der0(*args):
@@ -88,7 +75,7 @@ def decorator(func):
 
     return decorator
 
-HANDLED_FUNCTIONS = {}
+HANDLED_FUNCTIONS = {} # noqa
 HANDLED_UFUNCS = {}
 
 
@@ -100,11 +87,15 @@ def apply_func_elementwise(func, inputs, kwargs, result_dtype="object"):
         for index, x in np.ndenumerate(inputs[0]):
             inputs_ = [x if i == 0 else inputs[i] for i in range(len(inputs))]
             result[index] = func(*inputs_, **kwargs)
+        if inputs[0].ndim == 0:
+            result = result.item()
     elif isinstance(inputs[1], np.ndarray):
         result = np.empty_like(inputs[1], dtype=result_dtype)
         for index, x in np.ndenumerate(inputs[1]):
             inputs_ = [x if i == 1 else inputs[i] for i in range(len(inputs))]
             result[index] = func(*inputs_, **kwargs)
+        if inputs[1].ndim == 0:
+            result = result.item()
     else:
         result = func(*inputs, **kwargs)
     return result
@@ -185,7 +176,7 @@ def implement_ufunc(func_str, derivatives):
             @implements(func_str, "ufunc")
             def implementation(*inputs, **kwargs):
                 return apply_func_elementwise(
-                    cls.wrap(func, derivatives), inputs, kwargs)
+                    ops._wrap(cls, func, derivatives), inputs, kwargs)
             return implementation
 
         ufunc_derivatives = {

From 15d2d07ffdd546bfe6485054aa5d22860fcef1a3 Mon Sep 17 00:00:00 2001
From: Andrew <andrewgsavage@gmail.com>
Date: Fri, 12 Jul 2024 10:15:23 +0100
Subject: [PATCH 4/7] pass tests

---
 .coveragerc                  |  2 --
 uncertainties/ufloatnumpy.py | 12 +++++++++++-
 2 files changed, 11 insertions(+), 3 deletions(-)
 delete mode 100644 .coveragerc

diff --git a/.coveragerc b/.coveragerc
deleted file mode 100644
index 5852b6c5..00000000
--- a/.coveragerc
+++ /dev/null
@@ -1,2 +0,0 @@
-[run]
-omit = uncertainties_pandas/testsuite/*
diff --git a/uncertainties/ufloatnumpy.py b/uncertainties/ufloatnumpy.py
index 3c6be355..6c4b3321 100644
--- a/uncertainties/ufloatnumpy.py
+++ b/uncertainties/ufloatnumpy.py
@@ -87,6 +87,7 @@ def apply_func_elementwise(func, inputs, kwargs, result_dtype="object"):
         for index, x in np.ndenumerate(inputs[0]):
             inputs_ = [x if i == 0 else inputs[i] for i in range(len(inputs))]
             result[index] = func(*inputs_, **kwargs)
+        # unpack the result of operations with ndim=0 arrays
         if inputs[0].ndim == 0:
             result = result.item()
     elif isinstance(inputs[1], np.ndarray):
@@ -94,6 +95,7 @@ def apply_func_elementwise(func, inputs, kwargs, result_dtype="object"):
         for index, x in np.ndenumerate(inputs[1]):
             inputs_ = [x if i == 1 else inputs[i] for i in range(len(inputs))]
             result[index] = func(*inputs_, **kwargs)
+        # unpack the result of operations with ndim=0 arrays
         if inputs[1].ndim == 0:
             result = result.item()
     else:
@@ -221,10 +223,18 @@ def implementation(*inputs, **kwargs):
 
     @classmethod
     def _add_numpy_comparative_ufuncs(cls):
+        def recursive_to_affine_scalar(arr):
+            print(arr)
+            if isinstance(arr, (list, tuple)):
+                return type(arr)([recursive_to_affine_scalar(i) for i in arr])
+            if isinstance(arr, np.ndarray):
+                return np.array([recursive_to_affine_scalar(i) for i in arr], "object")
+            return cls._to_affine_scalar(arr)
+
         def implement_ufunc(func_str, func):
             @implements(func_str, "ufunc")
             def implementation(*inputs, **kwargs):
-                inputs = [cls._to_affine_scalar(i) for i in inputs]
+                inputs = recursive_to_affine_scalar(inputs)
                 return apply_func_elementwise(func, inputs, kwargs, result_dtype=bool)
 
             return implementation

From 9ae0a0bbd107df47d7d7cae42d8901ad3c449cec Mon Sep 17 00:00:00 2001
From: Andrew <andrewgsavage@gmail.com>
Date: Fri, 12 Jul 2024 10:20:14 +0100
Subject: [PATCH 5/7] lint

---
 tests/helpers.py             |   9 +--
 tests/test_ufloatnumpy.py    | 122 +++++++++++++++++++----------------
 uncertainties/core.py        |   1 +
 uncertainties/ufloatnumpy.py | 121 ++++++++++++++++++----------------
 4 files changed, 138 insertions(+), 115 deletions(-)

diff --git a/tests/helpers.py b/tests/helpers.py
index 79d1b6fc..27155442 100644
--- a/tests/helpers.py
+++ b/tests/helpers.py
@@ -190,17 +190,18 @@ def numbers_close(x, y, tolerance=1e-6):
     else:  # Either x or y is zero
         return abs(x or y) < tolerance
 
+
 def nominal_and_std_dev_close(x, y, tolerance=1e-6):
-    '''
+    """
     Tests if two numbers with uncertainties are close, NOT as random
     variables. Checks whether the magnitude of the nominal
     values and standard deviations are close.
 
     The tolerance is applied to both the nominal value and the
     standard deviation of the difference between the numbers.
-    '''
-    return (numbers_close(x.n, y.n, tolerance)
-            and numbers_close(x.s, y.s, tolerance))
+    """
+    return numbers_close(x.n, y.n, tolerance) and numbers_close(x.s, y.s, tolerance)
+
 
 def ufloats_close(x, y, tolerance=1e-6):
     """
diff --git a/tests/test_ufloatnumpy.py b/tests/test_ufloatnumpy.py
index 64b86f2e..75eafd71 100644
--- a/tests/test_ufloatnumpy.py
+++ b/tests/test_ufloatnumpy.py
@@ -1,14 +1,13 @@
-from uncertainties import unumpy, umath, ufloat, UFloat
-from helpers import (power_special_cases, power_all_cases, power_wrt_ref, numbers_close,
-    ufloats_close, compare_derivatives, uarrays_close, nominal_and_std_dev_close)
+from uncertainties import umath, ufloat
+from helpers import nominal_and_std_dev_close
 import numpy as np
-import pytest 
+import pytest
 
 a = ufloat(1, 0.1)
 b = ufloat(2, 0.2)
 
 
-class TestArithmetic():
+class TestArithmetic:
     @pytest.mark.parametrize(
         "first, second, expected",
         [
@@ -68,7 +67,6 @@ def test_divide(self, first, second, expected):
         result = np.true_divide(first, second)
         assert nominal_and_std_dev_close(result, expected)
 
-        
     @pytest.mark.parametrize(
         "first, second, expected",
         [
@@ -82,9 +80,9 @@ def test_floor_divide(self, first, second, expected):
 
         result = np.floor_divide(first, second)
         assert nominal_and_std_dev_close(result, expected)
-    
 
-class TestComparative():
+
+class TestComparative:
     @pytest.mark.parametrize(
         "first, second, expected",
         [
@@ -170,61 +168,72 @@ def test_greater_equal(self, first, second, expected):
         assert result == expected
 
 
-class TestUfuncs():
+class TestUfuncs:
     zero = ufloat(0.0, 0.1)
     one = ufloat(1.0, 0.1)
     pi_4 = ufloat(0.7853981633974483, 0.1)  # pi/4
     pi_2 = ufloat(1.5707963267948966, 0.1)  # pi/2
+
     @pytest.mark.parametrize(
         "numpy_func, umath_func, arg, expected",
         [
-            ('cos', 'cos',  zero, ufloat(1.0, 0.0)),
-            ('cos', 'cos',  pi_4, ufloat(0.7071067811865476, 0.07071067811865477)),
-            ('cos', 'cos',  pi_2, ufloat(6.123233995736766e-17, 0.1)),
-            ('cosh', 'cosh',  zero, ufloat(1.0, 0.0)),
-            ('cosh', 'cosh',  pi_4, ufloat(1.324609089252006, 0.08686709614860096)),
-            ('cosh', 'cosh',  pi_2, ufloat(2.5091784786580567, 0.2301298902307295)),
-            ('sin', 'sin',  zero, ufloat(0.0, 0.1)),
-            ('sin', 'sin',  pi_4, ufloat(0.7071067811865476, 0.07071067811865477)),
-            ('sin', 'sin',  pi_2, ufloat(1.0, 6.123233995736766e-18)),
-            ('sinh', 'sinh',  zero, ufloat(0.0, 0.1)),
-            ('sinh', 'sinh',  pi_4, ufloat(0.8686709614860095, 0.1324609089252006)),
-            ('sinh', 'sinh',  pi_2, ufloat(2.3012989023072947, 0.2509178478658057)),
-            ('tan', 'tan',  zero, ufloat(0.0, 0.1)),
-            ('tan', 'tan',  pi_4, ufloat(0.9999999999999999, 0.19999999999999998)),
-            ('tan', 'tan',  pi_2, ufloat(1.633123935319537e+16, 2.6670937881135717e+31)),
-            ('tanh', 'tanh',  zero, ufloat(0.0, 0.1)),
-            ('tanh', 'tanh',  pi_4, ufloat(0.6557942026326724, 0.05699339637933774)),
-            ('tanh', 'tanh',  pi_2, ufloat(0.9171523356672744, 0.015883159318006324)),
-            ('arccos', 'acos',  zero, ufloat(1.5707963267948966, 0.1)),
-            ('arccos', 'acos',  one, ufloat(0.0, float("nan"))),
-            ('arccosh', 'acosh',  one, ufloat(0.0, float("nan"))),
-            ('arcsin', 'asin',  zero, ufloat(0.0, 0.1)),
-            ('arcsin', 'asin',  one, ufloat(1.5707963267948966, float("nan"))),
-            ('arcsinh', 'asinh',  zero, ufloat(0.0, 0.1)),
-            ('arcsinh', 'asinh',  one, ufloat(0.8813735870195429, 0.07071067811865475)),
-            ('arctan', 'atan',  zero, ufloat(0.0, 0.1)),
-            ('arctan', 'atan',  one, ufloat(0.7853981633974483, 0.05)),
-            ('arctanh', 'atanh',  zero, ufloat(0.0, 0.1)),
-            ('exp', 'exp',  zero, ufloat(1.0, 0.1)),
-            ('exp', 'exp',  one, ufloat(2.718281828459045, 0.27182818284590454)),
-            ('exp2', None,  zero, ufloat(1.0, 0.06931471805599453)),
-            ('exp2', None,  one, ufloat(2.0, 0.13862943611198905)),
-            ('expm1', 'expm1',  zero, ufloat(0.0, 0.1)),
-            ('expm1', 'expm1',  one, ufloat(1.718281828459045, 0.27182818284590454)),
-            ('log10', 'log10',  one, ufloat(0.0, 0.04342944819032518)),
-            ('log1p', 'log1p',  zero, ufloat(0.0, 0.1)),
-            ('log1p', 'log1p',  one, ufloat(0.6931471805599453, 0.05)),
-            ('degrees', 'degrees',  zero, ufloat(0.0, 5.729577951308233)),
-            ('degrees', 'degrees',  one, ufloat(57.29577951308232, 5.729577951308233)),
-            ('radians', 'radians',  zero, ufloat(0.0, 0.0017453292519943296)),
-            ('radians', 'radians',  one, ufloat(0.017453292519943295, 0.0017453292519943296)),
-            ('rad2deg', 'degrees',  zero, ufloat(0.0, 5.729577951308233)),
-            ('rad2deg', 'degrees',  one, ufloat(57.29577951308232, 5.729577951308233)),
-            ('deg2rad', 'radians',  zero, ufloat(0.0, 0.0017453292519943296)),
-            ('deg2rad', 'radians',  one, ufloat(0.017453292519943295, 0.0017453292519943296)),
-            ('sqrt', 'sqrt',  zero, ufloat(0.0, float("nan"))),
-            ('sqrt', 'sqrt',  one, ufloat(1.0, 0.05)),
+            ("cos", "cos", zero, ufloat(1.0, 0.0)),
+            ("cos", "cos", pi_4, ufloat(0.7071067811865476, 0.07071067811865477)),
+            ("cos", "cos", pi_2, ufloat(6.123233995736766e-17, 0.1)),
+            ("cosh", "cosh", zero, ufloat(1.0, 0.0)),
+            ("cosh", "cosh", pi_4, ufloat(1.324609089252006, 0.08686709614860096)),
+            ("cosh", "cosh", pi_2, ufloat(2.5091784786580567, 0.2301298902307295)),
+            ("sin", "sin", zero, ufloat(0.0, 0.1)),
+            ("sin", "sin", pi_4, ufloat(0.7071067811865476, 0.07071067811865477)),
+            ("sin", "sin", pi_2, ufloat(1.0, 6.123233995736766e-18)),
+            ("sinh", "sinh", zero, ufloat(0.0, 0.1)),
+            ("sinh", "sinh", pi_4, ufloat(0.8686709614860095, 0.1324609089252006)),
+            ("sinh", "sinh", pi_2, ufloat(2.3012989023072947, 0.2509178478658057)),
+            ("tan", "tan", zero, ufloat(0.0, 0.1)),
+            ("tan", "tan", pi_4, ufloat(0.9999999999999999, 0.19999999999999998)),
+            ("tan", "tan", pi_2, ufloat(1.633123935319537e16, 2.6670937881135717e31)),
+            ("tanh", "tanh", zero, ufloat(0.0, 0.1)),
+            ("tanh", "tanh", pi_4, ufloat(0.6557942026326724, 0.05699339637933774)),
+            ("tanh", "tanh", pi_2, ufloat(0.9171523356672744, 0.015883159318006324)),
+            ("arccos", "acos", zero, ufloat(1.5707963267948966, 0.1)),
+            ("arccos", "acos", one, ufloat(0.0, float("nan"))),
+            ("arccosh", "acosh", one, ufloat(0.0, float("nan"))),
+            ("arcsin", "asin", zero, ufloat(0.0, 0.1)),
+            ("arcsin", "asin", one, ufloat(1.5707963267948966, float("nan"))),
+            ("arcsinh", "asinh", zero, ufloat(0.0, 0.1)),
+            ("arcsinh", "asinh", one, ufloat(0.8813735870195429, 0.07071067811865475)),
+            ("arctan", "atan", zero, ufloat(0.0, 0.1)),
+            ("arctan", "atan", one, ufloat(0.7853981633974483, 0.05)),
+            ("arctanh", "atanh", zero, ufloat(0.0, 0.1)),
+            ("exp", "exp", zero, ufloat(1.0, 0.1)),
+            ("exp", "exp", one, ufloat(2.718281828459045, 0.27182818284590454)),
+            ("exp2", None, zero, ufloat(1.0, 0.06931471805599453)),
+            ("exp2", None, one, ufloat(2.0, 0.13862943611198905)),
+            ("expm1", "expm1", zero, ufloat(0.0, 0.1)),
+            ("expm1", "expm1", one, ufloat(1.718281828459045, 0.27182818284590454)),
+            ("log10", "log10", one, ufloat(0.0, 0.04342944819032518)),
+            ("log1p", "log1p", zero, ufloat(0.0, 0.1)),
+            ("log1p", "log1p", one, ufloat(0.6931471805599453, 0.05)),
+            ("degrees", "degrees", zero, ufloat(0.0, 5.729577951308233)),
+            ("degrees", "degrees", one, ufloat(57.29577951308232, 5.729577951308233)),
+            ("radians", "radians", zero, ufloat(0.0, 0.0017453292519943296)),
+            (
+                "radians",
+                "radians",
+                one,
+                ufloat(0.017453292519943295, 0.0017453292519943296),
+            ),
+            ("rad2deg", "degrees", zero, ufloat(0.0, 5.729577951308233)),
+            ("rad2deg", "degrees", one, ufloat(57.29577951308232, 5.729577951308233)),
+            ("deg2rad", "radians", zero, ufloat(0.0, 0.0017453292519943296)),
+            (
+                "deg2rad",
+                "radians",
+                one,
+                ufloat(0.017453292519943295, 0.0017453292519943296),
+            ),
+            ("sqrt", "sqrt", zero, ufloat(0.0, float("nan"))),
+            ("sqrt", "sqrt", one, ufloat(1.0, 0.05)),
         ],
     )
     def test_single_arg(self, numpy_func, umath_func, arg, expected):
@@ -236,4 +245,3 @@ def test_single_arg(self, numpy_func, umath_func, arg, expected):
             func = getattr(umath, umath_func)
             result = func(arg)
             assert nominal_and_std_dev_close(result, expected)
-        
\ No newline at end of file
diff --git a/uncertainties/core.py b/uncertainties/core.py
index 37436aec..9882e4b5 100644
--- a/uncertainties/core.py
+++ b/uncertainties/core.py
@@ -333,6 +333,7 @@ def __getstate__(self):
     def __setstate__(self, state):
         (self.linear_combo,) = state
 
+
 class AffineScalarFunc(UFloatNumpy):
     """
     Affine functions that support basic mathematical operations
diff --git a/uncertainties/ufloatnumpy.py b/uncertainties/ufloatnumpy.py
index 6c4b3321..4abd40f7 100644
--- a/uncertainties/ufloatnumpy.py
+++ b/uncertainties/ufloatnumpy.py
@@ -1,5 +1,6 @@
 import numpy as np
 import math
+
 # ufuncs are listed at https://numpy.org/doc/stable/reference/ufuncs.html
 from . import ops
 # from .umath_core import log_der0,_deriv_copysign, _deriv_fabs, _deriv_pow_0, _deriv_pow_1
@@ -14,50 +15,56 @@ def log_der0(*args):
     Works whether 1 or 2 arguments are given.
     """
     if len(args) == 1:
-        return 1/args[0]
+        return 1 / args[0]
     else:
-        return 1/args[0]/math.log(args[1])  # 2-argument form
+        return 1 / args[0] / math.log(args[1])  # 2-argument form
 
     # The following version goes about as fast:
 
     ## A 'try' is used for the most common case because it is fast when no
     ## exception is raised:
-    #try:
+    # try:
     #    return log_1arg_der(*args)  # Argument number check
-    #except TypeError:
+    # except TypeError:
     #    return 1/args[0]/math.log(args[1])  # 2-argument form
 
-def _deriv_copysign(x,y):
+
+def _deriv_copysign(x, y):
     if x >= 0:
         return math.copysign(1, y)
     else:
         return -math.copysign(1, y)
 
+
 def _deriv_fabs(x):
     if x >= 0:
         return 1
     else:
         return -1
 
+
 def _deriv_pow_0(x, y):
     if y == 0:
-        return  0.
+        return 0.0
     elif x != 0 or y % 1 == 0:
-        return y*math.pow(x, y-1)
+        return y * math.pow(x, y - 1)
     else:
-        return float('nan')
+        return float("nan")
+
 
 def _deriv_pow_1(x, y):
     if x == 0 and y > 0:
-        return 0.
+        return 0.0
     else:
         return math.log(x) * math.pow(x, y)
 
+
 def is_upcast_type(t):
     # This can be used to allow downstream modules to overide operations; see pint
     # TODO add upcast_type list or dict to a public interface
     return False
 
+
 def implements(numpy_func_string, func_type):
     """Register an __array_function__/__array_ufunc__ implementation for UArray
     objects.
@@ -75,7 +82,8 @@ def decorator(func):
 
     return decorator
 
-HANDLED_FUNCTIONS = {} # noqa
+
+HANDLED_FUNCTIONS = {}  # noqa
 HANDLED_UFUNCS = {}
 
 
@@ -102,6 +110,7 @@ def apply_func_elementwise(func, inputs, kwargs, result_dtype="object"):
         result = func(*inputs, **kwargs)
     return result
 
+
 def numpy_wrap(func_type, func, args, kwargs, types):
     """Return the result from a NumPy function/ufunc as wrapped by uncertainties."""
 
@@ -122,6 +131,7 @@ def numpy_wrap(func_type, func, args, kwargs, types):
         return NotImplemented
     return handled[name](*args, **kwargs)
 
+
 class UFloatNumpy(object):
     # NumPy function/ufunc support
     __array_priority__ = 17
@@ -170,53 +180,58 @@ def __array_function__(self, func, types, args, kwargs):
     #     for func_str, derivatives in ufunc_derivatives.items():
     #         implement_ufunc(func_str, derivatives)
 
-
     @classmethod
     def _add_numpy_arithmetic_ufuncs(cls):
         def implement_ufunc(func_str, derivatives):
             func = getattr(np, func_str)
+
             @implements(func_str, "ufunc")
             def implementation(*inputs, **kwargs):
                 return apply_func_elementwise(
-                    ops._wrap(cls, func, derivatives), inputs, kwargs)
+                    ops._wrap(cls, func, derivatives), inputs, kwargs
+                )
+
             return implementation
 
         ufunc_derivatives = {
-            'add': [lambda x, y: 1., lambda x, y: 1.],
-            'subtract': [lambda x, y: 1., lambda x, y: -1.],
-            'multiply': [lambda x, y: y, lambda x, y: x],
-            'divide': [lambda x, y: 1./y, lambda x, y: -x/y**2],
-            'true_divide': [lambda x, y: 1./y, lambda x, y: -x/y**2],
-            'floor_divide': [lambda x, y: 0., lambda x, y: 0.],
-            
-            'arccos': [nan_if_exception(lambda x: -1/math.sqrt(1-x**2))],
-            'arccosh': [nan_if_exception(lambda x: 1/math.sqrt(x**2-1))],
-            'arcsin': [nan_if_exception(lambda x: 1/math.sqrt(1-x**2))],
-            'arcsinh': [nan_if_exception(lambda x: 1/math.sqrt(1+x**2))],
-            'arctan': [nan_if_exception(lambda x: 1/(1+x**2))],
-            'arctan2': [nan_if_exception(lambda y, x: x/(x**2+y**2)),  # Correct for x == 0
-                        nan_if_exception(lambda y, x: -y/(x**2+y**2))],  # Correct for x == 0
-            'arctanh': [nan_if_exception(lambda x: 1/(1-x**2))],
-            'cos': [lambda x: -math.sin(x)],
-            'cosh': [math.sinh],
-            'sin': [math.cos],
-            'sinh': [math.cosh],
-            'tan': [nan_if_exception(lambda x: 1+math.tan(x)**2)],
-            'tanh': [nan_if_exception(lambda x: 1-math.tanh(x)**2)],
-            'exp': [math.exp],
+            "add": [lambda x, y: 1.0, lambda x, y: 1.0],
+            "subtract": [lambda x, y: 1.0, lambda x, y: -1.0],
+            "multiply": [lambda x, y: y, lambda x, y: x],
+            "divide": [lambda x, y: 1.0 / y, lambda x, y: -x / y**2],
+            "true_divide": [lambda x, y: 1.0 / y, lambda x, y: -x / y**2],
+            "floor_divide": [lambda x, y: 0.0, lambda x, y: 0.0],
+            "arccos": [nan_if_exception(lambda x: -1 / math.sqrt(1 - x**2))],
+            "arccosh": [nan_if_exception(lambda x: 1 / math.sqrt(x**2 - 1))],
+            "arcsin": [nan_if_exception(lambda x: 1 / math.sqrt(1 - x**2))],
+            "arcsinh": [nan_if_exception(lambda x: 1 / math.sqrt(1 + x**2))],
+            "arctan": [nan_if_exception(lambda x: 1 / (1 + x**2))],
+            "arctan2": [
+                nan_if_exception(lambda y, x: x / (x**2 + y**2)),  # Correct for x == 0
+                nan_if_exception(lambda y, x: -y / (x**2 + y**2)),
+            ],  # Correct for x == 0
+            "arctanh": [nan_if_exception(lambda x: 1 / (1 - x**2))],
+            "cos": [lambda x: -math.sin(x)],
+            "cosh": [math.sinh],
+            "sin": [math.cos],
+            "sinh": [math.cosh],
+            "tan": [nan_if_exception(lambda x: 1 + math.tan(x) ** 2)],
+            "tanh": [nan_if_exception(lambda x: 1 - math.tanh(x) ** 2)],
+            "exp": [math.exp],
             "exp2": [lambda y: _deriv_pow_1(2, y)],
-            'expm1': [math.exp],
-            'log10': [nan_if_exception(lambda x: 1/x/math.log(10))],
-            'log1p': [nan_if_exception(lambda x: 1/(1+x))],
-            'degrees': [lambda x: math.degrees(1)],
-            'rad2deg': [lambda x: math.degrees(1)],
-            'radians': [lambda x: math.radians(1)],
-            'deg2rad': [lambda x: math.radians(1)],
-            'power': [_deriv_pow_0, _deriv_pow_1],
-            'sqrt': [nan_if_exception(lambda x: 0.5/math.sqrt(x))],
-            'hypot': [lambda x, y: x/math.hypot(x, y),
-                      lambda x, y: y/math.hypot(x, y)],
-            }
+            "expm1": [math.exp],
+            "log10": [nan_if_exception(lambda x: 1 / x / math.log(10))],
+            "log1p": [nan_if_exception(lambda x: 1 / (1 + x))],
+            "degrees": [lambda x: math.degrees(1)],
+            "rad2deg": [lambda x: math.degrees(1)],
+            "radians": [lambda x: math.radians(1)],
+            "deg2rad": [lambda x: math.radians(1)],
+            "power": [_deriv_pow_0, _deriv_pow_1],
+            "sqrt": [nan_if_exception(lambda x: 0.5 / math.sqrt(x))],
+            "hypot": [
+                lambda x, y: x / math.hypot(x, y),
+                lambda x, y: y / math.hypot(x, y),
+            ],
+        }
         # TODO: test arctan2, power, hypot
         for func_str, derivatives in ufunc_derivatives.items():
             implement_ufunc(func_str, derivatives)
@@ -240,12 +255,12 @@ def implementation(*inputs, **kwargs):
             return implementation
 
         ufunc_comparatives = {
-            'equal': ops.eq_on_aff_funcs,
-            'not_equal': ops.ne_on_aff_funcs,
-            'less': ops.lt_on_aff_funcs,
-            'less_equal': ops.le_on_aff_funcs,
-            'greater': ops.gt_on_aff_funcs,
-            'greater_equal': ops.ge_on_aff_funcs,
+            "equal": ops.eq_on_aff_funcs,
+            "not_equal": ops.ne_on_aff_funcs,
+            "less": ops.lt_on_aff_funcs,
+            "less_equal": ops.le_on_aff_funcs,
+            "greater": ops.gt_on_aff_funcs,
+            "greater_equal": ops.ge_on_aff_funcs,
         }
         for func_str, func in ufunc_comparatives.items():
             implement_ufunc(func_str, func)
@@ -291,5 +306,3 @@ def implementation(*inputs, **kwargs):
     #         "rad2deg": ("radian", "degree"),
     #         "logaddexp": ("", ""),
     #         "logaddexp2": ("", ""),
-
-            
\ No newline at end of file

From e928cbbd09a9b479480272761ea51b1fe591e33c Mon Sep 17 00:00:00 2001
From: Andrew <andrewgsavage@gmail.com>
Date: Sat, 13 Jul 2024 09:41:29 +0100
Subject: [PATCH 6/7] make AffineScalarFunc not require np

---
 uncertainties/core.py        | 10 ++++-
 uncertainties/ufloatnumpy.py | 73 ------------------------------------
 2 files changed, 9 insertions(+), 74 deletions(-)

diff --git a/uncertainties/core.py b/uncertainties/core.py
index 9882e4b5..82e9e849 100644
--- a/uncertainties/core.py
+++ b/uncertainties/core.py
@@ -334,7 +334,15 @@ def __setstate__(self, state):
         (self.linear_combo,) = state
 
 
-class AffineScalarFunc(UFloatNumpy):
+try:
+    import numpy
+except ImportError:
+    parent_classes = []
+else:
+    parent_classes = [UFloatNumpy]
+
+
+class AffineScalarFunc(*parent_classes):
     """
     Affine functions that support basic mathematical operations
     (addition, etc.).  Such functions can for instance be used for
diff --git a/uncertainties/ufloatnumpy.py b/uncertainties/ufloatnumpy.py
index 4abd40f7..3fcf16b3 100644
--- a/uncertainties/ufloatnumpy.py
+++ b/uncertainties/ufloatnumpy.py
@@ -83,7 +83,6 @@ def decorator(func):
     return decorator
 
 
-HANDLED_FUNCTIONS = {}  # noqa
 HANDLED_UFUNCS = {}
 
 
@@ -126,8 +125,6 @@ def numpy_wrap(func_type, func, args, kwargs, types):
         raise ValueError(f"Invalid func_type {func_type}")
 
     if name not in handled or any(is_upcast_type(t) for t in types):
-        print("NotImplemented L54")
-        raise TypeError
         return NotImplemented
     return handled[name](*args, **kwargs)
 
@@ -153,33 +150,6 @@ def __array_ufunc__(self, ufunc, method, *inputs, **kwargs):
     def __array_function__(self, func, types, args, kwargs):
         return numpy_wrap("function", func, args, kwargs, types)
 
-    # original code for _add_numpy_ufuncs. may be helpful for writing a generic wraps
-    # this can be deleted:
-    # @classmethod
-    # def _add_numpy_ufuncs(cls):
-    #     def implement_ufunc(func_str, derivatives):
-    #         func = getattr(np, func_str)
-    #         @implements(func_str, "ufunc")
-    #         def implementation(*inputs, **kwargs):
-    #             if isinstance(inputs[0], np.ndarray):
-    #                 result = np.empty_like(inputs[0], dtype=object)
-    #                 for index, x in np.ndenumerate(inputs[0]):
-    #                     inputs_ = (x if i == 0 else inputs[i] for i in range(len(inputs)))
-    #                     result[index] = cls.wrap(func, derivatives)(*inputs_, **kwargs)
-    #             elif isinstance(inputs[1], np.ndarray):
-    #                 result = np.empty_like(inputs[1], dtype=object)
-    #                 for index, x in np.ndenumerate(inputs[1]):
-    #                     inputs_ = (x if i == 1 else inputs[i] for i in range(len(inputs)))
-    #                     result[index] = cls.wrap(func, derivatives)(*inputs_, **kwargs)
-    #             else:
-    #                 result = cls.wrap(func, derivatives)(*inputs, **kwargs)
-    #             return result
-
-    #         return implementation
-
-    #     for func_str, derivatives in ufunc_derivatives.items():
-    #         implement_ufunc(func_str, derivatives)
-
     @classmethod
     def _add_numpy_arithmetic_ufuncs(cls):
         def implement_ufunc(func_str, derivatives):
@@ -239,7 +209,6 @@ def implementation(*inputs, **kwargs):
     @classmethod
     def _add_numpy_comparative_ufuncs(cls):
         def recursive_to_affine_scalar(arr):
-            print(arr)
             if isinstance(arr, (list, tuple)):
                 return type(arr)([recursive_to_affine_scalar(i) for i in arr])
             if isinstance(arr, np.ndarray):
@@ -264,45 +233,3 @@ def implementation(*inputs, **kwargs):
         }
         for func_str, func in ufunc_comparatives.items():
             implement_ufunc(func_str, func)
-
-    # 'copysign': [_deriv_copysign,
-    #              lambda x, y: 0],
-    # 'degrees': [lambda x: math.degrees(1)],
-    # 'erf': [lambda x: math.exp(-x**2)*erf_coef],
-    # 'erfc': [lambda x: -math.exp(-x**2)*erf_coef],
-    # 'fabs': [_deriv_fabs],
-    # 'hypot': [lambda x, y: x/math.hypot(x, y),
-    #           lambda x, y: y/math.hypot(x, y)],
-    # 'log': [log_der0,
-    #         lambda x, y: -math.log(x, y)/y/math.log(y)],
-    # 'pow': [_deriv_pow_0, _deriv_pow_1],
-    # 'radians': [lambda x: math.radians(1)],
-    # 'sqrt': [lambda x: 0.5/math.sqrt(x)],
-    # def _add_trig_ufuncs(cls):
-
-    #         "cumprod": ("", ""),
-    #         "arccos": ("", "radian"),
-    #         "arcsin": ("", "radian"),
-    #         "arctan": ("", "radian"),
-    #         "arccosh": ("", "radian"),
-    #         "arcsinh": ("", "radian"),
-    #         "arctanh": ("", "radian"),
-    #         "exp": ("", ""),
-    #         "expm1": ("", ""),
-    #         "exp2": ("", ""),
-    #         "log": ("", ""),
-    #         "log10": ("", ""),
-    #         "log1p": ("", ""),
-    #         "log2": ("", ""),
-    #         "sin": ("radian", ""),
-    #         "cos": ("radian", ""),
-    #         "tan": ("radian", ""),
-    #         "sinh": ("radian", ""),
-    #         "cosh": ("radian", ""),
-    #         "tanh": ("radian", ""),
-    #         "radians": ("degree", "radian"),
-    #         "degrees": ("radian", "degree"),
-    #         "deg2rad": ("degree", "radian"),
-    #         "rad2deg": ("radian", "degree"),
-    #         "logaddexp": ("", ""),
-    #         "logaddexp2": ("", ""),

From f45193ac38b09a2e5d44b616f375433ace38b3ea Mon Sep 17 00:00:00 2001
From: Andrew <andrewgsavage@gmail.com>
Date: Sat, 13 Jul 2024 09:43:50 +0100
Subject: [PATCH 7/7] make AffineScalarFunc not require np

---
 uncertainties/ufloatnumpy.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/uncertainties/ufloatnumpy.py b/uncertainties/ufloatnumpy.py
index 3fcf16b3..5373d535 100644
--- a/uncertainties/ufloatnumpy.py
+++ b/uncertainties/ufloatnumpy.py
@@ -84,6 +84,7 @@ def decorator(func):
 
 
 HANDLED_UFUNCS = {}
+HANDLED_FUNCTIONS = {}
 
 
 def apply_func_elementwise(func, inputs, kwargs, result_dtype="object"):