def float_test(self, dtype, significant=None):
     colname = 'col_%s' % dtype.__name__
     self.writeread(dtype)
     before, after = self.table_orig.data[colname], self.table_new.data[colname]
     self.assertEqual(before.shape, after.shape)
     self.assertEqual(before.dtype.type, after.dtype.type)
     if before.ndim == 1:
         for i in range(before.shape[0]):
             if(np.isnan(before[i])):
                 self.failUnless(np.isnan(after[i]))
             elif(np.isinf(before[i])):
                 self.failUnless(np.isinf(after[i]))
             else:
                 if significant:
                     self.assertAlmostEqualSig(before[i], after[i], significant=significant)
                 else:
                     self.assertEqual(before[i], after[i])
     else:
         for i in range(before.shape[0]):
             for j in range(before.shape[1]):
                 if(np.isnan(before[i, j])):
                     self.failUnless(np.isnan(after[i, j]))
                 elif(np.isinf(before[i, j])):
                     self.failUnless(np.isinf(after[i, j]))
                 else:
                     if significant:
                         self.assertAlmostEqualSig(before[i, j], after[i, j], significant=significant)
                     else:
                         self.assertEqual(before[i, j], after[i, j])

    def test_nan_inf(self):
        # Not-a-number
        q = u.Quantity('nan', unit='cm')
        assert np.isnan(q.value)

        q = u.Quantity('NaN', unit='cm')
        assert np.isnan(q.value)

        q = u.Quantity('-nan', unit='cm')  # float() allows this
        assert np.isnan(q.value)

        q = u.Quantity('nan cm')
        assert np.isnan(q.value)
        assert q.unit == u.cm

        # Infinity
        q = u.Quantity('inf', unit='cm')
        assert np.isinf(q.value)

        q = u.Quantity('-inf', unit='cm')
        assert np.isinf(q.value)

        q = u.Quantity('inf cm')
        assert np.isinf(q.value)
        assert q.unit == u.cm

        q = u.Quantity('Infinity', unit='cm')  # float() allows this
        assert np.isinf(q.value)

        # make sure these strings don't parse...
        with pytest.raises(TypeError):
            q = u.Quantity('', unit='cm')

        with pytest.raises(TypeError):
            q = u.Quantity('spam', unit='cm')

 def _get_sum(self):
     """Compute sum of non NaN / Inf values in the array."""
     try:
         return self._sum
     except AttributeError:
         self._sum = self.no_nan.sum()
         # The following 2 lines are needede as in Python 3.3 with NumPy
         # 1.7.1, numpy.ndarray and numpy.memmap aren't hashable.
         if type(self._sum) is numpy.memmap:
             self._sum = numpy.asarray(self._sum).item()
         if self.has_nan and self.no_nan.mask.all():
             # In this case the sum is not properly computed by numpy.
             self._sum = 0
         if numpy.isinf(self._sum) or numpy.isnan(self._sum):
             # NaN may happen when there are both -inf and +inf values.
             if self.has_nan:
                 # Filter both NaN and Inf values.
                 mask = self.no_nan.mask + numpy.isinf(self[1])
             else:
                 # Filter only Inf values.
                 mask = numpy.isinf(self[1])
             if mask.all():
                 self._sum = 0
             else:
                 self._sum = numpy.ma.masked_array(self[1], mask).sum()
             # At this point there should be no more NaN.
             assert not numpy.isnan(self._sum)
     return self._sum

