Added polynomial class and its tests

[jsmitty-here]

• (Not related to any known open issue)
• Executed pre-commit run --all-files with no errors
• The change is fully covered by automated unit tests
• Documented in docs/ as appropriate
• Added an entry to the CHANGES file

Explanation

I've been using pint for a project I am working on for fluid calculations. I needed a polynomial to model certain equations where the input and outputs of the polynomial have their own units; such as a pump curve which takes a flowrate unit and outputs a pressure unit, or vice versa.

I built a Polynomial class inheriting from the numpy Polynomial class that incorporates x & y unit variables for my own project and thought this could be a helpful class for the pint library as well. I have implemented unittests for this class in the test_polynomial.py file.

This is my first commit to this repo, so please let me know if I need to make any changes before this PR can be merged.

Example Uses

1. Example modeling a pump curve

flow_unit, pressure_unit = ureg.gallon / ureg.minute,  ureg.psi

pump_curve = Polynomial([46.029, -0.1546, -0.0053], flow_unit,  pressure_unit)
46.029 - 0.1546 * x - 0.0053 * x ** 2

pump_curve.solve(0 * flow_unit)  # Returns: 46.029 pound_force_per_square_inch
pump_curve.y_intercept         # Returns: 46.029 pound_force_per_square_inch
pump_curve.x_intercept          # Returns: [-108.91112924   79.74131792] gallon / minute
pump_curve.solve_for_x(10 *pressure_unit , min_value=0))      # Returns: 69.144 pound_force_per_square_inch
pump_curve.derivative_at(30 * flow_unit)   # Returns: -0.4726 minute * pound_force_per_square_inch / gallon

2. Example modeling a distance vs time polynomial

distance_unit, time_unit = ureg.meter,  ureg.second

distance_curve = Polynomial([0, 1, 2], time_unit ,  distance_unit)
0.0 + 1.0 * x + 2.0 * x ** 2

distance_curve.solve(0 * time_unit)  # Returns: 0.0 meter
distance_curve.y_intercept     # Returns: 0.0 meter
distance_curve.x_intercept     # Returns: -0.5 second
distance_curve.solve_for_x(10 * distance_unit, min_value=0))       # Returns: 1.9999999999999996 second
distance_curve.derivative_at(2 * distance_unit)    # Returns: 9.0 meter / second

velocity_curve = distance_curve.derivative  # Returns: Polynomial that is the derivative of the original
1.0 + 4.0 * x

velocity_curve.solve_for_x(5 * distance_unit / time_unit ) # Returns: 1.0 second

acceleration_curve = velocity_curve.derivative   # Returns: Polynomial that is the derivative of the original
4.0

acceleration_curve.solve(2 * time_unit)        # Returns: 4.0 meter / second ** 2

[hgrecco]

While I do agree that it is desirable to improve the numpy/scipy support, I would like to discuss a few things:

1. Is there a way to do it in a more direct manner (e.g. using wraps to avoid rewriting the internal code).
2. How should we organize it? Because it seems that this might explode easily. Maybe it is better to put it in a numpy support module?

[jsmitty-here]

Thanks for the response. Regarding your discussion points:

1. Do you have a vision on how this would be done? I may be misunderstanding what you mean. By a 'wrap to avoid rewriting the internal code' I assume you are referring to the dunder methods? If so I will certainly need to investigate how I could change this. Maybe something like this:

from numpy.polynomial import polynomial as p

def numpy_polynomial(cls):
      
    class Wrapper:
         """ Handles Integration with numpy Polynomial base class """

    return Wrapper
  
@numpy_polynomial
class Polynomial(p.Polynomial):
    """ Handles Initialization and any new polynomial functions. """

    def __init__(self, coef: list[float], x_unit: Unit, y_unit: Unit):
        super().__init__(coef)
        self.x_unit= x_unit
        self.y_unit= y_unit

Please let me know if that is not what you had in mind.

2. Agreed, I could see mathematical implementations incorporating pint growing beyond polynomials. Are you suggesting the 'polynomial' module (and modules like it) should be organized into a numpy support directory within pint or be segregated from the main pint library?

