Python pylab 模块，axis() 实例源码

def disp(iimg, label = "", gray=False):
    """ Display an image using pylab
    """
    try:
        import pylab
        dimage = iimg.copy()
        if iimg.ndim==3:
            dimage[...,0] = iimg[...,2]
            dimage[...,2] = iimg[...,0]

        pylab.imshow(dimage, interpolation='none')
        if gray: pylab.gray()
        #pylab.gca().format_coord = format_coord
        pylab.text(1500, -30, label)
        pylab.axis('off')
        pylab.show()
    except ImportError:
        print "Module pylab not available"

def edgescatter(self, ps):
        for ei,X in enumerate(self.edges):
            i,j = X[:2]
            matchdRA, matchdDec = X[10:12]
            mu = X[9]
            A = self.alignments[ei]

            plt.clf()
            if len(matchdRA) > 1000:
                plothist(matchdRA, matchdDec, 101)
            else:
                plt.plot(matchdRA, matchdDec, 'k.', alpha=0.5)
            plt.axvline(0, color='0.5')
            plt.axhline(0, color='0.5')
            plt.axvline(mu[0], color='b')
            plt.axhline(mu[1], color='b')
            for nsig in [1,2]:
                X,Y = A.getContours(nsigma=nsig)
                plt.plot(X, Y, 'b-')
            plt.xlabel('delta-RA (arcsec)')
            plt.ylabel('delta-Dec (arcsec)')
            plt.axis('scaled')
            ps.savefig()

def plotaffine(aff, RR, DD, exag=1000, affineOnly=False, doclf=True, **kwargs):
    import pylab as plt
    if doclf:
        plt.clf()
    if affineOnly:
        dr,dd = aff.getAffineOffset(RR, DD)
    else:
        rr,dd = aff.apply(RR, DD)
        dr = rr - RR
        dd = dd - DD
    #plt.plot(RR, DD, 'r.')
    #plt.plot(RR + dr*exag, DD + dd*exag, 'bx')
    plt.quiver(RR, DD, exag*dr, exag*dd,
               angles='xy', scale_units='xy', scale=1,
               pivot='middle', color='b', **kwargs)
               #pivot='tail'
    ax = plt.axis()
    plt.plot([aff.getReferenceRa()], [aff.getReferenceDec()], 'r+', mew=2, ms=5)
    plt.axis(ax)
    esuf = ''
    if exag != 1.:
        esuf = ' (x %g)' % exag
    plt.title('Affine transformation found' + esuf)

def mahalanobis_distsq(X, mu, C):
    #print 'X', X.shape
    #print 'mu', mu.shape
    #print 'C', C.shape
    assert(C.shape == (2,2))
    det = C[0,0]*C[1,1] - C[0,1]*C[1,0]
    Cinv = np.array([[ C[1,1]/det, -C[0,1]/det],
             [-C[1,0]/det,  C[0,0]/det]])
    #print 'Cinv', Cinv.shape
    #print 'Cinv * C:', np.dot(Cinv, C)
    #print 'C * Cinv:', np.dot(C, Cinv)
    d = X - mu
    #print 'd', d.shape
    Cinvd = np.dot(Cinv, d.T)
    #print 'Cinvd', Cinvd.shape
    M = np.sum(d * Cinvd.T, axis=1)
    #print 'M', M.shape
    return M

def plot_rectified(self):
        import pylab
        pylab.title('rectified')
        pylab.imshow(self.rectified)

        for line in self.vlines:
            p0, p1 = line
            p0 = self.inv_transform(p0)
            p1 = self.inv_transform(p1)
            pylab.plot((p0[0], p1[0]), (p0[1], p1[1]), c='green')

        for line in self.hlines:
            p0, p1 = line
            p0 = self.inv_transform(p0)
            p1 = self.inv_transform(p1)
            pylab.plot((p0[0], p1[0]), (p0[1], p1[1]), c='red')

        pylab.axis('image');
        pylab.grid(c='yellow', lw=1)
        pylab.plt.yticks(np.arange(0, self.l, 100.0));
        pylab.xlim(0, self.w)
        pylab.ylim(self.l, 0)

def plot_original(self):
        import pylab
        pylab.title('original')
        pylab.imshow(self.data)

        for line in self.lines:
            p0, p1 = line
            pylab.plot((p0[0], p1[0]), (p0[1], p1[1]), c='blue', alpha=0.3)

        for line in self.vlines:
            p0, p1 = line
            pylab.plot((p0[0], p1[0]), (p0[1], p1[1]), c='green')

        for line in self.hlines:
            p0, p1 = line
            pylab.plot((p0[0], p1[0]), (p0[1], p1[1]), c='red')

        pylab.axis('image');
        pylab.grid(c='yellow', lw=1)
        pylab.plt.yticks(np.arange(0, self.l, 100.0));
        pylab.xlim(0, self.w)
        pylab.ylim(self.l, 0)

