def run_div_test(fld, exact, show=False, ignore_inexact=False):
    t0 = time()
    result_numexpr = viscid.div(fld, preferred="numexpr", only=True)
    t1 = time()
    logger.info("numexpr magnitude runtime: %g", t1 - t0)

    result_diff = viscid.diff(result_numexpr, exact[1:-1, 1:-1, 1:-1])
    if not ignore_inexact and not (result_diff.data < 5e-5).all():
        logger.warn("numexpr result is far from the exact result")
    logger.info("min/max(abs(numexpr - exact)): %g / %g",
                np.min(result_diff.data), np.max(result_diff.data))

    planes = ["y=0f", "z=0f"]
    nrows = 2
    ncols = len(planes)
    ax = plt.subplot2grid((nrows, ncols), (0, 0))
    ax.axis("equal")

    for i, p in enumerate(planes):
        plt.subplot2grid((nrows, ncols), (0, i), sharex=ax, sharey=ax)
        mpl.plot(result_numexpr, p, show=False)
        plt.subplot2grid((nrows, ncols), (1, i), sharex=ax, sharey=ax)
        mpl.plot(result_diff, p, show=False)

    if show:
        mpl.mplshow()

def fit_and_plot(cand, spd):
    data = cand.profile
    n = len(data)
    rms = np.std(data[(n/2):])
    xs = np.linspace(0.0, 1.0, n, endpoint=False)
    G = gauss._compute_data(cand)
    print "    Reduced chi-squared: %f" % (G.get_chisqr(data) / G.get_dof(n))
    print "    Baseline rms: %f" % rms
    print "    %s" % G.components[0]

    fig1 = plt.figure(figsize=(10,10))
    plt.subplots_adjust(wspace=0, hspace=0)

    # upper
    ax1 = plt.subplot2grid((3,1), (0,0), rowspan=2, colspan=1)
    ax1.plot(xs, data/rms, color="black", label="data")
    ax1.plot(xs, G.components[0].make_gaussian(n), color="red", label="best fit")

    # lower
    ax2 = plt.subplot2grid((3,1), (2,0), sharex=ax1)
    ax2.plot(xs, data/rms - G.components[0].make_gaussian(n), color="black", label="residuals")
    ax2.set_xlabel("Fraction of pulse window")

    plt.figure()
    plt.pcolormesh(xs, spd.waterfall_freq_axis(), spd.data_zerodm_dedisp, cmap=Greys)
    plt.xlabel("Fraction of pulse window")
    plt.ylabel("Frequency (MHz)")
    plt.xlim(0, 1)
    plt.ylim(spd.min_freq, spd.max_freq)

    plt.show()

    def __init__(self):

        self.fps = 20

        self.N = 10
        self.L = 15

        self.x=np.array([])
        self.y=np.array([])

        self.fig1 = plt.figure(figsize = (6,16))
        self.ax1 = plt.subplot2grid((3,2),(1,0),colspan=2,rowspan=2)
        self.ax1.set_xlim([0,self.L])
        self.ax1.set_ylim([0,self.L])
        self.ax1.set_xlabel('Simulation')

        self.ax2 = plt.subplot2grid((3,2),(0,0),colspan=2)
        self.ax2.set_xlim([0,10])
        self.ax2.set_ylim([0,1.1])
        self.ax2.set_xlabel('Time',labelpad=-15)
        self.ax2.set_ylabel('$\phi$',fontsize=20)
        self.line, = self.ax2.plot([],[])

        #To do the animation
        self.ani1 = animation.FuncAnimation(self.fig1,self.update_scatt,interval =self.fps ,init_func = self.setup_plot,blit=True)

def createExamples() :
    numbersWeHave = range (0, 1)
    versionsWeHave = range(1, 5)
    path = '/Users/ahmadbarakat/Downloads/CCEP Term 7/Pattern Recognition/ass1/ass1/data/'
    
    arr = []    
    for eachNum in numbersWeHave:
        for eachVer in versionsWeHave:
            imgFilePath = path + str(eachNum) + '-' + str(eachVer) + '.jpg'
            ei = Image.open(imgFilePath)
            eiar = np.asarray(ei)
            print str(eachNum) + '-' + str(eachVer)
            eiar = threshold(eiar)
            arr.append(eiar)

        fog = plt.figure()
        ax1 = plt.subplot2grid((2,2), (0,0))
        ax2 = plt.subplot2grid((2,2), (0,1))
        ax3 = plt.subplot2grid((2,2), (1,0))
        ax4 = plt.subplot2grid((2,2), (1,1))

        ax1.imshow(graphArr(arr[0]))
        ax2.imshow(graphArr(arr[1]))
        ax3.imshow(graphArr(arr[2]))
        ax4.imshow(graphArr(arr[3]))
        plt.show()

