def plot_funnel(pi_true, delta_str):
    delta = float(delta_str)
    n = pl.exp(mc.rnormal(10, 2**-2, size=10000))
    p = pi_true*pl.ones_like(n)

    # old way:
    #delta = delta * p * n

    nb = rate_model.neg_binom_model('funnel', pi_true, delta, p, n)
    r = nb['p_pred'].value

    pl.vlines([pi_true], .1*n.min(), 10*n.max(),
              linewidth=5, linestyle='--', color='black', zorder=10)
    pl.plot(r, n, 'o', color=colors[0], ms=10,
            mew=0, alpha=.25)

    pl.semilogy(schiz['r'], schiz['n'], 's', mew=1, mec='white', ms=15,
                color=colors[1],
                label='Observed Values')

    pl.xlabel('Rate (Per 1000 PY)', size=32)
    pl.ylabel('Study Size (PY)', size=32)
    pl.axis([-.0001, .0101, 50., 15000000])
    pl.title(r'$\delta = %s$'%delta_str, size=48)
    pl.xticks([0, .005, .01], [0, 5, 10], size=30)
    pl.yticks(size=30)

def plot_profile(experiment, step, out_file):
    from pylab import figure, subplot, hold, plot, xlabel, ylabel, text, title, axis, vlines, savefig

    if out_file is None:
        out_file = "MISMIP_%s_A%d.pdf" % (experiment, step)

    xg = x_g(experiment, step)

    x = np.linspace(0, L(), N(2) + 1)
    thk = thickness(experiment, step, x)
    x_grounded, thk_grounded = x[x < xg],  thk[x < xg]
    x_floating, thk_floating = x[x >= xg], thk[x >= xg]

    figure(1)
    ax = subplot(111)
    hold(True)
    plot(x / 1e3, np.zeros_like(x), ls='dotted', color='red')
    plot(x / 1e3, -b(experiment, x), color='black')
    plot(x / 1e3, np.r_[thk_grounded - b(experiment, x_grounded),
                        thk_floating * (1 - rho_i() / rho_w())],
         color='blue')
    plot(x_floating / 1e3, -thk_floating * (rho_i() / rho_w()), color='blue')
    _, _, ymin, ymax = axis(xmin=0, xmax=x.max() / 1e3)
    vlines(xg / 1e3, ymin, ymax, linestyles='dashed', color='black')

    xlabel('distance from the summit, km')
    ylabel('elevation, m')
    text(0.6, 0.9, "$x_g$ (theory) = %4.0f km" % (xg / 1e3),
         color='black', transform=ax.transAxes)
    title("MISMIP experiment %s, step %d" % (experiment, step))
    savefig(out_file)

def plotT(x, y, plt):
    plt.scatter(x, y)
    plt.vlines(x, [0], y)
    plt.ylim((min(y)-abs(min(y)*0.1)),max(y)+max(y)*0.1)
    plt.hlines(0, x[0]-1, x[x.shape[0]-1]+1)
    plt.xlim(x[0]-1,x[x.shape[0]-1]+1)
    plt.grid()

def draw_fit(rl, pct):
    """Draw sigmoid for psychometric

    rl: x values
    pct: y values

    Fxn draws the curve
    """
    def sig(x, A, x0, k, y0):
        return A / (1 + np.exp(-k*(x-x0))) + y0
    def sig2(x, x0, k):
        return 1. / (1+np.exp(-k*(x-x0)))

    pl.xlabel('R-L stimuli')
    pl.ylabel('p(choose R)')
    pl.xlim([rl.min()-1, rl.max()+1])
    pl.ylim([-0.05, 1.05])

    popt,pcov = curve_fit(sig, rl, pct) # stretch and yshift are free params
    popt2,pcov2 = curve_fit(sig2, rl, pct) # stretch and yshift are fixed
    x = np.linspace(rl.min(), rl.max(), 200)
    y = sig(x, *popt)
    y2 = sig2(x, *popt2)
    pl.vlines(0,0,1,linestyles='--')
    pl.hlines(0.5,rl.min(),rl.max(),linestyles='--')
    pl.plot(x,y)
    #pl.plot(x,y2)
    return popt

