def marking_init(enhanced_img, mep, mbp):
# mark initially extracted minutiae points

    img_thin = np.array(enhanced_img[:])  # convert image to array for marking points

    fig = plt.figure(figsize=(10,8),dpi=30000)

    num1 = len(mep)
    num2 = len(mbp)

    figure, imshow(img_thin, cmap = cm.Greys_r)
    title('mark extracted points')
    plt.hold(True)

    for i in range(num1):
        xy = mep[i,:]
        u = xy[0]
        v = xy[1]
        if (u != 0.0) & (v != 0.0):
            plt.plot(v, u, 'r.', markersize = 7)

    plt.hold(True)

    for i in range(num2):
        xy = mbp[i,:]
        u = xy[0]
        v = xy[1]
        if (u != 0.0) & (v != 0.0):
            plt.plot(v, u, 'c+', markersize = 7)

    plt.show()
    cv2.imwrite("initial_extraction.png", img_thin)

def plotFit(min_x, max_x, mu, sigma, theta, p):
    #PLOTFIT Plots a learned polynomial regression fit over an existing figure.
    #Also works with linear regression.
    #   PLOTFIT(min_x, max_x, mu, sigma, theta, p) plots the learned polynomial
    #   fit with power p and feature normalization (mu, sigma).

    # Hold on to the current figure
    plt.hold(True)

    # We plot a range slightly bigger than the min and max values to get
    # an idea of how the fit will vary outside the range of the data points
    x = np.array(np.arange(min_x - 15, max_x + 25, 0.05)) # 1D vector

    # Map the X values 
    X_poly = pf.polyFeatures(x, p)
    X_poly = X_poly - mu
    X_poly = X_poly/sigma

    # Add ones
    X_poly = np.column_stack((np.ones((x.shape[0],1)), X_poly))

    # Plot
    plt.plot(x, np.dot(X_poly, theta), '--', linewidth=2)

    # Hold off to the current figure
    plt.hold(False)

def visualize(u1, t1, u2, t2, U, omega):
    plt.figure(1)
    plt.plot(t1, u1, 'r--o')
    t_fine = np.linspace(0, t1[-1], 1001)  # мелкая сетка для точного решения
    u_e = u_exact(t_fine, U, omega)
    plt.hold('on')
    plt.plot(t_fine, u_e, 'b-')
    plt.legend([u'приближенное', u'точное'], loc='upper left')
    plt.xlabel('$t$')
    plt.ylabel('$u$')
    tau = t1[1] - t1[0]
    plt.title('$\\tau = $ %g' % tau)
    umin = 1.2*u1.min();
    umax = -umin
    plt.axis([t1[0], t1[-1], umin, umax])
    plt.savefig('tmp1.png');  plt.savefig('tmp1.pdf')
    plt.figure(2)
    plt.plot(t2, u2, 'r--o')
    t_fine = np.linspace(0, t2[-1], 1001)  # мелкая сетка для точного решения
    u_e = u_exact(t_fine, U, omega)
    plt.hold('on')
    plt.plot(t_fine, u_e, 'b-')
    plt.legend([u'приближенное', u'точное'], loc='upper left')
    plt.xlabel('$t$')
    plt.ylabel('$u$')
    tau = t2[1] - t2[0]
    plt.title('$\\tau = $ %g' % tau)
    umin = 1.2 * u2.min();
    umax = -umin
    plt.axis([t2[0], t2[-1], umin, umax])
    plt.savefig('tmp2.png');
    plt.savefig('tmp2.pdf')