def fempipeWidget(alpha, pipedepth):
    respEW, respNS, X, Y = fempipe(alpha, pipedepth)
    fig = plt.figure(figsize = (8, 6))
    ax0 = plt.subplot2grid((2,2), (0,0))
    ax1 = plt.subplot2grid((2,2), (0,1))
    ax2 = plt.subplot2grid((2,2), (1,0), colspan=2)

    dat0 = ax0.imshow(respEW.real*100, extent=[X.min(),X.max(),Y.min(),Y.max()])
    dat1 = ax1.imshow(respNS.real*100, extent=[X.min(),X.max(),Y.min(),Y.max()])
    cb0 = plt.colorbar(dat0, ax = ax0)
    cb1 = plt.colorbar(dat1, ax = ax1)
    ax0.set_title("In-phase EW boom (%)", fontsize = 12)
    ax1.set_title("In-phase NS boom (%)", fontsize = 12)
    ax0.set_xlabel("Easting (m)", fontsize = 12)
    ax1.set_xlabel("Easting (m)", fontsize = 12)
    ax0.set_ylabel("Northing (m)", fontsize = 12)
    ax1.set_ylabel("Northing (m)", fontsize = 12)
    ax0.plot(np.r_[0., 0.], np.r_[-10., 10.], 'k--', lw=2)
    ax1.plot(np.r_[0., 0.], np.r_[-10., 10.], 'k--', lw=2)

    ax2.plot(Y[:,20],respEW[:, 20].real, 'k.-')
    ax2.plot(Y[:,20],respEW[:, 20].imag, 'k--')
    ax2.plot(Y[:,20],respNS[:, 20].real, 'r.-')
    ax2.plot(Y[:,20],respNS[:, 20].imag, 'r--')
    ax2.legend(('In-phase EW boom', 'Out-of-phase EW boom', 'In-phase NS boom', 'Out-of-phase NS boom'),loc=4)
    ax2.grid(True)
    ax2.set_ylabel('Hs/Hp (%)', fontsize = 16)
    ax2.set_xlabel('Northing (m)', fontsize = 16)
    ax2.set_title('Northing profile line at Easting 0 m', fontsize = 16)

    plt.tight_layout()
    plt.show()

def _onclick_help(event, params):
    """Function for drawing help window"""
    import matplotlib.pyplot as plt
    text, text2 = _get_help_text(params)

    width = 6
    height = 5

    fig_help = figure_nobar(figsize=(width, height), dpi=80)
    fig_help.canvas.set_window_title('Help')
    ax = plt.subplot2grid((8, 5), (0, 0), colspan=5)
    ax.set_title('Keyboard shortcuts')
    plt.axis('off')
    ax1 = plt.subplot2grid((8, 5), (1, 0), rowspan=7, colspan=2)
    ax1.set_yticklabels(list())
    plt.text(0.99, 1, text, fontname='STIXGeneral', va='top', weight='bold',
             ha='right')
    plt.axis('off')

    ax2 = plt.subplot2grid((8, 5), (1, 2), rowspan=7, colspan=3)
    ax2.set_yticklabels(list())
    plt.text(0, 1, text2, fontname='STIXGeneral', va='top')
    plt.axis('off')

    tight_layout(fig=fig_help)
    # this should work for non-test cases
    try:
        fig_help.canvas.draw()
        fig_help.show(warn=False)
    except Exception:
        pass