def clean_invalid(x,y,min_x=-numpy.inf,min_y=-numpy.inf,max_x=numpy.inf,max_y=numpy.inf):
    """Remove corresponding values from x and y when one or both of those is `nan` or `inf`,
    and optionally truncate values to minima and maxima

    Parameters
    ----------
    x, y : :class:`numpy.ndarray` or list
        Pair arrays or lists of corresponding numbers

    min_x, min_y, max_x, max_y : number, optional
        If supplied, set values below `min_x` to `min_x`, values larger
        than `max_x` to `max_x` and so for `min_y` and `max_y`

    Returns
    -------
    :class:`numpy.ndarray`
        A shortened version of `x`, excluding invalid values

    :class:`numpy.ndarray`
        A shortened version of `y`, excluding invalid values
    """
    x = numpy.array(x).astype(float)
    y = numpy.array(y).astype(float)

    x[x < min_x] = min_x
    x[x > max_x] = max_x
    y[y < min_y] = min_y
    y[y > max_y] = max_y
    
    newmask = numpy.isinf(x) | numpy.isnan(x) | numpy.isinf(y) | numpy.isnan(y) 
    x = x[~newmask]
    y = y[~newmask]


    return x,y 

    def __call__(self, value, clip=None):
        if clip is None:
            clip = self.clip

        if cbook.iterable(value):
            vtype = 'array'
            val = np.ma.asarray(value).astype(np.float)
        else:
            vtype = 'scalar'
            val = np.ma.array([value]).astype(np.float)

        val = np.ma.masked_where(np.isinf(val.data),val)

        self.autoscale_None(val)
        vmin, vmax = float(self.vmin), float(self.vmax)
        if vmin > vmax:
            raise ValueError("minvalue must be less than or equal to maxvalue")
        elif vmin<=0:
            raise ValueError("values must all be positive")
        elif vmin==vmax:
            return type(value)(0.0 * np.asarray(value))
        else:
            if clip:
                mask = np.ma.getmask(val)
                val = np.ma.array(np.clip(val.filled(vmax), vmin, vmax),
                                   mask=mask)
            result = (np.ma.log(val)-np.log(vmin))/(np.log(vmax)-np.log(vmin))
            result.data[result.data<0]=0.0
            result.data[result.data>1]=1.0
            result[np.isinf(val.data)] = -np.inf
            if result.mask is not np.ma.nomask:
                result.mask[np.isinf(val.data)] = False
        if vtype == 'scalar':
            result = result[0]
        return result

    def sample_representer_points(self):
        # Sample representer points only in the
        # configuration space by setting all environmental
        # variables to 1
        D = np.where(self.is_env == 0)[0].shape[0]

        lower = self.lower[np.where(self.is_env == 0)]
        upper = self.upper[np.where(self.is_env == 0)]

        self.sampling_acquisition.update(self.model)

        for i in range(5):
            restarts = np.random.uniform(low=lower,
                                         high=upper,
                                         size=(self.Nb, D))
            sampler = emcee.EnsembleSampler(self.Nb, D,
                                        self.sampling_acquisition_wrapper)

            self.zb, self.lmb, _ = sampler.run_mcmc(restarts, 50)
            if not np.any(np.isinf(self.lmb)):
                break
            else:
                print("Infinity")
        if np.any(np.isinf(self.lmb)):
            raise ValueError("Could not sample valid representer points! LogEI is -infinity")
        if len(self.zb.shape) == 1:
            self.zb = self.zb[:, None]
        if len(self.lmb.shape) == 1:
            self.lmb = self.lmb[:, None]

        # Project representer points to subspace
        proj = np.ones([self.zb.shape[0],
                    self.upper[self.is_env == 1].shape[0]])
        proj *= self.upper[self.is_env == 1].shape[0]
        self.zb = np.concatenate((self.zb, proj), axis=1)

