def inverse_magnification_tensors(self, imagepositions):
     tiny = 1e-12 # to deal with lens and image both exactly at origin
     ipos = np.atleast_2d(imagepositions)
     mag = np.zeros((len(ipos), 2, 2))
     mag[:,0,0] = 1.
     mag[:,1,1] = 1.
     # print "inverse_magnification_tensors: ipos = ",ipos
     dpos = ipos - self.position
     rcubed = np.sum(dpos * dpos, axis=1)**1.5 + tiny # tiny little hack
     if np.min(rcubed) <= 0.0:
         print "image positions = ",ipos
         print "lens position = ",self.position
         print "differences = ",dpos
         print "rcubed = ",rcubed
         print self
     else:
         mag[:,0,0] -= self.einsteinradius * dpos[:,1] * dpos[:,1] / rcubed
         mag[:,0,1] += self.einsteinradius * dpos[:,1] * dpos[:,0] / rcubed
         mag[:,1,0] += self.einsteinradius * dpos[:,0] * dpos[:,1] / rcubed
         mag[:,1,1] -= self.einsteinradius * dpos[:,0] * dpos[:,0] / rcubed
         mag[:,0,0] -= self.gammacos2phi
         mag[:,0,1] -= self.gammasin2phi
         mag[:,1,0] -= self.gammasin2phi
         mag[:,1,1] += self.gammacos2phi
     assert(np.min(rcubed) > 0.0)
     return mag

def cfl_superbee_theta(r,cfl,theta=0.95):
    r"""
    CFL-Superbee (Roe's Ultrabee) with theta parameter
    """
    a = np.empty((2,len(r)))
    b = np.zeros((2,len(r)))
    
    a[0,:] = 0.001
    a[1,:] = cfl
    cfmod1 = np.max(a,axis=0)
    a[0,:] = 0.999
    cfmod2 = np.min(a,axis=0)

    s1 = theta * 2.0 / cfmod1
    phimax = theta * 2.0 / (1.0 - cfmod2)

    a[0,:] = s1*r
    a[1,:] = phimax
    b[1,:] = np.min(a,axis=0)
    ultra = np.max(b,axis=0)
    
    a[0,:] = ultra
    b[0,:] = 1.0
    b[1,:] = r
    a[1,:] = np.max(b,axis=0)
    return np.min(a,axis=0)

def cada_torrilhon_limiter(r,cfl,epsilon=1.0e-3):
    r"""
    Cada-Torrilhon modified
    
    Additional Input:
     - *epsilon* = 
    """
    a = np.ones((2,len(r))) * 0.95
    b = np.empty((3,len(r)))

    a[0,:] = cfl
    cfl = np.min(a)
    a[1,:] = 0.05
    cfl = np.max(a)
    
    # Multiply all parts except b[0,:] by (1.0 - epsilon) as well
    b[0,:] = 1.0 + (1+cfl) / 3.0 * (r - 1)
    b[1,:] = 2.0 * np.abs(r) / (cfl + epsilon)
    b[2,:] = (8.0 - 2.0 * cfl) / (np.abs(r) * (cfl - 1.0 - epsilon)**2)
    b[1,::2] *= (1.0 - epsilon)
    a[0,:] = np.min(b)
    a[1,:] = (-2.0 * (cfl**2 - 3.0 * cfl + 8.0) * (1.0-epsilon)
                    / (np.abs(r) * (cfl**3 - cfl**2 - cfl + 1.0 + epsilon)))
    
    return np.max(a)

    def __init__(self, shape, successes, 
                 trials=None, coef=1., offset=None,
                 quadratic=None,
                 initial=None):

        smooth_atom.__init__(self,
                             shape,
                             offset=offset,
                             quadratic=quadratic,
                             initial=initial,
                             coef=coef)

        if sparse.issparse(successes):
            #Convert sparse success vector to an array
            self.successes = successes.toarray().flatten()
        else:
            self.successes = np.asarray(successes)

        if trials is None:
            if not set([0,1]).issuperset(np.unique(self.successes)):
                raise ValueError("Number of successes is not binary - must specify number of trials")
            self.trials = np.ones(self.successes.shape, np.float)
        else:
            if np.min(trials-self.successes) < 0:
                raise ValueError("Number of successes greater than number of trials")
            if np.min(self.successes) < 0:
                raise ValueError("Response coded as negative number - should be non-negative number of successes")
            self.trials = trials * 1.

        saturated = self.successes / self.trials
        deviance_terms = np.log(saturated) * self.successes + np.log(1-saturated) * (self.trials - self.successes)
        deviance_constant = -2 * coef * deviance_terms[~np.isnan(deviance_terms)].sum()

        devq = identity_quadratic(0,0,0,-deviance_constant)
        self.quadratic += devq

    def _crinfo_from_specific_data (self, data, margin):
