def entrymean(self, m, axis=None):
        """Average a matrix over the given axis. If the axis is None,
        average over both rows and columns, returning a scalar.
        (via some SciPy function)
        """
        # Mimic numpy's casting.  The int32/int64 check works around numpy
        # 1.5.x behavior of np.issubdtype, see gh-2677.
        if (np.issubdtype(m.dtype, np.float_) or
                np.issubdtype(m.dtype, np.int_) or
                    m.dtype in [np.dtype('int32'), np.dtype('int64')] or
                np.issubdtype(m.dtype, np.bool_)):
            res_dtype = np.float_
        elif np.issubdtype(m.dtype, np.complex_):
            res_dtype = np.complex_
        else:
            res_dtype = m.dtype

        m = m.astype(res_dtype)
        mu = m.sum(None) / m.getnnz()
        # if user or item has no ratings (stripped from training data), set to 0
        b_i = m.sum(0)
        b_u = m.sum(1)
        with np.errstate(invalid='ignore'):
            b_i = (b_i / m.getnnz(axis=0)) - mu
            b_u = (b_u.T / m.getnnz(axis=1)) - mu
        b_i[np.isnan(b_i)] = 0
        b_u[np.isnan(b_u)] = 0

        return mu, np.array(b_i)[0], np.array(b_u)[0]

    def mean(self, axis=None):
        """Average the matrix over the given axis.  If the axis is None,
        average over both rows and columns, returning a scalar.
        """
        # Mimic numpy's casting.  The int32/int64 check works around numpy
        # 1.5.x behavior of np.issubdtype, see gh-2677.
        if (np.issubdtype(self.dtype, np.float_) or
            np.issubdtype(self.dtype, np.int_) or
            self.dtype in [np.dtype('int32'), np.dtype('int64')] or
            np.issubdtype(self.dtype, np.bool_)):
                res_dtype = np.float_
        elif np.issubdtype(self.dtype, np.complex_):
            res_dtype = np.complex_
        else:
            res_dtype = self.dtype

        if axis is None:
            return self.sum(None) * 1.0 / (self.shape[0]*self.shape[1])

        if axis < 0:
            axis += 2
        if axis == 0:
            mean = self.astype(res_dtype).sum(0)
            mean *= 1.0 / self.shape[0]
            return mean
        elif axis == 1:
            mean = self.astype(res_dtype).sum(1)
            mean *= 1.0 / self.shape[1]
            return mean
        else:
            raise ValueError("axis out of bounds")

def assert_image_equal(actual, expected):
    if np.issubdtype(actual.dtype, np.integer):
        assert_equal(actual, expected)
    else:
        if np.issubdtype(expected.dtype, np.integer):
            expected = expected/float(np.iinfo(expected.dtype).max)
        assert_allclose(actual, expected, atol=1/256.)

 def transform_python_types(self, obj):
     """handle special scalars, default to default json encoder
     """
     # Pandas Timestamp
     if is_pandas and isinstance(obj, pd.tslib.Timestamp):
         return obj.value / 10**6.0  #nanosecond to millisecond
     elif np.issubdtype(type(obj), np.float):
         return float(obj)
     elif np.issubdtype(type(obj), np.int):
         return int(obj)
     elif np.issubdtype(type(obj), np.bool_):
         return bool(obj)
     # Datetime
     # datetime is a subclass of date.
     elif isinstance(obj, dt.datetime):
         return calendar.timegm(obj.timetuple()) * 1000. + obj.microsecond / 1000.
     # Date
     elif isinstance(obj, dt.date):
         return calendar.timegm(obj.timetuple()) * 1000.
     # Numpy datetime64
     elif isinstance(obj, np.datetime64):
         epoch_delta = obj - np.datetime64('1970-01-01T00:00:00Z')
         return (epoch_delta / np.timedelta64(1, 'ms'))
     # Time
     elif isinstance(obj, dt.time):
         return (obj.hour * 3600 + obj.minute * 60 + obj.second) * 1000 + obj.microsecond / 1000.
     elif is_dateutil and isinstance(obj, relativedelta):
         return dict(years=obj.years, months=obj.months, days=obj.days, hours=obj.hours,
             minutes=obj.minutes, seconds=obj.seconds, microseconds=obj.microseconds)
     # Decimal
     elif isinstance(obj, decimal.Decimal):
         return float(obj)
     else:
         return super(BokehJSONEncoder, self).default(obj)

 def _set_dtype(self, dtype, union=False):
     if np.issubdtype(dtype, np.complexfloating) \
            or np.issubdtype(self.dtype, np.complexfloating):
         self.dtype = np.complex_
     else:
         if not union or self.dtype != np.complex_:
             self.dtype = np.float_

