def visibility(vis,pol):
    """
    Helper function for plotMP. Defines the visibility/transparency of scatter
    symbols.

    @param vis: keyword among 'Visible', 'Invisible', and 'Transparent'.
    @type vis: string
    """

    if vis == "Visible":
        c = [pol]
        cmap = P.get_cmap('jet')
        v = True
        col = "0.0"
    elif vis == "Invisible": 
        c = [0.0]
        cmap = P.get_cmap('gist_yarg')
        v = False
        col = "1.0"
    elif vis == "Transparent": 
        c = [0.20]
        cmap = P.get_cmap('gist_yarg')
        v = True
        col = "0.80"
    else:
        c = None
        cmap = None
        v = None
        col = None
    
    return [c, cmap, v, col]

def show_matrix_with_names(matrix, vert_names, horiz_names, colormap='b_jet', overlay=None, force_range=False):

    def full_rename(namelist, subtypes):

        def rename(list, subtype, character):
            ret_list = []
            for bearer in list:
                modded = False
                for mod_decider in subtype:
                    if mod_decider in bearer:
                        modded = True
                        ret_list.append(bearer+' '+str(character))
                        break
                if not modded:
                    ret_list.append(bearer)
            return ret_list

        new_vert_names = namelist
        for idux, subtype in enumerate(subtypes, 1):
            new_vert_names = rename(new_vert_names, subtype, ''.join(['*']*idux))

        return new_vert_names

    if colormap == 'b_jet':
        prism_cmap = get_cmap('jet')
        prism_vals = prism_cmap(np.arange(0, 1, 0.01))
        prism_vals[99] = [0, 0, 0, 1]
        costum_cmap = colors.LinearSegmentedColormap.from_list('my_colormap', prism_vals)

    else:
        costum_cmap = get_cmap(colormap)

    if force_range:
        plt.imshow(matrix, interpolation='nearest', cmap=costum_cmap, vmin=force_range[0], vmax=force_range[1])
    else:
        plt.imshow(matrix, interpolation='nearest', cmap=costum_cmap)

    plt.colorbar()
    if overlay:
        print overlay[0]
        print overlay[1]
        overlay_y, overlay_x = np.nonzero(overlay[0])
        plt.scatter(overlay_x, overlay_y, c='k', marker='*', label=overlay[1])
        plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3, mode="expand", borderaxespad=0.)
    plt.tick_params(axis='both', labelsize=10)

    if type(vert_names) == list:
        plt.yticks(range(0, len(vert_names)), vert_names, rotation='horizontal')
    if type(horiz_names) == list:
        plt.xticks(range(0, len(horiz_names)), horiz_names, rotation='vertical')
    if type(vert_names) == tuple:
        vert_names = full_rename(vert_names[0], vert_names[1:])
        plt.yticks(range(0, len(vert_names)), vert_names, rotation='horizontal')
    if type(horiz_names) == tuple:
        horiz_names = full_rename(horiz_names[0], horiz_names[1:])
        plt.xticks(range(0, len(horiz_names)), horiz_names, rotation='vertical')


    plt.subplots_adjust(left=0.2, bottom=0.2)
    plt.show()

def compare(hdf5_in, hdf5_out, pdf_out, n_events):
    """Plots original and cropped img side by side."""
    _, original = load_dataset(hdf5_in)   # original dataset
    _, cropped = load_dataset(hdf5_out)  # cropped dataset

    # labels
    pid_ticks = [0, 1, 2, 3, 4, 5, 6, 7]
    pid_labels = ['nth', 'EM', 'mu', 'pi+', 'pi-', 'n', 'p', 'oth']

    with PdfPages(pdf_out) as pdf:
        for _ in range(n_events):
            # grab random event
            i = np.random.randint(0, cropped.shape[0])
            
            fig, axes = plt.subplots(1, 2)

            cmap = 'tab10'
            cbt = 'pid'

            fig.suptitle("Event: " + str(i))

            # plot pid image
            im = axes[0].imshow(original[i][0], cmap=pylab.get_cmap(cmap),
                                vmin=0, vmax=7)
            # plot pid image
            im = axes[1].imshow(cropped[i], cmap=pylab.get_cmap(cmap),
                                vmin=0, vmax=7)

            # just plot settings
            cbar = pylab.colorbar(im, fraction=0.04, ticks=pid_ticks)
            cbar.ax.set_yticklabels(pid_labels)
            cbar.set_label("pid", size=9)
            cbar.ax.tick_params(labelsize=6)

            pdf.savefig()

