def grafix1(VP,VPp,m,x,y,c):
    error = []
    for i in range(len(VP)):
        error.append(abs(VP[i]-VPp[i]))
    bins_s=60
    bins_vp = np.linspace(min(VP), max(VP), bins_s)
    scatterP= m+' \n $r=$'+str(round(np.corrcoef(VP,VPp)[0,1],2))
    label_hist_pl = '\n '+x+'\n $\overline{e} =$'+str(round(np.mean(error))) \
    +'\n $\sigma_e =$'+str(round(np.std(error)))
    #--------------------------------------------------------------------------------------------------#
    X_VP  = np.linspace(min(VP), max(VP),bins_s)
    dx_VP = np.histogram(VP ,bins=bins_vp)[1][1] - np.histogram(VP ,bins=bins_vp)[1][0]
    Y_VP  = mlab.normpdf(np.linspace(min(VP),max(VP),bins_s),np.mean(VP),np.sqrt(np.var(VP)))*len(VP)*dx_VP
    #-----------------------------------------------------------------------------------------------------#
    X_VPp  = np.linspace(min(VPp), max(VPp),bins_s)
    dx_VPp = np.histogram(VPp ,bins=bins_vp)[1][1] - np.histogram(VPp ,bins=bins_vp)[1][0]
    Y_VPp  = mlab.normpdf(np.linspace(min(VPp),max(VPp),bins_s),np.mean(VPp),np.sqrt(np.var(VPp)))*len(VPp)*dx_VPp
    #-----------------------------------------------------------------------------------------------------#
    fig = plt.figure(figsize= (12,12))

    ax1 = plt.subplot(222)
    ax1.hist(VP,bins_vp,histtype='bar',stacked=True,color='k',alpha=0.5,label='Valores $VP$')
    ax1.plot(X_VP,Y_VP,linewidth=2,color='k')
    ax1.hist(VPp , bins_vp, histtype='bar', stacked=True, color=c, alpha=0.3,label=label_hist_pl)
    ax1.plot(X_VPp,Y_VPp,linewidth = 2, color=c)
    plt.xlabel('Velocidades $(m / s)$');plt.ylabel('Distribuição');plt.grid();plt.xlim(xmax=max(VP),xmin=min(VP));
    plt.ylim(ymax=180,ymin=0);legend = ax1.legend(loc=1, shadow=True)

    ax2=plt.subplot(221);ax2.plot(VP,VP,'+k');ax2.plot(VP,VPp,'+'+c,label=scatterP);legend=ax2.legend(loc=4)
    plt.xlim(xmax=max(VP),xmin=min(VP));plt.ylim(ymax=max(VP),ymin=min(VP));
    plt.xlabel('Velocidade Original $VP$ em $m/s$')
    plt.ylabel('Velocidade Estimada '+x+' em $m/s$');plt.grid()

    plt.show()

def classify_2d(data_a, data_b, x):
    x1 = x[0]
    x2 = x[1]

    probability_a = data_a.shape[1] / (data_a.shape[1] + data_b.shape[1])
    probability_b = data_b.shape[1] / (data_a.shape[1] + data_b.shape[1])

    mean_x1_a = np.mean(data_a[0,:])
    mean_x2_a = np.mean(data_a[1,:])

    mean_x1_b = np.mean(data_b[0,:])
    mean_x2_b = np.mean(data_b[1,:])

    variance_x1_a = np.var(data_a[0,:])
    variance_x2_a = np.var(data_a[1,:])

    variance_x1_b = np.var(data_b[0,:])
    variance_x2_b = np.var(data_b[1,:])

    pd_x1_given_a = mlab.normpdf(x1, mean_x1_a, variance_x1_a)
    pd_x2_given_a = mlab.normpdf(x2, mean_x2_a, variance_x2_a)
    pd_x1_given_b = mlab.normpdf(x1, mean_x1_b, variance_x1_b)
    pd_x2_given_b = mlab.normpdf(x2, mean_x2_b, variance_x2_b)

    posterior_numerator_a = probability_a * pd_x1_given_a * pd_x2_given_a
    posterior_numerator_b = probability_b * pd_x1_given_b * pd_x2_given_b

    posterior_numerators = { 'A': posterior_numerator_a, 'B': posterior_numerator_b }

    return max(posterior_numerators.iterkeys(), key=(lambda k: posterior_numerators[k]))

