def plot_earth_moon_obliquity():
	from numpy import linspace, array, arange
	from pylab import rc, plot, xlim, xticks, xlabel, ylabel
	import matplotlib.pyplot as plt
	import mpmath
	rc('text', usetex=True)
	fig = plt.figure()
	# Initial params.
	#params = {'m2':1.9e27,'r2':70000.e3,'rot2':0.00017585125448028,\
	#	'r1':1,'rot1':6.28*100,'i_a':0.5,'ht': 0,\
	#	'a':421700E3 / 4,'e':0.01,'i':0.8,'h':2}
	rot1 = 7.978282018665196e-08
	rot2 = 2.972e-06
	params = {'m2':1.98892e30,'r2':695500000.0,'rot2':rot2,\
		'r1':384748.e3,'rot1':rot1,'i_a':10. * 2*mpmath.pi()/360,'ht': 0,\
		'a':149597870700. * 0.1,'e':0.017,'i':7. * 2*mpmath.pi()/360,'h':mpmath.pi()/2}
	# Set parameters.
	sp.parameters = params
	# Period.
	period = sp.wp_period
	ob = sp.obliquity_time
	lspace = linspace(0,5*period,250)
	xlim(0,float(lspace[-1]))
	xlabel(r'$\textnormal{Time (Ma)}$')
	ylabel(r'$\textnormal{Obliquity }\left( ^\circ \right)$')
	ob_series = array([ob(t).real*(360/(2*mpmath.pi())) for t in lspace])
	plot(lspace,ob_series,'k-',linewidth=2)
	xticks(arange(xlim()[0],xlim()[1],365*24*3600*100000000),[r'$'+str(int(_)*100)+r'$' for _ in [t[0] for t in enumerate(arange(xlim()[0],xlim()[1],365*24*3600*100000000))]])

def getOptCandGamma(cv_train, cv_label):
    print "Finding optimal C and gamma for SVM with RBF Kernel"
    C_range = 10.0 ** np.arange(-2, 9)
    gamma_range = 10.0 ** np.arange(-5, 4)
    param_grid = dict(gamma=gamma_range, C=C_range)
    cv = StratifiedKFold(y=cv_label, n_folds=40)

    # Use the svm.SVC() as the cost function to evaluate parameter choices
    # NOTE: Perhaps we should run computations in parallel if needed. Does it
    # do that already within the class?
    grid = GridSearchCV(svm.SVC(), param_grid=param_grid, cv=cv)
    grid.fit(cv_train, cv_label)

    score_dict = grid.grid_scores_
    scores = [x[1] for x in score_dict]
    scores = np.array(scores).reshape(len(C_range), len(gamma_range))
    pl.figure(figsize=(8,6))
    pl.subplots_adjust(left=0.05, right=0.95, bottom=0.15, top=0.95)
    pl.imshow(scores, interpolation='nearest', cmap=pl.cm.spectral)
    pl.xlabel('gamma')
    pl.ylabel('C')
    pl.colorbar()
    pl.xticks(np.arange(len(gamma_range)), gamma_range, rotation=45)
    pl.yticks(np.arange(len(C_range)), C_range)
    pl.show()

    print "The best classifier is: ", grid.best_estimator_

def createPlot(dataY, dataX, ticksX, annotations, axisY, axisX, dostep, doannotate):
    if not ticksX:
        ticksX = dataX
    
    if dostep:
        py.step(dataX, dataY, where='post', linestyle='-', label=axisY) # where=post steps after point
    else:
        py.plot(dataX, dataY, marker='o', ms=5.0, linestyle='-', label=axisY)
    
    if annotations and doannotate:
        for note, x, y in zip(annotations, dataX, dataY):
            py.annotate(note, (x, y), xytext=(2,2), xycoords='data', textcoords='offset points')

    py.xticks(np.arange(1, len(dataX)+1), ticksX, horizontalalignment='left', rotation=30)
    leg = py.legend()
    leg.draggable()
    py.xlabel(axisX)
    py.ylabel('time (s)')

    # Set X axis tick labels as rungs
    #print zip(dataX, dataY)
  
    py.draw()
    py.show()
    
    return

