def stem_hist(img, plotTitle):
  import matplotlib.pyplot as plt
  from matplotlib import cm
  imgarray = img.calc_hist()
  plt.figure()
  plt.stem(arange(img.grayLevel), imgarray)
  plt.title(plotTitle)

def impz(b, a=1):
    """Plot step and impulse response of an FIR filter.

    b : float
        Forward terms of the FIR filter.
    a : float
        Feedback terms of the FIR filter. (Default value = 1)

    From http://mpastell.com/2010/01/18/fir-with-scipy/

    Returns
    -------
    None

    """
    l = len(b)
    impulse = np.repeat(0., l)
    impulse[0] = 1.
    x = np.arange(0, l)
    response = sp.lfilter(b, a, impulse)
    plt.subplot(211)
    plt.stem(x, response)
    plt.ylabel('Amplitude')
    plt.xlabel(r'n (samples)')
    plt.title(r'Impulse response')
    plt.subplot(212)
    step = sp.cumsum(response)
    plt.stem(x, step)
    plt.ylabel('Amplitude')
    plt.xlabel(r'n (samples)')
    plt.title(r'Step response')
    plt.subplots_adjust(hspace=0.5)

def plotstft(sxx, Fs=100):
    winlen = int(len(sxx[0]))

    # with plt.():
    fig1 = plt.figure()
    ax = fig1.add_subplot(1,1,1)
    ctr = int(winlen / 2)
    faxis = np.multiply(Fs / 2, np.linspace(0, 1, ctr))*60
    ratio = []
    for dft in sxx:
        mag = abs(dft[0:ctr])
        max_idx, max_val = max(enumerate(mag), key = lambda p: p[1])
        ptotal = np.square(np.linalg.norm(mag,2))
        pmax = np.square(np.absolute(max_val))
        # print('max power: {}'.format(max_val))
        frac = pmax/ptotal
        # print(frac)
        ratio.append(frac)

        ax.plot(faxis, mag, linewidth=3.0)
        ax.set_xlabel('Frequency (RPM)')
        ax.set_ylabel('|H(f)|')
        ax.set_title('User 4 STFT Spectrum')
        font = {'family' : 'sans-serif ',
                'weight' : 'bold',
                'size'   : 30}
        rc('font', **font)

    fig1.savefig('STFTPlot.png')
    fig2 = plt.figure()
    plt.stem(np.linspace(1,len(ratio),num=len(ratio)),ratio, linewidth=3.0)
    plt.xlabel('Window (10s)', fontsize = 60)
    plt.ylabel('Symmetry in Pedaling', fontsize = 60)
    fig2.savefig('RatioPlot.png')

def run_OMP(plot=False, **options):
    """Recover one signal using OMP."""
    n = options.pop('n', 128)
    k = options.pop('k', 5)
    m = options.pop('m', 20)
    dist = options.pop('dist', 'uniform')
    seed = options.pop('seed', None)
    return_locals = options.pop('return_locals', True)

    print ('Recovering signal wiht n=%(n)i, k=%(k)i and m=%(m)i using OMP' %
           locals())

    if seed:
        np.random.seed(seed)
    x = get_sparse_x(n, k, dist=dist)
    if seed:
        np.random.seed(seed + 198)
    A = random_dict(m, n)
    y = np.dot(A, x)

    x_hat, residues, scores, Delta = omp(A, y, save_data=True)

    if plot:
        plt.figure()
        plt.stem(range(n), x, 'r-', 'ro', 'k:')
        plt.stem(range(n), x_hat, 'b:', 'bx', 'k:')
        plt.show()

    print 'error', norm(x - x_hat.reshape(n, 1))

    if return_locals:
        return locals()

