New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Considerations on the performance of astropy.units #7438

Closed
Juanlu001 opened this Issue May 4, 2018 · 6 comments

Comments

Projects
3 participants
@Juanlu001
Copy link
Member

Juanlu001 commented May 4, 2018

(This should have come before #7434, which I just closed)

A couple of years ago I was trying to use Quantity-enabled functions in poliastro that would be called in a for loop, and I was struck by some performance issues (see commit and corresponding pull request for a history of my attempts):

poliastro/poliastro@ed931ad

Yesterday at #pyastro18 @manodeep @eteq @adrn @duncanmmacleod, Bas Swinkles and myself gave a closer look on that and did some benchmarking. The first thing we tried was running a very typical poliastro program with line_profiler, you can check out the code here:

https://gist.github.com/Juanlu001/25388c43b4609c78d28e3265e2b8b2b7

(this of course requires poliastro to run)

This is an sample of the profiling results that we got, using latest stable version of astropy:

Wrote profile results to ex1.py.lprof
Timer unit: 1e-06 s

Total time: 2.62104 s
File: /home/juanlu/Development/Astropy/astropy-latest/astropy/units/core.py
Function: __div__ at line 642

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   642                                               @profile
   643                                               def __div__(self, m):
   644     17044      19208.0      1.1      0.7          if isinstance(m, (bytes, str)):
   645                                                       m = Unit(m)
   646                                           
   647     17044      14429.0      0.8      0.6          if isinstance(m, UnitBase):
   648     17036      66578.0      3.9      2.5              if m.is_unity():
   649         8          5.0      0.6      0.0                  return self
   650     17028    2519956.0    148.0     96.1              return CompositeUnit(1, [self, m], [1, -1], _error_check=False)
   651                                           
   652         8          2.0      0.2      0.0          try:
   653                                                       # Cannot handle this as Unit, re-try as Quantity
   654         8         54.0      6.8      0.0              from .quantity import Quantity
   655         8        805.0    100.6      0.0              return Quantity(1, self) / m
   656                                                   except TypeError:
   657                                                       return NotImplemented

Total time: 0.000511 s
File: /home/juanlu/Development/Astropy/astropy-latest/astropy/units/core.py
Function: __rdiv__ at line 659

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   659                                               @profile
   660                                               def __rdiv__(self, m):
   661         2          4.0      2.0      0.8          if isinstance(m, (bytes, str)):
   662                                                       return Unit(m) / self
   663                                           
   664         2          2.0      1.0      0.4          try:
   665                                                       # Cannot handle this as Unit.  Here, m cannot be a Quantity,
   666                                                       # so we make it into one, fasttracking when it does not have a
   667                                                       # unit, for the common case of <array> / <unit>.
   668         2         14.0      7.0      2.7              from .quantity import Quantity
   669         2          3.0      1.5      0.6              if hasattr(m, 'unit'):
   670                                                           result = Quantity(m)
   671                                                           result /= self
   672                                                           return result
   673                                                       else:
   674         2        488.0    244.0     95.5                  return Quantity(m, self**(-1))
   675                                                   except TypeError:
   676                                                       return NotImplemented

Total time: 0.032247 s
File: /home/juanlu/Development/Astropy/astropy-latest/astropy/units/core.py
Function: __mul__ at line 682

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   682                                               @profile
   683                                               def __mul__(self, m):
   684       106        127.0      1.2      0.4          if isinstance(m, (bytes, str)):
   685                                                       m = Unit(m)
   686                                           
   687       106         95.0      0.9      0.3          if isinstance(m, UnitBase):
   688       105       6463.0     61.6     20.0              if m.is_unity():
   689         2          1.0      0.5      0.0                  return self
   690       103       5210.0     50.6     16.2              elif self.is_unity():
   691                                                           return m
   692       103      20203.0    196.1     62.7              return CompositeUnit(1, [self, m], [1, 1], _error_check=False)
   693                                           
   694                                                   # Cannot handle this as Unit, re-try as Quantity.
   695         1          0.0      0.0      0.0          try:
   696         1          9.0      9.0      0.0              from .quantity import Quantity
   697         1        139.0    139.0      0.4              return Quantity(1, self) * m
   698                                                   except TypeError:
   699                                                       return NotImplemented

Total time: 3.84527 s
File: /home/juanlu/Development/Astropy/astropy-latest/astropy/units/core.py
Function: __rmul__ at line 701

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   701                                               @profile
   702                                               def __rmul__(self, m):
   703     50183      49549.0      1.0      1.3          if isinstance(m, (bytes, str)):
   704                                                       return Unit(m) * self
   705                                           
   706                                                   # Cannot handle this as Unit.  Here, m cannot be a Quantity,
   707                                                   # so we make it into one, fasttracking when it does not have a unit
   708                                                   # for the common case of <array> * <unit>.
   709     50183      24840.0      0.5      0.6          try:
   710     50183     272694.0      5.4      7.1              from .quantity import Quantity
   711     50183      49806.0      1.0      1.3              if hasattr(m, 'unit'):
   712                                                           result = Quantity(m)
   713                                                           result *= self
   714                                                           return result
   715                                                       else:
   716     50183    3448380.0     68.7     89.7                  return Quantity(m, self)
   717                                                   except TypeError:
   718                                                       return NotImplemented

Total time: 4.37844 s
File: /home/juanlu/Development/Astropy/astropy-latest/astropy/units/core.py
Function: __init__ at line 1973

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
  1973                                               @profile
  1974                                               def __init__(self, scale, bases, powers, decompose=False,
  1975                                                            decompose_bases=set(), _error_check=True):
  1976                                                   # There are many cases internal to astropy.units where we
  1977                                                   # already know that all the bases are Unit objects, and the
  1978                                                   # powers have been validated.  In those cases, we can skip the
  1979                                                   # error checking for performance reasons.  When the private
  1980                                                   # kwarg `_error_check` is False, the error checking is turned
  1981                                                   # off.
  1982     36326      27007.0      0.7      0.6          if _error_check:
  1983     17943      50541.0      2.8      1.2              scale = sanitize_scale(scale)
  1984     35980      26243.0      0.7      0.6              for base in bases:
  1985     18037      20410.0      1.1      0.5                  if not isinstance(base, UnitBase):
  1986                                                               raise TypeError(
  1987                                                                   "bases must be sequence of UnitBase instances")
  1988     17943     288092.0     16.1      6.6              powers = [validate_power(p) for p in powers]
  1989                                           
  1990     36326      29810.0      0.8      0.7          self._scale = scale
  1991     36326      23310.0      0.6      0.5          self._bases = bases
  1992     36326      22771.0      0.6      0.5          self._powers = powers
  1993     36326      23957.0      0.7      0.5          self._decomposed_cache = None
  1994     36326    3839278.0    105.7     87.7          self._expand_and_gather(decompose=decompose, bases=decompose_bases)
  1995     36326      27022.0      0.7      0.6          self._hash = None

Total time: 2.71706 s
File: /home/juanlu/Development/Astropy/astropy-latest/astropy/units/core.py
Function: _expand_and_gather at line 2036

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
  2036                                               @profile
  2037                                               def _expand_and_gather(self, decompose=False, bases=set()):
  2038     36326      56151.0      1.5      2.1          def add_unit(unit, power, scale):
  2039                                                       if unit not in bases:
  2040                                                           for base in bases:
  2041                                                               try:
  2042                                                                   scale *= unit._to(base) ** power
  2043                                                               except UnitsError:
  2044                                                                   pass
  2045                                                               else:
  2046                                                                   unit = base
  2047                                                                   break
  2048                                           
  2049                                                       if unit in new_parts:
  2050                                                           a, b = resolve_fractions(new_parts[unit], power)
  2051                                                           new_parts[unit] = a + b
  2052                                                       else:
  2053                                                           new_parts[unit] = power
  2054                                                       return scale
  2055                                           
  2056     36326      44891.0      1.2      1.7          new_parts = {}
  2057     36326      61148.0      1.7      2.3          scale = self.scale
  2058                                           
  2059     89876     166541.0      1.9      6.1          for b, p in zip(self.bases, self.powers):
  2060     53550      63961.0      1.2      2.4              if decompose and b not in bases:
  2061      1137       5277.0      4.6      0.2                  b = b.decompose(bases=bases)
  2062                                           
  2063     53550      79059.0      1.5      2.9              if isinstance(b, CompositeUnit):
  2064     18086      30420.0      1.7      1.1                  scale *= b._scale ** p
  2065     37606      56004.0      1.5      2.1                  for b_sub, p_sub in zip(b._bases, b._powers):
  2066     19520     327556.0     16.8     12.1                      a, b = resolve_fractions(p_sub, p)
  2067     19520     181461.0      9.3      6.7                      scale = add_unit(b_sub, a * b, scale)
  2068                                                       else:
  2069     35464     346411.0      9.8     12.7                  scale = add_unit(b, p, scale)
  2070                                           
  2071     36326     106591.0      2.9      3.9          new_parts = [x for x in new_parts.items() if x[1] != 0]
  2072     36326     145916.0      4.0      5.4          new_parts.sort(key=lambda x: (-x[1], getattr(x[0], 'name', '')))
  2073                                           
  2074     36326      83529.0      2.3      3.1          self._bases = [x[0] for x in new_parts]
  2075     36326     831713.0     22.9     30.6          self._powers = [validate_power(x[1]) for x in new_parts]
  2076     36326     130427.0      3.6      4.8          self._scale = sanitize_scale(scale)

Total time: 1.79216 s
File: /home/juanlu/Development/Astropy/astropy-latest/astropy/units/quantity.py
Function: __new__ at line 274

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   274                                               @profile
   275                                               def __new__(cls, value, unit=None, dtype=None, copy=True, order=None,
   276                                                           subok=False, ndmin=0):
   277                                           
   278     50384      64913.0      1.3      3.6          if unit is not None:
   279                                                       # convert unit first, to avoid multiple string->unit conversions
   280     50236     163725.0      3.3      9.1              unit = Unit(unit)
   281                                                       # if we allow subclasses, allow a class from the unit.
   282     50236      55267.0      1.1      3.1              if subok:
   283                                                           qcls = getattr(unit, '_quantity_class', cls)
   284                                                           if issubclass(qcls, cls):
   285                                                               cls = qcls
   286                                           
   287                                                   # optimize speed for Quantity with no dtype given, copy=False
   288     50384      75685.0      1.5      4.2          if isinstance(value, Quantity):
   289       156        180.0      1.2      0.0              if unit is not None and unit is not value.unit:
   290                                                           value = value.to(unit)
   291                                                           # the above already makes a copy (with float dtype)
   292                                                           copy = False
   293                                           
   294       156        195.0      1.2      0.0              if type(value) is not cls and not (subok and
   295                                                                                          isinstance(value, cls)):
   296        10        540.0     54.0      0.0                  value = value.view(cls)
   297                                           
   298       156        178.0      1.1      0.0              if dtype is None:
   299       156        163.0      1.0      0.0                  if not copy:
   300        20         21.0      1.1      0.0                      return value
   301                                           
   302       136        264.0      1.9      0.0                  if not np.can_cast(np.float32, value.dtype):
   303                                                               dtype = float
   304                                           
   305       136        158.0      1.2      0.0              return np.array(value, dtype=dtype, copy=copy, order=order,
   306       136       1401.0     10.3      0.1                              subok=True, ndmin=ndmin)
   307                                           
   308                                                   # Maybe str, or list/tuple of Quantity? If so, this may set value_unit.
   309                                                   # To ensure array remains fast, we short-circuit it.
   310     50228      55928.0      1.1      3.1          value_unit = None
   311     50228      71444.0      1.4      4.0          if not isinstance(value, np.ndarray):
   312     16794      21856.0      1.3      1.2              if isinstance(value, str):
   313                                                           # The first part of the regex string matches any integer/float;
   314                                                           # the second parts adds possible trailing .+-, which will break
   315                                                           # the float function below and ensure things like 1.2.3deg
   316                                                           # will not work.
   317                                                           pattern = (r'\s*[+-]?'
   318                                                                      r'((\d+\.?\d*)|(\.\d+)|([nN][aA][nN])|'
   319                                                                      r'([iI][nN][fF]([iI][nN][iI][tT][yY]){0,1}))'
   320                                                                      r'([eE][+-]?\d+)?'
   321                                                                      r'[.+-]?')
   322                                           
   323                                                           v = re.match(pattern, value)
   324                                                           unit_string = None
   325                                                           try:
   326                                                               value = float(v.group())
   327                                           
   328                                                           except Exception:
   329                                                               raise TypeError('Cannot parse "{0}" as a {1}. It does not '
   330                                                                               'start with a number.'
   331                                                                               .format(value, cls.__name__))
   332                                           
   333                                                           unit_string = v.string[v.end():].strip()
   334                                                           if unit_string:
   335                                                               value_unit = Unit(unit_string)
   336                                                               if unit is None:
   337                                                                   unit = value_unit  # signal no conversion needed below.
   338                                           
   339     16794      68771.0      4.1      3.8              elif (isiterable(value) and len(value) > 0 and
   340         5         25.0      5.0      0.0                    all(isinstance(v, Quantity) for v in value)):
   341                                                           # Convert all quantities to the same unit.
   342                                                           if unit is None:
   343                                                               unit = value[0].unit
   344                                                           value = [q.to_value(unit) for q in value]
   345                                                           value_unit = unit  # signal below that conversion has been done
   346                                           
   347     50228      59891.0      1.2      3.3          if value_unit is None:
   348                                                       # If the value has a `unit` attribute and if not None
   349                                                       # (for Columns with uninitialized unit), treat it like a quantity.
   350     50228      84637.0      1.7      4.7              value_unit = getattr(value, 'unit', None)
   351     50228      57686.0      1.1      3.2              if value_unit is None:
   352                                                           # Default to dimensionless for no (initialized) unit attribute.
   353     50228      58737.0      1.2      3.3                  if unit is None:
   354         1          2.0      2.0      0.0                      unit = cls._default_unit
   355     50228      57463.0      1.1      3.2                  value_unit = unit  # signal below that no conversion is needed
   356                                                       else:
   357                                                           try:
   358                                                               value_unit = Unit(value_unit)
   359                                                           except Exception as exc:
   360                                                               raise TypeError("The unit attribute {0!r} of the input could "
   361                                                                               "not be parsed as an astropy Unit, raising "
   362                                                                               "the following exception:\n{1}"
   363                                                                               .format(value.unit, exc))
   364                                           
   365                                                           if unit is None:
   366                                                               unit = value_unit
   367                                                           elif unit is not value_unit:
   368                                                               copy = False  # copy will be made in conversion at end
   369                                           
   370     50228      63914.0      1.3      3.6          value = np.array(value, dtype=dtype, copy=copy, order=order,
   371     50228     192727.0      3.8     10.8                           subok=False, ndmin=ndmin)
   372                                           
   373                                                   # check that array contains numbers or long int objects
   374     50228      75752.0      1.5      4.2          if (value.dtype.kind in 'OSU' and
   375                                                       not (value.dtype.kind == 'O' and
   376                                                            isinstance(value.item(() if value.ndim == 0 else 0),
   377                                                                       numbers.Number))):
   378                                                       raise TypeError("The value must be a valid Python or "
   379                                                                       "Numpy numeric type.")
   380                                           
   381                                                   # by default, cast any integer, boolean, etc., to float
   382     50228     103574.0      2.1      5.8          if dtype is None and (not np.can_cast(np.float32, value.dtype)
   383     50181      66032.0      1.3      3.7                                or value.dtype.kind == 'O'):
   384        47        218.0      4.6      0.0              value = value.astype(float)
   385                                           
   386     50228     147187.0      2.9      8.2          value = value.view(cls)
   387     50228     131282.0      2.6      7.3          value._set_unit(value_unit)
   388     50228      60293.0      1.2      3.4          if unit is value_unit:
   389     50228      52054.0      1.0      2.9              return value
   390                                                   else:
   391                                                       # here we had non-Quantity input that had a "unit" attribute
   392                                                       # with a unit different from the desired one.  So, convert.
   393                                                       return value.to(unit)

Total time: 0.52432 s
File: /home/juanlu/Development/Astropy/astropy-latest/astropy/units/utils.py
Function: validate_power at line 191

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   191                                           @profile
   192                                           def validate_power(p, support_tuples=False):
   193                                               """Convert a power to a floating point value, an integer, or a Fraction.
   194                                           
   195                                               If a fractional power can be represented exactly as a floating point
   196                                               number, convert it to a float, to make the math much faster; otherwise,
   197                                               retain it as a `fractions.Fraction` object to avoid losing precision.
   198                                               Conversely, if the value is indistinguishable from a rational number with a
   199                                               low-numbered denominator, convert to a Fraction object.
   200                                           
   201                                               Parameters
   202                                               ----------
   203                                               p : float, int, Rational, Fraction
   204                                                   Power to be converted
   205                                               """
   206     72810     292997.0      4.0     55.9      if isinstance(p, (numbers.Rational, Fraction)):
   207     72767      58848.0      0.8     11.2          denom = p.denominator
   208     72767      52143.0      0.7      9.9          if denom == 1:
   209     72763      72233.0      1.0     13.8              p = int(p.numerator)
   210                                                   # This is bit-twiddling hack to see if the integer is a
   211                                                   # power of two
   212         4          4.0      1.0      0.0          elif (denom & (denom - 1)) == 0:
   213         4         13.0      3.2      0.0              p = float(p)
   214                                               else:
   215        43         33.0      0.8      0.0          try:
   216        43         38.0      0.9      0.0              p = float(p)
   217                                                   except Exception:
   218                                                       if not np.isscalar(p):
   219                                                           raise ValueError("Quantities and Units may only be raised "
   220                                                                            "to a scalar power")
   221                                                       else:
   222                                                           raise
   223                                           
   224        43         42.0      1.0      0.0          if (p % 1.0) == 0.0:
   225                                                       # Denominators of 1 can just be integers.
   226         3          3.0      1.0      0.0              p = int(p)
   227        40         36.0      0.9      0.0          elif (p * 8.0) % 1.0 == 0.0:
   228                                                       # Leave alone if the denominator is exactly 2, 4 or 8, since this
   229                                                       # can be perfectly represented as a float, which means subsequent
   230                                                       # operations are much faster.
   231        40         26.0      0.7      0.0              pass
   232                                                   else:
   233                                                       # Convert floats indistinguishable from a rational to Fraction.
   234                                                       # Here, we do not need to test values that are divisors of a higher
   235                                                       # number, such as 3, since it is already addressed by 6.
   236                                                       for i in (10, 9, 7, 6):
   237                                                           scaled = p * float(i)
   238                                                           if((scaled + 4. * _float_finfo.eps) % 1.0 <
   239                                                              8. * _float_finfo.eps):
   240                                                               p = Fraction(int(round(scaled)), i)
   241                                                               break
   242                                           
   243     72810      47904.0      0.7      9.1      return p

Total time: 0.228974 s
File: /home/juanlu/Development/Astropy/astropy-latest/astropy/units/utils.py
Function: resolve_fractions at line 246

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   246                                           @profile
   247                                           def resolve_fractions(a, b):
   248                                               """
   249                                               If either input is a Fraction, convert the other to a Fraction.
   250                                               This ensures that any operation involving a Fraction will use
   251                                               rational arithmetic and preserve precision.
   252                                               """
   253     19659     114254.0      5.8     49.9      a_is_fraction = isinstance(a, Fraction)
   254     19659      89888.0      4.6     39.3      b_is_fraction = isinstance(b, Fraction)
   255     19659       7723.0      0.4      3.4      if a_is_fraction and not b_is_fraction:
   256                                                   b = Fraction(b)
   257     19659       9270.0      0.5      4.0      elif not a_is_fraction and b_is_fraction:
   258                                                   a = Fraction(a)
   259     19659       7839.0      0.4      3.4      return a, b

Total time: 0.077689 s
File: ex1.py
Function: accel at line 15

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
    15                                           @profile
    16                                           def accel(t0, state, k):
    17      8343       7838.0      0.9     10.1      v_vec = state[3:]
    18      8343      40963.0      4.9     52.7      norm_v = (v_vec * v_vec).sum() ** .5
    19      8343      28888.0      3.5     37.2      return 1e-5 * v_vec / norm_v

Total time: 7.09469 s
File: ex1.py
Function: accel_slow at line 21

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
    21                                           @profile
    22                                           def accel_slow(t0, state, k):
    23      8343      16870.0      2.0      0.2      r_vec, v_vec = state[:3], state[3:]
    24      8343    4201150.0    503.6     59.2      _k = k * (u.km ** 3 / u.s ** 2)
    25      8343      29063.0      3.5      0.4      body = Body(None, _k, "_Dummy")
    26      8343     689695.0     82.7      9.7      _r = r_vec * u.km
    27      8343    1964005.0    235.4     27.7      _v = v_vec * u.km / u.s
    28      8343      34757.0      4.2      0.5      ss = RVState(body, _r, _v)
    29      8343     102100.0     12.2      1.4      norm_v = (v_vec * v_vec).sum() ** .5
    30      8343      57051.0      6.8      0.8      return 1e-5 * v_vec / norm_v

Total time: 2.23195 s
File: ex1.py
Function: accel_so_so at line 36

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
    36                                           @profile
    37                                           def accel_so_so(t0, state, k):
    38      8343      13523.0      1.6      0.6      r_vec, v_vec = state[:3], state[3:]
    39      8343     735756.0     88.2     33.0      _k = k * km3s2
    40      8343      25122.0      3.0      1.1      body = Body(None, _k, "_Dummy")
    41      8343     663202.0     79.5     29.7      _r = r_vec * u.km
    42      8343     648396.0     77.7     29.1      _v = v_vec * kms
    43      8343      26060.0      3.1      1.2      ss = RVState(body, _r, _v)
    44      8343      76216.0      9.1      3.4      norm_v = (v_vec * v_vec).sum() ** .5
    45      8343      43672.0      5.2      2.0      return 1e-5 * v_vec / norm_v

Total time: 9.93957 s
File: ex1.py
Function: main at line 49

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
    49                                           @profile
    50                                           def main():
    51         1     188966.0 188966.0      1.9      initial.propagate(3 * u.day, method=cowell, ad=accel)
    52         1    7327508.0 7327508.0     73.7      initial.propagate(3 * u.day, method=cowell, ad=accel_slow)
    53         1    2423094.0 2423094.0     24.4      initial.propagate(3 * u.day, method=cowell, ad=accel_so_so)

So here are our main takeaways:

  • With this poliastro typical use case, the naïve user-defined function with units is ~100x slower than the unitless one
    • Taking into account the whole main program, the factor is ~40x
  • From our experiments, the Quantity and CompositeUnit initializers are the biggest contributors to the slowdown, as they have significant overhead. Specifically, we're talking about Quantity.__new__ and CompositeUnit.__init__.
  • Pre-computing the composite units outside of the function so they would not have to be re-computed on each iteration (such as km3s2 = u.km ** 3 / u.s ** 2) made the naïve version ~3x faster.
  • Most of the time of CompositeUnit.__init__ is spent in CompositeUnit._expand_and_gather. The functions called by the latter (mainly resolve_fractions and validate_power) make a few isinstance calls that are extremely slow and represent a significant part of the total running time.
  • Quantity.__new__, on the other hand, has more balanced running time throughout all the lines, although it should be noted that it's considerably long and it also makes heavy use of isinstance (slow). The reason is that many different ways to initialize a Quantity are supported.
  • As the creation of Quantity objects is slow, all the operations are affected. This can be noted in the performance analysis of UnitBase.__rmul__, for example.

Proposals and ideas:

  • In #7434 we explored the possibility of doing a extreme simplification of Quantity.__new__ and drop some APIs. This lead to an improvement of 2x in total running time of the benchmark above, both the naïve and the smart version. Notice though that this hack broke many tests and that it still requires lots of effort before getting to a state where it can be reasonably discussed.
  • There is some low-hanging fruit that was explored in these commits, but that only amounted to a small contribution on the total running time.
  • We also tried Pint, and it's not faster than astropy.units. In fact, it was much slower.
  • Pint has these tips to improve performance http://pint.readthedocs.io/en/latest/performance.html but they do not solve all the problems, since the conversion has still to be done on input and output.
@mhvk

This comment has been minimized.

Copy link
Contributor

mhvk commented May 4, 2018

@Juanlu001 - thanks for the good summary. One suggestion off-the-bat for faster initialization: rather than construct a new class, I would suggest class methods for specific types of input (say, from_quantity, from_ndarray, from_array_like, from_string, and from_quantity_sequence). Indeed, __new__ could then use these.

Also one comment on the test code: I think proper use of quantities means that all routines already get their input as quantities, so in that sense the slow case is not really representative. (Of course, that is also why you call it a "naive" version, but at the very least one should use u.Quantity(input, desired_unit), not input * desired_unit which gives the wrong unit for quantity input).

@Juanlu001

This comment has been minimized.

Copy link
Member

Juanlu001 commented May 4, 2018

One suggestion off-the-bat for faster initialization: ...

Those are great ideas!

Also one comment on the test code: ...

For this particular case, scipy.integrate.ode is calling the user-defined function on each iteration with a numpy array as an input, so if I want a quantity I have to create it (which I usually do with a decorator). So there's not much I can do about it.

@mhvk

This comment has been minimized.

Copy link
Contributor

mhvk commented May 4, 2018

OK.

For the particular case you mention, I think you know you're getting ndarray, so explicit u.Quantity(input, unit, copy=False) with predefined units should be about as fast as can be done via the standard initializer. Fastest possible would be

q = input.view(u.Quantity)
q._unit = unit

I mention this in part to give some benchmark, and since I think an important part of this will be to write more general documentation about how to use units with other packages, what the pitfalls are, and how to optimize performance (and fitting is definitely a major one! though perhaps also a case where, if at all possible, one would want to do unit conversion up front, and then simply works with arrays). That would of course be needed also if we have more class methods for initializing a Quantity.

@pllim pllim added this to Needs decision in 3.1 Feature Planning May 5, 2018

@pllim

This comment has been minimized.

Copy link
Member

pllim commented May 8, 2018

an important part of this will be to write more general documentation about how to use units with other packages, what the pitfalls are, and how to optimize performance

💯

@Juanlu001

This comment has been minimized.

Copy link
Member

Juanlu001 commented Aug 7, 2018

Noting that these results should be re evaluated after #7649 lands (also including the work of recent pull requests focusing on performance)

@Juanlu001

This comment has been minimized.

Copy link
Member

Juanlu001 commented Oct 30, 2018

@mhvk 's #7649 changed many things, so at the very least we should repeat the analysis in case this is still a bottleneck. According to the comments there, the numbers improved a lot 😁 Closing!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment