File: qto.py

package info (click to toggle)
python-pint 0.24.4-2
  • links: PTS, VCS
  • area: main
  • in suites: trixie
  • size: 1,924 kB
  • sloc: python: 19,983; makefile: 149
file content (424 lines) | stat: -rw-r--r-- 14,690 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
from __future__ import annotations

import bisect
import math
import numbers
import warnings
from typing import TYPE_CHECKING

from ...compat import (
    mip_INF,
    mip_INTEGER,
    mip_Model,
    mip_model,
    mip_OptimizationStatus,
    mip_xsum,
)
from ...errors import UndefinedBehavior
from ...util import infer_base_unit

if TYPE_CHECKING:
    from ..._typing import UnitLike
    from ...util import UnitsContainer
    from .quantity import PlainQuantity


def _get_reduced_units(
    quantity: PlainQuantity, units: UnitsContainer
) -> UnitsContainer:
    # loop through individual units and compare to each other unit
    # can we do better than a nested loop here?
    for unit1, exp in units.items():
        # make sure it wasn't already reduced to zero exponent on prior pass
        if unit1 not in units:
            continue
        for unit2 in units:
            # get exponent after reduction
            exp = units[unit1]
            if unit1 != unit2:
                power = quantity._REGISTRY._get_dimensionality_ratio(unit1, unit2)
                if power:
                    units = units.add(unit2, exp / power).remove([unit1])
                    break
    return units


def ito_reduced_units(quantity: PlainQuantity) -> None:
    """Return PlainQuantity scaled in place to reduced units, i.e. one unit per
    dimension. This will not reduce compound units (e.g., 'J/kg' will not
    be reduced to m**2/s**2), nor can it make use of contexts at this time.
    """

    # shortcuts in case we're dimensionless or only a single unit
    if quantity.dimensionless:
        return quantity.ito({})
    if len(quantity._units) == 1:
        return None

    units = quantity._units.copy()
    new_units = _get_reduced_units(quantity, units)

    return quantity.ito(new_units)


def to_reduced_units(
    quantity: PlainQuantity,
) -> PlainQuantity:
    """Return PlainQuantity scaled in place to reduced units, i.e. one unit per
    dimension. This will not reduce compound units (intentionally), nor
    can it make use of contexts at this time.
    """

    # shortcuts in case we're dimensionless or only a single unit
    if quantity.dimensionless:
        return quantity.to({})
    if len(quantity._units) == 1:
        return quantity

    units = quantity._units.copy()
    new_units = _get_reduced_units(quantity, units)

    return quantity.to(new_units)


def to_compact(
    quantity: PlainQuantity, unit: UnitsContainer | None = None
) -> PlainQuantity:
    """ "Return PlainQuantity rescaled to compact, human-readable units.

    To get output in terms of a different unit, use the unit parameter.


    Examples
    --------

    >>> import pint
    >>> ureg = pint.UnitRegistry()
    >>> (200e-9*ureg.s).to_compact()
    <Quantity(200.0, 'nanosecond')>
    >>> (1e-2*ureg('kg m/s^2')).to_compact('N')
    <Quantity(10.0, 'millinewton')>
    """

    if not isinstance(quantity.magnitude, numbers.Number) and not hasattr(
        quantity.magnitude, "nominal_value"
    ):
        warnings.warn(
            "to_compact applied to non numerical types has an undefined behavior.",
            UndefinedBehavior,
            stacklevel=2,
        )
        return quantity

    if (
        quantity.unitless
        or quantity.magnitude == 0
        or math.isnan(quantity.magnitude)
        or math.isinf(quantity.magnitude)
    ):
        return quantity

    SI_prefixes: dict[int, str] = {}
    for prefix in quantity._REGISTRY._prefixes.values():
        try:
            scale = prefix.converter.scale
            # Kludgy way to check if this is an SI prefix
            log10_scale = int(math.log10(scale))
            if log10_scale == math.log10(scale):
                SI_prefixes[log10_scale] = prefix.name
        except Exception:
            SI_prefixes[0] = ""

    SI_prefixes_list = sorted(SI_prefixes.items())
    SI_powers = [item[0] for item in SI_prefixes_list]
    SI_bases = [item[1] for item in SI_prefixes_list]

    if unit is None:
        unit = infer_base_unit(quantity, registry=quantity._REGISTRY)
    else:
        unit = infer_base_unit(quantity.__class__(1, unit), registry=quantity._REGISTRY)

    q_base = quantity.to(unit)

    magnitude = q_base.magnitude
    # Support uncertainties
    if hasattr(magnitude, "nominal_value"):
        magnitude = magnitude.nominal_value

    units = list(q_base._units.items())
    units_numerator = [a for a in units if a[1] > 0]

    if len(units_numerator) > 0:
        unit_str, unit_power = units_numerator[0]
    else:
        unit_str, unit_power = units[0]

    if unit_power > 0:
        power = math.floor(math.log10(abs(magnitude)) / float(unit_power) / 3) * 3
    else:
        power = math.ceil(math.log10(abs(magnitude)) / float(unit_power) / 3) * 3

    index = bisect.bisect_left(SI_powers, power)

    if index >= len(SI_bases):
        index = -1

    prefix_str = SI_bases[index]

    new_unit_str = prefix_str + unit_str
    new_unit_container = q_base._units.rename(unit_str, new_unit_str)

    return quantity.to(new_unit_container)


def to_preferred(
    quantity: PlainQuantity, preferred_units: list[UnitLike] | None = None
) -> PlainQuantity:
    """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, 'watt * second')>
    """

    units = _get_preferred(quantity, preferred_units)
    return quantity.to(units)


def ito_preferred(
    quantity: PlainQuantity, preferred_units: list[UnitLike] | None = None
) -> PlainQuantity:
    """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, 'watt * second')>
    """

    units = _get_preferred(quantity, preferred_units)
    return quantity.ito(units)


def _get_preferred(
    quantity: PlainQuantity, preferred_units: list[UnitLike] | None = None
) -> PlainQuantity:
    if preferred_units is None:
        preferred_units = quantity._REGISTRY.default_preferred_units

    if not quantity.dimensionality:
        return quantity._units.copy()

    # 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(quantity.dimensionality)
        self_exps = [quantity.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 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(quantity._REGISTRY.Unit, quantity._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(quantity._REGISTRY.Unit, quantity._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(quantity._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(quantity.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(
        {u for d, u in unit_selections.items() if d in preferred_dims}
    )
    preferred_units.sort(key=str)  # for determinism

    # and unpreferred_units are the selected units that weren't originally preferred
    unpreferred_units = list(
        {u for d, u in unit_selections.items() if d not in preferred_dims}
    )
    unpreferred_units.sort(key=str)  # 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 = [quantity.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 (quantity)
    # 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 = quantity._REGISTRY.Unit("")
            for var in vars:
                if var.xi(i):
                    temp_unit *= quantity._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 result_unit

    # for whatever reason, a solution wasn't found
    # return the original quantity
    return quantity._units.copy()