def plot_scaling():
    X, y = make_blobs(n_samples=50, centers=2, random_state=4, cluster_std=1)
    X += 3

    plt.figure(figsize=(15, 8))
    main_ax = plt.subplot2grid((2, 4), (0, 0), rowspan=2, colspan=2)

    main_ax.scatter(X[:, 0], X[:, 1], c=y, cmap=cm2, s=60)
    maxx = np.abs(X[:, 0]).max()
    maxy = np.abs(X[:, 1]).max()

    main_ax.set_xlim(-maxx + 1, maxx + 1)
    main_ax.set_ylim(-maxy + 1, maxy + 1)
    main_ax.set_title("Original Data")
    other_axes = [plt.subplot2grid((2, 4), (i, j)) for j in range(2, 4) for i in range(2)]

    for ax, scaler in zip(other_axes, [StandardScaler(), RobustScaler(),
                                       MinMaxScaler(), Normalizer(norm='l2')]):
        X_ = scaler.fit_transform(X)
        ax.scatter(X_[:, 0], X_[:, 1], c=y, cmap=cm2, s=60)
        ax.set_xlim(-2, 2)
        ax.set_ylim(-2, 2)
        ax.set_title(type(scaler).__name__)

    other_axes.append(main_ax)

    for ax in other_axes:
        ax.spines['left'].set_position('center')
        ax.spines['right'].set_color('none')
        ax.spines['bottom'].set_position('center')
        ax.spines['top'].set_color('none')
        ax.xaxis.set_ticks_position('bottom')
        ax.yaxis.set_ticks_position('left')

def visualize_tss(flare_class, prediction_threshold, flare_threshold,tss):
    hits_x = []
    hits_y = []
    miss_x = []
    miss_y = []

    for x,y in plot_xy:
        if (x>=prediction_threshold) == (y>=flare_threshold):
            hits_x.append(x)
            hits_y.append(y)
        else:
            miss_x.append(x)
            miss_y.append(y)

    plt.rcParams['figure.figsize'] = (12.8,9.6)
    plt.subplot2grid((1,1),(0,0), colspan=1, rowspan=1)
    plt.plot(miss_x, miss_y, 'mo',markersize=1.0, markeredgecolor='r')
    plt.plot(hits_x, hits_y, 'mo',markersize=1.0, markeredgecolor='b')
    plt.gca().set_xscale('log')
    plt.gca().set_yscale('log')
    plt.gca().set_xlabel("AIA 193nm thresholded sum ({})".format(threshold_value))
    plt.gca().set_ylabel("GOES 1-8A 24hour future max")
    filename = "Flarepredict-{}-{}.png".format(flare_class, threshold_value)
    plt.title("{}-class flare prediction with TI({}) : TSS = {}".format(flare_class, threshold_value, tss))

    plt.savefig(filename, dpi=100)
    plt.close("all")

