def main():
    SAMPLE_NUM = 10
    degree = 9
    x, y = sin_wgn_sample(SAMPLE_NUM)
    fig = pylab.figure(1)
    pylab.grid(True)
    pylab.xlabel('x')
    pylab.ylabel('y')
    pylab.axis([-0.1,1.1,-1.5,1.5])

    # sin(x) + noise
    # markeredgewidth mew
    # markeredgecolor mec
    # markerfacecolor mfc

    # markersize      ms
    # linewidth       lw
    # linestyle       ls
    pylab.plot(x, y,'bo',mew=2,mec='b',mfc='none',ms=8)

    # sin(x)
    x2 = linspace(0, 1, 1000)
    pylab.plot(x2,sin(2*x2*pi),'#00FF00',lw=2,label='$y = \sin(x)$')

    # polynomial fit
    reg = exp(-18)
    w = curve_poly_fit(x, y, degree,reg) #w = polyfit(x, y, 3)
    po = poly1d(w)      
    xx = linspace(0, 1, 1000)
    pylab.plot(xx, po(xx),'-r',label='$M = 9, \ln\lambda = -18$',lw=2)
    
    pylab.legend()
    pylab.show()
    fig.savefig("poly_fit9_10_reg.pdf")

def TestOverIntDim():
    nDim = 10
    numOfParticles = 20
    maxIteration = 200
    minX = array([-100.0]*nDim)
    maxX = array([100.0]*nDim)
    maxV = 0.2*(maxX - minX)
    minV = -1.0*maxV
    numOfTrial = 20
    alpha = 0.0
    for intDim in xrange(0,11,2):
        gBest = array([0.0]*maxIteration)
        for i in xrange(numOfTrial):
            p1 = AUPSO.PSOProblem(nDim, numOfParticles, maxIteration, minX, maxX, minV, maxV, AUPSO.Griewank,intDim,alpha)
            p1.run()
            gBest = gBest + p1.gBestArray[:maxIteration]
        gBest = gBest / numOfTrial
        pylab.plot(range(maxIteration), log10(gBest),label='intDim='+str(intDim))
    pylab.title('$G_{best}$ over 20 trials'+' alpha='+str(alpha))
    pylab.xlabel('The $N^{th}$ Iteratioin')
    pylab.ylabel('Average gBest over '+str(numOfTrial)+' runs')
    pylab.grid(True)
#    pylab.yscale('log')
    ylim = [-6, 1]
    ystep = 1.0
#    pylab.ylim(ylim[0], ylim[1])
#    yticks = linspace(ylim[0], ylim[1], int((ylim[1]-ylim[0])/ystep+1))
#    pylab.yticks(tuple(yticks), tuple(map(str,yticks)))
    pylab.legend(loc='lower left')
    pylab.show()

def param_set_averages_plot(results):
    averages_ocr = [
        a[1] for a in sorted(
            param_set_averages(results, metric='ocr').items(),
            key=lambda x: int(x[0].split('-')[1]))
    ]
    averages_q = [
        a[1] for a in sorted(
            param_set_averages(results, metric='q').items(),
            key=lambda x: int(x[0].split('-')[1]))
    ]
    averages_mse = [
        a[1] for a in sorted(
            param_set_averages(results, metric='mse').items(),
            key=lambda x: int(x[0].split('-')[1]))
    ]
    fig = plt.figure(figsize=(6, 4))
    # plt.tight_layout()
    plt.plot(averages_ocr, label='OCR', linewidth=2.0)
    plt.plot(averages_q, label='Q', linewidth=2.0)
    plt.plot(averages_mse, label='MSE', linewidth=2.0)
    plt.ylim([0, 1])
    plt.xlabel(u'Paslėptų neuronų skaičius')
    plt.ylabel(u'Vidurinė Q įverčio pokyčio reikšmė')
    plt.grid(True)
    plt.tight_layout()
    plt.legend(loc='lower right')
    plt.show()