def plotMatrix(matrix, color_scheme, row_headers, col_headers, vmin, vmax, options):

    pylab.imshow(matrix,
                 cmap=color_scheme,
                 origin='lower',
                 vmax=vmax,
                 vmin=vmin,
                 interpolation='nearest')

    # offset=0: x=center,y=center
    # offset=0.5: y=top/x=right
    offset = 0.0

    if options.xticks:
        pylab.xticks([offset + x for x in range(len(options.xticks))],
                     options.xticks,
                     rotation="vertical",
                     fontsize="8")
    else:
        if col_headers and len(col_headers) < 100:
            pylab.xticks([offset + x for x in range(len(col_headers))],
                         col_headers,
                         rotation="vertical",
                         fontsize="8")

    if options.yticks:
        pylab.yticks([offset + y for y in range(len(options.yticks))],
                     options.yticks,
                     fontsize="8")
    else:
        if row_headers and len(row_headers) < 100:
            pylab.yticks([offset + y for y in range(len(row_headers))],
                         row_headers,
                         fontsize="8")

    def plot_cost(self):
        if self.show_cost not in self.train_outputs[0][0]:
            raise ShowNetError("Cost function with name '%s' not defined by given convnet." % self.show_cost)
        train_errors = [o[0][self.show_cost][self.cost_idx] for o in self.train_outputs]
        test_errors = [o[0][self.show_cost][self.cost_idx] for o in self.test_outputs]

        numbatches = len(self.train_batch_range)
        test_errors = numpy.row_stack(test_errors)
        test_errors = numpy.tile(test_errors, (1, self.testing_freq))
        test_errors = list(test_errors.flatten())
        test_errors += [test_errors[-1]] * max(0,len(train_errors) - len(test_errors))
        test_errors = test_errors[:len(train_errors)]

        numepochs = len(train_errors) / float(numbatches)
        pl.figure(1)
        x = range(0, len(train_errors))
        pl.plot(x, train_errors, 'k-', label='Training set')
        pl.plot(x, test_errors, 'r-', label='Test set')
        pl.legend()
        ticklocs = range(numbatches, len(train_errors) - len(train_errors) % numbatches + 1, numbatches)
        epoch_label_gran = int(ceil(numepochs / 20.)) # aim for about 20 labels
        epoch_label_gran = int(ceil(float(epoch_label_gran) / 10) * 10) # but round to nearest 10
        ticklabels = map(lambda x: str((x[1] / numbatches)) if x[0] % epoch_label_gran == epoch_label_gran-1 else '', enumerate(ticklocs))

        pl.xticks(ticklocs, ticklabels)
        pl.xlabel('Epoch')
#        pl.ylabel(self.show_cost)
        pl.title(self.show_cost)