def _cut_windows_vertically(self, door_top, roof_top, sky_sig, win_strip):
        win_sig = np.percentile(win_strip, 85, axis=1)
        win_sig[sky_sig > 0.5] = 0
        if win_sig.max() > 0:
            win_sig /= win_sig.max()
        win_sig[:roof_top] = 0
        win_sig[door_top:] = 0
        runs, starts, values = run_length_encode(win_sig > 0.5)
        win_heights = runs[values]
        win_tops = starts[values]
        if len(win_heights) > 0:
            win_bottom = win_tops[-1] + win_heights[-1]
            win_top = win_tops[0]
            win_vertical_spacing = np.diff(win_tops).mean() if len(win_tops) > 1 else 0
        else:
            win_bottom = win_top = win_vertical_spacing = -1

        self.top = int(win_top)
        self.bottom = int(win_bottom)
        self.vertical_spacing = int(win_vertical_spacing)
        self.vertical_scores = make_list(win_sig)
        self.heights = np.array(win_heights)
        self.tops = np.array(win_tops)

def ransac_guess_color(colors, n_iter=50, std=2):
    colors = rgb2lab(colors)
    colors = colors.reshape(-1, 3)
    masked = colors[:, 0] < 0.1
    colors = colors[~masked]
    assert len(colors) > 0, "Must have at least one color"

    best_mu = np.array([0, 0, 0])
    best_n = 0
    for k in range(n_iter):
        subset = colors[np.random.choice(np.arange(len(colors)), 1)]

        mu = subset.mean(0)
        #inliers = (((colors - mu) ** 2 / std) < 1).all(1)
        inliers = ((np.sqrt(np.sum((colors - mu)**2, axis=1))  / std) < 1)

        mu = colors[inliers].mean(0)
        n = len(colors[inliers])
        if n > best_n:
            best_n = n
            best_mu = mu
    #import ipdb; ipdb.set_trace()
    best_mu = np.squeeze(lab2rgb(np.array([[best_mu]])))
    return best_mu

def joint_density(X, Y, bounds=None):
    """
Plots joint distribution of variables.
Inherited from method in src/graphics.py module in project 
git://github.com/aflaxman/pymc-example-tfr-hdi.git
    """
    if bounds:
        X_min, X_max, Y_min, Y_max = bounds
    else:
        X_min = X.min()
        X_max = X.max()
        Y_min = Y.min()
        Y_max = Y.max()

    pylab.plot(X, Y, linestyle='none', marker='o', color='green', mec='green', alpha=.2, zorder=-99)

    gkde = scipy.stats.gaussian_kde([X, Y])
    x,y = pylab.mgrid[X_min:X_max:(X_max-X_min)/25.,Y_min:Y_max:(Y_max-Y_min)/25.]
    z = pylab.array(gkde.evaluate([x.flatten(), y.flatten()])).reshape(x.shape)
    pylab.contour(x, y, z, linewidths=2)

    pylab.axis([X_min, X_max, Y_min, Y_max])

def plot_profiles(self):
        """
        Plot TOPP profiles, e.g. for debugging.
        """
        import pylab
        pylab.ion()
        self.topp.WriteProfilesList()
        self.topp.WriteSwitchPointsList()
        profileslist = TOPP.TOPPpy.ProfilesFromString(
            self.topp.resprofilesliststring)
        switchpointslist = TOPP.TOPPpy.SwitchPointsFromString(
            self.topp.switchpointsliststring)
        TOPP.TOPPpy.PlotProfiles(profileslist, switchpointslist)
        TOPP.TOPPpy.PlotAlphaBeta(self.topp)
        pylab.title("%s phase profile" % type(self).__name__)
        pylab.axis([0, 1, 0, 10])

def drop_variable_before_preparation(self, df):
        # Acceptable limit of NaN in features
        limit_of_nans = 0.3*df.shape[0]
        # limit_of_nans = 800
        for feature in self.features_with_missing_values_in_dataframe(df).index:
            if df[feature].isnull().sum() > limit_of_nans:
                df = df.drop([feature], axis=1)

        # df = df.drop(['Alley'], axis=1)
        # df = df.drop(['MasVnrType'], axis=1)
        # df = df.drop(["Utilities", "LotFrontage", "Alley", "MasVnrType", "MasVnrArea", "BsmtQual",
        #               "BsmtCond", "BsmtExposure", "BsmtFinType1", "BsmtFinType2",
        #               "Electrical", "FireplaceQu", "GarageType", "GarageYrBlt",
        #               "GarageFinish", "GarageQual", "GarageCond", "PoolQC",
        #               "Fence", "MiscFeature"], axis=1)
        return df

def drop_variable(self, df):
        # if HousePrices._is_one_hot_encoder:
            # Drop all categorical feature helping columns ('Num')
            # Todo: is it defined when importing data set? _feature_names_num
            # for feature_name in HousePrices._feature_names_num:
            #     df = df.drop([feature_name], axis=1)

        # is_with_feature_agglomeration = 0
        # if is_with_feature_agglomeration:
        #     print(df.shape)
        #     df = HousePrices.feature_agglomeration(df)
        #     print(df.shape)

        # df = df.drop(['Fireplaces'], axis=1)
        df = df.drop(['Id'], axis=1)

        if not any(tuple(df.columns == 'SalePrice')):
            # All feature var names occuring in test data is assigned the public varaible df_test_all_feature_var_names.
            self.df_test_all_feature_var_names = df.columns
        return df