def plotEventFlop(library, num, eventNames, sizes, times, events, filename = None):
  from pylab import legend, plot, savefig, semilogy, show, title, xlabel, ylabel
  import numpy as np

  arches = sizes.keys()
  bs     = events[arches[0]].keys()[0]
  data   = []
  names  = []
  for event, color in zip(eventNames, ['b', 'g', 'r', 'y']):
    for arch, style in zip(arches, ['-', ':']):
      if event in events[arch][bs]:
        names.append(arch+'-'+str(bs)+' '+event)
        data.append(sizes[arch][bs])
        data.append(1e-3*np.array(events[arch][bs][event])[:,1])
        data.append(color+style)
      else:
        print 'Could not find %s in %s-%d events' % (event, arch, bs)
  semilogy(*data)
  title('Performance on '+library+' Example '+str(num))
  xlabel('Number of Dof')
  ylabel('Computation Rate (GF/s)')
  legend(names, 'upper left', shadow = True)
  if filename is None:
    show()
  else:
    savefig(filename)
  return

    def test_mask_LUT(self):
        """
        The masked image has a masked ring around 1.5deg with value -10
        without mask the pixels should be at -10 ; with mask they are at 0
        """
        x1 = self.ai.xrpd_LUT(self.data, 1000)
#        print self.ai._lut_integrator.lut_checksum
        x2 = self.ai.xrpd_LUT(self.data, 1000, mask=self.mask)
#        print self.ai._lut_integrator.lut_checksum
        x3 = self.ai.xrpd_LUT(self.data, 1000, mask=numpy.zeros(shape=self.mask.shape, dtype="uint8"), dummy= -20.0, delta_dummy=19.5)
#        print self.ai._lut_integrator.lut_checksum
        res1 = numpy.interp(1.5, *x1)
        res2 = numpy.interp(1.5, *x2)
        res3 = numpy.interp(1.5, *x3)
        if logger.getEffectiveLevel() == logging.DEBUG:
            pylab.plot(*x1, label="nomask")
            pylab.plot(*x2, label="mask")
            pylab.plot(*x3, label="dummy")
            pylab.legend()
            pylab.show()
            raw_input()

        self.assertAlmostEqual(res1, -10., 1, msg="Without mask the bad pixels are around -10 (got %.4f)" % res1)
        self.assertAlmostEqual(res2, 0., 4, msg="With mask the bad pixels are actually at 0 (got %.4f)" % res2)
        self.assertAlmostEqual(res3, -20., 4, msg="Without mask but dummy=-20 the dummy pixels are actually at -20 (got % .4f)" % res3)

def Doplots_monthly(mypathforResults,PlottingDF,variable_to_fill, Site_ID,units,item):   
    ANN_label=str(item+"_NN")     #Do Monthly Plots
    print "Doing MOnthly  plot"
    #t = arange(1, 54, 1)
    NN_label='Fc'
    Plottemp = PlottingDF[[NN_label,item]][PlottingDF['day_night']!=1]
    #Plottemp = PlottingDF[[NN_label,item]].dropna(how='any')
    figure(1)
    pl.title('Nightime ANN v Tower by year-month for '+item+' at '+Site_ID)

    try:
	xdata1a=Plottemp[item].groupby([lambda x: x.year,lambda x: x.month]).mean()
	plotxdata1a=True
    except:
	plotxdata1a=False
    try:
	xdata1b=Plottemp[NN_label].groupby([lambda x: x.year,lambda x: x.month]).mean()
	plotxdata1b=True
    except:
	plotxdata1b=False 
    if plotxdata1a==True:
	pl.plot(xdata1a,'r',label=item) 
    if plotxdata1b==True:
	pl.plot(xdata1b,'b',label=NN_label)
    pl.ylabel('Flux')    
    pl.xlabel('Year - Month')       
    pl.legend()
    pl.savefig(mypathforResults+'/ANN and Tower plots by year and month for variable '+item+' at '+Site_ID)
    #pl.show()
    pl.close()
    time.sleep(1)

