def standard_process():
    """
    standard process from scikit learn
    """
    X = np.atleast_2d([1., 3., 5., 6., 7., 8.]).T
    y = f(X).ravel()
    x = np.atleast_2d(np.linspace(0, 10, 1000)).T
    gp = gaussian_process.GaussianProcess(theta0=1e-2, thetaL=1e-4, thetaU=1e-1)
    gp.fit(X, y)
    y_pred, MSE = gp.predict(x, eval_MSE=True)
    sigma = np.sqrt(MSE)

    # Plot the function, the prediction and the 95% confidence interval based on
    # the MSE
    fig = plt.figure()
    plt.plot(x, f(x), 'r:', label=u'$f(x) = x\,\sin(x)$')
    plt.plot(X, y, 'r.', markersize=10, label=u'Observations')
    plt.plot(x, y_pred, 'b-', label=u'Prediction')
    plt.fill(np.concatenate([x, x[::-1]]),
            np.concatenate([y_pred - 1.9600 * sigma,
                           (y_pred + 1.9600 * sigma)[::-1]]),
            alpha=.5, fc='b', ec='None', label='95% confidence interval')
    plt.xlabel('$x$')
    plt.ylabel('$f(x)$')
    plt.ylim(-10, 20)
    plt.legend(loc='upper left')
    plt.show()

def plot_data_and_line(w1,w2):
    w1,w2 = float(w1),float(w2)
    if w2 != 0 :
        y1,y2 = (-w1*(xl1))/w2, (-w1*(xl2))/w2
        vx1,vy1 = [xl1,xl2,xl2,xl1,xl1], [y1,y2,yl2,yl2,y1]
        vx2,vy2 = [xl1,xl2,xl2,xl1,xl1], [y1,y2,yl1,yl1,y1]
    elif w1 != 0:
        vx1,vy1 = [xl2,0,0,xl2,xl2], [yl1,yl1,yl2,yl2,yl1]
        vx2,vy2 = [xl1,0,0,xl1,xl1], [yl1,yl1,yl2,yl2,yl1]
    else:
        print "ERROR, Invalid w1 and w2."
        return;
    if  w2 > 0 or ( w2 == 0 and w1 > 0):
        c1,c2 = 'b','r'
    else:
        c1,c2 = 'r','b'
    fig = plt.figure()
    ax = fig.add_subplot(111)
    plt.scatter(X1[Y > 0], X2[Y > 0], s = 80, c = 'b', marker = "o")
    plt.scatter(X1[Y<= 0], X2[Y<= 0], s = 80, c = 'r', marker = "^")
    plt.fill(vx1, vy1, c1, alpha = 0.25)
    plt.fill(vx2, vy2, c2, alpha = 0.25)
    ax.set_title(("w1 = %s, w2 = %s")%( w1, w2))
    ax.set_xlim(xl1, xl2)
    ax.set_ylim(yl1, yl2)
    fig.set_size_inches(6, 6)
    plt.show()    

def plot_triangle(x, y, color=None, draw_lines=True, fill=False):
    """plot an rgb triangle in xy

    Args:
        x,y (numpy.array): [r, g, b] coords

    kwargs:
        color (str): if none, 1st point--> red, 2nd --> green, 3rd -> blue

        drawLines (bool): draw outline

        fill (bool): fill triangle

    """
    if fill:
        plt.fill(x, y, color='grey', alpha='0.5')
    if draw_lines:
        indexVal = np.hstack([np.arange(x.size), 0])
        plt.plot(x[indexVal], y[indexVal], '-k')

    if color:
        plt.plot(x[0], y[0], 'o', x[1], y[1], 'o', x[2], y[2], 'o',
                 color=color)
    else:
        plt.plot(x[0], y[0], 'or', x[1], y[1], 'og', x[2], y[2], 'ob')