def outlier_identification(self, model, x_train, y_train):
        # Split the training data into an extra set of test
        x_train_split, x_test_split, y_train_split, y_test_split = train_test_split(x_train, y_train)
        print('\nOutlier shapes')
        print(np.shape(x_train_split), np.shape(x_test_split), np.shape(y_train_split), np.shape(y_test_split))
        model.fit(x_train_split, y_train_split)
        y_predicted = model.predict(x_test_split)
        residuals = np.absolute(y_predicted - y_test_split)
        rmse_pred_vs_actual = self.rmse(y_predicted, y_test_split)
        outliers_mask = residuals >= rmse_pred_vs_actual
        # outliers_mask = np.insert(np.zeros((np.shape(y_train_split)[0],), dtype=np.int), np.shape(y_train_split)[0],
        #                           outliers_mask)
        outliers_mask = np.concatenate([np.zeros((np.shape(y_train_split)[0],), dtype=bool), outliers_mask])
        not_an_outlier = outliers_mask == 0
        # Resample the training set from split, since the set was randomly split
        x_out = np.insert(x_train_split, np.shape(x_train_split)[0], x_test_split, axis=0)
        y_out = np.insert(y_train_split, np.shape(y_train_split)[0], y_test_split, axis=0)
        return x_out[not_an_outlier, ], y_out[not_an_outlier, ]

def tile_images(image_batch, image_width=28, image_height=28, image_channel=1, dir=None, filename="images"):
    if dir is None:
        raise Exception()
    try:
        os.mkdir(dir)
    except:
        pass
    fig = pylab.gcf()
    fig.set_size_inches(16.0, 16.0)
    pylab.clf()
    pylab.gray()
    for m in range(100):
        pylab.subplot(10, 10, m + 1)
        pylab.imshow(image_batch[m].reshape((image_width, image_height)), interpolation="none")
        pylab.axis("off")
    pylab.savefig("{}/{}.png".format(dir, filename))

def view_(_pred,_lable):

    fname = ['Captcha/lv3/%i.jpg' %i for i in range(20)]
    img = []
    for fn in fname:
        img.append(Image.open(open(fn)))
        #img.append(misc.imread(fn).astype(np.float))
    for i in range(len(img)):
        pylab.subplot(4,5,i+1); pylab.axis('off')

        pylab.imshow(img[i])
        #pylab.imshow( np.dot(np.array(img[i])[...,:3],[0.299,0.587,0.114]) , cmap=plt.get_cmap("gray"))
        #pylab.text(40,60,_pred[i],color = 'b')
        if ( _pred[i] == _lable[i] ):
            pylab.text(40,65,_pred[i],color = 'b',size = 15)
        else:
            pylab.text(40,65,_pred[i],color = 'r',size = 15)

        pylab.text(40,92,_lable[i],color = 'g',size = 15)

    pylab.show()

def visualize_reconstruction(xp, model, x, visualization_dir, epoch, gpu=False):
    x_variable = chainer.Variable(xp.asarray(x))
    _x = model.decode(model.encode(x_variable), test=True)
    _x.to_cpu()
    _x = _x.data

    fig = pylab.gcf()
    fig.set_size_inches(8.0, 8.0)
    pylab.clf()
    pylab.gray()
    for m in range(50):
        i = m / 10
        j = m % 10
        pylab.subplot(10, 10, 20 * i + j + 1, xticks=[], yticks=[])
        pylab.imshow(x[m].reshape((28, 28)), interpolation="none")
        pylab.subplot(10, 10, 20 * i + j + 10 + 1, xticks=[], yticks=[])
        pylab.imshow(_x[m].reshape((28, 28)), interpolation="none")
        # pylab.imshow(np.clip((_x_batch.data[m] + 1.0) / 2.0, 0.0, 1.0).reshape(
        # (config.img_channel, config.img_width, config.img_width)), interpolation="none")
        pylab.axis("off")
    pylab.savefig("{}/reconstruction_{}.png".format(visualization_dir, epoch))
    # pylab.show()

def tile_binary_images(x, dir=None, filename="x", row=10, col=10):
    if dir is None:
        raise Exception()
    try:
        os.mkdir(dir)
    except:
        pass
    fig = pylab.gcf()
    fig.set_size_inches(col * 2, row * 2)
    pylab.clf()
    pylab.gray()
    for m in range(row * col):
        pylab.subplot(row, col, m + 1)
        pylab.imshow(np.clip(x[m], 0, 1), interpolation="none")
        pylab.axis("off")
    pylab.savefig("{}/{}.png".format(dir, filename))

def get_sqe_above(self,high_p,rich_return=False):
        p = dict()
        p['g_exc'] = high_p[0] # nS
        p['sen2dec_exc'] = high_p[1] # pA/Hz
        if len(high_p)>2:
            p['frac_proj'] = high_p[2]
        else:
            p['frac_proj'] = 0.1

        win1_all = np.zeros((6,6))
        for i_m in xrange(len(self.m_plot)):
            for i_c in xrange(len(self.c_plot)):
                m = self.m_plot[i_m]
                c = self.c_plot[i_c]
                win1_all[i_m,i_c] = self.quick_Newsome_model(p,m,c)

        sqe = 0
        sqe += sum((self.p_m_m-win1_all.mean(axis=1))**2)
        sqe += sum((self.p_m_c-win1_all.mean(axis=0))**2)

        if rich_return:
            return p
        else:
            return sqe