def plot_stats(X,Y,model,costs):
	#two plots, the decision fcn and points and the cost over time
	y_onehot = Trainer.class_to_onehot(Y)
	f,(p1,p2) = plot.subplots(1,2)
	p2.plot(range(len(costs)),costs)
	p2.set_title("Cost over time")
	
	#plot points/centroids/decision fcn
	cls_ct = y_onehot.shape[1]
	y_cls = Trainer.onehot_to_int(y_onehot)
	colors = get_cmap("RdYlGn")(np.linspace(0,1,cls_ct))
	
	#model_cents = model.c.get_value()
	#p1.scatter(model_cents[:,0], model_cents[:,1], c='black', s=81)
	for curclass,curcolor in zip(range(cls_ct),colors):
		inds = [i for i,yi in enumerate(y_cls) if yi==curclass]
		p1.scatter(X[inds,0], X[inds,1], c=curcolor)
		
	nx,ny = 200, 200
	x = np.linspace(X[:,0].min()-1,X[:,0].max()+1,nx)
	y = np.linspace(X[:,1].min()-1,X[:,1].max()+1,ny)
	xv,yv = np.meshgrid(x,y)
	
	Z = np.array([z for z in np.c_[xv.ravel(), yv.ravel()]])
	Zp = Trainer.onehot_to_int(np.array(model.probability(Z)))
	Zp = Zp.reshape(xv.shape)
	p1.imshow(Zp, interpolation='nearest', 
				extent=(xv.min(), xv.max(), yv.min(), yv.max()),
				origin = 'lower', cmap=get_cmap("Set1"))
	
	p1.set_title("Decision boundaries and centroids")
	f.tight_layout()
	plot.show()					