def plot_calibration_curve(est, name, fig_index):
    """Plot calibration curve for est w/o and with calibration. """
    # Calibrated with isotonic calibration
    isotonic = CalibratedClassifierCV(est, cv=2, method='isotonic')

    # Calibrated with sigmoid calibration
    sigmoid = CalibratedClassifierCV(est, cv=2, method='sigmoid')

    # Calibrated with ROC convex hull calibration
    rocch = CalibratedClassifierCV(est, cv=2, method='rocch')

    # Logistic regression with no calibration as baseline
    lr = LogisticRegression(C=1., solver='lbfgs')

    fig = plt.figure(fig_index, figsize=(10, 10))
    ax1 = plt.subplot2grid((3, 1), (0, 0), rowspan=2)
    ax2 = plt.subplot2grid((3, 1), (2, 0))

    ax1.plot([0, 1], [0, 1], "k:", label="Perfectly calibrated")
    for clf, name in [(lr, 'Logistic'),
                      (est, name),
                      (isotonic, name + ' + Isotonic'),
                      (sigmoid, name + ' + Sigmoid'),
                      (rocch, name + ' + ROCConvexHull')]:
        clf.fit(X_train, y_train)
        y_pred = clf.predict(X_test)
        if hasattr(clf, "predict_proba"):
            prob_pos = clf.predict_proba(X_test)[:, 1]
        else:  # use decision function
            prob_pos = clf.decision_function(X_test)
            prob_pos = \
                (prob_pos - prob_pos.min()) / (prob_pos.max() - prob_pos.min())

        clf_score = brier_score_loss(y_test, prob_pos, pos_label=y.max())
        print("%s:" % name)
        print("\tBrier: %1.4f" % (clf_score))
        print("\tPrecision: %1.3f" % precision_score(y_test, y_pred))
        print("\tRecall: %1.3f" % recall_score(y_test, y_pred))
        print("\tF1: %1.3f" % f1_score(y_test, y_pred))
        print("\tAuc: %1.4f\n" % roc_auc_score(y_test, prob_pos))

        fraction_of_positives, mean_predicted_value = \
            calibration_curve(y_test, prob_pos, n_bins=10)

        ax1.plot(mean_predicted_value, fraction_of_positives, "s-",
                 label="%s (%1.4f)" % (name, clf_score))

        ax2.hist(prob_pos, range=(0, 1), bins=10, label=name,
                 histtype="step", lw=2)

    ax1.set_ylabel("Fraction of positives")
    ax1.set_ylim([-0.05, 1.05])
    ax1.legend(loc="lower right")
    ax1.set_title('Calibration plots  (reliability curve)')

    ax2.set_xlabel("Mean predicted value")
    ax2.set_ylabel("Count")
    ax2.legend(loc="upper center", ncol=2)

    plt.tight_layout()

    def plot(self, weights=True, assets=True, portfolio_label='PORTFOLIO', **kwargs):
        """ Plot equity of all assets plus our strategy.
        :param weights: Plot weights as a subplot.
        :param assets: Plot asset prices.
        :return: List of axes.
        """
        res = ListResult([self], [portfolio_label])
        if not weights:
            ax1 = res.plot(assets=assets, **kwargs)
            return [ax1]
        else:
            plt.figure(1)
            ax1 = plt.subplot2grid((3, 1), (0, 0), rowspan=2)
            res.plot(assets=assets, ax=ax1, **kwargs)
            ax2 = plt.subplot2grid((3, 1), (2, 0), sharex=ax1)

            # plot weights as lines
            if self.B.values.min() < -0.01:
                self.B.plot(ax=ax2, ylim=(min(0., self.B.values.min()), max(1., self.B.sum(1).max())),
                            legend=False, colormap=plt.get_cmap('jet'))
            else:
                # fix rounding errors near zero
                if self.B.values.min() < 0:
                    B = self.B - self.B.values.min()
                else:
                    B = self.B
                B.plot(ax=ax2, ylim=(0., max(1., B.sum(1).max())),
                       legend=False, colormap=plt.get_cmap('jet'), kind='area', stacked=True)
            plt.ylabel('weights')
            return [ax1, ax2]

def calibration_plot(clf, xtest, ytest):
	prob = clf.predict_proba(xtest)[:, 1]
	outcome = ytest
	data = pd.DataFrame(dict(prob=prob, outcome=outcome))

	# group outcomes into bins of similar probability
	bins = np.linspace(0, 1, 20)
	cuts = pd.cut(prob, bins)
	binwidth = bins[1] - bins[0]

	# freshness ratio and number of examples in each bin
	cal = data.groupby(cuts).outcome.agg(['mean', 'count'])
	cal['pmid'] = (bins[:-1] + bins[1:]) / 2
	cal['sig'] = np.sqrt(cal.pmid * ( 1- cal.pmid) / cal['count'])

	# the calibration plot
	ax = plt.subplot2grid((3, 1), (0, 0), rowspan=2)
	p = plt.errorbar(cal.pmid, cal['mean'], cal['sig'])
	plt.plot(cal.pmid, cal.pmid, linestyle='--', lw=l, color='k')
	plt.ylabel("Empirical P(Fresh)")
	remove_border(ax)

	# the distirubtion of P(fresh)
	ax = plt.subplot2grid((3, 1), (2, 0), sharex=ax)

	plt.bar(left=cal.pmid - binwidth / 2, height=cal['count'],
			width=.95 * (bins[1] - bins[0]),
			fc=p[0].get_color())

	plt.xlabel("Predicted P(Fresh)")
	remove_border()
	plt.ylabel("Number")

def main():

  print "loading data.."
  # Open hdf5 file
  df = pd.read_csv('train.shuffled', nrows = 100);
  df = df.dropna();
  # specifies the parameters of our graphs
  fig = plt.figure(figsize=(18,6)); 
  alpha=alpha_scatterplot = 0.2 
  alpha_bar_chart = 0.55
   # lets us plot many diffrent shaped graphs together 
  ax1 = plt.subplot2grid((2,3),(0,0))