def main():
    x1 = 0.0
    N = 5000
    x1Values = np.zeros(N)
    x2Values = np.zeros(N)

    for i in range(N):
        x2 = updateX2(x1)
        x2Values[i] = x2
        x1 = updateX1(x2)
        x1Values[i] = x1

    #plot
    plt.hist(x1Values, bins=20, alpha=0.6, label='calc. marginal', normed=True)
    x = np.linspace(-3,5)
    plt.plot(x,mlab.normpdf(x,1.0,1.0), lw=3)
    plt.xlabel("Value")
    plt.ylabel("Frequency")
    plt.legend(loc='upper right')
    plt.savefig('px1')
    plt.clf()

    plt.hist(x2Values, bins=20, alpha=0.6, label='calc. marginal', normed=True)
    x = np.linspace(-3,5)
    plt.plot(x,mlab.normpdf(x,1.0,1.0), lw=3)
    plt.xlabel("Value")
    plt.ylabel("Frequency")
    plt.legend(loc='upper right')
    plt.savefig('px2')

def plot_score_distributions(threshold, neg_devel, pos_devel, neg_test, pos_test, filename='score_dist.png'):

    plt.clf()
    plt.figure(1)
    plt.subplot(211)
    plt.title("Score distributions (Deve set)")
    n, bins, patches = plt.hist(neg_devel, bins=25, normed=1, histtype='bar', label='Negative class')
    na, bins_a, patches_a = plt.hist(pos_devel, bins=25, normed=1, histtype='bar', label='Positive class')

    # add a line showing the expected distribution
    y = mlab.normpdf(bins, np.mean(neg_devel), np.std(neg_devel))
    plt.plot(bins, y, 'k--', linewidth=1.5)
    y = mlab.normpdf(bins_a, np.mean(pos_devel), np.std(pos_devel))
    plt.plot(bins_a, y, 'k--', linewidth=1.5)
    plt.axvline(x=threshold, linewidth=2, color='blue')
    plt.legend()

    plt.subplot(212)
    plt.title("Score distributions (Test set)")
    n, bins, patches = plt.hist(neg_test, bins=25, normed=1, facecolor='green', alpha=0.5, histtype='bar',
                                label='Negative class')
    na, bins_a, patches_a = plt.hist(pos_test, bins=25, normed=1, facecolor='red', alpha=0.5, histtype='bar',
                                     label='Positive class')

    # add a line showing the expected distribution
    y = mlab.normpdf(bins, np.mean(neg_test), np.std(neg_test))
    plt.plot(bins, y, 'k--', linewidth=1.5)
    y = mlab.normpdf(bins_a, np.mean(pos_test), np.std(pos_test))
    plt.plot(bins_a, y, 'k--', linewidth=1.5)
    plt.axvline(x=threshold, linewidth=2, color='blue')
    plt.legend()

    current_dir = os.getcwd()
    output = '{0}/{1}.png'.format(current_dir, filename)
    plt.savefig(output)