This is my first attempt at contributing to an open source library so I apologize for any informalities. Additionally, If the polynomial feature does not align with the overall goal of the pint library just let me know and I will take no offense. Like I said, I had already wrote this class and it's tests for a personal project I am working on.

[hgrecco]

(This is a great contribution, first or otherwise. Thanks for submitting it)

Regarding (1), I was wondering if using this could be useful.

In relation to (2), indeed. I think a numpy support module might be useful. What do you think @jules-ch

[jsmitty-here]

I may need some assistance with using wraps of this type. My initial attempt was to add the wraps into the class methods:
Converting a function like this:

class Polynomial(p.Polynomial)
    ...
    def solve(self, value: Quantity, min_value: float = -inf) -> Quantity:
        return max(self(value), min_value) * self.y_unit

to this:

class Polynomial(p.Polynomial)
    ...
    @staticmethod
    @ureg.wraps(self.y_unit, (self.x_unit, self.x_unit))
    def solve(value: float, min_value: float = -inf) -> float:
        return max(self(value), min_value)

however, I can not access the x_unit, and y_unit instance variables from self for the ureg.wraps nor can I use self inside the functions as the method now needs to be a static method to be compatible with the ureg.wraps decorator function.

Additionally, I do not see how these ureg unit conversion wraps could be used on any of the functions or dunder methods which return 'Polynomial' as they return a new instance of the polynomial class, not a float or Quantity.

I do agree that there is a lot of repeating code considering the class primarily contains just a numpy polynomial, and it's unit values. I just have not discovered the best way to implement that yet.

[jsmitty-here]

Ok, after experimenting, I have discovered a way to make this work for the dunder methods:

def wrap_numpy_poly(dunder_method):

    @wraps(dunder_method)
    def wrap_numpy_poly_inner(pint_polynomial_class, other):
        """ Call the method of the numpy polynomial super class by method name and provide it the 'other' argument.
              Adapt the units if needed (multiplication, division, etc.)  """
        x_unit, y_unit = pint_polynomial_class.x_unit, pint_polynomial_class.y_unit
        if isinstance(other, pint_polynomial_class.__class__):
          try:
            x_unit = getattr(x_unit, dunder_method.__name__)(other.x_unit)
            y_unit = getattr(y_unit, dunder_method.__name__)(other.y_unit)
          except AttributeError:
            pass
          finally:
            new_np_poly = getattr(pint_polynomial_class._poly, dunder_method.__name__)(other._poly)
        else:
          new_np_poly = getattr(pint_polynomial_class._poly, dunder_method.__name__)(other)
        
        return pint_polynomial_class.__class__(new_np_poly.coef, x_unit, y_unit)

    return wrap_numpy_poly_inner

Then in the Polynomial Class, the dunder methods only require the wrapper. As long as the super class can handle the arithmetic, then the subclass will too.

class Polynomial(p.Polynomial):

    def __init__(self, coef: List[float], x_unit: Unit = ureg.dimensionless, y_unit: Unit = ureg.dimensionless):
        super().__init__(coef)
        self._poly = p.Polynomial(coef)
        self.x_unit = x_unit
        self.y_unit = y_unit

    @wrap_numpy_poly
    def __mul__(self, other) -> 'Polynomial':
        pass

    @wrap_numpy_poly
    def __add__(self, other) -> 'Polynomial':
        pass

    @wrap_numpy_poly
    def __sub__(self, other) -> 'Polynomial':
        pass

Now, I am curious if there is an even better way to add these wrappers to each necessary dunder method. Possibly a wrapper to the Polynomial class or by using a MetaClass. I will continue playing around with this.

[andrewgsavage]

Now, I am curious if there is an even better way to add these wrappers to each necessary dunder method. Possibly a wrapper to the Polynomial class or by using a MetaClass. I will continue playing around with this.

Take a look at how pandas does it. I'm not sure how you'd do it with wraps though.
https://github.com/pandas-dev/pandas/blob/01d0874a70a9d2957a8f96501e72a31b5da63717/pandas/core/arrays/base.py#L1705

