def SVM(X, Y, ktype): 

	svmScores = []
	kf = KFold(n=len(Y), n_folds = nfolds)
	for train_index, test_index in kf:
		X_train, X_test = X[train_index], X[test_index]
		y_train, y_test = Y[train_index], Y[test_index]

		# SVC fit
		clf = svm.SVC(C=1.0, kernel= ktype)
		clf.fit(X_train, y_train)
		svmScores.append(clf.score(X_test, y_test))

		print "scores" , svmScores
		
		xx, yy = np.meshgrid(np.linspace(-10, 10, 500), np.linspace(-10, 10, 500))

		Z = clf.decision_function(np.c_[xx.ravel(), yy.ravel()])
		Z = Z.reshape(xx.shape)

		plt.contourf(xx, yy, Z, levels=np.linspace(Z.min(), 0, 7), cmap=pl.cm.Blues_r)
		a = pl.contour(xx, yy, Z, levels=[0], linewidths=2, colors='red')
		pl.contourf(xx, yy, Z, levels=[0, Z.max()], colors='orange')

		colors = ['w' if i == 0 else 'k' for i in Y]
		plt.scatter(X[:,0], X[:,1], color = colors, alpha=1.0)

		# Plt that plot yoz. and make the xylim not shitty
		plt.xlim([np.min(X)-5,np.max(X)+5])
		plt.ylim([np.min(X)-5,np.max(X)+5])
		plt.show()

def Contour(X,Y,Z, label='', levels=None, cmapidx=0, colors=None, fmt='%g', lwd=1, fsz=10,
        inline=0, wire=True, cbar=True, zorder=None, markZero='', clabels=True):
    """
    Plot contour
    ============
    """
    L = None
    if levels != None:
        if not hasattr(levels, "__iter__"): # not a list or array...
            levels = linspace(Z.min(), Z.max(), levels)
    if colors==None:
        c1 = contourf (X,Y,Z, cmap=Cmap(cmapidx), levels=levels, zorder=None)
    else:
        c1 = contourf (X,Y,Z, colors=colors, levels=levels, zorder=None)
    if wire:
        c2 = contour (X,Y,Z, colors=('k'), levels=levels, linewidths=[lwd], zorder=None)
        if clabels:
            clabel (c2, inline=inline, fontsize=fsz)
    if cbar:
        cb = colorbar (c1, format=fmt)
        cb.ax.set_ylabel (label)
    if markZero:
        c3 = contour(X,Y,Z, levels=[0], colors=[markZero], linewidths=[2])
        if clabels:
            clabel(c3, inline=inline, fontsize=fsz)

def ScatterWithLinearSVC(X_,class_,C_,TrainedSVC,Labels_):
    ##"Train" the data, aka fit with vectorial distance to a line, since we are using a linear kernel
    ## Here we essentially minimize the total distance that each point is to a decision boundary
    ## Based on which side of the boundary the point lies (if distance is "positive" or "negative"), the decision is made
    #TrainedSVC = SVC(C = C_, kernel = 'linear').fit(X_,class_)
    if np.shape(X_)[1] >2:
        print "X is larger than 2!"
    Xvalues=X_[:,0]
    Yvalues=X_[:,1]
    xmax, ymax, xmin,ymin =Xvalues.max(), Yvalues.max(), Xvalues.min(), Yvalues.min()
    binning=100.
    ### Make a big grid of coordinates, as a matrix of X positions and a matrix of y positions
    X, Y = np.meshgrid(np.arange(xmin,xmax, (xmax-xmin)/binning),
                       np.arange(ymin,ymax, (ymax-ymin)/binning))
    ### Ravel the X and Y matrices up into vectors, put them together, and feed them to the predict function
    GridPredictions=TrainedSVC.predict(np.c_[X.ravel(), Y.ravel()])
    ### Re-form a matrix of Grid predictions
    GridPredictions=np.reshape(GridPredictions,np.shape(X))
    ### Plot a contour, with coloring each area, fading to 1/3 of the color "strength"
    fig=plt.figure()
    plt.contourf(X,Y,GridPredictions,alpha=0.33)
    plt.scatter(Xvalues, Yvalues,c=class_)   
    plt.scatter(Xvalues, Yvalues,c=class_)   
    plt.ylabel(Labels_[0])
    plt.xlabel(Labels_[1])
    #plt.legend([GridPredictions], ["Training accuracy"])

    fig.savefig('result.png')