def plotDensityMap(bins,binnedData,xscale='log',yscale='log',normaliseX=True,logScale=True):
    if logScale:
        l, b, r, t = 0.1, 0.12, 1.0, 0.970
    else:
        l, b, r, t = 0.1, 0.12, 1.05, 0.970
    axes_rect = [l,b,r-l,t-b]
    fig=p.figure()
    fig.subplots_adjust(left=0.01, bottom=.05, right=.985, top=.95, wspace=.005, hspace=.05)
    ax = fig.add_axes(axes_rect)
    ax.set_xscale(xscale)
    ax.set_yscale(yscale)
    X,Y=bins.edge_grids
    if normaliseX:
        ySum=p.sum(binnedData,1)
        ySum=p.ma.masked_array(ySum,ySum==0)
        Z=binnedData.transpose()/ySum
    else:
        Z=binnedData.transpose()

    if logScale:
        mappable=ax.pcolor(X,Y,p.ma.array(Z,mask=Z==0),cmap=p.get_cmap("jet"),norm=matplotlib.colors.LogNorm())
    else:
        mappable=ax.pcolor(X,Y,p.ma.array(Z,mask=Z==0),cmap=p.get_cmap("jet"))

    markersize=5.0
    linewidth=2.0

    fig.colorbar(mappable)

    #ax.set_ylim((T/2.0/bins2D.centers[0][-1],T/2.0*2))
    #ax.set_xlim(bins2D.centers[0][0]*day,bins2D.centers[0][-1]*day)
    return fig,ax

	def create_plot_2d_speeches(self, withLabels = True):
		if withLabels:
			font = { 'fontname':'Tahoma', 'fontsize':0.5, 'verticalalignment': 'top', 'horizontalalignment':'center' }
			labels= [ (i['who'], i['date']) for i in self.data ]
			pylab.subplots_adjust(bottom =0.1)
			pylab.scatter(self.speech_2d[:,0], self.speech_2d[:,1], marker = '.' ,cmap = pylab.get_cmap('Spectral'))
			for label, x, y in zip(labels, self.speech_2d[:, 0], self.speech_2d[:, 1]):
			    pylab.annotate(
			        label, 
			        xy = (x, y), xytext = None,
			        ha = 'right', va = 'bottom', **font)
			        #,textcoords = 'offset points',bbox = dict(boxstyle = 'round,pad=0.5', fc = 'yellow', alpha = 0.5),
			        #arrowprops = dict(arrowstyle = '->', connectionstyle = 'arc3,rad=0'))

			pylab.title('U.S. Presidential Speeches(1790-2006)')
			pylab.xlabel('X')	
			pylab.ylabel('Y')

			pylab.savefig('plot_with_labels', bbox_inches ='tight', dpi = 1000, orientation = 'landscape', papertype = 'a0')
		else:
			pylab.subplots_adjust(bottom =0.1)
			pylab.scatter(self.speech_2d[:,0], self.speech_2d[:,1], marker = 'o' ,cmap = pylab.get_cmap('Spectral'))
			pylab.title('U.S. Presidential Speeches(1790-2006)')
			pylab.xlabel('X')	
			pylab.ylabel('Y')
			pylab.savefig('plot_without_labels', bbox_inches ='tight', dpi = 1000, orientation = 'landscape', papertype = 'a0')
		pylab.close()

    def __call__(self, inputs):
        from numpy import outer, arange, ones
        from pylab import figure, axis, imshow, title, show, subplot, text, clf, subplots_adjust
        maps = self.get_input('colormap')
        a=outer(arange(0,1,0.01),ones(10))

        if self.get_input('showall') is True:
            figure(figsize=(10,5))
            clf()
            l = len(tools.cmaps)
            subplots_adjust(top=0.9,bottom=0.05,left=0.01,right=0.99)
            for index, m in enumerate(tools.cmaps):
                #print index
                subplot(int(l/2)+l%2+1, 2, index+1)
                #print int(l/2)+l%2, 2, (index+1)/2+(index+1)%2+1
                axis("off")
                imshow(a.transpose(),aspect='auto',cmap=tools.get_cmap(m),origin="lower")
                #title(m,rotation=0,fontsize=10)
                text(0.5,0.5, m)
            show()
        elif self.get_input('show') is True:
            figure(figsize=(10,5))
            clf()
            axis("off")
            imshow(a.transpose(),aspect='auto',cmap=get_cmap(maps),origin="lower")
            title(maps,rotation=0,fontsize=10)
            show()
        from pylab import get_cmap
        res = get_cmap(maps)
        return res

def run_main(options):
    try:
        inp = open(options.input_file,'rb')
        geo_input = pickle.load(inp)
        inp.close()
    except:
        print 'Problems in opening the %s. Check the file.' % options.input_file
        exit()

    width,height = map(int,options.size.split(','))
    mode_colormap = options.mode
    n_connections = options.n_connections

    #Create image
    im = Image.new('RGBA', (width, height), (0, 0, 0, 0))
    draw = ImageDraw.Draw(im)
    geo_final = {}

    #Remove the same start == end and repeated cities (start,end) == (end,start).
    for origem,destino in geo_input:
        if origem != destino:
            if not geo_final.has_key((origem,destino)):
                if not geo_final.has_key((destino,origem)):
                    geo_final.setdefault((origem,destino),0)
                    geo_final[(origem,destino)] += geo_input[(origem,destino)]
                else:
                    geo_final.setdefault((destino,origem),0)
                    geo_final[(destino,origem)] += geo_input[(origem,destino)]

    matrix_geo = {}
    #Find all valid cities and calculate the distances
    max_distance = -10000000
    for origem,destino in geo_final:
        if geo_final[(origem,destino)] > n_connections:
            x1,y1 = geo2pixel(origem[0],origem[1],width,height)
            x2,y2 = geo2pixel(destino[0],destino[1],width,height)
            distance = sqrt((y2-y1)**2 + (x2-x1)**2)*40000/360.0
            if distance > max_distance:  max_distance = distance
            matrix_geo[(x1,y1),(x2,y2)] = distance

    #Sort by distance.
    sorted_matrix = sorted(matrix_geo.items(),key= lambda x: x[1], reverse=True)

    cmap =  get_cmap('jet')  if mode_colormap == 'all' else get_cmap('jet')
    dist_bins = logspace(1,log(max_distance)/log(10),255)


    #Plot using the correct colormap
    for ((x1,y1),(x2,y2)),distance in sorted_matrix:
        for i in range(len(dist_bins)-1):
            if distance > dist_bins[i] and distance < dist_bins[i+1]:
                break
        p = rgb2hex(converter.to_rgb(cmap(255-i)))
        draw.line([(y1,x1),(y2,x2)], p )

    #Save the image.
    im.save(options.output_file)