def smooth_demo():

    Smoother = SmoothClass()



    xfile = ArrayClass("DatasetX")
    #print xfile.array
    xn = np.array(xfile.array)

    plt.subplot(211)
    #plt.plot(np.ones(ws))

    windows=['blackman']
    #windows=['flat', 'hanning', 'hamming', 'bartlett', 'blackman']

    plt.hold(True)


    plt.axis([0,30,0,1.1])

    plt.legend(windows)
    plt.title("The smoothing windows")
    plt.subplot(212)
    #plt.plot(xn)
    plt.plot(xn)
    plt.plot(Smoother.smooth(xn,10,'blackman'))

 def test_prop(self):
     N = 800.0
     V = linspace(5.0,51.0,50)
     rho = 1.2255
     beta = 45.0
     J    = list()
     CT   = list()
     CP   = list()
     effy = list()
     for v in V:
         data = self.analyze_prop(beta,N,v,rho)
         J.append(data[2])
         CT.append(data[3])
         CP.append(data[4])
         effy.append(data[5])
     plt.figure(1)
     plt.grid(True)
     plt.hold(True)
     plt.plot(J,CT,'o-')
     plt.xlabel('J')
     plt.plot(J,CP,'ro-')
     plt.axis([0,2.5,0,0.15])
     plt.figure(2)
     plt.plot(J,effy,'gs-')
     plt.hold(True)
     plt.grid(True)
     plt.axis([0,2.5,0,1.0])
     plt.xlabel('advance ratio')
     plt.ylabel('efficiency')
     plt.show()

def plot_categorical_scatter_with_mean(vals, categoryLabels, jitter=True, colours=None, xlabel=None, ylabel=None, title=None):
    import matplotlib.colors
    import scipy.stats
    import pdb
    numCategories = len(vals)
    plt.hold(True)
    if colours is None:
        colours = plt.cm.gist_rainbow(np.linspace(0,1,numCategories))
    for category in range(numCategories):
        edgeColour = matplotlib.colors.colorConverter.to_rgba(colours[category], alpha=0.5)
        xval = (category+1)*np.ones(len(vals[category]))
        if jitter:
            jitterAmt = np.random.random(len(xval))
            xval = xval + (0.3 * jitterAmt) - 0.15
        #pdb.set_trace()
        plt.plot(xval, vals[category], 'o', mec=edgeColour, mew = 4, mfc='none', ms=16)
        mean = np.mean(vals[category])
        sem = scipy.stats.sem(vals[category])
        print mean, sem
        plt.plot(category+1, mean, 'o', color='k', mec=colours[category], ms=20)
        plt.errorbar(category+1, mean, yerr = sem, color=colours[category])
    plt.xlim(0,numCategories+1)
    plt.ylim(0,1)
    ax = plt.gca()
    ax.set_xticks(range(1,numCategories+1))
    ax.set_xticklabels(categoryLabels, fontsize=16)
    if xlabel is not None:
        plt.xlabel(xlabel, fontsize=20)
    if ylabel is not None:
        plt.ylabel(ylabel, fontsize=20)
    if title is not None:
        plt.title(title)
    plt.show()

def createResponsePlot(dataframe,plotdir):
    mag = dataframe['MAGPDE'].as_matrix()
    response = (dataframe['TFIRSTPUB'].as_matrix())/60.0
    response[response > 60] = 60 #anything over 60 minutes capped at 6 minutes
    imag5 = (mag >= 5.0).nonzero()[0]
    imag55 = (mag >= 5.5).nonzero()[0]
    fig = plt.figure(figsize=(8,6))
    n,bins,patches = plt.hist(response[imag5],color='g',bins=60,range=(0,60))
    plt.hold(True)
    plt.hist(response[imag55],color='b',bins=60,range=(0,60))
    plt.xlabel('Response Time (min)')
    plt.ylabel('Number of earthquakes')
    plt.xticks(np.arange(0,65,5))
    ymax = text.ceilToNearest(max(n),10)
    yinc = ymax/10
    plt.yticks(np.arange(0,ymax+yinc,yinc))
    plt.grid(True,which='both')
    plt.hold(True)
    x = [20,20]
    y = [0,ymax]
    plt.plot(x,y,'r',linewidth=2,zorder=10)
    s1 = 'Magnitude 5.0, Events = %i' % (len(imag5))
    s2 = 'Magnitude 5.5, Events = %i' % (len(imag55))
    plt.text(35,.85*ymax,s1,color='g')
    plt.text(35,.75*ymax,s2,color='b')
    plt.savefig(os.path.join(plotdir,'response.pdf'))
    plt.savefig(os.path.join(plotdir,'response.png'))
    plt.close()
    print 'Saving response.pdf'

