def plot_likelihood_overspace(info, session, task_labels, colours, filepath=None):
    for task_label in task_labels:
        zones = getattr(session, task_label).zones
        likelihood = np.nanmean(np.array(getattr(session, task_label).likelihoods), axis=(0, 1))

        likelihood[np.isnan(likelihood)] = 0

        xx, yy = np.meshgrid(info.xedges, info.yedges)
        xcenters, ycenters = get_bin_centers(info)
        xxx, yyy = np.meshgrid(xcenters, ycenters)

        maze_highlight = "#fed976"
        plt.plot(session.position.x, session.position.y, ".", color=maze_highlight, ms=1, alpha=0.2)
        pp = plt.pcolormesh(xx, yy, likelihood, cmap='bone_r')
        for label in ["u", "shortcut", "novel"]:
            plt.contour(xxx, yyy, zones[label], levels=0, linewidths=2, colors=colours[label])
        plt.colorbar(pp)
        plt.axis('off')

        plt.tight_layout()
        if filepath is not None:
            filename = info.session_id + "_" + task_label + "_likelihoods-overspace.png"
            plt.savefig(os.path.join(filepath, filename))
            plt.close()
        else:
            plt.show()

    def region_images(self, save_name=None, percentile=80.):
      '''

      Creates saved PDF plots of several quantities/images.

      '''

      if self.verbose:
        pass
      else:
          threshold = scoreatpercentile(self.image[~np.isnan(self.image)], percentile)
          p.imshow(self.image, vmax=threshold, origin="lower", interpolation="nearest")
          p.contour(self.mask)
          p.title("".join([save_name," Contours at ", str(round(threshold))]))
          p.savefig("".join([save_name,"_filaments.pdf"]))
          p.close()

          ## Skeletons
          masked_image = self.image * self.mask
          skel_points = np.where(self.skeleton==1)
          for i in range(len(skel_points[0])):
              masked_image[skel_points[0][i],skel_points[1][i]] = np.NaN
          p.imshow(masked_image, vmax=threshold, interpolation=None, origin="lower")
          p.savefig("".join([save_name,"_skeletons.pdf"]))
          p.close()

      return self

def plot_haxby(activation, title):
    z = 25

    fig = plt.figure(figsize=(4, 5.4))
    fig.subplots_adjust(bottom=0., top=1., left=0., right=1.)
    plt.axis('off')
    # pl.title('SVM vectors')
    plt.imshow(mean_img[:, 4:58, z].T, cmap=pl.cm.gray,
              interpolation='nearest', origin='lower')
    plt.imshow(activation[:, 4:58, z].T, cmap=pl.cm.hot,
              interpolation='nearest', origin='lower')

    mask_house = nib.load(h.mask_house[0]).get_data()
    mask_face = nib.load(h.mask_face[0]).get_data()

    plt.contour(mask_house[:, 4:58, z].astype(np.bool).T, contours=1,
            antialiased=False, linewidths=4., levels=[0],
            interpolation='nearest', colors=['blue'], origin='lower')

    plt.contour(mask_face[:, 4:58, z].astype(np.bool).T, contours=1,
            antialiased=False, linewidths=4., levels=[0],
            interpolation='nearest', colors=['limegreen'], origin='lower')

    p_h = Rectangle((0, 0), 1, 1, fc="blue")
    p_f = Rectangle((0, 0), 1, 1, fc="limegreen")
    plt.legend([p_h, p_f], ["house", "face"])
    plt.title(title, x=.05, ha='left', y=.90, color='w', size=28)