def transform_to_2d(data, max_axis):
    """
    Projects 3d data cube along one axis using maximum intensity with
    preservation of the signs. Adapted from nilearn.
    """
    import numpy as np
    # get the shape of the array we are projecting to
    new_shape = list(data.shape)
    del new_shape[max_axis]

    # generate a 3D indexing array that points to max abs value in the
    # current projection
    a1, a2 = np.indices(new_shape)
    inds = [a1, a2]
    inds.insert(max_axis, np.abs(data).argmax(axis=max_axis))

    # take the values where the absolute value of the projection
    # is the highest
    maximum_intensity_data = data[inds]

    return np.rot90(maximum_intensity_data)

def tile_binary_images(x, dir=None, filename="x", row=10, col=10):
    if dir is None:
        raise Exception()
    try:
        os.mkdir(dir)
    except:
        pass
    fig = pylab.gcf()
    fig.set_size_inches(col * 2, row * 2)
    pylab.clf()
    pylab.gray()
    for m in range(row * col):
        pylab.subplot(row, col, m + 1)
        pylab.imshow(np.clip(x[m], 0, 1), interpolation="none")
        pylab.axis("off")
    pylab.savefig("{}/{}.png".format(dir, filename))

def tile_binary_images(x, dir=None, filename="x"):
    if dir is None:
        raise Exception()
    try:
        os.mkdir(dir)
    except:
        pass
    fig = pylab.gcf()
    fig.set_size_inches(16.0, 16.0)
    pylab.clf()
    pylab.gray()
    for m in range(100):
        pylab.subplot(10, 10, m + 1)
        pylab.imshow(np.clip(x[m], 0, 1), interpolation="none")
        pylab.axis("off")
    pylab.savefig("{}/{}.png".format(dir, filename))

def plot_signals(input_data, filename=None,downsamplefactor=1,n_columns=1):
    import pylab as pl
    n_rows = input_data.signal.n_channels()
    n_rows = int(n_rows/n_columns)
    print str(n_rows) + ' ' + str(n_columns)
    for row in range(n_rows):
        for col in range(n_columns):
            print (row)*n_columns+col+1
            pl.subplot(n_rows, n_columns, row*n_columns+col+1)
            if downsamplefactor==1:
                pl.plot(input_data.timebase, input_data.signal.get_channel(row*n_columns+col))
                pl.axis([-0.01,0.1,-5, 5])
            else:
                plotdata=input_data.signal.get_channel(row*n_columns+col)
                timedata=input_data.timebase
                pl.plot(timedata[0:len(timedata):downsamplefactor], plotdata[0:len(timedata):downsamplefactor])
                pl.axis([-0.01,0.1,-5,5])
    if filename != None:
        pl.savefig(filename)
    else:
        pl.show()

def visualize_x(reconstructed_x_batch, image_width=28, image_height=28, image_channel=1, dir=None):
    if dir is None:
        raise Exception()
    try:
        os.mkdir(dir)
    except:
        pass
    fig = pylab.gcf()
    fig.set_size_inches(16.0, 16.0)
    pylab.clf()
    if image_channel == 1:
        pylab.gray()
    for m in range(100):
        pylab.subplot(10, 10, m + 1)
        if image_channel == 1:
            pylab.imshow(reconstructed_x_batch[m].reshape((image_width, image_height)), interpolation="none")
        elif image_channel == 3:
            pylab.imshow(reconstructed_x_batch[m].reshape((image_channel, image_width, image_height)), interpolation="none")
        pylab.axis("off")
    pylab.savefig("%s/reconstructed_x.png" % dir)

def removeIllumination2(self, size, title = ''):
        out = ndimage.filters.gaussian_filter(self.image, size)
        pylab.figure()
        pylab.subplot(2,2,1)
        pylab.axis('off')
        pylab.imshow(self.image)
        pylab.subplot(2,2,2)
        pylab.axis('off')
        pylab.imshow(out)
        pylab.subplot(2,2,3)
        pylab.axis('off')
        pylab.imshow(self.image - out)
        pylab.subplot(2,2,4)
        pylab.axis('off')
        pylab.imshow(self.smooth - out)
        if title != '':
            pylab.savefig(title)
            pylab.close()
        else:
            pylab.show()
        self.smooth -= out
        return self.image - out

def draw2D(self, title, image=[]):
        pylab.figure()
        if image == []:
            pylab.imshow(self.image, 'gray')
        else:
            pylab.imshow(image, 'gray')
        pylab.axis('off')
        pylab.autoscale(False)
        for i in xrange(self.nComponents):
            xeq = lambda t: self.params[6 * i + 3] * np.cos(t) * np.cos(self.params[6 * i + 5]) + self.params[
                                                                                                      6 * i + 4] * np.sin(
                t) * np.sin(self.params[6 * i + 5]) + self.params[6 * i + 1]
            yeq = lambda t: - self.params[6 * i + 3] * np.cos(t) * np.sin(self.params[6 * i + 5]) + self.params[
                                                                                                        6 * i + 4] * np.sin(
                t) * np.cos(self.params[6 * i + 5]) + self.params[6 * i + 2]
            t = np.linspace(0, 2 * np.pi, 100)
            x = xeq(t)
            y = yeq(t)
            pylab.scatter(self.params[6 * i + 2], self.params[6 * i + 1], color='k')
            pylab.plot(y.astype(int), x.astype(int), self.colors[i] + '-')
        pylab.savefig(title)
        pylab.close()