def param_sweeping(clf, obj, X, y, param_dist, metric, param, clfName):
	'''Plot a parameter sweeping (ps) curve with the param_dist as a axis, and the scoring based on metric as y axis.
	Keyword arguments:
	clf - - classifier
	X - - feature matrix
	y - - target array
	param - - a parameter of the classifier
	param_dist - - the parameter distribution of param
	clfName - - the name of the classifier
	metric - - the metric we use to evaluate the performance of the classifiers
	obj - - the name of the dataset we are using'''
	scores = []
	for i in param_dist:
		y_true = []
		y_pred = []
		# new classifer each iteration
		newclf = eval("clf.set_params("+ param + "= i)")
		y_pred, y_true, gs_score_list, amp = testAlgo(newclf, X, y, clfName)
		mean_fpr, mean_tpr, mean_auc = plot_unit_prep(y_pred, y_true, metric)
		scores.append(mean_auc)
		print("Area under the ROC curve : %f" % mean_auc)
	fig = pl.figure(figsize=(8,6),dpi=150)
	paramdist_len = len(param_dist)
	pl.plot(range(paramdist_len), scores, label = 'Parameter Sweeping Curve')
	pl.xticks(range(paramdist_len), param_dist, size = 15, rotation = 45)
	pl.xlabel(param.upper(),fontsize=30)
	pl.ylabel(metric.upper(),fontsize=30)
	pl.title('Parameter Sweeping Curve',fontsize=25)
	pl.legend(loc='lower right')
	pl.tight_layout()
	fig.savefig('plots/'+obj+'/'+ clfName +'_' + param +'_'+'ps.pdf')
	pl.show()

    def plot(self, filesuffix=('.png',)):

        pylab.figure()

        kPluses = 10**numpy.linspace(0, 3, 100)
        kMinuses = 10**numpy.linspace(6, 9, 100)
        figureOfMerits = numpy.zeros((len(kPluses), len(kMinuses), 4), 'd')
        for i, kPlus in enumerate(kPluses):
            for j, kMinus in enumerate(kMinuses):
                figureOfMerits[i, j, :] = self.figureOfMerit(self.generateData({'kPlus' : kPlus, 'kMinus' : kMinus}))
        data = self.generateData({'kPlus' : kPluses[0], 'kMinus' : kMinuses[0]})

        self._contourf(kMinuses, kPluses, figureOfMerits, data)

        pylab.xticks((10**6, 10**7, 10**8, 10**9), fontsize=self.fontsize)
        pylab.yticks((10**0, 10**1, 10**2, 10**3), fontsize=self.fontsize)
        pylab.xlabel(r'$k^-$ $\left(1\per\metre\right)$', fontsize=self.fontsize)
        pylab.ylabel(r'$k^+$ $\left(\power{\metre}{3}\per\mole\cdot\second\right)$', fontsize=self.fontsize)

        pylab.text(2 * 10**6, 7 * 10**2, r'I', fontsize=self.fontsize)
        pylab.text(3 * 10**7, 7 * 10**2, r'II', fontsize=self.fontsize)
        pylab.text(6 * 10**8, 7 * 10**2, r'III', fontsize=self.fontsize)
        pylab.text(6 * 10**8, 7 * 10**1, r'IV', fontsize=self.fontsize)

        for fP, kPlus, paxeslabel in ((1.143,3.51e+00, False), (0.975, 9.33e+00, False), (0.916, 3.51e+01, False), (0.89, 9.33e+01, False), (0.87, 3.51e+02, True)):
            for fM, kMinus, maxeslabel in ((1.4, 2.48e+06, False), (1.07, 7.05e+6, False), (0.96, 2.48e+07, False), (0.91, 7.05e+7, False), (0.88, 2.48e+08, True)):
                xpos = (numpy.log10(kMinus) - 6.) / 3. *  fM
                ypos = numpy.log10(kPlus) / 3. * fP        
                self.makeBackGroundPlot({'kPlus' : kPlus, 'kMinus' : kMinus}, xpos, ypos, axeslabel=paxeslabel and maxeslabel)

        for fs in filesuffix:
            pylab.savefig('kPlusVkMinus' + fs)
            pylab.close('all')

    def plot_predictions(self):
        data = self.get_next_batch(train=False)[2] # get a test batch
        num_classes = self.test_data_provider.get_num_classes()
        NUM_ROWS = 2
        NUM_COLS = 4
        NUM_IMGS = NUM_ROWS * NUM_COLS
        NUM_TOP_CLASSES = min(num_classes, 4) # show this many top labels

        label_names = self.test_data_provider.batch_meta['label_names']
        if self.only_errors:
            preds = n.zeros((data[0].shape[1], num_classes), dtype=n.single)
        else:
            preds = n.zeros((NUM_IMGS, num_classes), dtype=n.single)
            rand_idx = nr.randint(0, data[0].shape[1], NUM_IMGS)
            print rand_idx
            data[0] = n.require(data[0][:,rand_idx], requirements='C')
            data[1] = n.require(data[1][:,rand_idx], requirements='C')
        data += [preds]
        temp = data[0]
        print data
        print temp.ndim,temp.shape,temp.size
        # Run the model
        self.libmodel.startFeatureWriter(data, self.sotmax_idx)
        self.finish_batch()

        fig = pl.figure(3)
        fig.text(.4, .95, '%s test case predictions' % ('Mistaken' if self.only_errors else 'Random'))
        if self.only_errors:
            err_idx = nr.permutation(n.where(preds.argmax(axis=1) != data[1][0,:])[0])[:NUM_IMGS] # what the net got wrong
            data[0], data[1], preds = data[0][:,err_idx], data[1][:,err_idx], preds[err_idx,:]

        data[0] = self.test_data_provider.get_plottable_data(data[0])
        for r in xrange(NUM_ROWS):
            for c in xrange(NUM_COLS):
                img_idx = r * NUM_COLS + c
                if data[0].shape[0] <= img_idx:
                    break
                pl.subplot(NUM_ROWS*2, NUM_COLS, r * 2 * NUM_COLS + c + 1)
                pl.xticks([])
                pl.yticks([])
                try:
                    img = data[0][img_idx,:,:,:]
                except IndexError:
                    # maybe greyscale?
                    img = data[0][img_idx,:,:]
                pl.imshow(img, interpolation='nearest')
                true_label = int(data[1][0,img_idx])

                img_labels = sorted(zip(preds[img_idx,:], label_names), key=lambda x: x[0])[-NUM_TOP_CLASSES:]
                pl.subplot(NUM_ROWS*2, NUM_COLS, (r * 2 + 1) * NUM_COLS + c + 1, aspect='equal')

                ylocs = n.array(range(NUM_TOP_CLASSES)) + 0.5
                height = 0.5
                width = max(ylocs)
                pl.barh(ylocs, [l[0]*width for l in img_labels], height=height, \
                        color=['r' if l[1] == label_names[true_label] else 'b' for l in img_labels])
                pl.title(label_names[true_label])
                pl.yticks(ylocs + height/2, [l[1] for l in img_labels])
                pl.xticks([width/2.0, width], ['50%', ''])
                pl.ylim(0, ylocs[-1] + height*2)