def plot_meth_and_twobeds(coverage, methylated, mod):
    l1 = len(mod.bed_list_gt)
    l2 = len(mod.bed_list_h)
    n_cells = np.shape(coverage)[0]
    plt.figure()
    # Get current size

    fig_size_temp = plt.rcParams["figure.figsize"]
    fig_size = fig_size_temp
    fig_size[0] = 500
    fig_size[1] = 40
    plt.rcParams["figure.figsize"] = fig_size

    fig, axarr = plt.subplots(n_cells+l1+l2+1, 1, sharex=True)
    plt.hold(True)
    plot_meth(axarr[:n_cells], coverage, methylated)

    for i in range(0, l1):
        axn = n_cells+i
        plot_bed([axarr[axn]], [mod.bed_list_gt[i]])
        axarr[axn].set_ylabel(mod.state_name_gt[i])

    for i in range(0, l2):
        axn = n_cells+l1+1+i
        plot_bed([axarr[axn]], [mod.bed_list_h[i]])
        axarr[axn].set_ylabel(mod.state_name_h[i])

    fig.savefig(mod.path_name + mod.bed_title + 'l1 = ' + str(l1) + 'l2 = ' + str(l2) + 'n_cells = ' + str(n_cells)+'_l='+str(mod.l)+'_l_test='+str(mod.l_test) + '.pdf')
    plt.hold(False)
    plt.rcParams["figure.figsize"] = fig_size_temp
    plt.close(fig)

    def plot_filter_characteristics(self):
        w, h = freqz(self.freq_filter.num, self.freq_filter.denom)
        plt.figure(1)
        plt.subplot(2,1,1)
        plt.hold(True)
        powa = plt.plot((self.filter_parameters.sample_rate*0.5/pi)*w, abs(h),'b-', label = 'Char. amplitudowa')
        plt.title('Charakterystyki filtru')
        plt.xlabel('Czestotliwosc [Hz]')
        plt.ylabel('Amplituda')


        plt.twinx(ax=None)
        angles = unwrap(angle(h))
        plt.znie = plot((self.filter_parameters.sample_rate*0.5/pi)*w,angles, 'g-', label = 'Char. fazowa')
        plt.ylabel('Faza')

        plt.grid()
        tekst = powa + znie
        wybierz = [l.get_label() for l in tekst]

        plt.legend(tekst, wybierz, loc='best')
    ########################################################################################################################
        plt.subplot(2,1,2)

        w2, gd = group_delay((num, denom))

        plt.plot((sample_rate*0.5/pi)*w2, gd)
        plt.grid()
        plt.xlabel('Czestotliwosc [Hz]')
        plt.ylabel('Opoznienie grupowe [probki]')
        plt.title('Opoznienie grupowe filtru')

        plt.show()

def plot(x, dict_res, plot_func):
    colors = ['#e41a1c', '#377eb8', '#4daf4a', '#984ea3']

    plt.hold(True)
    k = 0
    for (alg, dtype) in dict_res.keys():
        for comp in dict_res[(alg, dtype)]:
            if len(alg) < 4:
                label = '{0:4s}'.format(alg.upper()) + ' - '
            else:
                label = '{0:4s}'.format(alg.upper()) + ' - '
            if comp:
                label += '{0:5s} - '.format('comp.')
                linestyle = '-'
            else:
                label += '{0:5s} - '.format('QR')
                linestyle = '--'
            label += dtype

            y = dict_res[(alg, dtype)][comp]
            valid = np.isfinite(y).flatten()
            if not np.any(valid):
                continue
            plot_func(x[valid], y[valid], label=label,
                      linestyle=linestyle, linewidth=2,
                      marker='o', markeredgecolor='none', color=colors[k])
        k += 1

    plt.hold(False)

