quantity.py 31.6 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
#!/bin/env python
"""
Module simtk.unit.quantity

Physical quantities with units, intended to produce similar functionality
to Boost.Units package in C++ (but with a runtime cost).
Uses similar API as Scientific.Physics.PhysicalQuantities
but different internals to satisfy our local requirements.
In particular, there is no underlying set of 'canonical' base
units, whereas in Scientific.Physics.PhysicalQuantities all
units are secretly in terms of SI units.  Also, it is easier
to add new fundamental dimensions to simtk.dimensions.  You
Justin MacCallum's avatar
Justin MacCallum committed
13
might want to make new dimensions for, say, "currency" or
14
15
16
17
18
19
20
21
22
23
24
25
"information".

Some features of this implementation:
  * Quantities are a combination of a value and a unit.  The value
    part can be any python type, including numbers, lists, numpy
    arrays, and anything else.  The unit part must be a simtk.unit.Unit.
  * Operations like adding incompatible units raises an error.
  * Multiplying or dividing units/quantities creates new units.
  * Users can create new Units and Dimensions, but most of the useful
    ones are predefined.
  * Conversion factors between units are applied transitively, so all
    possible conversions are available.
Justin MacCallum's avatar
Justin MacCallum committed
26
27
28
  * I want dimensioned Quantities that are compatible with numpy arrays,
    but do not necessarily require the python numpy package. In other
    words, Quantities can be based on either numpy arrays or on built in
29
    python types.
Justin MacCallum's avatar
Justin MacCallum committed
30
31
32
33
34
  * Units are NOT necessarily stored in terms of SI units internally.
    This is very important for me, because one important application
    area for us is at the molecular scale. Using SI units internally
    can lead to exponent overflow in commonly used molecular force
    calculations. Internally, all unit systems are equally fundamental
35
36
37
38
39
    in SimTK.

Two possible enhancements that have not been implemented are
  1) Include uncertainties with propagation of errors
  2) Incorporate offsets for celsius <-> kelvin conversion
40
41
42
43
44
45
46
47
48
49
50
51



This is part of the OpenMM molecular simulation toolkit originating from
Simbios, the NIH National Center for Physics-Based Simulation of
Biological Structures at Stanford, funded under the NIH Roadmap for
Medical Research, grant U54 GM072970. See https://simtk.org.

Portions copyright (c) 2012 Stanford University and the Authors.
Authors: Christopher M. Bruns
Contributors: Peter Eastman

Justin MacCallum's avatar
Justin MacCallum committed
52
Permission is hereby granted, free of charge, to any person obtaining a
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
copy of this software and associated documentation files (the "Software"),
to deal in the Software without restriction, including without limitation
the rights to use, copy, modify, merge, publish, distribute, sublicense,
and/or sell copies of the Software, and to permit persons to whom the
Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
USE OR OTHER DEALINGS IN THE SOFTWARE.
69
70
"""

Peter Eastman's avatar
Peter Eastman committed
71
72
from __future__ import division

73
74
75
76
77
78
79
80
81
82
83
__author__ = "Christopher M. Bruns"
__version__ = "0.5"


import math
import copy
from standard_dimensions import *
from unit import Unit, is_unit, dimensionless

class Quantity(object):
    """Physical quantity, such as 1.3 meters per second.
Justin MacCallum's avatar
Justin MacCallum committed
84

85
86
87
88
89
90
91
92
93
94
    Quantities contain both a value, such as 1.3; and a unit,
    such as 'meters per second'.

    Supported value types include:
      1 - numbers (float, int, long)
      2 - lists of numbers, e.g. [1,2,3]
      3 - tuples of numbers, e.g. (1,2,3)
            Note - unit conversions will cause tuples to be converted to lists
      4 - lists of tuples of numbers, lists of lists of ... etc. of numbers
      5 - numpy.arrays
Justin MacCallum's avatar
Justin MacCallum committed
95

96
97
98
99
100
101
102
103
104
105
106
    Create numpy.arrays with units using the Quantity constructor, not the
    multiply operator.  e.g.

      Quantity(numpy.array([1,2,3]), centimeters) # correct

        *NOT*

      numpy.array([1,2,3]) * centimeters # won't work

    because numpy.arrays already overload the multiply operator for EVERYTHING.
    """