def VisualizeFit(Xval, pval, epsilon, mu, sigma):
    """
    Visualize the fitter data
    :param Xval: the validation data set (only the first two columns are used)
    :param pval: A vector containing probabilities for example data in Xval
    :param mu: Estimate for the mean, using the training data
    :param sigma: Estimate for the variance, using the training data
    :return:
    """
    np.seterr(over='ignore')
    x1, x2 = np.meshgrid(np.arange(0, 35, 0.5), np.arange(0, 35, 0.5))
    z1 = np.asarray(x1).reshape(-1)

    mat = np.zeros([len(z1), 2])
    mat[:, 0] = np.asarray(x1).reshape(-1)
    mat[:, 1] = np.asarray(x2).reshape(-1)

    Z = MultivariateGaussian(mat, mu, sigma)
    Z = np.reshape(Z, np.shape(x1))

    x = [10 ** x for x in np.arange(-20, 0, 3)]

    plt.figure(1)

    plt.scatter(Xval[:, 0], Xval[:, 1], c=None, s=25, alpha=None, marker="+")

    points = np.where(pval < epsilon)
    plt.scatter(Xval[:, 0][points], Xval[:, 1][points], s=50, marker='+', color='red')
    plt.contour(x1, x2, Z, x)

    plt.show()

def steepVsConj():
    Q = np.array([[1.,0], [0,10.]])
    b = np.array([0.,0.])
    def f(pt):
        return .5*(pt*Q.dot(pt)).sum() - (b*pt).sum()
    #pts = [np.array([1.5,1.5]), np.array([1.,1.]), np.array([.5,.5])]
    x0 = np.array([5.,.5])
    pts = steepestDescent(Q, b, x0, niter=20)
    Q = np.array([[1.,0], [0,10.]])
    b = np.array([0.,0.])
    x0 = np.array([5.,.5])
    pts2 = conjugateGradient(Q, b, x0)
    dom = np.linspace(-5, 5, 100)
    X, Y = np.meshgrid(dom, dom)
    Z = np.empty(X.shape)
    for i in xrange(X.shape[0]):
        for j in xrange(X.shape[1]):
            pt = np.array([X[i,j], Y[i,j]])
            Z[i,j] = f(pt)
    vals = np.empty(len(pts))
    for i in xrange(len(pts)):
        vals[i] = f(pts[i])
    plt.contour(X,Y,Z, vals[:5], colors='gray')
    plt.plot(np.array(pts)[:,0],np.array(pts)[:,1], '*-')
    plt.plot(np.array(pts2)[:,0],np.array(pts2)[:,1], '*-')
    plt.savefig('steepVsConj.pdf')
    plt.clf()

def PlotDecisionBoundary(predictor, train_data, test_data):
  x1 = np.arange(0, 1, 0.01)
  x2 = np.arange(0, 1, 0.01)
  x1,x2 = np.meshgrid(x1, x2)
  y = np.zeros_like(x1)
  m,n = y.shape
  for i in range(m):
    for j in range(n):
      y[i,j] = predictor.Classify([x1[i,j], x2[i,j]])
  # now do plots
  train_pos = train_data[:,2] == 1
  train_neg = train_data[:,2] == -1
  test_pos = test_data[:,2] == 1
  test_neg = test_data[:,2] == -1
  plt.figure(figsize=(12,5))
  plt.subplot(121)
  plt.plot(train_data[train_pos,0], train_data[train_pos,1], 'r+')
  plt.plot(train_data[train_neg,0], train_data[train_neg,1], 'b+')
  plt.contour(x1, x2, y, levels=[0])
  plt.xlim(0, 1)
  plt.xlabel(r'$x_1$', fontsize=20)
  plt.ylim(0, 1)
  plt.ylabel(r'$x_2$', fontsize=20)
  plt.title('Decision boundary overimposed on training data', fontsize=15)
  plt.subplot(122)
  plt.plot(test_data[test_pos,0], test_data[test_pos,1], 'r+')
  plt.plot(test_data[test_neg,0], test_data[test_neg,1], 'b+')
  plt.contour(x1, x2, y, levels=[0])
  plt.xlim(0, 1)
  plt.xlabel(r'$x_1$', fontsize=20)
  plt.ylim(0, 1)
  plt.ylabel(r'$x_2$', fontsize=20)
  plt.title('Decision boundary overimposed on test data', fontsize=15)
  plt.tight_layout()
  plt.show()