def plot_C_gamma_grid_search(grid, C_range, gamma_range, score):
    '''
    Plots the scores computed on a grid. 
    
    Arguments: 
        grid - the grid search object created using GridSearchCV()
        C_range - the C parameter range 
        gamma_range - the gamma parameter range 
        score - the scoring function  
        
    
    '''

    # grid_scores_ contains parameter settings and scores
    # We extract just the scores
    scores = [x[1] for x in grid.grid_scores_]
    scores = np.array(scores).reshape(len(C_range), len(gamma_range))
    
    # draw heatmap of accuracy as a function of gamma and C
    pl.figure(figsize=(8, 6))
    pl.subplots_adjust(left=0.05, right=0.95, bottom=0.15, top=0.95)
    pl.imshow(scores, interpolation='nearest', cmap=pl.cm.spectral)
    pl.title("Grid search on C and gamma for best %s" % score)
    pl.xlabel('gamma')
    pl.ylabel('C')
    pl.colorbar()
    pl.xticks(np.arange(len(gamma_range)), gamma_range, rotation=45)
    pl.yticks(np.arange(len(C_range)), C_range)
    
    pl.show()

def paint1(data):
    order_data = sorted(data.items(), key=lambda x: x[0])
    print order_data

    x = [i[0] for i in order_data]

    y0 = [i[1][0] for i in order_data]
    y1 = [i[1][1] for i in order_data]
    y2 = [i[1][2] for i in order_data]
    y3 = [i[1][3] for i in order_data]

    #print "x", x
    #print "y0", y0
    #print "y1", y1


    pl.plot(x, y0, 'r', label="NONE")
    pl.plot(x, y1, 'b', label="MPI")
    pl.plot(x, y2, 'g', label="OMP")
    pl.plot(x, y3, 'y', label="OMP_MPI")

    pl.legend(loc="upper left")
    pl.xticks(list(xrange(3, 11)), [str(i) for i in xrange(3, 11)])
	
    pl.xlabel('String Length')
    pl.ylabel('Time')

    pl.show()