Justin MacCallum's avatar
Justin MacCallum committed
107

108
109
110
    def __init__(self, value=None, unit=None):
        """
        Create a new Quantity from a value and a unit.
Justin MacCallum's avatar
Justin MacCallum committed
111

112
113
114
115
116
        Parameters
         - value: (any type, usually a number) Measure of this quantity
         - unit: (Unit) the physical unit, e.g. simtk.unit.meters.
        """
        # When no unit is specified, bend over backwards to handle all one-argument possibilities
117
        if unit is None: # one argument version, copied from UList
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
            if is_unit(value):
                # Unit argument creates an empty list with that unit attached
                unit = value
                value = []
            elif is_quantity(value):
                # Ulist of a Quantity is just the Quantity itself
                unit = value.unit
                value = value._value
            elif _is_string(value):
                unit = dimensionless
            else:
                # Is value a container?
                is_container = True
                try:
                    i = iter(value)
                except TypeError:
                    is_container = False
                if is_container:
                    if len(value) < 1:
                        unit = dimensionless
                    else:
                        first_item = iter(value).next()
                        # Avoid infinite recursion for string, because a one-character
                        # string is its own first element
142
143
144
145
146
147
148
149
150
                        try:
                            isstr = bool(value == first_item)
                        except ValueError:
                            # For numpy, value == first_item returns a numpy
                            # array of booleans, which cannot be evaluated for
                            # truthiness (a ValueError is raised). So in this
                            # case, we don't have a string
                            isstr = False
                        if isstr:
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
                            unit = dimensionless
                        else:
                            unit = Quantity(first_item).unit
                     # Notice that tuples, lists, and numpy.arrays can all be initialized with a list
                    new_container = Quantity([], unit)
                    for item in value:
                        new_container.append(Quantity(item)) # Strips off units into list new_container._value
                    # __class__ trick does not work for numpy.arrays
                    try:
                        import numpy
                        if isinstance(value, numpy.ndarray):
                            value = numpy.array(new_container._value)
                        else:
                            # delegate contruction to container class from list
                            value = value.__class__(new_container._value)
                    except ImportError:
                        # delegate contruction to container class from list
                        value = value.__class__(new_container._value)
                else:
                    # Non-Quantity, non container
                    # Wrap in a dimensionless Quantity
Justin MacCallum's avatar
Justin MacCallum committed
172
                    unit = dimensionless
173
174
175
176
177
        # Accept simple scalar quantities as units
        if is_quantity(unit):
            value = value * unit._value
            unit = unit.unit
        # Use empty list for unspecified values
178
        if value is None:
Justin MacCallum's avatar
Justin MacCallum committed
179
180
            value = []

181
182
        self._value = value
        self.unit = unit
Justin MacCallum's avatar
Justin MacCallum committed
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
    def __getstate__(self):
        state = dict()
        state['_value'] = self._value
        state['unit'] = self.unit
        return state

    def __setstate__(self, state):
        self._value = state['_value']
        self.unit = state['unit']
        return

    def __copy__(self):
        """
        Shallow copy produces a new Quantity with the shallow copy of value and the same unit.
        Because we want copy operations to work just the same way they would on the underlying value.
        """
        return Quantity(copy.copy(self._value), self.unit)

    def __deepcopy__(self, memo):
        """
        Deep copy produces a new Quantity with a deep copy of the value, and the same unit.
        Because we want copy operations to work just the same way they would on the underlying value.
        """
        return Quantity(copy.deepcopy(self._value, memo), self.unit)

    def __getattr__(self, attribute):
        """
        Delegate unrecognized attribute calls to the underlying value type.
        """
        ret_val = getattr(self._value, attribute)
        return ret_val
Justin MacCallum's avatar
Justin MacCallum committed
215

216
217
    def __str__(self):
        """Printable string version of this Quantity.
Justin MacCallum's avatar
Justin MacCallum committed
218

219
220
221
        Returns a string consisting of quantity number followed by unit abbreviation.
        """
        return str(self._value) + ' ' + str(self.unit.get_symbol())