def decorate(mean):
    pl.legend(loc='upper center', bbox_to_anchor=(.5,-.3))
    xmin, xmax, ymin, ymax = pl.axis()
    pl.vlines([mean], -ymax, ymax*10, linestyle='dashed', zorder=20)
    pl.xticks([0, .5, 1])
    pl.yticks([])
    pl.axis([-.01, 1.01, -ymax*.05, ymax*1.01])

def noisyAbsorption(recipe, curvNo, noIndicators):
    worker = recipe.getWorker("Slicing")
    noisyAbsorption = worker.plugExtract.getResult()
    worker = recipe.getWorker("ThicknessModeller")[0]
    simulation = worker.plugCalcAbsorption.getResult()
    worker = recipe.getWorker("MRA Exp")
    minimaPos = worker.plugMra.getResult()[r'\lambda_{min}'].inUnitsOf(simulation.dimensions[1])
    maximaPos = worker.plugMra.getResult()[r'\lambda_{max}'].inUnitsOf(simulation.dimensions[1])
    worker = recipe.getWorker("AddColumn")
    table = worker.plugCompute.getResult(subscriber=TextSubscriber("Add Column"))
    xPos = table[u"x-position"]
    yPos = table[u"y-position"]
    thickness = table[u"thickness"]
    index = curvNo2Index(table[u"pixel"], curvNo)
    result = "$%s_{%s}$(%s %s,%s %s)=%s %s" % (thickness[index].shortname,curvNo,
                                               xPos.data[index],xPos.unit.unit.name(),
                                               yPos.data[index],yPos.unit.unit.name(),
                                               thickness.data[index],
                                               thickness.unit.unit.name())
    pylab.plot(noisyAbsorption.dimensions[1].inUnitsOf(simulation.dimensions[1]).data,
               noisyAbsorption.data[index,:],label="$%s$"%noisyAbsorption.shortname)
    if not noIndicators:
        pylab.vlines(minimaPos.data[:,index],0.1,1.0,
                     label ="$%s$"%minimaPos.shortname)
        pylab.vlines(maximaPos.data[:,index],0.1,1.0,
                     label ="$%s$"%maximaPos.shortname)
    pylab.title(result)
    pylab.xlabel(simulation.dimensions[1].label)
    pylab.ylabel(simulation.label)

 def view_simple( self, stats, thetas ):
   # plotting params
   nbins       = 20
   alpha       = 0.5
   label_size  = 8
   linewidth   = 3
   linecolor   = "r"
   
   # extract from states
   #thetas = states_object.get_thetas()[burnin:,:]
   #stats  = states_object.get_statistics()[burnin:,:]
   #nsims  = states_object.get_sim_calls()[burnin:]
   
   # plot sample distribution of thetas, add vertical line for true theta, theta_star
   f = pp.figure()
   sp = f.add_subplot(111)
   pp.plot( self.fine_theta_range, self.posterior, linecolor+"-", lw = 1)
   ax = pp.axis()
   pp.hist( thetas, self.nbins_coarse, range=self.range,normed = True, alpha = alpha )
   
   pp.fill_between( self.fine_theta_range, self.posterior, color="m", alpha=0.5)
   
   pp.plot( self.posterior_bars_range, self.posterior_bars, 'ro')
   pp.vlines( thetas.mean(), ax[2], ax[3], color="b", linewidths=linewidth)
   #pp.vlines( self.theta_star, ax[2], ax[3], color=linecolor, linewidths=linewidth )
   pp.vlines( self.posterior_mode, ax[2], ax[3], color=linecolor, linewidths=linewidth )
   
   pp.xlabel( "theta" )
   pp.ylabel( "P(theta)" )
   pp.axis([self.range[0],self.range[1],ax[2],ax[3]])
   set_label_fonsize( sp, label_size )
   pp.show()

    def plotcdf(self, x=None, xmin=None, alpha=None, pointcolor='k',
            pointmarker='+', **kwargs):
        """
        Plots CDF and powerlaw
        """
        if x is None: x=self.data
        if xmin is None: xmin=self._xmin
        if alpha is None: alpha=self._alpha

        x=numpy.sort(x)
        n=len(x)
        xcdf = numpy.arange(n,0,-1,dtype='float')/float(n)

        q = x[x>=xmin]
        fcdf = (q/xmin)**(1-alpha)
        nc = xcdf[argmax(x>=xmin)]
        fcdf_norm = nc*fcdf

        D_location = argmax(xcdf[x>=xmin]-fcdf_norm)
        pylab.vlines(q[D_location],xcdf[x>=xmin][D_location],fcdf_norm[D_location],color='m',linewidth=2)

        #plotx = pylab.linspace(q.min(),q.max(),1000)
        #ploty = (plotx/xmin)**(1-alpha) * nc

        pylab.loglog(x,xcdf,marker=pointmarker,color=pointcolor,**kwargs)
        #pylab.loglog(plotx,ploty,'r',**kwargs)
        pylab.loglog(q,fcdf_norm,'r',**kwargs)

