diff options
Diffstat (limited to 'pint/facets/plain/quantity.py')
| -rw-r--r-- | pint/facets/plain/quantity.py | 228 |
1 files changed, 228 insertions, 0 deletions
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: |