def band_select(spikeTimeStamps, eventOnsetTimes, amplitudes, bandwidths, timeRange, fullRange = [0.0, 2.0]):
    numBands = np.unique(bandwidths)
    numAmps = np.unique(amplitudes)
    spikeArray = np.zeros((len(numBands), len(numAmps)))
    errorArray = np.zeros_like(spikeArray)
    trialsEachCond = behavioranalysis.find_trials_each_combination(bandwidths, 
                                                                   numBands, 
                                                                   amplitudes, 
                                                                   numAmps)
    spikeTimesFromEventOnset, trialIndexForEachSpike, indexLimitsEachTrial = spikesanalysis.eventlocked_spiketimes(
                                                                                                        spikeTimeStamps, 
                                                                                                        eventOnsetTimes,
                                                                                                        fullRange)
    spikeCountMat = spikesanalysis.spiketimes_to_spikecounts(spikeTimesFromEventOnset, indexLimitsEachTrial, timeRange)
    baseTimeRange = [timeRange[1]+0.5, fullRange[1]]
    baseSpikeCountMat = spikesanalysis.spiketimes_to_spikecounts(spikeTimesFromEventOnset, indexLimitsEachTrial, baseTimeRange)
    baselineSpikeRate = np.mean(baseSpikeCountMat)/(baseTimeRange[1]-baseTimeRange[0])
    plt.hold(True)
    for amp in range(len(numAmps)):
        trialsThisAmp = trialsEachCond[:,:,amp]
        for band in range(len(numBands)):
            trialsThisBand = trialsThisAmp[:,band]
            if spikeCountMat.shape[0] != len(trialsThisBand):
                spikeCountMat = spikeCountMat[:-1,:]
                print "FIXME: Using bad hack to make event onset times equal number of trials"
            thisBandCounts = spikeCountMat[trialsThisBand].flatten()
            spikeArray[band, amp] = np.mean(thisBandCounts)
            errorArray[band,amp] = stats.sem(thisBandCounts)
    return spikeArray, errorArray, baselineSpikeRate

def showResults(challenger_data, model):
    ''' Show the original data, and the resulting logit-fit'''
    
    temperature = challenger_data[:,0]
    failures = challenger_data[:,1]
    
    # First plot the original data
    plt.figure()
    setFonts()
    sns.set_style('darkgrid')
    np.set_printoptions(precision=3, suppress=True)
    
    plt.scatter(temperature, failures, s=200, color="k", alpha=0.5)
    plt.yticks([0, 1])
    plt.ylabel("Damage Incident?")
    plt.xlabel("Outside Temperature [F]")
    plt.title("Defects of the Space Shuttle O-Rings vs temperature")
    plt.tight_layout
    
    # Plot the fit
    x = np.arange(50, 85)
    alpha = model.params[0]
    beta = model.params[1]
    y = logistic(x, beta, alpha)
    
    plt.hold(True)
    plt.plot(x,y,'r')
    plt.xlim([50, 85])
    
    outFile = 'ChallengerPlain.png'
    showData(outFile)

def iterer_los_plot(x0,max_error,max_iter,f):
	# For aa visualisere iterasjonsprosessen tar vi vare paa punktene underveis
	X=array(x0)
	F=array(0.0)
	# Setter startverdi for feilen lik 1.0 slik at while-lokken garantert starter
	error = 1.0
	j = 0
	while error > max_error and j < max_iter:
		x1=f(x0)
		X=append(X,[x0,x0])
		F=append(F,[x0,x1])
		error = abs(x1-x0)
		x0 = x1
		j = j+1
	# Gi beskjed dersom vi har brukt opp max antall iterasjoner
	# uten at feilen har blitt mindre enn max_error
	if j == max_iter:
		print('Max antall iterasjoner brukt opp.\n')
	else:
		print('Losning funnet: x=%3.12f' % (x0))
	# Plotter x, f(x) og iterasjonsprosessen(stiplet)
	x = arange(0,1,0.001)
	figure('Iterasjoner')
	hold(True)
	plot(x,f(x),'b',label=r'$f(x)$')
	plot(X,F,'k--',label='iterasjoner')
	plot(x,x,'g',label=r'$y=x$')
	legend(loc='best')
	show()
	# Returner losningen
	return x0