def find_linear_scale(data):
    scale = []
    scale_name = []
    linear_scale = False
    longest = None
    if type(data.columns) == pd.MultiIndex:
        for n, l in enumerate(data.columns.levels):
            if l.dtype == np.dtype('O'):  # Object; maybe str?
                if longest is None or len(l) > longest:
                    longest = len(l)

            elif np.issubdtype(l.dtype, np.integer) or np.issubdtype(l.dtype, np.float):
                linear_scale = True
                scale = [v[n] for v in data.columns.values]
                scale_name = data.columns.names[n]

                if np.issubdtype(l.dtype, np.float):
                    # Prefer float scales, assume more accurate
                    break
    else:
        scale = []
        linear_scale = True
        for x in data.columns.values:
            try:
                scale.append(float(x))
            except:
                linear_scale = False
                break

    return scale, linear_scale, scale_name

    def to_json(o, level=0):
        ''' format JSON with no line break between list items '''
        INDENT = 3
        SPACE = " "
        NEWLINE = "\n"

        ret = ""
        if isinstance(o, dict):
            ret += "{" + NEWLINE
            comma = ""
            for k,v in o.iteritems():
                ret += comma
                comma = ",\n"
                ret += SPACE * INDENT * (level+1)
                ret += '"' + str(k) + '":' + SPACE
                ret += Utils.to_json(v, level + 1)

            ret += NEWLINE + SPACE * INDENT * level + "}"
        elif isinstance(o, basestring):
            ret += '"' + o + '"'
        elif isinstance(o, list):
            ret += "[" + ",".join([Utils.to_json(e, level+1) for e in o]) + "]"
        elif isinstance(o, bool):
            ret += "true" if o else "false"
        elif isinstance(o, int):
            ret += str(o)
        elif isinstance(o, float):
            ret += '%.7g' % o
        elif isinstance(o, numpy.ndarray) and numpy.issubdtype(o.dtype, numpy.integer):
            ret += "[" + ','.join(map(str, o.flatten().tolist())) + "]"
        elif isinstance(o, numpy.ndarray) and numpy.issubdtype(o.dtype, numpy.inexact):
            ret += "[" + ','.join(map(lambda x: '%.7g' % x, o.flatten().tolist())) + "]"
        else:
            raise TypeError("Unknown type '%s' for json serialization" % str(type(o)))
        return ret

  def DtypeToNumberConverter(self, dtype):
    """Converts a Numpy dtype to a converter method if applicable.

      The converter method takes in a numpy array of objects of the provided
      dtype
      and returns a numpy array of the numbers backing that object for
      statistical
      analysis. Returns None if no converter is necessary.

    Args:
      dtype: The numpy dtype to make a converter for.

    Returns:
      The converter method or None.
    """
    if np.issubdtype(dtype, np.datetime64):

      def DatetimesToNumbers(dt_list):
        return np.array([pd.Timestamp(dt).value for dt in dt_list])

      return DatetimesToNumbers
    elif np.issubdtype(dtype, np.timedelta64):

      def TimedetlasToNumbers(td_list):
        return np.array([pd.Timedelta(td).value for td in td_list])

      return TimedetlasToNumbers
    else:
      return None