def _periodogram_plot(title, column, data, trend, peaks):
    """display periodogram results using matplotlib"""

    periods, power = periodogram(data)
    plt.figure(1)
    plt.subplot(311)
    plt.title(title)
    plt.plot(data, label=column)
    if trend is not None:
        plt.plot(trend, linewidth=3, label="broad trend")
        plt.legend()
        plt.subplot(312)
        plt.title("detrended")
        plt.plot(data - trend)
    else:
        plt.legend()
        plt.subplot(312)
        plt.title("(no detrending specified)")
    plt.subplot(313)
    plt.title("periodogram")
    plt.stem(periods, power)
    for peak in peaks:
        period, score, pmin, pmax = peak
        plt.axvline(period, linestyle='dashed', linewidth=2)
        plt.axvspan(pmin, pmax, alpha=0.2, color='b')
        plt.annotate("{}".format(period), (period, score * 0.8))
        plt.annotate("{}...{}".format(pmin, pmax), (pmin, score * 0.5))
    plt.tight_layout()
    plt.show()

    def plotImpulseResponse(self, xmin=None, xmax=None, ymin_imp=None, ymax_imp=None, ymin_step=None, ymax_step=None):
        """Plot the frequency and phase response of the filter object.

      :param xmin: Minimum value for x-axis.
      :param xmax: Maximum value for x-axis.
      :param ymin_imp: Minimum value for y-axis for the impulse response plot.
      :param ymax_imp: Maximum value for y-axis for the impulse response plot.
      :param ymin_step: Minimum value for y-axis for the step response plot.
      :param ymax_step: Maximum value for y-axis for the step response plot.
      """
        #  def plotImpulseResponse(b,a=1):
        l = len(self.ir)
        impulse = np.repeat(0.0, l)
        impulse[0] = 1.0
        x = np.arange(0, l)
        response = sp.signal.lfilter(self.ir, 1, impulse)
        mp.subplot(211)
        mp.stem(x, response)
        mp.ylabel("Amplitude")
        mp.xlabel(r"n (samples)")
        mp.title(r"Impulse response")
        mp.subplot(212)

        step = np.cumsum(response)
        mp.stem(x, step)
        mp.ylabel("Amplitude")
        mp.xlabel(r"n (samples)")
        mp.title(r"Step response")
        mp.subplots_adjust(hspace=0.5)
        mp.show()

def moment_ss_shear_bending(L,Pin,ain):
    '''
    Shear Bending plot of moment loads of a simply supported beam
    L = 4 # total length of beam
    Pin = [5]  # point moment load
    ain = [2]  # location of point load

    # or more multiple point moments
    L = 10
    Pin = [3,-15]
    ain = [2,6]
    '''

    import numpy as np
    import matplotlib.pyplot as plt
    
    x = np.arange(0,L,L*0.02)
    V = np.zeros(len(x))
    M = np.zeros(len(x))
    
    for a, P in zip(ain, Pin):
        V += -P/L
        M[x<=a] += -P*x[x<=a]/L
        M[x>a] += P*(1-x[x>a]/L)    
    
    plt.figure()
    plt.title('Point Moment Loads')
    plt.subplot(2,1,1)    
    plt.stem(x,V)
    plt.ylabel('V,shear')
    plt.subplot(2,1,2)
    plt.stem(x,M)
    plt.ylabel('M,moment')    

