Add to_preferred

5 files changed, 269 insertions, 0 deletions

diff --git a/pint/compat.py b/pint/compat.py
index f5b03e3..e324d62 100644
--- a/pint/compat.py
+++ b/pint/compat.py
@@ -130,6 +130,13 @@ try:
 except ImportError:
     HAS_BABEL = False
 
+try:
+    import mip
+
+    HAS_MIP = True
+except ImportError:
+    HAS_MIP = False
+
 # Defines Logarithm and Exponential for Logarithmic Converter
 if HAS_NUMPY:
     from numpy import exp  # noqa: F401
@@ -141,6 +148,13 @@ else:
 if not HAS_BABEL:
     babel_parse = babel_units = missing_dependency("Babel")  # noqa: F811
 
+if not HAS_MIP:
+    class MissingDependency:
+        def __getattribute__(self, name: str):
+            missing_dependency("mip")()
+
+    mip = MissingDependency()
+
 # Define location of pint.Quantity in NEP-13 type cast hierarchy by defining upcast
 # types using guarded imports
 upcast_types = []
diff --git a/pint/facets/plain/quantity.py b/pint/facets/plain/quantity.py
index 5e4f33d..913504e 100644
--- a/pint/facets/plain/quantity.py
+++ b/pint/facets/plain/quantity.py
@@ -40,6 +40,7 @@ from ...compat import (
     eq,
     is_duck_array_type,
     is_upcast_type,
+    mip,
     np,
     zero_or_nan,
 )