def test_apply_loop_invariant_optimisation_integer():
    variables = {'v': Variable('v', scalar=False),
                 'N': Constant('N', 10),
                 'b': Variable('b', scalar=True, dtype=int),
                 'c': Variable('c', scalar=True, dtype=int),
                 'd': Variable('d', scalar=True, dtype=int),
                 'y': Variable('y', scalar=True, dtype=float),
                 'z': Variable('z', scalar=True, dtype=float),
                 'w': Variable('w', scalar=True, dtype=float),
                 }
    statements = [Statement('v', '=', 'v % (2*3*N)', '', np.float32),
                  # integer version doesn't get rewritten but float version does
                  Statement('a', ':=', 'b//(c//d)', '', int),
                  Statement('x', ':=', 'y/(z/w)', '', float),
                  ]
    scalar, vector = optimise_statements([], statements, variables)
    assert len(scalar) == 3
    assert np.issubdtype(scalar[0].dtype, np.signedinteger)
    assert scalar[0].var == '_lio_1'
    expr = scalar[0].expr.replace(' ', '')
    assert expr=='6*N' or expr=='N*6'
    assert np.issubdtype(scalar[1].dtype, np.signedinteger)
    assert scalar[1].var == '_lio_2'
    expr = scalar[1].expr.replace(' ', '')
    assert expr=='b//(c//d)'
    assert np.issubdtype(scalar[2].dtype, np.floating)
    assert scalar[2].var == '_lio_3'
    expr = scalar[2].expr.replace(' ', '')
    assert expr=='(y*w)/z' or expr=='(w*y)/z'

def load_image(image_file):
    """
    Loads an analyze/nifti image, generally for 3D images. 
    Casts input as 32-bit float (if float) or 32-bit uint (if int).
    
    Paramaters
    ----------
    image_file: str
        path to data
    
    Returns
    -------
    dat: nparray
        Image as numpy array (i.e., 3D array)
    """
    img = nib.load(image_file)
    dat = img.get_data()
    # Ensure that data is cast as at least 32-bit
    if np.issubdtype(dat.dtype, float):
        dat = dat.astype("float32")
        # Check for negative values
        if (dat < 0).any():
            print "found negative values, setting to zero (see file: %s)" % image_file
            dat[dat < 0] = 0
    elif np.issubdtype(dat.dtype, int):
        dat = dat.astype("uint32")
    else:
        msg = "Error: Unknown datatype %s" % dat.dtype
        print msg
        raise Exception(msg)
    return dat