def write_final(newmep, newmbp, enhanced_img):

    eminutiae_array = removezero(newmep)
    bminutiae_array = removezero(newmbp)

    num1 = len(eminutiae_array)
    num2 = len(bminutiae_array)

    img_thin = np.array(enhanced_img[:])
    fig = plt.figure(figsize=(15,12),dpi=30000)
    figure, imshow(img_thin, cmap = cm.Greys_r)
    title('minutiae marking')
    plt.hold(True)

    for i in range(num1):
        xy = eminutiae_array[i,:]
        u = xy[0]
        v = xy[1]
        plt.plot(v, u, 'r.', markersize = 10)

    plt.hold(True)

    for i in range(num2):
        xy = bminutiae_array[i,:]
        u = xy[0]
        v = xy[1]
        plt.plot(v, u, 'c+', markersize = 15)

    plt.show()
    cv2.imwrite("final_minutiae_image.png", img_thin)

 def internal_surf_afteraxes(cd):
     plt.hold(True)
     plt.title('')
     plt.ylabel('m')
     plt.subplots_adjust(hspace=0.05)        
     plt.plot([multilayer_data.bathy_location,multilayer_data.bathy_location],bottom_surf_zoomed,'--k')
     plt.hold(False)

def plot_conv(all_JSD,all_JSDs,rest_type):
    fold = len(all_JSD)
    rounds = len(all_JSDs[0])
    n_rest = len(rest_type)
    new_JSD = [[] for i in range(n_rest)]
    for i in range(len(all_JSD)):
        for j in range(n_rest):
            new_JSD[j].append(all_JSD[i][j])
    JSD_dist = [[] for i in range(n_rest)]
    JSD_std = [[] for i in range(n_rest)]
    for rest in range(n_rest):
        for f in range(fold):
            temp_JSD = all_JSDs[f][:,rest]
            JSD_dist[rest].append(np.mean(temp_JSD))
            JSD_std[rest].append(np.std(temp_JSD))
    plt.figure(figsize=(10,5*n_rest))
    x = np.arange(100./fold,101.,fold)
    colors = ['red','blue','green','black','magenta','gold','navy']
    for i in range(n_rest):
        plt.subplot(n_rest,1,i+1)
        plt.plot(x,new_JSD[i],'o-',color=colors[i],label=rest_type[i])
        plt.hold(True)
        plt.plot(x,JSD_dist[i],'o',color=colors[i],label=rest_type[i])
        plt.fill_between(x,np.array(JSD_dist[i])+np.array(JSD_std[i]),np.array(JSD_dist[i])-np.array(JSD_std[i]),color=colors[i],alpha=0.2)
        plt.xlabel('dataset (%)')
        plt.ylabel('JSD')
        plt.legend(loc='best')
    plt.tight_layout()
    plt.savefig('convergence.pdf')

    def update_figures(self):
        plt.figure(self.figure.number)
        x = np.arange(0, 256, 0.1)  # artificial x-axis
        # self.figure.gca().cla()  # clearing the figure, just to be sure

        # plt.subplot(411)
        plt.plot(self.bins, self.hist, 'k')
        plt.hold(True)
        if self.rv_healthy and self.rv_hypo and self.rv_hyper:
            healthy_y = self.rv_healthy.pdf(x)
            if self.win.params['unaries_as_cdf']:
                hypo_y = (1 - self.rv_hypo.cdf(x)) * self.rv_healthy.pdf(self.rv_healthy.mean())
                hyper_y = self.rv_hyper.cdf(x) * self.rv_healthy.pdf(self.rv_healthy.mean())
            else:
                hypo_y = self.rv_hypo.pdf(x)
                hyper_y = self.rv_hyper.pdf(x)
            y_max = max(healthy_y.max(), hypo_y.max(), hyper_y.max())
            fac = self.hist.max() / y_max

            plt.plot(x, fac * healthy_y, 'g', linewidth=2)
            plt.plot(x, fac * hypo_y, 'b', linewidth=2)
            plt.plot(x, fac * hyper_y, 'r', linewidth=2)
            plt.title('all PDFs')
        ax = plt.axis()
        plt.axis([0, 256, ax[2], ax[3]])
        plt.hold(False)

        self.canvas.draw()

