def histogram(A, B, nameA, nameB):
	plt.hist(A, bins=255, alpha=0.5, color='b', label = nameA)
	plt.hist(B, bins=255, alpha=0.5, color='r', label = nameB)
	plt.xlabel('Intensity')
	plt.ylabel('Number of occurrencies')
	plt.legend()
	plt.show()

def hist(fname, data, bins, xlabel, ylabel, title, facecolor='green', alpha=0.5, transparent=True, **kwargs):
    plt.clf()
    plt.xlabel(xlabel)
    plt.ylabel(ylabel)
    plt.title(title)
    plt.hist(x=data, bins=bins, facecolor=facecolor, alpha=alpha, **kwargs)
    plt.savefig(fname, transparent=transparent)

def predicted_probabilities(y_true, y_pred, n_groups=30):
    """Plots the distribution of predicted probabilities.

    Parameters
    ----------
    y_true : array_like
        Observed labels, either 0 or 1.
    y_pred : array_like
        Predicted probabilities, floats on [0, 1].
    n_groups : int, optional
        The number of groups to create. The default value is 30.

    Notes
    -----
    .. plot:: pyplots/predicted_probabilities.py
    """
    plt.hist(y_pred, n_groups)
    plt.xlim([0, 1])
    plt.xlabel('Predicted Probability')
    plt.ylabel('Count')

    title = 'Distribution of Predicted Probabilities (n = {})'
    plt.title(title.format(len(y_pred)))

    plt.tight_layout()

def plot_net_distribution(net_mat, n_bins):
    """Plot the network distribution.

    Parameters
    ----------
    net_mat: np.ndarray
        the net represented in a matrix way.
    n_bins: int
        the number of intervals we want to use to plot the distribution.

    Returns
    -------
    fig: matplotlib.pyplot.figure
        the figure of the distribution required of the relations between
        elements defined by the `net_mat`.

    """
    net_mat = net_mat.reshape(-1)

    fig = plt.figure()
    plt.hist(net_mat, n_bins)

    l1 = plt.axvline(net_mat.mean(), linewidth=2, color='k', label='Mean',
                     linestyle='--')
    plt.legend([l1], ['Mean'])

    return fig

def plot_ekf_vs_mc():

    def fx(x):
        return x**3

    def dfx(x):
        return 3*x**2

    mean = 1
    var = .1
    std = math.sqrt(var)

    data = normal(loc=mean, scale=std, size=50000)
    d_t = fx(data)

    mean_ekf = fx(mean)

    slope = dfx(mean)
    std_ekf = abs(slope*std)


    norm = scipy.stats.norm(mean_ekf, std_ekf)
    xs = np.linspace(-3, 5, 200)
    plt.plot(xs, norm.pdf(xs), lw=2, ls='--', color='b')
    plt.hist(d_t, bins=200, normed=True, histtype='step', lw=2, color='g')

    actual_mean = d_t.mean()
    plt.axvline(actual_mean, lw=2, color='g', label='Monte Carlo')
    plt.axvline(mean_ekf, lw=2, ls='--', color='b', label='EKF')
    plt.legend()
    plt.show()

    print('actual mean={:.2f}, std={:.2f}'.format(d_t.mean(), d_t.std()))
    print('EKF    mean={:.2f}, std={:.2f}'.format(mean_ekf, std_ekf))

def plot_scatter_with_histograms(xvals, yvals, colour='k', oneToOneLine=True, xlabel=None, ylabel=None, title=None):
    gs = gridspec.GridSpec(5, 5)
    xmin = np.floor(min(xvals))
    xmax = np.ceil(max(xvals))
    ymin = np.floor(min(yvals))
    ymax = np.ceil(max(yvals))
    plt.subplot(gs[1:, 0:4])
    plt.plot(xvals, yvals, 'o', color=colour)
    if xlabel is not None:
        plt.xlabel(xlabel)
    if ylabel is not None:
        plt.ylabel(ylabel)
    if oneToOneLine:
        oneToOneMax = max([max(xvals),max(yvals)])
        plt.plot([0,oneToOneMax],[0,oneToOneMax],'b--')
    plt.xlim(xmin,xmax)
    plt.ylim(ymin,ymax)
    plt.subplot(gs[0, 0:4])
    plt.hist(xvals, np.linspace(xmin,xmax,50))
    plt.axis('off')
    plt.subplot(gs[1:,4])
    plt.hist(yvals, np.linspace(ymin,ymax,50), orientation='horizontal')
    plt.axis('off')
    if title is not None:
        plt.suptitle(title)