def from_bcolz(x, chunksize=None, categorize=True, index=None, **kwargs):
    """ Read dask Dataframe from bcolz.ctable

    Parameters
    ----------

    x : bcolz.ctable
        Input data
    chunksize : int (optional)
        The size of blocks to pull out from ctable.  Ideally as large as can
        comfortably fit in memory
    categorize : bool (defaults to True)
        Automatically categorize all string dtypes
    index : string (optional)
        Column to make the index

    See Also
    --------

    from_array: more generic function not optimized for bcolz
    """
    import dask.array as da
    import bcolz
    if isinstance(x, (str, unicode)):
        x = bcolz.ctable(rootdir=x)
    bc_chunklen = max(x[name].chunklen for name in x.names)
    if chunksize is None and bc_chunklen > 10000:
        chunksize = bc_chunklen

    categories = dict()
    if categorize:
        for name in x.names:
            if (np.issubdtype(x.dtype[name], np.string_) or
                    np.issubdtype(x.dtype[name], np.unicode_) or
                    np.issubdtype(x.dtype[name], np.object_)):
                a = da.from_array(x[name], chunks=(chunksize * len(x.names),))
                categories[name] = da.unique(a)

    columns = tuple(x.dtype.names)
    divisions = (0,) + tuple(range(-1, len(x), chunksize))[1:]
    if divisions[-1] != len(x) - 1:
        divisions = divisions + (len(x) - 1,)
    new_name = 'from_bcolz' + next(tokens)
    dsk = dict(((new_name, i),
                (dataframe_from_ctable,
                 x,
                 (slice(i * chunksize, (i + 1) * chunksize),),
                 None, categories))
               for i in range(0, int(ceil(len(x) / chunksize))))

    result = DataFrame(dsk, new_name, columns, divisions)

    if index:
        assert index in x.names
        a = da.from_array(x[index], chunks=(chunksize * len(x.names),))
        q = np.linspace(0, 100, len(x) // chunksize + 2)
        divisions = da.percentile(a, q).compute()
        return set_partition(result, index, divisions, **kwargs)
    else:
        return result

    def load_onto_vtk(self, vtk_data):
        """ Load the stored information onto a vtk data container.

        Parameters
        ----------
        vtk_data : vtkPointData or vtkCellData
            The vtk container to load the value onto.

        Data are loaded onto the vtk container based on their data
        type. The name of the added array is the name of the CUBA key
        (i.e. :samp:`{CUBA}.name`). Currently only scalars and three
        dimensional vectors are supported.

        """
        def replacer(data):
            return nan if data is None else data

        for cuba in self.keys:
            default = dummy_cuba_value(cuba)
            if (numpy.issubdtype(type(default), numpy.float) or
                    numpy.issubdtype(type(default), numpy.int)):
                data = numpy.array(self._data[cuba], dtype=float)
                index = vtk_data.add_array(data)
                vtk_data.get_array(index).name = cuba.name
            elif isinstance(default, numpy.ndarray) and default.size == 3:
                nan = numpy.array([None, None, None], dtype=float)
                data = numpy.array(
                    tuple(replacer(data) for data in self._data[cuba]),
                    dtype=numpy.float)
                index = vtk_data.add_array(data)
                vtk_data.get_array(index).name = cuba.name
            else:
                message = 'property {!r} is currently ignored'
                warnings.warn(message.format(cuba))

def has_inf_or_nan(datum, tensor):
  """A predicate for whether a tensor consists of any bad numerical values.

  This predicate is common enough to merit definition in this module.
  Bad numerical values include nans and infs.
  The signature of this function follows the requiremnet of DebugDumpDir's
  find() method.

  Args:
    datum: (DebugTensorDatum) Datum metadata.
    tensor: (numpy.ndarray or None) Value of the tensor. None represents
      an uninitialized tensor.

  Returns:
    (bool) True if and only if tensor consists of any nan or inf values.
  """

  _ = datum  # Datum metadata is unused in this predicate.
  if tensor is None:
    # Uninitialized tensor doesn't have bad numerical values.
    return False
  elif (np.issubdtype(tensor.dtype, np.float) or
        np.issubdtype(tensor.dtype, np.complex) or
        np.issubdtype(tensor.dtype, np.integer)):
    return np.any(np.isnan(tensor)) or np.any(np.isinf(tensor))
  else:
    return False

    def mean(self, axis=None):
        """Average the matrix over the given axis.  If the axis is None,
        average over both rows and columns, returning a scalar.
        """
        # Mimic numpy's casting.
        if (np.issubdtype(self.dtype, np.float_) or
                np.issubdtype(self.dtype, np.integer) or
                np.issubdtype(self.dtype, np.bool_)):
            res_dtype = np.float_
        elif np.issubdtype(self.dtype, np.complex_):
            res_dtype = np.complex_
        else:
            res_dtype = self.dtype

        if axis is None:
            return self.sum(None) * 1.0 / (self.shape[0]*self.shape[1])

        if axis < 0:
            axis += 2
        if axis == 0:
            mean = self.astype(res_dtype).sum(0)
            mean *= 1.0 / self.shape[0]
            return mean
        elif axis == 1:
            mean = self.astype(res_dtype).sum(1)
            mean *= 1.0 / self.shape[1]
            return mean
        else:
            raise ValueError("axis out of bounds")

    def _diff(self):
        if self.a.shape != self.b.shape:
            self.diff_dimensions = (self.a.shape, self.b.shape)
            # Don't do any further comparison if the dimensions differ
            # TODO: Perhaps we could, however, diff just the intersection
            # between the two images
            return

        # Find the indices where the values are not equal
        # If neither a nor b are floating point, ignore self.tolerance
        if not ((np.issubdtype(self.a.dtype, float) or
                 np.issubdtype(self.a.dtype, complex)) or
                (np.issubdtype(self.b.dtype, float) or
                 np.issubdtype(self.b.dtype, complex))):
            tolerance = 0
        else:
            tolerance = self.tolerance

        diffs = where_not_allclose(self.a, self.b, atol=0.0, rtol=tolerance)

        self.diff_total = len(diffs[0])

        if self.diff_total == 0:
            # Then we're done
            return

        if self.numdiffs < 0:
            numdiffs = self.diff_total
        else:
            numdiffs = self.numdiffs

        self.diff_pixels = [(idx, (self.a[idx], self.b[idx]))
                            for idx in islice(izip(*diffs), 0, numdiffs)]
        self.diff_ratio = float(self.diff_total) / float(len(self.a.flat))

def half_fft_convolve(in1, in2, size, mode = 'full', return_type='real'):
    """
    Rewrite of fftconvolve from scipy.signal ((c) Travis Oliphant 1999-2002)
    to deal with fft convolution where one signal is not fft transformed
    and the other one is.  Application is, for example, in a loop where
    convolution happens repeatedly with different kernels over the same
    signal.  First input is not transformed, second input is.
    """
    s1 = np.array(in1.shape)
    s2 = size - s1 + 1
    complex_result = (np.issubdtype( in1.dtype, np.complex) or
                      np.issubdtype( in2.dtype, np.complex) )

    # Always use 2**n-sized FFT
    fsize = 2 **np.ceil( np.log2( size) )
    IN1 = fftn(in1, fsize)
    IN1 *= in2
    fslice = tuple( [slice( 0, int(sz)) for sz in size] )
    ret = ifftn(IN1)[fslice].copy()
    del IN1
    if not complex_result:
        ret = ret.real
    if return_type == 'real':
        ret = ret.real
    if mode == 'full':
        return ret
    elif mode == 'same':
        if np.product(s1, axis=0) > np.product(s2, axis=0):
            osize = s1
        else:
            osize = s2
        return _centered(ret, osize)
    elif mode == 'valid':
        return _centered(ret, abs(s2 - s1) + 1)

def into(a, b, **kwargs):
    dialect = b.dialect.copy()
    del dialect['lineterminator']
    dates = [i for i, typ in enumerate(b.schema[0].types)
               if 'date' in str(typ)]
    schema = b.schema
    if '?' in str(schema):
        schema = dshape(str(schema).replace('?', ''))

    dtypes = valmap(to_numpy_dtype, schema[0].dict)

    datenames = [name for name in dtypes
                      if np.issubdtype(dtypes[name], np.datetime64)]

    dtypes = dict((k, v) for k, v in dtypes.items()
                         if not np.issubdtype(v, np.datetime64))

    if 'strict' in dialect:
        del dialect['strict']

    # Pass only keyword arguments appropriate for read_csv
    kws = keywords(pd.read_csv)
    options = toolz.merge(dialect, kwargs)
    options = toolz.keyfilter(lambda k: k in kws, options)

    if b.open == gzip.open:
        options['compression'] = 'gzip'

    return pd.read_csv(b.path,
                       skiprows=1 if b.header else 0,
                       dtype=dtypes,
                       parse_dates=datenames,
                       names=b.columns,
                       **options)

    def _argsortData(self, data, order):
        if data.ndim == 1:
            indices = np.argsort(data, kind='mergesort')
            if order == Qt.DescendingOrder:
                indices = indices[::-1]
            # Always sort NaNs last
            if np.issubdtype(data.dtype, np.number):
                indices = np.roll(indices, -np.isnan(data).sum())
        else:
            assert np.issubdtype(data.dtype, np.number), \
                'We do not deal with non numeric values in sorting by ' \
                'multiple values'
            if order == Qt.DescendingOrder:
                data[:, -1] = -data[:, -1]

            # In order to make sure NaNs always appear at the end, insert a
            # indicator whether NaN or not. Note that the data array must
            # contain an empty column of zeros at index -2 since inserting an
            # extra column after the fact can result in a MemoryError for data
            # with a large amount of variables
            assert np.all(data[:, -2] == 0), \
                'Add an empty column of zeros at index -2 to accomodate NaNs'
            np.isnan(data[:, -1], out=data[:, -2])

            indices = np.lexsort(np.flip(data.T, axis=0))

        return indices

    def from_numpy_dtype(self, dt):
        """
        From Numpy dtype.

        >>> from datashape import CType
        >>> from numpy import dtype
        >>> CType.from_numpy_dtype(dtype('int32'))
        ctype("int32")
        >>> CType.from_numpy_dtype(dtype('i8'))
        ctype("int64")
        >>> CType.from_numpy_dtype(dtype('M8'))
        DateTime(None)
        >>> CType.from_numpy_dtype(dtype('U30'))
        ctype("string[30, 'U32']")
        """
        try:
            return Type.lookup_type(dt.name)
        except KeyError:
            pass
        if np.issubdtype(dt, np.datetime64):
            unit, _ = np.datetime_data(dt)
            defaults = {'D': date_, 'Y': date_, 'M': date_, 'W': date_}
            return defaults.get(unit, datetime_)
        elif np.issubdtype(dt, np.timedelta64):
            unit, _ = np.datetime_data(dt)
            return TimeDelta(unit=unit)
        elif np.issubdtype(dt, np.unicode_):
            return String(dt.itemsize // 4, 'U32')
        elif np.issubdtype(dt, np.str_) or np.issubdtype(dt, np.bytes_):
            return String(dt.itemsize, 'ascii')
        raise NotImplementedError("NumPy datatype %s not supported" % dt)

def isEqual(left, right, eps=None, masked_equal=True):
  ''' This function checks if two numpy arrays or scalars are equal within machine precision, and returns a scalar logical. '''
  diff_type = "Both arguments to function 'isEqual' must be of the same class!"
  if isinstance(left,np.ndarray):
    # ndarray
    if not isinstance(right,np.ndarray): raise TypeError(diff_type)
    if not left.dtype==right.dtype:
      right = right.astype(left.dtype) # casting='same_kind' doesn't work...
    if np.issubdtype(left.dtype, np.inexact): # also catch float32 etc
      if eps is None: return ma.allclose(left, right, masked_equal=masked_equal)
      else: return ma.allclose(left, right, masked_equal=masked_equal, atol=eps)
    elif np.issubdtype(left.dtype, np.integer) or np.issubdtype(left.dtype, np.bool):
      return np.all( left == right ) # need to use numpy's all()
  elif isinstance(left,(float,np.inexact)):
    # numbers
    if not isinstance(right,(float,np.inexact)): raise TypeError(diff_type)
    if eps is None: eps = 100.*floateps # default
    if ( isinstance(right,float) or isinstance(right,float) ) or left.dtype.itemsize == right.dtype.itemsize: 
      return np.absolute(left-right) <= eps
    else:
      if left.dtype.itemsize < right.dtype.itemsize: right = left.dtype.type(right)
      else: left = right.dtype.type(left)
      return np.absolute(left-right) <= eps  
  elif isinstance(left,(int,bool,np.integer,np.bool)):
    # logicals
    if not isinstance(right,(int,bool,np.integer,np.bool)): raise TypeError(diff_type)
    return left == right
  else: raise TypeError(left)