def plot_histogram():
    import numpy as np
    import pylab as P

    # The hist() function now has a lot more options
    #

    #
    # first create a single histogram
    #
    P.figure()
    mu, sigma = 40, 35


    x = abs(np.random.normal(mu, sigma, 1000000))

    # the histogram of the data with histtype='step'
    n, bins, patches = P.hist(x, 100, normed=1, histtype='stepfilled')
    P.setp(patches, 'facecolor', 'g', 'alpha', 0.50)
    P.vlines(np.mean(x), 0, max(n))
    P.vlines(np.median(x), 0, max(n))
    # add a line showing the expected distribution
    y = np.abs(P.normpdf( bins, mu, sigma))
    l = P.plot(bins, y, 'k--', linewidth=1.5)

    P.show()

    def OnCalcShift(self, event):
        if (len(self.PSFLocs) > 0):
            import pylab
            
            x,y,z = self.PSFLocs[0]
            
            z_ = numpy.arange(self.image.data.shape[2])*self.image.mdh['voxelsize.z']*1.e3
            z_ -= z_.mean()
            
            pylab.figure()
            p_0 = 1.0*self.image.data[x,y,:,0].squeeze()
            p_0  -= p_0.min()
            p_0 /= p_0.max()

            #print (p_0*z_).sum()/p_0.sum()
            
            p0b = numpy.maximum(p_0 - 0.5, 0)
            z0 = (p0b*z_).sum()/p0b.sum()
            
            p_1 = 1.0*self.image.data[x,y,:,1].squeeze()
            p_1 -= p_1.min()
            p_1 /= p_1.max()
            
            p1b = numpy.maximum(p_1 - 0.5, 0)
            z1 = (p1b*z_).sum()/p1b.sum()
            
            dz = z1 - z0

            print(('z0: %f, z1: %f, dz: %f' % (z0,z1,dz)))
            
            pylab.plot(z_, p_0)
            pylab.plot(z_, p_1)
            pylab.vlines(z0, 0, 1)
            pylab.vlines(z1, 0, 1)
            pylab.figtext(.7,.7, 'dz = %3.2f' % dz)