def main():
    train = pd.DataFrame.from_csv('train.csv')
    places_index = train['place_id'].values

    places_loc_sqr_wei = []
    for i, place_id in enumerate(train['place_id'].unique()):
        if not i % 100:
            print(i)
        place_df = train.iloc[places_index == place_id]
        place_weights_acc_sqred = 1 / (place_df['accuracy'].values ** 2)

        places_loc_sqr_wei.append([place_id,
                                   np.average(place_df['x'].values, weights=place_weights_acc_sqred),
                                   np.std(place_df['x'].values),
                                   np.average(place_df['y'].values, weights=place_weights_acc_sqred),
                                   np.std(place_df['y'].values),
                                   np.average(np.log(place_df['accuracy'].values)),
                                   np.std(np.log(place_df['accuracy'].values)),
                                   place_df.shape[0]])

        # print(places_loc_sqr_wei[-1])
        # plt.hist2d(place_df['x'].values, place_df['y'].values, bins=100)
        # plt.show()
        plt.hist(np.log(place_df['accuracy'].values), bins=20)
        plt.show()
    places_loc_sqr_wei = np.array(places_loc_sqr_wei)
    column_names = ['x_mean', 'x_sd', 'y_mean', 'y_sd', 'accuracy_mean', 'accuracy_sd', 'n_persons']
    places_loc_sqr_wei = pd.DataFrame(data=places_loc_sqr_wei[:, 1:], index=places_loc_sqr_wei[:, 0],
                                      columns=column_names)

    now = str(datetime.datetime.now().strftime("%Y-%m-%d-%H-%M"))
    places_loc_sqr_wei.to_csv('places_loc_sqr_weights_%s.csv' % now)

def createHistogram(df, pic, bins=45, rates=False):
    data=mergeMatrix(df, pic)
    matrix=sortMatrix(df, pic)


    density = gaussian_kde(data)
    xs = np.linspace(min(data), max(data), max(data))
    density.covariance_factor = lambda : .25
    density._compute_covariance()
    #xs = np.linspace(min(data), max(data), 1000)

    fig,ax1 = plt.subplots()
    #plt.xlim([0, 4000])
    plt.hist(data, bins=bins, range=[-500, 4000], histtype='stepfilled', color='grey', alpha=0.5)
    lims = plt.ylim()
    height=lims[1]-2
    for i in range(0,len(matrix)):
        currentRow = matrix[i][np.nonzero(matrix[i])]
        plt.plot(currentRow, np.ones(len(currentRow))*height, '|', color='black')
        height -= 2

    plt.axvline(x=0, color='red', linestyle='dashed')
    #plt.axvline(x=1000, color='black', linestyle='dashed')
    #plt.axvline(x=2000, color='black', linestyle='dashed')
    #plt.axvline(x=3000, color='black', linestyle='dashed')

    if rates:
        rates = get_rate(df, pic)
        ax1.text(-250, 4, str(rates[0]), size=15, ha='center', va='center', color='green')
        ax1.text(500, 4, str(rates[1]), size=15, ha='center', va='center', color='green')
        ax1.text(1500, 4, str(rates[2]), size=15, ha='center', va='center', color='green')
        ax1.text(2500, 4, str(rates[3]), size=15, ha='center', va='center', color='green')
        ax1.text(3500, 4, str(rates[4])+ r' $\frac{\mathsf{Spikes}}{\mathsf{s}}$', size=15, ha='center', va='center', color='green')
    plt.ylim([0,lims[1]+5])
    plt.xlim([0, 4000])
    plt.title('Histogram for ' + str(pic))
    ax1.set_xticklabels([-500, 'Start\nStimulus', 500, 1000, 1500, 2000, 2500, 3000, 3500, 4000])
    plt.xlabel('Time (ms)')
    plt.ylabel('Counts (Spikes)')


    print lims
    arr_hand = getPic(pic)
    imagebox = OffsetImage(arr_hand, zoom=.3)
    xy = [3200, lims[1]+5]               # coordinates to position this image

    ab = AnnotationBbox(imagebox, xy, xybox=(30., -30.), xycoords='data',boxcoords="offset points")
    ax1.add_artist(ab)

    ax2 = ax1.twinx() #Necessary for multiple y-axes

    #Use ax2.plot to draw the hypnogram.  Be sure your x values are in seconds
    ax2.plot(xs, density(xs) , 'g', drawstyle='steps')
    plt.ylim([0,0.001])
    plt.yticks([0.0001,0.0002, 0.0003, 0.0004, 0.0005, 0.0006, 0.0007, 0.0008, 0.0009])
    ax2.set_yticklabels([1,2,3,4, 5, 6, 7, 8, 9])
    plt.ylabel(r'Density ($\cdot \mathsf{10^{-4}}$)', color='green')
    plt.gcf().subplots_adjust(right=0.89)
    plt.gcf().subplots_adjust(bottom=0.2)
    plt.savefig(pic, dpi=150)