# hledáme automatický ořez, nonzero dá indexy
        nzi = np.nonzero(data)

        x1 = np.min(nzi[0]) - margin[0]
        x2 = np.max(nzi[0]) + margin[0] + 1
        y1 = np.min(nzi[1]) - margin[0]
        y2 = np.max(nzi[1]) + margin[0] + 1
        z1 = np.min(nzi[2]) - margin[0]
        z2 = np.max(nzi[2]) + margin[0] + 1 

# ošetření mezí polí
        if x1 < 0:
            x1 = 0
        if y1 < 0:
            y1 = 0
        if z1 < 0:
            z1 = 0

        if x2 > data.shape[0]:
            x2 = data.shape[0]-1
        if y2 > data.shape[1]:
            y2 = data.shape[1]-1
        if z2 > data.shape[2]:
            z2 = data.shape[2]-1

# ořez
        crinfo = [[x1, x2],[y1,y2],[z1,z2]]
        #dataout = self._crop(data,crinfo)
        #dataout = data[x1:x2, y1:y2, z1:z2]
        return crinfo

  def testEncodeAdjacentPositions(self, verbose=False):
    repetitions = 100
    n = 999
    w = 25
    radius = 10
    minThreshold = 0.75
    avgThreshold = 0.90
    allOverlaps = np.empty(repetitions)

    for i in range(repetitions):
      overlaps = overlapsForRelativeAreas(n, w,
                                          np.array([i * 10, i * 10]), radius,
                                          dPosition=np.array([0, 1]),
                                          num=1)
      allOverlaps[i] = overlaps[0]

    self.assertGreater(np.min(allOverlaps), minThreshold)
    self.assertGreater(np.average(allOverlaps), avgThreshold)

    if verbose:
      print ("===== Adjacent positions overlap "
             "(n = {0}, w = {1}, radius = {2}) ===").format(n, w, radius)
      print "Max: {0}".format(np.max(allOverlaps))
      print "Min: {0}".format(np.min(allOverlaps))
      print "Average: {0}".format(np.average(allOverlaps))

def grid_xyz(xyz, n_x, n_y, **kwargs):
    """ Grid data as a list of X,Y,Z coords into a 2D array

    Parameters
    ----------
    xyz: np.array
        Numpy array of X,Y,Z values, with shape (n_points, 3)
    n_x: int
        Number of points in x direction (fastest varying!)
    n_y: int
        Number of points in y direction

    Returns
    -------
    gridded_data: np.array
        2D array of gridded data, with shape (n_x, n_y)

    Notes
    -----
    'x' is the inner dimension, i.e. image dimensions are (n_y, n_x). This is
    counterintuitive (to me at least) but in line with numpy definitions.
    """
    x, y, z = xyz[:, 0], xyz[:, 1], xyz[:, 2]
    x_ax = np.linspace(np.min(x), np.max(x), n_x)
    y_ax = np.linspace(np.min(y), np.max(y), n_y)

    xg, yg = np.meshgrid(x_ax, y_ax)

    data = griddata(xyz[:, :2], z, (xg, yg), **kwargs)
    return data    

    def function1D(self, t):
        A  = self.getParamValue(0)
        B  = self.getParamValue(1)
        R  = self.getParamValue(2)
        T0 = self.getParamValue(3)
        Scale = self.getParamValue(4)
        HatWidth  = self.getParamValue(5)
        KConv  = self.getParamValue(6)

        # A/2 Scale factor has been removed to make A and Scale independent
        f_int = Scale*((1-R)*np.power((A*(t-T0)),2)*
                       np.exp(-A*(t-T0))+2*R*A**2*B/np.power((A-B),3) *
                       (np.exp(-B*(t-T0))-np.exp(-A*(t-T0))*(1+(A-B)*(t-T0)+0.5*np.power((A-B),2)*np.power((t-T0),2))))
        f_int[t<T0] = 0

        mid_point_hat = len(f_int)//2
        gc_x = np.array(range(len(f_int))).astype(float)
        ppd = 0.0*gc_x
        lowIDX  = int(np.floor(np.max([mid_point_hat-np.abs(HatWidth),0])))
        highIDX = int(np.ceil(np.min([mid_point_hat+np.abs(HatWidth),len(gc_x)])))

        ppd[lowIDX:highIDX] = 1.0
        ppd = ppd/sum(ppd)

        gc_x = np.array(range(len(f_int))).astype(float)
        gc_x = 2*(gc_x-np.min(gc_x))/(np.max(gc_x)-np.min(gc_x))-1
        gc_f = np.exp(-KConv*np.power(gc_x,2))
        gc_f = gc_f/np.sum(gc_f)

        npad = len(f_int) - 1
        first = npad - npad//2
        f_int = np.convolve(f_int,ppd,'full')[first:first+len(f_int)]
        f_int = np.convolve(f_int,gc_f,'full')[first:first+len(f_int)]

        return f_int

