def simFlips(numFlips, numTrials):      # performs and displays the simulation result
    diffs = []                          # diffs to know if there was a fair Trial. It has the absolute differences of heads and tails in each trial
    for i in xrange(0, numTrials):           
        heads, tails = flipTrial(numFlips)
        diffs.append(abs(heads - tails))
                
    diffs = pylab.array(diffs)          # create an array of diffs
    diffMean = sum(diffs)/len(diffs)    # average of absolute differences of heads and tails from each trial
    diffPercent = (diffs/float(numFlips)) * 100     # create an array of percentage of each diffs from its no. of flips.
    percentMean = sum(diffPercent)/len(diffPercent)     # create a percent mean of all diffPercents in the array
    
    pylab.hist(diffs)                   # displays the distribution of elements in diffs array
    pylab.axvline(diffMean, color = 'r', label = 'Mean')
    pylab.legend()
    titleString = str(numFlips) + ' Flips, ' + str(numTrials) + ' Trials'
    pylab.title(titleString)
    pylab.xlabel('Difference between heads and tails')
    pylab.ylabel('Number of Trials')
    
    pylab.figure()
    pylab.plot(diffPercent)
    pylab.axhline(percentMean, color = 'r', label = 'Mean')
    pylab.legend()
    pylab.title(titleString)
    pylab.xlabel('Trial Number')
    pylab.ylabel('Percent Difference between heads and tails')

def simulationDelayedTreatment(numTrials, condition=75):
    """
    Runs simulations and make histograms for problem 1.

    Runs numTrials simulations to show the relationship between delayed
    treatment and patient outcome using a histogram.

    Histograms of final total virus populations are displayed for delays of 300,
    150, 75, 0 timesteps (followed by an additional 150 timesteps of
    simulation).

    numTrials: number of simulation runs to execute (an integer)
    """
    
    trialResults = {trialNum: 0 for trialNum in range(numTrials)}
    for trial in range(numTrials):
        viruses = [ResistantVirus(0.1, 0.05, {'guttagonol': False}, 0.005) for x in range(100)]
        treatedPatient = TreatedPatient(viruses, 1000)
        for timeStep in range(0,condition+150):
            treatedPatient.update()
            if timeStep == condition:
                treatedPatient.addPrescription('guttagonol')
        print str(trial) + " Completed"
        trialResults[trial] = treatedPatient.update()
    print trialResults
    pylab.hist(trialResults.values(), bins=20)
    pylab.title("Final Resistant Population - Prescription Given After " + str(condition) + " Time Steps for " + str(numTrials) + " Trials")
    pylab.xlabel("Final Total Virus Population")
    pylab.ylabel("Number of Trials")
    pylab.legend(loc='best')
    pylab.show()

def plotHistogram(data, preTime):
    pylab.figure(1)
    pylab.hist(data, bins=10)
    pylab.xlabel("Virus Population At End of Simulation")
    pylab.ylabel("Number of Trials")
    pylab.title("{0} Time Steps Before Treatment Simulation".format(preTime))
    pylab.show()

def stage_plots2(sdsscoimgs=None, coimgs=None,
                 comods=None, comods2=None, resids=None,
                 bands=None,
                 **kwargs):

    for band,co in zip(bands, sdsscoimgs):
        print 'co', co.shape
        plt.clf()
        plt.hist(co.ravel(), range=(-0.1, 0.1), bins=100)
        plt.title('SDSS %s band' % band)
        plt.savefig('sdss-%s.png' % band)

        print band, 'band 16th and 84th pcts:', np.percentile(co.ravel(), [16,84])
        
    kwa = dict(mnmx=(-2,10), scales=dict(g=(2,0.02), r=(1,0.03),
    z=(0,0.1)))
    #z=(0,0.22)))
    

    plt.clf()
    dimshow(get_rgb(sdsscoimgs, bands, **kwa), ticks=False)
    plt.savefig('sdss2.png')

    plt.clf()
    dimshow(get_rgb(coimgs, bands, **kwa), ticks=False)
    plt.savefig('img2.png')

