AffineScalarFunction and related refactoring.

[jagerber48]

There have been many discussions targeting a refactoring of AffineScalarFunction and the related Variable, LinearCombination etc.

Here's a bit of code I wrote that re-implements a lot of the core uncertainties functionality but in ways which I feel are more clear.

Here I laid out some goals for this refactoring

• Resolve a bug or two we discovered in how correlations between variables were tracked
• Make the code dramatically easier to understand (On one level, this just appeals to a desire I have for code to be clean/readable. On a more practical level readable code improves maintainability and improvability).
• Ensure there is a strategy to save uncertainties objects in such a way that their correlations can be saved and reloaded
• Ensure performance is not heavily impacted by any changes we make (maybe this one ought to be higher in my list)

from __future__ import annotations

from collections import defaultdict
from copy import deepcopy
from dataclasses import dataclass, field
from functools import cache, wraps
import inspect
from math import sqrt
import sys
from typing import Callable, Collection, Tuple, Union, TYPE_CHECKING
import uuid

if TYPE_CHECKING:
    from inspect import Signature


@dataclass(frozen=True)
class UncertaintyAtom:
    """
    Custom class to keep track of "atoms" of uncertainty. Note e.g. that
    UncertaintyAtom(3) is UncertaintyAtom(3)
    returns False.
    """
    unc: float
    uuid: uuid.uuid4 = field(default_factory=uuid.uuid4, init=False)

UncertaintyCombo = Tuple[
    Tuple[
        Union[UncertaintyAtom, "UncertaintyCombo"],
        float
    ],
    ...
]
UncertaintyComboExpanded = Tuple[
    Tuple[
        UncertaintyAtom,
        float
    ],
    ...
]


@cache
def get_expanded_combo(combo: UncertaintyCombo) -> UncertaintyComboExpanded:
    """
    Recursively expand a linear combination of uncertainties out into the base Atoms.
    It is a performance optimization to sometimes store unexpanded linear combinations.
    For example, there may be a long calculation involving many layers of UFloat
    manipulations. We need not expand the linear combination until the end when a
    calculation of a standard deviation on a UFloat is requested.
    """
    expanded_dict = defaultdict(float)
    for unc_combo_1, weight_1 in combo:
        if isinstance(unc_combo_1, UncertaintyAtom):
            expanded_dict[unc_combo_1] += weight_1
        else:
            expanded_sub_combo = get_expanded_combo(unc_combo_1)
            for unc_atom, weight_2 in expanded_sub_combo:
                expanded_dict[unc_atom] += weight_2 * weight_1

    return tuple((unc_atom, weight) for unc_atom, weight in expanded_dict.items())


Value = Union["UValue", float]


class UFloat:
    """
    Core class. Stores a mean value (val, nominal_value, n) and an uncertainty stored
    as a (possibly unexpanded) linear combination of uncertainty atoms. Two UFloat's
    which share non-zero weight for a certain uncertainty atom are correlated.

    UFloats can be combined using arithmetic and more sophisticated mathematical
    operations. The uncertainty is propagation using rules of linear uncertainty
    propagation.
    """
    def __init__(
            self,
            /,
            val: float,
            unc: Union[UncertaintyCombo, float] = ()):
        self._val = float(val)
        if isinstance(unc, (float, int)):
            unc_atom = UncertaintyAtom(unc)
            unc_combo = ((unc_atom, 1.0),)
            self.unc_linear_combo = unc_combo
        else:
            self.unc_linear_combo = unc

    @property
    def val(self: "UFloat") -> float:
        return self._val

    @property
    def unc(self: "UFloat") -> float:
        expanded_combo = get_expanded_combo(self.unc_linear_combo)
        return sqrt(sum([(weight * unc_atom.unc)**2 for unc_atom, weight in expanded_combo]))

    @property
    def nominal_value(self: "UFloat") -> float:
        return self.val

    @property
    def n(self: "UFloat") -> float:
        return self.val

    @property
    def std_dev(self: "UFloat") -> float:
        return self.unc

    @property
    def s(self: "UFloat") -> float:
        return self.unc

    def __repr__(self) -> str:
        return f'{self.__class__.__name__}({self.val}, {self.unc})'


STEP_SIZE = sqrt(sys.float_info.epsilon)


def get_param_name(sig: Signature, param: Union[int, str]):
    if isinstance(param, int):
        param_name = list(sig.parameters.keys())[param]
    else:
        param_name = param
    return param_name


def partial_derivative(
        f: Callable[..., float],
        target_param: Union[str, int],
        *args,
        **kwargs
):
    """
    Calculate the partial derivative of a function f with respect to the target_param
    (string name of the float variable to f to be varied) holding all other arguments of
    *args and **kwargs constant.
    """
    sig = inspect.signature(f)
    lower_bound_sig = sig.bind(*args, **kwargs)

    for arg, val in lower_bound_sig.arguments.items():
        if isinstance(val, UFloat):
            lower_bound_sig.arguments[arg] = val.val

    upper_bound_sig = deepcopy(lower_bound_sig)

    target_param_name = get_param_name(sig, target_param)

    x = lower_bound_sig.arguments[target_param_name]
    dx = abs(x) * STEP_SIZE  # Step size scaled by scale of x
    # dx = STEP_SIZE  # Or keep step size fixed?

    # Inject x - dx into target_param and evaluate f
    lower_bound_sig.arguments[target_param_name] = x - dx
    lower_y = f(*lower_bound_sig.args, **lower_bound_sig.kwargs)

    # Inject x + dx into target_param and evaluate f
    upper_bound_sig.arguments[target_param_name] = x + dx
    upper_y = f(*upper_bound_sig.args, **upper_bound_sig.kwargs)

    derivative = (upper_y - lower_y) / (2 * dx)
    return derivative


class ToUFunc:
    """
    Decorator to convert a function which typically accepts float inputs into a function
    which accepts UFloat inputs.

    e.g.
    > @ToUFunc(('x', 'y'))
    > def multiply(x, y, print_str='print this string!', do_print=False):
    >     if do_print:
    >         print(print_str)
    >     return x * y

    Pass in a list of parameter names which correspond to float inputs that should now
    accept UFloat inputs.

    To calculate the output nominal value the decorator replaces all float inputs with
    their respective nominal values and evaluates the function directly.

    To calculate the output uncertainty linear combination the decorator calculates the
    partial derivative of the function with respect to each UFloat entry and appends the
    uncertainty linear combination corresponding to that UFloat, weighted by the
    corresponding partial derivative.
    """
    def __init__(self, u_float_params: Collection[Union[str, int]]):
        self.u_float_params = u_float_params

    def __call__(self, f: Callable[..., float]):
        sig = inspect.signature(f)

        @wraps(f)
        def wrapped(*args, **kwargs):
            """
            Calculate the
            """
            unc_linear_combo = []
            bound = sig.bind(*args, **kwargs)
            return_u_val = False

            for u_float_param in self.u_float_params:
                target_param_name = get_param_name(sig, u_float_param)
                arg = bound.arguments[target_param_name]
                if isinstance(arg, UFloat):
                    return_u_val = True
                    bound.arguments[target_param_name] = arg.val
                    sub_unc_linear_combo = arg.unc_linear_combo
                    derivative = partial_derivative(f, target_param_name, *args, **kwargs)
                    unc_linear_combo.append((sub_unc_linear_combo, derivative))

            new_val = f(*bound.args, **bound.kwargs)

            if not return_u_val:
                return new_val

            unc_linear_combo = tuple(unc_linear_combo)
            return UFloat(new_val, unc_linear_combo)

        return wrapped