def plotter(fromdat,filename):
    
    plt.figure() 
    bins = fromdat.bins
    plt.hist(fromdat.all_val, bins=bins, color=(0, 0, 0, 1 ),
                 histtype='step',label = 'All Hits' )
    plt.ylabel('Counts' )
    plt.xlabel('Energy kev' )
    plt.title('All Detectors Spectrum\n'+ filename )
    plt.legend(loc='upper right' ) 
    plt.show() 

    plt.figure() 
    his_det1 = plt.hist(fromdat.det1_val, bins=bins, color=(0, 0, 0, 0.7),
                 histtype='step', label = fromdat.detector1 )
   
    his_det2 = plt.hist(fromdat.det2_val, bins=bins, color=(0, 1, 0, 0.7 ),
                 histtype='step', label = fromdat.detector2 )
    plt.ylabel('Counts' )
    plt.xlabel('Energy kev' )
    plt.title('Overlay Plot Both Spectrum \n ' + filename)
    plt.legend(loc='upper right' ) 
    plt.show()

    his_det3 = plt.hist(fromdat.det3_val, bins=bins, color=(0, 0, 0, 0.5 ),
             histtype='step',label = fromdat.detector3 )
    plt.ylabel('Counts' )
    plt.xlabel('Energy kev' )
    plt.title( fromdat.detector3)
    plt.legend(loc='upper right' ) 
    plt.show() 

    def show(self):
        figure = plt.figure(self.figure_num)
        num_histograms = len(self.histograms)
        num_subplots = len(self.subplots)
        y_dim = 4.0
        x_dim = math.ceil((num_subplots + num_histograms)/y_dim)

        for i in range(len(self.subplots)):
            title, img = self.subplots[i]

            print "plotting: " + str(title)
            print img.shape

            ax = plt.subplot(x_dim, y_dim, i + 1)
            format_subplot(ax, img)
            plt.title(title)
            plt.imshow(img)

        for i in range(len(self.histograms)):
            title, img = self.histograms[i]

            print "plotting: " + str(title)
            print img.shape

            plt.subplot(x_dim,y_dim, num_subplots + i + 1)
            plt.title(title)
            plt.hist(img, bins=10, alpha=0.5)

def CNS(directory):
    print directory
    MASegDict = defaultdict(list)
    seqCount = Counter()
    numFeatures = defaultdict(list)
    speciesDistributionMaster = defaultdict(list)
    for species in [file for file in os.listdir(directory) if file.endswith('.bed')]:
        try:
            print directory+species
            seqCount[species] = 0
            speciesDistribution = Counter()
            with open(directory+species,'r') as f:
                lines = f.readlines()
                numFeatures[species] = [len(lines)]
                if species.endswith('ConservedElements.bed'):
                    for line in lines:
                        if line:
                            lineList = line.split('\t')
                            lineList2 = lineList[-1].split(';')
                            lineList3 = lineList2[1].split(',')
                            tempDict = {word.split(':')[0]:int(word.split(':')[1] != '0') for word in lineList3}
                            MASegDict[lineList2[2].replace('SegmentID=','')] = sum(tempDict.values())
                            seqCount[species] += int(lineList[2])-int(lineList[1])
                            for species2 in tempDict.keys():
                                if species2 not in speciesDistribution.keys():
                                    speciesDistribution[species2] = 0
                                else:
                                    speciesDistribution[species2] += tempDict[species2]
                else:
                    for line in lines:
                        if line:
                            lineList = line.split('\t')
                            lineList2 = lineList[-1].split(';')
                            lineList3 = lineList2[1].split(',')
                            tempDict = {word.split(':')[0]:int(word.split(':')[1] != '0') for word in lineList3}
                            seqCount[species] += int(lineList[2])-int(lineList[1])
                            for species2 in tempDict.keys():
                                if species2 not in speciesDistribution.keys():
                                    speciesDistribution[species2] = 0
                                else:
                                    speciesDistribution[species2] += tempDict[species2]
                speciesDistributionMaster[species] = speciesDistribution
                #print speciesDistributionMaster
                #print numFeatures
                #print ','.join('%s:%d'%(key,speciesDistributionMaster[species][key]) for key in speciesDistributionMaster[species].keys())
        except:
            print 'Error with ' + species
    with open(directory+'CNSStatistics.txt','w') as f:
        for species in sorted(numFeatures.keys()):
            if species:
                try:
                    f.write(species+'\nTotalSequenceAmount=%dbps\nNumberOfElements=%d\n%s\n\n'%(seqCount[species],numFeatures[species][0],'SpeciesDistribution='+','.join('%s:%d'%(key,speciesDistributionMaster[species][key]) for key in speciesDistributionMaster[species].keys())))#FIXME Add species number and graph
                except:
                    print 'Error writing ' + species
    plt.figure()
    plt.hist(MASegDict.values(),bins=np.arange(0,int(np.max(MASegDict.values()))) + 0.5)
    plt.title('Distribution of Number of Species for Conserved Segments')
    plt.ylabel('Count')
    plt.xlabel('Number of species in Conserved Segment')
    plt.savefig(directory+'SpeciesNumberDistribution.png')

