Fixed pr/1574 conflicts

5 files changed, 294 insertions, 0 deletions

diff --git a/pint/compat.py b/pint/compat.py
index cbb60c7..de149ac 100644
--- a/pint/compat.py
+++ b/pint/compat.py
@@ -130,6 +130,20 @@ try:
 except ImportError:
     HAS_BABEL = False
 
+try:
+    import mip
+
+    mip_model = mip.model
+    mip_Model = mip.Model
+    mip_INF = mip.INF
+    mip_INTEGER = mip.INTEGER
+    mip_xsum = mip.xsum
+    mip_OptimizationStatus = mip.OptimizationStatus
+
+    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 +155,14 @@ else:
 if not HAS_BABEL:
     babel_parse = babel_units = missing_dependency("Babel")  # noqa: F811
 
+if not HAS_MIP:
+    mip_missing = missing_dependency("mip")
+    mip_model = mip_missing
+    mip_Model = mip_missing
+    mip_INF = mip_missing
+    mip_INTEGER = mip_missing
+    mip_xsum = mip_missing
+    mip_OptimizationStatus = mip_missing
 
 # Define location of pint.Quantity in NEP-13 type cast hierarchy by defining upcast
 # types using guarded imports
diff --git a/pint/facets/plain/quantity.py b/pint/facets/plain/quantity.py
index 4667c99..f4608c7 100644
--- a/pint/facets/plain/quantity.py
+++ b/pint/facets/plain/quantity.py
@@ -40,6 +40,12 @@ from ...compat import (
     eq,
     is_duck_array_type,
     is_upcast_type,
+    mip_INF,
+    mip_INTEGER,
+    mip_model,
+    mip_Model,
+    mip_OptimizationStatus,
+    mip_xsum,
     np,
     zero_or_nan,
 )
@@ -695,6 +701,228 @@ 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()
+        model.verbose = 0
+
+        # 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 628701e..94cc5ff 100644
--- a/pint/facets/plain/registry.py
+++ b/pint/facets/plain/registry.py
@@ -266,6 +266,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.
@@ -501,6 +504,7 @@ class PlainRegistry(metaclass=RegistryMeta):
 
     def _add_unit(self, definition: UnitDefinition):
         if definition.is_base:
+            self._base_units.append(definition.name)
             for dim_name in definition.reference.keys():
                 if dim_name not in self._dimensions:
                     self._add_dimension(DimensionDefinition(dim_name))
diff --git a/pint/testsuite/helpers.py b/pint/testsuite/helpers.py
index 0348a45..4c560fb 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,7 @@ 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 b1f15e5..8fb712a 100644
--- a/pint/testsuite/test_quantity.py
+++ b/pint/testsuite/test_quantity.py
@@ -369,6 +369,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):
         # Conversions with single units take a different codepath than