def visualCost(X,y,theta):
    theta0_vals = np.linspace(-10, 10, 500)
    theta1_vals = np.linspace(-1, 4, 500)
    J_vals = np.zeros((len(theta0_vals),len(theta1_vals)))
    for i, elem0 in enumerate(theta0_vals):
        for j, elem1 in enumerate(theta1_vals):
            theta_vals = np.array([[elem0], [elem1]])
            J_vals[i, j] = computeCost(X, y, theta_vals)       
    theta0_vals, theta1_vals = np.meshgrid(theta0_vals,theta1_vals)

    J_vals = J_vals.T
    #surface plot
    fig4 = plt.figure()
    ax = fig4.gca(projection='3d')
    ax.plot_surface(theta0_vals, theta1_vals, J_vals, cmap=cm.jet)
    ax.view_init(azim = 180+40,elev = 25)
    plt.xlabel("theta_0")
    plt.ylabel("theta_1")
    plt.savefig("cost3D.pdf")
    
    #contour plot
    plt.figure(5)
    plt.contour(theta0_vals, theta1_vals, J_vals, levels=np.logspace(-2, 3, 20))
    plt.scatter(theta[0,0], theta[1,0],  c='r', marker='x', s=20, linewidths=1)  # mark converged minimum
    plt.xlim([-10, 10])
    plt.ylim([-1, 4])
    plt.xlabel("theta_0")
    plt.ylabel("theta_1")
    plt.savefig("costContour.pdf")

def main(argv):
    # Load Station Data
    conn = sqlite3.connect(_BABS_DATABASE_PATH)
    station_data = pd.read_sql(_BABS_QUERY, conn)
    station_data['Usage'] = 0.5*(station_data['Start Count'] + station_data['End Count'])

    print("\n\nLoading model to file ... ")
    with open(os.path.join(os.path.dirname(os.path.realpath(__file__)), os.path.pardir, 'models', 'babs_usage_model.dill'), 'r') as f:
            babs_usage_model = dill.load(f)

    # Interpolate all features from feature models
    lats, longs = np.meshgrid(np.linspace(37.7,37.82,30), np.linspace(-122.55,-122.35,30))
    transformed_features = babs_usage_model.named_steps['features_from_lat_long'].transform(pd.DataFrame({'Latitude': lats.reshape(1,-1).squeeze(), 'Longitude': longs.reshape(1,-1).squeeze()}))

    prediction_features = pd.DataFrame({'Latitude': lats.reshape(1,-1).squeeze(), 'Longitude': longs.reshape(1,-1).squeeze()})
    usage_predictions = babs_usage_model.predict(prediction_features)
    usage_predictions[np.array(transformed_features['Elevation']<0)] = np.nan
    usage_predictions = np.lib.scimath.logn(100, usage_predictions - np.nanmin(usage_predictions) + 1)
    usage_predictions[np.where(np.isnan(usage_predictions))] = 0

    plt.contourf(longs, lats, usage_predictions.reshape(30,-1),
                 norm=colors.Normalize(np.mean(usage_predictions)-(1*np.std(usage_predictions)), np.mean(usage_predictions)+(1*np.std(usage_predictions)), clip=True),
                 levels=np.linspace(0.01,max(usage_predictions),300))
    plt.contour(longs, lats, (transformed_features['Elevation']).reshape(30,-1), linewidth=0.2, colors='white')
    plt.scatter(station_data[station_data['Landmark']=='San Francisco']['Longitude'], station_data[station_data['Landmark']=='San Francisco']['Latitude'], s=2, )

    # plt.scatter(longs,
    #             lats,
    #             #s=(usage_predictions<0)*10,
    #             s=(transformed_features['Elevation']>0)*10,
    #             cmap=matplotlib.cm.Reds)
    plt.show()