def signmag_plot(a, b, z, ref):
    imdata1 = np.sign(ref)
    cmap1 = plt.cm.RdBu
    cmap1.set_bad('k', 1)

    imdata2 = np.log10(np.abs(ref))
    cmap2 = plt.cm.YlOrRd
    cmap2.set_bad('k', 1)

    fig, axarr = plt.subplots(ncols=2, figsize=(12, 6))
    axarr[0].pcolormesh(a, b, imdata1, cmap=cmap1, vmin=-1, vmax=1)
    im = axarr[1].pcolormesh(a, b, imdata2, cmap=cmap2,
                             vmin=np.percentile(imdata2,  5),
                             vmax=np.percentile(imdata2, 95))

    for ax in axarr:
        ax.set_xlim((np.min(a), np.max(a)))
        ax.set_ylim((np.min(b), np.max(b)))
        ax.set_xlabel("a")
        ax.set_ylabel("b")
        ax.set(adjustable='box-forced', aspect='equal')

    fig.subplots_adjust(right=0.8)
    cbar_ax = fig.add_axes([0.85, 0.15, 0.03, 0.7])
    fig.colorbar(im, cax=cbar_ax)

    axarr[0].set_title("Sign of hyp1f1")
    axarr[1].set_title("Magnitude of hyp1f1")
    plt.suptitle("z = {:.2e}".format(np.float64(z)))

    return fig

def explore_city_data(city_data):
    """Calculate the Boston housing statistics."""

    # Get the labels and features from the housing data
    housing_prices = city_data.target
    housing_features = city_data.data

    ###################################
    ### Step 1. YOUR CODE GOES HERE ###
    ###################################

    # Please calculate the following values using the Numpy library
    print "Size of data (number of houses)"
    print np.size(housing_prices)
    print "Number of features"
    print np.size(housing_features, 1)
    print "Minimum price"
    print np.min(housing_prices)
    print "Maximum price"
    print np.max(housing_prices)
    print "Calculate mean price"
    print np.mean(housing_prices)
    print "Calculate median price"
    print np.median(housing_prices)
    print "Calculate standard deviation"
    print np.std(housing_prices)

 def coverage_string(self):
     """Coverage of reader to be reported as string for debug output"""
     corners = self.xy2lonlat([self.xmin, self.xmin, self.xmax, self.xmax],
                              [self.ymax, self.ymin, self.ymax, self.ymin])
     return '%.2f-%.2fE, %.2f-%.2fN' % (
                 np.min(corners[0]), np.max(corners[0]),
                 np.min(corners[1]), np.max(corners[1]))

def hausdorffnorm(A, B):
    '''
    Finds the hausdorff norm between two matrices A and B.
    INPUTS:
    A: numpy array
    B : numpy array
    OUTPUTS:
    Housdorff norm between matrices A and B
    '''
    # ensure matrices are 3 dimensional, and shaped conformably
    if len(A.shape) == 1:
        A = np.atleast_2d(A)

    if len(B.shape) == 1:
        B = np.atleast_2d(B)

    A = np.atleast_3d(A)
    B = np.atleast_3d(B)

    x, y, z = B.shape
    A = np.reshape(A, (z, x, y))
    B = np.reshape(B, (z, x, y))

    # find hausdorff norm: starting from A to B
    z, x, y = B.shape
    temp1 = np.tile(np.reshape(B.T, (y, z, x)), (max(A.shape), 1))
    temp2 = np.tile(np.reshape(A.T, (y, x, z)), (1, max(B.shape)))
    D1 = np.min(np.sqrt(np.sum((temp1-temp2)**2, 0)), axis=0)

    # starting from B to A
    temp1 = np.tile(np.reshape(A.T, (y, z, x)), (max(B.shape), 1))
    temp2 = np.tile(np.reshape(B.T, (y, x, z)), (1, max(A.shape)))
    D2 = np.min(np.sqrt(np.sum((temp1-temp2)**2, 0)), axis=0)

    return np.max([D1, D2])