def plot_heatingrate(data_dict, filename, do_show=True):
    pl.figure(201)
    color_list = ['b','r','g','k','y','r','g','b','k','y','r',]
    fmtlist = ['s','d','o','s','d','o','s','d','o','s','d','o']
    result_dict = {}
    for key in data_dict.keys():
        x = data_dict[key][0]
        y = data_dict[key][1][:,0]
        y_err = data_dict[key][1][:,1]

        p0 = np.polyfit(x,y,1)
        fit = LinFit(np.array([x,y,y_err]).transpose(), show_graph=False)
        p1 = [0,0]
        p1[0] = fit.param_dict[0]['Slope'][0]
        p1[1] = fit.param_dict[0]['Offset'][0]
        print fit
        x0 = np.linspace(0,max(x))
        cstr = color_list.pop(0)
        fstr = fmtlist.pop(0)
        lstr = key + " heating: {0:.2f} ph/ms".format((p1[0]*1e3)) 
        pl.errorbar(x/1e3,y,y_err,fmt=fstr + cstr,label=lstr)
        pl.plot(x0/1e3,np.polyval(p0,x0),cstr)
        pl.plot(x0/1e3,np.polyval(p1,x0),cstr)
        result_dict[key] = 1e3*np.array(fit.param_dict[0]['Slope'])
    pl.xlabel('Heating time (ms)')
    pl.ylabel('nbar')
    if do_show:
        pl.legend()
        pl.show()
    if filename != None:
        pl.savefig(filename)
    return result_dict

def plot_dat(ax, file_name):

    with open(file_name, 'rb') as datfile:
        l=[]
        for row in datfile:
            if len(row.split('|')[-1].split()):
                l.append(row.split('|')[-1].split())
#                 print  row
    lengend_names=l[1]
    l=l[2:]
    data=[]
    for row in l:
        for i in range(len(row)):
            try:
                type=row[i][-1]
                row[i]=float(row[i][:-1])
                if type=='G':
                    row[i]*=1000.0
            except:
#                 print i
                row[i]=0.
        data.append([row[0]])
    data=zip(*data)
    data=numpy.array(data)
    shape=data.transpose().shape

    
    ax.plot(numpy.mgrid[0:shape[0]*10:10,0:1][0],
               100*(data.transpose()-data.transpose()[0,0])/(1533.0+59900.0))
    pylab.legend([lengend_names[0]])
    pylab.ylabel('Memory (MB)')
    pylab.xlabel('Time (sec)')
    pylab.show()

def PlotNCodonMuts(allmutations, plotfile, title):
    """Plots number of nucleotide changes per codon mutation.

    allmutations -> list of all mutations as tuples (wtcodon, r, mutcodon)
    plotfile -> name of the plot file we create.
    title -> string giving the plot title.
    """
    pylab.figure(figsize=(3.5, 2.25))
    (lmargin, rmargin, bmargin, tmargin) = (0.16, 0.01, 0.21, 0.07)
    pylab.axes([lmargin, bmargin, 1.0 - lmargin - rmargin, 1.0 - bmargin - tmargin])
    nchanges = {1:0, 2:0, 3:0}
    nmuts = len(allmutations)
    for (wtcodon, r, mutcodon) in allmutations:
        assert 3 == len(wtcodon) == len(mutcodon)
        diffs = len([i for i in range(3) if wtcodon[i] != mutcodon[i]])
        nchanges[diffs] += 1
    barwidth = 0.6
    xs = [1, 2, 3]
    nactual = [nchanges[x] for x in xs]
    nexpected = [nmuts * 9. / 63., nmuts * 27. / 63., nmuts * 27. / 63.]
    bar = pylab.bar([x - barwidth / 2.0 for x in xs], nactual, width=barwidth)
    pred = pylab.plot(xs, nexpected, 'rx', markersize=6, mew=3)
    pylab.gca().set_xlim([0.5, 3.5])
    pylab.gca().set_ylim([0, max(nactual + nexpected) * 1.1])
    pylab.gca().xaxis.set_major_locator(matplotlib.ticker.MaxNLocator(4))
    pylab.gca().yaxis.set_major_locator(matplotlib.ticker.MaxNLocator(5))
    pylab.xlabel('nucleotide changes in codon')
    pylab.ylabel('number of mutations')
    pylab.legend((bar[0], pred[0]), ('actual', 'expected'), loc='upper left', numpoints=1, handlelength=0.9, borderaxespad=0, handletextpad=0.4)
    pylab.title(title, fontsize=12)
    pylab.savefig(plotfile)
    time.sleep(0.5)
    pylab.show()