def plot_x_and_psd_with_estimated_omega(x, sample_step=1, dt=1.0):
    y = x[::sample_step]
    F = freq(x, sample_step, dt)
    T = 1.0 / F

    pylab.clf()

    # plot PSD
    pylab.subplot(211)
    pylab.psd(y, Fs=1.0 / sample_step / dt)
    ylim = pylab.ylim()
    pylab.vlines(F, *ylim)
    pylab.ylim(ylim)

    # plot time series
    pylab.subplot(223)
    pylab.plot(x)

    # plot time series (three cycles)
    n = int(T / dt) * 3
    m = n // sample_step
    pylab.subplot(224)
    pylab.plot(x[:n])
    pylab.plot(numpy.arange(0, n, sample_step)[:m], y[:m], 'o')

    pylab.suptitle('F = %s' % F)

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 plot_matrices(cov, prec, title, subject_n=0):
    """Plot covariance and precision matrices, for a given processing. """

    # Put zeros on the diagonal, for graph clarity.
    size = prec.shape[0]
    prec[range(size), range(size)] = 0

    span = max(abs(prec.min()), abs(prec.max()))
    title = "{0:d} {1}".format(subject_n, title)

    # Display covariance matrix
    pl.figure()
    pl.imshow(cov, interpolation="nearest",
              vmin=-1, vmax=1, cmap=pl.cm.get_cmap("bwr"))
    pl.hlines([(pl.ylim()[0] + pl.ylim()[1]) / 2],
              pl.xlim()[0], pl.xlim()[1])
    pl.vlines([(pl.xlim()[0] + pl.xlim()[1]) / 2],
              pl.ylim()[0], pl.ylim()[1])
    pl.colorbar()
    pl.title(title + " / covariance")

    # Display precision matrix
    pl.figure()
    pl.imshow(prec, interpolation="nearest",
              vmin=-span, vmax=span,
              cmap=pl.cm.get_cmap("bwr"))
    pl.hlines([(pl.ylim()[0] + pl.ylim()[1]) / 2],
              pl.xlim()[0], pl.xlim()[1])
    pl.vlines([(pl.xlim()[0] + pl.xlim()[1]) / 2],
              pl.ylim()[0], pl.ylim()[1])
    pl.colorbar()
    pl.title(title + " / precision")

    def update_loop(self):
        import pylab
        import numpy
        pylab.ion()
        self.fig = pylab.figure()
        vmin, vmax = -1, 1
        self.img = pylab.imshow(self.retina.image, vmin=vmin, vmax=vmax,
                               cmap='gray', interpolation='none')
        self.line_y = pylab.vlines(self.retina.p_y, 0, 128)
        self.line_x = pylab.vlines(self.retina.p_x, 0, 128)
        while True:
            #self.fig.clear()
            #pylab.imshow(self.retina.image, vmin=vmin, vmax=vmax,
            #                   cmap='gray', interpolation='none')
            self.img.set_data(self.retina.image)
            #self.scatter[0].set_data([self.retina.p_x], [self.retina.p_y])
            #if len(self.retina.deltas) > 0:
            #    pylab.hist(self.retina.deltas, 20, range=(0.0, 0.1))
            #pylab.ylim(0, 100)
            #pylab.vlines(self.retina.p_y, 0, 128)
            #pylab.hlines(self.retina.p_x, 0, 128)

            y = self.retina.p_y
            self.line_y.set_segments(numpy.array([[(y,0), (y,128)]]))
            x = self.retina.p_x
            self.line_x.set_segments(numpy.array([[(0,x), (128,x)]]))

            self.retina.clear_image()
            pylab.draw()
            time.sleep(0.04)

def plot_hist_rmsRho_Levs012(rho_L0, rho_L1, rho_L2, tSim, fname):
    ''''''
    import pylab as py
    #make 6 hists
    rms = rho_L0
    lab = r"$\rho$, Lev 0"
    rms1d = np.reshape(rms, np.product(rms.shape))
    logrms= np.log10(rms1d)
    cnt,bins,patches = py.hist(logrms, bins=100, color="blue", alpha=0.5, label=lab)
    rms = rho_L1
    lab = r"$\rho$, Lev 1"
    rms1d = np.reshape(rms, np.product(rms.shape))
    logrms= np.log10(rms1d)
    cnt,bins,patches = py.hist(logrms, bins=100, color="red", alpha=0.5, label=lab)
    rms = rho_L2
    lab = r"$\rho$, Lev 2"
    rms1d = np.reshape(rms, np.product(rms.shape))
    logrms= np.log10(rms1d)
    #plot quantities
    Tratio = tSim / ic.tCr
    #plot
    cnt,bins,patches = py.hist(logrms, bins=100, color="green", alpha=0.5,label=lab)
    py.vlines(np.log10( ic.rho0 ), 0, cnt.max(), colors="black", linestyles='dashed',label=r"$\rho_{0}$ = %2.2g [g/cm^3]" % ic.rho0)
    #py.xlim([-13.,-9.])
    py.xlabel("Log10 Density [g/cm^3]")
    py.ylabel("count")
    py.title(r"$T/T_{\rmCross}$ = %g" % Tratio)
    py.legend(loc=0, fontsize="small")
    py.savefig(fname,format="pdf")
    py.close()