def contains_inf(arr):
    """
    Test whether a numpy.ndarray contains any `np.inf` values.

    Parameters
    ----------
    arr : np.ndarray

    Returns
    -------
    contains_inf : bool
        `True` if the array contains any `np.inf` values, `False` otherwise.

    Notes
    -----
    Tests for the presence of `np.inf`'s by determining whether the
    values returned by `np.nanmin(arr)` and `np.nanmax(arr)` are finite.
    This approach is more memory efficient than the obvious alternative,
    calling `np.any(np.isinf(ndarray))`, which requires the construction of a
    boolean array with the same shape as the input array.
    """
    if isinstance(arr, theano.gof.type.CDataType._cdata_type):
        return False
    elif isinstance(arr, np.random.mtrand.RandomState):
        return False
    return np.isinf(np.nanmax(arr)) or np.isinf(np.nanmin(arr))

    def common_limits(datasets, default_min=0, default_max=0):
        """Find the global maxima and minima of a list of datasets.

        Parameters
        ----------
        datasets : `iterable`
            list (or any other iterable) of data arrays to analyse.

        default_min : `float`, optional
            fall-back minimum value if datasets are all empty.

        default_max : `float`, optional
            fall-back maximum value if datasets are all empty.

        Returns
        -------
        (min, max) : `float`
            2-tuple of common minimum and maximum over all datasets.
        """
        from glue import iterutils
        if isinstance(datasets, numpy.ndarray) or not iterable(datasets[0]):
            datasets = [datasets]
        max_stat = max(list(iterutils.flatten(datasets)) + [-numpy.inf])
        min_stat = min(list(iterutils.flatten(datasets)) + [numpy.inf])
        if numpy.isinf(-max_stat):
            max_stat = default_max
        if numpy.isinf(min_stat):
            min_stat = default_min
        return min_stat, max_stat

 def _check_for_infinities(self, tif):
     try:
         if np.any(np.isinf(tif)):
             tif[np.isinf(tif)] = 0
             g.alert('Some array values were inf. Setting those values to 0')
     except MemoryError:
         pass

def circumcircle(P1,P2,P3):
    ''' 
    Adapted from:
    http://local.wasp.uwa.edu.au/~pbourke/geometry/circlefrom3/Circle.cpp
    '''
    delta_a = P2 - P1
    delta_b = P3 - P2
    if np.abs(delta_a[0]) <= 0.000000001 and np.abs(delta_b[1]) <= 0.000000001:
        center_x = 0.5*(P2[0] + P3[0])
        center_y = 0.5*(P1[1] + P2[1])
    else:
        aSlope = delta_a[1]/delta_a[0]
        bSlope = delta_b[1]/delta_b[0]

        if aSlope == 0.0:
            aSlope = 1E-6

        if bSlope == 0.0:
            bSlope = 1E-6

        if np.isinf(aSlope):
            aSlope = 1E6

        if np.isinf(bSlope):
            bSlope = 1E6

        if np.abs(aSlope-bSlope) <= 0.000000001:
            return None
        center_x= (aSlope*bSlope*(P1[1] - P3[1]) + bSlope*(P1[0] + P2 [0]) \
                        - aSlope*(P2[0]+P3[0]) )/(2* (bSlope-aSlope) )
        center_y = -1*(center_x - (P1[0]+P2[0])/2)/aSlope +  (P1[1]+P2[1])/2;
    return center_x, center_y

def __compare( verify_obj, obj, nsig, ndec ):
  if isinstance(verify_obj,tuple):
    if len(verify_obj) == len(obj):
      return all([ __compare( vo, o, nsig, ndec, title='#%d' % i )
        for i, (vo,o) in enumerate( zip( verify_obj, obj ) ) ])
    log.error( 'non matching lenghts: %d != %d' % ( len(verify_obj), len(obj) ) )
  elif not isinstance(verify_obj,float):
    if verify_obj == obj:
      return True
    log.error( 'non equal: %s != %d' % ( obj, verify_obj ) )
  elif numpy.isnan(verify_obj):
    if numpy.isnan(obj):
      return True
    log.error( 'expected nan: %s' % obj )
  elif numpy.isinf(verify_obj):
    if numpy.isinf(obj):
      return True
    log.error( 'expected inf: %s' % obj )
  else:
    if verify_obj:
      n = numeric.floor( numpy.log10( abs(verify_obj) ) )
      N = max( n-(nsig-1), -ndec )
    else:
      N = -ndec
    maxerr = .5 * 10.**N
    if abs(verify_obj-obj) <= maxerr:
      return True
    log.error( 'non equal to %s digits: %e != %e' % ( nsig, obj, verify_obj ) )
  return False