def rootSpI(img, list_remove=[], sc=None, lut_range = False, verbose=False):
    """
    case where the data is a spatialimage
    """
    # -- cells are positionned inside a structure, the polydata, and assigned a scalar value.
    polydata,polydata2 = img2polydata_complexe(img, list_remove=list_remove, sc=sc, verbose=verbose)
    m = tvtk.PolyDataMapper(input=polydata.output)
    m2 = tvtk.PolyDataMapper(input=polydata2.output)
    
    # -- definition of the scalar range (default : min to max of the scalar value).
    if sc:
        ran=[sc[i] for i in sc.keys() if i not in list_remove]
        if (lut_range != None) and (lut_range != False):
            print lut_range
            m.scalar_range = lut_range[0],lut_range[1]
        else:
            m.scalar_range = np.min(ran), np.max(ran)
    else:
        m.scalar_range=np.min(img), np.max(img)
    
    # -- actor that manage changes of view if memory is short.
    a = tvtk.QuadricLODActor(mapper=m)
    a.property.point_size=8
    a2 = tvtk.QuadricLODActor(mapper=m2)
    a2.property.point_size=8
    #scalebar
    if lut_range != None:
        sc=tvtk.ScalarBarActor(orientation='vertical',lookup_table=m.lookup_table)
    return a, a2, sc, m, m2

    def _makewindows(self, indices, window):
        """
        Make masks used by windowing functions

        Given a list of indices specifying window centers,
        and a window size, construct a list of index arrays,
        one per window, that index into the target array

        Parameters
        ----------
        indices : array-like
            List of times specifying window centers

        window : int
            Window size
        """
        div = divmod(window, 2)
        before = div[0]
        after = div[0] + div[1]
        index = asarray(self.index)
        indices = asarray(indices)
        if where(index == max(indices))[0][0] + after > len(index):
            raise ValueError("Maximum requested index %g, with window %g, exceeds length %g"
                             % (max(indices), window, len(index)))
        if where(index == min(indices))[0][0] - before < 0:
            raise ValueError("Minimum requested index %g, with window %g, is less than 0"
                             % (min(indices), window))
        masks = [arange(where(index == i)[0][0]-before, where(index == i)[0][0]+after, dtype='int') for i in indices]
        return masks

    def quantify(self):
        """Quantify shape of the contours."""
        four_pi = 4. * np.pi
        for edge in self.edges:
            # Positions
            x = edge['x']
            y = edge['y']

            A, perimeter, x_center, y_center, distances = \
                self.get_shape_factor(x, y)

            # Set values.
            edge['area'] = A
            edge['perimeter'] = perimeter
            edge['x_center'] = x_center
            edge['y_center'] = y_center
            # Circle is 1. Rectangle is 0.78. Thread-like is close to zero.
            edge['shape_factor'] = four_pi * edge['area'] / \
                                   edge['perimeter'] ** 2.

            # We assume that the radius of the edge
            # as the median value of the distances from the center.
            radius = np.median(distances)
            edge['radius_deviation'] = np.std(distances - radius) / radius

            edge['x_min'] = np.min(x)
            edge['x_max'] = np.max(x)
            edge['y_min'] = np.min(y)
            edge['y_max'] = np.max(y)

def allclose_with_out(x, y, atol=0.0, rtol=1.0e-5):
    # run the np.allclose on x and y
    # if it fails print some stats
    # before returning
    ac = np.allclose(x, y, rtol=rtol, atol=atol)
    if not ac:
        dd = np.abs(x - y)
        neon_logger.display('abs errors: %e [%e, %e] Abs Thresh = %e'
                            % (np.median(dd), np.min(dd), np.max(dd), atol))
        amax = np.argmax(dd)

        if np.isscalar(x):
            neon_logger.display('worst case: %e %e' % (x, y.flat[amax]))
        elif np.isscalar(y):
            neon_logger.display('worst case: %e %e' % (x.flat[amax], y))
        else:
            neon_logger.display('worst case: %e %e' % (x.flat[amax], y.flat[amax]))

        dd = np.abs(dd - atol) / np.abs(y)
        neon_logger.display('rel errors: %e [%e, %e] Rel Thresh = %e'
                            % (np.median(dd), np.min(dd), np.max(dd), rtol))
        amax = np.argmax(dd)
        if np.isscalar(x):
            neon_logger.display('worst case: %e %e' % (x, y.flat[amax]))
        elif np.isscalar(y):
            neon_logger.display('worst case: %e %e' % (x.flat[amax], y))
        else:
            neon_logger.display('worst case: %e %e' % (x.flat[amax], y.flat[amax]))
    return ac

 def min(self, axis=None, out=None, keepdims=False):
     self._prepare_out(out=out)
     try:
         value = np.min(self.value, axis=axis, out=out, keepdims=keepdims)
     except:  # numpy < 1.7
         value = np.min(self.value, axis=axis, out=out)
     return self._new_view(value)

