util.py 9.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256
  1. """
  2. Several methods to simplify expressions involving unit objects.
  3. """
  4. from functools import reduce
  5. from collections.abc import Iterable
  6. from typing import Optional
  7. from sympy import default_sort_key
  8. from sympy.core.add import Add
  9. from sympy.core.containers import Tuple
  10. from sympy.core.mul import Mul
  11. from sympy.core.power import Pow
  12. from sympy.core.sorting import ordered
  13. from sympy.core.sympify import sympify
  14. from sympy.matrices.common import NonInvertibleMatrixError
  15. from sympy.physics.units.dimensions import Dimension, DimensionSystem
  16. from sympy.physics.units.prefixes import Prefix
  17. from sympy.physics.units.quantities import Quantity
  18. from sympy.physics.units.unitsystem import UnitSystem
  19. from sympy.utilities.iterables import sift
  20. def _get_conversion_matrix_for_expr(expr, target_units, unit_system):
  21. from sympy.matrices.dense import Matrix
  22. dimension_system = unit_system.get_dimension_system()
  23. expr_dim = Dimension(unit_system.get_dimensional_expr(expr))
  24. dim_dependencies = dimension_system.get_dimensional_dependencies(expr_dim, mark_dimensionless=True)
  25. target_dims = [Dimension(unit_system.get_dimensional_expr(x)) for x in target_units]
  26. canon_dim_units = [i for x in target_dims for i in dimension_system.get_dimensional_dependencies(x, mark_dimensionless=True)]
  27. canon_expr_units = set(dim_dependencies)
  28. if not canon_expr_units.issubset(set(canon_dim_units)):
  29. return None
  30. seen = set()
  31. canon_dim_units = [i for i in canon_dim_units if not (i in seen or seen.add(i))]
  32. camat = Matrix([[dimension_system.get_dimensional_dependencies(i, mark_dimensionless=True).get(j, 0) for i in target_dims] for j in canon_dim_units])
  33. exprmat = Matrix([dim_dependencies.get(k, 0) for k in canon_dim_units])
  34. try:
  35. res_exponents = camat.solve(exprmat)
  36. except NonInvertibleMatrixError:
  37. return None
  38. return res_exponents
  39. def convert_to(expr, target_units, unit_system="SI"):
  40. """
  41. Convert ``expr`` to the same expression with all of its units and quantities
  42. represented as factors of ``target_units``, whenever the dimension is compatible.
  43. ``target_units`` may be a single unit/quantity, or a collection of
  44. units/quantities.
  45. Examples
  46. ========
  47. >>> from sympy.physics.units import speed_of_light, meter, gram, second, day
  48. >>> from sympy.physics.units import mile, newton, kilogram, atomic_mass_constant
  49. >>> from sympy.physics.units import kilometer, centimeter
  50. >>> from sympy.physics.units import gravitational_constant, hbar
  51. >>> from sympy.physics.units import convert_to
  52. >>> convert_to(mile, kilometer)
  53. 25146*kilometer/15625
  54. >>> convert_to(mile, kilometer).n()
  55. 1.609344*kilometer
  56. >>> convert_to(speed_of_light, meter/second)
  57. 299792458*meter/second
  58. >>> convert_to(day, second)
  59. 86400*second
  60. >>> 3*newton
  61. 3*newton
  62. >>> convert_to(3*newton, kilogram*meter/second**2)
  63. 3*kilogram*meter/second**2
  64. >>> convert_to(atomic_mass_constant, gram)
  65. 1.660539060e-24*gram
  66. Conversion to multiple units:
  67. >>> convert_to(speed_of_light, [meter, second])
  68. 299792458*meter/second
  69. >>> convert_to(3*newton, [centimeter, gram, second])
  70. 300000*centimeter*gram/second**2
  71. Conversion to Planck units:
  72. >>> convert_to(atomic_mass_constant, [gravitational_constant, speed_of_light, hbar]).n()
  73. 7.62963087839509e-20*hbar**0.5*speed_of_light**0.5/gravitational_constant**0.5
  74. """
  75. from sympy.physics.units import UnitSystem
  76. unit_system = UnitSystem.get_unit_system(unit_system)
  77. if not isinstance(target_units, (Iterable, Tuple)):
  78. target_units = [target_units]
  79. if isinstance(expr, Add):
  80. return Add.fromiter(convert_to(i, target_units, unit_system)
  81. for i in expr.args)
  82. expr = sympify(expr)
  83. target_units = sympify(target_units)
  84. if not isinstance(expr, Quantity) and expr.has(Quantity):
  85. expr = expr.replace(lambda x: isinstance(x, Quantity),
  86. lambda x: x.convert_to(target_units, unit_system))
  87. def get_total_scale_factor(expr):
  88. if isinstance(expr, Mul):
  89. return reduce(lambda x, y: x * y,
  90. [get_total_scale_factor(i) for i in expr.args])
  91. elif isinstance(expr, Pow):
  92. return get_total_scale_factor(expr.base) ** expr.exp
  93. elif isinstance(expr, Quantity):
  94. return unit_system.get_quantity_scale_factor(expr)
  95. return expr
  96. depmat = _get_conversion_matrix_for_expr(expr, target_units, unit_system)
  97. if depmat is None:
  98. return expr
  99. expr_scale_factor = get_total_scale_factor(expr)
  100. return expr_scale_factor * Mul.fromiter(
  101. (1/get_total_scale_factor(u)*u)**p for u, p in
  102. zip(target_units, depmat))
  103. def quantity_simplify(expr, across_dimensions: bool=False, unit_system=None):
  104. """Return an equivalent expression in which prefixes are replaced
  105. with numerical values and all units of a given dimension are the
  106. unified in a canonical manner by default. `across_dimensions` allows
  107. for units of different dimensions to be simplified together.
  108. `unit_system` must be specified if `across_dimensions` is True.
  109. Examples
  110. ========
  111. >>> from sympy.physics.units.util import quantity_simplify
  112. >>> from sympy.physics.units.prefixes import kilo
  113. >>> from sympy.physics.units import foot, inch, joule, coulomb
  114. >>> quantity_simplify(kilo*foot*inch)
  115. 250*foot**2/3
  116. >>> quantity_simplify(foot - 6*inch)
  117. foot/2
  118. >>> quantity_simplify(5*joule/coulomb, across_dimensions=True, unit_system="SI")
  119. 5*volt
  120. """
  121. if expr.is_Atom or not expr.has(Prefix, Quantity):
  122. return expr
  123. # replace all prefixes with numerical values
  124. p = expr.atoms(Prefix)
  125. expr = expr.xreplace({p: p.scale_factor for p in p})
  126. # replace all quantities of given dimension with a canonical
  127. # quantity, chosen from those in the expression
  128. d = sift(expr.atoms(Quantity), lambda i: i.dimension)
  129. for k in d:
  130. if len(d[k]) == 1:
  131. continue
  132. v = list(ordered(d[k]))
  133. ref = v[0]/v[0].scale_factor
  134. expr = expr.xreplace({vi: ref*vi.scale_factor for vi in v[1:]})
  135. if across_dimensions:
  136. # combine quantities of different dimensions into a single
  137. # quantity that is equivalent to the original expression
  138. if unit_system is None:
  139. raise ValueError("unit_system must be specified if across_dimensions is True")
  140. unit_system = UnitSystem.get_unit_system(unit_system)
  141. dimension_system: DimensionSystem = unit_system.get_dimension_system()
  142. dim_expr = unit_system.get_dimensional_expr(expr)
  143. dim_deps = dimension_system.get_dimensional_dependencies(dim_expr, mark_dimensionless=True)
  144. target_dimension: Optional[Dimension] = None
  145. for ds_dim, ds_dim_deps in dimension_system.dimensional_dependencies.items():
  146. if ds_dim_deps == dim_deps:
  147. target_dimension = ds_dim
  148. break
  149. if target_dimension is None:
  150. # if we can't find a target dimension, we can't do anything. unsure how to handle this case.
  151. return expr
  152. target_unit = unit_system.derived_units.get(target_dimension)
  153. if target_unit:
  154. expr = convert_to(expr, target_unit, unit_system)
  155. return expr
  156. def check_dimensions(expr, unit_system="SI"):
  157. """Return expr if units in addends have the same
  158. base dimensions, else raise a ValueError."""
  159. # the case of adding a number to a dimensional quantity
  160. # is ignored for the sake of SymPy core routines, so this
  161. # function will raise an error now if such an addend is
  162. # found.
  163. # Also, when doing substitutions, multiplicative constants
  164. # might be introduced, so remove those now
  165. from sympy.physics.units import UnitSystem
  166. unit_system = UnitSystem.get_unit_system(unit_system)
  167. def addDict(dict1, dict2):
  168. """Merge dictionaries by adding values of common keys and
  169. removing keys with value of 0."""
  170. dict3 = {**dict1, **dict2}
  171. for key, value in dict3.items():
  172. if key in dict1 and key in dict2:
  173. dict3[key] = value + dict1[key]
  174. return {key:val for key, val in dict3.items() if val != 0}
  175. adds = expr.atoms(Add)
  176. DIM_OF = unit_system.get_dimension_system().get_dimensional_dependencies
  177. for a in adds:
  178. deset = set()
  179. for ai in a.args:
  180. if ai.is_number:
  181. deset.add(())
  182. continue
  183. dims = []
  184. skip = False
  185. dimdict = {}
  186. for i in Mul.make_args(ai):
  187. if i.has(Quantity):
  188. i = Dimension(unit_system.get_dimensional_expr(i))
  189. if i.has(Dimension):
  190. dimdict = addDict(dimdict, DIM_OF(i))
  191. elif i.free_symbols:
  192. skip = True
  193. break
  194. dims.extend(dimdict.items())
  195. if not skip:
  196. deset.add(tuple(sorted(dims, key=default_sort_key)))
  197. if len(deset) > 1:
  198. raise ValueError(
  199. "addends have incompatible dimensions: {}".format(deset))
  200. # clear multiplicative constants on Dimensions which may be
  201. # left after substitution
  202. reps = {}
  203. for m in expr.atoms(Mul):
  204. if any(isinstance(i, Dimension) for i in m.args):
  205. reps[m] = m.func(*[
  206. i for i in m.args if not i.is_number])
  207. return expr.xreplace(reps)