"""
YTArray class.



"""
from __future__ import print_function
#-----------------------------------------------------------------------------
# Copyright (c) 2013, yt Development Team.
#
# Distributed under the terms of the Modified BSD License.
#
# The full license is in the file COPYING.txt, distributed with this software.
#-----------------------------------------------------------------------------

import copy
import numpy as np

from distutils.version import LooseVersion
from functools import wraps
from numpy import \
    add, subtract, multiply, divide, logaddexp, logaddexp2, true_divide, \
    floor_divide, negative, power, remainder, mod, absolute, rint, \
    sign, conj, exp, exp2, log, log2, log10, expm1, log1p, sqrt, square, \
    reciprocal, sin, cos, tan, arcsin, arccos, arctan, arctan2, \
    hypot, sinh, cosh, tanh, arcsinh, arccosh, arctanh, deg2rad, rad2deg, \
    bitwise_and, bitwise_or, bitwise_xor, invert, left_shift, right_shift, \
    greater, greater_equal, less, less_equal, not_equal, equal, logical_and, \
    logical_or, logical_xor, logical_not, maximum, minimum, fmax, fmin, \
    isreal, iscomplex, isfinite, isinf, isnan, signbit, copysign, nextafter, \
    modf, ldexp, frexp, fmod, floor, ceil, trunc, fabs, spacing

try:
    # numpy 1.13 or newer
    from numpy import positive, divmod as divmod_, isnat, heaviside
except ImportError:
    positive, divmod_, isnat, heaviside = (None,)*4

from yt.units.unit_object import Unit, UnitParseError
from yt.units.unit_registry import UnitRegistry
from yt.units.dimensions import \
    angle, \
    current_mks, \
    dimensionless, \
    em_dimensions
from yt.utilities.exceptions import \
    YTUnitOperationError, YTUnitConversionError, \
    YTUfuncUnitError, YTIterableUnitCoercionError, \
    YTInvalidUnitEquivalence, YTEquivalentDimsError
from yt.utilities.lru_cache import lru_cache
from numbers import Number as numeric_type
from yt.utilities.on_demand_imports import _astropy
from sympy import Rational
from yt.units.unit_lookup_table import \
    default_unit_symbol_lut
from yt.units.equivalencies import equivalence_registry
from yt.utilities.logger import ytLogger as mylog
from .pint_conversions import convert_pint_units

NULL_UNIT = Unit()
POWER_SIGN_MAPPING = {multiply: 1, divide: -1}

# redefine this here to avoid a circular import from yt.funcs
def iterable(obj):
    try: len(obj)
    except: return False
    return True

def return_arr(func):
    @wraps(func)
    def wrapped(*args, **kwargs):
        ret, units = func(*args, **kwargs)
        if ret.shape == ():
            return YTQuantity(ret, units)
        else:
            # This could be a subclass, so don't call YTArray directly.
            return type(args[0])(ret, units)
    return wrapped

@lru_cache(maxsize=128, typed=False)
def sqrt_unit(unit):
    return unit**0.5

@lru_cache(maxsize=128, typed=False)
def multiply_units(unit1, unit2):
    return unit1 * unit2

def preserve_units(unit1, unit2=None):
    return unit1

@lru_cache(maxsize=128, typed=False)
def power_unit(unit, power):
    return unit**power

@lru_cache(maxsize=128, typed=False)
def square_unit(unit):
    return unit*unit

@lru_cache(maxsize=128, typed=False)
def divide_units(unit1, unit2):
    return unit1/unit2

@lru_cache(maxsize=128, typed=False)
def reciprocal_unit(unit):
    return unit**-1

def passthrough_unit(unit, unit2=None):
    return unit

def return_without_unit(unit, unit2=None):
    return None

def arctan2_unit(unit1, unit2):
    return NULL_UNIT

def comparison_unit(unit1, unit2=None):
    return None

def invert_units(unit):
    raise TypeError(
        "Bit-twiddling operators are not defined for YTArray instances")

def bitop_units(unit1, unit2):
    raise TypeError(
        "Bit-twiddling operators are not defined for YTArray instances")

def get_inp_u_unary(ufunc, inputs, out_arr=None):
    inp = inputs[0]
    u = getattr(inp, 'units', None)
    if u is None:
        u = NULL_UNIT
    if u.dimensions is angle and ufunc in trigonometric_operators:
        inp = inp.in_units('radian').v
        if out_arr is not None:
            out_arr = ufunc(inp).view(np.ndarray)
    return out_arr, inp, u