def hist_Ne(sel_coinc_ids, coords, data, n):
    global histNe2
    c_index = data.root.coincidences.c_index
    observ = data.root.coincidences.observables
    core_rec = data.root.core_reconstructions.reconstructions
    
    histNe = core_rec.readCoordinates(coords, field='reconstructed_shower_size')
    #histNe = [x for x in histNe if x > 0]  # for showersize smaller than 0 
    
    d = 10**10.2
    
    histNe2 = [x*d for x in histNe]
    #histNe *= 10 ** 10.2
    
      
    pylab.hist(np.log10(histNe2), 100, log=True) # histtype="step"
       
    pylab.xlabel('Showerenergy log(eV)')
    pylab.ylabel('count')
    pylab.title('showersize bij N==%s' %n)
    pylab.ylim(ymin=1)
    pylab.grid(True)
    pylab.show()
    
    return histNe2

def plotHist(result, title, xLabel, yLabel):
    pylab.hist(result)
    pylab.title(title)
    pylab.xlabel(xLabel)
    pylab.ylabel(yLabel)
    pylab.legend(loc = 1)
    pylab.show()

def make_cdf(ham, spam, similarity_name, file_name=None):
    hist(
        ham,
        alpha=0.5,
        bins=100,
        normed=True,
        cumulative=True,
        histtype='stepfilled',
        label='Ham (' + str(len(ham)) + ')',
        color='b')
    hist(
        spam,
        alpha=0.5,
        bins=100,
        normed=True,
        cumulative=True,
        histtype='stepfilled',
        label='Spam (' + str(len(spam)) + ')',
        color='r')
    legend(loc=2)
    title("CDF of %s similarity for spam and ham" % (similarity_name))
    xlabel("%s similarity" % (similarity_name))
    ylabel("proportion")
    if file_name is None:
        savefig("%s_cdf.png" % (similarity_name))
    else:
        savefig(file_name)
    show()
    clf()

def plotVowelProportionHistogram(wordList, numBins=15):
    """
    Plots a histogram of the proportion of vowels in each word in wordList
    using the specified number of bins in numBins
    """
    vowels = 'aeiou'
    vowelProportions = []

    for word in wordList:
        vowelsCount = 0.0

        for letter in word:
            if letter in vowels:
                vowelsCount += 1

        vowelProportions.append(vowelsCount / len(word))

    meanProportions = sum(vowelProportions) / len(vowelProportions)
    print "Mean proportions: ", meanProportions

    pylab.figure(1)
    pylab.hist(vowelProportions, bins=15)
    pylab.title("Histogram of Proportions of Vowels in Each Word")
    pylab.ylabel("Count of Words in Each Bucket")
    pylab.xlabel("Proportions of Vowels in Each Word")

    ymin, ymax = pylab.ylim()
    ymid = (ymax - ymin) / 2
    pylab.text(0.03, ymid, "Mean = {0}".format(
        str(round(meanProportions, 4))))
    pylab.vlines(0.5, 0, ymax)
    pylab.text(0.51, ymax - 0.01 * ymax, "0.5", verticalalignment = 'top')

    pylab.show()

