def i5():
    # griddata
    r, dr = 10, 4
    dth = dr/(1.0*r)
    
    rr = linspace(r,r+dr,15)
    th = linspace(-0.5*dth,0.5*dth,16)
    R,TH = make_grid(rr,th)
    X,Y = R*cos(TH),R*sin(TH)
    Z = ff(X,Y)
    points = grid_to_points(X,Y)
    values = grid_to_points(Z)
    
    xfine = linspace(r,r+dr,50)
    yfine = linspace(-0.5*dr,0.5*dr,51)
    XF, YF = make_grid(xfine, yfine)
    desired_points = grid_to_points(XF, YF)

    # Only 'nearest' seems to work
    # 'linear' and 'cubic' just give fill value
    desired_values = scipy.interpolate.griddata(points, values, desired_points, method='nearest', fill_value=1.0)
    ZF = desired_values.reshape(XF.shape)
    
    pl.clf()
    pl.pcolormesh(XF, YF, ZF)
    pl.colorbar()
    return ZF

def i6():
    # Rbf: this works great except that somewhat weird stuff happens
    # when you're off the grid.
    r, dr = 10, 4
    dth = dr/(1.0*r)
    
    rr = linspace(r,r+dr,25)
    th = linspace(-0.5*dth,0.5*dth,26)
    R,TH = make_grid(rr,th)
    X,Y = R*cos(TH),R*sin(TH)
    Z = ff(X,Y)
    xlist, ylist, zlist = X.ravel(), Y.ravel(), Z.ravel()
    
    # function, epsilon, smooth, norm, 
    fint = scipy.interpolate.Rbf(xlist, ylist, zlist )

    xfine = linspace(-10+r,r+dr+10,99)
    yfine = linspace(-10-0.5*dr,10+0.5*dr,101)
    XF, YF = make_grid(xfine, yfine)
    ZF = fint(XF, YF)
    
    pl.clf()
    pl.pcolormesh(XF, YF, ZF)
    pl.colorbar()
    return ZF

def i1():
    # interp2d
    # This does not seem to work.  It also does not seem to honor
    # bounds_error or fill_value

    #xx = linspace(0,2*pi,39)
    #yy = linspace(-2,2,41)
    xx = linspace(0,2*pi,12)
    yy = linspace(-2,2,11)
    X,Y = make_grid(xx,yy)
    Z = ff(X,Y)

    # Note that in this use, xx,yy and 1d, Z is 2d.  
    # The output shapes are 
    # fint: ()   ()   => (1,)
    # fint: (n,) ()   => (n,)
    # fint: ()   (m,) => (1,m)
    # fint: (n,) (m,) => (n,m)

    fint = scipy.interpolate.interp2d(xx,yy, Z)

    xfine = linspace(0.0,2*pi,99)
    yfine = linspace(-2,2,101)
    XF, YF = make_grid(xfine, yfine)
    ZF = fint(xfine,yfine).transpose()

    pl.clf()
    pl.pcolormesh(XF,YF,ZF)
    pl.colorbar()

def i2():
    # interp2d -- do it on a structured gird, but use calling
    # convention for unstructured grid.

    xx = linspace(0,2*pi,15)
    yy = linspace(-2,2,14)
    X,Y = make_grid(xx,yy)
    Z = ff(X,Y)

    # The output shapes are 
    # fint: ()   ()   => (1,)
    # fint: (n,) ()   => (n,)
    # fint: ()   (m,) => (1,m)
    # fint: (n,) (m,) => (n,m)
    # 
    # Linear interpolation is all messed up. Cant get sensible results.
    # Cubic seems extremely sensitive to the exact number of data points.
    # Doesn't respect bounds_error
    # Doesn't respect fill_value
    # Doesn't respect copy
    fint = scipy.interpolate.interp2d(X,Y,Z, kind='quintic')
    
    xfine = linspace(-2,2*pi+2,99)
    yfine = linspace(-4,4,101)
    XF, YF = make_grid(xfine, yfine)
    # NB -- TRANSPOSE HERE!
    ZF = fint(xfine,yfine).transpose()

    pl.clf()
    pl.pcolormesh(XF, YF, ZF)
    pl.colorbar()