@@ -689,6 +690,213 @@ class PlainQuantity(PrettyIPython, SharedRegistryObject, Generic[_MagnitudeType]
 
         return self.to(new_unit_container)
 
+    def to_preferred(self, preferred_units: List[UnitLike]) -> PlainQuantity[_MagnitudeType]:
+        """Return Quantity converted to a unit composed of the preferred units.
+
+        Examples
+        --------
+
+        >>> import pint
+        >>> ureg = pint.UnitRegistry()
+        >>> (1*ureg.acre).to_preferred([ureg.meters])
+        <Quantity(4046.87261, 'meter ** 2')>
+        >>> (1*(ureg.force_pound*ureg.m)).to_preferred([ureg.W])
+        <Quantity(4.44822162, 'second * watt')>
+        """
+
+        if not self.dimensionality:
+            return self
+
+        # The optimizer isn't perfect, and will sometimes miss obvious solutions.
+        # This sub-algorithm is less powerful, but always finds the very simple solutions.
+        def find_simple():
+            best_ratio = None
+            best_unit = None
+            self_dims = sorted(self.dimensionality)
+            self_exps = [self.dimensionality[d] for d in self_dims]
+            s_exps_head, *s_exps_tail = self_exps
+            n = len(s_exps_tail)
+            for preferred_unit in preferred_units:
+                dims = sorted(preferred_unit.dimensionality)
+                if dims == self_dims:
+                    p_exps_head, *p_exps_tail = [preferred_unit.dimensionality[d] for d in dims]
+                    if all(
+                        s_exps_tail[i] * p_exps_head == p_exps_tail[i] ** s_exps_head
+                        for i in range(n)
+                    ):
+                        ratio = p_exps_head / s_exps_head
+                        ratio = max(ratio, 1 / ratio)
+                        if best_ratio is None or ratio < best_ratio:
+                            best_ratio = ratio
+                            best_unit = preferred_unit ** (s_exps_head / p_exps_head)
+            return best_unit
+
+        simple = find_simple()
+        if simple is not None:
+            return self.to(simple)
+
+        # For each dimension (e.g. T(ime), L(ength), M(ass)), assign a default base unit from
+        # the collection of base units
+
+        unit_selections = {
+            base_unit.dimensionality: base_unit
+            for base_unit in map(self._REGISTRY.Unit, self._REGISTRY._base_units)
+        }
+
+        # Override the default unit of each dimension with the 1D-units used in this Quantity
+        unit_selections.update({
+            unit.dimensionality: unit
+            for unit in map(self._REGISTRY.Unit, self._units.keys())
+        })
+
+        # Determine the preferred unit for each dimensionality from the preferred_units
+        # (A prefered unit doesn't have to be only one dimensional, e.g. Watts)
+        preferred_dims = {
+            preferred_unit.dimensionality: preferred_unit
+            for preferred_unit in map(self._REGISTRY.Unit, preferred_units)
+        }
+
+        # Combine the defaults and preferred, favoring the preferred
+        unit_selections.update(preferred_dims)
+
+        # This algorithm has poor asymptotic time complexity, so first reduce the considered
+        # dimensions and units to only those that are useful to the problem
+
+        ## The dimensions (without powers) of this Quantity
+        dimension_set = set(self.dimensionality)
+
+        ## Getting zero exponents in dimensions not in dimension_set can be facilitated
+        ## by units that interact with that dimension and one or more dimension_set members.
+        ## For example MT^1 * LT^-1 lets you get MLT^0 when T is not in dimension_set.
+        ## For each candidate unit that interacts with a dimension_set member, add the
+        ## candidate unit's other dimensions to dimension_set, and repeat until no more
+        ## dimensions are selected.
+
+        discovery_done = False
+        while not discovery_done:
+            discovery_done = True
+            for d in unit_selections:
+                unit_dimensions = set(d)
+                intersection = unit_dimensions.intersection(dimension_set)
+                if 0 < len(intersection) < len(unit_dimensions):
+                # there are dimensions in this unit that are in dimension set
+                # and others that are not in dimension set
+                    dimension_set = dimension_set.union(unit_dimensions)
+                    discovery_done = False
+                    break
+
+        ## filter out dimensions and their unit selections that don't interact with any
+        ## dimension_set members
+        unit_selections = {
+            dimensionality: unit
+            for dimensionality, unit in unit_selections.items()
+            if set(dimensionality).intersection(dimension_set)
+        }
+
+        ## update preferred_units with the selected units that were originally preferred
+        preferred_units = list(set(u for d, u in unit_selections.items() if d in preferred_dims))
+        preferred_units.sort(key=lambda unit: str(unit))  # for determinism
+
+        ## and unpreferred_units are the selected units that weren't originally preferred
+        unpreferred_units = list(set(
+            u for d, u in unit_selections.items() if d not in preferred_dims))
+        unpreferred_units.sort(key=lambda unit: str(unit))  # for determinism
+
+        # for indexability
+        dimensions = list(dimension_set)
+        dimensions.sort()  # for determinism
+
+        # the powers for each elemet of dimensions (the list) for this Quantity
+        dimensionality = [self.dimensionality[dimension] for dimension in dimensions]
+
+        # Now that the input data is minimized, setup the optimization problem
+
+        ## use mip to select units from preferred units
+
+        model = mip.Model()
+
+        ## Make one variable for each candidate unit
+
+        vars = [
+            model.add_var(str(unit), lb=-mip.INF, ub=mip.INF, var_type=mip.INTEGER)
+            for unit in (preferred_units + unpreferred_units)
+        ]
+
+        ## where [u1 ... uN] are powers of N candidate units (vars)
+        ## and [d1(uI) ... dK(uI)] are the K dimensional exponents of candidate unit I
+        ## and [t1 ... tK] are the dimensional exponents of the quantity (self)
+        ## create the following constraints
+        ##
+        ##                ⎡ d1(u1) ⋯ dK(u1) ⎤
+        ## [ u1 ⋯ uN ] * ⎢    ⋮    ⋱         ⎢ = [ t1 ⋯ tK ]
+        ##                ⎣ d1(uN)    dK(uN) ⎦
+        ##
+        ## in English, the units we choose, and their exponents, when combined, must have the
+        ## target dimensionality
+
+        matrix = [
+            [preferred_unit.dimensionality[dimension] for dimension in dimensions]
+            for preferred_unit in (preferred_units + unpreferred_units)
+        ]
+
+        ## Do the matrix multiplication with mip.model.xsum for performance and create constraints
+        for i in range(len(dimensions)):
+            dot = mip.model.xsum([var * vector[i] for var, vector in zip(vars, matrix)])
+            # add constraint to the model
+            model += dot == dimensionality[i]
+
+        ## where [c1 ... cN] are costs, 1 when a preferred variable, and a large value when not
+        ## minimize sum(abs(u1) * c1 ... abs(uN) * cN)
+
+        ## linearize the optimization variable via a proxy
+        objective = model.add_var("objective", lb=0, ub=mip.INF, var_type=mip.INTEGER)
+
+        ## Constrain the objective to be equal to the sums of the absolute values of the preferred
+        ## unit powers. Do this by making a separate constraint for each permutation of signedness.
+        ## Also apply the cost coefficient, which causes the output to prefer the preferred units
+
+        ### prefer units that interact with fewer dimensions
+        cost = [len(p.dimensionality) for p in preferred_units]
+
+        ### set the cost for non preferred units to a higher number
+        bias = max(map(abs, dimensionality)) * max((1, *cost)) * 10  # arbitrary, just needs to be larger
+        cost.extend([bias] * len(unpreferred_units))
+
+        for i in range(1 << len(vars)):
+            sum = mip.xsum(
+                [(-1 if i & 1 << (len(vars) - j - 1) else 1) * cost[j] * var for j, var in enumerate(vars)]
+            )
+            model += objective >= sum
+
+        model.objective = objective
+
+        # run the mips minimizer and extract the result if successful
+        if model.optimize() == mip.OptimizationStatus.OPTIMAL:
+            optimal_units = []
+            min_objective = float("inf")
+            for i in range(model.num_solutions):
+                if model.objective_values[i] < min_objective:
+                    min_objective = model.objective_values[i]
+                    optimal_units.clear()
+                elif model.objective_values[i] > min_objective:
+                    continue
+
+                temp_unit = self._REGISTRY.Unit("")
+                for var in vars:
+                    if var.xi(i):
+                        temp_unit *= self._REGISTRY.Unit(var.name) ** var.xi(i)
+                optimal_units.append(temp_unit)
+
+            sorting_keys = {tuple(sorted(unit._units)): unit for unit in optimal_units}
+            min_key = sorted(sorting_keys)[0]
+            result_unit = sorting_keys[min_key]
+
+            return self.to(result_unit)
+        else:
+            # for whatever reason, a solution wasn't found
+            # return the original quantity
+            return self
+
     # Mathematical operations
     def __int__(self) -> int:
         if self.dimensionless:
diff --git a/pint/facets/plain/registry.py b/pint/facets/plain/registry.py
index 8572fec..7f2c9a7 100644
--- a/pint/facets/plain/registry.py
+++ b/pint/facets/plain/registry.py
@@ -272,6 +272,9 @@ class PlainRegistry(metaclass=RegistryMeta):
         #: Might contain prefixed units.
         self._units: Dict[str, UnitDefinition] = {}
 
+        #: List base unit names
+        self._base_units: List[str] = []
+
         #: Map unit name in lower case (string) to a set of unit names with the right
         #: case.
         #: Does not contain prefixed units.
@@ -453,6 +456,8 @@ class PlainRegistry(metaclass=RegistryMeta):
             # For a plain units, we need to define the related dimension
             # (making sure there is only one to define)
             if definition.is_base:
+                self._base_units.append(definition.name)
+
                 for dimension in definition.reference.keys():
                     if dimension in self._dimensions:
                         if dimension != "[]":
diff --git a/pint/testsuite/helpers.py b/pint/testsuite/helpers.py
index d72b5a3..ca0e978 100644
--- a/pint/testsuite/helpers.py
+++ b/pint/testsuite/helpers.py
@@ -10,6 +10,7 @@ from pint.testing import assert_equal as assert_quantity_equal  # noqa: F401
 
 from ..compat import (
     HAS_BABEL,
+    HAS_MIP,
     HAS_NUMPY,
     HAS_NUMPY_ARRAY_FUNCTION,
     HAS_UNCERTAINTIES,
@@ -139,6 +140,9 @@ requires_uncertainties = pytest.mark.skipif(
 requires_not_uncertainties = pytest.mark.skipif(
     HAS_UNCERTAINTIES, reason="Requires Uncertainties not to be installed."
 )
+requires_mip = pytest.mark.skipif(
+    not HAS_MIP, reason="Requires MIP"
+)
 
 # Parametrization
 
diff --git a/pint/testsuite/test_quantity.py b/pint/testsuite/test_quantity.py
index c9b3fe9..08b9a70 100644
--- a/pint/testsuite/test_quantity.py
+++ b/pint/testsuite/test_quantity.py
@@ -370,6 +370,44 @@ class TestQuantity(QuantityTestCase):
             round(abs(self.Q_("2 second").to("millisecond").magnitude - 2000), 7) == 0
         )
 
+    @helpers.requires_mip
+    def test_to_preferred(self):
+        ureg = UnitRegistry()
+        Q_ = ureg.Quantity
+
+        ureg.define("pound_force_per_square_foot = 47.8803 pascals = psf")
+        ureg.define("pound_mass = 0.45359237 kg = lbm")
+
+        preferred_units = [
+            ureg.ft,                    # distance      L
+            ureg.slug,                  # mass          M
+            ureg.s,                     # duration      T
+            ureg.rankine,               # temperature   Θ
+            ureg.lbf,                   # force         L M T^-2
+            ureg.psf,                   # pressure      M L^−1 T^−2
+            ureg.lbm * ureg.ft**-3,     # density       M L^-3
+            ureg.W                      # power         L^2 M T^-3
+        ]
+
+        temp = (Q_("1 lbf") * Q_("1 m/s")).to_preferred(preferred_units)
+        assert temp.units == ureg.W
+
+        temp = (Q_(" 1 lbf*m")).to_preferred(preferred_units)
+        # would prefer this to be repeatable, but mip doesn't guarantee that currently
+        assert temp.units in [ureg.W * ureg.s, ureg.ft * ureg.lbf]
+
+        temp = Q_("1 kg").to_preferred(preferred_units)
+        assert temp.units == ureg.slug
+
+        result = Q_("1 slug/m**3").to_preferred(preferred_units)
+        assert result.units == ureg.lbm * ureg.ft**-3
+
+        result = Q_("1 amp").to_preferred(preferred_units)
+        assert result.units == ureg.amp
+
+        result = Q_("1 volt").to_preferred(preferred_units)
+        assert result.units == ureg.volts
+
     @helpers.requires_numpy
     def test_convert_numpy(self):