def add_float_funcs_to_uvalue():
    """
    Monkey-patch common float instance methods over to UFloat
    """
    float_funcs_dict = {
        '__abs__': (0,),
        '__add__': (0, 1),
        '__radd__': (0, 1),
        '__sub__': (0, 1),
        '__rsub__': (0, 1),
        '__mul__': (0, 1),
        '__rmul__': (0, 1),
        '__truediv__': (0, 1),
        '__rtruediv__': (0, 1),
        '__pow__': (0, 1),
        '__rpow__': (0, 1),
    }

    for func_name, u_val_params in float_funcs_dict.items():
        float_func = getattr(float, func_name)
        ufloat_ufunc = ToUFunc(u_val_params)(float_func)
        setattr(UFloat, func_name, ufloat_ufunc)


add_float_funcs_to_uvalue()

Examples:

from math import sin

usin = ToUFunc((0,))(sin)

x = UFloat(10, 1)

y = UFloat(10, 1)

z = UFloat(20, 2)

print(f'{x=}')
print(f'{3*x=}')
print(f'{x-x=}  # A UFloat is correlated with itself')

print(f'{y=}')
print(f'{x-y=}  # Two distinct UFloats are not correlated unless they have the same Uncertainty Atoms')

print(f'{z=}')

print(f'{x*z=}')
print(f'{x/z=}')
print(f'{x**z=}')

print(f'{usin(x)=}  # We can UFloat-ify complex functions')

# x=UFloat(10.0, 1.0)
# 3*x=UFloat(30.0, 3.0)
# x-x=UFloat(0.0, 0.0)  # A UFloat is correlated with itself
# y=UFloat(10.0, 1.0)
# x-y=UFloat(0.0, 1.4142135623730951)  # Two distinct UFloats are not correlated unless they have the same Uncertainty Atoms
# z=UFloat(20.0, 2.0)
# x*z=UFloat(200.0, 28.284271247461902)
# x/z=UFloat(0.5, 0.07071067811865477)
# x**z=UFloat(1e+20, 5.0207163276303525e+20)
# usin(x)=UFloat(-0.5440211108893698, 0.8390715289860964)  # We can UFloat-ify complex functions

I don't know if this code resolves the stated bugs. I need to go back through the discussion to find what those bugs are.

---

I think this code is much more understandable/readable than the current source code. One way it is greatly helped is by using the newer inspect.signature methods for function signature introspection and manipulation. This is usually much easier to work with and understand than the old argspec approaches. But I also think some of the other technical details (linear combination expansion, partial derivative calculation, uncertainty propagation) are also more clear.

Importantly this code also gets around some of the circularity that seems to be present in the existing source code via the introduction of the simple UncertaintyAtom dataclass. There is one point of self-referencing, but that only appears in the type alias for UncertaintyLinearCombination.

---

I think this code will be easier to make saving/loading work. The only sticking point I see is that it is possible there could be two UFloats that are equal, meaning they have the same value and their uncertainty linear combinations expand to be equal, but the unexpanded forms of their linear combinations differ from each other. The question is, should two such UFloats be equal to each other? I think yes. Should they have the same hash? I'm not sure.

It's also pretty clear in this code that all the major classes could be made/interpreted to be immutable.

---

I don't know how this code compares in performance to the current implementation. I tried to follow a similar strategy so I hope it's not too far off. It does look like there was some caching that might be implemented in the current source code which I didn't implement here.

Originally I was considering schemes that looked something like:

• Give a UID to each variable.
• If any variable X is correlated with other variables Y, Z then store the UIDs Y, Z in X along with a covariance matrix capturing the correlations. @lebigot had a comment elsewhere that this would be inefficient and there are comments in the code to a similar effect.

So to this end I tried to stick with the "linear combination of uncertainty atoms" approach. I hope this means performance is similar.

---

Obviously this code doesn't capture the full functionality of uncertainties. But I do think it captures a lot of the core functionality in many fewer lines and in a way that is more readable. I'm curious what bits of this implementation might be ported into the main source code

---

Edit: Updated code so that you can pass either function parameter names that should allow UFloat, or parameter numbers. The latter is very useful when dealing with built-in math functions that typically expose positional arguments whose actual names are hard to track down. The latter may be useful for users to UFloat-ify their own math functions.

---

Edit2: It was necessary to add a uuid field to UncertaintyAtom to ensure UFloat(1, 1) - UFloat(1, 1) did NOT give UFloat(0, 0). Not that if x=UFloat(1, 0) then x-x DOES give UFloat(0, 0).

[newville]

@jagerber48 I like this in general. I might have a few small nits:
1. somehow 'unc' really does not work for me, could that be 'std_dev'?;
2. I really do not like black formatting. I can work around it, but refuse to contribute to projects that require it (I'm a privileged
old man, and I can say no to sh*t, and black is pure evil).

But, also: I think this is really interesting.

In

    float_funcs_dict = { '__abs__': (0,),   '__add__': (0, 1), ...}

and the setting of (and use of) u_float_params, could that just be

    float_funcs_dict = { '__abs__': 1,   '__add__': 2, ...}

to indicate how many positional arguments there were? I think that will always be (0, 1, ..., n-1), so why not just say n.
You could probably just have a list of methods, and use len(inspect.signature(f)._parameters), but then to make that more general (beyond "float methods"), you might want to have to check for which are positional only and which are keyword arguments.

And also:

could "float" be replaced with "float or complex or ndarray"? Where is "float and Python math" required here, and could that be expanded to "available numeric representations"?

[jagerber48]

@newville thanks for the feedback.

1. Yes, can do unc -> std_dev no problem. I like the word uncertainty, because that is semantically what the parameter is, and I like that unc is fewer keystrokes than std_dev. But it is also semantically the standard deviation of the random variable. So yes, like I said, no problem to do unc -> std_dev.
2. I didn't run a formatter on this :p it's just the formatting I've gotten used to. It's according to the rules of PyCharm and Ruff. Depending on which bit of formatting you don't like I may or may not have strong opinions about it, but this is a side-topic.

---

Yeah, the float_funcs_dict is pretty confusing in this code. I'm working this code into a branch on my fork of uncertainties and I'm rehashing this bit pretty heavily. https://github.com/jagerber48/uncertainties/blob/major/core_rewrite/uncertainties/core_new.py The equivalent code in master is also a bit complicated. But, in short, in cases where we know the function being replaced uses positional arguments, yes, we can simplify the specification of the custom derivative function specification. The way I've implemented it so far is I've kept the implementation above for the general case: the users supplies a list of ufloat params, specified by number or name, together with a dict that optionally provides custom derivative calculation functions. But then I've included a helper decorator as another layer of abstraction for use in the "we're assuming positional arguments" case that is what we need most of the time. In this case the user passes a list of functions or functions-as-strings, like -x/y**2.

That said, I'm tempted to do away with all of this (confusing, imo) monkey-patching business and just directly define the dunder methods (__mul__, etc.) on the UFloat class. It will be easier to read, even though it's more code, and it will give code inspection tools an easier time. We would have something like