def naive_bayes(w1train,w2train,test):
	# prior
	n = w1train.shape[0]+w2train.shape[0]
	w_1 = w1train.shape[0] / float(n)
	w_2 = w2train.shape[0] / float(n)
	print 'prior w1:', w_1
	print 'prior w2:', w_2
	# likelihood
	mu_1, s_1 = gauss_mle_1d(w1train)
	mu_2, s_2 = gauss_mle_1d(w2train)
	post_1 = gaussian(test,mu_1,s_1)
	post_2 = gaussian(test,mu_2,s_2)
	print 'p(w1|x)=',post_1
	print 'p(w2|x)=',post_2
	p_1 = post_1*w_1
	p_2 = post_2*w_2
	print 'class 1',p_1
	print 'class 2',p_2
	print 'bla_1', p_1 / (p_1+p_2)
	print 'bla_2', p_2 / (p_1+p_2)
	x1 = np.linspace(-3,8,100)
	plt.title('Epic Info')
	plt.ylabel('Y axis')
	plt.xlabel('X axis')
	plt.plot(x1,mlab.normpdf(x1,mu_1,s_1),label='estimate class1')
	plt.plot(x1,mlab.normpdf(x1,mu_2,s_2),label='estimate class2')
	plt.legend()
	plt.text(-2,0.7,'class 1:%s\nclass2: %s'%(p_1,p_2))
	plt.plot(test,0,'o',label='test point')
	plt.show()

    def test_kernel_smoothing(self):
        # Qualitatively view kernel smoothed noisy Gaussian landscape

        # Make Normal distribution, and secondary smaller normal.
        realNorm = np.array([mlab.normpdf(i,40,10) for i in range(100)])
        realNorm = realNorm / np.max(realNorm)
        gaussBlip = np.array([mlab.normpdf(i,80,3) for i in range(100)])
        gaussBlip = gaussBlip / np.max(gaussBlip)
        signal = realNorm + (gaussBlip * 0.5)

        # Add noise.
        noise = np.random.random(100)        
        signal = (signal * noise) + (0.2 * noise)
        signal = signal / np.sum(signal)
        signal = np.concatenate((signal, signal))
        
        # Apply kernel smoothing.
        smoothGauss = analysis.force._kernel_smoothing(signal, 0.25)
        smootherGauss = analysis.force._kernel_smoothing(signal, 0.75)
        smoothestGauss = analysis.force._kernel_smoothing(signal, 0.98)

        # Compare kernel smoothed plot to noisy data plot.
        plot.plot(signal, 'k')
        plot.hold(True)
        plot.plot(smoothGauss, 'm--')
        plot.plot(smootherGauss, 'c--')
        plot.plot(smoothestGauss, 'r--')
        plot.hold(False)
        plot.show()

def get_prob_for_distributions(p):
    """
    Based on the integral of the three normal distributions,
    the likelihood from which of the three distributions a distance is to be drawn 
    is calculated here.
    Returns the three probabilities for the three distributions.
    """
    w1 = p[0]
    mu1 = p[1]
    sigma1 = p[2]
    w2 = p[3]
    mu2 = p[4]
    sigma2 = p[5]
    w3 = p[6]
    mu3 = p[7]
    sigma3 = p[8]
    dist_range = (0, 4.330310991999920844e+01)
    x = np.linspace(dist_range[0], dist_range[1], 1000)
    A1 = np.array(w1 * mlab.normpdf(x, mu1, sigma1)).sum()
    A2 = np.array(w2 * mlab.normpdf(x, mu2, sigma2)).sum()
    A3 = np.array(w3 * mlab.normpdf(x, mu3, sigma3)).sum()
    p1 = A1 / (A1 + A2 + A3)
    p2 = A2 / (A1 + A2 + A3)
    p3 = A3 / (A1 + A2 + A3)
    return p1, p2, p3

def kldistancecluster(planets):
    nlist = nall(knownplanets, 'earth')[1]
    ntrue = nall(knownplanets, 'earth')[0]
    difference = variance(ntrue,nlist)
    uniformdist = np.asarray(np.random.uniform(0.0,0.5,len(difference)))
    difference = np.asarray(difference)

    plt.hist(nlist, bins = 25, color = 'blue', alpha = 0.7, normed = True)
    plt.hist(ntrue, bins = 25, color = 'green', alpha = 0.5, normed = True)
    
    # Find best fit
    x = np.linspace(0.0, 8, 25)
    best_fit_uniform = mlab.normpdf(x, np.mean(nlist), np.std(nlist))
    best_fit_dif = mlab.normpdf(x, np.mean(ntrue), np.std(ntrue))
    plt.plot(x, best_fit_uniform, label = 'unif')
    plt.plot(x, best_fit_dif, label = 'dif')
    plt.xlabel('Distribution Value')
    plt.ylabel('Frequency')
    
    blue = mpatches.Patch(color='blue', label = 'Normed PDF for Integer Distribution')
    green = mpatches.Patch(color = 'green', label = 'Normed PDF for Calculate Rank')
    plt.legend(handles = [blue, green])
   
    plt.text(3.4, 1.0, 'KL Divergence: \n  ( rank, integer distribution) = 0.0112', style='italic', bbox={'facecolor':'red', 'alpha':0.5, 'pad':10})
    plt.show()
    #kldiv = stats.entropy(difference, qk=uniformdist, base=None)
    kldiv = stats.entropy(nlist, qk=ntrue, base=None)
    return(kldiv)