def logRegress(X,Y):

	scores = []
	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]
		logModel = linear_model.LogisticRegression()
		logModel.fit(X_train,y_train)
		scores.append(logModel.score(X_test, y_test))
		
		print "Scores" , scores

		x_min, x_max = X[:, 0].min() - .5, X[:, 0].max() + .5
		y_min, y_max = X[:, 1].min() - .5, X[:, 1].max() + .5
		xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))
		Z = logreg.predict(np.c_[xx.ravel(), yy.ravel()])

		# Put the result into a color plot
		Z = Z.reshape(xx.shape)
		pl.figure(1, figsize=(4, 3))
		pl.pcolormesh(xx, yy, Z, cmap=pl.cm.Paired)

		# Plot also the training points
		pl.scatter(X[:, 0], X[:, 1], c=Y, edgecolors='k', cmap=pl.cm.Paired)
		pl.xlabel('Sepal length')
		pl.ylabel('Sepal width')

		pl.xlim(xx.min(), xx.max())
		pl.ylim(yy.min(), yy.max())
		pl.xticks(())
		pl.yticks(())

		pl.show()

def the_plot():
    x,y = linspace(0,2*pi,100), linspace(-2,2,100)
    X,Y = make_grid(x,y)
    Z = ff(X,Y)
    pl.clf()
    pl.pcolormesh(X,Y,Z)
    pl.colorbar()

def main(plot=True):

    # Setup grid
    g = 4
    nx, ny = 200, 50
    Lx, Ly = 26, 26/4

    (x, y), (dx, dy) = ghosted_grid([nx, ny], [Lx, Ly], 0)

    # monkey patch the velocity
    uc = MultiFab(sizes=[nx, ny], n_ghost=4, dof=2)
    uc.validview[0] = (y > Ly / 3) * (2 * Ly / 3 > y)
    uc.validview[1] = np.sin(2 * pi * x / Lx) * .3 / (2 * pi / Lx)
    # state.u = np.random.rand(*x.shape)

    tad = BarotropicSolver()
    tad.geom.dx = dx
    tad.geom.dy = dx

    dt = min(dx, dy) / 4

    if plot:
        import pylab as pl
        pl.ion()

    for i, (t, uc) in enumerate(steps(tad.onestep, uc, dt, [0.0, 10000 * dt])):
        if i % 100 == 0:
            if plot:
                pl.clf()
                pl.pcolormesh(uc.validview[0])
                pl.colorbar()
                pl.pause(.01)

def plot_results_with_hyperplane(clf, clf_name, df, plt_nmbr):
    x_min, x_max = df.x.min() - .5, df.x.max() + .5
    y_min, y_max = df.y.min() - .5, df.y.max() + .5

    # step between points. i.e. [0, 0.02, 0.04, ...]
    step = .02
    # to plot the boundary, we're going to create a matrix of every possible point
    # then label each point as a wolf or cow using our classifier
    xx, yy = np.meshgrid(np.arange(x_min, x_max, step), np.arange(y_min, y_max, step))
    Z = clf.predict(np.c_[xx.ravel(), yy.ravel()])
    # this gets our predictions back into a matrix
    Z = Z.reshape(xx.shape)

    # create a subplot (we're going to have more than 1 plot on a given image)
    pl.subplot(2, 2, plt_nmbr)
    # plot the boundaries
    pl.pcolormesh(xx, yy, Z, cmap=pl.cm.Paired)

    # plot the wolves and cows
    for animal in df.animal.unique():
        pl.scatter(df[df.animal==animal].x,
                   df[df.animal==animal].y,
                   marker=animal,
                   label="cows" if animal=="x" else "wolves",
                   color='black',
                   c=df.animal_type, cmap=pl.cm.Paired)
    pl.title(clf_name)
    pl.legend(loc="best")