def get_inp_u_binary(ufunc, inputs):
    inp1 = coerce_iterable_units(inputs[0])
    inp2 = coerce_iterable_units(inputs[1])
    unit1 = getattr(inp1, 'units', None)
    unit2 = getattr(inp2, 'units', None)
    ret_class = get_binary_op_return_class(type(inp1), type(inp2))
    if unit1 is None:
        unit1 = Unit(registry=getattr(unit2, 'registry', None))
    if unit2 is None and ufunc is not power:
        unit2 = Unit(registry=getattr(unit1, 'registry', None))
    elif ufunc is power:
        unit2 = inp2
        if isinstance(unit2, np.ndarray):
            if isinstance(unit2, YTArray):
                if unit2.units.is_dimensionless:
                    pass
                else:
                    raise YTUnitOperationError(ufunc, unit1, unit2)
            unit2 = 1.0
    return (inp1, inp2), (unit1, unit2), ret_class

def handle_preserve_units(inps, units, ufunc, ret_class):
    if units[0] != units[1]:
        any_nonzero = [np.any(inps[0]), np.any(inps[1])]
        if any_nonzero[0] == np.bool_(False):
            units = (units[1], units[1])
        elif any_nonzero[1] == np.bool_(False):
            units = (units[0], units[0])
        else:
            if not units[0].same_dimensions_as(units[1]):
                raise YTUnitOperationError(ufunc, *units)
            inps = (inps[0], ret_class(inps[1]).to(
                ret_class(inps[0]).units))
    return inps, units

def handle_comparison_units(inps, units, ufunc, ret_class, raise_error=False):
    if units[0] != units[1]:
        u1d = units[0].is_dimensionless
        u2d = units[1].is_dimensionless
        any_nonzero = [np.any(inps[0]), np.any(inps[1])]
        if any_nonzero[0] == np.bool_(False):
            units = (units[1], units[1])
        elif any_nonzero[1] == np.bool_(False):
            units = (units[0], units[0])
        elif not any([u1d, u2d]):
            if not units[0].same_dimensions_as(units[1]):
                raise YTUnitOperationError(ufunc, *units)
            else:
                if raise_error:
                    raise YTUfuncUnitError(ufunc, *units)
                inps = (inps[0], ret_class(inps[1]).to(
                    ret_class(inps[0]).units))
    return inps, units

def handle_multiply_divide_units(unit, units, out, out_arr):
    if unit.is_dimensionless and unit.base_value != 1.0:
        if not units[0].is_dimensionless:
            if units[0].dimensions == units[1].dimensions:
                out_arr = np.multiply(out_arr.view(np.ndarray),
                                      unit.base_value, out=out)
                unit = Unit(registry=unit.registry)
    return out, out_arr, unit

def coerce_iterable_units(input_object):
    if isinstance(input_object, np.ndarray):
        return input_object
    if iterable(input_object):
        if any([isinstance(o, YTArray) for o in input_object]):
            ff = getattr(input_object[0], 'units', NULL_UNIT, )
            if any([ff != getattr(_, 'units', NULL_UNIT) for _ in input_object]):
                raise YTIterableUnitCoercionError(input_object)
            # This will create a copy of the data in the iterable.
            return YTArray(input_object)
        return input_object
    else:
        return input_object

def sanitize_units_mul(this_object, other_object):
    inp = coerce_iterable_units(this_object)
    ret = coerce_iterable_units(other_object)
    # If the other object is a YTArray and has the same dimensions as the object
    # under consideration, convert so we don't mix units with the same
    # dimensions.
    if isinstance(ret, YTArray):
        if inp.units.same_dimensions_as(ret.units):
            ret.in_units(inp.units)
    return ret

def sanitize_units_add(this_object, other_object, op_string):
    inp = coerce_iterable_units(this_object)
    ret = coerce_iterable_units(other_object)
    # Make sure the other object is a YTArray before we use the `units`
    # attribute.
    if isinstance(ret, YTArray):
        if not inp.units.same_dimensions_as(ret.units):
            # handle special case of adding or subtracting with zero or
            # array filled with zero
            if not np.any(other_object):
                return ret.view(np.ndarray)
            elif not np.any(this_object):
                return ret
            raise YTUnitOperationError(op_string, inp.units, ret.units)
        ret = ret.in_units(inp.units)
    else:
        # If the other object is not a YTArray, then one of the arrays must be
        # dimensionless or filled with zeros
        if not inp.units.is_dimensionless and np.any(ret):
            raise YTUnitOperationError(op_string, inp.units, dimensionless)
    return ret