def toBin(img, thr=200):
    #return (img < thr).astype(np.uint8)
    return (img < thr).max(axis=2).astype(np.uint8)

def GetLineSegments(itrim):
    """
    Segment an image in lines and interline spaces
    Returns lists of both (position width)
    """
        # sum along pixel lines
    asum = np.sum(itrim, axis=1)
    abin = asum > 0
    sp = []
    tx = []
    lastval=-1
    lastpos=-1
    for i in range(0, abin.size):
        if abin[i] != lastval:
            lastval = abin[i]
            if lastval:
                tx.append(np.array((i,0)))
                if i>0: 
                    sp[-1][1] = i-sp[-1][0]
            else:
                sp.append(np.array((i,0)))
                if i>0: 
                    tx[-1][1] = i-tx[-1][0]
    # set the last segment lenght
    if tx[-1][1] == 0: tx[-1][1] = itrim.shape[0] - tx[-1][0]
    if sp==[]:# empy if there is just one line in the image
        sp.append(np.array((0,0)))
    else:
        if sp[-1][1] == 0: sp[-1][1] = itrim.shape[0] - sp[-1][0]
    return tx, sp

def align(l1, l2, axis):
    if axis == 1: #horizontal alignment, we do not care about the right line end
    #cw = min(l2.shape[1],l1.shape[1])
    #l1 = l1[:,:cw]
    #l2 = l2[:,:cw]
        #compute correlation
        sc1 = np.sum(l1, axis=1-axis)
        sc2 = np.sum(l2, axis=1-axis)
        cor = np.correlate(sc1,sc2,"same")
        posErr =  np.argmax(cor)-sc1.shape[0]/2
        #place at right position
        if posErr > 0:
            l2c = l2.copy()
            l2c[:]=0
            l2c[:,posErr:] = l2[:,:-posErr]
            l2 = l2c
        elif posErr < 0:
            l1c = l1.copy()
            l1c[:]=0
            l1c[:,-posErr:] = l1[:,:posErr]
            l1=l1c
    else: #vertical alignment, we cate about both ends
        #compute correlation
        sc1 = np.sum(l1, axis=1-axis)
        sc2 = np.sum(l2, axis=1-axis)
        cor = np.correlate(sc1,sc2,"same")
        posErr =  np.argmax(cor)-sc1.shape[0]/2
        #place at right position
        if posErr > 0:
            l2c=l2.copy()
            l2c[:]=0
            l2c[posErr:,:] = l2[:-posErr,:]
            l2 = l2c
        elif posErr < 0:
            l1c=l1.copy()
            l1c[:]=0
            l1c[-posErr:,:]=l1[:posErr,:]
            l1 = l1c
    return posErr, l1, l2

def edgeplot(self, TT, ps):
        for ei,X in enumerate(self.edges):
            (i, j) = X[:2]
            Ta = TT[i]
            Tb = TT[j]
            plt.clf()
            if len(Ta) > 1000:
                nbins = 101
                ra = np.hstack((Ta.ra, Tb.ra))
                dec = np.hstack((Ta.dec, Tb.dec))
                H,xe,ye = np.histogram2d(ra, dec, bins=nbins)
                (matchRA, matchDec, dr,dd) = self.edge_matches(ei, goodonly=True)
                G,xe,ye = np.histogram2d(matchRA, matchDec, bins=(xe,ye))
                assert(G.shape == H.shape)
                img = antigray(H / H.max())
                img[G>0,:] = matplotlib.cm.hot(G[G>0] / H[G>0])
                ax = setRadecAxes(xe[0], xe[-1], ye[0], ye[-1])
                plt.imshow(img, extent=(min(xe), max(xe), min(ye), max(ye)),
                           aspect='auto', origin='lower', interpolation='nearest')
                plt.axis(ax)

            else:
                self.plotallstars([Ta,Tb])
                self.plotmatchedstars(ei)
                plt.xlabel('RA (deg)')
                plt.ylabel('Dec (deg)')
            ps.savefig()

    # one plot per edge

def hsvoffsets(self, TT, rad, apply=False):
        print 'hsv offsets plot'
        plt.clf()

        for ix,X in enumerate(self.edges):
            X = self.get_edge_dradec_arcsec(ix, corrected=apply, goodonly=True)
            (matchRA, matchDec, dra, ddec) = X

            print 'matchRA,Dec:', len(matchRA), len(matchDec)
            print 'dra,dec:', len(dra), len(ddec)

            for ra,dec,dr,dd in zip(matchRA, matchDec, dra, ddec):
                angle = arctan2(dd, dr) / (2.*pi)
                angle = fmod(angle + 1, 1.)
                mag = hypot(dd, dr)
                mag = min(1, mag/(0.5*rad))
                rgb = colorsys.hsv_to_rgb(angle, mag, 0.5)
                plt.plot([ra], [dec], '.', color=rgb, alpha=0.5)

        # legend in top-right corner.
        ax=plt.axis()
        xlo,xhi = plt.gca().get_xlim()
        ylo,yhi = plt.gca().get_ylim()
        # fraction
        keycx = xlo + 0.90 * (xhi-xlo)
        keycy = ylo + 0.90 * (yhi-ylo)
        keyrx = 0.1 * (xhi-xlo) / 1.4 # HACK
        keyry = 0.1 * (yhi-ylo)
        nrings = 5
        for i,(rx,ry) in enumerate(zip(np.linspace(0, keyrx, nrings), np.linspace(0, keyry, nrings))):
            nspokes = ceil(i / float(nrings-1) * 30)
            angles = np.linspace(0, 2.*pi, nspokes, endpoint=False)
            for a in angles:
                rgb = colorsys.hsv_to_rgb(a/(2.*pi), float(i)/(nrings-1), 0.5)
                plt.plot([keycx + rx*sin(a)], [keycy + ry*cos(a)], '.', color=rgb, alpha=1.)
        plt.axis(ax)
        plt.xlabel('RA (deg)')
        plt.ylabel('Dec (deg)')