def paint2(data):
    order_data = sorted(data.items(), key=lambda x: x[0])
    print order_data
    dic = ['', "NONE", "MPI", "OMP", "OMP_MPI", '']

    x = [i for i in xrange(4)]
    y0 = order_data[0][1]
    y1 = order_data[1][1]
    y2 = order_data[2][1]
    y3 = order_data[3][1]
    y4 = order_data[4][1]
    #y4 = order_data[4][1]

    #print "x", x
    #print "y0", y0
    #print "y1", y1


    pl.plot(x, y0, 'r', label="length=4")
    pl.plot(x, y1, 'b', label="length=5")
    pl.plot(x, y2, 'g', label="length=6")
    pl.plot(x, y3, 'y', label="length=7")
    pl.plot(x, y4, 'c', label="length=8")
    #pl.plot(x, y4, 'k', label="length=9")

    pl.xticks(list(xrange(-1, 5)), list(dic))

    pl.legend(loc="upper right")
	
    pl.xlabel('Compiling Modes')
    pl.ylabel('Time')
	
    pl.show()

    def plot_bases(self, autoscale=True, stampsize=None):
        import pylab as plt

        N = len(self.psfbases)
        cols = int(np.ceil(np.sqrt(N)))
        rows = int(np.ceil(N / float(cols)))
        plt.clf()
        plt.subplots_adjust(hspace=0, wspace=0)

        cut = 0
        if stampsize is not None:
            H, W = self.shape
            assert H == W
            cut = max(0, (H - stampsize) / 2)

        ima = dict(interpolation="nearest", origin="lower")
        if autoscale:
            mx = self.psfbases.max()
            ima.update(vmin=-mx, vmax=mx)
        nil, xpows, ypows = self.polynomials(0.0, 0.0, powers=True)
        for i, (xp, yp, b) in enumerate(zip(xpows, ypows, self.psfbases)):
            plt.subplot(rows, cols, i + 1)

            if cut > 0:
                b = b[cut:-cut, cut:-cut]
            if autoscale:
                plt.imshow(b, **ima)
            else:
                mx = np.abs(b).max()
                plt.imshow(b, vmin=-mx, vmax=mx, **ima)
            plt.xticks([])
            plt.yticks([])
            plt.title("x^%i y^%i" % (xp, yp))
        plt.suptitle("PsfEx eigen-bases")

 def postProcess(self):
     coreName = self.study.cacheDir + os.sep + "hz-" + self.name
     f = open(coreName)
     vals = {}
     for l in f:
         toks = l.split("\t")
         vals[toks[0], toks[1]] = toks[2:]
     f.close()
     for pop in self.pops:
         pylab.clf()
         pylab.title(pop)
         labels = []
         ehzs = []
         ohzs = []
         ns = []
         for indiv in self.study.pops.getIndivs(pop):
             labels.append("\n".join(list(indiv)))
             o, e, n, f = vals[indiv]
             ehz = 1 - float(e) / int(n)
             ohz = 1 - float(o) / int(n)
             ehzs.append(ehz)
             ohzs.append(ohz)
             ns.append(int(n))
         pylab.xticks(rotation=90, fontsize="x-small")
         pylab.legend()
         pylab.plot(ehzs, "+", label="ExpHe")
         pylab.plot(ohzs, "+", label="ObsHe")
         a2 = pylab.twinx()
         a2.plot(ns, ".")
         pylab.xticks(list(range(len(labels))), labels)
         xmin, xmax = pylab.xlim()
         pylab.xlim(xmin - 1, xmax + 1)
         pylab.savefig(coreName + "-" + pop + ".png")

def plot_io_obliquity_2():
	from numpy import linspace, array, arange
	from pylab import rc, plot, xlim, xticks, xlabel, ylabel
	import matplotlib.pyplot as plt
	import mpmath
	rc('text', usetex=True)
	fig = plt.figure()
	# Initial params.
	#params = {'m2':1.9e27,'r2':70000.e3,'rot2':0.00017585125448028,\
	#	'r1':1,'rot1':6.28*100,'i_a':0.5,'ht': 0,\
	#	'a':421700E3 / 4,'e':0.01,'i':0.8,'h':2}
	params = {'m2':1.9e27,'r2':70000.e3,'rot2':0.00017585125448028,\
		'r1':1821E3,'rot1':6.28/(21 * 3600),'i_a':0.17453,'ht': 0,\
		'a':421700E3,'e':0.05,'i':0.349,'h':0.35}
	# Set parameters.
	sp.parameters = params
	# Period.
	period = sp.wp_period
	ob = sp.obliquity_time
	lspace = linspace(0,5*period,250)
	xlim(0,float(lspace[-1]))
	xlabel(r'$\textnormal{Time (Ma)}$')
	ylabel(r'$\textnormal{Obliquity }\left( ^\circ \right)$')
	ob_series = array([ob(t).real*(360/(2*mpmath.pi())) for t in lspace])
	plot(lspace,ob_series,'k-',linewidth=2)
	xticks(arange(xlim()[0],xlim()[1],365*24*3600*1000000),[r'$'+str(int(_)*1)+r'$' for _ in [t[0] for t in enumerate(arange(xlim()[0],xlim()[1],365*24*3600*1000000))]])