def weighted_mean(_line):
    max_weight = 50
    
    # print _line.shape
    
    median_2d = bottleneck.nanmedian(_line, axis=1).reshape(_line.shape[0],1).repeat(_line.shape[1], axis=1)
    std = bottleneck.nanstd(_line, axis=1)
    std_2d = std.reshape(_line.shape[0],1).repeat(_line.shape[1], axis=1)
    
    weight_2d = numpy.fabs(std_2d / (_line - median_2d))
#    weight_2d[weight_2d > max_weight] = max_weight
    weight_2d[numpy.isinf(weight_2d)] = max_weight
    
    for i in range(3):
        avg = bottleneck.nansum(_line*weight_2d, axis=1)/bottleneck.nansum(weight_2d, axis=1)
        avg_2d = avg.reshape(_line.shape[0],1).repeat(_line.shape[1], axis=1)
        
        std = numpy.sqrt(bottleneck.nansum(((_line - avg_2d)**2 * weight_2d), axis=1)/bottleneck.nansum(weight_2d, axis=1))
        std_2d = std.reshape(_line.shape[0],1).repeat(_line.shape[1], axis=1)
        
        weight_2d = numpy.fabs(std_2d / (_line - avg_2d))
        #weight_2d[weight_2d > max_weight] = max_weight
        weight_2d[numpy.isinf(weight_2d)] = max_weight
    
    return bottleneck.nansum(_line*weight_2d, axis=1)/bottleneck.nansum(weight_2d, axis=1)

    def test_nans_infs(self):
        oldsettings = np.seterr(all='ignore')
        try:
            # Check some of the ufuncs
            assert_equal(np.isnan(self.all_f16), np.isnan(self.all_f32))
            assert_equal(np.isinf(self.all_f16), np.isinf(self.all_f32))
            assert_equal(np.isfinite(self.all_f16), np.isfinite(self.all_f32))
            assert_equal(np.signbit(self.all_f16), np.signbit(self.all_f32))
            assert_equal(np.spacing(float16(65504)), np.inf)

            # Check comparisons of all values with NaN
            nan = float16(np.nan)

            assert_(not (self.all_f16 == nan).any())
            assert_(not (nan == self.all_f16).any())

            assert_((self.all_f16 != nan).all())
            assert_((nan != self.all_f16).all())

            assert_(not (self.all_f16 < nan).any())
            assert_(not (nan < self.all_f16).any())

            assert_(not (self.all_f16 <= nan).any())
            assert_(not (nan <= self.all_f16).any())

            assert_(not (self.all_f16 > nan).any())
            assert_(not (nan > self.all_f16).any())

            assert_(not (self.all_f16 >= nan).any())
            assert_(not (nan >= self.all_f16).any())
        finally:
            np.seterr(**oldsettings)

def lscsum0(lx):
  """
  Accepts log-values as input, exponentiates them, sums down the rows
  (first dimension), then converts the sum back to log-space and returns the result.
  Handles underflow by rescaling so that the largest values is exactly 1.0.
  """
  # rows = lx.shape[0]
  # columns = numpy.prod(lx.shape[1:])
  # lx = lx.reshape(rows, columns)
  # bases = lx.max(1).reshape(rows, 1)
  # bases = lx.max(0).reshape((1,) + lx.shape[1:])
  lx = numpy.asarray(lx)
  bases = lx.max(0) # Don't need to reshape in the case of 0.
  x = numpy.exp(lx - bases)
  ssum = x.sum(0)

  result = numpy.log(ssum) + bases
  try:
    conventional = numpy.log(numpy.exp(lx).sum(0))

    if not similar(result, conventional):
      if numpy.isinf(conventional).any() and not numpy.isinf(result).any():
        # print "Scaled log sum down axis 0 avoided underflow or overflow."
        pass
      else:
        import sys
        print >>sys.stderr, "Warning: scaled log sum down axis 0 did not match."
        print >>sys.stderr, "Scaled log result:"
        print >>sys.stderr, result
        print >>sys.stderr, "Conventional result:"
        print >>sys.stderr, conventional
  except FloatingPointError, e:
    # print "Scaled log sum down axis 0 avoided underflow or overflow."
    pass