def plotalignment(A, nbins=200, M=None, rng=None, doclf=True, docolorbar=True,
                  docutcircle=True, docontours=True, dologhist=False,
                  doaxlines=False, imshowargs={}):
    import pylab as plt
    from astrometry.util.plotutils import plothist, loghist
    if doclf:
        plt.clf()
    if M is None:
        M = A.match
    if dologhist:
        f = loghist
    else:
        f = plothist
    H,xe,ye = f(M.dra_arcsec*1000., M.ddec_arcsec*1000., nbins,
                range=rng, doclf=doclf, docolorbar=docolorbar,
                imshowargs=imshowargs)
    ax = plt.axis()
    if A is not None:
        # The EM fit is based on a subset of the matches;
        # draw the subset cut circle.
        if docutcircle:
            angle = np.linspace(0, 2.*pi, 360)
            plt.plot((A.cutcenter[0] + A.cutrange * np.cos(angle))*1000.,
                     (A.cutcenter[1] + A.cutrange * np.sin(angle))*1000., 'r-')
        if docontours:
            for i,c in enumerate(['b','c','g']*2):
                if i == A.ngauss:
                    break
                for nsig in [1,2]:
                    XY = A.getContours(nsig, c=i)
                    if XY is None:
                        break
                    X,Y = XY
                    plt.plot(X*1000., Y*1000., '-', color=c)#, alpha=0.5)
    if doaxlines:
        plt.axhline(0., color='b', alpha=0.5)
        plt.axvline(0., color='b', alpha=0.5)
    plt.axis(ax)
    plt.xlabel('dRA (mas)')
    plt.ylabel('dDec (mas)')
    return H,xe,ye

def _solve_lr(vlines, w, l, opt_options=OPTIMIZATION_OPTIONS, opt_method=OPTIMIZATION_METHOD, limit=0.3):
    """ Solve for the left and right edge displacement.
    This routine estimates the amount to move the upper left and right cornders of the image
     in a horizontal direction in order to make the given lines parallel and vertical.

    :param vlines: Lines that we want to map to vertical lines.
    :param w:  The width of the image
    :param l:   The height of the image
    :param opt_options: Optimization options passed into `minimize`
    :param opt_method: The optimization method.
    :param limit:  A limit on the amount of displacement -- beyond this and we will assume failure.
    :return: (dl, dr),   the horizontal displacement of the left and right corners.
    """
    if len(vlines) == 0:
        return 0, 0

    a = np.append(vlines[:, 0, :], np.ones((len(vlines), 1)), axis=1)
    b = np.append(vlines[:, 1, :], np.ones((len(vlines), 1)), axis=1)

    def objective(x):
        dl, dr = x
        Hv = np.linalg.inv(H_v(dl, dr, w, l))
        return np.sum(np.abs(Hv[0, :].dot(a.T) / Hv[2, :].dot(a.T) - Hv[0, :].dot(b.T) / Hv[2, :].dot(b.T)))

    res = minimize(objective, (0., 0.),
                   options=opt_options,
                   method=opt_method)
    dl, dr = res.x

    # Give up if the solution is not plausible (this indicates that the 'vlines' are too noisy
    if abs(dl) > limit * w:
        dl = 0
    if abs(dr) > limit * w:
        dr = 0
    return dl, dr

def _solve_ud(hlines, dl, dr, w, l, opt_options=OPTIMIZATION_OPTIONS, opt_method=OPTIMIZATION_METHOD, limit=0.3):
    """ Solve for the left top and bottom edge displacement.
    This routine estimates the amount to move the upper left and lower left corners of the image
    in a vertical direction in order to make the given lines parallel and horizontal.

    :param hlines: Lines that we want to map to horizontal lines.
    :param w:  The width of the image
    :param l:   The height of the image
    :param opt_options: Optimization options passed into `minimize`
    :param opt_method: The optimization method.
    :param limit:  A limit on the amount of displacement -- beyond this and we will assume failure.
                   It is expressed as a fraction of the image height.
    :return: (dl, dr),   the horizontal displacement of the left and right corners.
    """
    if len(hlines) == 0:
        return 0, 0

    a = np.append(hlines[:, 0, :], np.ones((len(hlines), 1)), axis=1)
    b = np.append(hlines[:, 1, :], np.ones((len(hlines), 1)), axis=1)

    Hv = np.linalg.inv(H_v(dl, dr, w, l))
    a = Hv.dot(a.T).T
    b = Hv.dot(b.T).T

    def objective(x):
        du, dd = x
        Hh = np.linalg.inv(H_h(du, dd, w, l))
        return np.sum(np.abs(Hh[1, :].dot(a.T) / Hh[2, :].dot(a.T) - Hh[1, :].dot(b.T) / Hh[2, :].dot(b.T)))

    res = minimize(objective, (0., 0.),
                   options=opt_options,
                   method=opt_method)
    du, dd = res.x

    # Give up if the result is not plausible. We are better off nor warping.
    if abs(du) > limit * l:
        du = 0
    if abs(dd) > limit * l:
        dd = 0
    return du, dd