def makeboxplot(filteredclusts, dblibrary, figname, pool=False):
    '''takes a filtered dict of clusts worth keeping and creates a boxplot of either by lane (default) or pool'''
    indiv_cluster_count = defaultdict(int) 
    for clust, inddict in filteredclusts.items():
        for ind, reads in inddict.items():
            if ind in indiv_cluster_count.keys():
                indiv_cluster_count[ind]+=1
            else:
                indiv_cluster_count[ind]+=1 
    
    t = gdata_tools.get_table_as_dict(dblibrary)
    db_ind_countd = Util.countdict([d['sampleid'] for d in t if d['sampleid'] in indiv_cluster_count.keys()[3]]) #creates a table of individual dicts from google spreadsheet
    indiv_by_group = defaultdict(list)
    for d in t:
        if 'pool' in d:
            indkey = (d.get('flowcell',None),d.get('lane',None),d.get('index',None),d.get('sampleid',None))
            if indkey in indiv_cluster_count:
                if pool == True:
                    indiv_by_group[(d['flowcell'],d['lane'],d.get('index',None),d['pool'])].append(indiv_cluster_count[indkey]) 
                else:
                    indiv_by_group[(d['flowcell'],d['lane'],d.get('index',None))].append(indiv_cluster_count[indkey])
    
    boxes = []
    labels = []
    for group,indcounts in indiv_by_group.items():
        boxes.append(indcounts)
        labels.append(group)
    boxplt = pylab.figure(1)
    pylab.boxplot(boxes)
    pylab.xticks(arange(1,(len(labels)+1)),labels,fontsize='small') #legend with best location (0) if pools
    boxplt.savefig(figname)

def align_subplots(N,M,xlim=None, ylim=None):
    """make all of the subplots have the same limits, turn off unnecessary ticks"""
    #find sensible xlim,ylim
    if xlim is None:
        xlim = [np.inf,-np.inf]
        for i in range(N*M):
            pb.subplot(N,M,i+1)
            xlim[0] = min(xlim[0],pb.xlim()[0])
            xlim[1] = max(xlim[1],pb.xlim()[1])
    if ylim is None:
        ylim = [np.inf,-np.inf]
        for i in range(N*M):
            pb.subplot(N,M,i+1)
            ylim[0] = min(ylim[0],pb.ylim()[0])
            ylim[1] = max(ylim[1],pb.ylim()[1])

    for i in range(N*M):
        pb.subplot(N,M,i+1)
        pb.xlim(xlim)
        pb.ylim(ylim)
        if (i)%M:
            pb.yticks([])
        else:
            removeRightTicks()
        if i<(M*(N-1)):
            pb.xticks([])
        else:
            removeUpperTicks()

def plot_importances(imp, clfName, obj):
    imp=np.vstack(imp)
    print imp
    mean_importance = np.mean(imp,axis=0)
    std_importance = np.std(imp,axis=0)
    indices = np.argsort(mean_importance)[::-1]
    print indices
    print featureNames
    featureList = []
    # num_features = len(featureNames)
    print("Feature ranking:")
    for f in range(num_features):
        featureList.append(featureNames[indices[f]])
        print("%d. feature %s (%.2f)" % (f, featureNames[indices[f]], mean_importance[indices[f]]))
    fig = pl.figure(figsize=(8,6),dpi=150)
    pl.title("Feature importances",fontsize=30)
    pl.bar(range(num_features), mean_importance[indices],
            yerr = std_importance[indices], color=paired[0], align="center",
            edgecolor=paired[0],ecolor=paired[1])
    pl.xticks(range(num_features), featureList, size=15,rotation=90)
    pl.ylabel("Importance",size=30)
    pl.yticks(size=20)
    pl.xlim([-1, num_features])
    # fix_axes()
    pl.tight_layout()
    save_path = 'plots/'+obj+'/'+clfName+'_feature_importances.pdf'
    fig.savefig(save_path)