def simulationDelayedTreatment():
    """
    Runs simulations and make histograms for problem 5.
    Runs multiple simulations to show the relationship between delayed treatment
    and patient outcome.
    Histograms of final total virus populations are displayed for delays of 300,
    150, 75, 0 timesteps (followed by an additional 150 timesteps of
    simulation).    
    """
    histBins = [i*50 for i in range(11)]
    delays = [0, 75, 150, 300]
    subplot = 1
    for d in delays:
        results = []
        for t in xrange(NUM_TRIALS):
            results.append(runSimulation(d))
        pylab.subplot(2, 2, subplot)
        subplot += 1
        pylab.hist(results, bins=histBins, label='delayed ' + str(d))
        pylab.xlim(0, 500)
        pylab.ylim(0, NUM_TRIALS)
        pylab.ylabel('number of patients')
        pylab.xlabel('total virus population')
        pylab.title(str(d) + ' time step delay')
        popStd = numpy.std(results)
        popMean = numpy.mean(results)
        ## print str(d)+' step delay standard deviation: '+str(popStd)
        ## print str(d)+' step mean: '+str(popMean)
        ## print str(d)+' step CV: '+str(popStd / popMean)
        ## print str(d) + ' step delay: ' + str(results)
    pylab.suptitle('Patient virus populations after 150 time steps when ' +\
                   'prescription\n' +\
                   'is applied after delays of 0, 75, 150, 300 time steps')
    pylab.show()    

    def test_wald_sample(self):
        acc=ShiftedWaldAccumulator(.2, .2, 2.0)
        nsamples=100000
        x=np.linspace(0,10, nsamples)
        
        import pylab as pl
        samp=acc.sample(nsamples)
        #dens=scipy.stats.gaussian_kde(samp[samp<10])

        pl.hist(acc.sample(nsamples),200, normed=True)
        h,hx=np.histogram(samp, density=True, bins=1000)
        hx=hx[:-1]+(hx[1]-hx[0])/2.
        #assert np.all(np.abs(h-acc.pdf(hx))<1.5)

        # kolmogoroff smirnov tests whether samples come from CDF
        D,pv=scipy.stats.kstest(samp, acc.cdf)
        print D,pv
        assert pv>.05, "D=%f,p=%f"%(D,pv)
        if True:
            pl.clf()
            #pl.subplot(2,1,1)
            #pl.hist(samp[samp<10],300, normed=True, alpha=.3)


            #pl.subplot(2,1,2)
            pl.bar(hx, h, alpha=.3, width=hx[1]-hx[0])
            pl.plot(x,acc.pdf(x), color='red', label='analytical')
            #pl.plot(x,dens(x),    color='green', label='kde')
            pl.xlim(0,3)
            pl.legend()
            self.savefig()

def distribution_plot(item_array, index_num):
    vals = [x[index_num] for x in item_array]
    #vals = []
    #for db_domain in cur_db.keys():
    #    vals.append(cur_db[db_domain][item_key])
    pylab.hist(vals, 20)
    pylab.show()

def simulationDelayedTreatment(numTrials):
    """
    Runs simulations and make histograms for problem 1.

    Runs numTrials simulations to show the relationship between delayed
    treatment and patient outcome using a histogram.

    Histograms of final total virus populations are displayed for delays of 300,
    150, 75, 0 timesteps (followed by an additional 150 timesteps of
    simulation).

    numTrials: number of simulation runs to execute (an integer)
    """
    # TODO
    patients =  [TreatedPatient(viruses=[ResistantVirus(maxBirthProb=0.1,
                                                        clearProb=0.05,
                                                        resistances={'guttagonol': False},
                                                        mutProb=0.005) 
                                        for x in xrange(100)],
                                maxPop=1000) for i in xrange(numTrials)]
    delay = 75
    viruspops = []
    for p in patients:
        for i in xrange(delay):
            p.update() 
        p.addPrescription('guttagonol')
        for i in xrange(150):
            p.update()
        viruspops.append(p.getTotalPop())

    pylab.hist(viruspops,10)
    pylab.show()

    def plot_function_data_scientific(self):
        """plot the crap in a more scientific way!
           more information here:
           http://matplotlib.sourceforge.net/api/pyplot_api.html#matplotlib.pyplot.plot
        """

        print '[main] plot data'

        function_name_to_plot_regex = re.compile('^.*(flip|Critical).*$')

        #size in inches, come on...
        pylab.figure(figsize = (32, 18))

        #plot data
        for element in self.__function_name_value_tuple_list:
            if function_name_to_plot_regex.match(element[0]):
                #for time_value in element[1]:
                    #counter_dict = collections.Counter(element[1])

                pylab.hist(element[1], label = element[0], histtype = 'bar')
                #pylab.plot(counter_dict.keys(), counter_dict.values(), 'O', label = element[0])

        #decorate plot
        pylab.xlabel('used time in function [ns]')
        pylab.ylabel('value count')
        pylab.title('JNI ByteFlipper Performance Test (size of byte array: %d)' % (self.__byte_array_size))
        pylab.grid(True)
        pylab.legend(loc = 'best')

        plot_file_name = self.log_file_name.strip('.log') + '_scientific.png'
        print '[main] write simple plot to image (%s)' % plot_file_name
        pylab.savefig(plot_file_name)