Justin MacCallum's avatar
Justin MacCallum committed
222

223
224
225
226
227
    def __repr__(self):
        """
        """
        return (Quantity.__name__ + '(value=' + repr(self._value) + ', unit=' +
                str(self.unit) + ')')
Justin MacCallum's avatar
Justin MacCallum committed
228

229
230
231
232
233
    def format(self, format_spec):
        return format_spec % self._value + ' ' + str(self.unit.get_symbol())

    def __add__(self, other):
        """Add two Quantities.
Justin MacCallum's avatar
Justin MacCallum committed
234

235
236
        Only Quantities with the same dimensions (e.g. length)
        can be added.  Raises TypeError otherwise.
Justin MacCallum's avatar
Justin MacCallum committed
237

238
239
240
        Parameters
         - self: left hand member of sum
         - other: right hand member of sum
Justin MacCallum's avatar
Justin MacCallum committed
241

242
243
244
245
246
247
248
249
250
251
252
        Returns a new Quantity that is the sum of the two arguments.
        """
        # can only add using like units
        if not self.unit.is_compatible(other.unit):
            raise TypeError('Cannot add two quantities with incompatible units "%s" and "%s".' % (self.unit, other.unit))
        value = self._value + other.value_in_unit(self.unit)
        unit = self.unit
        return Quantity(value, unit)

    def __sub__(self, other):
        """Subtract two Quantities.
Justin MacCallum's avatar
Justin MacCallum committed
253

254
255
        Only Quantities with the same dimensions (e.g. length)
        can be subtracted.  Raises TypeError otherwise.
Justin MacCallum's avatar
Justin MacCallum committed
256

257
258
259
        Parameters
         - self: left hand member (a) of a - b.
         - other: right hand member (b) of a - b.
Justin MacCallum's avatar
Justin MacCallum committed
260

261
262
263
264
265
266
267
        Returns a new Quantity that is the difference of the two arguments.
        """
        if not self.unit.is_compatible(other.unit):
            raise TypeError('Cannot subtract two quantities with incompatible units "%s" and "%s".' % (self.unit, other.unit))
        value = self._value - other.value_in_unit(self.unit)
        unit = self.unit
        return Quantity(value, unit)
Justin MacCallum's avatar
Justin MacCallum committed
268

269
270
271
272
273
    def __eq__(self, other):
        """
        """
        if not is_quantity(other):
            return False
Peter Eastman's avatar
Peter Eastman committed
274
275
276
        if not self.unit.is_compatible(other.unit):
            return False
        return self.value_in_unit(other.unit) == other._value
Justin MacCallum's avatar
Justin MacCallum committed
277

278
279
280
    def __ne__(self, other):
        """
        """
Peter Eastman's avatar
Peter Eastman committed
281
        return not self.__eq__(other)
282

Peter Eastman's avatar
Peter Eastman committed
283
    def __lt__(self, other):
284
        """Compares two quantities.
Justin MacCallum's avatar
Justin MacCallum committed
285

286
        Raises TypeError if the Quantities are of different dimension (e.g. length vs. mass)
Justin MacCallum's avatar
Justin MacCallum committed
287

Peter Eastman's avatar
Peter Eastman committed
288
        Returns True if self < other, False otherwise.
289
        """
Peter Eastman's avatar
Peter Eastman committed
290
        return self._value < other.value_in_unit(self.unit)
291
292

    def __ge__(self, other):
Justin MacCallum's avatar
Justin MacCallum committed
293
        return self._value >= (other.value_in_unit(self.unit))
294
    def __gt__(self, other):
Justin MacCallum's avatar
Justin MacCallum committed
295
        return self._value > (other.value_in_unit(self.unit))
296
    def __le__(self, other):
Justin MacCallum's avatar
Justin MacCallum committed
297
        return self._value <= (other.value_in_unit(self.unit))
298
    def __lt__(self, other):
Justin MacCallum's avatar
Justin MacCallum committed
299
        return self._value < (other.value_in_unit(self.unit))