def _make_logp_histogram(points, logp, nbins, rangeci, weights, cbar):
    from numpy import (ones_like, searchsorted, linspace, cumsum, diff, 
        sort, argsort, array, hstack, log10, exp)
    if weights == None: weights = ones_like(logp)
    edges = linspace(rangeci[0],rangeci[1],nbins+1)
    idx = searchsorted(points, edges)
    weightsum = cumsum(weights)
    heights = diff(weightsum[idx])/weightsum[-1]  # normalized weights

    import pylab
    vmin,vmax,cmap = cbar
    cmap_steps = linspace(vmin,vmax,cmap.N+1)
    bins = [] # marginalized maximum likelihood
    for h,s,e,xlo,xhi in zip(heights,idx[:-1],idx[1:],edges[:-1],edges[1:]):
        if s == e: continue
        pv = -logp[s:e]
        pidx = argsort(pv)
        pw = weights[s:e][pidx]
        x = array([xlo,xhi],'d')
        y = hstack((0,cumsum(pw)))  
        z = pv[pidx][:,None]
        # centerpoint, histogram height, maximum likelihood for each bin
        bins.append(((xlo+xhi)/2,y[-1],exp(vmin-z[0])))
        if len(z) > cmap.N:
           # downsample histogram bar according to number of colors
           pidx = searchsorted(z[1:-1].flatten(), cmap_steps)
           if pidx[-1] < len(z)-1: pidx = hstack((pidx,-1))
           y,z = y[pidx],z[pidx]
        pylab.pcolormesh(x,y,z,vmin=vmin,vmax=vmax,hold=True,cmap=cmap)
        # Draw bars around each histogram bin
        #pylab.plot([xlo,xlo,xhi,xhi],[y[0],y[-1],y[-1],y[0]],'-k',linewidth=0.1,hold=True)
    centers,height,maxlikelihood = array(bins).T
    pylab.plot(centers, maxlikelihood*max(height), '-g', hold=True)

    def load_data(self):
        self.input_file = self.get_input()

        if self.input_file is None:
            print "No input file selected. Exiting..."
            import sys
            sys.exit(0)

        self.nc = NC(self.input_file)
        nc = self.nc

        self.x = np.array(nc.variables['x'][:], dtype=np.double)
        self.y = np.array(nc.variables['y'][:], dtype=np.double)
        self.z = np.array(np.squeeze(nc.variables['usurf'][:]), dtype=np.double)
        self.thk = np.array(np.squeeze(nc.variables['thk'][:]), dtype=np.double)

        self.mask = dbg.initialize_mask(self.thk)
        print "Mask initialization: done"

        plt.figure(1)
        plt.pcolormesh(self.x, self.y, self.mask)
        plt.contour(self.x, self.y, self.z, colors='black')
        plt.axis('tight')
        plt.axes().set_aspect('equal')
        plt.show()

    def compute_mask(self):
        import matplotlib.nxutils as nx

        if self.pts is not None:
            def correct_mask(mask, x, y, pts):
                    for j in range(y.size):
                        for i in range(x.size):
                            if mask[j,i] > 0:
                                if nx.pnpoly(x[i], y[j], pts):
                                    mask[j,i] = 2
                                else:
                                    mask[j,i] = 1

            correct_mask(self.mask, self.x, self.y, self.pts)

        dbg.upslope_area(self.x, self.y, self.z, self.mask)
        print "Drainage basin computation: done"
        self.mask_computed = True

        plt.figure(1)
        plt.pcolormesh(self.x, self.y, self.mask)
        plt.contour(self.x, self.y, self.z, colors='black')
        plt.axis('tight')
        plt.axes().set_aspect('equal')
        plt.show()

def plotRadTilt(plot_data, plot_cmaps, plot_titles, grid, file_name, base_ref=None):
    subplot_base = 220
    n_plots = 4

    xs, ys, gs_x, gs_y, map = grid

    pylab.figure(figsize=(12,12))
    pylab.subplots_adjust(left=0.02, right=0.98, top=0.98, bottom=0.02, wspace=0.04)

    for plot in range(n_plots):
        pylab.subplot(subplot_base + plot + 1)
        cmap, vmin, vmax = plot_cmaps[plot]
        cmap.set_under("#ffffff", alpha=0.0)

        pylab.pcolormesh(xs - gs_x / 2, ys - gs_y / 2, plot_data[plot] >= -90, vmin=0, vmax=1, cmap=mask_cmap)

        pylab.pcolormesh(xs - gs_x / 2, ys - gs_y / 2, plot_data[plot], vmin=vmin, vmax=vmax, cmap=cmap)
        pylab.colorbar()

        if base_ref is not None:
            pylab.contour(xs, ys, base_ref, levels=[10.], colors='k', lw=0.5)