def wiggle (traces, skipt=1,scale=1.,lwidth=.1,offsets=None,redvel=0., manthifts=None, tshift=0.,sampr=1.,clip=10., dx=1., color='black',fill=True,line=True, ax=None):

  ns = traces.shape[1]
  ntr = traces.shape[0]
  t = np.arange(ns)*sampr
  timereduce = lambda offsets, redvel, shift: [float(offset) / redvel + shift for offset in offsets]

  if (offsets is not None):
    shifts = timereduce(offsets, redvel, tshift)
  elif (manthifts is not None):
    shifts = manthifts
  else:
    shifts = np.zeros((ntr,))

  for i in range(0, ntr, skipt):
    trace = traces[i].copy()
    trace[0] = 0
    trace[-1] = 0

    if ax == None:
      if (line):
        plt.plot(i*dx + clipsign(trace / scale, clip), t - shifts[i], color=color, linewidth=lwidth)
      if (fill):
        for j in range(ns):
          if (trace[j] < 0):
            trace[j] = 0
        plt.fill(i*dx + clipsign(trace / scale, clip), t - shifts[i], color=color, linewidth=0)
    else:
      if (line):
        ax.plot(i*dx + clipsign(trace / scale, clip), t - shifts[i], color=color, linewidth=lwidth)
      if (fill):
        for j in range(ns):
          if (trace[j] < 0):
            trace[j] = 0
        ax.fill(i*dx + clipsign(trace / scale, clip), t - shifts[i], color=color, linewidth=0)

 def dist_load():
     
             # add distributed load
     plt.subplot(3,1,1)
     for k in np.linspace(0,L,20):
         ax1.arrow(k, 11+L/10, 0, -3, head_width=L*0.01, head_length=L*0.1, fc='k', ec='k')
     plt.title('Free Body Diagram')
     plt.axis('off') # removes axis and labels
     #ax1.set_yticklabels('') 
     
     # dist load shear
     x = [0,0,L,L]
     y = [0,5,-5,0]
     plt.subplot(3,1,2)
     plt.ylabel('Shear, V')
     plt.title('Shear Diagram')
     plt.fill(x, y, 'b', alpha=0.25)
     plt.grid(True)
     plt.xlim([-1, L+1])
     
     # dist load bending
     x = np.linspace(-L/2,L/2,100)
     y = -(x**2)+(np.max(x**2))
     x = np.linspace(0,L,100)
     plt.subplot(3,1,3)
     plt.title('Bending Diagram')
     plt.ylabel('Moment, M')
     plt.fill(x, y, 'b', alpha=0.25)
     plt.grid(True)
     plt.xlim([-1, L+1])

def radar_plot():
    """
    radar plot
    """
    # 生成测试数据
    labels = np.array(["A组", "B组", "C组", "D组", "E组", "F组"])
    data = np.array([68, 83, 90, 77, 89, 73])
    theta = np.linspace(0, 2*np.pi, len(data), endpoint=False)

    # 数据预处理
    data = np.concatenate((data, [data[0]]))
    theta = np.concatenate((theta, [theta[0]]))

    # 画图方式
    plt.subplot(111, polar=True)
    plt.title("雷达图", fontproperties=myfont)

    # 设置"theta grid"/"radar grid"
    plt.thetagrids(theta*(180/np.pi), labels=labels, fontproperties=myfont)
    plt.rgrids(np.arange(20, 100, 20), labels=np.arange(20, 100, 20), angle=0)
    plt.ylim(0, 100)

    # 画雷达图,并填充雷达图内部区域
    plt.plot(theta, data, "bo-", linewidth=2)
    plt.fill(theta, data, color="red", alpha=0.25)

    # 图形显示
    plt.show()
    return