300

301
302
    _reduce_cache = {}

303
304
305
306
    def reduce_unit(self, guide_unit=None):
        """
        Combine similar component units and scale, to form an
        equal Quantity in simpler units.
Justin MacCallum's avatar
Justin MacCallum committed
307

308
309
        Returns underlying value type if unit is dimensionless.
        """
310
311
312
313
314
315
316
        key = (self.unit, guide_unit)
        if key in Quantity._reduce_cache:
            (unit, value_factor) = Quantity._reduce_cache[key]
        else:
            value_factor = 1.0
            canonical_units = {} # dict of dimensionTuple: (Base/ScaledUnit, exponent)
            # Bias result toward guide units
317
            if guide_unit is not None:
318
319
320
321
322
                for u, exponent in guide_unit.iter_base_or_scaled_units():
                    d = u.get_dimension_tuple()
                    if d not in canonical_units:
                        canonical_units[d] = [u, 0]
            for u, exponent in self.unit.iter_base_or_scaled_units():
323
                d = u.get_dimension_tuple()
324
                # Take first unit found in a dimension as canonical
325
                if d not in canonical_units:
326
327
328
329
330
331
332
333
334
335
336
337
338
                    canonical_units[d] = [u, exponent]
                else:
                    value_factor *= (u.conversion_factor_to(canonical_units[d][0])**exponent)
                    canonical_units[d][1] += exponent
            new_base_units = {}
            for d in canonical_units:
                u, exponent = canonical_units[d]
                if exponent != 0:
                    assert u not in new_base_units
                    new_base_units[u] = exponent
            # Create new unit
            if len(new_base_units) == 0:
                unit = dimensionless
339
            else:
340
341
342
343
344
345
346
347
348
349
                unit = Unit(new_base_units)
            # There might be a factor due to unit conversion, even though unit is dimensionless
            # e.g. suppose unit is meter/centimeter
            if unit.is_dimensionless():
                unit_factor = unit.conversion_factor_to(dimensionless)
                if unit_factor != 1.0:
                    value_factor *= unit_factor
                    # print "value_factor = %s" % value_factor
                unit = dimensionless
            Quantity._reduce_cache[key] = (unit, value_factor)
350
351
352
353
354
355
356
357
358
        # Create Quantity, then scale (in case value is a container)
        # That's why we don't just scale the value.
        result = Quantity(self._value, unit)
        if value_factor != 1.0:
            # __mul__ strips off dimensionless, if appropriate
            result = result * value_factor
        if unit.is_dimensionless():
            assert unit is dimensionless # should have been set earlier in this method
            if is_quantity(result):
359
                result = copy.deepcopy(result._value)
360
361
362
363
        return result

    def __mul__(self, other):
        """Multiply a quantity by another object
Justin MacCallum's avatar
Justin MacCallum committed
364

365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
        Returns a new Quantity that is the product of the self * other,
        unless the resulting unit is dimensionless, in which case the
        underlying value type is returned, instead of a Quantity.
        """
        if is_unit(other):
            # print "quantity * unit"
            # Many other mul/div operations delegate to here because I was debugging
            # a dimensionless unit conversion problem, which I ended up fixing within
            # the reduce_unit() method.
            unit = self.unit * other
            return Quantity(self._value, unit).reduce_unit(self.unit)
        elif is_quantity(other):
            # print "quantity * quantity"
            # Situations where the units cancel can result in scale factors from the unit cancellation.
            # To simplify things, delegate Quantity * Quantity to (Quantity * scalar) * unit
            return (self * other._value) * other.unit
        else:
            # print "quantity * scalar"
            return self._change_units_with_factor(self.unit, other, post_multiply=False)
Justin MacCallum's avatar
Justin MacCallum committed
384