def visualize_fit(X, mu_vec, var_vec):
    """ Plots dataset and estimated Gaussian distribution.

    Args:
      X: Matrix of features.
      mu_vec: Vector of mean values of features.
      var_vec: Vector of variances of features.

    Returns:
      None.
    """
    pyplot.scatter(X[:, 0], X[:, 1], s=80, marker='x', color='b')
    pyplot.ylabel('Throughput (mb/s)', fontsize=18)
    pyplot.xlabel('Latency (ms)', fontsize=18)
    pyplot.hold(True)
    u_vals = numpy.linspace(0, 35, num=71)
    v_vals = numpy.linspace(0, 35, num=71)
    z_vals = numpy.zeros((u_vals.shape[0], v_vals.shape[0]))
    for u_index in range(0, u_vals.shape[0]):
        for v_index in range(0, v_vals.shape[0]):
            z_vals[u_index, v_index] = (
                multivariate_gaussian(numpy.c_[u_vals[u_index],
                                               v_vals[v_index]], mu_vec,
                                      var_vec))
    z_vals_trans = numpy.transpose(z_vals)
    u_vals_x, v_vals_y = numpy.meshgrid(u_vals, v_vals)
    z_vals_reshape = z_vals_trans.reshape(u_vals_x.shape)
    exp_seq = numpy.linspace(-20, 1, num=8)
    pow_exp_seq = numpy.power(10, exp_seq)
    pyplot.contour(u_vals, v_vals, z_vals_reshape, pow_exp_seq)
    pyplot.hold(False)
    return None

def plot2(X,Y,beta):
	#function to plot non linear decison boundary
	for i in range(1,X.shape[0]):
		if Y[i] == 1:
			plt.plot(X[i][1],X[i][2],'rs')
		else:
			plt.plot(X[i][1],X[i][2],'bs')
	#Here is the grid range
	#50 values btw -1 and 1.5
    	u = np.linspace(-1, 1.5, 50);
    	v = np.linspace(-1, 1.5, 50);

    	z = np.zeros((len(u), len(v)));
    	#Evaluate z = theta*x over the grid
    	#decision boundary is given by theta*x = 0
    	for i in range(0,len(u)):
    		for j in range(1,len(v)):
    			temp = np.matrix([u[i],v[j]])
    	        	z[i,j] = mapFeatures(temp,6)*beta;
    	 
    	z = z.T; # important to transpose z before calling contour
	
	#Plot z = 0
	#Notice you need to specify the range [0, 0]
	plt.contour(u, v, z, [0, 0])
	
	plt.show()	

def show_segmented_result(original_img, closed_heathly, closed_cancer):
    # We load with opencv and plot with pyplot (BGR->RGB)
    original_img = cv2.cvtColor(original_img, cv2.COLOR_BGR2RGB)
    plt.imshow(original_img)
    plt.contour(closed_heathly, [0.5], linewidths=0.7, colors='g')
    plt.contour(closed_cancer, [0.5], linewidths=0.7, colors='r')
    plt.savefig(sys.argv[1]+"_out.png", dpi=700)

def gaussFit(data, inix=0, iniy=0, amp=1.):
   
    xv = np.linspace(np.max(data.xaxis), np.min(data.xaxis), data.xaxis.shape[0])
    yv = np.linspace(np.min(data.yaxis), np.max(data.yaxis), data.yaxis.shape[0])

    x, y = np.meshgrid(xv, yv)
    
    p_init = models.Gaussian2D(amplitude=amp, x_mean=inix, y_mean=iniy, x_stddev=.01 , y_stddev=.1,theta=0.)
    fit_p = fitting.LevMarLSQFitter()
    
#    for i in range(256):
#        for j in range(256):
#            if (data.yaxis[j]>-51. ): #12*data.xaxis[i]-373.8):
#                data.image[j,i]=0.
    
    
    p = fit_p(p_init, x, y, data.image)
    print p
    
    th=p.theta.value
    a=(np.cos(th)**2)/(2*p.x_stddev.value**2) + (np.sin(th)**2)/(2*p.y_stddev.value**2)
    b=-(np.sin(2*th))/(4*p.x_stddev.value**2) + (np.sin(2*th)) /(4*p.y_stddev.value**2)
    c=(np.sin(th)**2)/(2*p.x_stddev.value**2) + (np.cos(th)**2)/(2*p.y_stddev.value**2)

    z=p.amplitude.value*np.exp(-(a*(x-p.x_mean.value)**2 - 2*b*(x-p.x_mean.value)*(y-p.y_mean.value) + c*(y-p.y_mean.value)**2  ))
    plt.contour(x,y,z, linewidths=2)