def hinton(W, max_weight=None):
    """ Draws a Hinton diagram for visualizing a weight matrix. 
    """
    if max_weight is None:
        max_weight = 2 ** np.ceil(np.log2(np.max(np.abs(W))))

    # Temporarily disable matplotlib interactive mode if it is on,
    # otherwise this takes forever.
    isinteractive = plt.isinteractive()
    if isinteractive:
        plt.ioff()

    height, width = W.shape
    plt.clf()
    plt.fill(np.array([0, width, width, 0]), np.array([0, 0, height, height]), "gray")

    plt.axis("off")
    plt.axis("equal")
    for (y, x), w in np.ndenumerate(W):
        if w != 0:
            color = "white" if w > 0 else "black"
            area = min(1, np.abs(w) / max_weight)
            _blob(x + 0.5, height - y - 0.5, area, color)

    if isinteractive:
        plt.ion()

def hist_with_peak(x, bins=None, range=None, ax=None, orientation='vertical',
                   histtype=None, **kwargs):
    if ax is None:
        ax = plt.gca()
    hist, bin_edges = np.histogram(x, bins=bins, range=range)
    hist_n = hist * 1.0/hist.max()
    width = bin_edges[1] - bin_edges[0]
    x = np.ravel(zip(bin_edges[:-1], bin_edges[:-1]+width))
    y = np.ravel(zip(hist_n, hist_n))
    if histtype == 'step':
        if orientation == 'vertical':
            plt.plot(x, y, **kwargs)
        elif orientation == 'horizontal':
            plt.plot(y, x, **kwargs)
        else:
            raise ValueError
    elif histtype == 'stepfilled':
        if orientation == 'vertical':
            plt.fill(x, y, **kwargs)
        elif orientation == 'horizontal':
            plt.fill(y, x, **kwargs)
        else:
            raise ValueError
    else:
        raise ValueError

def noiseless():
  X = np.atleast_2d([1.,3.,5.,6.,7.,8.,]).T
  # observations
  y = f(X).ravel() # reshape to 1-D
  # mesh the input space to cover all points to predict f(x) and the MSE
  x = np.atleast_2d(np.linspace(0,10,1000)).T # and flatten
  # instantiate gauss process
  gp = GaussianProcess(corr='cubic', theta0=1e-2, thetaL=1e-4,
                       thetaU=1e-1, random_start=100)
  # fit to data using maximum likelihood estimation of params
  gp.fit(X,y)
  # make predictions on meshed x-axis, also return MSE
  y_pred, MSE = gp.predict(x, eval_MSE=True)
  sigma = np.sqrt(MSE)
  
  # plot
  fig = plt.figure()
  plt.plot(x, f(x), 'r:', label=u'$f(x) = x\,\sin(x)$')
  plt.plot(X, y, 'r.', markersize=10, label=u'Observations')
  plt.plot(x, y_pred, 'b-', label=u'Prediction')
  # fill the space between the +/-MSE
  plt.fill(np.concatenate([x, x[::-1]]), # reverse order of x
           np.concatenate([y_pred - 1.9600 * sigma,
                          (y_pred + 1.9600 * sigma)[::-1]]),
           alpha=0.5, fc='b', ec='None', label='95% confidence interval')
           # shade, fill color, edge color
  plt.title('Noiseless case')
  plt.xlabel('$x$')
  plt.ylabel('$f(x)$')
  plt.ylim(-10, 20)
  plt.legend(loc='upper left')
  return

	def GPVisualize1D(self, locations, measurements, predict_range = (0, 1), num_samples = 1000):
		"""
		Visualize posterior in graphical form
		NOTE: very ineffecient since we are using the weight space view to vizualize this
		"""	

		# Grid points
		x = np.atleast_2d(np.linspace(predict_range[0], predict_range[1], num_samples, endpoint=False)).T

		# Compute predictions - very inefficient because we are using the weight space view
		predicted_mean = [0.0] * num_samples
		predicted_variance = [0.0] * num_samples
		for i in xrange(num_samples):
			predicted_mean[i] = self.GPMean(locations, measurements, x[i])[0]
			predicted_variance[i] = self.GPVariance2(locations, x[i])[0]

		# Plot posterior mean and variances
		pl.plot(x, self.GPRegressionTestEnvironment(x), 'r:', label=u'$f(x)$')
		pl.plot(locations, measurements, 'r.', markersize=10, label=u'Observations')
		pl.plot(x, predicted_mean, 'b-', label=u'Prediction')
		pl.fill(np.concatenate([x, x[::-1]]),
        np.concatenate([predicted_mean - 1.9600 * np.sqrt(predicted_variance),
                       (predicted_mean + 1.9600 * np.sqrt(predicted_variance))[::-1]]),
        				alpha=.5, fc='b', ec='None', label='95% confidence interval')
		pl.xlabel('$x$')
		pl.ylabel('$f(x)$')
		pl.legend(loc='upper left')

		pl.show()