def flamePlot(mapZ, noNanZ, markRegions, flatten, colorCeil):
	print '\nPlotting flameview'

	ratios=[4,1,4]
	gs = gridspec.GridSpec(len(ratios), 1, height_ratios=ratios)
	plt.figure(figsize=(24,9))


	# Direct
	plt.subplot(gs[0])

	plt.imshow(mapZ, cmap=get_cmap(fColorMap), interpolation='nearest', aspect='auto')
	plt.clim(-colorCeil,colorCeil)

	plt.title('Overview')
	plt.xlabel('Probe number')
	plt.ylabel('Window Step')


	# Markers
	plt.subplot(gs[1])

	flatMarkers = [0]*len(markRegions[0])
	for i, val in enumerate(markRegions[0]):
		if abs(val) > flatten:
			flatMarkers[i] = val
	markRegions.append(flatMarkers)

	plt.imshow(markRegions, cmap=get_cmap(fColorMap), interpolation='nearest', aspect='auto')
	plt.clim(-colorCeil,colorCeil)

	plt.title('Maximum values')
	plt.xlabel('Probe number')
	plt.ylabel('Post|Pre')

	# Flattened
	for i in range(len(mapZ)):
		for j in range(len(noNanZ)):
			stouff=mapZ[i][j]
			if abs(stouff) < flatten:
				stouff = 0
			mapZ[i][j]=stouff

	plt.subplot(gs[2])
	plt.imshow(mapZ, cmap=get_cmap(fColorMap), interpolation='nearest', aspect='auto')
	plt.clim(-colorCeil,colorCeil)

	plt.title('Filtered at abs(z) > ' + str(flatten))
	plt.xlabel('Probe number')
	plt.ylabel('Window Step')

	plt.tight_layout()
	return plt

 def _get__pl(self):
     '''return the LinearSegmentedColormap describing this CustomColormap'''
     if self.cmap=='file' and self.fname:
         colors = lut_manager.parse_lut_file(self.fname)
         if self.reverse:
             colors.reverse()
         return LinearSegmentedColormap.from_list('file',colors)
     elif self.cmap=='custom_heat':
         return gen_heatmap(t=self.threshold,reverse=self.reverse)
     elif self.reverse:
         return get_cmap(self.cmap+'_r')
     else:
         return get_cmap(self.cmap)

def getColorMapPlotSettings(trackFile, colorPropertyName, settings):
	if 'colorSettings' in settings:
		colorSettings = settings['colorSettings']
		vmin = colorSettings['vmin']
		vmax = colorSettings['vmax']
		separation = colorSettings['separation']
		colormap = P.get_cmap(colorSettings['colorMap'])
	else:
		vmin = (trackFile.axisLimits[colorPropertyName])[0]
		vmax = (trackFile.axisLimits[colorPropertyName])[1]
		separation = 35
		colormap = P.get_cmap('jet')

	return vmin, vmax, separation, colormap