def validate_comparison_units(this, other, op_string):
    # Check that other is a YTArray.
    if hasattr(other, 'units'):
        if this.units.expr is other.units.expr:
            if this.units.base_value == other.units.base_value:
                return other
        if not this.units.same_dimensions_as(other.units):
            raise YTUnitOperationError(op_string, this.units, other.units)
        return other.in_units(this.units)

    return other

@lru_cache(maxsize=128, typed=False)
def _unit_repr_check_same(my_units, other_units):
    """
    Takes a Unit object, or string of known unit symbol, and check that it
    is compatible with this quantity. Returns Unit object.

    """
    # let Unit() handle units arg if it's not already a Unit obj.
    if not isinstance(other_units, Unit):
        other_units = Unit(other_units, registry=my_units.registry)

    equiv_dims = em_dimensions.get(my_units.dimensions, None)
    if equiv_dims == other_units.dimensions:
        if current_mks in equiv_dims.free_symbols:
            base = "SI"
        else:
            base = "CGS"
        raise YTEquivalentDimsError(my_units, other_units, base)

    if not my_units.same_dimensions_as(other_units):
        raise YTUnitConversionError(
            my_units, my_units.dimensions, other_units, other_units.dimensions)

    return other_units

unary_operators = (
    negative, absolute, rint, sign, conj, exp, exp2, log, log2,
    log10, expm1, log1p, sqrt, square, reciprocal, sin, cos, tan, arcsin,
    arccos, arctan, sinh, cosh, tanh, arcsinh, arccosh, arctanh, deg2rad,
    rad2deg, invert, logical_not, isreal, iscomplex, isfinite, isinf, isnan,
    signbit, floor, ceil, trunc, modf, frexp, fabs, spacing, positive, isnat,
)

binary_operators = (
    add, subtract, multiply, divide, logaddexp, logaddexp2, true_divide, power,
    remainder, mod, arctan2, hypot, bitwise_and, bitwise_or, bitwise_xor,
    left_shift, right_shift, greater, greater_equal, less, less_equal,
    not_equal, equal, logical_and, logical_or, logical_xor, maximum, minimum,
    fmax, fmin, copysign, nextafter, ldexp, fmod, divmod_, heaviside
)

trigonometric_operators = (
    sin, cos, tan,
)