def hinton(W, maxWeight=None):
    """
    Draws a Hinton diagram for visualizing a weight matrix.
    Temporarily disables matplotlib interactive mode if it is on,
    otherwise this takes forever.
    """
    reenable = False
    if plt.isinteractive():
        plt.ioff()
    plt.clf()
    height, width = W.shape
    if not maxWeight:
        maxWeight = 2**np.ceil(np.log(np.max(np.abs(W)))/np.log(2))

    plt.fill(np.array([0, width, width, 0]), np.array([0, 0, height, height]),
             'gray')
    plt.axis('off')
    plt.axis('equal')
    for x in xrange(width):
        for y in xrange(height):
            _x = x+1
            _y = y+1
            w = W[y, x]
            if w > 0:
                _blob(_x - 0.5, height - _y + 0.5,
                      min(1, w/maxWeight), 'white')
            elif w < 0:
                _blob(_x - 0.5, height - _y + 0.5,
                      min(1, -w/maxWeight), 'black')
    if reenable:
        plt.ion()
    plt.show()

def main():
    x = linspace(0, 1, 200)
    y = sin(4 * pi * x) * exp(-5 * x)

    fill(x, y, 'r')
    grid(True)
    show()

def convexview(linearcombination):
    # draw K
    plt.plot([1.0, 2.0, 20.0, 3.0, 1.0],[2.0, 4.0, 0.0, 1.0, 2.0],'b',lw=2.0)
    plt.fill([1.0, 2.0, 20.0, 3.0, 1.0],[2.0, 4.0, 0.0, 1.0, 2.0],
             facecolor='b', alpha=0.2, edgecolor='none')
    plt.text(3.5,1.2,r'$K$',fontsize=28,color='b')
    # draw axes
    plt.plot([0.5, 0.5],[-0.5, 4.5],'k',lw=1.0)
    plt.plot([0.0, 4.0],[0.5, 0.5],'k',lw=1.0)
    # plot u, v
    u = np.array([1.6,1.7])
    v = np.array([3.0,3.0])
    plt.plot(u[0],u[1],'.',markersize=14,color='k')
    plt.text(u[0]-0.3,u[1]-0.3,r'$u$',fontsize=24,color='k')
    plt.plot(v[0],v[1],'.',markersize=14,color='k')
    plt.text(v[0]+0.1,v[1]+0.1,r'$v$',fontsize=24,color='k')
    # show either interpretation
    eps = 0.6
    z = (1.0 - eps) * u + eps * v
    plt.plot(z[0],z[1],'.',markersize=14,color='k')
    if linearcombination:
        # lin-comb : (1 - eps) u + eps v
        plt.plot([u[0],v[0]], [u[1],v[1]],'k',lw=2.0)
        plt.text(z[0],z[1]-0.3,r'$(1-\varepsilon) u + \varepsilon v$',fontsize=24,color='k')
    else:
        # intoK : u + eps (v - u)
        ax = plt.axes()
        w = eps * (v - u)
        ax.arrow(u[0],u[1],w[0],w[1],
                 head_width=0.1, head_length=0.2, length_includes_head=True, fc='k', ec='k')
        plt.text(z[0],z[1]-0.3,r'$u + \varepsilon (v - u)$',fontsize=24,color='k')
    plt.axis([-1.0,4.0,0.0,5.0],'k',lw=2.0)
    plt.axis('off')