385
386
387
    # value type might not be commutative for multiplication
    def __rmul__(self, other):
        """Multiply a scalar by a Quantity
Justin MacCallum's avatar
Justin MacCallum committed
388
389

        Returns a new Quantity with the same units as self, but with the value
390
391
392
393
394
395
396
397
398
399
400
401
        multiplied by other.
        """
        if is_unit(other):
            raise NotImplementedError('programmer is surprised __rmul__ was called instead of __mul__')
            # print "R unit * quantity"
        elif is_quantity(other):
            # print "R quantity * quantity"
            raise NotImplementedError('programmer is surprised __rmul__ was called instead of __mul__')
        else:
            # print "scalar * quantity"
            return self._change_units_with_factor(self.unit, other, post_multiply=True)
            # return Quantity(other * self._value, self.unit)
Justin MacCallum's avatar
Justin MacCallum committed
402

Peter Eastman's avatar
Peter Eastman committed
403
    def __truediv__(self, other):
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
        """Divide a Quantity by another object

        Returns a new Quantity, unless the resulting unit type is dimensionless,
        in which case the underlying value type is returned.
        """
        if is_unit(other):
            # print "quantity / unit"
            return self * pow(other, -1.0)
            # unit = self.unit / other
            # return Quantity(self._value, unit).reduce_unit(self.unit)
        elif is_quantity(other):
            # print "quantity / quantity"
            # Delegate quantity/quantity to (quantity/scalar)/unit
            return (self/other._value) / other.unit
        else:
            # print "quantity / scalar"
            return self * pow(other, -1.0)
            # return Quantity(self._value / other, self.unit)

423
424
    __div__ = __truediv__

Peter Eastman's avatar
Peter Eastman committed
425
    def __rtruediv__(self, other):
426
        """Divide a scalar by a quantity.
Justin MacCallum's avatar
Justin MacCallum committed
427

428
429
430
431
        Returns a new Quantity.  The resulting units are the inverse of the self argument units.
        """
        if is_unit(other):
            # print "R unit / quantity"
Peter Eastman's avatar
Peter Eastman committed
432
            raise NotImplementedError('programmer is surprised __rtruediv__ was called instead of __truediv__')
433
        elif is_quantity(other):
Peter Eastman's avatar
Peter Eastman committed
434
            raise NotImplementedError('programmer is surprised __rtruediv__ was called instead of __truediv__')
435
436
437
438
439
        else:
            # print "R scalar / quantity"
            return other * pow(self, -1.0)
            # return Quantity(other / self._value, pow(self.unit, -1.0))

440
441
    __rdiv__ = __rtruediv__

442
443
    def __pow__(self, exponent):
        """Raise a Quantity to a power.
Justin MacCallum's avatar
Justin MacCallum committed
444

445
446
447
448
449
        Generally both the value and the unit of the Quantity are affected by this operation.

        Returns a new Quantity equal to self**exponent.
        """
        return Quantity(pow(self._value, exponent), pow(self.unit, exponent))
Justin MacCallum's avatar
Justin MacCallum committed
450

451
452
453
    def sqrt(self):
        """
        Returns square root of a Quantity.
Justin MacCallum's avatar
Justin MacCallum committed
454

455
456
457
458
459
460
461
462
463
464
465
        Raises ArithmeticError if component exponents are not even.
        This behavior can be changed if you present a reasonable real life case to me.
        """
        # There might be a conversion factor from taking the square root of the unit
        new_value = math.sqrt(self._value)
        new_unit = self.unit.sqrt()
        unit_factor = self.unit.conversion_factor_to(new_unit*new_unit)
        if unit_factor != 1.0:
            new_value *= math.sqrt(unit_factor)
        return Quantity(value=new_value, unit=new_unit)

466
    def sum(self, *args, **kwargs):
467
468
469
470
471
472
        """
        Computes the sum of a sequence, with the result having the same unit as
        the current sequence.

        If the value is not iterable, it raises a TypeError (same behavior as if
        you tried to iterate over, for instance, an integer).
473
474
475
476

        This function can take as arguments any arguments recognized by
        `numpy.sum`. If arguments are passed to a non-numpy array, a TypeError
        is raised
477
478
479
        """
        try:
            # This will be much faster for numpy arrays
480
            mysum = self._value.sum(*args, **kwargs)
481
        except AttributeError:
482
483
            if args or kwargs:
                raise TypeError('Unsupported arguments for Quantity.sum')