def _cut_windows_horizontally(self, s, win_strip):
        win_horizontal_scores = []
        if len(self.heights) > 0:
            win_horizontal_scores = np.percentile(win_strip[self.top:self.bottom], 85, axis=0)
            runs, starts, values = run_length_encode(win_horizontal_scores > 0.5)
            starts += s
            win_widths = runs[np.atleast_1d(values)]
            win_widths = np.atleast_1d(win_widths)
            win_lefts = np.atleast_1d(starts[values])
            if len(win_widths) > 0:
                win_left = win_lefts[0]
                win_right = win_lefts[-1] + win_widths[-1]
                win_horizontal_spacing = np.diff(win_lefts).mean() if len(win_lefts) > 1 else 0
                # win_width = win_widths.mean()
            else:
                win_left = win_right = win_horizontal_spacing = -1  # win_width = -1
        else:
            win_widths = win_lefts = []
            win_left = win_right = win_horizontal_spacing = -1

        self.horizontal_spacing = int(win_horizontal_spacing)
        self.left = int(win_left)
        self.right = int(win_right)
        self.horizontal_scores = win_horizontal_scores
        self.lefts = np.array(win_lefts)
        self.widths = np.array(win_widths)

def _create_mini_facade(self, left, right, wall_colors):
        door_strip = i12.door(self.facade_layers)[:, left:right].copy()
        shop_strip = i12.shop(self.facade_layers)[:, left:right]
        door_strip = np.max((door_strip, shop_strip), axis=0)
        win_strip = self.window_scores[:, left:right].copy()
        sky_strip = self._sky_mask[:, left:right].copy()
        rgb_strip = wall_colors[:, left:right]
        win_strip[:, :1] = win_strip[:, -1:] = 0  # edge effects
        sky_strip[:, :1] = sky_strip[:, -1:] = 0  # edge effects

        facade = FacadeCandidate(self, left, right, sky_strip, door_strip, win_strip, rgb_strip)
        facade.find_regions(self.facade_layers)
        return facade

def plot_facade_cuts(self):

        facade_sig = self.facade_edge_scores.sum(0)
        facade_cuts = find_facade_cuts(facade_sig, dilation_amount=self.facade_merge_amount)
        mu = np.mean(facade_sig)
        sigma = np.std(facade_sig)

        w = self.rectified.shape[1]
        pad=10

        gs1 = pl.GridSpec(5, 5)
        gs1.update(wspace=0.5, hspace=0.0)  # set the spacing between axes.

        pl.subplot(gs1[:3, :])
        pl.imshow(self.rectified)
        pl.vlines(facade_cuts, *pl.ylim(), lw=2, color='black')
        pl.axis('off')
        pl.xlim(-pad, w+pad)

        pl.subplot(gs1[3:, :], sharex=pl.gca())
        pl.fill_between(np.arange(w), 0, facade_sig, lw=0, color='red')
        pl.fill_between(np.arange(w), 0, np.clip(facade_sig, 0, mu+sigma), color='blue')
        pl.plot(np.arange(w), facade_sig, color='blue')

        pl.vlines(facade_cuts, facade_sig[facade_cuts], pl.xlim()[1], lw=2, color='black')
        pl.scatter(facade_cuts, facade_sig[facade_cuts])

        pl.axis('off')

        pl.hlines(mu, 0, w, linestyle='dashed', color='black')
        pl.text(0, mu, '$\mu$ ', ha='right')

        pl.hlines(mu + sigma, 0, w, linestyle='dashed', color='gray',)
        pl.text(0, mu + sigma, '$\mu+\sigma$ ', ha='right')
        pl.xlim(-pad, w+pad)

def __call__(self, X, Y=None):
        XX = np.sum(X * X, axis=1)[:,np.newaxis]
        if Y is None:
            Y = X
            YY = XX.T
        else:
            YY = np.sum(Y * Y, axis=1)[np.newaxis,:]
        distances = XX + YY # Using broadcasting
        distances -= 2 * np.dot(X, Y.T)
        distances = np.maximum(distances, 0)
        return np.exp(- self.gamma * distances)

def drop_num_features(df):
        # Drop all categorical feature helping columns ('Num')
        # Todo: is it defined when importing data set? _feature_names_num
        for feature_name in HousePrices._feature_names_num:
            df = df.drop([feature_name], axis=1)
        return df

def features_with_null_logical(df, axis=1):
        row_length = len(df._get_axis(0))
        # Axis to count non null values in. aggregate_axis=0 implies counting for every feature
        aggregate_axis = 1 - axis
        features_non_null_series = df.count(axis=aggregate_axis)
        # Whenever count() differs from row_length it implies a null value exists in feature column and a False in mask
        mask = row_length == features_non_null_series
        return mask