#       if plot == 0:
#           pylab.contour(xs, ys, refl_88D[:, :, 0], levels=np.arange(10, 80, 10), colors='#808080', lw=0.5)

#       pylab.plot(radar_x, radar_y, 'ko')

        pylab.title(plot_titles[plot])

        drawPolitical(map)

    pylab.savefig(file_name)
    pylab.close()

    return

def plotGridPref(gridscore, clfName, obj , metric = 'roc_auc'):
    ''' Plot Grid Performance
    '''
    # data_path = data_path+obj+'/'+clfName+'_opt.h5'
    # f=hp.File(data_path, 'r')
    # gridscore = f['grids_score'].value
    # Get numblocks
    CV = np.unique(gridscore["i_CV"])
    folds = np.unique(gridscore["i_fold"])
    numblocks = len(CV) * len(folds)
    paramNames = gridscore.dtype.fields.keys()
    paramNames.remove("mean_validation_score")
    paramNames.remove("std")
    paramNames.remove("i_CV")
    paramNames.remove("i_fold")
    score = gridscore["mean_validation_score"]
    std = gridscore["std"]
    newgridscore = gridscore[paramNames]
    num_params = len(paramNames)
    ### get index of hit ###
    hitindex = []
    n_iter = len(score)/numblocks
    for k in range(numblocks):
        hit0index = np.argmax(score[k*n_iter: (k+1)*n_iter])
        hitindex.append(k*n_iter+hit0index )

    for m in range(num_params-1):
        i = paramNames[m]
        x = newgridscore[i]
        for n in range(m+1, num_params):
        # for j in list(set(paramNames)- set([i])):
            j = paramNames[n]
            y = newgridscore[j]
            compound = [x,y]
            # Only plot heat map if dtype of all elements of x, y are int or float
            if [True]* len(compound)== map(lambda t: np.issubdtype(t.dtype,  np.float) or np.issubdtype(t.dtype, np.int), compound):
                gridsize = 50
                fig = pl.figure()
                points = np.vstack([x,y]).T
                #####Construct MeshGrids##########
                xnew = np.linspace(max(x), min(x), gridsize)
                ynew = np.linspace(max(y), min(y), gridsize)
                X, Y = np.meshgrid(xnew, ynew)
                #####Interpolate Z on top of MeshGrids#######
                Z = griddata(points, score, (X, Y), method = "cubic", tol = 1e-2)
                z_min = min(score)
                z_max = max(score)
                pl.pcolormesh(X,Y,Z, cmap='RdBu', vmin=z_min, vmax=z_max)
                pl.axis([x.min(), x.max(), y.min(), y.max()])
                pl.xlabel(i, fontsize = 30)
                pl.ylabel(j, fontsize = 30)
                cb = pl.colorbar()
                cb.set_label(metric, fontsize = 30)
                ##### Mark the "hit" points #######
                hitx = x[hitindex]
                hity = y[hitindex]
                pl.plot(hitx, hity, 'rx')
                # Save the plot
                save_path = plot_path +obj+'/'+ clfName +'_' +metric+'_'+ i +'_'+ j+'.pdf'
                fig.savefig(save_path)

    def __call__(self, n):
        if len(self.f.shape) == 3:
            # f = f[x,v,t], 2 dim in phase space
            ft = self.f[n,:,:]
            pylab.pcolormesh(self.X, self.V, ft.T, cmap = 'jet')
            pylab.colorbar()
            pylab.clim(0,0.38) # for Landau test case
            pylab.grid()
            pylab.axis([self.xmin, self.xmax, self.ymin, self.ymax])
            pylab.xlabel('$x$', fontsize = 18)
            pylab.ylabel('$v$', fontsize = 18)
            pylab.title('$N_x$ = %d, $N_v$ = %d, $t$ = %2.1f' % (self.x.N, self.v.N, self.it*self.t.width))
            pylab.savefig(self.path + self.filename)
            pylab.clf()
            return None

        if len(self.f.shape) == 2:
            # f = f[x], 1 dim in phase space
            ft = self.f[n,:]
            pylab.plot(self.x.gridvalues,ft,'ob')
            pylab.grid()
            pylab.axis([self.xmin, self.xmax, self.ymin, self.ymax])
            pylab.xlabel('$x$', fontsize = 18)
            pylab.ylabel('$f(x)$', fontsize = 18)
            pylab.savefig(self.path + self.filename)
            return None