def show_ae(autoencoder):
    encoder = autoencoder.Encoder()
    decoder = autoencoder.Decoder()

    encoded_imgs = encoder.predict(X_test)
    decoded_imgs = decoder.predict(encoded_imgs)

    n = 10
    plt.figure(figsize=(20, 6))
    for i in range(n):

        ax = plt.subplot(3, n, i + 1)
        plt.imshow(X_test[i].reshape(28, 28))
        plt.gray()
        ax.get_xaxis().set_visible(False)
        ax.get_yaxis().set_visible(False)

        ax = plt.subplot(3, n, i + 1 + n)
        plt.stem(encoded_imgs[i].reshape(-1))
        plt.gray()
        ax.get_xaxis().set_visible(False)
        ax.get_yaxis().set_visible(False)

        ax = plt.subplot(3, n, i + 1 + n + n)
        plt.imshow(decoded_imgs[i].reshape(28, 28))
        plt.gray()
        ax.get_xaxis().set_visible(False)
        ax.get_yaxis().set_visible(False)

    plt.show()

    def plotVibEpoch(self, epochTimes, signal=None, points=False):

        # indStart = int(epochTimes[0] * qu.s * self.entireVibrationSignal.sampling_rate + self.recordingStartIndex)
        # indEnd = int(epochTimes[1] * qu.s * self.entireVibrationSignal.sampling_rate + self.recordingStartIndex)

        # epochTVec = self.entireVibrationSignal.t_start + np.arange(indStart, indEnd) * self.entireVibrationSignal.sampling_period

        # plt.plot(epochTVec, self.entireVibrationSignal[indStart:indEnd], 'g' + extra)

        indStart = int(epochTimes[0] * qu.s * self.vibrationSignalDown.sampling_rate)
        indEnd = int(epochTimes[1] * qu.s * self.vibrationSignalDown.sampling_rate)

        epochTVec = self.vibrationSignalDown.t_start + np.arange(indStart, indEnd) * self.vibrationSignalDown.sampling_period

        stimEnds = (np.array(self.stimEndInds)) / self.downSamplingFactor
        stimStarts = (np.array(self.stimStartInds)) / self.downSamplingFactor

        stimEndsPresent = [x * self.vibrationSignalDown.sampling_period + self.vibrationSignalDown.t_start
                           for x in stimEnds if indStart <= x <= indEnd]
        stimStartsPresent = [x * self.vibrationSignalDown.sampling_period + self.vibrationSignalDown.t_start
                             for x in stimStarts if indStart <= x <= indEnd]

        extra = ''
        if points:
            extra = '*-'

        plt.plot(epochTVec, self.vibrationSignalDown[indStart:indEnd], 'g' + extra)
        plt.stem(stimStartsPresent, np.ones(np.shape(stimStartsPresent)), 'k')
        plt.stem(stimEndsPresent, np.ones(np.shape(stimEndsPresent)), 'm')
        if not signal is None:
            plt.plot(epochTVec, signal[indStart:indEnd], 'r' + extra)

        plt.plot(epochTVec, 2 * self.vibrationSignalDownStdDev * np.ones(epochTVec.shape), 'y')

def stemplotf(v, p):
	fig2=plt.figure()
	fig2.suptitle('Problem 2b: Data Set %s' % p)
	plt.stem(v)
	plt.xlabel('k')
	plt.ylabel('1/E[alpha_k]')
	fig2.savefig('Problem2b_data%s' % p)

def plot_station_res(sta_res,ires,perr=False,scale=1):
    
    titles = ('Tangent Force','Normal Force','Tangent Moment','Normal moment','Inflow Angle','Angle of Attack','Reynolds number',
              'Local vel.','Axial induction','Radial Induction','Effective Velocity','Lift Coeff','Drag Coeff','Loss Factor',
              'CT','CQ','CP','dP','radius','azim','height','ind_vel_a','ind_vel_r')

    axis = get_station_res(sta_res,18)
    res = get_station_res(sta_res,ires)

    axis = np.array(axis) * scale

    #plt.figure()
    if perr is True:
        plt.subplot(211)
        plt.plot(axis,res,'-o')
        plt.subplot(212)
        err = get_station_res(sta_res,-1)
        plt.stem(axis, err, 'r-')
    else:
        plt.figure()
        plt.plot(axis,res,'-o')
    plt.grid()
    plt.title(titles[ires])
    plt.show()

    return 0

def stepplot(x, y, labels, plot_titles):

        """Generates Correlation Graph. 
        With the x,y coordinates, labels, and plot titles
        established, the step-plots can be generated. Output 
        is PDF file format"""

        plt.figure()      #makes new image for each plot  
        
        #plot x & y stemplot. format plot points
        plt.stem(x, y, linefmt='k--', markerfmt='ro', basefmt='k-')               
        
        #set x-axis labels and set them vertical. size 10 font.
        plt.xticks(x, labels, rotation='vertical', fontsize = 10)
              
        #set titles for graph and axes
        plt.title(plot_titles[0])
        plt.xlabel("Biomarkers")
        plt.ylabel("Correlation Values")
        
        # slightly move axis away from plot. prevents clipping the labels 
        plt.margins(0.2)
        
        # Tweak spacing to prevent clipping of tick-labels
        plt.subplots_adjust(bottom=0.15)
        
        plt.tight_layout()   #prevents labels from being clipped
                  

        with PdfPages(plot_titles[0]+'.pdf') as pdf: #creates new file for each figure
            pdf.savefig()