Peter Eastman's avatar
Peter Eastman committed
484
485
486
487
488
489
            if len(self._value) == 0:
                mysum = 0
            else:
                mysum = self._value[0]
                for i in range(1, len(self._value)):
                    mysum += self._value[i]
490
491
        return Quantity(mysum, self.unit)

492
    def mean(self, *args, **kwargs):
493
494
495
496
497
        """
        Computes the mean of a sequence, with the result having the same unit as
        the current sequence.

        If the value is not iterable, it raises a TypeError
498
499
500
501

        This function can take as arguments any arguments recognized by
        `numpy.mean`. If arguments are passed to a non-numpy array, a TypeError
        is raised
502
503
504
        """
        try:
            # Faster for numpy arrays
505
            mean = self._value.mean(*args, **kwargs)
506
        except AttributeError:
507
508
            if args or kwargs:
                raise TypeError('Unsupported arguments for Quantity.mean')
509
510
511
            mean = self.sum() / len(self._value)
        return Quantity(mean, self.unit)

512
    def std(self, *args, **kwargs):
513
514
515
516
517
        """
        Computes the square root of the variance of a sequence, with the result
        having the same unit as the current sequence.

        If the value is not iterable, it raises a TypeError
518
519
520
521

        This function can take as arguments any arguments recognized by
        `numpy.std`. If arguments are passed to a non-numpy array, a TypeError
        is raised
522
523
524
        """
        try:
            # Faster for numpy arrays
525
            std = self._value.std(*args, **kwargs)
526
        except AttributeError:
527
528
            if args or kwargs:
                raise TypeError('Unsupported arguments for Quantity.std')
Jason Swails's avatar
Jason Swails committed
529
            mean = self.mean()
530
            for val in self._value:
Jason Swails's avatar
Jason Swails committed
531
532
533
534
                res = mean - val
                var += res * res
            var /= len(self._value)
            std = math.sqrt(var)
535
536
        return Quantity(std, self.unit)

537
    def max(self, *args, **kwargs):
538
539
540
541
542
        """
        Computes the maximum value of the sequence, with the result having the
        same unit as the current sequence.

        If the value is not iterable, it raises a TypeError
543
544
545
546

        This function can take as arguments any arguments recognized by
        `numpy.max`. If arguments are passed to a non-numpy array, a TypeError
        is raised
547
548
549
        """
        try:
            # Faster for numpy arrays
550
            mymax = self._value.max(*args, **kwargs)
551
        except AttributeError:
552
553
            if args or kwargs:
                raise TypeError('Unsupported arguments for Quantity.max')
554
555
556
            mymax = max(self._value)
        return Quantity(mymax, self.unit)

557
    def min(self, *args, **kwargs):
558
559
560
561
562
        """
        Computes the minimum value of the sequence, with the result having the
        same unit as the current sequence.

        If the value is not iterable, it raises a TypeError
563
564
565
566

        This function can take as arguments any arguments recognized by
        `numpy.min`. If arguments are passed to a non-numpy array, a TypeError
        is raised
567
568
569
        """
        try:
            # Faster for numpy arrays
570
            mymin = self._value.min(*args, **kwargs)
571
        except AttributeError:
572
573
            if args or kwargs:
                raise TypeError('Unsupported arguments for Quantity.min')
574
575
576
            mymin = min(self._value)
        return Quantity(mymin, self.unit)

577
578
579
580
581
582
    def reshape(self, shape, order='C'):
        """
        Same as numpy.ndarray.reshape, except the result is a Quantity with the
        same units as the current object rather than a plain numpy.ndarray
        """
        try:
583
            return Quantity(self._value.reshape(shape, order=order), self.unit)
584
585
586
587
        except AttributeError:
            raise AttributeError('Only numpy array Quantity objects can be '
                                 'reshaped')

588
589
590
    def __abs__(self):
        """
        Return absolute value of a Quantity.
Justin MacCallum's avatar
Justin MacCallum committed
591

592
593
594
595
        The unit is unchanged.  A negative value of self will result in a positive value
        in the result.
        """
        return Quantity(abs(self._value), self.unit)