def get_region_boxes(sp, reg2sp):
  x = np.arange(0, sp.shape[1])
  y = np.arange(0, sp.shape[0])
  xv, yv = np.meshgrid(x, y)
  maxsp = np.max(sp)
  sp1=sp.reshape(-1)-1
  xv = xv.reshape(-1)
  yv = yv.reshape(-1)
  spxmin = accum.my_accumarray(sp1,xv, maxsp, 'min')
  spymin = accum.my_accumarray(sp1,yv, maxsp, 'min')
  spxmax = accum.my_accumarray(sp1,xv, maxsp, 'max')
  spymax = accum.my_accumarray(sp1,yv, maxsp, 'max')
  
  Z = reg2sp.astype(float, copy=True)
  Z[reg2sp==0] = np.inf
  xmin = np.nanmin(np.multiply(spxmin.reshape(-1,1), Z),0)
  ymin = np.nanmin(np.multiply(spymin.reshape(-1,1), Z),0)
  xmax = np.amax(np.multiply(spxmax.reshape(-1,1), reg2sp),0)
  ymax = np.amax(np.multiply(spymax.reshape(-1,1), reg2sp), 0)
  xmin[np.isinf(xmin)]=0
  ymin[np.isinf(ymin)]=0
  

  boxes = np.hstack((xmin.reshape(-1,1), ymin.reshape(-1,1), xmax.reshape(-1,1), ymax.reshape(-1,1)))
  return boxes 

def find_reasonable_epsilon(theta0, grad0, logp0, f):
    """ Heuristic for choosing an initial value of epsilon """
    epsilon = 1.
    r0 = np.random.normal(0., 1., len(theta0))

    # Figure out what direction we should be moving epsilon.
    _, rprime, gradprime, logpprime = leapfrog(theta0, r0, grad0, epsilon, f)
    # brutal! This trick make sure the step is not huge leading to infinite
    # values of the likelihood. This could also help to make sure theta stays
    # within the prior domain (if any)
    k = 1.
    while np.isinf(logpprime) or np.isinf(gradprime).any():
        k *= 0.5
        _, rprime, _, logpprime = leapfrog(theta0, r0, grad0, epsilon * k, f)

    epsilon = 0.5 * k * epsilon

    # acceptprob = np.exp(logpprime - logp0 - 0.5 * (np.dot(rprime, rprime.T) - np.dot(r0, r0.T)))
    # a = 2. * float((acceptprob > 0.5)) - 1.
    logacceptprob = logpprime-logp0-0.5*(np.dot(rprime, rprime)-np.dot(r0,r0))
    a = 1. if logacceptprob > np.log(0.5) else -1.
    # Keep moving epsilon in that direction until acceptprob crosses 0.5.
    # while ( (acceptprob ** a) > (2. ** (-a))):
    while a * logacceptprob > -a * np.log(2):
        epsilon = epsilon * (2. ** a)
        _, rprime, _, logpprime = leapfrog(theta0, r0, grad0, epsilon, f)
        # acceptprob = np.exp(logpprime - logp0 - 0.5 * ( np.dot(rprime, rprime.T) - np.dot(r0, r0.T)))
        logacceptprob = logpprime-logp0-0.5*(np.dot(rprime, rprime)-np.dot(r0,r0))

    print("find_reasonable_epsilon=", epsilon)

    return epsilon

    def _set_spinbox_limits(self, bottom_val, top_val):
        # turn off signals on the spin boxes
        reset_state = [(sb, sb.blockSignals(True)) for sb in
                       (self._spin_max,
                        self._spin_min)]
        try:
            # set the top and bottom limits on the spinboxs to be in bounds
            self._spin_max.setMinimum(bottom_val)
            self._spin_min.setMinimum(bottom_val)

            self._spin_max.setMaximum(top_val)
            self._spin_min.setMaximum(top_val)
            # don't let the step be bigger than the total allowed range
            self._spin_step.setMaximum(top_val - bottom_val)

            if not np.isinf(bottom_val) or not np.isinf(top_val):
                # set the current values
                self._spin_min.setValue(bottom_val)
                self._spin_max.setValue(top_val)

                # this will trigger via the call-back updating everything else
                self._spin_step.setValue(
                    (top_val - bottom_val) / 100)
        finally:
            # un-wrap the signal blocking
            [sb.blockSignals(state) for sb, state in reset_state]