def show_spec_in_graph(graph, vertex, spec, pos, weight, file_name):
    dist = 1.0 - squareform(pdist(spec.T, 'cosine'))

    plt.figure()
    plt.stem(dist[vertex, :], markerfmt=' ')

    rim = graph.new_vertex_property('vector<double>')
    rim.set_2d_array(np.array([0, 0, 0, 1]))
    rim[graph.vertex(vertex)] = [0.8941176471, 0.1019607843, 0.1098039216, 1]
    rim_width = graph.new_vertex_property('float', vals=0.5)
    rim_width.a[vertex] = 2
    shape = graph.new_vertex_property('int', vals=0)
    shape[graph.vertex(vertex)] = 2
    size = graph.new_vertex_property('double', vals=10)
    size.a[vertex] = 15
    correlation = graph.new_vertex_property('double', vals=2)
    correlation.a = dist[vertex, :]
    vorder = graph.new_vertex_property('int', vals=0)
    vorder.a[vertex] = 1

    palette = sns.cubehelix_palette(256)
    cmap = colors.ListedColormap(palette)
    gt_draw.graph_draw(graph, pos=pos, vertex_color=rim, vorder=vorder,
                       vertex_pen_width=rim_width,
                       vertex_shape=shape, vertex_fill_color=correlation,
                       vcmap=cmap, vertex_size=size, edge_color=[0, 0, 0, 0.7],
                       edge_pen_width=weight, output=file_name + '.png',
                       output_size=(1200, 1200))

    plt.figure()
    utils.plot_colorbar(cmap, np.arange(0, 1.01, 0.2), file_name)

	def graph(b,a=1):
			#make a graph
			w,h = signal.freqz(b,a)
			h_dB = 20 * np.log10 (abs(h))
			plt.figure()
			#plt.subplot(311)
			plt.plot(w/max(w),h_dB)
			plt.ylim(-150, 5)
			plt.ylabel('Magnitude (db)')
			plt.xlabel(r'Normalized Frequency (x$\pi$rad/sample)')
			plt.title(r'Frequency response')
			plt.show()
			plt.figure()
			l = len(b)
			impulse = np.repeat(0.,l); impulse[0] =1.
			x = arange(0,l)
			response = signal.lfilter(b,a,impulse)
			#plt.subplot(312)
			plt.stem(x, response)
			plt.ylabel('Amplitude')
			plt.xlabel(r'n (samples)')
			plt.title(r'Impulse response')
			plt.show()
			#plt.figure()
			#plt.subplot(313)
			#step = np.cumsum(response)
			#plt.stem(x, step)
			#plt.ylabel('Amplitude')
			#plt.xlabel(r'n (samples)')
			#plt.title(r'Step response')
			#plt.subplots_adjust(hspace=0.5)
			#plt.show()
			return 1

def viz_rank(Oi,k=None):
    if k is None:
        U, s, VT = np.linalg.svd(Oi)
    else:
        U, s, VT = thin_svd_randomized(Oi,k)

    plt.figure()
    plt.stem(np.arange(s.shape[0]),s)