def render(
    src: array,
    height: int,
    width: int,
    isostep: int,
    dpi: int,
    output_file_path: Path,
) -> None:

    height = height // MAP_SCALE
    width = width // MAP_SCALE

    image_size = (
        (width / dpi),
        (height / dpi),
    )

    data = np.array(src).reshape((height, width))
    isohypses = list(range(0, data.max(), isostep))

    plt.clf()
    plt.axis('off')
    plt.figure(dpi=dpi)

    fig = plt.figure(figsize=image_size, frameon=False)
    fig.add_axes([0, 0, 1, 1])

    contourf(data, 256, cmap=CMAP)
    contour(data, isohypses, colors=ISOHYPSE_COLOR, linewidths=ISOHYPSE_WIDTH)

    output_file_path.parent.parent.mkdir(parents=True, exist_ok=True)
    plt.savefig(str(output_file_path), bbox_inches=0, dpi=dpi)

        def doSubplot(multiplier=1.0, layout=(-1, -1)):
            exp_vort = cPickle.load(open("vort_time_height_%s.pkl" % exp[5:], 'r'))
            vort_data[exp] = exp_vort
            print "Max vorticity:", exp_vort.max()

            cmap = cm.get_cmap('RdYlBu_r')
            cmap.set_under('#ffffff')

#           pylab.pcolormesh(boundCoordinate(np.array(temp.getTimes())), boundCoordinate(z_column / 1000.), exp_vort.T, cmap=cmap, vmin=min_vorticity, vmax=max_vorticity)
            pylab.contourf(temp.getTimes(), z_column / 1000., exp_vort.T, cmap=cmap, levels=np.arange(min_vorticity, max_vorticity + 0.0024, 0.0025))

            pylab.xlim(temp.getTimes()[0], temp.getTimes()[-1])
            pylab.ylim([z_column[0] / 1000, 10])

            layout_r, layout_c = layout
            if layout_r == 2 and layout_c == 1 or layout_r == -1 and layout_c == -1:
                pylab.xlabel("Time (UTC)", size='large')
                pylab.ylabel("Height (km MSL)", size='large')

#           pylab.axvline(14400, color='k', linestyle=':')

                pylab.xticks(temp.getTimes(), temp.getStrings('%H%M', aslist=True), size='large')
                pylab.yticks(size='large')
            else:
                pylab.xticks(temp.getTimes(), [ '' for t in temp ])
                pylab.yticks(pylab.yticks()[0], [ '' for z in pylab.yticks()[0] ])

            pylab.text(0.05, 0.95, experiments[exp], ha='left', va='top', transform=pylab.gca().transAxes, size=18 * multiplier)
            pylab.xlim(temp.getTimes()[0], temp.getTimes()[3])
            pylab.ylim([z_column[1] / 1000, 6])
#           pylab.title(exp)
            return

        def doSubplot(multiplier=1.0, layout=(-1, -1)):
            data = cPickle.load(open("cold_pool_%s.pkl" % exp, 'r'))
            wdt = temporal.getTimes().index(time_sec)

            try:
                mo = ARPSModelObsFile("%s/%s/KCYSan%06d" % (base_path, exp, time_sec))
            except AssertionError:
                mo = ARPSModelObsFile("%s/%s/KCYSan%06d" % (base_path, exp, time_sec), mpi_config=(2, 12))
            except:
                print "Can't load reflectivity ..."
                mo = {'Z':np.zeros((1, 255, 255), dtype=np.float32)}

            cmap = matplotlib.cm.get_cmap('Blues_r')
            cmap.set_over('#ffffff')
#           cmap.set_under(tuple( cmap._segmentdata[c][0][-1] for c in ['red', 'green', 'blue'] ))
#           cmap.set_under(cmap[0])

            xs, ys = grid.getXY()
#           pylab.pcolormesh(xs, ys, data['t'][2][domain_bounds], cmap=cmap, vmin=288., vmax=295.)
            pylab.contour(xs, ys, mo['Z'][0][domain_bounds], levels=np.arange(10, 80, 10), colors='k', zorder=10)
            pylab.contourf(xs, ys, data['t'][wdt][domain_bounds], levels=range(289, 296), cmap=cmap)

            grid.drawPolitical(scale_len=10)