def plotLL(fname='out4.npy'):
    plt.figure()
    h= np.linspace(0,1,21)
    g= np.linspace(0,1,21)
    m=np.linspace(0,2,21)
    d=np.linspace(0,2,21)
    out=np.load(fname)
    print np.nanmax(out),np.nanmin(out)
    rang=np.nanmax(out)-np.nanmin(out)
    maxloc= np.squeeze(np.array((np.nanmax(out)==out).nonzero()))
    H,G=np.meshgrid(h,g)
    print maxloc
    for mm in range(m.size/2):
        for dd in range(d.size/2):
            plt.subplot(10,10,(9-mm)*10+dd+1)
            plt.pcolormesh(h,g,out[:,:,mm*2,dd*2].T,
                           vmax=np.nanmax(out),vmin=np.nanmax(out)-rang/4.)
            plt.gca().set_xticks([])
            plt.gca().set_yticks([])
            if mm==maxloc[2]/2 and dd==maxloc[3]/2:
                plt.plot(h[maxloc[0]],g[maxloc[1]],'ow',ms=8)
            if dd==0:
                print mm,dd
                plt.ylabel('%.1f'%m[mm*2])
            if mm==0: plt.xlabel('%.1f'%d[dd*2])
    plt.title(fname[:6])

def plot_iris_knn():
    iris = datasets.load_iris()
    X = iris.data[:, :2]  # we only take the first two features. We could
                        # avoid this ugly slicing by using a two-dim dataset
    y = iris.target

    knn = neighbors.KNeighborsClassifier(n_neighbors=3)
    knn.fit(X, y)

    x_min, x_max = X[:, 0].min() - .1, X[:, 0].max() + .1
    y_min, y_max = X[:, 1].min() - .1, X[:, 1].max() + .1
    xx, yy = np.meshgrid(np.linspace(x_min, x_max, 100),
                         np.linspace(y_min, y_max, 100))
    Z = knn.predict(np.c_[xx.ravel(), yy.ravel()])

    # Put the result into a color plot
    Z = Z.reshape(xx.shape)
    pl.figure()
    pl.pcolormesh(xx, yy, Z, cmap=cmap_light)

    # Plot also the training points
    pl.scatter(X[:, 0], X[:, 1], c=y, cmap=cmap_bold)
    pl.xlabel('sepal length (cm)')
    pl.ylabel('sepal width (cm)')
    pl.axis('tight')

 def outside_juristiction(self, data):
     data = copy.deepcopy(data)
     xynames = ['X', 'Y']
     data['Prediction'] = pandas.Series(self.clf.predict(data[xynames]))
     data = data[data.PdDistrictInt != data.Prediction]
     data = data.reset_index(drop=True)
     data['Prob'] = pandas.Series(
         [max(x) for x in self.clf.predict_proba(data[xynames])]
     )
     data = data[data.Prob > 0.95].reset_index(drop=True)
     fig = pl.figure()
     x_min, x_max = data.X.min() - 0.01, data.X.max() + 0.01
     y_min, y_max = data.Y.min() - 0.01, data.Y.max() + 0.01
     xx, yy = np.meshgrid(
         np.arange(x_min, x_max, self.step),
         np.arange(y_min, y_max, self.step)
     )
     Z = self.clf.predict(np.c_[xx.ravel(), yy.ravel()])
     Z = Z.reshape(xx.shape)
     pl.xlim(xx.min(), xx.max())
     pl.ylim(yy.min(), yy.max())
     pl.pcolormesh(xx, yy, Z, cmap='jet', alpha=1.0)
     pl.scatter(data.X, data.Y, c=data.PdDistrictInt, cmap='jet', alpha=1.0)
     pl.savefig('plots/knn_PD_2.png')
     pl.close(fig)
     return data

 def plotter(view=None, **kw):
     x,y = kw[args[0]],kw[args[1]]
     r = linspace(range[0],range[1],200)
     X,Y = meshgrid(x+r,y+r)
     kw['x'],kw['y'] = X,Y
     pylab.pcolormesh(x+r,y+r,nllf(**kw))
     pylab.plot(x,y,'o',hold=True, markersize=6,markerfacecolor='red',markeredgecolor='black',markeredgewidth=1, alpha=0.7)