def svm_figure_generate(w, b, support_vectors, X):
    k = - w[0]/w[1]
    x = np.linspace(-5, 5)
    y = k*x - b/w[1]
    sv_1 = support_vectors[0]
    yy_down = k*x + (sv_1[1]-k*sv_1[0])
    sv_2 = support_vectors[-1]
    yy_up = k*x + (sv_2[1] - k*sv_2[0])
    pl.plot(x, y, 'k-')
    pl.plot(x, yy_up, 'k--')
    pl.plot(x, yy_down, 'k--')
    pl.scatter(support_vectors[:, 0], support_vectors[:, 1], s=80, facecolor='none')
    pl.scatter(X[:, 0], X[:, 1], c='Y', cmap=pl.cm.Paired)
    pl.axis('tight')
    pl.show()

def adjust_layout(self) :
        x0, x1, y0, y1 = plt.axis()
        plot_margin_x = 0.01 * float(x1)
        plot_margin_y = 0.01 * float(y1)
        plt.axis((x0 - plot_margin_x, x1 + plot_margin_x, y0, y1 + plot_margin_y))
        plt.tight_layout()

def cmPlot(targ_ra, targ_dec, data, iso, g_radius, nbhd, type):
    """Color-magnitude plot"""

    angsep = ugali.utils.projector.angsep(targ_ra, targ_dec, data['RA'], data['DEC'])
    annulus = (angsep > g_radius) & (angsep < 1.)

    mag_g = data[mag_g_dred_flag]
    mag_r = data[mag_r_dred_flag]

    if type == 'stars':
        filter = star_filter(data)
        plt.title('Stellar Color-Magnitude')
    elif type == 'galaxies':
        filter = galaxy_filter(data)
        plt.title('Galactic Color-Magnitude')

    iso_filter = (iso.separation(mag_g, mag_r) < 0.1)

    # Plot background objects
    plt.scatter(mag_g[filter & annulus] - mag_r[filter & annulus], mag_g[filter & annulus], c='k', alpha=0.1, edgecolor='none', s=1)

    # Plot isochrone
    ugali.utils.plotting.drawIsochrone(iso, lw=2, label='{} Gyr, z = {}'.format(iso.age, iso.metallicity))

    # Plot objects in nbhd
    plt.scatter(mag_g[filter & nbhd] - mag_r[filter & nbhd], mag_g[filter & nbhd], c='g', s=5, label='r < {:.3f}$^\circ$'.format(g_radius))

    # Plot objects in nbhd and near isochrone
    plt.scatter(mag_g[filter & nbhd & iso_filter] - mag_r[filter & nbhd & iso_filter], mag_g[filter & nbhd & iso_filter], c='r', s=5, label='$\Delta$CM < 0.1')

    plt.axis([-0.5, 1, 16, 24])
    plt.gca().invert_yaxis()
    plt.gca().set_aspect(1./4.)
    plt.legend(loc='upper left')
    plt.xlabel('g-r (mag)')
    plt.ylabel('g (mag)')

def nullspace(A, eps=1e-15):
    u, s, vh = sp.linalg.svd(A,full_matrices=1,compute_uv=1)
    # Pad so that we get the nullspace of a wide matrix. 
    N = A.shape[1]
    K = s.shape[0]
    if K < N:
        s[K+1:N] = 0
        s2 = np.zeros((N))
        s2[0:K] = s
        s = s2
    null_mask = (s <= eps)
    null_space = sp.compress(null_mask, vh, axis=0)
    return sp.transpose(null_space)

# return smallest singular vector of A (or the nullspace if A is 2x3)

def analytic_twopath(self,p,rate1,rate2):
        '''
        Population of neurons receive input from two pathways, the first path is gated-on
        rate1 and rate2 are the input rate of each pathway
        First we need to convert the input rate into conductance,
        the dend_IO(exc, inh) function takes total excitatory and inhibitory
        conductances as inputs
        '''
        # number of synapses
        num_syn = 15
        g_exc = p['g_exc']*num_syn
        # gating variable
        s1 = MCM.meansNMDA(rate1)
        s2 = MCM.meansNMDA(rate2)

        # Total conductance input
        Exc1 = self.Exc1_raw*s1*g_exc # nS
        Exc2 = self.Exc2_raw*s2*g_exc # nS

        Exc = Exc1+Exc2

        #frac_proj = 0.1 # fraction projection
        N_proj = p['frac_proj']*self.params['n_pyr']
        N_proj0 = np.floor(N_proj)
        N_proj0 = min((N_proj0,self.params['n_pyr']-1))
        N_proj0 = max((N_proj0,0))

        DendV = dend_IO(Exc[:(N_proj0+1)*self.params['n_dend_each']],
                            self.Inh1[:(N_proj0+1)*self.params['n_dend_each']])
        meanDendV = DendV.reshape(N_proj0+1,self.params['n_dend_each']).mean(axis=1)
        SomaR = soma_fv(meanDendV)

        # Make sure firing rate depend smoothly on frac_proj
        rboth = (SomaR[:N_proj0].sum()+SomaR[N_proj0]*(N_proj-N_proj0))/N_proj

        return rboth