def plot_decision_kernel_boundary(X,y,scaler, sigma, clf,  xlabel, ylabel, legend):

    ax = plot_twoclass_data(X,y,xlabel,ylabel,legend)
    ax.autoscale(False)

    # create a mesh to plot in
    h = 0.05
    x1_min, x1_max = X[:, 0].min() - 1, X[:, 0].max() + 1
    x2_min, x2_max = X[:, 1].min() - 1, X[:, 1].max() + 1
    xx1, xx2 = np.meshgrid(np.arange(x1_min, x1_max, h),
                         np.arange(x2_min, x2_max, h))

    ZZ = np.array(np.c_[xx1.ravel(), xx2.ravel()])
    K = np.array([gaussian_kernel(x1,x2,sigma) for x1 in ZZ for x2 in X]).reshape((ZZ.shape[0],X.shape[0]))

    # need to scale it
    scaleK = scaler.transform(K)

    # and add the intercept column of ones
    KK = np.vstack([np.ones((scaleK.shape[0],)),scaleK.T]).T

    # make predictions on this mesh
    Z = clf.predict(KK)

    # Put the result into a contour plot
    Z = Z.reshape(xx1.shape)
    plt.contour(xx1,xx2,Z,cmap=plt.cm.gray,levels=[0.5])

def contourRegrid(l,S,limits,time,before,lvls=[],vlims=[]):
    plt.close()
    def makeplot():
#        plt.clabel(ph, inline=1, fontsize=10)        
#        plt.axis('equal')
        plt.axis(limits)
        if before:
            plt.title(str+' before regridding, time = %03f' % time)
            fname = os.path.expanduser('~/VEsims/') + str +'Regridding%02dBefore.pdf' % int(round(time))
        else:
            plt.title(str+' after regridding, time = %03f' % time)
            fname = os.path.expanduser('~/VEsims/') + str +'Regridding%02dRegrid.pdf' % int(round(time))
        plt.savefig(fname)
    ph=plt.contour(l[:,:,0],l[:,:,1],S[:,:,0,0],levels=lvls[0])
#    ph = plt.pcolor(l[:,:,0],l[:,:,1],S[:,:,0,0],vmin=vlims[0],vmax=vlims[1],cmap=cm.cool)
    plt.colorbar(ph)
    str='S11'
    makeplot()
    plt.clf()
    ph=plt.contour(l[:,:,0],l[:,:,1],S[:,:,1,1],levels=lvls[1])
#    ph = plt.pcolor(l[:,:,0],l[:,:,1],S[:,:,1,1],vmin=vlims[2],vmax=vlims[3],cmap=cm.cool)
    plt.colorbar(ph)
    str='S22'
    makeplot()
    plt.close()

def test_contour_heatmap():
    from scipy.interpolate import griddata
    from matplotlib import pyplot as plt

    mesh3D = mesh(200)
    mesh2D = proj_to_2D(mesh3D)

    data = np.zeros((3,3))
    data[0,1] += 2

    vals = np.exp(log_dirichlet_density(mesh3D,2.,data=data.sum(0)))
    temp = log_censored_dirichlet_density(mesh3D,2.,data=data)
    censored_vals = np.exp(temp - temp.max())

    xi = np.linspace(-1,1,1000)
    yi = np.linspace(-0.5,1,1000)

    plt.figure()
    plt.contour(griddata((mesh2D[:,0],mesh2D[:,1]),vals,(xi[None,:],yi[:,None]), method='cubic'))
    plt.axis('off')
    plt.title('uncensored likelihood')

    plt.figure()
    plt.contour(griddata((mesh2D[:,0],mesh2D[:,1]),censored_vals,(xi[None,:],yi[:,None]),method='cubic'))
    plt.axis('off')
    plt.title('censored likelihood')