def plot_bourgdata(N1,N2):
	A=TRICLAIRModele()
	Tb15 = A.get_data_triathlon(link='/triathlon-bourg-resultats-1996.htm',year=2015)
	Tb14 = A.get_data_triathlon(link='/triathlon-bourg-resultats-1715.htm',year=2014)	
	S15_ = map(lambda x: x.total_seconds()/60,Tb15['Scratch'].dropna())
	S14_ = map(lambda x: x.total_seconds()/60,Tb14['Scratch'].dropna())
	
	S15 = S15_[N1:N2]
	S14 = S14_[N1:N2]

	(mu14, sigma14) = norm.fit(S14)
	(mu15, sigma15) = norm.fit(S15)

	N_BINS = 50
	fig = plt.figure()
	ax = fig.add_subplot(1, 1, 1)
	n, bins, patches = ax.hist(S14, N_BINS,normed=1, facecolor='red', alpha=0.5,label=r'$\mathrm{2014:}\ \mu=%.3f,\ \sigma=%.3f$' %(mu14, sigma14))
	y = mlab.normpdf( bins, mu14, sigma14)
	l = ax.plot(bins, y, 'r-', linewidth=4)
	n, bins, patches = ax.hist(S15, N_BINS, normed=1, facecolor='green', alpha=0.5,label=r'$\mathrm{ 2015:}\ \mu=%.3f,\ \sigma=%.3f$' %(mu15, sigma15))
	y = mlab.normpdf( bins, mu15, sigma15)
	l = ax.plot(bins, y, 'g-', linewidth=4)

	fig.tight_layout()
	ax.set_xlabel('Scratch Time (minutes)')
	ax.set_ylabel('Number of athletes per scratch time (normalized)')
	ax.legend(loc='best', fancybox=True, framealpha=0.5)
	ax.set_title(r'$\mathrm{Athletes\ from\ rank\ } %d \mathrm{\ to\ } %d$' %(N1, N2))

	plt.show()
	

def avg_score_distribution(n, m, data_points=10e4, bins=100, visualize=False):
    '''
    Returns estimated mean and standard deviation of score distribution
    for randomized amino acid recognition result
     
    n := sum of all fragment lengths
    m := length of sequence
    '''
    
    assert n <= m
    
    scores = []
    for i in range(int(data_points)):
    
        p = random(n)
        avg = 1 - ((m - n + p.sum()) / m)
        scores.append(avg)
        
    data = array(scores)
    mu = mean(data) ## mean value
    sigma = std(data) ## standard deviation
    
    if visualize:
        n, bins, patches = plt.hist(data, bins, normed=1, alpha=.3)
        y = mlab.normpdf(bins, mu, sigma)
        plt.plot(bins, y, 'r-', linewidth=1)
        plt.vlines(mu, 0, mlab.normpdf([mu], mu, sigma), colors='r')
        plt.show()
    
    return mu, sigma

def compute_costed_threshold(weight, thresholds, meanS1, sdS1, meanS2, sdS2):
    """Compute the costed threshold of two spike responses to
two independent stimuli.

    Args:
        weight: costed threshold multiplier
        thresholds: values to test for suitability
        meanS1: mean of firing rate distribution triggered by stimulus 1
        sdS1: standard deviation of firing rate distribution triggered by stimulus 1
        meanS2: mean of firing rate distribution triggered by stimulus 2
        sdS1: standard deviation of firing rate distribution triggered by stimulus 2

    Returns:
        neuronal firing tate as which to set optimum costed threshold"""

    opt_thresh = 0.0

    for threshold in thresholds:
        Ps1 = mlab.normpdf(threshold, meanS1, sdS1)
        Ps2 = mlab.normpdf(threshold, meanS2, sdS2)

        ratio_raw = Ps2 / Ps1  # make likelihood twice for Ps2
        ratio = round(ratio_raw, 1)

        print "Ps1 = %s, Ps2 = %s. Threshold = %s. Ratio raw = %s Weight = %s\n" % (Ps1, Ps2, threshold, ratio_raw, ratio)
        if ratio == weight:
            opt_thresh = threshold

    return opt_thresh

    def extract_coarse_coding_features_absolute(self, phone_duration):
        dur = int(phone_duration)

        cc_feat_matrix = numpy.zeros((dur, 3))

        npoints1 = (dur*2)*10+1
        npoints2 = (dur-1)*10+1
        npoints3 = (2*dur-1)*10+1

        x1 = numpy.linspace(-dur, dur, npoints1)
        x2 = numpy.linspace(1, dur, npoints2)
        x3 = numpy.linspace(1, 2*dur-1, npoints3)

        mu1 = 0
        mu2 = (1+dur)/2
        mu3 = dur
        variance = 1
        sigma = variance*((dur/10)+2)
        sigma1 = sigma
        sigma2 = sigma-1
        sigma3 = sigma

        y1 = mlab.normpdf(x1, mu1, sigma1)
        y2 = mlab.normpdf(x2, mu2, sigma2)
        y3 = mlab.normpdf(x3, mu3, sigma3)

        for i in range(dur):
            cc_feat_matrix[i,0] = y1[(dur+1+i)*10]
            cc_feat_matrix[i,1] = y2[i*10]
            cc_feat_matrix[i,2] = y3[i*10]

        for i in range(3):
            cc_feat_matrix[:,i] = cc_feat_matrix[:,i]/max(cc_feat_matrix[:,i])

        return cc_feat_matrix