class UFloat:
    ...
    @ToUFuncPositional(('y', 'x'))
    def __mul__(self, other):
        return self * other

    @ToUFuncPositional(('x', 'y'))
    def __rmul__(self, other):
        return self * other

or maybe

class UFloat:
    ...
    def __mul__(self: "UFloat", other):
        return ToUFuncPositional(('x', 'y'))(float.__mul__)(self, other)

    def __rmul__(self: "UFloat", other):
        return ToUFuncPositional(('y', 'x'))(float.__rmul__)(self, other)

---

could "float" be replaced with "float or complex or ndarray"? Where is "float and Python math" required here, and could that be expanded to "available numeric representations"?

For extending to complex I'd need to think a little bit about what it means to have a complex random variable. I worry the partial derivative calculations won't be straightforward..

For ndarray I'm not totally clear what the goal what the goal would be. There nothing that is specifically requiring float or python math. I'm researching #253 trying to understand how to get a numpy array that follows the rules of a new scalar type but having trouble coming to a conclusion.

[andrewgsavage]

I'm finding it harder than it needs to be to follow since you've used new class and attributes names. It'd be easier to compare current behaviour to this if less has changed.

This looks like a good future state. It'd be good to see how you'd handle analytical derivatives. I don't think numerical derivatives should be used when the current behaviour uses analytical for many functions. I found the wrapping to return an AffineScalarType and dealing with nans felt more difficult than it should be.

A partial derivative decorator that can handle ndarrays would be very helpful for numpy support. Matrix multiplication is an example of a function that'd be time consuming to evaluate numerically when the sizes of matrices gets big.

Recursively expand a linear combination of uncertainties out into the base Atoms.
It is a performance optimization to sometimes store unexpanded linear combinations.
For example, there may be a long calculation involving many layers of UFloat
manipulations. We need not expand the linear combination until the end when a
calculation of a standard deviation on a UFloat is requested.

I still question the value of delaying the expansion. This doesn't save any time unless the user does not need the std_dev - in which case, why are they using uncertainties?

[andrewgsavage]

There's quite a lot of changes in here so I'm not sure how much can be done alone. Is it possible to make a PR that only changes the linear combinations, or does this change how operations need coding?

[jagerber48]

I'm finding it harder than it needs to be to follow since you've used new class and attributes names. It'd be easier to compare current behaviour to this if less has changed.

Here is where I'm working right now: https://github.com/jagerber48/uncertainties/blob/major/core_rewrite/uncertainties/core_new.py
This is a total rewrite. Initially it was for my own sake of understanding how the current code works. I think it may help others understand also. I open it as a discussion rather than PR because I'm not, yet, proposing these changes actually be pulled into the source. But I am interested to discuss this.

UFloat is the main class here. Variable and AffineScalarFunc have been collapsed into UFloat. There is no LinearCombination class, rather the uncertainty_lin_combo attribute on UFloat is a tuple of pairs of uncertainties and weights. Each "uncertainty" can be an UncertaintyAtom (like a single dx in a calculation that has a value for std_dev and a unique identifier) or itself a linear combination of uncertainties.

It'd be good to see how you'd handle analytical derivatives. I don't think numerical derivatives should be used when the current behaviour uses analytical for many functions. I found the wrapping to return an AffineScalarType and dealing with nans felt more difficult than it should be.

Oh I thought the code above did handle analytical derivatives. I've been updating this code, largely to handle analytic derivatives. I'm curious to know how much of a speedup it is to use analytic compared to numerical derivatives. I guess it's pretty large but it would be interesting to know how much. See https://github.com/jagerber48/uncertainties/blob/major/core_rewrite/uncertainties/core_new.py#L392

I think the wrapping here is a little easier to understand that what is going on in ops.py now but there is a bit of technical machinery/abstraction needed that trades off readability for easy human readable specifications of analytical partial derivative functions (like passing "-x/y**2" instead of lambda x, y: -x/y**2).

A partial derivative decorator that can handle ndarrays would be very helpful for numpy support. Matrix multiplication is an example of a function that'd be time consuming to evaluate numerically when the sizes of matrices gets big.

Unclear what exactly is desired here. Can't all partial derivative calculations be deferred to the scalar operations? But anyways, it would be pretty easy to extend the partial derivative calculation to handle variables that are numpy arrays I think. But anyways, you have ideas about numpy arrays that I'm not up to speed yet on. Hopefully I can get up to speed in the numpy thread then this request will make more sense to me.

I still question the value of delaying the expansion. This doesn't save any time unless the user does not need the std_dev - in which case, why are they using uncertainties?

I think the theory is that it saves time if the user is calculating a lot of intermediate UFloats but is only interested in the std_dev of the final result. Whether that use case occurs regularly in a performance-demanding way, and whether this optimization actually provides a speed up is unknown (to me). But maybe Lebigot did some testing in the past?

There's quite a lot of changes in here so I'm not sure how much can be done alone. Is it possible to make a PR that only changes the linear combinations, or does this change how operations need coding?

My goals in such a PR would be to

1. remove some of the self-referential issues in the current implementation.
2. after having done 1. I think it would be more clear how to make AffineScalarFunc immutable and hashable
3. Make UFloat the real class and AffineScalarFunc the alias
4. Totally remove Variable

The change to linear combination would probably require at least a minor rework of how operations need to get coded. Specifically the bit where derivatives are being calculated and set as weights of new linear combinations. Rather than continue reverse-engineering how that works maybe it would be easier to just take the new implementation which is being forward engineered?

[jagerber48]

Here's a PR on my fork to more easily track the changes. I also made an actual test script to more clearly demonstrate implemented functionality:

PR: #256
tests: https://github.com/lmfit/uncertainties/pull/256/files#diff-52cc195c675039d1dd06f95855eb13ba74b5644cfc3bf5fc90f847e63ff9c040

[jagerber48]

Closed PR on this main repo (that was a mistake) and re-opened a PR on my fork: jagerber/uncertainties/2.

• More filling out and cleanup of core_new.
• Added umath_new
• Some test cases

[jagerber48]

refactored into a new sub-package and new tests to keep code separate. The uncertainties.new sub-package has all of my new code at the moment.

• ucombo.py: Has the UAtom and UCombo class which are how uncertainties are represented. UAtom is an independent unit of uncertainty with some std_dev. UCombo represents nested linear combinations of UAtom. std_dev can be calculated on UCombo. Some caching is done to speed up expansion of combos and calculations of std_dev. These are all frozen dataclasses so they're immutable and hashable.
• ufloat.py: Holds the UFloat class which keeps a float value together with an uncertainty UCombo. UFloat class will be monkey patched to include math functions with error propagation. Because of an unfortunate quirk of dataclass and typing it wasn't possible to make UFloat a frozen dataclass. Instead __hash__ is implemented directly and immutability is accomplished with properties. This is plenty good for our purposes.
• func_conversion.py: This is sort of the "core" file in a way where the action is happening. This module exposes the function wrappers that implement uncertainty propagation. Supports numerical partial derivates where necessary, or "analytic" partial derivatives where possible.
• umath.py: Here is the repository of functions which become error-propagation aware. There are some standard float operations like __mul__ that get monkey patched into UFloat. There are a number of functions in math that get monkey patched into umath and finally there are a number of mathematical functions that support the numpy ufunc protocol which also get monkey patched into UFloat. The monkey patching functions live in this module and are imported into __init__ where they are executed at import time.
• UArray holds a thin ndarray wrapper that mostly allows a few convenience methods you'd want on an array of UFloat. It provides this convenience while providing all the power of numpy for UFloat arrays.