def embed_dat_matrix_two_dimensions(low_dimension_data_matrix,
                                    y=None,
                                    labels=None,
                                    density_colormap='Blues',
                                    instance_colormap='YlOrRd'):
    from sklearn.preprocessing import scale
    low_dimension_data_matrix = scale(low_dimension_data_matrix)
    # make mesh
    x_min, x_max = low_dimension_data_matrix[:, 0].min(), low_dimension_data_matrix[:, 0].max()
    y_min, y_max = low_dimension_data_matrix[:, 1].min(), low_dimension_data_matrix[:, 1].max()
    step_num = 50
    h = min((x_max - x_min) / step_num, (y_max - y_min) / step_num)  # step size in the mesh
    b = h * 10  # border size
    x_min, x_max = low_dimension_data_matrix[:, 0].min() - b, low_dimension_data_matrix[:, 0].max() + b
    y_min, y_max = low_dimension_data_matrix[:, 1].min() - b, low_dimension_data_matrix[:, 1].max() + b
    xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))

    # induce a one class model to estimate densities
    from sklearn.svm import OneClassSVM
    gamma = max(x_max - x_min, y_max - y_min)
    clf = OneClassSVM(gamma=gamma, nu=0.1)
    clf.fit(low_dimension_data_matrix)

    # Plot the decision boundary. For that, we will assign a color to each
    # point in the mesh [x_min, m_max] . [y_min, y_max].
    if hasattr(clf, "decision_function"):
        score_matrix = clf.decision_function(np.c_[xx.ravel(), yy.ravel()])
    else:
        score_matrix = clf.predict_proba(np.c_[xx.ravel(), yy.ravel()])[:, 1]
    # Put the result into a color plot
    levels = np.linspace(min(score_matrix), max(score_matrix), 40)
    score_matrix = score_matrix.reshape(xx.shape)

    if y is None:
        y = 'white'

    plt.contourf(xx, yy, score_matrix, cmap=plt.get_cmap(density_colormap), alpha=0.9, levels=levels)
    plt.scatter(low_dimension_data_matrix[:, 0], low_dimension_data_matrix[:, 1],
                alpha=.5,
                s=70,
                edgecolors='gray',
                c=y,
                cmap=plt.get_cmap(instance_colormap))
    # labels
    if labels is not None:
        for id in range(low_dimension_data_matrix.shape[0]):
            label = labels[id]
            x = low_dimension_data_matrix[id, 0]
            y = low_dimension_data_matrix[id, 1]
            plt.annotate(label, xy=(x, y), xytext=(0, 0), textcoords='offset points')

    def plot_vapor(self,imagename='vapor.png', title='',timestep=0):
	"""plot  qvapor
	Optional input:
	     --------
	     imagename = 'your file name.png', defaults to vapor.png'
	     title = 'your title' defaults to none
	     timestep = integer. If there are multiple time periods per file then
	     choose your timestep, defaults to 0'

	    Returns:
	    --------
		contour plot of qvapor
	
	"""

	pl.figure(2)
	pl.clf()
	try:
	    vapor_full=self.variable_dict['QVAPOR']
	except:
	    vapor_full=self.get_var('QVAPOR')

	qvapor=vapor_full[timestep,0,:,:].copy()


	custom_map=pl.get_cmap('jet',lut=22)
	try:
	    m=self.map_proj
	    x=self.map_x
	    y=self.map_y
	except:
	    self.bmap()
	    m=self.map_proj
	    x=self.map_x
	    y=self.map_y

	
	custom_map=pl.get_cmap('jet',lut=10)
	contour_levs=n.arange(0.012,0.025, 0.001)
	vapor_plot=m.contourf(x,y,qvapor[:,:],contour_levs,cmap=custom_map, title=title)#,cmap=custom_map)
	vapor_plot.set_clim((0.012,0.025))

	self.map_lines()

	pl.colorbar(orientation='horizontal')
	pl.title(title)

	pl.savefig(self.plot_directory+'/'+imagename)

def ST_pi_pure(r,steps = 99):
	"""
	Makes a heat map over ST space for a given value of r, using the pure strategy model, plotting
	fitness.

	Inputs
	======
	r : float
		Value of relatedness
	steps : int {99}
		Number of points to sample S and T at.

	"""

	Ss = np.linspace(-1,2,steps)
	Ts = np.linspace(0,3,steps)
	dataFull = np.zeros( (steps,steps) )
	for i,S in enumerate(Ss):
		for j,T in enumerate(Ts):
			x = ESS_pure(r,S,T)
			dataFull[i,j] = fit_pure(x,r,S,T)/maximal_possible_fitness(S,T)

	#pl.figure()
	cmap = pl.get_cmap('Reds')
	pl.imshow(dataFull, origin = [Ts[0], Ss[0]], interpolation = 'nearest',\
			extent = [ Ts[0],Ts[-1],Ss[0],Ss[-1] ], cmap = cmap, vmin = 0, vmax = .5*(Ss[-1]+Ts[-1]))
	pl.plot( [ Ts[0], Ts[-1] ],[ 0,0 ],color='black', linewidth = 2.5 )
	pl.plot( [ 1, 1 ],[ Ss[0],Ss[-1] ],'--',color='black',  linewidth = 2.5 )
	pl.plot( [ 0, 3 ],[ 2, -1 ],'--',color='black',  linewidth = 2.5 )