def main():
    amps = [0.167e-9,
            0.25e-9,
            0.333e-9]
    model_dict = setup_model()
    for ii, a in enumerate(amps):
        do_sim(model_dict['stimulus'], a)   
        config.logger.info('##### %d' % (model_dict['tab_vm'].size))
        vm = model_dict['tab_vm'].vector
        inject = model_dict['tab_stim'].vector.copy()
        t = np.linspace(0, simtime, len(vm))
        fname = 'data_fig_a3_%s.txt' % (chr(ord('A')+ii))
        np.savetxt(fname,
                   np.vstack((t, inject, vm)).transpose())
        msg = 'Saved data for %g A current pulse in %s' % (a, fname)
        config.logger.info(msg)
        print(msg)
        pylab.subplot(3,1,ii+1)
        pylab.title('%g nA' % (a*1e9))
        pylab.plot(t, vm, label='soma-Vm (mV)')
        stim_boundary = np.flatnonzero(np.diff(inject))
        pylab.plot((t[stim_boundary[0]]), (vm.min()), 'r^', label='stimulus start')                           
        pylab.plot((t[stim_boundary[-1]]), (vm.min()), 'gv', label='stimulus end')        
        pylab.legend()    
    pylab.savefig('fig_a3.png')
    pylab.show()

def simulationWithoutDrug(numViruses, maxPop, maxBirthProb, clearProb,
                          numTrials):
    """
    Run the simulation and plot the graph for problem 3 (no drugs are used,
    viruses do not have any drug resistance).    
    For each of numTrials trial, instantiates a patient, runs a simulation
    for 300 timesteps, and plots the average virus population size as a
    function of time.

    numViruses: number of SimpleVirus to create for patient (an integer)
    maxPop: maximum virus population for patient (an integer)
    maxBirthProb: Maximum reproduction probability (a float between 0-1)        
    clearProb: Maximum clearance probability (a float between 0-1)
    numTrials: number of simulation runs to execute (an integer)
    """

    # TODO
    steps = 300
    trialResults = [[] for s in range(steps)]

    for i in range(numTrials): 
        viruses = [SimpleVirus(maxBirthProb,clearProb) for v in range(numViruses)]
        patient = Patient(viruses, maxPop)
        for step in range(300): 
            trialResults[step].append(patient.update())
    resultsSummary = [sum(l) / float(numTrials) for l in trialResults]
    pylab.plot(resultsSummary, label="Total Virus Population")
    pylab.title("SimpleVirus simulation")
    pylab.xlabel("Time Steps")
    pylab.ylabel("Average Virus Population")
    pylab.legend()
    pylab.show()

def plot_datasets(dataset_ids, title=None, legend=True, labels=True):
    """
    Plots one or more dataset.

    :param dataset_ids: list of datasets to plot
    :type dataset_ids: list of integers
    :param title: title of the plot
    :type title: string
    :param legend: whether or not to show legend
    :type legend: boolean
    :param labels: whether or not to plot point labels
    :type labels: boolean
    """
    title = title if title else "Datasets " + ",".join(
        [str(d) for d in dataset_ids])
    pl.title(title)

    data = {k: v for k, v in npoints.items() if k in dataset_ids}

    lines = [pl.plot(zip(*p)[0], zip(*p)[1], 'o-')[0] for p in data.values()]

    if legend:
        pl.legend(lines, data.keys())

    if labels:
        for x, y, l in [i for s in data.values() for i in s]:
            pl.annotate(str(l), xy=(x, y), xytext=(x, y + 0.1))

    pl.grid(True)

    return pl

def plotMonthlyTrend(keywords, title, monthList):
    db = mysql(host, user, passwd, dbName)
    db.connect()
    allKeywordTrend = []
    for k in keywords:
        allCount = []
        for m in monthList:
            rows = db.getMonthlyKeywordCount(k, m)
            print rows
            count = 0
            for r in rows:
                count += r[0]
            persent = count*1.0
            cc = db.getMonthlyTweetCount(m)
            if cc == 0:
                persent = 0.0
            else:
                persent /= cc
            allCount.append(persent)
        allKeywordTrend.append(allCount)
    db.close()

    for p in allKeywordTrend:
        pylab.plot(range(1, len(p)+1), p)
    pylab.title(title)
    pylab.legend(keywords)
    pylab.xlabel("month")
    pylab.ylabel("frequency of occurrence")
    pylab.show()