class YTArray(np.ndarray):
    """
    An ndarray subclass that attaches a symbolic unit object to the array data.

    Parameters
    ----------

    input_array : :obj:`!iterable`
        A tuple, list, or array to attach units to
    input_units : String unit specification, unit symbol object, or astropy units
        The units of the array. Powers must be specified using python
        syntax (cm**3, not cm^3).
    registry : ~yt.units.unit_registry.UnitRegistry
        The registry to create units from. If input_units is already associated
        with a unit registry and this is specified, this will be used instead of
        the registry associated with the unit object.
    dtype : data-type
        The dtype of the array data. Defaults to the dtype of the input data,
        or, if none is found, uses np.float64
    bypass_validation : boolean
        If True, all input validation is skipped. Using this option may produce
        corrupted, invalid units or array data, but can lead to significant
        speedups in the input validation logic adds significant overhead. If set,
        input_units *must* be a valid unit object. Defaults to False.

    Examples
    --------

    >>> from yt import YTArray
    >>> a = YTArray([1, 2, 3], 'cm')
    >>> b = YTArray([4, 5, 6], 'm')
    >>> a + b
    YTArray([ 401.,  502.,  603.]) cm
    >>> b + a
    YTArray([ 4.01,  5.02,  6.03]) m

    NumPy ufuncs will pass through units where appropriate.

    >>> import numpy as np
    >>> a = YTArray(np.arange(8) - 4, 'g/cm**3')
    >>> np.abs(a)
    YTArray([4, 3, 2, 1, 0, 1, 2, 3]) g/cm**3

    and strip them when it would be annoying to deal with them.

    >>> np.log10(a)
    array([       -inf,  0.        ,  0.30103   ,  0.47712125,  0.60205999,
            0.69897   ,  0.77815125,  0.84509804])

    YTArray is tightly integrated with yt datasets:

    >>> import yt
    >>> ds = yt.load('IsolatedGalaxy/galaxy0030/galaxy0030')
    >>> a = ds.arr(np.ones(5), 'code_length')
    >>> a.in_cgs()
    YTArray([  3.08600000e+24,   3.08600000e+24,   3.08600000e+24,
             3.08600000e+24,   3.08600000e+24]) cm

    This is equivalent to:

    >>> b = YTArray(np.ones(5), 'code_length', registry=ds.unit_registry)
    >>> np.all(a == b)
    True

    """
    _ufunc_registry = {
        add: preserve_units,
        subtract: preserve_units,
        multiply: multiply_units,
        divide: divide_units,
        logaddexp: return_without_unit,
        logaddexp2: return_without_unit,
        true_divide: divide_units,
        floor_divide: divide_units,
        negative: passthrough_unit,
        power: power_unit,
        remainder: preserve_units,
        mod: preserve_units,
        fmod: preserve_units,
        absolute: passthrough_unit,
        fabs: passthrough_unit,
        rint: return_without_unit,
        sign: return_without_unit,
        conj: passthrough_unit,
        exp: return_without_unit,
        exp2: return_without_unit,
        log: return_without_unit,
        log2: return_without_unit,
        log10: return_without_unit,
        expm1: return_without_unit,
        log1p: return_without_unit,
        sqrt: sqrt_unit,
        square: square_unit,
        reciprocal: reciprocal_unit,
        sin: return_without_unit,
        cos: return_without_unit,
        tan: return_without_unit,
        sinh: return_without_unit,
        cosh: return_without_unit,
        tanh: return_without_unit,
        arcsin: return_without_unit,
        arccos: return_without_unit,
        arctan: return_without_unit,
        arctan2: arctan2_unit,
        arcsinh: return_without_unit,
        arccosh: return_without_unit,
        arctanh: return_without_unit,
        hypot: preserve_units,
        deg2rad: return_without_unit,
        rad2deg: return_without_unit,
        bitwise_and: bitop_units,
        bitwise_or: bitop_units,
        bitwise_xor: bitop_units,
        invert: invert_units,
        left_shift: bitop_units,
        right_shift: bitop_units,
        greater: comparison_unit,
        greater_equal: comparison_unit,
        less: comparison_unit,
        less_equal: comparison_unit,
        not_equal: comparison_unit,
        equal: comparison_unit,
        logical_and: comparison_unit,
        logical_or: comparison_unit,
        logical_xor: comparison_unit,
        logical_not: return_without_unit,
        maximum: preserve_units,
        minimum: preserve_units,
        fmax: preserve_units,
        fmin: preserve_units,
        isreal: return_without_unit,
        iscomplex: return_without_unit,
        isfinite: return_without_unit,
        isinf: return_without_unit,
        isnan: return_without_unit,
        signbit: return_without_unit,
        copysign: passthrough_unit,
        nextafter: preserve_units,
        modf: passthrough_unit,
        ldexp: bitop_units,
        frexp: return_without_unit,
        floor: passthrough_unit,
        ceil: passthrough_unit,
        trunc: passthrough_unit,
        spacing: passthrough_unit,
        positive: passthrough_unit,
        divmod_: passthrough_unit,
        isnat: return_without_unit,
        heaviside: preserve_units,
    }

    __array_priority__ = 2.0

    def __new__(cls, input_array, input_units=None, registry=None, dtype=None,
                bypass_validation=False):
        if dtype is None:
            dtype = getattr(input_array, 'dtype', np.float64)
        if bypass_validation is True:
            obj = np.asarray(input_array, dtype=dtype).view(cls)
            obj.units = input_units
            if registry is not None:
                obj.units.registry = registry
            return obj
        if input_array is NotImplemented:
            return input_array.view(cls)
        if registry is None and isinstance(input_units, (str, bytes)):
            if input_units.startswith('code_'):
                raise UnitParseError(
                    "Code units used without referring to a dataset. \n"
                    "Perhaps you meant to do something like this instead: \n"
                    "ds.arr(%s, \"%s\")" % (input_array, input_units)
                    )
        if isinstance(input_array, YTArray):
            ret = input_array.view(cls)
            if input_units is None:
                if registry is None:
                    ret.units = input_array.units
                else:
                    units = Unit(str(input_array.units), registry=registry)
                    ret.units = units
            elif isinstance(input_units, Unit):
                ret.units = input_units
            else:
                ret.units = Unit(input_units, registry=registry)
            return ret
        elif isinstance(input_array, np.ndarray):
            pass
        elif iterable(input_array) and input_array:
            if isinstance(input_array[0], YTArray):
                return YTArray(np.array(input_array, dtype=dtype),
                               input_array[0].units, registry=registry)

        # Input array is an already formed ndarray instance
        # We first cast to be our class type

        obj = np.asarray(input_array, dtype=dtype).view(cls)

        # Check units type
        if input_units is None:
            # Nothing provided. Make dimensionless...
            units = Unit()
        elif isinstance(input_units, Unit):
            if registry and registry is not input_units.registry:
                units = Unit(str(input_units), registry=registry)
            else:
                units = input_units
        else:
            # units kwarg set, but it's not a Unit object.
            # don't handle all the cases here, let the Unit class handle if
            # it's a str.
            units = Unit(input_units, registry=registry)

        # Attach the units
        obj.units = units

        return obj

    def __repr__(self):
        """

        """
        return super(YTArray, self).__repr__()+' '+self.units.__repr__()

    def __str__(self):
        """

        """
        return str(self.view(np.ndarray)) + ' ' + str(self.units)

    #
    # Start unit conversion methods
    #

    def convert_to_units(self, units):
        """
        Convert the array and units to the given units.

        Parameters
        ----------
        units : Unit object or str
            The units you want to convert to.

        """
        new_units = _unit_repr_check_same(self.units, units)
        (conversion_factor, offset) = self.units.get_conversion_factor(new_units)

        self.units = new_units
        values = self.d
        values *= conversion_factor

        if offset:
            np.subtract(self, offset*self.uq, self)

        return self

    def convert_to_base(self, unit_system="cgs"):
        """
        Convert the array and units to the equivalent base units in
        the specified unit system.

        Parameters
        ----------
        unit_system : string, optional
            The unit system to be used in the conversion. If not specified,
            the default base units of cgs are used.

        Examples
        --------
        >>> E = YTQuantity(2.5, "erg/s")
        >>> E.convert_to_base(unit_system="galactic")
        """
        return self.convert_to_units(self.units.get_base_equivalent(unit_system))

    def convert_to_cgs(self):
        """
        Convert the array and units to the equivalent cgs units.

        """
        return self.convert_to_units(self.units.get_cgs_equivalent())

    def convert_to_mks(self):
        """
        Convert the array and units to the equivalent mks units.

        """
        return self.convert_to_units(self.units.get_mks_equivalent())

    def in_units(self, units, equivalence=None, **kwargs):
        """
        Creates a copy of this array with the data in the supplied 
        units, and returns it.

        Optionally, an equivalence can be specified to convert to an 
        equivalent quantity which is not in the same dimensions. 

        .. note::

            All additional keyword arguments are passed to the 
            equivalency, which should be used if that particular 
            equivalency requires them.

        Parameters
        ----------
        units : Unit object or string
            The units you want to get a new quantity in.
        equivalence : string, optional
            The equivalence you wish to use. To see which 
            equivalencies are supported for this unitful 
            quantity, try the :meth:`list_equivalencies` 
            method. Default: None

        Returns
        -------
        YTArray
        """
        if equivalence is None:
            new_units = _unit_repr_check_same(self.units, units)
            (conversion_factor, offset) = self.units.get_conversion_factor(new_units)

            new_array = type(self)(self.ndview * conversion_factor, new_units)

            if offset:
                np.subtract(new_array, offset*new_array.uq, new_array)

            return new_array
        else:
            return self.to_equivalent(units, equivalence, **kwargs)

    def to(self, units, equivalence=None, **kwargs):
        """
        An alias for YTArray.in_units().

        See the docstrings of that function for details.
        """
        return self.in_units(units, equivalence=equivalence, **kwargs)

    def to_value(self, units=None, equivalence=None, **kwargs):
        """
        Creates a copy of this array with the data in the supplied
        units, and returns it without units. Output is therefore a 
        bare NumPy array.

        Optionally, an equivalence can be specified to convert to an 
        equivalent quantity which is not in the same dimensions. 

        .. note::

            All additional keyword arguments are passed to the 
            equivalency, which should be used if that particular 
            equivalency requires them.

        Parameters
        ----------
        units : Unit object or string, optional
            The units you want to get the bare quantity in. If not
            specified, the value will be returned in the current units.

        equivalence : string, optional
            The equivalence you wish to use. To see which 
            equivalencies are supported for this unitful 
            quantity, try the :meth:`list_equivalencies` 
            method. Default: None

        Returns
        -------
        NumPy array
        """
        if units is None:
            v = self.value
        else:
            v = self.in_units(units, equivalence=equivalence, **kwargs).value
        if isinstance(self, YTQuantity):
            return float(v)
        else:
            return v

    def in_base(self, unit_system="cgs"):
        """
        Creates a copy of this array with the data in the specified unit system,
        and returns it in that system's base units.

        Parameters
        ----------
        unit_system : string, optional
            The unit system to be used in the conversion. If not specified,
            the default base units of cgs are used.

        Examples
        --------
        >>> E = YTQuantity(2.5, "erg/s")
        >>> E_new = E.in_base(unit_system="galactic")
        """
        return self.in_units(self.units.get_base_equivalent(unit_system))

    def in_cgs(self):
        """
        Creates a copy of this array with the data in the equivalent cgs units,
        and returns it.

        Returns
        -------
        Quantity object with data converted to cgs units.

        """
        return self.in_units(self.units.get_cgs_equivalent())

    def in_mks(self):
        """
        Creates a copy of this array with the data in the equivalent mks units,
        and returns it.

        Returns
        -------
        Quantity object with data converted to mks units.

        """
        return self.in_units(self.units.get_mks_equivalent())

    def to_equivalent(self, unit, equiv, **kwargs):
        """
        Convert a YTArray or YTQuantity to an equivalent, e.g., something that is
        related by only a constant factor but not in the same units.

        Parameters
        ----------
        unit : string
            The unit that you wish to convert to.
        equiv : string
            The equivalence you wish to use. To see which equivalencies are
            supported for this unitful quantity, try the
            :meth:`list_equivalencies` method.

        Examples
        --------
        >>> a = yt.YTArray(1.0e7,"K")
        >>> a.to_equivalent("keV", "thermal")
        """
        conv_unit = Unit(unit, registry=self.units.registry)
        if self.units.same_dimensions_as(conv_unit):
            return self.in_units(conv_unit)
        this_equiv = equivalence_registry[equiv]()
        oneway_or_equivalent = (
            conv_unit.has_equivalent(equiv) or this_equiv._one_way)
        if self.has_equivalent(equiv) and oneway_or_equivalent:
            new_arr = this_equiv.convert(
                self, conv_unit.dimensions, **kwargs)
            if isinstance(new_arr, tuple):
                try:
                    return type(self)(new_arr[0], new_arr[1]).in_units(unit)
                except YTUnitConversionError:
                    raise YTInvalidUnitEquivalence(equiv, self.units, unit)
            else:
                return new_arr.in_units(unit)
        else:
            raise YTInvalidUnitEquivalence(equiv, self.units, unit)

    def list_equivalencies(self):
        """
        Lists the possible equivalencies associated with this YTArray or
        YTQuantity.
        """
        self.units.list_equivalencies()

    def has_equivalent(self, equiv):
        """
        Check to see if this YTArray or YTQuantity has an equivalent unit in
        *equiv*.
        """
        return self.units.has_equivalent(equiv)

    def ndarray_view(self):
        """
        Returns a view into the array, but as an ndarray rather than ytarray.

        Returns
        -------
        View of this array's data.
        """
        return self.view(np.ndarray)

    def to_ndarray(self):
        """
        Creates a copy of this array with the unit information stripped

        """
        return np.array(self)

    @classmethod
    def from_astropy(cls, arr, unit_registry=None):
        """
        Convert an AstroPy "Quantity" to a YTArray or YTQuantity.

        Parameters
        ----------
        arr : AstroPy Quantity
            The Quantity to convert from.
        unit_registry : yt UnitRegistry, optional
            A yt unit registry to use in the conversion. If one is not
            supplied, the default one will be used.
        """
        # Converting from AstroPy Quantity
        u = arr.unit
        ap_units = []
        for base, exponent in zip(u.bases, u.powers):
            unit_str = base.to_string()
            # we have to do this because AstroPy is silly and defines
            # hour as "h"
            if unit_str == "h": unit_str = "hr"
            ap_units.append("%s**(%s)" % (unit_str, Rational(exponent)))
        ap_units = "*".join(ap_units)
        if isinstance(arr.value, np.ndarray):
            return YTArray(arr.value, ap_units, registry=unit_registry)
        else:
            return YTQuantity(arr.value, ap_units, registry=unit_registry)


    def to_astropy(self, **kwargs):
        """
        Creates a new AstroPy quantity with the same unit information.
        """
        if _astropy.units is None:
            raise ImportError("You don't have AstroPy installed, so you can't convert to " +
                              "an AstroPy quantity.")
        return self.value*_astropy.units.Unit(str(self.units), **kwargs)

    @classmethod
    def from_pint(cls, arr, unit_registry=None):
        """
        Convert a Pint "Quantity" to a YTArray or YTQuantity.

        Parameters
        ----------
        arr : Pint Quantity
            The Quantity to convert from.
        unit_registry : yt UnitRegistry, optional
            A yt unit registry to use in the conversion. If one is not
            supplied, the default one will be used.

        Examples
        --------
        >>> from pint import UnitRegistry
        >>> import numpy as np
        >>> ureg = UnitRegistry()
        >>> a = np.random.random(10)
        >>> b = ureg.Quantity(a, "erg/cm**3")
        >>> c = yt.YTArray.from_pint(b)
        """
        p_units = []
        for base, exponent in arr._units.items():
            bs = convert_pint_units(base)
            p_units.append("%s**(%s)" % (bs, Rational(exponent)))
        p_units = "*".join(p_units)
        if isinstance(arr.magnitude, np.ndarray):
            return YTArray(arr.magnitude, p_units, registry=unit_registry)
        else:
            return YTQuantity(arr.magnitude, p_units, registry=unit_registry)

    def to_pint(self, unit_registry=None):
        """
        Convert a YTArray or YTQuantity to a Pint Quantity.

        Parameters
        ----------
        arr : YTArray or YTQuantity
            The unitful quantity to convert from.
        unit_registry : Pint UnitRegistry, optional
            The Pint UnitRegistry to use in the conversion. If one is not
            supplied, the default one will be used. NOTE: This is not
            the same as a yt UnitRegistry object.

        Examples
        --------
        >>> a = YTQuantity(4.0, "cm**2/s")
        >>> b = a.to_pint()
        """
        from pint import UnitRegistry
        if unit_registry is None:
            unit_registry = UnitRegistry()
        powers_dict = self.units.expr.as_powers_dict()
        units = []
        for unit, pow in powers_dict.items():
            # we have to do this because Pint doesn't recognize
            # "yr" as "year"
            if str(unit).endswith("yr") and len(str(unit)) in [2,3]:
                unit = str(unit).replace("yr","year")
            units.append("%s**(%s)" % (unit, Rational(pow)))
        units = "*".join(units)
        return unit_registry.Quantity(self.value, units)

    #
    # End unit conversion methods
    #

    def write_hdf5(self, filename, dataset_name=None, info=None, group_name=None):
        r"""Writes a YTArray to hdf5 file.

        Parameters
        ----------
        filename: string
            The filename to create and write a dataset to

        dataset_name: string
            The name of the dataset to create in the file.

        info: dictionary
            A dictionary of supplementary info to write to append as attributes
            to the dataset.

        group_name: string
            An optional group to write the arrays to. If not specified, the arrays
            are datasets at the top level by default.

        Examples
        --------
        >>> a = YTArray([1,2,3], 'cm')
        >>> myinfo = {'field':'dinosaurs', 'type':'field_data'}
        >>> a.write_hdf5('test_array_data.h5', dataset_name='dinosaurs',
        ...              info=myinfo)
        """
        from yt.utilities.on_demand_imports import _h5py as h5py
        from yt.extern.six.moves import cPickle as pickle
        if info is None:
            info = {}

        info['units'] = str(self.units)
        info['unit_registry'] = np.void(pickle.dumps(self.units.registry.lut))

        if dataset_name is None:
            dataset_name = 'array_data'

        f = h5py.File(filename)
        if group_name is not None:
            if group_name in f:
                g = f[group_name]
            else:
                g = f.create_group(group_name)
        else:
            g = f
        if dataset_name in g.keys():
            d = g[dataset_name]
            # Overwrite without deleting if we can get away with it.
            if d.shape == self.shape and d.dtype == self.dtype:
                d[...] = self
                for k in d.attrs.keys():
                    del d.attrs[k]
            else:
                del f[dataset_name]
                d = g.create_dataset(dataset_name, data=self)
        else:
            d = g.create_dataset(dataset_name, data=self)

        for k, v in info.items():
            d.attrs[k] = v
        f.close()

    @classmethod
    def from_hdf5(cls, filename, dataset_name=None, group_name=None):
        r"""Attempts read in and convert a dataset in an hdf5 file into a 
        YTArray.

        Parameters
        ----------
        filename: string
        The filename to of the hdf5 file.

        dataset_name: string
            The name of the dataset to read from.  If the dataset has a units
            attribute, attempt to infer units as well.

        group_name: string
            An optional group to read the arrays from. If not specified, the 
            arrays are datasets at the top level by default.

        """
        import h5py
        from yt.extern.six.moves import cPickle as pickle

        if dataset_name is None:
            dataset_name = 'array_data'

        f = h5py.File(filename)
        if group_name is not None:
            g = f[group_name]
        else:
            g = f
        dataset = g[dataset_name]
        data = dataset[:]
        units = dataset.attrs.get('units', '')
        if 'unit_registry' in dataset.attrs.keys():
            unit_lut = pickle.loads(dataset.attrs['unit_registry'].tostring())
        else:
            unit_lut = None
        f.close()
        registry = UnitRegistry(lut=unit_lut, add_default_symbols=False)
        return cls(data, units, registry=registry)

    #
    # Start convenience methods
    #

    @property
    def value(self):
        """Get a copy of the array data as a numpy ndarray"""
        return np.array(self)

    v = value

    @property
    def ndview(self):
        """Get a view of the array data."""
        return self.ndarray_view()

    d = ndview

    @property
    def unit_quantity(self):
        """Get a YTQuantity with the same unit as this array and a value of 
        1.0"""
        return YTQuantity(1.0, self.units)

    uq = unit_quantity

    @property
    def unit_array(self):
        """Get a YTArray filled with ones with the same unit and shape as this 
        array"""
        return np.ones_like(self)

    ua = unit_array

    def __getitem__(self, item):
        ret = super(YTArray, self).__getitem__(item)
        if ret.shape == ():
            return YTQuantity(ret, self.units, bypass_validation=True)
        else:
            if hasattr(self, 'units'):
                ret.units = self.units
            return ret

    #
    # Start operation methods
    #

    if LooseVersion(np.__version__) < LooseVersion('1.13.0'):

        def __add__(self, right_object):
            """
            Add this ytarray to the object on the right of the `+` operator. 
            Must check for the correct (same dimension) units.

            """
            ro = sanitize_units_add(self, right_object, "addition")
            return super(YTArray, self).__add__(ro)

        def __radd__(self, left_object):
            """ See __add__. """
            lo = sanitize_units_add(self, left_object, "addition")
            return super(YTArray, self).__radd__(lo)

        def __iadd__(self, other):
            """ See __add__. """
            oth = sanitize_units_add(self, other, "addition")
            np.add(self, oth, out=self)
            return self

        def __sub__(self, right_object):
            """
            Subtract the object on the right of the `-` from this ytarray. Must
            check for the correct (same dimension) units.

            """
            ro = sanitize_units_add(self, right_object, "subtraction")
            return super(YTArray, self).__sub__(ro)

        def __rsub__(self, left_object):
            """ See __sub__. """
            lo = sanitize_units_add(self, left_object, "subtraction")
            return super(YTArray, self).__rsub__(lo)

        def __isub__(self, other):
            """ See __sub__. """
            oth = sanitize_units_add(self, other, "subtraction")
            np.subtract(self, oth, out=self)
            return self

        def __neg__(self):
            """ Negate the data. """
            return super(YTArray, self).__neg__()

        def __mul__(self, right_object):
            """
            Multiply this YTArray by the object on the right of the `*` 
            operator. The unit objects handle being multiplied.

            """
            ro = sanitize_units_mul(self, right_object)
            return super(YTArray, self).__mul__(ro)

        def __rmul__(self, left_object):
            """ See __mul__. """
            lo = sanitize_units_mul(self, left_object)
            return super(YTArray, self).__rmul__(lo)

        def __imul__(self, other):
            """ See __mul__. """
            oth = sanitize_units_mul(self, other)
            np.multiply(self, oth, out=self)
            return self

        def __div__(self, right_object):
            """
            Divide this YTArray by the object on the right of the `/` operator.

            """
            ro = sanitize_units_mul(self, right_object)
            return super(YTArray, self).__div__(ro)

        def __rdiv__(self, left_object):
            """ See __div__. """
            lo = sanitize_units_mul(self, left_object)
            return super(YTArray, self).__rdiv__(lo)

        def __idiv__(self, other):
            """ See __div__. """
            oth = sanitize_units_mul(self, other)
            
np.divide(self, oth, out=self)