def __main__():
    '''test code'''
    # make up data points
    #np.random.seed(1234)
    points = np.random.rand(3, 2)

    # compute Voronoi tesselation
    vor = Voronoi(points)

    # plot
    regions, vertices = voronoi_finite_polygons_2d(vor,radius=100000)
    print("--")
    print(regions)
    print("--")
    print(vertices)

    # colorize
    for region in regions:
        polygon = vertices[region]
        plt.fill(*zip(*polygon), alpha=0.4)

    plt.plot(points[:,0], points[:,1], 'ko')
    plt.axis('equal')
    plt.xlim(vor.min_bound[0] - 0.1, vor.max_bound[0] + 0.1)
    plt.ylim(vor.min_bound[1] - 0.1, vor.max_bound[1] + 0.1)

    #plt.savefig('voro.png')
    plt.show()

    def test_voronoi(self):
        # make up data points
        np.random.seed(1234)
        points = np.random.rand(15, 2)

        # compute Voronoi tesselation
        vor = Voronoi(points)

        # plot
        regions, vertices = self.voronoi_finite_polygons_2d(vor)
        print "--"
        print regions
        print "--"
        print vertices

        # colorize
        for region in regions:
            polygon = vertices[region]
            plt.fill(*zip(*polygon), alpha=0.4)

        plt.plot(points[:, 0], points[:, 1], "ko")
        plt.xlim(vor.min_bound[0] - 0.1, vor.max_bound[0] + 0.1)
        plt.ylim(vor.min_bound[1] - 0.1, vor.max_bound[1] + 0.1)

        plt.show()

def histogram_gdal_direct(bag_file, patch_name,verbose):
    # This is the one 
    'Plot the depth histogram.  Use the tif until gdal 1.7.0 comes out'

    patch = osgeo.gdal.Open(bag_file)
    band = patch.GetRasterBand(1)
    bandmin,bandmax = band.ComputeRasterMinMax()
    hist = band.GetDefaultHistogram()

    #print hist
    #print '----------------------------------------------------------------------'
    #print hist[3]
    hist_vals = hist[3][:-1]
    while hist_vals[-1] == 0:
        hist_vals.pop()

    xticks = ['%02.1f' % depth for depth in  np.arange(hist[0],hist[1],(hist[1]-hist[0])/5)]
    xticks.append('%.1f' % hist[1])
    if verbose:
        print (xticks)  # Why is the 0.0 tick not showing?
    plt.xticks([val * len(hist_vals)/5 for val in range(len(xticks))],xticks) # Yuck!
    plt.fill(range(len(hist_vals)),hist_vals)
    plt.grid(True)
    plt.savefig(patch_name+'-hist.png') #,dpi=50)
    with file(patch_name+'.hist','w') as o:
        o.write('\n'.join([str(v) for v in hist_vals]))

    def plot_kde_overlap(self, terms, color1='#0067a2', color2='#e8a945', overlap_color='#dddddd', **kwargs):

        term1 = terms[0]
        term2 = terms[1]

        t1 = self.stem(term1)
        t2 = self.stem(term2)

        bc = self.score_braycurtis(t1, t2, **kwargs)

        kde1 = self.kde(t1, **kwargs)
        kde2 = self.kde(t2, **kwargs)
        plt.plot(kde1, color=color1, label=term1)
        plt.plot(kde2, color=color2, label=term2)

        overlap = np.minimum(kde1, kde2)
        plt.fill(overlap, color=overlap_color)
        plt.title(term1+', '+term2+' - '+str(round(bc, 4)))

        plt.xlabel('Word Offset')
        plt.ylabel('Number of Occurrences')
        plt.legend(loc='upper right')

        fig = plt.gcf()
        fig.set_size_inches(10, 4)
        fig.tight_layout()

        return plt

	def display(self):
		length = self.length
		width = self.width
		base_x = self.x - .5*length
		base_y = self.y - .5*width
		xs = [base_x, base_x + length, base_x+length, base_x, base_x]
		ys = [base_y, base_y, base_y+width, base_y+width, base_y]
		plt.plot(xs, ys, 'k')
		if not self.flyable:
			plt.fill(xs, ys, 'r')
		for letter in self.exitnodelist:
			if letter == 'C':
				x = base_x + .5*length
				y = base_y +.5*width
			if letter == 'N':
				x = base_x + .5*length
				y = base_y + width
			if letter == 'S':
				x = base_x + .5*length
				y = base_y
			if letter == 'E':
				x = base_x + length
				y = base_y + .5*width
			if letter == 'W':
				x = base_x
				y = base_y + .5*width
			plt.plot([x, base_x+.5*length],[y, base_y +.5*width], 'g')
			plt.plot([x, base_x+.5*length],[y, base_y +.5*width], 'go')