def simulationTwoDrugsDelayedTreatment(numTrials):
    """
    Runs simulations and make histograms for problem 2.

    Runs numTrials simulations to show the relationship between administration
    of multiple drugs and patient outcome.

    Histograms of final total virus populations are displayed for lag times of
    300, 150, 75, 0 timesteps between adding drugs (followed by an additional
    150 timesteps of simulation).

    numTrials: number of simulation runs to execute (an integer)
    """
    # simulation incorrect...
    
    wait_time = (300, 150, 75, 0)
    finals = [[] for i in wait_time]
    for wait in wait_time:
        results = []
        print "\ntrial (wait time = " + str(wait) + ") working",
        for i in range(numTrials):
            tmp = simulationWithDrug(100, 1000, 0.1, 0.05, {'guttagonol': False},
                                     0.005, 1, before_drug=150, after_drug=wait, plot=False)
            sec_tmp = simulationWithDrug(int(tmp[-1]), 1000, 0.1, 0.05, {'grimpex': False},
                                         0.005, 1, before_drug=0, after_drug=150, plot=False)
            results.append(tmp[-1])
            if i%2: print ".",
        finals[wait_time.index(wait)] += results
    for results in finals:    
        pylab.hist(results, [x for x in range(1000) if x % 50 == 0])
        pylab.title("Wait Time = " + str(wait_time[finals.index(results)]))
        pylab.xlabel("Final Virus Populations")
        pylab.ylabel("Trials")
        pylab.ylim(0, numTrials)
        pylab.show()

def plot_distribution_true_false(prediction_df):
    """

    :param prediction_df:
    :return:
    """
    mask_well_classified = prediction_df.expected == prediction_df.class_predicted
    for group_name, g in prediction_df.groupby('expected'):
        print("group_name %s" % group_name)
        pylab.figure()
        try:
            v = g.confidence[mask_well_classified]
            pylab.hist(list(v), color='g', alpha=0.3, normed=0, range=(0,1), bins=10)
            #sns.distplot(g.confidence[mask_well_classified], color='g', bins=11)  # TRUE POSITIVE
        except Exception as e:
            print(e)
        mask_wrong = (prediction_df.class_predicted == group_name) & (prediction_df.expected != group_name)  # FALSE POSITIVE
        #v = g.confidence[~mask_well_classified]
        try:
            v = prediction_df.confidence[mask_wrong]
            pylab.hist(list(v), color='r', alpha=0.3, normed=0, range=(0,1), bins=10)
            #sns.distplot(v, color='r', bins=11)
        except Exception as e:
            print(e)
            #print(len(v))
            pass
        print("FIN figure %s" % group_name)
        pylab.show()
        print("")