Justin MacCallum's avatar
Justin MacCallum committed
596

597
598
599
600
601
    def __pos__(self):
        """
        Returns a reference to self.
        """
        return Quantity(+(self._value), self.unit)
Justin MacCallum's avatar
Justin MacCallum committed
602

603
604
    def __neg__(self):
        """Negate a Quantity.
Justin MacCallum's avatar
Justin MacCallum committed
605

606
607
608
        Returns a new Quantity with a different sign on the value.
        """
        return Quantity(-(self._value), self.unit)
Justin MacCallum's avatar
Justin MacCallum committed
609

610
611
612
613
614
615
616
617
618
619
620
621
622
    def __nonzero__(self):
        """Returns True if value underlying Quantity is zero, False otherwise.
        """
        return bool(self._value)

    def __complex__(self):
        return Quantity(complex(self._value), self.unit)
    def __float__(self):
        return Quantity(float(self._value), self.unit)
    def __int__(self):
        return Quantity(int(self._value), self.unit)
    def __long__(self):
        return Quantity(int(self._value), self.unit)
Justin MacCallum's avatar
Justin MacCallum committed
623

624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
    def value_in_unit(self, unit):
        """
        Returns underlying value, in the specified units.
        """
        val = self.in_units_of(unit)
        if is_quantity(val):
            return val._value
        else: # naked dimensionless
            return val

    def value_in_unit_system(self, system):
        """
        Returns the underlying value type, after conversion to a particular unit system.
        """
        result = self.in_unit_system(system)
        if is_quantity(result):
            return result._value
        else:
            return result # dimensionless
Justin MacCallum's avatar
Justin MacCallum committed
643

644
645
646
647
    def in_unit_system(self, system):
        """
        Returns a new Quantity equal to this one, expressed in a particular unit system.
        """
648
649
650
        new_units = system.express_unit(self.unit)
        f = self.unit.conversion_factor_to(new_units)
        return self._change_units_with_factor(new_units, f)
651
652
653
654
655
656
657

    def in_units_of(self, other_unit):
        """
        Returns an equal Quantity expressed in different units.

        If the units are the same as those in self, a reference to self is returned.
        Raises a TypeError if the new unit is not compatible with the original unit.
Justin MacCallum's avatar
Justin MacCallum committed
658

659
660
661
662
663
664
665
666
        The post_multiply argument is used in case the multiplication operation is not commutative.
          i.e. result = factor * value when post_multiply is False
          and  result = value * factor when post_multiply is True
        """
        if not self.unit.is_compatible(other_unit):
            raise TypeError('Unit "%s" is not compatible with Unit "%s".' % (self.unit, other_unit))
        f = self.unit.conversion_factor_to(other_unit)
        return self._change_units_with_factor(other_unit, f)
Justin MacCallum's avatar
Justin MacCallum committed
667

668
669
670
671
672
673
674
675
676
677
    def _change_units_with_factor(self, new_unit, factor, post_multiply=True):
        # numpy arrays cannot be compared with 1.0, so just "try"
        factor_is_identity = False
        try:
            if (factor == 1.0):
                factor_is_identity = True
        except ValueError:
            pass
        if factor_is_identity:
            # No multiplication required
678
            result = Quantity(copy.deepcopy(self._value), new_unit)
679
680
681
682
683
684
        else:
            try:
                # multiply operator, if it exists, is preferred
                if post_multiply:
                    value = self._value * factor # works for number, numpy.array, or vec3, e.g.
                else:
Justin MacCallum's avatar
Justin MacCallum committed
685
                    value = factor * self._value # works for number, numpy.array, or vec3, e.g.
686
687
                result = Quantity(value, new_unit)
            except TypeError:
688
                value = copy.deepcopy(self._value)
689
                result = Quantity(self._scale_sequence(value, factor, post_multiply), new_unit)
690
691
692
693
694
        if (new_unit.is_dimensionless()):
            return result._value
        else:
            return result

695
696
697
    def _scale_sequence(self, value, factor, post_multiply):
        try:
            if post_multiply:
698
                value = value*factor
699
            else:
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
                value = factor*value
        except TypeError:
            try:
                if post_multiply:
                    if isinstance(value, tuple):
                        value = tuple([x*factor for x in value])
                    else:
                        for i in range(len(value)):
                            value[i] = value[i]*factor
                else:
                    if isinstance(value, tuple):
                        value = tuple([factor*x for x in value])
                    else:
                        for i in range(len(value)):
                            value[i] = factor*value[i]
            except TypeError as ex:
                if isinstance(value, tuple):
                    value = tuple([self._scale_sequence(x, factor, post_multiply) for x in value])
718
719
                else:
                    for i in range(len(value)):
720
                        value[i] = self._scale_sequence(value[i], factor, post_multiply)
721
722
723
        return value


724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741

    ####################################
    ### Sequence methods of Quantity ###
    ###  in case value is a sequence ###
    ####################################

    def __len__(self):
        """
        Return size of internal value type.
        """
        return len(self._value)

    def __getitem__(self, key):
        """
        Keep the same units on contained elements.
        """
        assert not is_quantity(self._value[key])
        return Quantity(self._value[key], self.unit)
Justin MacCallum's avatar
Justin MacCallum committed
742

743
744
745
746
747
748
749
750
751
752
753
754
755
756
    def __setitem__(self, key, value):
        # Delegate slices to one-at-a time ___setitem___
        if isinstance(key, slice): # slice
            indices = key.indices(len(self))
            for i in range(*indices):
                self[i] = value[i]
        else: # single index
            # Check unit compatibility
            if self.unit.is_dimensionless() and is_dimensionless(value):
                pass # OK
            elif not self.unit.is_compatible(value.unit):
                raise TypeError('Unit "%s" is not compatible with Unit "%s".' % (self.unit, value.unit))
            self._value[key] = value / self.unit
            assert not is_quantity(self._value[key])
Justin MacCallum's avatar
Justin MacCallum committed
757

758
759
760
761
762
    def __delitem__(self, key):
        del(self._value[key])

    def __contains__(self, item):
        return self._value.__contains__(item.value_in_unit(self.unit))
Justin MacCallum's avatar
Justin MacCallum committed
763

764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
    def __iter__(self):
        for item in self._value:
            yield Quantity(item, self.unit)

    def count(self, item):
        return self._value.count(item.value_in_unit(self.unit))
    def index(self, item):
        return self._value.index(item.value_in_unit(self.unit))
    def append(self, item):
        if is_quantity(item):
            return self._value.append(item.value_in_unit(self.unit))
        elif is_dimensionless(self.unit):
            return self._value.append(item)
        else:
            raise TypeError("Cannot append item without units into list with units")
    def extend(self, rhs):
        self._value.extend(rhs.value_in_unit(self.unit))
    def insert(self, index, item):
        self._value.insert(index, item.value_in_unit(self.unit))
    def remove(self, item):
        self._value.remove(item)
    def pop(self, *args):
        return self._value.pop(*args) * self.unit
    # list.reverse will automatically delegate correctly
    # list.sort with no arguments will delegate correctly
    # list.sort with a comparison function cannot be done correctly


def is_quantity(x):
    """
    Returns True if x is a Quantity, False otherwise.
    """
    return isinstance(x, Quantity)

def is_dimensionless(x):
    """
    """
    if is_unit(x):
        return x.is_dimensionless()
    elif is_quantity(x):
        return x.unit.is_dimensionless()
    else:
        # everything else in the universe is dimensionless
        return True

# Strings can cause trouble
# as can any container that has infinite levels of containment
def _is_string(x):
     # step 1) String is always a container
     # and its contents are themselves containers.
     if isinstance(x, str):
         return True
     try:
         first_item = iter(x).next()
         inner_item = iter(first_item).next()
         if first_item is inner_item:
             return True
         else:
             return False
     except TypeError:
         return False
     except StopIteration:
         return False


# run module directly for testing
if __name__=='__main__':
    # Test the examples in the docstrings
    import doctest, sys
    doctest.testmod(sys.modules[__name__])