def plot_genome(out_file, data_file, samples, dpi=300, screen=False):
    if screen: PL.rcParams.update(PLOT_PARAMS_SCREEN)
    LOG.info("plot_genome - out_file=%s, data_file=%s, samples=%s, dpi=%d"%(out_file, data_file, str(samples), dpi))
    colors = 'bgryckbgryck'

    data = read_posterior(data_file)
    if samples is None or len(samples) == 0: samples = data.keys()
    if len(samples) == 0: return

    PL.figure(None, [14, 4])
    right_end = 0 # rightmost plotted base pair
    for chrm in sort_chrms(data.values()[0]): # for chromosomes in ascending order
        max_site = max(data[samples[0]][chrm]['L']) # length of chromosome
        for s, sample in enumerate(samples): # plot all samples
            I = SP.where(SP.array(data[sample][chrm]['SD']) < 0.3)[0] # at sites that have confident posteriors
            PL.plot(SP.array(data[sample][chrm]['L'])[I] + right_end, SP.array(data[sample][chrm]['AF'])[I], alpha=0.4, color=colors[s], lw=2) # offset by the end of last chromosome
        if right_end > 0: PL.plot([right_end, right_end], [0,1], 'k--', lw=0.4, alpha=0.2) # plot separators between chromosomes
        new_right = right_end + max(data[sample][chrm]['L'])
        PL.text(right_end + 0.5*(new_right - right_end), 0.9, str(chrm), horizontalalignment='center')
        right_end = new_right # update rightmost end
    PL.plot([0,right_end], [0.5,0.5], 'k--', alpha=0.3)
    PL.xlim(0,right_end)
    xrange = SP.arange(0,right_end, 1000000)
    PL.xticks(xrange, ["%d"%(int(x/1000000)) for x in xrange])
    PL.xlabel("Genome (Mb)"), PL.ylabel("Reference allele frequency")
    PL.savefig(out_file, dpi=dpi)

def plotCoeff(X, y, obj, featureNames, whichReg):
    """ Plot Regression's Coeff
    """
    clf = classifiers[whichReg]
    clf,_,_ = fitAlgo(clf, X,y, opt= True, param_dict = param_dist_dict[whichReg])
    if whichReg == "LogisticRegression":
    	coeff = np.absolute(clf.coef_[0])
    else:
    	coeff = np.absolute(clf.coef_)
    print coeff
    indices = np.argsort(coeff)[::-1]
    print indices
    print featureNames
    featureList = []
    # num_features = len(featureNames)
    print("Feature ranking:")
    for f in range(num_features):
        featureList.append(featureNames[indices[f]])
        print("%d. feature %s (%.2f)" % (f, featureNames[indices[f]], coeff[indices[f]]))
    fig = pl.figure(figsize=(8,6),dpi=150)
    pl.title("Feature importances",fontsize=30)
    # pl.bar(range(num_features), coeff[indices],
    #         yerr = std_importance[indices], color=paired[0], align="center",
    #         edgecolor=paired[0],ecolor=paired[1])
    pl.bar(range(num_features), coeff[indices], color=paired[0], align="center",
            edgecolor=paired[0],ecolor=paired[1])
    pl.xticks(range(num_features), featureList, size=15,rotation=90)
    pl.ylabel("Importance",size=30)
    pl.yticks(size=20)
    pl.xlim([-1, num_features])
    # fix_axes()
    pl.tight_layout()
    save_path = 'plots/'+obj+'/'+whichReg+'_feature_importances.pdf'
    fig.savefig(save_path)

def display_coeff(data=None):
    betaAll,betaErrAll, R2adjAll = measure_stamp_coeff(data = data, zernike_max_order=20)
    ind = np.arange(len(betaAll[0]))
    momname = ('M20','M22.Real','M22.imag','M31.real','M31.imag','M33.real','M33.imag')
    fmtarr = ['bo-','ro-','go-','co-','mo-','yo-','ko-']
    pl.figure(figsize=(17,13))
    for i in range(7):
        pl.subplot(7,1,i+1)
        pl.errorbar(ind,betaAll[i],yerr = betaErrAll[i],fmt=fmtarr[i])
        pl.grid()
        pl.xlim(-1,21)
        if i ==0:
            pl.ylim(-10,65)
        elif i ==1:
            pl.ylim(-5,6)
        elif i ==2:
            pl.ylim(-5,6)
        elif i == 3:
            pl.ylim(-0.1,0.1)
        elif i == 4:
            pl.ylim(-0.1,0.1)
        elif i ==5:
            pl.ylim(-100,100)
        elif i == 6:
            pl.ylim(-100,100)
        pl.xticks(ind,('','','','','','','','','','','','','','','','','','','',''))
        pl.ylabel(momname[i])
    pl.xticks(ind,('0','1','2','3','4','5','6','7','8','9','10','11','12','13','14','15','16','17','18','19'))
    pl.xlabel('Zernike Coefficients')
    return '--- done ! ----'