def makeimg(wav):
	global callpath
	global imgpath

	fs, frames = wavfile.read(os.path.join(callpath, wav))
	
	pylab.ion()

	# generate specgram
	pylab.figure(1)
	
	# generate specgram
	pylab.specgram(
		frames,
		NFFT=256, 
		Fs=22050, 
		detrend=pylab.detrend_none,
		window=numpy.hamming(256),
		noverlap=192,
		cmap=pylab.get_cmap('Greys'))
	
	x_width = len(frames)/fs
	
	pylab.ylim([0,11025])
	pylab.xlim([0,round(x_width,3)-0.006])
	
	img_path = os.path.join(imgpath, wav.replace(".wav",".png"))

	pylab.savefig(img_path)
	
	return img_path

def ST_beta(r,steps = 99):
	"""
	Makes a heat map of beta over ST space, for a given value of r.
	Uses the full model.
	Steps : int {99}
		Number of steps to sample S and T at
	Returns
	=======
	None
	"""

	Ss = np.linspace(-1,3,steps)
	Ts = np.linspace(0,4,steps)
	dataFull = np.zeros( (steps,steps) )
	for i,S in enumerate(Ss):
		for j,T in enumerate(Ts):
			mod = ESS_class(r,S,T)
			state = mod.getESS()
			dataFull[i,j] = mod.beta(state,r)

	#pl.figure()
	cmap = pl.get_cmap('coolwarm')
	pl.imshow(dataFull, origin = [Ts[0], Ss[0]],\
			extent = [ Ts[0],Ts[-1],Ss[0],Ss[-1] ], vmin = -1,vmax = 1, cmap = cmap)
	pl.plot( [ Ts[0], Ts[-1] ],[ 0,0 ],color='black', linewidth = 2.5 )
	pl.plot( [ 1, 1 ],[ Ss[0],Ss[-1] ],'--',color='black',  linewidth = 2.5 )
	pl.plot( [ 0, 3 ],[ 2, -1 ],'--',color='black',  linewidth = 2.5 )

def ST_xA(r,steps = 99):
	"""
	Makes a heat map for x over ST space for a given r using the full model.
	inputs
	======
	r : float
		Value of relatedness
	steps: int {99}
		Number of decrete points at which to sample S and T

	"""


	Ss = np.linspace(-1,2,steps)
	Ts = np.linspace(0,3,steps)
	dataFull = np.zeros( (steps,steps) )
	for i,S in enumerate(Ss):
		for j,T in enumerate(Ts):
			mod = ESS_class(r,S,T)
			state = mod.getESS()
			dataFull[i,j] = mod.xA(state,r)

	#pl.figure()
	cmap = pl.get_cmap('Reds')
	pl.imshow(dataFull, origin = [Ts[0], Ss[0]],\
			extent = [ Ts[0],Ts[-1],Ss[0],Ss[-1] ], vmin = 0,vmax = 1, cmap = cmap)
	pl.plot( [ Ts[0], Ts[-1] ],[ 0,0 ],color='black', linewidth = 2.5 )
	pl.plot( [ 1, 1 ],[ Ss[0],Ss[-1] ],'--',color='black',  linewidth = 2.5 )
	pl.plot( [ 0, 3 ],[ 2, -1 ],'--',color='black',  linewidth = 2.5 )