Some effort has been made to give slick __repr__ and __str__ to UFloat, UCombo and UAtom, more work can be done.

Some notes/thoughts I've had during this re-implementation process:

• I don't think >, >=, <, <= should be implemented on UFloat. As explained in the docs, they don't satisfy trichotomy because == relies on the value and uncertainty whereas the comparison operators only depend on the value. I say if users want to compare nominal values they can just do unum.n rather than borking the UFloat data model. If folks want to discuss this point I say we open a dedicated discussion to it, though it's likely controversial, it's tangential to the overall rewrite.
• I've left a number of math operations unpopulated in the new code so far. It would be straightforward to implement them but I wonder, "is it worth it"? The functions I've left off are things having to do with finer detail of float arithmetic. I'm very much based in the physical sciences so I think: would someone in physical sciences need this computer-science oriented function. Is it worth implementing everything for exhaustiveness? Or can we just throw in what we think people will need and add more esoteric things as they're requested?
• With new.ufloat, new.umath, and new.uarray I think I've implemented a lot of the existing uncertainties functionality in a way that is, in my opinion, much easier to follow. Probably many fewer lines of source code. The organization is more modular.
• Notably missing is formatting, but this could be included rapidly since the formatting code is totally orthogonal to the error propagation code.
• Also missing are some things like hooks for copy or pickle but these could also be added easily.
• A hash function on UFloat is included helping resolve Should `Variable.__eq__` be based on `id()` or a separate identifier? #218 and Implement hash for AffineScalarFunc #219.
• Functions to generate a covariance matrix from a sequence of UFloat and to generate a sequence of UFloat from a covariance matrix and sequence of floats are not yet implemented. Though the implementation should be straightforward and can likely address Use cholesky decomposition in correlated_values and correlated_values_norm #103.
• I would argue the numpy support is a solid candidate worth consideration for resolving implement numpy ufuncs #249 and Numpy support #253
• I've also include type "stubs" for the float functions in UFloat and the math functions in umath to help type checkers and code inspectors/completers to be more useful tools despite the monkey patching in use. Type "stubs" haven't been included for the numpy ufuncs on UFloat yet, but they could be easily (just a lot of lines of boilerplate). Perhaps a proper type stub file would be a better way to go? I don't have experience with this.

There are likely many more detailed features I'm missing. Continuing on this path next steps would be:

• Get agreement on what features need to be covered
• Port tests from the old implementation to the new implementation. This will be the easiest way to confirm functionality is preserved.
• Time performance benchmarking?
• This is a total rewrite of the existing code. This would definitely be the major version bump 4.0.0. I would say if we're going to pay backwards incompatible changes (like removing comparison operators) to buy improved usability or maintainability now would be a good time to do it.
• This is a total rewrite of the existing code. I realize this might rub some people the wrong way. Not trying to offend anyone. Obviously all credit goes to Lebigot for the architecture and inspiration for this rewrite. The new code does the same stuff just using the tools and conventions of much more modern python than when uncertainties was first authored. I also hope I'm not stepping on toes of any large swaths of work that has gone on recently since the new maintainers stepped in. Probably implement numpy ufuncs #249 is the biggest risk here, so apologies for there. That PR also inspired a lot of research and the implementation I came up with for UArray. I do think there's an argument that extending this code to cover uncertainties functionality may be easier than porting improvements in bit-by-bit into the existing source code. Especially for stuff like ops.py or core.py. Besides all of this, it was fun for me to do and a nice learning experience.
• The MAJOR DANGER with new code like this is that it's not battle tested for years compared to old code. Adoption will require extensive testing. >= the amount of existing testing.

[jagerber48]

I think a better approach would be to work with the heavily/fully re-written branch but to get all of the old tests (and more) to work with minimal modification and then move forward guided by whatever diff ends up being required on the tests.

Yea that's an option, I'd run it by other maintainers first before spending time on it for discussion and to confirm they're happy with it.

Thanks for at least acknowledging this is an option. My guess is it won't be very well received, but I appreciate at least having it discussed for a moment.

Yes, I understand how it's no fun fiddling around with this code. What have you started changing in that branch? Again I'd suggest trying to make small incremental changes - I'd hope it'd be possible to replace LinearCombination with your approach without changing too much

I'll have a stab at just replacing LinearCombination on its own. This will be just a baby step towards solving the hash problem, but maybe we can get there incrementally.

[jagerber48]

@andrewgsavage ok, I tried to minimally replace linear_combination. I got somewhere, but this is already going to require API/test changes that I need to discuss with folks before moving on.

Currently one can do

x = ufloat(1, 1)
y = ufloat(2, 2)
z = 2 * x + 4 * y

ec = z.error_components()

assert x in ec
assert ec[x] == 2

That is the Variable x is a key to the error_components dictionary for z. The intuition here is that z is a variable which depends on x. In this sense Variables depend on other Variables. Some Variables are "atomic" in that their error components dictionary only has one component. But other Variables are composite in that their error components dictionary may have multiple components.

In the new framework the concept of Variable is essentially removed. Now we just have UFloat where a UFloat has a UCombo which consists of multiple UAtom. If variable z depends on variable x then this means that z and x share some UAtom. It would be that they share all UAtom, but it is possible that the weight of one or more of x's UAtoms in z is zero.

What we might have in this case (though not the general case of more complex variable dependencies) is

x = ufloat(1, 1)
y = ufloat(2, 2)
z = 2 * x + 4 * y

z_ec = z.error_components()
x_ec = z.error_components()

assert set(x_ec.keys()).issubset(z_ec.keys())
for x_uatom, x_weight in x_ec.items():
    assert z_ec[x_uatom] == 2 * x_weight

So in my proposal of using UAtoms independent from Variable, we lose the semantic idea that Variables can depend directly on Variables. Is this something we're ok with? It is still possible e.g. to calculate correlations between two UFloat, and, in the end, that's all we should really care about. That is, we don't care about the functional depedence between z and x, that is something that can be warped in strange ways as calculations go on. But we do care about the correlation between them, and that is something that is still retained in the new model.

[lebigot]

I think a better approach would be to work with the heavily/fully re-written branch but to get all of the old tests (and more) to work with minimal modification and then move forward guided by whatever diff ends up being required on the tests.

There is a feature of the original uncertainties that may not appear in existing tests: the speed of calculations. In particular, quite a few people need to sum many numbers with uncertainties, and the existing code is optimized for this. Non-regression tests that check the speed of summations are important for such users.

[jagerber48]

Ah, thank you very much for linking the release with information about the quadratic -> linear performance optimization. This is extremely helpful. My understanding is that the performance optimization is to essentially ensure that linear combinations of uncertainties are evaluated lazily. The new UCombo/UAtom has included this optimization. But with the detail found in your link I should be able to put together an actual performance test seeing if performance is impacted.

Note that either the total refactor or the incremental change strategy dramatically modifies the specific code responsible for this optimization, it is the main target of this change. So in either case I think it's important we have some performance testing in the test suite. Before we had been wondering if this performance optimization was necessary or if it was somehow a case of premature optimization. Your link makes it very clear that it is a very important optimization.

I don't know how to properly do performance testing in the test suite (I have guesses but they feel naive) but I can go ahead and learn how to do that. Or someone can share if they already know how.