def main():
    axes = lsst.afw.geom.ellipses.Axes(4, 3, 1)
    ellipse = lsst.afw.geom.ellipses.Ellipse(
        axes, lsst.geom.Point2D(0.25338, 0.76032))
    region = lsst.afw.geom.ellipses.PixelRegion(ellipse)
    for bbox in [ellipse.computeBBox(), lsst.geom.Box2D(region.getBBox())]:
        corners = bbox.getCorners()
        pyplot.fill([p.getX() for p in corners], [p.getY()
                                                  for p in corners], alpha=0.2)
    envelope = region.getBBox()
    envelope.grow(2)
    ellipse.plot(alpha=0.2)
    ellX = []
    ellY = []
    allX, allY = numpy.meshgrid(
        numpy.arange(envelope.getBeginX(), envelope.getEndX()),
        numpy.arange(envelope.getBeginY(), envelope.getEndY())
    )
    gt = ellipse.getGridTransform()
    mgt = gt.getMatrix()
    transX = mgt[0, 0] * allX + mgt[0, 1] * allY + mgt[0, 2]
    transY = mgt[1, 0] * allX + mgt[1, 1] * allY + mgt[1, 2]
    allR = (transX**2 + transY**2)**0.5
    pyplot.plot(ellX, ellY, 'ro', markeredgewidth=0, alpha=0.5)
    pyplot.plot(allX[allR < 1], allY[allR < 1], '+')
    pyplot.plot(allX[allR > 1], allY[allR > 1], 'x')

def load_world():
	shpRecords = shpUtils.loadShapefile('world_borders/world_borders.shp')['features']
	colors = load_color_map()
	plt.figure(figsize=(16, 9))

	last = ''
	for i in range(0,len(shpRecords)):
		if shpRecords[i]['info']['CNTRY_NAME'] != last:
			print shpRecords[i]['info']['CNTRY_NAME']
			last = shpRecords[i]['info']['CNTRY_NAME']

		# x and y are empty lists to be populated with the coords of each geometry.
		x = []
		y = []
		#print shpRecords[i]
		for j in range(0,len(shpRecords[i]['shape']['parts'][0]['points'])):
			tempx = float(shpRecords[i]['shape']['parts'][0]['points'][j][0])
			tempy = float(shpRecords[i]['shape']['parts'][0]['points'][j][1])
			x.append(tempx)
			y.append(tempy) # Populate the lists  

		# Creates a polygon in matplotlib for each geometry in the shapefile
		if shpRecords[i]['info']['CNTRY_NAME'] in colors:
			if shpRecords[i]['info']['CNTRY_NAME'] == 'Congo' or shpRecords[i]['info']['CNTRY_NAME'] == 'Zaire':
				plt.fill(x,y, facecolor=colors['Democratic Republic of the Congo'])
			plt.fill(x,y, facecolor=colors[shpRecords[i]['info']['CNTRY_NAME']])

	plt.axis('equal')
	plt.savefig('world_calls.png', dpi=100, format='png')
	plt.show()