def RosenbrockTest():
    nDim = 3
    numOfParticles = 20
    maxIteration = 200
    minX = array([-5.0]*nDim)
    maxX = array([5.0]*nDim)
    maxV = 0.2*(maxX - minX)
    minV = -1.0*maxV
    numOfTrial = 20
    gBest = array([0.0]*maxIteration)
    for i in xrange(numOfTrial):
        p1 = RPSO.PSOProblem(nDim, numOfParticles, maxIteration, minX, maxX, minV, maxV, RPSO.Rosenbrock)
        p1.run()
        gBest = gBest + p1.gBestArray[:maxIteration]
    gBest = gBest / numOfTrial
    pylab.title('$G_{best}$ over 20 trials')
    pylab.xlabel('The $N^{th}$ Iteratioin')
    pylab.ylabel('Average gBest over '+str(numOfTrial)+' runs (logscale)')
    pylab.grid(True)
#    pylab.yscale('log')
    ymin, ymax = -1.5, 2.5
    ystep = 0.5
    pylab.ylim(ymin, ymax)
    yticks = linspace(ymin, ymax, (ymax-ymin)/ystep+1)
    pylab.yticks(tuple(yticks),tuple(map(str,yticks)))
    pylab.plot(range(maxIteration), log10(gBest),'-', label='Global best')
    pylab.legend()
    pylab.show()

    def rmsdSpreadSubplot(multiplier=1.0, layout=(-1, -1)):
        rmsd_data   = dict( (e, rad_data[e]['innov'][quant])  for e in rad_data.iterkeys() )
        spread_data = dict( (e, rad_data[e]['spread'][quant]) for e in rad_data.iterkeys() )

        times = temp.getTimes()
        n_t = len(times)

        for exp, exp_name in exp_names.iteritems():
            pylab.plot(sawtooth(times, times)[:(n_t + 1)], rmsd_data[exp][:(n_t + 1)], color=colors[exp], linestyle='-')
            pylab.plot(times[(n_t / 2):], rmsd_data[exp][n_t::2], color=colors[exp], linestyle='-')
 
        for exp, exp_name in exp_names.iteritems():
            pylab.plot(sawtooth(times, times)[:(n_t + 1)], spread_data[exp][:(n_t + 1)], color=colors[exp], linestyle='--')
            pylab.plot(times[(n_t / 2):], spread_data[exp][n_t::2], color=colors[exp], linestyle='--')

        ylim = pylab.ylim()
        pylab.plot(times, -1 * np.ones((len(times),)), color='#999999', linestyle='-', label="RMS Innovation")
        pylab.plot(times, -1 * np.ones((len(times),)), color='#999999', linestyle='--', label="Spread")

        pylab.axhline(y=7, color='k', linestyle=':')
        pylab.axvline(x=14400, color='k', linestyle=':')

        pylab.ylabel("RMS Innovation/Spread (dBZ)", size='large')

        pylab.xlim(times[0], times[-1])
        pylab.ylim(ylim)

        pylab.legend(loc=4)

        pylab.xticks(times[::2], [ "" for t in times[::2] ])
        pylab.yticks(size='x-large')
        return

def plotForce():
    figure(size=3,aspect=0.5)
    subplot(1,2,1)
    from EvalTraj import plotFF
    plotFF(vp=351,t=28,f=900,cm=0.6,foffset=8)
    subplot_annotate()
    
    subplot(1,2,2)
    for i in [1,2,3,4]:
        R=np.squeeze(np.load('Rdpse%d.npy'%i))
        R=stats.nanmedian(R,axis=2)[:,1:,:]
        dps=np.linspace(-1,1,201)[1:]
        plt.plot(dps,R[:,:,2].mean(0));
    plt.legend([0,0.1,0.2,0.3],loc=3) 
    i=2
    R=np.squeeze(np.load('Rdpse%d.npy'%i))
    R=stats.nanmedian(R,axis=2)[:,1:,:]
    mn=np.argmin(R,axis=1)
    y=np.random.randn(mn.shape[0])*0.00002+0.0438
    plt.plot(np.sort(dps[mn[:,2]]),y,'+',mew=1,ms=6,mec=[ 0.39  ,  0.76,  0.64])
    plt.xlabel('Displacement of Force Origin')
    plt.ylabel('Average Net Force Magnitude')
    hh=dps[mn[:,2]]
    err=np.std(hh)/np.sqrt(hh.shape[0])*stats.t.ppf(0.975,hh.shape[0])
    err2=np.std(hh)/np.sqrt(hh.shape[0])*stats.t.ppf(0.75,hh.shape[0])
    m=np.mean(hh)
    print m, m-err,m+err
    np.save('force',[m, m-err,m+err,m-err2,m+err2])
    plt.xlim([-0.5,0.5])
    plt.ylim([0.0435,0.046])
    plt.grid(b=True,axis='x')
    subplot_annotate()