def collect_statistics_from_sigma_bins(sigma, bins_start, bins_end, burnin=0, smooth=True, area_fraction=0.68, numpoints=1000):
    sigma = flatten_commander_chain(sigma, burnin)
    lmax = sigma.shape[2] - 1
    means = []
    stds = []
    mls = []
    uppers = []
    lowers = []
    for lstart, lend in zip(bins_start, bins_end):
        vars = []
        sigmas = []
        for l in range(lstart, lend+1):
            print l
            vars.append(np.var(sigma[:, 0, l]))
            print vars[-1]
            sigmas.append(sigma[:, 0, l] / vars[-1])
        vars = np.array(vars)
        sigmas = np.array(sigmas)
        if lstart == lend:
            samps = sigmas
            samps = samps * vars
        else:
            samps = np.sum(sigmas, axis=0)
            samps = samps / np.sum(1 / vars)
        print np.min(samps)
        print np.max(samps)
        x = np.linspace(np.min(samps), np.max(samps), numpoints)
        mean, std, ml, upper, lower = collect_statistics(samps, x, area_fraction=area_fraction, smooth=smooth)
        means.append(mean)
        stds.append(std)
        mls.append(ml)
        uppers.append(upper)
        lowers.append(lower)
    return means, stds, mls, uppers, lowers

def plot_embedding(X, title=None):
    x_min, x_max = np.min(X, 0), np.max(X, 0)
    X = (X - x_min) / (x_max - x_min)

    pl.figure()
    ax = pl.subplot(111)
    for i in range(digits.data.shape[0]):
        pl.text(
            X[i, 0],
            X[i, 1],
            str(digits.target[i]),
            color=pl.cm.Set1(digits.target[i] / 10.0),
            fontdict={"weight": "bold", "size": 9},
        )

    if hasattr(offsetbox, "AnnotationBbox"):
        # only print thumbnails with matplotlib > 1.0
        shown_images = np.array([[1.0, 1.0]])  # just something big
        for i in range(digits.data.shape[0]):
            dist = np.sum((X[i] - shown_images) ** 2, 1)
            if np.min(dist) < 4e-3:
                # don't show points that are too close
                continue
            shown_images = np.r_[shown_images, [X[i]]]
            imagebox = offsetbox.AnnotationBbox(offsetbox.OffsetImage(digits.images[i], cmap=pl.cm.gray_r), X[i])
            ax.add_artist(imagebox)
    pl.xticks([]), pl.yticks([])
    if title is not None:
        pl.title(title)

def check_min_samples_split(name):
    X, y = hastie_X, hastie_y
    ForestEstimator = FOREST_ESTIMATORS[name]

    # test boundary value
    assert_raises(ValueError,
                  ForestEstimator(min_samples_split=-1).fit, X, y)
    assert_raises(ValueError,
                  ForestEstimator(min_samples_split=0).fit, X, y)
    assert_raises(ValueError,
                  ForestEstimator(min_samples_split=1.1).fit, X, y)

    est = ForestEstimator(min_samples_split=10, n_estimators=1, random_state=0)
    est.fit(X, y)
    node_idx = est.estimators_[0].tree_.children_left != -1
    node_samples = est.estimators_[0].tree_.n_node_samples[node_idx]

    assert_greater(np.min(node_samples), len(X) * 0.5 - 1,
                   "Failed with {0}".format(name))

    est = ForestEstimator(min_samples_split=0.5, n_estimators=1, random_state=0)
    est.fit(X, y)
    node_idx = est.estimators_[0].tree_.children_left != -1
    node_samples = est.estimators_[0].tree_.n_node_samples[node_idx]

    assert_greater(np.min(node_samples), len(X) * 0.5 - 1,
                   "Failed with {0}".format(name))