def oneRoundDelay(step = 10000):
    injectRate = 0.9
    factory = GridNetworkFactory(makeSimpleNode(),Queues)
    factory.constructNetwork(6,6)\
        .setFlow(Flow((0,0),(0,5),injectRate))\
        .setFlow(Flow((5,0),(5,5),injectRate))\
        .setFlow(Flow((2,0),(3,5),injectRate))

    network = factory.getNetwork()

    packetFactory = PacketFactory()
    simulator = \
        Simulator(network,step,ConstLinkRateGenerator(1),packetFactory)
    simulator.run()
    #simulator.printNetwork()
    stat = simulator.getStaticsInfo()
    print stat['aveDelay']
    packetPool = sorted(stat['packetPool'],key=lambda p: p.getID)


    
    py.subplot(211)
    py.vlines([p.getCreateTime() for p in packetPool],[1],[p.getDelay() for p in packetPool],'r')
    py.xlabel('Packet create time(bp with $\lambda$ = 0.9)')
    py.ylabel('delay')
    py.grid(True)




    injectRate = 0.9
    factory = GridNetworkFactory(makeMNode(2),Queues)
    factory.constructNetwork(6,6)\
        .setFlow(Flow((0,0),(0,5),injectRate))\
        .setFlow(Flow((5,0),(5,5),injectRate))\
        .setFlow(Flow((2,0),(3,5),injectRate))

    network = factory.getNetwork()

    packetFactory = PacketFactory()
    simulator = \
        Simulator(network,step,ConstLinkRateGenerator(1),packetFactory)
    simulator.run()
    #simulator.printNetwork()

    stat = simulator.getStaticsInfo()
    print stat['aveDelay']
    packetPool = sorted(stat['packetPool'],key=lambda p: p.getID)

    py.subplot(212)
    py.vlines([p.getCreateTime() for p in packetPool],[1],[p.getDelay() for p in packetPool],'b')
    py.xlabel('Packet create time (m=2 with $\lambda$ = 0.9)')
    py.ylabel('delay')
    py.grid(True)
    py.savefig('packetDelayInOneRound_09')
    py.show()

def myPlot(X,fun,X0,funp,c1,c2):
    p.figure(figsize=(10,6))
    F=do(f,X)
    p.plot(X0,funp,c1+'o-',lw=2,)
    p.vlines(X0,0,funp,linestyle='dotted')
    p.hlines(0,0,1)
    p.xticks(X0,['$x_{%s}$'%i for i in range(len(X0))],fontsize='large')
    p.yticks([])
    p.ylabel('$f(x)$',fontsize='large')
    p.axis([X0[0],X0[-1],0,F.max()+0.1])