# plots a bar graph of those who surived vs those who did not.               
  df.click.value_counts().plot(kind='bar', alpha=alpha_bar_chart)
  ax1.set_xlim(-1, 2)
# puts a title on our graph
  plt.title("Distribution of click, (1 = clicked)")    

  plt.subplot2grid((2,3),(0,1))
  plt.scatter(df.click, df.banner_pos, alpha=alpha_scatterplot)
# sets the y axis lable
  plt.ylabel("banner pos")
# formats the grid line style of our graphs                          
  plt.grid(b=True, which='major', axis='y')  
  plt.title("Survial by banner pos,  (1 = clicked)")

  df["date_time"] = df['hour'].apply(string_to_date);
  dfts = df.set_index("date_time");
  time_mean =dfts.groupby(lambda x: x.time()).mean();
  print "saving data to hdf5"
  print time_mean
  test = pd.read_csv('test');

def paint_clustering(results, clusters, num, chrom, tad_names):
    dendros = []
    axes = []
    prev = 0
    xlim = [-100, 100]
    tmp = []
    for i, result in enumerate(results):
        if axes:
            axes[-1].set_xticklabels([], visible=False)
        clust = linkage(result, method='ward')
        tmp = dendrogram(clust, orientation='right', no_plot=True)['leaves']
        dendros += reversed(list([clusters[i][n] for n in tmp]))
        axes.append(plt.subplot2grid((num, 9),(prev, 0), rowspan=len(result),
                                     colspan=4))
        dendrogram(clust, orientation='right',
                   labels=[tad_names[c] for c in clusters[i]])
        if xlim[0] < axes[-1].get_xlim()[0]:
            xlim[0] = axes[-1].get_xlim()[0]
        if xlim[1] > axes[-1].get_xlim()[1]:
            xlim[1] = axes[-1].get_xlim()[1]
        prev += len(result)
    for ax in axes:
        ax.set_xlim(left=xlim[0], right=xlim[1])
    axes = []
    for i, j in enumerate(dendros):
        axes.append(plt.subplot2grid((num, 9),(i, 4)))#gs1[i]))
        chrom.visualize('exp1',
                        tad=chrom.get_experiment('exp1').tads[tad_names[j]],
                        axe=axes[-1], show=False)
        axes[-1].set_axis_off()
    ax4 = plt.subplot2grid((num, 9),(0, 5), rowspan=num, colspan=4)
    chrom.visualize('exp1', paint_tads=True, axe=ax4)
    plt.draw()

def plot_reg_hist(datafile):
    AVE_A = 0.42128425
    STD_A = 0.28714299
    AVE_Z = 0.09973424
    STD_Z = 0.05802179

    ID, real, res = ReadData(datafile)

    uAn = res[:, 0]
    rAn = real[:, 0]
    uZn = res[:, 1]
    rZn = real[:, 1]

    dAn = uAn - rAn
    dZn = uZn - rZn

    uA = 3.0 * STD_A * uAn + AVE_A
    uZ = 3.0 * STD_Z * uZn + AVE_Z

    rA = 3.0 * STD_A * rAn + AVE_A
    rZ = 3.0 * STD_Z * rZn + AVE_Z

    dA = uA - rA
    dZ = uZ - rZ

    fig = plt.figure(figsize=(17,8))
    ax1 = plt.subplot2grid((1,2), (0,0))
    ax2 = plt.subplot2grid((1,2), (0,1))
    ax1.set_title('Extinction SD Hist')
    ax1.hist(dA)
    ax2.set_title('Redshift SD Hist')
    ax2.hist(dZ)