def knn(x_train, y_train, x_valid):
    x_train=np.log(x_train+1)
    x_valid=np.log(x_valid+1)

    where_are_nan = np.isnan(x_train)
    where_are_inf = np.isinf(x_train)
    x_train[where_are_nan] = 0
    x_train[where_are_inf] = 0
    where_are_nan = np.isnan(x_valid)
    where_are_inf = np.isinf(x_valid)
    x_valid[where_are_nan] = 0
    x_valid[where_are_inf] = 0

    scale=StandardScaler()
    scale.fit(x_train)
    x_train=scale.transform(x_train)
    x_valid=scale.transform(x_valid)

    #pca = PCA(n_components=10)
    #pca.fit(x_train)
    #x_train = pca.transform(x_train)
    #x_valid = pca.transform(x_valid)

    kneighbors=KNeighborsClassifier(n_neighbors=200,n_jobs=-1)
    knn_train, knn_test = stacking(kneighbors, x_train, y_train, x_valid, "knn")
    return knn_train, knn_test, "knn"

 def __init__(self, pt1, pt2, imageSize=None):
   if pt1[0] <= pt2[0]:  # ensure pt1 is to the left of pt2 for easier computations later on
     self.pt1 = pt1
     self.pt2 = pt2
   else:
     self.pt1 = pt2
     self.pt2 = pt1
   self.delta = np.subtract(self.pt2, self.pt1)
   self.length = sqrt(self.delta[0]**2 + self.delta[1]**2)
   self.m = float(self.delta[1]) / float(self.delta[0]) if self.delta[0] != 0.0 else (np.inf if self.delta[1] >=0 else -np.inf)
   self.c = self.pt1[1] - self.m * self.pt1[0]
   #print "delta = {0}, m = {1}, c = {2}".format(self.delta, self.m, self.c)
   
   # Check for validity/stability
   if np.isinf(self.m) or np.isinf(self.c):
     self.angle = 0.0
     self.valid = False
     return
   
   # Compute angle in degrees
   self.angle = degrees(atan2(self.delta[1], self.delta[0]))
   self.valid = True
   
   # Compute points on left and right edges, if an imageSize is given
   if imageSize is None:
     self.ptLeft = self.pt1
     self.ptRight = self.pt2
   else:
     self.ptLeft = (0, int(self.c))
     self.ptRight = (imageSize[0] - 1, int(self.m * (imageSize[0] - 1) + self.c))

def equal(a, b, exact):
    if array_equal(a, b):
        return True

    if hasattr(a, 'dtype') and a.dtype in ['f4','f8']:
        nnans = isnan(a).sum()
        if nnans > 0:
            # For results containing NaNs, just check that the number
            # of NaNs is the same in both arrays.  This check could be
            # made more exhaustive, but checking element by element in
            # python space is very expensive in general.
            return nnans == isnan(b).sum()
        ninfs = isinf(a).sum()
        if ninfs > 0:
            # Ditto for Inf's
            return ninfs == isinf(b).sum()
    if exact:
        return (shape(a) == shape(b)) and alltrue(ravel(a) == ravel(b), axis=0)
    else:
        if hasattr(a, 'dtype') and a.dtype == 'f4':
            atol = 1e-5   # Relax precission for special opcodes, like fmod
        else:
            atol = 1e-8
        return (shape(a) == shape(b) and
                allclose(ravel(a), ravel(b), atol=atol))