#           pylab.plot(track_xs, track_ys, 'mv-', lw=2.5, mfc='k', ms=8)

            pylab.text(0.05, 0.95, experiments[exp], ha='left', va='top', transform=pylab.gca().transAxes, size=14 * multiplier)
            return

def plotComparison(ens_mean, ens_ob_mean, ens_ob_std, obs, ob_locations, refl, grid, levels, cmap, title, file_name):
    max_std = 5.0

    xs, ys = grid.getXY()
    obs_xs, obs_ys = grid(*ob_locations)

    clip_box = Bbox([[0, 0], [1, 1]])

    pylab.figure()

    pylab.contourf(xs, ys, ens_mean, cmap=cmap, levels=levels)
    pylab.colorbar()

    pylab.contour(xs, ys, refl, colors='k', levels=np.arange(20, 80, 20))

    for ob_x, ob_y, ob, ob_mean, ob_std in zip(obs_xs, obs_ys, obs, ens_ob_mean, ens_ob_std):
        color_bin = np.argmin(np.abs(ob - levels))
        if ob > levels[color_bin]: color_bin += 1
        color_level = float(color_bin) / len(levels)

        ob_z_score = (ob - ob_mean) / ob_std

        print "Ob z-score:", ob_z_score, "Ob:", ob, "Ob mean:", ob_mean, "Ob std:", ob_std

        pylab.plot(ob_x, ob_y, 'ko', markerfacecolor=cmap(color_level), markeredgecolor=std_cmap(ob_z_score / max_std), markersize=4, markeredgewidth=1)
#       pylab.text(ob_x - 1000, ob_y + 1000, "%5.1f" % temp_K, ha='right', va='bottom', size='xx-small', clip_box=clip_box, clip_on=True)

    grid.drawPolitical()

    pylab.suptitle(title)
    pylab.savefig(file_name)
    pylab.close()
    return

def main():
    zbar = 0.7
    vbar = 500  # km/s comoving distance
    ncc = 1024
    npp = 200000  # total mass = 1e9*2e5 M_sun/h
    bsz = 3.0  # Mpc/h comoving distance

    file_particles = sys.argv[1]

    sdens = call_sph_sdens(file_particles, bsz, ncc, npp)
    sdens.astype(np.float32).tofile("./output_files/cnfw_sdens.bin")

    ksz_map = ksz_based_on_sdens(zbar, vbar, sdens)
    # ---------------------------
    # Save 2d array to a binary file.
    ksz_map.astype(np.float32).tofile("./output_files/cnfw_ksz.bin")
    # ---------------------------

    # ---------------------------
    # Read the output binary files.
    a0 = np.fromfile("./output_files/cnfw_ksz.bin", dtype=np.float32)
    a0 = a0.reshape((ncc, ncc))
    # ---------------------------

    # ---------------------------
    # Plot the contours
    import pylab as pl

    pl.contourf(np.log10(a0))
    pl.colorbar()
    pl.show()

    return 0

def plot_function_contour(x, y, function, **kwargs):
    """Make a contour plot of a function of two variables.

    Parameters
    ----------
    x, y : array_like of float
        The positions of the nodes of a plotting grid.
    function : function
        The function to plot.
    filling : bool
        Fill contours if True (default).
    num_contours : int
        The number of contours to plot, 50 by default.
    xlabel, ylabel : str, optional
        The axes labels. Empty by default.
    title : str, optional
        The title. Empty by default.

    """
    X, Y = numpy.meshgrid(x, y)
    Z = []
    for y_value in y:
        Z.append([])
        for x_value in x:
            Z[-1].append(function(x_value, y_value))
    Z = numpy.array(Z)
    num_contours = kwargs.get('num_contours', 50)
    if kwargs.get('filling', True):
        pylab.contourf(X, Y, Z, num_contours, cmap=pylab.cm.jet)
    else:
        pylab.contour(X, Y, Z, num_contours, cmap=pylab.cm.jet)
    pylab.xlabel(kwargs.get('xlabel', ''))
    pylab.ylabel(kwargs.get('ylabel', ''))
    pylab.title(kwargs.get('title', ''))

        def doSubplot(multiplier=1.0, layout=(-1, -1)):
            pylab.contour(xs, ys, vr_obs[0][bounds], colors='k', linestyles='-', levels=np.arange(-50, 60, 10))
            pylab.contour(xs, ys, model_obs['vr'][0][bounds], colors='k', linestyles='--', levels=np.arange(-50, 60, 10))
            pylab.contourf(xs, ys, (model_obs['vr'][0] - vr_obs[0])[bounds], cmap=matplotlib.cm.get_cmap('RdBu_r'), zorder=-1)

            grid.drawPolitical()
            return