def MonotoneRegion():
    matplotlib.rcParams['xtick.direction'] = 'out'
    matplotlib.rcParams['ytick.direction'] = 'out'

    delta = 0.025
    w2 = np.arange(-5.0, 5.0, delta)
    a = np.arange(0.0, 5.0, delta)
    W2, A = np.meshgrid(w2, a)
    LAM_cc = -1. + 5.*(A**2.) - np.sqrt(1. + A**4. + 2.*(A**2.)*(-1. + 8.*(W2**2.)))  # They're actually the same boundary!!!!
    # solve for when square root above is less than b^2 and you get same as below, LAM_eg
    LAM_eg = -1. + 3.*(A**2.) - 2.*(W2**2.)
    # You can force all the contours to be the same color.
    plt.figure()
    CS = plt.contour(W2, A, LAM_cc, 10,
                     colors='k',  # negative contours will be dashed by default
                     )
    plt.clabel(CS, fontsize=9, inline=1)
    plt.title('Monotone Region: Crossing the Curl')
    plt.savefig('monregion_cc.png')
    plt.close()

    plt.figure()
    CS = plt.contour(W2, A, LAM_eg, 10,
                     colors='k',  # negative contours will be dashed by default
                     )
    plt.clabel(CS, fontsize=9, inline=1)
    plt.title('Monotone Region: Extragradient')
    plt.savefig('monregion_eg.png')
    plt.close()

def multiple_optima(gene_number=937, resolution=80, model_restarts=10, seed=10000, max_iters=300, optimize=True, plot=True):
    """
    Show an example of a multimodal error surface for Gaussian process
    regression. Gene 939 has bimodal behaviour where the noisy mode is
    higher.
    """

    # Contour over a range of length scales and signal/noise ratios.
    length_scales = np.linspace(0.1, 60., resolution)
    log_SNRs = np.linspace(-3., 4., resolution)

    try:import pods
    except ImportError:
        print('pods unavailable, see https://github.com/sods/ods for example datasets')
        return
    data = pods.datasets.della_gatta_TRP63_gene_expression(data_set='della_gatta',gene_number=gene_number)
    # data['Y'] = data['Y'][0::2, :]
    # data['X'] = data['X'][0::2, :]

    data['Y'] = data['Y'] - np.mean(data['Y'])

    lls = GPy.examples.regression._contour_data(data, length_scales, log_SNRs, GPy.kern.RBF)
    if plot:
        pb.contour(length_scales, log_SNRs, np.exp(lls), 20, cmap=pb.cm.jet)
        ax = pb.gca()
        pb.xlabel('length scale')
        pb.ylabel('log_10 SNR')

        xlim = ax.get_xlim()
        ylim = ax.get_ylim()

    # Now run a few optimizations
    models = []
    optim_point_x = np.empty(2)
    optim_point_y = np.empty(2)
    np.random.seed(seed=seed)
    for i in range(0, model_restarts):
        # kern = GPy.kern.RBF(1, variance=np.random.exponential(1.), lengthscale=np.random.exponential(50.))
        kern = GPy.kern.RBF(1, variance=np.random.uniform(1e-3, 1), lengthscale=np.random.uniform(5, 50))

        m = GPy.models.GPRegression(data['X'], data['Y'], kernel=kern)
        m.likelihood.variance = np.random.uniform(1e-3, 1)
        optim_point_x[0] = m.rbf.lengthscale
        optim_point_y[0] = np.log10(m.rbf.variance) - np.log10(m.likelihood.variance);

        # optimize
        if optimize:
            m.optimize('scg', xtol=1e-6, ftol=1e-6, max_iters=max_iters)

        optim_point_x[1] = m.rbf.lengthscale
        optim_point_y[1] = np.log10(m.rbf.variance) - np.log10(m.likelihood.variance);

        if plot:
            pb.arrow(optim_point_x[0], optim_point_y[0], optim_point_x[1] - optim_point_x[0], optim_point_y[1] - optim_point_y[0], label=str(i), head_length=1, head_width=0.5, fc='k', ec='k')
        models.append(m)

    if plot:
        ax.set_xlim(xlim)
        ax.set_ylim(ylim)
    return m # (models, lls)