def plotHist(data, bins=None, figsize=(7,7), title="", **kwargs):
    if (bins==None):
        bins=len(data)
    plt.figure(figsize=figsize);
    plt.hist(data,bins=bins, **kwargs)
    plt.title(title)
    plt.show()

def create_random_sample_from_beta(success, total, sample_size=10000, plot=False):
    """ Create random sample from the Beta distribution """

    failures = total - success
    data = stats.beta.rvs(success, failures, size=sample_size)
    if plot: hist(data, 100); show()
    return data

    def fluence_dist(self):
        """ Plots the fluence distribution and gives the mean and median fluence
        values of the sample """
        fluences = []
        for i in range(0,len(self.fluences),1):
            try:
                fluences.append(float(self.fluences[i]))

            except ValueError:
                continue

        fluences = np.array(fluences)
        mean_fluence = np.mean(fluences)
        median_fluence = np.median(fluences)
        print('Mean Fluence =',mean_fluence,'(15-150 keV) [10^-7 erg cm^-2]')
        print('Median Fluence =',median_fluence,'(15-150 keV) [10^-7 erg cm^-2]')

        plt.figure()
        plt.xlabel('Fluence (15-150 keV) [$10^{-7}$ erg cm$^{-2}$]')
        plt.ylabel('Number of GRBs')
        plt.xscale('log')
        minimum, maximum = min(fluences), max(fluences)
        plt.axvline(mean_fluence,color='red',linestyle='-')
        plt.axvline(median_fluence,color='blue',linestyle='-')
        plt.hist(fluences,bins= 10**np.linspace(np.log10(minimum),np.log10(maximum),20),color='grey',alpha=0.5)
        plt.show()

def test_power():
    a = 5.  # shape
    samples = 10000
    s1 = np.random.power(a, samples)
    s2 = common.rand_pow_array(a, samples)

    plt.figure('power test')
    count1, bins1, ignored1 = plt.hist(s1,
                                       bins=30,
                                       label='numpy',
                                       histtype='step')
    x = np.linspace(0, 1, 100)
    y = a * x**(a - 1.0)
    normed_y1 = samples * np.diff(bins1)[0] * y
    plt.plot(x, normed_y1, label='numpy.random.power fit')

    count2, bins2, ignored2 = plt.hist(s2,
                                       bins=30,
                                       label='joinmarket',
                                       histtype='step')
    normed_y2 = samples * np.diff(bins2)[0] * y
    plt.plot(x, normed_y2, label='common.rand_pow_array fit')
    plt.title('testing power distribution')
    plt.legend(loc='upper left')
    plt.show()