def showExamplePolyFit(xs,ys,fitDegree1 = 1,fitDegree2 = 2):
    pylab.figure()    
    pylab.plot(xs,ys,'r.',ms=2.0,label = "measured")

    # poly fit to noise
    coeeff = numpy.polyfit(xs, ys, fitDegree1)

    # Predict the curve
    pys = numpy.polyval(numpy.poly1d(coeeff), xs)

    se = mse(ys, pys)
    r2 = rSquared(ys, pys)

    pylab.plot(xs,pys, 'g--', lw=5,label="%d degree fit, SE = %0.10f, R2 = %0.10f" %(fitDegree1,se,r2))

    # Poly fit to noise
    coeeffs = numpy.polyfit(xs, ys, fitDegree2)

    # Predict the curve
    pys = numpy.polyval(numpy.poly1d(coeeffs), xs)

    se = mse(ys, pys)
    r2 = rSquared(ys, pys)

    pylab.plot(xs,pys, 'b--', lw=5,label="%d degree fit, SE = %0.10f, R2 = %0.10f" %(fitDegree2,se,r2))

    pylab.legend()

    def plot_cost(self):
        if self.show_cost not in self.train_outputs[0][0]:
            raise ShowNetError("Cost function with name '%s' not defined by given convnet." % self.show_cost)
        train_errors = [o[0][self.show_cost][self.cost_idx] for o in self.train_outputs]
        test_errors = [o[0][self.show_cost][self.cost_idx] for o in self.test_outputs]

        numbatches = len(self.train_batch_range)
        test_errors = numpy.row_stack(test_errors)
        test_errors = numpy.tile(test_errors, (1, self.testing_freq))
        test_errors = list(test_errors.flatten())
        test_errors += [test_errors[-1]] * max(0,len(train_errors) - len(test_errors))
        test_errors = test_errors[:len(train_errors)]

        numepochs = len(train_errors) / float(numbatches)
        pl.figure(1)
        x = range(0, len(train_errors))
        pl.plot(x, train_errors, 'k-', label='Training set')
        pl.plot(x, test_errors, 'r-', label='Test set')
        pl.legend()
        ticklocs = range(numbatches, len(train_errors) - len(train_errors) % numbatches + 1, numbatches)
        epoch_label_gran = int(ceil(numepochs / 20.)) # aim for about 20 labels
        epoch_label_gran = int(ceil(float(epoch_label_gran) / 10) * 10) # but round to nearest 10
        ticklabels = map(lambda x: str((x[1] / numbatches)) if x[0] % epoch_label_gran == epoch_label_gran-1 else '', enumerate(ticklocs))

        pl.xticks(ticklocs, ticklabels)
        pl.xlabel('Epoch')
#        pl.ylabel(self.show_cost)
        pl.title(self.show_cost)

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 plot_sphere_x( s, fname ):
  """ put plot of ionization fractions from sphere `s` into fname """

  plt.figure()
  s.Edges.units = 'kpc'
  s.r_c.units = 'kpc'
  xx = s.r_c
  L = s.Edges[-1]

  plt.plot( xx, np.log10( s.xHe1 ),
            color='green', ls='-', label = r'$x_{\rm HeI}$' )
  plt.plot( xx, np.log10( s.xHe2 ),
            color='green', ls='--', label = r'$x_{\rm HeII}$' )
  plt.plot( xx, np.log10( s.xHe3 ),
            color='green', ls=':', label = r'$x_{\rm HeIII}$' )

  plt.plot( xx, np.log10( s.xH1 ),
            color='red', ls='-', label = r'$x_{\rm HI}$' )
  plt.plot( xx, np.log10( s.xH2 ),
            color='red', ls='--', label = r'$x_{\rm HII}$' )

  plt.xlim( -L/20, L+L/20 )
  plt.xlabel( 'r_c [kpc]' )

  plt.ylim( -4.5, 0.2 )
  plt.ylabel( 'log 10 ( x )' )

  plt.grid()
  plt.legend(loc='best', ncol=2)
  plt.tight_layout()
  plt.savefig( 'doc/img/x_' + fname )