def overlap_plot(profile_A, profile_B, w, y_slice):
    """Make nice overlap plots.
    """
    import matplotlib.pyplot as plt
    plt.figure()
    plt.subplots_adjust(hspace=0.001)
    
    norm_product_AB = profile_A * profile_B 
    ax1 = plt.subplot2grid((4, 1), (0, 0), rowspan=3)
    ax1.plot(profile_A, label='Source A', linewidth=4, color='b')
    ax1.plot(profile_B, label='Source B', linewidth=4, color='g')
    ax1.tick_params(axis='both', which='major', labelsize=16)
    plt.yticks(np.arange(0, 1.2, 0.2))
    plt.legend()
    
    ax2 = plt.subplot2grid((4, 1), (3, 0), rowspan=1)
    ax2.plot(norm_product_AB, label='Product', linewidth=4, color='r')
    plt.yticks(np.arange(0, 0.5, 0.25))
    plt.ylim(0, 0.5)
    ax2.tick_params(axis='both', which='major', labelsize=16)
    
    pos = np.linspace(0, profile_A.size, 11)
    y_labels, x_labels = w.wcs_pix2world(pos, np.ones_like(pos) * y_slice, 0)
    plt.xticks(pos, x_labels)
    plt.setp(ax1.get_xticklabels(), visible=False)
    plt.legend()    

def viewQueryGoalProportions(distribution_1, distribution_2):
    
    proportion_1, goals_1, max_goal_1 = infoDistribution(distribution_1)
    proportion_2, goals_2, max_goal_2 = infoDistribution(distribution_2)

    width = 0.5 # gives histogram aspect to the bar diagram

    plt.figure(figsize=(16, 5))
    ax1 = plt.subplot2grid((1, 2), (0, 0))
    #plt.subplot(1,2,1)
    #ax = plt.axes()
    #ax.set_xticks(pos_1 + (width / 2))
    #ax.set_xticklabels(goals_1)
    ax1.bar(goals_1, proportion_1, width, color='r')
    plt.ylabel('Goals distribution for query')
    plt.xlabel('Goals Distribution 1')   
    plt.grid(True)  
    
    ax2= plt.subplot2grid((1, 2), (0, 1))
    ax2.bar(goals_2, proportion_2, width, color='blue')
    plt.ylabel('Goals distribution for query')
    plt.xlabel('Goals Distribution 2')        
    
    plt.grid(True)    
    plt.show()
    
    max_goals = []
    max_goals.append(max_goal_1)
    max_goals.append(max_goal_2)

    return max_goals