def figure_binary_offset(title, axes, f_parent, signature, n, i):
    f = f_parent.add_subplot(*axes, title=title)
    plt.stem(n, i, linefmt='b-', markerfmt='b.', basefmt='b|')
    plt.axis([n[0], n[-1], signature.min, signature.max],
        figure=f)
    plt.xlabel('Time')
    plt.ylabel('Magnitude')
    return f

 def simple_plot(self, arr, filename, type='normal'):
     plt.figure()
     if type == 'normal':
         plt.plot(arr)
     elif type == 'stem':
         plt.stem(arr)
     plt.savefig('../images/{}.jpg'.format(filename))
     plt.close()

 def _make_channel_plots(self):
     plt.figure()
     channel = self._system.get_block(self._system.MULTI_PATH_CHANNEL)
     plt.stem(channel._impulse_response, basefmt='')
     plt.title(self._build_title("Odpowiedź impulsowa kanału"))
     plt.grid('on')
     length = len(channel._impulse_response)
     plt.xlim((-length / 20, length))
     self._save_plt("channel_impulse")

def training(nfiltbank, orderLPC):
    nSpeaker = 8
    nCentroid = 16
    codebooks_mfcc = np.empty((nSpeaker,nfiltbank,nCentroid))
    codebooks_lpc = np.empty((nSpeaker, orderLPC, nCentroid))
    directory = os.getcwd() + '/train';
    fname = str()

    for i in range(nSpeaker):
        fname = '/s' + str(i+1) + '.wav'
        print('Now speaker ', str(i+1), 'features are being trained' )
        (fs,s) = read(directory + fname)
        mel_coeff = mfcc(s, fs, nfiltbank)
        lpc_coeff = lpc(s, fs, orderLPC)
        codebooks_mfcc[i,:,:] = lbg(mel_coeff, nCentroid)
        codebooks_lpc[i,:,:] = lbg(lpc_coeff, nCentroid)
        
        plt.figure(i)
        plt.title('Codebook for speaker ' + str(i+1) + ' with ' + str(nCentroid) +  ' centroids')
        for j in range(nCentroid):
            plt.subplot(211)
            plt.stem(codebooks_mfcc[i,:,j])
            plt.ylabel('MFCC')
            plt.subplot(212)
            markerline, stemlines, baseline = plt.stem(codebooks_lpc[i,:,j])
            plt.setp(markerline,'markerfacecolor','r')
            plt.setp(baseline,'color', 'k')
            plt.ylabel('LPC')
            plt.axis(ymin = -1, ymax = 1)
            plt.xlabel('Number of features')
    
    plt.show()
    print('Training complete')
    
    #plotting 5th and 6th dimension MFCC features on a 2D plane
    #comment lines 54 to 71 if you don't want to see codebook
    codebooks = np.empty((2, nfiltbank, nCentroid))
    mel_coeff = np.empty((2, nfiltbank, 68))
   
    for i in range(2):
        fname = '/s' + str(i+2) + '.wav'
        (fs,s) = read(directory + fname)
        mel_coeff[i,:,:] = mfcc(s, fs, nfiltbank)[:,0:68]
        codebooks[i,:,:] = lbg(mel_coeff[i,:,:], nCentroid)
        
    
    plt.figure(nSpeaker + 1)
    s1 = plt.scatter(mel_coeff[0,6,:], mel_coeff[0,4,:],s = 100,  color = 'r', marker = 'o')
    c1 = plt.scatter(codebooks[0,6,:], codebooks[0,4,:], s = 100, color = 'r', marker = '+')
    s2 = plt.scatter(mel_coeff[1,6,:], mel_coeff[1,4,:],s = 100,  color = 'b', marker = 'o')
    c2 = plt.scatter(codebooks[1,6,:], codebooks[1,4,:], s = 100, color = 'b', marker = '+')
    plt.grid()
    plt.legend((s1, s2, c1, c2), ('Sp1','Sp2','Sp1 centroids', 'Sp2 centroids'), scatterpoints = 1, loc = 'upper left')    
    plt.show()
   
    
    return (codebooks_mfcc, codebooks_lpc)

def plotHistogram(fg):
   nSubBands = fg.getNumberOfSubBands()
   histogram = []
   for i in range(0,nSubBands):
      histogram.append(fg.getHistogram(i))
   plt.stem(range(0,nSubBands), histogram)
   plt.ylabel('Number of FA in a given sub-band')
   plt.xlabel('sub-band')
   plt.draw()