[jagerber48]

str(sum(UFloat(1, 1) for _ in range(1000)))

This is a great test. My UCombo code fails on this test with a max recursion depth error during the combination expansion code. I'm trying to do memoization which should help with that but clearly I'm missing something.

[jagerber48]

This code "just works" for pickling and copying. This isn't so surprising since the objects are all immutable and there's not really anything much more complicated then a nested tuple.

[andrewgsavage]

pint uses codspeed https://github.com/hgrecco/pint/blob/master/.github/workflows/bench.yml I haven't looked into it myself

…

On Fri, Aug 9, 2024 at 2:36 PM Justin Gerber ***@***.***> wrote: Ah, thank you very much for linking the release with information about the quadratic -> linear performance optimization. This is extremely helpful. My understanding is that the performance optimization is to essentially ensure that linear combinations of uncertainties are evaluated *lazily*. The new UCombo/UAtom has included this optimization. But with the detail found in your link I should be able to put together an actual performance test seeing if performance is impacted. Note that either the total refactor or the incremental change strategy dramatically modifies the specific code responsible for this optimization, it is the main target of this change. So in either case I think it's important we have some performance testing in the test suite. Before we had been wondering if this performance optimization was necessary or if it was somehow a case of premature optimization. Your link makes it very clear that it is a very important optimization. I don't know how to properly do performance testing in the test suite (I have guesses but they feel naive) but I can go ahead and learn how to do that. Or someone can share if they already know how. — Reply to this email directly, view it on GitHub <#251 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADEMLEHQZEXKTTMFIJAARO3ZQTAV5AVCNFSM6AAAAABK2LBTYKVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTAMRYG4ZDOMY> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[jagerber48]

Ok, I feel like I'm finally up-to-speed on the whole "situation" we are having with respect to the data models for LinearCombination, AffineScalarFunc and Variable. I think this whole topic was spurred on my two things:

• First is the discussion about making the uncertainties objects hashable. e.g. Implement hash for AffineScalarFunc #219 or Add hash function for AffineScalarFunc class #189.
• Second is the general "confusingness" of the interplay between these three classes.

One of the centerpieces of the confusion is the fact that a Variable is and AffineScalarFunc but an AffineScalarFunc has a LinearCombination which is a mapping from Union[Variable, LinearCombination] to float (in code, prior to "expansion", it's represented as a list of tuples, but it is still a type of mapping even if it's not technically a dict). However, when we create a Variable the LinearCombination it has maps self to 1.0. This seemingly circular dependency is cause for concern whether or not the logic is sound.

A note on LinearCombination. LinearCombination has its main linear_combo attribute. This attribute is either a dict mapping dict[Variable, float] or a list of tuples mapping List[Tuple[float, Union[Variable, LinearCombination]]] This former dict is the expanded form while the latter list of tuples is the un-expanded form. The unexpanded form is advantageous because it makes operations between AffineScalarFunc O(1) in time because we simply combine the LinearCombinations of one or more AffineScalarFuncs that are being operated together without expansion into their base "Variable's".

Note that Variable has a nominal_value and std_dev.

---

To understand why things are the way they are it is important to understand AffineScalarFunc.derivatives and AffineScalarFunc.error_components. Every AffineScalarFunc has a LinearCombination. The expansion of this LinearCombination into its expanded dict[Variable, float] is the AffineScalarFunc.derivatives property. If y.derivatives = {x: 2.0} then this encodes that dy/dx = 2. So we can calculate the std_dev of an AffineScalarFunc using its derivatives according to

df = sqrt(((df/dx)*dx)**2 + ((df/dy)*dy)**2 + ((df/dz)*dz)**2)