def plotting(ls1, ls2, head):
    # ls1 normal, ls2 satire
    mu1 = np.mean(ls1)
    mu2 = np.mean(ls2)
    sigma1 = np.std(ls1) # standard deviation of distribution
    sigma2 = np.std(ls2)
    x = ls1
    y = ls2
    plt.figure(1)
    num_bins = 100
    # the histogram of the data
    n1, bins1, patches1 = plt.hist(x, num_bins, normed=1, facecolor='green', alpha=0.5, label='normal')
    plt.legend(loc=2)
    # add a 'best fit' line
    x1 = mlab.normpdf(bins1, mu1, sigma1)

    n2, bins2, patches2 = plt.hist(y, num_bins, normed=1, facecolor='red', alpha=0.5, label = 'satire')
    plt.legend(loc=2)
    # add a 'best fit' line
    y2 = mlab.normpdf(bins2, mu2, sigma2)

    plt.plot(bins1, x1, 'b--')
    plt.plot(bins2, y2, 'b--')
    plt.xlabel('Variance')
    plt.ylabel('Density')
    plt.title('Distribution of docs')
    

    # Tweak spacing to prevent clipping of ylabel
    plt.subplots_adjust(left=0.15)
    if head == True: filename = 'dist_head.png'
    elif head == False: filename = 'dist.png'
    plt.savefig(filename, format='png')

def gaussian_1d(data, Pi, means, sds, N, K):
	x = np.linspace(min(data), max(data), N);
	gmm = 0*mlab.normpdf(x,0,1);
	for i in range(len(Pi)):
		gmm += Pi[i]*mlab.normpdf(x,means[i],sds[i]);
	
	plt.plot(x, gmm);

def peval_binormal(x, p):
    # p[0] = w1
    # p[1] = mu1
    # p[2] = sigma1
    # p[3] = w2
    # p[4] = mu2
    # p[5] = sigma2
    return (p[0] * mlab.normpdf(x, p[1], p[2]) + p[3] * mlab.normpdf(x, p[4], p[5]))

 def max_score_func(m):
     mu = 40
     sigma = 8
     if m <= mu:
         return 1.
     else:
         top_val = normpdf(mu, mu, sigma)
         return normpdf(m, mu, sigma) / top_val

def kalman_plot(prediction, measurement, correction):
    """Helper to draw all curves in each filter step."""
    plot([normpdf(x, prediction.mu, sqrt(prediction.sigma2))
          for x in range(*arena)], color = 'b', linewidth=2)
    plot([normpdf(x, measurement.mu, sqrt(measurement.sigma2))
          for x in range(*arena)], color = 'g', linewidth=2)
    plot([normpdf(x, correction.mu, sqrt(correction.sigma2))
          for x in range(*arena)], color = 'r', linewidth=2)