def plot_knn_boundary():
    ## Training dataset preparation
    # use sklearn iris dataset
    iris_dataset = datasets.load_iris()


    # first two dimensions as the features
    # it's easy to plot boundary in 2D
    train_data = iris_dataset.data[:,:2] 
    print "init:",train_data

    # get labels
    labels = iris_dataset.target # labels
    print "init2:",labels

    
    ## Test dataset preparation
    h = 0.1
    
    x0_min = train_data[:,0].min() - 0.5
    x0_max = train_data[:,0].max() + 0.5
    
    x1_min = train_data[:,1].min() - 0.5
    x1_max = train_data[:,1].max() + 0.5
    
    x0_features, x1_features = np.meshgrid(np.arange(x0_min, x0_max, h), 
                                           np.arange(x1_min, x1_max, h))
    
    # test dataset are samples from the whole regions of feature domains
    test_data = np.c_[x0_features.ravel(), x1_features.ravel()]
    
    ## KNN classification
    p_labels = []   # prediction labels
    for test_sample in test_data:
        # knn prediction
        p_label = knn_predict(train_data, labels, test_sample, n_neighbors = 6)
        p_labels.append(p_label)
    
    # list to array
    p_labels = np.array(p_labels)
    p_labels = p_labels.reshape(x0_features.shape)
    
    ## Boundary plotting  边界策划
    pl.figure(1)
    pl.set_cmap(pl.cm.Paired)
    pl.pcolormesh(x0_features, x1_features, p_labels)
    
    pl.scatter(train_data[:,0], train_data[:,1], c = labels)
    # x y轴的名称
    pl.xlabel('feature 0')
    pl.ylabel('feature 1')

    # 设置x,y轴的上下限
    pl.xlim(x0_features.min(), x0_features.max())
    pl.ylim(x1_features.min(), x1_features.max())
    # 设置x,y轴记号
    pl.xticks(())
    pl.yticks(())
    
    pl.show()

def plotTiming(vortex_prob, vortex_times, times, map, grid_spacing, tornado_track, title, file_name, obs=None, centers=None, min_prob=0.1):
    nx, ny = vortex_prob.shape
    gs_x, gs_y = grid_spacing

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

    time_color_map = matplotlib.cm.Accent
    time_color_map.set_under('#ffffff')

    vortex_times = np.where(vortex_prob >= min_prob, vortex_times, -1)

    track_xs, track_ys = map(*reversed(tornado_track))

    pylab.figure(figsize=(10, 8))
    pylab.axes((0.025, 0.025, 0.95, 0.925))

    pylab.pcolormesh(xs, ys, vortex_times, cmap=time_color_map, vmin=times.min(), vmax=times.max())

    tick_labels = [ (datetime(2009, 6, 5, 18, 0, 0) + timedelta(seconds=int(t))).strftime("%H%M") for t in times ]
    bar = pylab.colorbar()#orientation='horizontal', aspect=40)
    bar.locator = FixedLocator(times)
    bar.formatter = FixedFormatter(tick_labels)
    bar.update_ticks()

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

    drawPolitical(map, scale_len=(xs[-1, -1] - xs[0, 0]) / 10.)

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