def cluster_show(_3D_chan1, _3D_chan2, _3D_labels, n_clusters_, means, std, hulls):
    """

    :param _3D_chan1:
    :param _3D_chan2:
    :param _3D_labels:
    """
    red = (1.0, 0.0, 0.0)
    green = (0.0, 1.0, 0.1)

    mlab.pipeline.volume(mlab.pipeline.scalar_field(_3D_chan1), color=green)
    mlab.pipeline.volume(mlab.pipeline.scalar_field(_3D_chan2), color=red)
    mlab.show()

    cm = get_cmap('gist_rainbow')

    for i in range(0, n_clusters_):
        re_img = np.zeros(_3D_chan2.shape)
        re_img[_3D_labels == i] = _3D_chan2[_3D_labels == i]
        mlab.pipeline.volume(mlab.pipeline.scalar_field(re_img), color=tuple(cm(1.*i/n_clusters_)[:-1]))

    # mean_arr = np.zeros((n_clusters_, 3))
    # std_arr = np.zeros((n_clusters_))
    # for i in range(0, n_clusters_):
    #     mean_arr[i, :] = means[i]
    #     std_arr[i] = std[i]
    # x,y,z = mean_arr.T
    # mlab.points3d(x,y,z, std_arr)

    for hull in hulls:
        x,y,z = hull.points.T
        triangles = hull.simplices
        mlab.triangular_mesh(x, y, z, triangles, representation='wireframe', color=(0, 0, 0))

    mlab.show()

def color_iterator(grouping, cmap='jet', n=None):
    '''
    Given a Matplotlib colormap, iterate through the color range in equal-sized
    increments based on the size of the group.

    Parameters
    ----------
    grouping : iterable
        Values to return on each iteration
    cmap : string
        Matplotlib colormap to use
    n : int
        Size of the group.  If not provided, will attempt to estimate the size
        from the iterable.
    '''

    # Attempt to get the length of the iterator first.  If we can't get the
    # length, then  we need to convert the iterator into a list so we can get
    # the number of elements.
    if n is None:
        try:
            n = len(grouping)
        except:
            grouping = list(grouping)
            n = len(grouping)

    if isinstance(cmap, basestring):
        cmap = pylab.get_cmap(cmap)
    for i, g in enumerate(grouping):
        if n == 1:
            yield cmap(1), g
        else:
            yield cmap(i/(n-1)), g

    def _plot_annotation_on_ax(self, signal, ax, autoscale=False, colourmap="flag"):

        if autoscale:
            xstart = 0
            xdelta = signal.duration.total_seconds()
        else:

            xstart,ystart,xdelta,ydelta = ax.viewLim.bounds

        if xstart <0:
            start_time = timedelta()
        else:
            start_time = timedelta(seconds=xstart)

        stop_time = timedelta(seconds=xdelta) +  timedelta(seconds=xstart)
        sub_sig = signal[start_time:stop_time]
        xs =np.linspace(0, sub_sig.duration.total_seconds() ,sub_sig.size) + start_time.total_seconds()
        ys = sub_sig.values
        probs = sub_sig.probas

        ys = ys.reshape((1,ys.size))

        zoom_f = float(self.max_point_amplitude_plot)/ sub_sig.size

        ys = zoom(ys,[1, zoom_f], order=0)


        ax.imshow(ys, extent=[np.min(xs), np.max(xs), 1.5, -0.5], aspect="auto",
                  cmap=colourmap, vmin=0, vmax=255, origin='lower')


        ax.plot(xs,probs,"-", color="k", linewidth=3)
        ax.plot(xs,probs,"-", color="y", linewidth=1,alpha=0.5)


        jet = cm = pl.get_cmap(colourmap)
        cNorm  = colors.Normalize(vmin=0, vmax=255)
        scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=jet)

        states = np.unique(ys)

        boxes = [pl.Rectangle((0, 0), 1, 1, fc=scalarMap.to_rgba(col)) for col in states]
        labels = [chr(s) for  s in states]
        pl.legend(boxes,labels, loc='lower right')

        n_labels = 8 #fixme magic number

        if len(xs) > n_labels:
            trimming = int(float(len(xs)) / float(n_labels))
            xs_trimmed = np.round(xs[::trimming])
        else:
            xs_trimmed = xs

        time_strings = [str(timedelta(seconds=s)) for s in xs_trimmed]


        ax.set_xticks(xs_trimmed)
        ax.set_xticklabels(time_strings, rotation=70)

        return