def three_d_contour(file,show):

    xloc = raw_input("xlocation(two_d_initial_x_value.txt):") or "two_d_initial_x_value.txt"
    yloc = raw_input("ylocation(two_d_initial_y_value.txt):") or "two_d_initial_y_value.txt"
    uloc = raw_input("ulocation(two_d_cavity_u_results.txt):") or "two_d_cavity_u_results.txt"
    vloc = raw_input("vlocation(two_d_cavity_v_results.txt):") or "two_d_cavity_v_results.txt"
    ploc = raw_input("plocation(two_d_cavity_p_results.txt):") or "two_d_cavity_p_results.txt"
    u=np.loadtxt(uloc)
    v=np.loadtxt(vloc)
    p=np.loadtxt(ploc)
    X=np.loadtxt(xloc)
    Y=np.loadtxt(yloc)
    X,Y=np.meshgrid(X,Y)
#    print X,Y
    fig = plt.figure()
    plt.contourf(X,Y,p,alpha=0.5)    ###plnttong the pressure field as a contour
    plt.colorbar()
    plt.contour(X,Y,p)               ###plotting the pressure field outlines
    plt.quiver(X[::3,::3],Y[::3,::3],u[::3,::3],v[::3,::3]) ##plotting velocity
    plt.xlabel('X')
    plt.ylabel('Y')
    plt.title(file)
    sep = '.'
    filename = file.split(sep, 1)[0]
    plt.savefig(filename+'.png')   # save figure as .png file
    if show:
        # print('test')
        plt.show()
    plt.close(fig)

def plot_mv_normal(mu, covarianceMatrix, xmin=-5, xmax=5, ymin=-5, ymax=5, resolution=50):
    X, Y = numpy.meshgrid(numpy.linspace(xmin, xmax, resolution),
                          numpy.linspace(ymin, ymax, resolution))

    # 2d array to store mv_normal pdf values
    pdfv = numpy.zeros(X.shape)

    # calculate pdf values at each X,Y location
    for i in range(X.shape[0]):
        for j in range(X.shape[1]):
            pdfv[i][j] = mv_normal(numpy.array([X[i][j], Y[i][j]]), mu, covarianceMatrix)

    # Contour plot
    plt.figure()
    plt.contour(X, Y, pdfv)
    plt.xlabel(r'$x$')
    plt.ylabel(r'$y$')

    # Surface plot
    fig = plt.figure()
    ax = fig.gca(projection='3d')
    surf = ax.plot_surface(X, Y, pdfv, rstride=1, cstride=1,
                           cmap=matplotlib.cm.jet,
                           linewidth=0, antialiased=False)
    plt.xlabel('$x$')
    plt.ylabel('$y$')
    fig.colorbar(surf)

    plt.show()

    def test(self):
        xlim=[-2.0,2.0]
        ylim=[-2.0,2.0]
        
        resol = 0.025
        sample_x = np.arange(-2.0,2.0,resol)
        sample_y = np.arange(-2.0,2.0,resol)
        sample_X, sample_Y = np.meshgrid(sample_x, sample_y)
        sample_XX = np.hstack([sample_X.reshape(sample_X.size,1),sample_Y.reshape(sample_Y.size,1)])

        sample_Z1 = np.exp(self.gmm_0.score(sample_XX))
        sample_Z2 = np.exp(self.gmm_1.score(sample_XX))

        plt.figure()
        ax1 = plt.subplot(121, aspect='equal')    
        CS = plt.contour(sample_X, sample_Y, sample_Z1.reshape(sample_X.shape))
        plt.clabel(CS, inline=1, fontsize=10)

        ax2 = plt.subplot(122, aspect='equal')    
        CS = plt.contour(sample_X, sample_Y, sample_Z2.reshape(sample_X.shape))
        plt.clabel(CS, inline=1, fontsize=10)

        ax1.set_xlim(xlim)
        ax1.set_ylim(ylim)

        ax2.set_xlim(xlim)
        ax2.set_ylim(ylim)

        plt.show()

        #print gmm_1.get_params(deep=True)
        print self.gmm_0.weights_
        print self.gmm_0.means_
        print self.gmm_0.covars_