def main():
    buffer_width=40000
    reports = loadReports("snowfall_2013022600_m.txt")
#   map = Basemap(projection='lcc', resolution='i', area_thresh=10000.,
#                 llcrnrlat=20., llcrnrlon=-120., urcrnrlat=48., urcrnrlon=-60.,
#                 lat_0=34.95, lon_0=-97.0, lat_1=30., lat_2=60.)

    map = Basemap(projection='lcc', resolution='i', area_thresh=10000.,
                  llcrnrlat=25., llcrnrlon=-110., urcrnrlat=40., urcrnrlon=-90.,
                  lat_0=34.95, lon_0=-97.0, lat_1=30., lat_2=60.)

    width, height = map(map.urcrnrlon, map.urcrnrlat)
    report_xs, report_ys = map(reports['longitude'], reports['latitude'])

    keep_idxs = np.where((report_xs >= -buffer_width) & (report_xs <= width + buffer_width) & (report_ys >= -buffer_width) & (report_ys <= height + buffer_width))

    spacing = findAvgStationSpacing(*map(reports['longitude'][keep_idxs], reports['latitude'][keep_idxs]))

    dx = dy = floor(spacing * 0.75 / 2000) * 2000
    print dx, dy
    xs, ys = np.meshgrid(np.arange(0, width, dx), np.arange(0, height, dy))

    grid = griddata((report_xs[keep_idxs], report_ys[keep_idxs]), reports['amount'][keep_idxs], (xs, ys))

    pylab.contourf(xs, ys, grid / 2.54, levels=np.arange(2, 22, 2), cmap=pylab.cool())
    pylab.colorbar(orientation='horizontal')

    map.drawcoastlines()
    map.drawcountries()
    map.drawstates()
    pylab.savefig("snow.png")
    return

 def view_pcem(self,PCEM):
   nbins = np.sqrt(len(PCEM))
   levels=[0.0,0.01,0.02,0.03,0.04,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5]
   #levels2=[0.0,0.01,0.05,0.1,0.2,0.3,0.4,0.5,1.0]
 
   T1 = PCEM[:,0].reshape( (nbins,nbins))
   T2 = PCEM[:,1].reshape( (nbins,nbins))
   E  = PCEM[:,2].reshape( (nbins,nbins))
   EC  = PCEM[:,3].reshape( (nbins,nbins))
   M  = PCEM[:,4].reshape( (nbins,nbins))
   pp.figure()
   pp.subplot(1,2,1)
   pp.contour( T1, T2, E, levels,linewidths=0.5,colors="k",linestyles='solid'  )
   pp.contourf( T1, T2, E, levels, alpha=0.75  )
   pp.xlabel( "proposal")
   pp.ylabel( "current")
   pp.title( "Error" )
   pp.colorbar()
   # pp.subplot(1,3,2)
   # pp.contour( T1, T2, EC, levels,linewidths=0.5,colors="k",linestyles='solid'  )
   # pp.contourf( T1, T2, EC, levels, alpha=0.75  )
   # pp.xlabel( "proposal")
   # pp.ylabel( "current")
   # pp.colorbar()
   # pp.title( "Corrected Error" )
   pp.subplot(1,2,2)
   pp.contour( T1, T2, M, 20,linewidths=1.0,colors="k",linestyles='solid'  )
   pp.contourf( T1, T2, M, 20  )
   pp.xlabel( "proposal")
   pp.ylabel( "current")
   pp.colorbar()
   pp.title("MU Z")