def test(ym):
    
    runtime = 1000.*100

    V, t = ym.simulate(runtime)
    df = ym.sim_data
    print "Sim done"

    Q = df.loc[0.0:,'Q10':'Q40']
    C = df.loc[0.0:,'C10':'C40']
    
    Qt = Q.T
    Ct = C.T
    
    fig = plt.figure()
    '''
    ax1 = plt.subplot(211)
    ax2 = plt.subplot(222)
    ax3 = plt.subplot(223)
    '''
    ax1 = plt.subplot2grid((2,2), (0,0), colspan=2)
    ax2 = plt.subplot2grid((2,2), (1,0))
    ax3 = plt.subplot2grid((2,2), (1, 1))


    ax1.plot(t,V)
    
    plot_states(ax2,Qt,t)
    plot_states(ax3,Ct,t)

    plt.tight_layout()
    return plt.gca()
    
    '''

    def plotDistanceInfo(self, title=None, printInfo=False):
        """
        Utility function for generating informative plots of the ensemble relative to the observation
        """
        if self.observation_type == dautils.ObservationType.UnderlyingFlow or \
           self.observation_type == dautils.ObservationType.DirectUnderlyingFlow:
            return self.plotVelocityInfo(title=title, printInfo=printInfo)

        fig = plt.figure(figsize=(10,6))
        gridspec.GridSpec(2, 3)

        # PLOT POSITIONS OF PARTICLES AND OBSERVATIONS
        ax0 = plt.subplot2grid((2,3), (0,0))
        plt.plot(self.observeParticles()[:,:,0].flatten(), \
                 self.observeParticles()[:,:,1].flatten(), 'b.')
        plt.plot(self.observeTrueState()[:,0], \
                 self.observeTrueState()[:,1], 'r.')
        ensembleMean = self.getEnsembleMean()
        if ensembleMean is not None:
            plt.plot(ensembleMean[0], ensembleMean[1], 'r+')
        plt.xlim(0, self.getDomainSizeX())
        plt.xlabel('x')
        plt.ylabel('y')
        plt.ylim(0, self.getDomainSizeY())
        plt.title("Particle positions")

        # PLOT DISCTRIBUTION OF PARTICLE DISTANCES AND THEORETIC OBSERVATION PDF
        ax0 = plt.subplot2grid((2,3), (0,1), colspan=2)
        innovations = self.getInnovationNorms()
        obs_var = self.getObservationVariance()
        plt.hist(innovations, bins=30, \
                 range=(0, max(min(self.getDomainSizeX(), self.getDomainSizeY()), np.max(innovations))),\
                 normed=True, label="particle innovations")

        # With observation 
        x = np.linspace(0, max(self.getDomainSizeX(), self.getDomainSizeY()), num=100)
        gauss_pdf = self.getGaussianWeight(x, normalize=False)
        plt.plot(x, gauss_pdf, 'g', label="pdf directly from innovations")
        plt.legend()
        plt.title("Distribution of particle innovations")

        # PLOT SORTED DISTANCES FROM OBSERVATION
        ax0 = plt.subplot2grid((2,3), (1,0), colspan=3)
        gaussWeights = self.getGaussianWeight()
        indices_sorted_by_observation = innovations.argsort()
        ax0.plot(gaussWeights[indices_sorted_by_observation]/np.max(gaussWeights), 'g', label="Gauss weight")
        ax0.set_ylabel('Relative weight')
        ax0.grid()
        ax0.set_ylim(0,1.4)
        plt.legend(loc=7)

        ax1 = ax0.twinx()
        ax1.plot(innovations[indices_sorted_by_observation], label="innovations")
        ax1.set_ylabel('Innovations from observation', color='b')

        plt.title("Sorted innovations from observation")

        if title is not None:
            plt.suptitle(title, fontsize=16)
        return fig

def calibration_curve_plotter(y_test, prob_pos, n_bins=10):

    brier = brier_score_loss(y_test, prob_pos, pos_label=1)

    fig = plt.figure(0, figsize=(10, 10))
    ax1 = plt.subplot2grid((3, 1), (0, 0), rowspan=2)
    ax2 = plt.subplot2grid((3, 1), (2, 0))

    df = pd.DataFrame({"true": y_test})
    bins = np.linspace(0.0, 1.0, n_bins + 1)
    binids = np.digitize(prob_pos, bins) - 1
    df["Bin center"] = bins[binids] + 0.5 / n_bins
    df[""] = "Model calibration: (%1.5f)" % brier
    o = bins + 0.5 / n_bins

    df2 = pd.DataFrame({"true": o, "Bin center": o})
    df2[""] = "Perfect calibration"

    df = pd.concat([df, df2])

    sns.pointplot(x="Bin center", y="true", data=df, order=o, hue="", ax=ax1)

    ax2.hist(prob_pos, range=(0, 1), bins=10, label="Model", histtype="step", lw=2)

    ax1.set_ylabel("Fraction of positives")
    ax1.set_ylim([-0.05, 1.05])
    # ax1.legend(loc="lower right")
    ax1.set_title("Calibration plots")

    ax2.set_xlabel("Predicted Probability")
    ax2.set_ylabel("Count")

    plt.tight_layout()

def plot_all_3(results, colors):
  x = np.arange(len(results))

  m_results = tuple(r[1] for r in results)
  m_results = tuple(zip(*m_results)) # transpose
  m_results = tuple((m,r) for m, r in zip(methods, m_results))
  print(m_results)

  ax1 = plt.subplot2grid((3, 5), (0, 1), colspan=4)
  ax2 = plt.subplot2grid((3, 5), (1, 1), colspan=4)
  ax3 = plt.subplot2grid((3, 5), (2, 1), colspan=4)
  ax1.set_title('all')
  ax2.set_ylabel('time / ms')
  ax2.set_ylim(0,m_results[-3][1][-1])
  ax3.set_ylim(0,0.05)

  ax1.set_xticks(x)
  ax1.set_xticklabels(['' for _ in range(5)])
  ax2.set_xticks(x)
  ax2.set_xticklabels(['' for _ in range(5)])
  ax3.set_xticks(x)
  ax3.set_xticklabels([r[0] for r in results])

  for idx, m_result in enumerate(m_results):
    title, times = m_result
    ax1.plot(x, times, color=colors[idx])
    ax2.plot(x, times, color=colors[idx])
    ax3.plot(x, times, color=colors[idx])
    
    # ax2.barh(x[:-2], times[:-2], color=colors[:-2])
    # plt.savefig(title+'.png')
  plt.show()