So, for each Variable x, y, z there are two factors that determine the contribution to the std_dev of f: The derivative with respect to the Variable and the std_dev of the variable. The derivatives are initially stored in the LinearCombination linear_part attribute on f and then they are stored in the derivatives attribute. The std_devsare stored in thestd_devof the composingVariable`s.

It is central to the current architecture of AffineScalarFunc, Variable and LinearCombination that it is possible to extract the derivative of AffineScalarFunc with respect to Variables that it depends on. This is also where we pick up the circular dependency. When we create a new Variable x it doesn't depend on any other variables than itself. So since x is also an AffineScalarFunc, it must contain a LinearCombination, and this LinearCombination must map {self: 1.0}.

---

Some downsides of this approach.

The "objects" that we work with are always AffineScalarFunc, but they are sometimes also Variable.

from uncertainties import ufloat

x = ufloat(1, 1)
print(type(x))
# <class 'uncertainties.core.Variable'>
print(type(1.0*x))
# <class 'uncertainties.core.AffineScalarFunc'>

As soon as any operation is done on a Variable it is converted to an AffineScalarFunc.
This dual nature feels a bit superfluous at first glance to the developer and would be pretty confusing to a user introspecting into these classes. I especially imagine a user type checking against these types. For example they might see that type(ufloat(1, 1)) is Variable, so they would type check if isinstance(x, Variable): ... but be disappointed when x**2 fails their check.

This dual nature also has a downside involving derivatives. Recall that any AffineScalarFunc has derivatives.

from uncertainties import ufloat

x = ufloat(1, 1)
y = ufloat(2, 2)
z = ufloat(3, 3)

combo = 3*x + 2*y + z
print(combo.derivatives)
# defaultdict(None, {3.0+/-3.0: 1.0, 2.0+/-2.0: 2.0, 1.0+/-1.0: 3.0})

This makes sense. However, we get the sad

print(x in combo.derivatives)
True
print(1*x in combo.derivatives)
False

The problem is that x is the variable on which combo depends, but 1*x, even though it is equal to x, is an AffineScalarFunc and is not the Variable on which combo depends.

---

My feeling on this is that the derivatives architecture is nice but that the pitfalls are problematic enough that I'd be happy to jettison the derivatives approach for another proposal like my own. However, this will involve an API change to the user so I'm interested to get feedback from the others.

My proposal is the following:

• There are UAtom objects which are uniquely identifiable (they hold a uuid) and is to be thought of as having a standard deviation of 1.0 and is statistically independent from all other UAtoms.
• AffineScalarFunc becomes the main object of uncertainties. It is renamed to UFloat with AffineScalarFunc remaining as a legacy name for some amount of time > 0.
• UFloat has a float nominal_value
• UFloat has an uncertainty represented by a UCombo object (very similar to the LinearCombination object) which maps Union[UAtom, UCombo] to float.

With respect to the old way of doing things what has happened is that the Variable object has been removed. Uncertainties are now fundamentally held in the UAtom objects. But UAtom objects aren't compatible with arithmetic the way Variable objects were. Also, now users directly instantiate UFloat/AffineScalarFunc objects when working with uncertainties. Either through the UFloat constructor or the familiar ufloat helper function (could be converted to an alias for the UFloat constructor FYI). Users will likely never directly instantiate UAtom objects (we have the option to either expose them to the public API or not).

What happens in this new framework is that when operations are performed with two UFloat objects we find that the resulting UFloat object's UCombo will contain UAtoms from both the constituents. That is we will have something like

from uncertainties import ufloat

x = ufloat(1, 1)
y = ufloat(2, 2)
z = ufloat(3, 3)

combo = 3*x + 2*y + z

print(set(x.uncertainty.uatoms).issubset(combo.uncertainty.uatoms))
# True
print(set(y.uncertainty.uatoms).issubset(combo.uncertainty.uatoms))
# True
print(set(z.uncertainty.uatoms).issubset(combo.uncertainty.uatoms))
# True

But the advantage here is that we also will have

print(set((1*x).uncertainty.uatoms).issubset(combo.uncertainty.uatoms))

In this new framework correlation can be detected by looking at the intersection of the sets of UAtoms on which two UFloats, x and y depend. In fact this is exactly how the covariance of two UFloat can be calculated.

The big thing that is lost here is that we don't explicitly keep track of the dependency between an AffineScalarFunc and the Variables which went into generating that AffineScalarFunc. That is, if we do

x = ufloat(1, 1)
y = 2*x

there is no explicit record that dy/dx = 2. For simple examples like this, we can infer something like that by looking at the relative weightings of the single shared UAtom but for more complex examples any relationship is lost.
However, I argue that this isn't really a major loss. In the current framework we have

x1 = ufloat(1, 1)
y1 = x1 / 2
z1 = 3 * y1
print(z1.derivatives[x1])
1.5

y2 = ufloat(0.5, 0.5)
x2 = 2 * y2
z2 = 3 * y2
print(z2.derivatives[y2])
3.0

It is basically incidental which exact expression a given AffineScalarFunc depends on. Programmatically any AffineScalarFunc depends on the Variables that went into creating it but mathematically there isn't really anything privileged for Variables compared to AffineScalarFunc.

The new framework eliminates Variable and therefore the unearned privilege of some expressions.

---

In the new framework instead of relationships between variables being captured by the derivatives attribute on AffineScalarFunc, the statistical relationship between variables can be quantified by the covariance between ANY two UFloats.

x = ufloat(1, 1)
y = 1*x
z = 3 * y

assert x.covariance(x) == (x.std_dev)**2 
assert x.covariance(y) == (x.std_dev) * (y.std_dev)
assert z.covariance(x) == (x.std_dev) * (z.std_dev)
assert z.covariance(y) == (y.std_dev) * (z.std_dev)

These relationships will of course be more complicated in the multi-variate case where the answer depends on the specific statistical covariances.

x = ufloat(1, 1)
y = ufloat(2, 2)
z = 10*x + 3*y

assert z.covariance(x) = (x.std_dev) * (10 * x.std_dev)
assert z.covariance(y) = (y.std_dev) * (3 * y.std_dev)

---

A note: above my proposal is to basically remove the derivatives attribute of AffineScalarFunc. If we have that UAtoms represent random variables with standard deviation of unity this has to be done. If we let UAtoms have arbitrary non-negative std_dev then we can retain the difference between UFloat.derivatives and UFloat.error_components, but, from a computational perspective, this is redundant and, in my view, I don't see what is gained. For example, it won't be immediately obvious that for z=3*x that dz/dx = 3 like it is in the current implementation. Rather, all that will be seen is that z's dependence on the UAtoms that x depends on will be three times that of x's dependence.

---

For now, I imagine users will construct UFloat objects py passing nominal_value and std_dev floats. In this case a UCombo with a corresponding UAtom will be created with a weight corresponding to the user-passed std_dev. Recall that UAtom always has a std_dev of 1.0 so the standard deviation of UFloats are always stored in the weights in the UCombo for the different UAtoms. But, importantly, note that there is no "special" relationship between the UAtom that is created during the UFloat construction and the UFloat that is constructed other than the fact the new UAtom is the unique UAtom on which the new UFloat depends.

We could open the possibility for users to create their own UAtom and UCombo objects and use these to construct UFloat but I'm not sure what is really gained by this

---

Note that UFloat.error_components() will still persist. But instead of a dict mapping dict[Variable, float] it will map dict[UAtom, float].

---

So question to the group, @lebigot, @newville , @andrewgsavage , @wshanks what are peoples thoughts on removing (1) the Variable object and (2) the AffineScalarFunc.derivatives property in light of the discussion above? Right now I have a branch that has gotten through some of the tests removing Variable and AffineScalarFunc.derivatives. But at this point, seeing the situation clearly, I think it's the right time to discuss the proposed changes to the interface.

Names and implementation details are of course up for discussion. The only comment I'll make about the implementation here is that the important "lazy linear combination expansion" performance optimization lebigot put in years ago is just as possible in the new framework as it is in the old.

[newville]

@jagerber48 +1 on dropping the name AffineScalarFunc. I might slightly prefer UVariable or Variable over UFloat. That is especially because I sort of think that "Float" (ie, scalar) is a limitation we may not need to impose. That is, are we sure that the "Value with Uncertainty" class can not accommodate Complex numbers or NDArrays?

But also: I get that ufloat() is the main interface, and we really cannot break that. And ufloat() returning an object of type UFloat makes sense.

I'm OK with keeping "lazy derivatives" but I am also skeptical about potential performance problems for scalar calculations. I would suggest focusing on an implementation we all can understand and maintain.

[lebigot]

Thank, you for the description of your method: I now see what you're doing, with UAtom and your "new proposal" with the UAtom → float mapping. I must say that this idea popped up a few times while I was working on uncertainties, so I feel that it may be a fine, alternate approach to the current approach.

I also found the self-referential nature of Variables a bit surprising but it only says that dx/dx = 1; this comes from the fact that a variable x can also be considered as a function of itself (id(x) = x). This keeps the code handle both functions (or "dependent variables") and variables streamlined (and fast, since there is no test on whether an object with uncertainty is a function or a variable). The UAtom approach should also give a streamlined code, so this is not a criterion for deciding between both approaches.

The main questions are apparently around the importance for users of (1) the variable / function distinction, (2) more importantly, an access to derivatives and (3) even more importantly, an access to sources of uncertainty through error_components. Let me go through these points in order.

Variable/function distinction

Independent variables and functions of these variables are fundamentally different objects, so I'm not sure why users should lose this distinction that they currently have (that they can either use or ignore). Thus, in the example y = 1*x, it is still (mathematically) true that y is a function of x (which happens to be the identity), so I don't find this to be a problem. The most important point, here, is that removing the concept of independent variable removes a feature. I guess it's not used much, but it might be used by some.

Derivatives

More people might use derivatives than check the class of a number with uncertainty (Variable or not), but again calculating is not the main goal of this package and I've never seen anybody use that.

Your reason for removing derivatives is centered on the idea that covariance is all that matters. I don't see the problem with z in y.derivatives returning false: z is indeed not a variable of function y. The distinction between independent variables and functions of these variables is mathematically fundamental.

That said, it's a valid question to consider whether all that matters to users may be covariances.

Sources of uncertainty

My answer to this question is that it's mostly true that covariances between quantities is all that matter, but users would lose a useful feature, if the distinction between variable and function is removed.

In fact, more than one scientist relies on error_components: error budgets are a typical feature in metrology papers, and scientific paper in general. With the current approach, they can take any value with uncertainty, and get an answer to the question: "where does my total uncertainty of my final result come from?". If I understand correctly, with the new proposal, they can do this, but at a high price: they need to remember all the variables that went into their final result, which is both cumbersome and error prone:

Conclusion

In conclusion:

1. I find handling uncertainties internally though atomic sources of uncertainty with zero nominal value and unit standard deviation reasonable in principle. This still allows covariances to be calculated.
2. However, so far, this eliminates the concept of independent variable, which removes user features: the ability to know which variables are used in a code (probably not used much), and more critically, the ability to get a direct and robust access to the sources of uncertainty of a final result through error_components (which is very useful for many scientists). If I try to think of a way of keeping this important user feature, I get back to the current choice of implementation, essentially. But maybe you could think of something else and put back this feature in your new proposal?

In order to know for sure what users need, maybe it would be reasonable to ask those who rely on the features that may disappear to raise their hand? Maybe with a Python warning that asks for that and says that the variable / non-variable distinction and, more importantly, the error_components and derivatives might be removed soon (and explains how to disable the warning, for the others)?

[newville]

@jagerber48 @lebigot (sorry, might be a bit late and/or out-of-order).

I agree that @jagerber48 covariance is more useful than derivative.
I like the UAtom approach. I would also call it elegant, and likely to work. So, I'm +1 on expecting this to be the solution, or at the very least the core of the solution. I'm OK with going slow and incrementally ;).

Removing the self-referential "Linear Combination" would almost certainly make all the challenges with "id" and pickling either go away or at least be much simpler to deal with. +1.

It would be helpful to be able to query a result of some calculation

x = ufloat(1, 0.2)
theta = np.sin(x / 2.0)
scale = ufloat(100, 0.8)
m = np.exp( - theta**2/scale)

and be able to discern that the uncertainty in m is related to that in x. I think m.covariance(x) != 0 is fine. But maybe also one would like a list of all this UVariables that m is related too (x, scale, .... maybe theta?).

I think I said this at some point, but it might be nice to be able to define a "context" -- effectively a namespace where a bunch of potentially related UAtoms and UVariables live. Currently, it is a little hard to reason about whether that is __local__ or __globals, and maybe it should be something settable by the user.
I think that this is neither a high priority, nor all that difficult to imagine with your approach.

[jagerber48]

@lebigot Thank you very much for that response. Yes, you're understanding the proposal now.

The main questions are apparently around the importance for users of (1) the variable / function distinction, (2) more importantly, an access to derivatives and (3) even more importantly, an access to sources of uncertainty through error_components. Let me go through these points in order.

Yes, these are the main questions.

Independent variables and functions of these variables are fundamentally different objects, so I'm not sure why users should lose this distinction that they currently have (that they can either use or ignore). Thus, in the example y = 1*x, it is still (mathematically) true that y is a function of x (which happens to be the identity), so I don't find this to be a problem. The most important point, here, is that removing the concept of independent variable removes a feature. I guess it's not used much, but it might be used by some.

The new proposal does drop the distinction between independent variables and functions of variables. Instead, it puts all of the user-facing uncertainties objects on equal footing. Another Variable/AffineScalarFunc example is the awkward (to me) fact that if y=1*x and z=2*y then z has a special relationship with x, but not y, i.e .x in z.derivatives but y not in z.derivatives despite z=2*y.

The distinction between independent variables and functions of these variables is mathematically fundamental.

I'm not so sure about this. The y=1*x and z=2*y example illustrates this. Mathematically when we have multi-variable functions, there is ambiguity about implicit and explicit total derivatives with respect to different variables. I would say it is just as reasonable for a user to be interested in dz/dy as it is for them to be interested in dz/dx, but the current uncertainties architecture privileges one of these quantities over the other. It's a philosophical dislike for this asymmetry that motivated me to eliminate the distinction and instead move towards "everything being equally treated random variables" architecture.

Sources of uncertainty

...

To be clear, in the new framework the UVariable class will still have an error_contributions attribute. But instead of returning a dict that maps IndependentVariable -> float it will return a dict that maps UAtom -> float. How the user can use this or some other feature to answer

"where does my total uncertainty of my final result come from?"

is a very good question that I want to explore further. One confession I have is that it's not really clear to me how the tag. There's at least one example in the docs:

syst_error = math.sqrt(
    sum(  # Error from *all* systematic errors
        error**2
        for (var, error) in result.error_components().items()
        if var.tag == "systematic"
    )
)

So one idea I've already toyed around with in code is that we could continue to have a tag attribute, but it would stay on the UAtom, not the UVariable. So we would have

x = UVariable(1, 0.1, tag="systematic")

with pytest.raises(AttributeError):
    x.tag

x_uatom = next(iter(x.error_components))
assert x_uatom.tag == "systematic"

and we could realize the example above with

syst_error = math.sqrt(
    sum(  # Error from *all* systematic errors
        weight**2
        for (uatom, weight) in result.error_components().items()
        if uatom.tag == "systematic"
    )
)

I think with some care this would allow users to acheive some error budget tracking functionality similar to what exists now. We could even supply helper filtering functions that allow users to filter the error_components UAtom -> float dictionary by tag or extract weights by tag. Maybe with the tag defaulting to the uuid in case no tag is supplied by the user.

@lebigot I directly agree with your statement that this new proposal drops a user feature. That's exactly why I'm trying to carefully have this conversation. I think the path forward on this question is to understand how this feature is used (i.e. flagging concrete use cases as code snippets) and figuring out what functionality can/should retained in the new proposal.

In order to know for sure what users need, maybe it would be reasonable to ask those who rely on the features that may disappear to raise their hand? Maybe with a Python warning that asks for that and says that the variable / non-variable distinction and, more importantly, the error_components and derivatives might be removed soon (and explains how to disable the warning, for the others)?

I'm open to this as an attempt at a data gathering tool. Just a reminder that error_components won't be removed but it will be changed. A proposed message for any users of either parameter could be

msg = 'uncertainties proposal: There is an active proposal under consideration to remove the `derivatives` property and to change the behavior of the `error_components` function. If these properties are important to you we invite you to share your use case and thoughts on the proposal by posting at https://github.com/lmfit/uncertainties/discussions/... . To suppress this warning please use `warnings.filterwarnings("ignore", msg="uncertainties proposal: ")`'
warnings.warn(msg, UserWarning)

I don't think DeprecationWarning would be appropriate until we actually decide that it is appropriate to move forward with the propsoal.

---

@newville

Thanks for sharing the points you agree with.

About this:

I think I said this at some point, but it might be nice to be able to define a "context" -- effectively a namespace where a bunch of potentially related UAtoms and UVariables live.

Yes, this is an interesting idea that would help with some of the issues. Perhaps an approach like https://stackoverflow.com/questions/328851/printing-all-instances-of-a-class or https://stackoverflow.com/questions/54000173/how-to-get-all-instances-of-a-class could be followed to register and track each UVariable that is created. Then there could be a helper function that takes one UVariable and returns all the other UVariable with which it is correlated.

The global nature worries me, and I also worry that this cache/registry could get quite large if users are dealing with lots of UVariable, e.g. the print(sum(ufloat(1, 0.1) for _ in range(1000000))) test. I also worry if maintaining this registry could slow down performance.

An idea that would ensure no performance hit to uncertainties would be a strategy where no active registry is maintained, but, at user-request, the __locals__ and/or __globals__ can be filtered for any matching UVariable objects. Again, I'm not sure how easy such an approach would be to implement.

[lebigot]

I see two main subjects in your last posts, @newville and @jagerber48: error budgets, and derivatives that are not with respect to independent variables, that I'll address in turn.

Error budget and UAtoms

@jagerber48: thank you for clarifying that the UVariable class still has an error_components (or similar) attribute. This addresses the error budget tracking that some users need.

I introduced the tag system only so that users easily identify the variables/sources of errors that their quantities depend on (otherwise, sources of uncertainties would be displayed as an "anonymous" … ± … quantity): (independent) variable can be named (thanks to the tags).

Your example with a UVariable works well: people can identify the contribution of each tagged source of uncertainty.

I agree that the right place for the tag is a UAtom (as they are the ultimate sources of uncertainty).

As for the API, I'm thinking that most current users of error_components would not need any access to its UAtoms (as they are all "the same"): error_components could simply return (tag, weight) pairs, or maybe, for a self-explanatory API, (tag, variance) pairs. I'm not sure if it's super useful to return the id() of the UAtom; doing so would also in principle constrain users to not use integers for their tags (which they might want to do, if they automatize the creation of many variables and want to tag them with integers).

In conclusion:

1. I now consider that the "unit standard errors" (UAtom) approach would allow most users to work mostly like before.
2. derivatives (i.e. derivatives with respect to variables) would not be available anymore, but I don't expect too much suffering from that. Warning users ahead of time and suggesting that they give feedback seems like a considerate way to proceed.
3. error_components would change its semantics, since Variables would be gone. I believe that its main application is error budgets, and @jagerber48 you showed that this feature can remain with UAtoms. Again, a warning that mentions that it will be phased out and replaced by another, similar method (I suggest providing (tag, variance) pairs, with a the current None default tag if none was given).
4. In all honesty, I'm not sure to see the benefits for users of an implementation with unit random variables, since they would be losing features—even if not too many. My argument that variables and functions are fundamentally different, for users, is very clear in the fact that in Python, expressions are isomorphic to trees, and like everybody knows, trees have something specific called leaves—which represent random variables, here—; so I see as perfectly acceptable to distinguish between both with two different classes (Variable and AffineScalarFunc), especially since they are usually hidden from users.

Derivatives with respect to "functions"

There are quite many points I could share on the subject of allowing (or not) users to calculate things like dz/dy in your example, @jagerber48, or dm/d_theta in @newville's example), but I can share them later. For the time being, I'll just say that this opens relatively complicated questions, in more general cases than y = 1*x, z = 2*y so that I could want dz/dy (that mostly have answers, I believe). Also I'm only seeing limited benefits (essentially, I'm not seeing use cases that cannot already be addressed). I'd be happy to keep discussing the subject, though!

I'll note one thing, because I'm afraid that there is some confusion (that could prove problematic): z.covariance(y) != 0 is very different from "y is a variable of z", so this test cannot be used for listing the "variables" that z depends on. A case like z = u + v and y = t + v is one examples among many.

Ah, and I share @jagerber48's unease with playing with some global list of variables (for the same reasons).

In any case, it seems to me that more general derivatives would be an additional feature, which is less of a priority than exploring the unit random variable implementation.

[newville]

@jagerber48 @lebigot um, the current uncertainties linear combinations mechanism uses globals.

>>> x = uncertainties.ufloat(10, 1)
>>> y = uncertainties.ufloat(20, 0.4)
>>> z = x + y
>>> z
30.0+/-1.077032961426901
>>> z._linear_part.linear_combo
defaultdict(<class 'float'>, {20.0+/-0.4: 1.0, 10.0+/-1.0: 1.0})

z has to reference x.

OTOH, using a context would allow a set of UVariables to share a correlation matrix, but have no risk of confusion between different contexts. The context would hold a set of variable (not much more than {name:nominal_value} pairs), and a correlation matrix. Any calculation using the variables in the same context would propagate uncertainties with the contexts correlation. The correlation could changed within a context. Variables in different contexts are uncorrelated.

Sure, the default context could be a "global state". Or a user could create a new one, set it as the current one, and return to another.

[jagerber48]

Ok, I've been busy for a long time, but now I have some limited time to look at this again.

I think we actually left the conversation off in a pretty good spot a few months ago. The conclusion from my perspective was that:

• We are agreed to move from AffineScalarFunc + Variable to UFloat + UAtom (in concept, not necessarily name).
• Notably this change will eliminate the concept of derivative from uncertainties. Instead of covariance being tracked by noting dy/dz = c it will be tracked by noting that both y and z depend on a shared UAtom. This information will be tracked in the error_components attribute. I'll comment that this is an appreciable change to how the package work, and it may even change usage for some users. But the result of the previous conversation is that we are in agreement that this is all ok.
• We should attempt to make this upgrade as incrementally as possible. That is, we shouldn't e.g. use the wholesale rewrite that I worked out as a concept. I think this will be possible but there will be one pretty large PR that makes most of the change. I think it should be possible to contain this change mostly to the API without modifying numerical internals. Downstream changes could address some of the hard-to-follow monkey patching going on in the internals.

I'll try to move forward on this as I have time.

[newville]

@jagerber48 thanks (and sorry for the delay in responding). #262 looks like it is coming along nicely.

Do you have a suggestion for the steps to do this incrementally?

Perhaps tagging this as 4.0.0.rc1 could be an okay approach, accompanied by an explicit statement of not pushing to PyPI until it reaches 4.0.1 or something. That might signal "please test and raise concerns" to the other developers and brave users.

[jagerber48]

I'm finding it challenging to make this change more incrementally than I am now. What I am leaving off compared to my wholesale rewrite is that I am keeping all of the underlying "uncertainty propagation" code. Basically I'm leaving ops.py relatively untouched.

I think all of these changes need to go together:

1. Change from AffineScalarFunc + Variable to UFloat + UCombo + UAtom
2. Remove the derivative property. Now error_components maps UAtom -> float.
3. UFloat objects don't have tags now, UAtom objects have tags now.
4. Copying and pickling/unpickling have different, but more sensible, behaviors now

(1) is obviously the goal for this whole push.

(2) requires the most substantial changes to the tests. In some cases, like if a test is looking at the relationship between UFloat x and UFloat y which was directly derived from x then the test can often be rewritten using get_single_uatom(x) and get_single_uatom(y) and applying the derivative test to the single UAtoms on which x and y depend. But in some cases tests need to be entirely removed. In some cases where entire tests are removed it is sometimes possibly to replace the tests with similar tests in the new framework, but there are a few tests in the old framework which I don't think have an analog in the new framework and are just obsolete.

(3) is minor and I think we're all ok with it

( 4) I need to open an issue discussing this. I find the old copy behavior to be very confusing. I don't think there will be objections to the change here but the discussion should be captured.

---

The status going forward is I need to continue to make changes to tests. I also haven't touch the "array" side of things at all. I'm not sure if that will be easy or hard. In any case, @newville and others, looking at #262 as it is right now, do these changes look appropriate for e.g. 4.0.0.rc1? Or are they too sweeping and we need to find a way to break them down more?

[andrewgsavage]

you could make a branch of lmfit that uses this refactor before a 4.0.0.rc1, that would help show it in action.

[newville]

@jagerber48 Thanks for that. I agree with all your steps and design considerations.

Doing steps 1, 2, and 3 in the current 3.2 (or 3.3) series seems fine to me. For step 4, "yes please" on simplifying pickle/unpickle.

I would also be OK with making a "4.0rc" branch (or just calling your branch for #262 that) and expecting that this is a larger, not-so-incremental rewrite. But if you think steps 1, 2, and 3 can be (and are better) done incrementally, I'm OK with that too.