def make_intergenerational_figure(data, lowerbound, upperbound, rows, title):
    plt.figure(figsize=(10,10))
    plt.suptitle(title,fontsize=20)
    for index in range(4):
        plt.subplot(2,2,index+1)    
        #simulation distribution
        plt.hist(accepted[:,rows[index]], normed=True, bins = range(0,100,5), color = col)
        #simulation values
        value = np.mean(accepted[:,rows[index]])
        std = 2*np.std(accepted[:,rows[index]])
        plt.errorbar((value,), (red_marker_location-0.02), xerr=((std,),(std,)),
                     color=col, fmt='o', linewidth=2, capsize=5, mec = col)
        #survey values
        value = data[index]
        lb = lowerbound[index]
        ub = upperbound[index]
        plt.errorbar((value,), (red_marker_location,), xerr=((value-lb,),(ub-value,)),
                     color='r', fmt='o', linewidth=2, capsize=5, mec = 'r')
        #labeling    
        plt.ylim(0,ylimit)
        plt.xlim(0,100)
    #make subplots pretty
    plt.subplot(2,2,1)
    plt.title("Males")
    plt.ylabel("'05\nFrequency")
    plt.subplot(2,2,2)
    plt.title("Females")
    plt.subplot(2,2,3)
    plt.ylabel("'08\nFrequency")
    plt.xlabel("Percent Responding Affirmatively")
    plt.subplot(2,2,4)
    plt.xlabel("Percent Responding Affirmatively")

    def testProbabilites(self):

        pmf=PMFList(10)
        count={}
        times=10000
        exp=re.compile('\d+')

        o=[]
        for i in range(0,times):
            
            sel=pmf.choose()
            o.append(int(exp.search(sel[0]).group()))

            try: 
                count[sel]=count[sel]+1
            except KeyError:
                count[sel]=1

        print("-------------------------------\n"+
              "Results table of the PMF test:\n"+
              "-------------------------------")

        for k,i in count.items():
            print("Item name: "+k[0]+" | Item probability: "+
                  k[1].__str__()+" | Ocurrences: "+i.__str__()+
                  " | Item empirical probability ("+
                  times.__str__()+" runs)"+
                  ": "+(float(i)/times).__str__())
            

        plt.hist(o)
        plt.xlabel('Item number (1-'+len(count).__str__()+')')
        plt.ylabel('Number of occurrences')
        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 main():
    # produce gaussian noise and show as standard color coded image
    ar = np.random.randn(100, 200)
    plt.imshow(ar)
    plt.colorbar(shrink=.5)
    plt.title('Gaussian noise color coded')
    plt.show()

    # show as grayscale coded image
    plt.imshow(ar, cmap=cm.gray)
    plt.colorbar(shrink=.5)
    plt.title('Gaussian noise grayscale coded')
    plt.show()

    # flatten and show as histogram
    ar_flat = ar.reshape(100 * 200)
    plt.hist(ar_flat)
    plt.title('Gaussian noise histogram')
    plt.show()

    # show sin(x)*sin(y) sample
    l = 100
    sin_x = np.sin(range(l))
    array = np.asarray([sin_x * x for x in sin_x])
    plt.imshow(array, interpolation = None)
    plt.colorbar()
    plt.title('sin(x)sin(y)')
    plt.show()

    # replace al negative values by 0
    array_pos = np.maximum(array, 0)
    plt.imshow(array_pos, interpolation = None)
    plt.colorbar()
    plt.title('sin(x)sin(y) only positive')
    plt.show()

def plot_tothist(infile, tot, maxy, binsize=3):
   """
   Plot the total-score histogram, where the total score (tot) has been
   previous calculated or read-in by the input functions
   """

   """ Calculate moments of the distribution """
   mn = tot.mean()
   med = np.median(tot)
   mp = tot.mean() + tot.std()
   mm = tot.mean() - tot.std()

   """ Report on the properties of the distibution """
   print ""
   print "Statistics for %s" % infile
   print "---------------------------------"
   print "  Mean:         %5.1f" % mn
   print "  Median:       %5.1f" % med
   print "  Sigma:        %5.1f" % tot.std()
   print "  Mean - 1 sig: %5.1f" % mm
   print "  Mean + 1 sig: %5.1f" % mp
   print ""

   """ Plot the distribution """
   binhist = range(int(tot.min())-1,int(tot.max())+3,binsize)
   plt.hist(tot,binhist,histtype='step',ec='k')
   plt.ylim(0,maxy)
   plt.axvline(x=mn, ymin=0, ymax=maxy, c='r', lw=3)
   plt.axvline(x=mm, ymin=0, ymax=maxy, c='b', lw=3)
   plt.axvline(x=mp, ymin=0, ymax=maxy, c='b', lw=3)
   plt.title("Distribution of scores for %s" % infile)
   plt.xlabel("Scores")
   plt.ylabel("N")
   plt.show()