def trainStep(fnn, trainer, trndata, tstdata):
    trainer.trainEpochs(1)
    trnresult = percentError(trainer.testOnClassData(), trndata["class"])
    tstresult = percentError(trainer.testOnClassData(dataset=tstdata), tstdata["class"])

    print "epoch: %4d" % trainer.totalepochs, "  train error: %5.2f%%" % trnresult, "  test error: %5.2f%%" % tstresult

    out = fnn.activateOnDataset(griddata)
    out = out.argmax(axis=1)  # the highest output activation gives the class
    out = out.reshape(X.shape)

    figure(1)
    ioff()  # interactive graphics off
    clf()  # clear the plot
    hold(True)  # overplot on
    for c in [0, 1, 2]:
        here, _ = where(trndata["class"] == c)
        plot(trndata["input"][here, 0], trndata["input"][here, 1], "o")
    if out.max() != out.min():  # safety check against flat field
        contourf(X, Y, out)  # plot the contour
    ion()  # interactive graphics on
    draw()  # update the plot

    figure(2)
    ioff()  # interactive graphics off
    clf()  # clear the plot
    hold(True)  # overplot on
    for c in [0, 1, 2]:
        here, _ = where(tstdata["class"] == c)
        plot(tstdata["input"][here, 0], tstdata["input"][here, 1], "o")
    if out.max() != out.min():  # safety check against flat field
        contourf(X, Y, out)  # plot the contour
    ion()  # interactive graphics on
    draw()  # update the plot

def main():

    nnn = 512
    boxsize = 4.0
    dsx = boxsize/nnn
    xi1 = np.linspace(-boxsize/2.0,boxsize/2.0-dsx,nnn)+0.5*dsx
    xi2 = np.linspace(-boxsize/2.0,boxsize/2.0-dsx,nnn)+0.5*dsx
    xi1,xi2 = np.meshgrid(xi1,xi2)

    #----------------------------------------------------
    # lens parameters for main halo
    xlc1 = 0.0
    xlc2 = 0.0
    ql0 = 0.799999999999
    rc0 = 0.100000000001
    re0 = 1.0
    phi0 = 0.0

    g_ycen = 0.0
    g_xcen = 0.0

    phi,td,ai1,ai2,kappa,mu,yi1,yi2 = nie_all(xi1,xi2,xlc1,xlc2,re0,rc0,ql0,phi0,g_ycen,g_xcen)

    pl.figure(figsize=(10,10))
    pl.contourf(td)
    pl.show()

 def plot(self,scatter=None,screen=True,out=None):
     """ Plot the VLEP landscape. """
     X = np.linspace(0.5,3.2,100)
     Y = np.linspace(-2,1.4,100)
     Z=np.zeros((100,100))
     for i,x in enumerate(X):
         for j,y in enumerate(Y):
             Z[j,i]=self._function(x,y)
             
     pl.contourf(X,Y,Z,100)
     pl.hot()
     
     if scatter!=None:
         pl.scatter(scatter[0,:],scatter[1,:],color='b')
     #if plot!=None:
         #pl.plot(plot[0,:],plot[1,:],color='y')
         #pl.scatter(plot[0,:],plot[1,:],color='y')
         
     if screen==True:
         pl.show()
     else:
         assert out!=None
         pl.savefig(out)
         pl.clf()
     self.it+=1

        def doSubplot(multiplier=1.0, layout=(-1, -1)):
            if time_sec < 16200:
                xs, ys = xs_1, ys_1
                domain_bounds = bounds_1sthalf
                grid = grid_1
            else:
                xs, ys = xs_2, ys_2
                domain_bounds = bounds_2ndhalf
                grid = grid_2

            try:
                mo = ARPSModelObsFile("%s/%s/KCYS%03dan%06d" % (base_path, exp, min_ens, time_sec))
            except AssertionError:
                mo = ARPSModelObsFile("%s/%s/KCYS%03dan%06d" % (base_path, exp, min_ens, time_sec), mpi_config=(2, 12))
            except:
                print "Can't load reflectivity ..."
                mo = {'Z':np.zeros((1, 255, 255), dtype=np.float32)}

            pylab.contour(xs, ys, wind[exp]['w'][wdt][domain_bounds], levels=np.arange(2, 102, 2), styles='-', colors='k')
            pylab.contour(xs, ys, wind[exp]['w'][wdt][domain_bounds], levels=np.arange(-100, 0, 2), styles='--', colors='k')

            pylab.quiver(xs[thin], ys[thin], wind[exp]['u'][wdt][domain_bounds][thin], wind[exp]['v'][wdt][domain_bounds][thin])

            pylab.contourf(xs, ys, mo['Z'][0][domain_bounds], levels=np.arange(10, 85, 5), cmap=NWSRef, zorder=-10)

            grid.drawPolitical(scale_len=10)

            row, col = layout
            if col == 1:
                pylab.text(-0.075, 0.5, exp_names[exp], transform=pylab.gca().transAxes, rotation=90, ha='center', va='center', size=12 * multiplier)