def hist_scores(scores, labels, weights=None, **kwargs):

    if hasattr(scores, "values"):
        scores = scores.values
    if hasattr(labels, "values"):
        labels = labels.values
    if hasattr(weights, "values"):
        weights = weights.values

    if weights != None:
        s_w = weights[labels == 1]
        b_w = weights[labels == 0]
    else:
        s_w = b_w = None

    kwargs["bins"] = kwargs.get("bins", 100)
    kwargs["normed"] = kwargs.get("normed", True)
    kwargs["histtype"] = kwargs.get("histtype", "stepfilled")

    if scores.ndim == 2:
        scores = scores[:, 1]

    plt.hist(scores[labels == 1], alpha=0.5, label="signal", weights=s_w, color="blue", **kwargs)
    plt.hist(scores[labels == 0], alpha=0.5, label="bkgd", weights=b_w, color="red", **kwargs)
    plt.legend(loc="best")

    def plot(self):
        rows = np.sqrt(self.M).round()
        cols = np.ceil(self.M/rows)

        if self.D==1:
            xmin = self.X.min()
            xmax = self.X.max()
            xmin,xmax = xmin-0.1*(xmax-xmin), xmax+0.1*(xmax-xmin)
            Xgrid = np.linspace(xmin,xmax,100)[:,None]
            zz = self.predict(Xgrid)

        if self.D==2:
            xmin,ymin = np.vstack(self.X).min(0)
            xmax,ymax = np.vstack(self.X).max(0)
            xmin,xmax = xmin-0.1*(xmax-xmin), xmax+0.1*(xmax-xmin)
            ymin,ymax = ymin-0.1*(ymax-ymin), ymax+0.1*(ymax-ymin)
            xx,yy = np.mgrid[xmin:xmax:100j,ymin:ymax:100j]
            Xgrid = np.vstack((xx.flatten(),yy.flatten())).T
            zz = self.predict(Xgrid).reshape(100,100,self.M)

        for m in range(self.M):
            pb.subplot(rows,cols,m+1)
            if self.D==1:
                pb.hist(self.X[m,:,0],self.N/10.,normed=True)
                pb.plot(Xgrid,zz[:,m],'r',linewidth=2)
            elif self.D==2:
                pb.plot(self.X[m,:,0],self.X[m,:,1],'rx',mew=2)
                zz_data = self.predict(self.X[m])[:,m]
                pb.contour(xx,yy,zz[:,:,m],[stats.scoreatpercentile(zz_data,5)],colors='r',linewidths=1.5)
                pb.imshow(zz[:,:,m].T,extent=[xmin,xmax,ymin,ymax],origin='lower',cmap=pb.cm.binary,vmin=0.,vmax=zz_data.max())

def plotQuizzes():
  scores = generateScores(10000)
  pylab.hist(scores, bins=7)
  pylab.xlabel("Final Score")
  pylab.ylabel("Number of Trials")
  pylab.title("Distribution of Scores")
  pylab.show()

def explore_city_data(city_data):
    """Calculate the Boston housing statistics."""

    # Get the labels and features from the housing data
    housing_prices = city_data.target
    housing_features = city_data.data

    ###################################
    ### Step 1. YOUR CODE GOES HERE ###
    ###################################

    # Please calculate the following values using the Numpy library
    print "Statistics of the Boston housing dataset."
    print
    # Size of data (number of houses)?
    print "Number of houses: ", housing_features.shape[0]
    # Number of features?
    print "Number of features: ", housing_features.shape[1]
    # Minimum price?
    print "Minimum price (in $1,000s): ", min(housing_prices)
    # Maximum price?
    print "Maximum price (in $1,000s): ", max(housing_prices)
    # Calculate mean price?
    print "Mean price (in $1,000s): ", housing_prices.mean() 
    # Calculate median price?
    print "Median price (in $1,000s): ", np.median(housing_prices)
    # Calculate standard deviation?
    print "Standard deviation (in $1,000s): ", housing_prices.std()

    # Make a plot?
    pl.figure()
    pl.hist(housing_prices, bins=15)
    pl.xlabel('Price in $1000s')
    pl.ylabel('Counts')
    pl.show()

def pair_plot(data, savefile=None, display=True, **kwargs):
    chan = data.channels
    l = len(chan)
    figure = pylab.figure()
    pylab.subplot(l, l, 1)
    for i in range(l):
        for j in range(i + 1):
            pylab.subplot(l, l, i * l + j + 1)
            if i == j:
                pylab.hist(data[:, i], bins=200, histtype='stepfilled')
            else:
                pylab.scatter(data[:, i], data[:, j], **kwargs)

            if j == 0:
                pylab.ylabel(chan[i])
            if i == l - 1:
                pylab.xlabel(chan[j])

    if display:
        pylab.show()

    if savefile:
        pylab.savefig(savefile)

    return figure