def plotPaper(appeared,
              citedBy,
              pubLabel):
    hold(True)

    yrange = range(appeared, curDate + 1)
    months = len(yrange)
    cites = [0] * months
    citeVenues = {}

    for citation in citedBy:
        
        (venue, date) = citation
        
        if venue in citeVenues:
            citeVenues[venue] += 1
        else:
            citeVenues[venue] = 1
            
        for i in range(date - appeared, months):
            cites[i] += 1
            
        for i in range(date - startDate, curDate - startDate + 1):
            citations[i] += 1

    citeAx.plot(yrange,
                cites,
                label = pubLabel)

    logCiteAx.semilogy(yrange,
                       cites,
                       label = pubLabel)

    venues[pubLabel] = citeVenues
    articleCitations[pubLabel] = (appeared, cites)

    def gauge_after_axes(cd):

        if cd.gaugeno in [1,2,3,4]:
            axes = plt.gca()
            # # Add Kennedy gauge data
            # kennedy_gauge = kennedy_gauges[gauge_name_trans[cd.gaugeno]]
            # axes.plot(kennedy_gauge['t'] - seconds2days(date2seconds(gauge_landfall[0])), 
            #          kennedy_gauge['mean_water'] + kennedy_gauge['depth'], 'k-', 
            #          label='Gauge Data')

            # Add GeoClaw gauge data
            geoclaw_gauge = cd.gaugesoln
            axes.plot(seconds2days(geoclaw_gauge.t - date2seconds(gauge_landfall[1])),
                  geoclaw_gauge.q[3,:] + gauge_surface_offset[0], 'b--', 
                  label="GeoClaw")

            # Add ADCIRC gauge data
            # ADCIRC_gauge = ADCIRC_gauges[kennedy_gauge['gauge_no']]
            # axes.plot(seconds2days(ADCIRC_gauge[:,0] - gauge_landfall[2]), 
            #          ADCIRC_gauge[:,1] + gauge_surface_offset[1], 'r-.', label="ADCIRC")

            # Fix up plot
            axes.set_title('Station %s' % cd.gaugeno)
            axes.set_xlabel('Days relative to landfall')
            axes.set_ylabel('Surface (m)')
            axes.set_xlim([-2,1])
            axes.set_ylim([-1,5])
            axes.set_xticks([-2,-1,0,1])
            axes.set_xticklabels([r"$-2$",r"$-1$",r"$0$",r"$1$"])
            axes.grid(True)
            axes.legend()

            plt.hold(False)

def PlotEDepSummary(gFiles,nFiles,figureName='EDepSummary.png',tParse=GetThickness,
  histKey='eDepHist'):
  """ PlotEDepSummary
  Plotss the energy deposition summary
  """
  # Extrating the average values
  gT = list()
  gDep = list()
  gDepError = list()
  nT = list()
  nDep = list()
  nDepError = list()
  for fname in gFiles:
    f = TFile(fname,'r')
    hist = f.Get(histKey)
    gT.append(GetThickness(fname))
    gDep.append(hist.GetMean())
    gDepError.append(hist.GetMeanError())
  for fname in nFiles:
    f = TFile(fname,'r')
    hist = f.Get(histKey)
    nT.append(GetThickness(fname))
    nDep.append(hist.GetMean())
    nDepError.append(hist.GetMeanError())
  # Plotting
  plt.errorbar(gT,gDep,yerr=gDepError,fmt='r+')
  plt.hold(True)
  plt.errorbar(nT,nDep,yerr=nDepError,fmt='go')
  plt.xlabel("Thickness (mm)")
  plt.ylabel("Average Energy Deposition (MeV)")
  plt.legend(["Co-60","Cf-252"])
  plt.xscale("log")
  plt.yscale("log")
  plt.grid(True)
  plt.savefig(figureName)