[andrewgsavage]

I've written something similar for water pump calculations too. My approach was to use pint-pandas to convert units to metric (any coherent unit system would work) and then drop units so I can work with floats. This allows me to use code from examples and documentation much easier than would be possible with the unit aware Polynomial class, which needs each method to be added and tested. I can convert units at the end if needed, although I've not needed to do this as metric has been sufficient.

Is this approach (using pint or pint-pandas to convert values to a coherent system) something you considered? I think curve fitting is a common use case for pint, and I don't think the pint documentation points people to convert then drop units. I wonder if it's something worth adding - this is how every FEA and MBD simulation package I've used deals with units.

That being said, the Polynomial class is a nice idea and will be great to use when it's in a mature state. I think a really cool example would be to have two pumps curves with different units. Create a Polynomial object for each, and then can add the two Polynomial objects to get a pump curve for the two pumps in series.

[jsmitty-here]

@andrewgsavage
I have not yet investigated pint-pandas and maybe this could solve some of my issues.

In the application I am building, I am hoping to not drop units at any point. Hence the creation of a polynomial class with the units caked in.
My main driver for making this polynomial class was for pump curves and generating system curves (fitting a polynomial to the pressure drop through a piping system at a range of flowrates), however, I see how it could be used for any polynomial physics equations where units matter (see my first comment where the derivative of the 'distance_curve', with units 'second' vs 'meter', generated a new polynomial with units 'second' vs 'meter / second') .

Additionally, I envision the units of the pump curve not actually being important thanks to Pint. The user could input any Quantity (of comparable units) and convert the returned Quantity as needed. For example:

# Define units for simplicity
GPH = ureg.gallon / ureg.hour
PSI = ureg.psi
LPM = ureg.liter/ ureg.minute
PASCAL = ureg.pascal

curve_1 = Polynomial([20, 0, -0.05], GPH, PSI)      # f( x GPH) = y PSI = 20.0 - 0.05 * x ** 2
curve_2 = Polynomial([30, -0.1, -0.02], LPM, PASCAL) # f( x LPM) = y PASCAL= 30.0 + -0.1* x - 0.02 * x ** 2

# The class handles the unit conversion for the input argument internally
curve_1.solve(1 * GPH)    # 19.95 pound_force_per_square_inch
curve_1.solve(1 * GPH).to(PASCAL)    # 137550.41 pascal
curve_1.solve(1 * LPM)    # 7.44 pound_force_per_square_inch
curve_1.solve(1 * LPM).to(PASCAL)    # 51285.71 pascal

curve_2.solve(1 * GPH)    # 29.99 pascal
curve_2.solve(1 * GPH).to(PSI)    # 0.0044 pound_force_per_square_inch
curve_2.solve(1 * LPM)    # 29.88 pascal
curve_2.solve(1 * LPM).to(PSI)    # 0.0043 pound_force_per_square_inch

Now if the user actually wants to convert the entire polynomial to a new unit system, I have a method like such:

curve_1 = Polynomial([20, 0, -0.05], GPH, PSI)      # f( x GPH) = y PSI = 20.0 - 0.05 * x ** 2
# May rename this method from "to_units()" to "to()" to match the Quantity's "to()" unit conversion method.
curve_1.to_units(x_unit=LPM, y_unit=PASCAL)      # f(x LPM) = y PASCAL = 137895.14586336722 − 86609.4395918349 * x ** 2

In the case of adding two (or more) pump curves:

curve_1 = Polynomial([20, 0, -0.05], GPH, PSI)      # f( x GPH) = y PSI = 20.0 - 0.05 * x ** 2
curve_2 = Polynomial([30, -0.1, -0.02], LPM, PASCAL) # f( x LPM) = y PASCAL= 30.0 + -0.1* x - 0.02 * x ** 2

combined_curve = curve_1 + curve_2    # + ... + curve_N
print(combined_curve)    # f( x GPH) = y PSI = 20.004351132131905 - 9.150459358810584e-07 * x - 0.05000001154608556 * x ** 2