def main():
    base_path = "/caps2/tsupinie/1kmf-control/"
    temp = goshen_1km_temporal(start=14400, end=14400)
    grid = goshen_1km_grid()
    n_ens_members = 40

    np.seterr(all='ignore')

    ens = loadEnsemble(base_path, [ 11 ], temp.getTimes(), ([ 'pt', 'p' ], computeDensity))
    ens = ens[0, 0]

    zs = decompressVariable(nio.open_file("%s/ena001.hdfgrdbas" % base_path, mode='r', format='hdf').variables['zp'])
    xs, ys = grid.getXY()
    xs = xs[np.newaxis, ...].repeat(zs.shape[0], axis=0)
    ys = ys[np.newaxis, ...].repeat(zs.shape[0], axis=0)

    eff_buoy = effectiveBuoyancy(ens, (zs, ys, xs), plane={'z':10})
    print eff_buoy

    pylab.figure()
    pylab.contourf(xs[0], ys[0], eff_buoy[0], cmap=matplotlib.cm.get_cmap('RdBu_r'))
    pylab.colorbar()

    grid.drawPolitical()

    pylab.suptitle("Effective Buoyancy")
    pylab.savefig("eff_buoy.png")
    pylab.close()
    return

def plotSurface(pt, td, winds, map, stride, title, file_name):
    pylab.figure()
    pylab.axes((0.05, 0.025, 0.9, 0.9))

    u, v = winds

    nx, ny = pt.shape
    gs_x, gs_y = goshen_3km_gs
    xs, ys = np.meshgrid(gs_x * np.arange(nx), gs_y * np.arange(ny))   

    data_thin = tuple([ slice(None, None, stride) ] * 2)

    td_cmap = matplotlib.cm.get_cmap('Greens')
    td_cmap.set_under('#ffffff')
    pylab.contourf(xs, ys, td, levels=np.arange(40, 80, 5), cmap=td_cmap)
    pylab.colorbar()
    CS = pylab.contour(xs, ys, pt, colors='r', linestyles='-', linewidths=1.5, levels=np.arange(288, 324, 4))
    pylab.clabel(CS, inline_spacing=0, fmt="%d K", fontsize='x-small')
    pylab.quiver(xs[data_thin], ys[data_thin], u[data_thin], v[data_thin])

    drawPolitical(map, scale_len=75)

    pylab.suptitle(title)
    pylab.savefig(file_name)
    pylab.close()
    return

  def __init__(self, name, nx, ng, np, opts={}):
    Callback.__init__(self)

    self.nx = nx
    self.ng = ng
    self.np = np

    opts['input_scheme'] = nlpsol_out()
    opts['output_scheme'] = ['ret']

    figure(1)
    subplot(111)
    
    x_,y_ = mgrid[-1:1.5:0.01,-1:1.5:0.01]
    z_ = DM.zeros(x_.shape)
    
    for i in range(x_.shape[0]):
      for j in range(x_.shape[1]):
        z_[i,j] = fcn(x_[i,j],y_[i,j])
    contourf(x_,y_,z_)
    colorbar()
    title('Iterations of Rosenbrock')
    draw()
    
    self.x_sols = []
    self.y_sols = []

    # Initialize internal objects
    self.construct(name, opts)

def postplot(x,y,arr,plottitle,showplot=True):
    """ Post process in 2D"""
    levels = linspace(min(arr.ravel())*0.9,max(arr.ravel())*1.1,50)
    contourf(x,y,arr.T,levels=levels);
    title(plottitle); colorbar();
    show()
    return