def plot_density(mean, class_prior):
    x = np.linspace(-5, 7, 100)
    pyplot.xticks([1*k for k in range(-5, 8)])
    pyplot.plot(x, (mlab.normpdf(x, mean[0], 1)*class_prior[0] + mlab.normpdf(x, mean[1], 1)*class_prior[1]))
    pyplot.title("Density plot for the data")
    pyplot.xlabel("Data")
    pyplot.ylabel("Density")
    pyplot.savefig('./density_plot.png')
    pyplot.close()

  def plot_PDFs(self):
    """
    Plot probability density functions
    """
    import matplotlib.mlab as mlab
    NB_class = len(self.opdict['types'])
    feat_1 = self.x.columns[0]
    feat_2 = self.x.columns[1]

    binwidth = .05
    lim_sup_1 = (int(np.max(self.x[feat_1])/binwidth)+2)*binwidth
    lim_inf_1 = (int(np.min(self.x[feat_1])/binwidth)-2)*binwidth
    bins_1 = np.arange(lim_inf_1, lim_sup_1 + binwidth, binwidth)
    lim_sup_2 = (int(np.max(self.x[feat_2])/binwidth)+2)*binwidth
    lim_inf_2 = (int(np.min(self.x[feat_2])/binwidth)-2)*binwidth
    bins_2 = np.arange(lim_inf_2, lim_sup_2 + binwidth, binwidth)

    x_hist, y_hist = [],[]
    g_x, g_y = {}, {}
    for i in range(NB_class):
      index = self.y[self.y.NumType.values==i].index
      x1 = self.x.reindex(columns=[feat_1],index=index).values
      x2 = self.x.reindex(columns=[feat_2],index=index).values
      g_x[i] = mlab.normpdf(bins_1, np.mean(x1), np.std(x1))
      g_y[i] = mlab.normpdf(bins_2, np.mean(x2), np.std(x2))
      x_hist.append(x1)
      y_hist.append(x2)

    if NB_class > 2:
      colors_g = ('y','orange','r')
      colors_h = ('k','gray','w')
    elif NB_class == 2:
      colors_g = ('y','r')
      colors_h = ('k','w')

    fig = plt.figure()
    fig.set_facecolor('white')
    plt.hist(x_hist,bins=bins_1,color=colors_h,normed=1,histtype='stepfilled',alpha=.5)
    for key in sorted(g_x):
      plt.plot(bins_1,g_x[key],color=colors_g[key],lw=2.,label='Class %s'%self.opdict['types'][key])
    plt.xlabel(feat_1)
    plt.legend(loc=2)
    plt.savefig('%s/histo_%s.png'%(self.opdict['fig_path'],feat_1))

    fig = plt.figure()
    fig.set_facecolor('white')
    plt.hist(y_hist,bins=bins_2,color=colors_h,normed=1,histtype='stepfilled',alpha=.5)
    for key in sorted(g_y):
      plt.plot(bins_2,g_y[key],color=colors_g[key],lw=2.,label='Class %s'%self.opdict['types'][key])
    plt.plot([.5,.5],[0,2],'g--',lw=2.)
    plt.figtext(.52,.7,'?',color='g',size=20)
    plt.xlabel(feat_2)
    plt.legend(loc=2)
    plt.savefig('%s/histo_%s.png'%(self.opdict['fig_path'],feat_2))
    plt.show()

def plot(muVar, sigmaVar, muSDTrue, sigmaSDTrue):
    xMin = min(muVar,muSDTrue)-2*max(sigmaVar,np.sqrt(varSDTrue))
    xMax = max(muVar,muSDTrue)+2*max(sigmaVar,np.sqrt(varSDTrue))
    x = np.linspace(xMin,xMax, 1000)
    #plt.title('m= %f, sigma_e=%f, weight=2, no bias' % (m, sigma_e)) 
    plt.plot(x,mlab.normpdf(x,muVar,sigmaVar),label='variational mu=%f, sd=%f' % (muVar, sigmaVar),color='blue')
    plt.plot(x,mlab.normpdf(x,muSDTrue,sigmaSDTrue), ls='--',color='red')
    #, label='true mu=%f, sd=%f' % (muSDTrue, sigmaSDTrue))
    #plt.legend(('variational mu=%f, sd=%f' % (muVar, sigmaVar), 'true posterior mu=%f, sd=%f' % (muSDTrue, sigmaSDTrue)))
    #plt.legend()
    plt.show()