def periodogram_scipy(data, frequencies, plot=False):
    freqs, pgram = lombscargle_scipy(data, frequencies)
    periods = map(lambda x: (2 * math.pi) / x, freqs)

    if plot:
        pylab.vlines(periods, 0, np.array(pgram), color='k', linestyles='solid')
        pylab.xlabel("period, days")
        pylab.ylabel("power")
        pylab.title("CenX-3 18-60 KeV Periodogram")
        pylab.show()

    def plotpdf(self,x=None,xmin=None,alpha=None,nbins=50,dolog=True,dnds=False,
            drawstyle='steps-post', histcolor='k', plcolor='r', **kwargs):
        """
        Plots PDF and powerlaw.

        kwargs is passed to pylab.hist and pylab.plot
        """
        if not(x): x=self.data
        if not(xmin): xmin=self._xmin
        if not(alpha): alpha=self._alpha

        x=numpy.sort(x)
        n=len(x)

        pylab.gca().set_xscale('log')
        pylab.gca().set_yscale('log')

        if dnds:
            hb = pylab.histogram(x,bins=numpy.logspace(log10(min(x)),log10(max(x)),nbins))
            h = hb[0]
            b = hb[1]
            db = hb[1][1:]-hb[1][:-1]
            h = h/db
            pylab.plot(b[:-1],h,drawstyle=drawstyle,color=histcolor,**kwargs)
            #alpha -= 1
        elif dolog:
            hb = pylab.hist(x,bins=numpy.logspace(log10(min(x)),log10(max(x)),nbins),log=True,fill=False,edgecolor=histcolor,**kwargs)
            alpha -= 1
            h,b=hb[0],hb[1]
        else:
            hb = pylab.hist(x,bins=numpy.linspace((min(x)),(max(x)),nbins),fill=False,edgecolor=histcolor,**kwargs)
            h,b=hb[0],hb[1]
        # plotting points are at the center of each bin
        b = (b[1:]+b[:-1])/2.0

        q = x[x>=xmin]
        px = (alpha-1)/xmin * (q/xmin)**(-alpha)

        # Normalize by the median ratio between the histogram and the power-law
        # The normalization is semi-arbitrary; an average is probably just as valid
        plotloc = (b>xmin)*(h>0)
        norm = numpy.median( h[plotloc] / ((alpha-1)/xmin * (b[plotloc]/xmin)**(-alpha))  )
        px = px*norm

        plotx = pylab.linspace(q.min(),q.max(),1000)
        ploty = (alpha-1)/xmin * (plotx/xmin)**(-alpha) * norm

        #pylab.loglog(q,px,'r',**kwargs)
        pylab.loglog(plotx,ploty,color=plcolor,**kwargs)

        axlims = pylab.axis()
        pylab.vlines(xmin,axlims[2],max(px),colors=plcolor,linestyle='dashed')

        pylab.gca().set_xlim(min(x),max(x))

 def view_results( self, states_object, burnin = 1 ):
   # plotting params
   nbins       = 20
   alpha       = 0.5
   label_size  = 8
   linewidth   = 3
   linecolor   = "r"
   
   # extract from states
   thetas = states_object.get_thetas()[burnin:,:]
   stats  = states_object.get_statistics()[burnin:,:]
   nsims  = states_object.get_sim_calls()[burnin:]
   
   f=pp.figure()
   for i in range(6):
     sp=f.add_subplot(2,10,i+1)
     pp.hist( thetas[:,i], 10, normed=True, alpha = 0.5)
     pp.title( self.theta_names[i])
     set_label_fonsize( sp, 6 )
     set_tick_fonsize( sp, 6 )
     set_title_fonsize( sp, 8 )
   for i in range(10):
     sp=f.add_subplot(2,10,10+i+1)
     pp.hist( stats[:,i], 10, normed=True, alpha = 0.5)
     ax=pp.axis()
     pp.vlines( self.obs_statistics[i], 0, ax[3], color="r", linewidths=2)
     # if self.obs_statistics[i] < ax[0]:
     #   ax[0] = self.obs_statistics[i]
     # elif self.obs_statistics[i] > ax[1]:
     #   ax[1] = self.obs_statistics[i]
     pp.axis( [ min(ax[0],self.obs_statistics[i]), max(ax[1],self.obs_statistics[i]), ax[2],ax[3]] )
     pp.title( self.stats_names[i])
     set_label_fonsize( sp, 6 )
     set_tick_fonsize( sp, 6 )
     set_title_fonsize( sp, 8 )
   pp.suptitle( "top: posterior, bottom: post pred with true")
   
   f = pp.figure()  
   I = np.random.permutation( len(thetas) )
   for i in range(16):
     sp=pp.subplot(4,4,i+1)
     theta = thetas[ I[i],:]
     test_obs = self.simulation_function( theta )
     test_stats = self.statistics_function( test_obs )
     err = np.sum( np.abs( self.obs_statistics - test_stats ) )
     pp.title( "%0.2f"%( err ))
     pp.plot( self.observations/1000.0 )
     pp.plot(test_obs/1000.0)
     pp.axis("off")
     set_label_fonsize( sp, 6 )
     set_tick_fonsize( sp, 6 )
     set_title_fonsize( sp, 8 )
   pp.suptitle( "time-series from random draws of posterior")