def random_test():
    rand_periods = np.zeros(1000)
    periods = np.zeros(1000)
    for i in range(1000):
        rand_periods[i] = np.random.uniform(low = 0.0, high = 2.0)
        periods[i] = 10**rand_periods[i]

    true_periods = np.zeros(1000)
    for i in range(1000):
        data = np.genfromtxt('/Users/angusr/angusr/ACF/star_spot_sim/tests/sim_period%s.txt' %(i+1))
        true_periods[i] = data

    p.close(4)
    p.figure(4)
    p.subplot(3,1,1)
    p.plot(rand_periods, 'k.')
    p.subplot(3,1,2)
    p.plot(periods, 'k.')
    p.subplot(3,1,3)
    p.plot(np.log10(true_periods) ,'k.')

    ''' Plotting as close to original periods as I can'''
    p.close(10)
    p.figure(10)
    p.subplot(1,2,1)
    orig_periods = np.zeros(100)
    for i in range(100):
        data = np.genfromtxt('/Users/angusr/angusr/ACF/star_spot_sim/sim_period%s.txt' %(i+1)).T
        p.axhline(np.log10(data[4]), color = 'k')
    p.subplot(1,2,2)
    for i in range(100):
        data = np.genfromtxt('/Users/angusr/angusr/ACF/star_spot_sim/grid/%sparams.txt' %(i+1)).T
        p.axhline(np.log10(data[7]), color = 'k')

def test():
    from pandas import DataFrame
    X = np.linspace(0.01, 1.0, 10)
    Y = np.log(X)
    Y -= Y.min()
    Y /= Y.max()
    Y *= 0.95

    #Y = X

    df = DataFrame({'X': X, 'Y': Y})
    P = Pareto(df, 'X', 'Y')

    data = []
    for val in np.linspace(0,1,15):
        data.append(dict(val=val, x=P.lookup_x(val), y=P.lookup_y(val)))
        pl.axvline(val, alpha=.5)
        pl.axhline(val, alpha=.5)
    dd = DataFrame(data)
    pl.scatter(dd.y, dd.val, lw=0, c='r')
    pl.scatter(dd.val, dd.x, lw=0, c='g')
    print dd

    #P.scatter(c='r', lw=0)
    P.show_frontier(c='r', lw=4)
    pl.show()

 def imshow_box(f,im, x,y,s):
     '''imshow_box(f,im, x,y,s)
     f: figure
     im: image
     x: center coordinate for box
     y: center coord
     s: box shape, (width, height)
     '''
     global coord
     P.figure(f.number)
     P.clf();
     P.imshow(im);
     P.axhline(y-s[1]/2.)
     P.axhline(y+s[1]/2.)
     P.axvline(x-s[0]/2.)
     P.axvline(x+s[0]/2.)
     xy=crop(m,s,y,x)
     coord=(0.5*(xy[2]+xy[3]), 0.5*(xy[0]+xy[1]))
     P.title(str('x: %d y: %d' % (x,y)));        
     P.figure(999);
     P.imshow(master[xy[0]:xy[1],xy[2]:xy[3]])
     P.title('Master');
     P.figure(998);
     df=(master[xy[0]:xy[1],xy[2]:xy[3]]-slave)
     P.imshow(np.abs(df))
     P.title(str('RMS: %0.6f' % np.sqrt((df**2.).mean()) ));        

def show_table(table_name,ls="none", fmt="o", legend=False, name="m", do_half=0):
	bt = fi.FITS(table_name)[1].read()
	rgpp = (np.unique(bt["rgp_lower"])+np.unique(bt["rgp_upper"]))/2
	nbins = rgpp.size

	plt.xscale("log")
	colours=["purple", "forestgreen", "steelblue", "pink", "darkred", "midnightblue", "gray", "sienna", "olive", "darkviolet"]
	pts = ["o", "D", "x", "^", ">", "<", "1", "s", "*", "+", "."]
	for i,r in enumerate(rgpp):
		sel = (bt["i"]==i)
		snr = 10** ((np.log10(bt["snr_lower"][sel]) + np.log10(bt["snr_upper"][sel]))/2)

		if do_half==1 and i>nbins/2:
			continue
		elif do_half==2 and i<nbins/2:
			continue
		if legend:
			plt.errorbar(snr, bt["%s"%name][i*snr.size:(i*snr.size)+snr.size], bt["err_%s"%name][i*snr.size:(i*snr.size)+snr.size], color=colours[i], ls=ls, fmt=pts[i], lw=2.5, label="$R_{gpp}/R_p = %1.2f-%1.2f$"%(np.unique(bt["rgp_lower"])[i],np.unique(bt["rgp_upper"])[i]))
		else:
			plt.errorbar(snr, bt["%s"%name][i*snr.size:(i*snr.size)+snr.size], bt["err_%s"%name][i*snr.size:(i*snr.size)+snr.size], color=colours[i], ls=ls, fmt=pts[i], lw=2.5)

	plt.xlim(10,300)
	plt.axhline(0, lw=2, color="k")
	
	plt.xlabel("Signal-to-Noise $SNR_w$")
	if name=="m":
		plt.ylim(-0.85,0.05)
		plt.ylabel("Multiplicative Bias $m \equiv (m_1 + m_2)/2$")
	elif name=="alpha":
		plt.ylabel(r"PSF Leakage $\alpha \equiv (\alpha _1 + \alpha _2)/2$")
		plt.ylim(-0.5,2)



	plt.legend(loc="lower right")

def _show_rates(rate, wo, wt, attenuator, tau_NP, tau_P):
    import pylab

    #pylab.figure()
    pylab.errorbar(rate, wt[0], yerr=wt[1], fmt='g.', label='attenuated')
    pylab.errorbar(rate, wo[0], yerr=wo[1], fmt='b.', label='unattenuated')

    pylab.xscale('log')
    pylab.yscale('log')
    pylab.xlabel('incident rate (counts/second)')
    pylab.ylabel('observed rate (counts/second)')
    pylab.legend(loc='best')
    pylab.grid(True)
    pylab.plot(rate, rate/attenuator, 'g-', label='target')
    pylab.plot(rate, rate, 'b-', label='target')

    Ipeak, Rpeak = peak_rate(tau_NP=tau_NP, tau_P=tau_P)
    if rate[0] <= Ipeak <= rate[-1]:
        pylab.axvline(x=Ipeak, ls='--', c='b')
        pylab.text(x=Ipeak, y=0.05, s=' %g'%Ipeak,
                   ha='left', va='bottom',
                   transform=pylab.gca().get_xaxis_transform())
    if False:
        pylab.axhline(y=Rpeak, ls='--', c='b')
        pylab.text(y=Rpeak, x=0.05, s=' %g\n'%Rpeak,
                   ha='left', va='bottom',
                   transform=pylab.gca().get_yaxis_transform())

def plot_resid(x, resid):
    import pylab
    pylab.plot(x, resid, '.')
    pylab.gca().locator_params(axis='y', tight=True, nbins=4)
    pylab.axhline(y=1, ls='dotted')
    pylab.axhline(y=-1, ls='dotted')
    pylab.ylabel("Residuals")

def gfe4():
  x2=plt.linspace(1e-20,.13,90000)
  xmin2=((4*np.pi*(SW.R)**3)/3)*1e-20
  xmax2=((4*np.pi*(SW.R)**3)/3)*.13
  xff2 = plt.linspace(xmin2,xmax2,90000)
  thigh=100
  plt.figure()
  plt.title('Grand free energy per volume vs ff @ T=%0.4f'%Tlist[thigh])
  plt.ylabel('Grand free energy per volume')
  plt.xlabel('filling fraction')  
  plt.plot(xff2,SW.phi(Tlist[thigh],x2,nR[thigh]),color='#f36118',linewidth=3)
  #plt.axvline(nL[thigh])
  #plt.axvline(nR[thigh])
  #plt.axhline(SW.phi(Tlist[thigh],nR[thigh]))
  #plt.plot(x2,x2-x2,'c')
  plt.plot(nL[thigh]*((4*np.pi*(SW.R)**3)/3),SW.phi(Tlist[thigh],nL[thigh],nR[thigh]),'ko')
  plt.plot(nR[thigh]*((4*np.pi*(SW.R)**3)/3),SW.phi(Tlist[thigh],nR[thigh],nR[thigh]),'ko')
  plt.axhline(SW.phi(Tlist[thigh],nR[thigh],nR[thigh]),color='c',linewidth=2)
  print(Tlist[100])
  print(nL[100],nR[100])
  plt.savefig('figs/gfe_cotangent.pdf')

  plt.figure()
  plt.plot(xff2,SW.phi(Tlist[thigh],x2,nR[thigh]),color='#f36118',linewidth=3)
  plt.plot(nL[thigh]*((4*np.pi*(SW.R)**3)/3),SW.phi(Tlist[thigh],nL[thigh],nR[thigh]),'ko')
  plt.plot(nR[thigh]*((4*np.pi*(SW.R)**3)/3),SW.phi(Tlist[thigh],nR[thigh],nR[thigh]),'ko')
  plt.axhline(SW.phi(Tlist[thigh],nR[thigh],nR[thigh]),color='c',linewidth=2)
  plt.xlim(0,0.0003)
  plt.ylim(-.000014,0.000006)
  print(Tlist[100])
  print(nL[100],nR[100])
  plt.savefig('figs/gfe_insert_cotangent.pdf')

def findSeriesLength(teamProb):
    numSeries = 200
    maxLen = 2500
    step = 10

    def fracWon(teamProb, numSeries, seriesLen):
        won = 0.0
        for series in range(numSeries):
            if playSeries(seriesLen, teamProb):
                won += 1
        return won/numSeries

    winFrac = []
    xVals = []
    for seriesLen in range(1, maxLen, step):
        xVals.append(seriesLen)
        winFrac.append(fracWon(teamProb, numSeries, seriesLen))

    pylab.plot(xVals, winFrac, linewidth=5)
    pylab.xlabel("Length of Series")
    pylab.ylabel("Probability of winning a series")
    pylab.title(str(round(teamProb, 4)) + ' probability of team winning a game')
    # draw horizontal line at 0.95
    pylab.axhline(0.95)
    pylab.show()

def doplot(data, title, ylim=None, yaxis = 'Quantity', meanval=False,
                                                      color='b',
                                                      label='quantity',
                                                      showfig = False):

    pl.title(title, fontsize = 15, color='k')
    fig = pl.plot(range(len(data)), data, color, label=label, linewidth=1.5)
    pl.ylabel(yaxis)
    #pl.xticks( np.arange(0, len(data)+len(data)*0.1, len(data)*0.1 ) )
    x_mean = np.mean(data)
    if meanval==True:
        pl.axhline( x_mean, 0, len(data), color='r', linewidth=1.3, label='mean')
    if ylim<>None:
        pl.ylim(ylim)
    pl.xlabel('Number of maps')
    pl.xticks(np.arange(0,len(data),1))
    pl.yticks(np.arange(0,1.1,0.1))
    pl.grid(False)
    pl.gca().yaxis.grid(True)
    pl.legend(loc='upper left', numpoints = 1)

    if showfig:
        pl.show()


    return fig

def plot_tm_rbf_decision_values(class_ids, dec_values, plot_title = '', plot_file = ''):
    
    import pylab as pl
    from collections import defaultdict
    
    print 
    true_class = defaultdict(list)
    for i, class_id in enumerate(class_ids):
        print '#%d true class: %d decision value: %.5f' % (i, class_id, dec_values[i])
        true_class[class_id] += [dec_values[i]]
    print 
    
    pl.clf()
    pl.plot(true_class[IRRELEVANT_CLASS_ID], 'bo', label='Irrelevant')
    x2 = range(len(true_class[IRRELEVANT_CLASS_ID]), len(class_ids))
    pl.plot(x2, true_class[RELEVANT_CLASS_ID], 'r+', label='Relevant')
    pl.axhline(0, color='black')
    pl.xlabel('Documents')
    pl.ylabel('Decision values')
    pl.title(plot_title)
    pl.legend(loc='lower right', prop={'size':9})
    pl.grid(True)

    if (plot_file == ''):
        pl.show()
    else: 
        pl.savefig(plot_file, dpi=300, bbox_inches='tight', pad_inches=0.1)
    pl.close()
    pl.clf()

def plotramp(sci, err, tsamp, nsamp, plotfit=False, **kwargs):
    """ plot the up-the-ramp sampling sequence 
    for the given pixel, given in _ima coordinates.
    kwargs are passed to pylab.errorbar().
    """
    from pylab import plot, errorbar, legend, xlabel, ylabel, axes, axhline

    # sci, err, tsamp, nsamp  = getrampdat(imafile, x, y )
    if plotfit:
        ax1 = axes([0.1, 0.35, 0.85, 0.6])
        m, b = fitramp(sci, err, tsamp, nsamp)
        fit = m * nsamp + b
        plot(nsamp, fit, ls="--", color="k", marker="")
        errorbar(nsamp, sci * tsamp, err * tsamp, **kwargs)
        ylabel("cumulative counts (sci*tsamp)")
        legend(loc="upper left")

        ax2 = axes([0.1, 0.1, 0.85, 0.25], sharex=ax1)
        kwargs["ls"] = " "
        errorbar(nsamp[1:], sci[1:] * tsamp[1:] - fit[1:], err[1:] * tsamp[1:], **kwargs)
        axhline(ls="--", color="k")

        xlabel("SAMP NUMBER")
    else:
        errorbar(nsamp, sci * tsamp, err * tsamp, **kwargs)
        xlabel("SAMP NUMBER")
        ylabel("cumulative counts (sci*tsamp)")
        legend(loc="upper left")

def PlotSlice(LogLikelihood,par,low,upp,par_in,func_args=(),plot_samp=100):
  """
  Plot the conditional distributions for each variable parameter. Used to visualise the
  conditional errors, and get sensible inputs to ConditionalErrors function.
  
  """
  
  i = par_in
  op_par = np.copy(par)
  max_loglik = LogLikelihood(op_par,*func_args)  
  
  par_range = np.linspace(low,upp,plot_samp)
  log_lik = np.zeros(plot_samp)
  temp_par = np.copy(op_par)
  for q,par_val in enumerate(par_range):
    temp_par[i] = par_val
    log_lik[q] = LogLikelihood(temp_par,*func_args)
  print np.exp(log_lik-max_loglik)
  pylab.clf()
  pylab.subplot(211)
  pylab.plot(par_range,log_lik)
  pylab.axhline(max_loglik-0.5,color='g',ls='--')
  pylab.xlabel("p[%s]" % str(i))
  pylab.ylabel("log Posterior")
  pylab.subplot(212)
  pylab.plot(par_range,np.exp(log_lik-max_loglik))
  pylab.axvline(op_par[i],color='r')
  pylab.axhline(0.6065,color='g',ls='--')
  pylab.xlabel("p[%s]" % str(i))
  pylab.ylabel("Posterior")

def PlotConditionals(LogLikelihood,par,err,low,upp,func_args=(),plot_samp=100,opt=False,par_in=None,wait=False):
  """
  Plot the conditional distributions for each variable parameter. Used to visualise the
  conditional errors, and get sensible inputs to ConditionalErrors function.
  
  """
  
  #first optimise the log likelihood?
  if opt: op_par = Optimise(LogLikelihood,par[:],func_args,fixed=(np.array(err) == 0)*1)
  else: op_par = np.copy(par)
  
  max_loglik = LogLikelihood(op_par,*func_args)
  
  if par_in == None: par_in = np.where(np.array(err) != 0.)[0]
  
  for i in par_in:
   
   par_range = np.linspace(low[i],upp[i],plot_samp)
   log_lik = np.zeros(plot_samp)
   temp_par = np.copy(op_par)
   for q,par_val in enumerate(par_range):
     temp_par[i] = par_val
     log_lik[q] = LogLikelihood(temp_par,*func_args)
   pylab.clf()
   pylab.plot(par_range,log_lik)
   pylab.plot(par_range,max_loglik-(par_range-op_par[i])**2/2./err[i]**2,'r--')
   pylab.axvline(op_par[i],color='r')
   pylab.axvline(op_par[i]+err[i],color='g')
   pylab.axvline(op_par[i]-err[i],color='g')
   pylab.axhline(max_loglik-0.5,color='g',ls='--')
   pylab.xlabel("p[%s]" % str(i))
   pylab.ylabel("log Posterior")
   #pylab.xlims(low[i],upp[i])
   if wait: raw_input("")  

def plot_one_PSTH(arr, n_trial, ax, rng=DEF_PSTH_RANGE_MS, xrng=DEF_PSTH_RANGE_MS, \
        aggregate=10, ymax=250, color='#999999', nticks=5, visible=False, txt=None, nondraw=False):
    arr = np.array(arr)
    # plot PSTH
    if n_trial > 0:
        arr = arr/1000.                       # to ms
        interval = rng[1] - rng[0]
        bins = int(interval / aggregate)
        weight = 1000. / aggregate / n_trial  # to convert into Spiks/s.
        weights = np.array([weight] * len(arr))
        if nondraw:
            rtn = np.histogram(arr, range=rng, bins=bins, weights=weights)
            pl.axhline(y=0, color='#333333')
        else:
            rtn = pl.hist(arr, range=rng, bins=bins, weights=weights, fc=color, ec=color)
    # beutify axes
    pl.xlim(xrng)
    pl.ylim([0,ymax])
    ax.xaxis.set_major_locator(pl.MaxNLocator(nticks))
    ax.yaxis.set_major_locator(pl.MaxNLocator(nticks))
    if not visible:
        ax.set_xticklabels([''] * nticks)
        ax.set_yticklabels([''] * nticks)
    pl.axvline(x=0, ymin=0, ymax=1, lw=0.5, color='r')
    if txt != None:
        if nondraw:
            pl.text(xrng[1] - 20, ymax - 70, txt, size=6, va='top', ha='right')
        else:
            pl.text(xrng[1] - 20, ymax - 20, txt, size=6, va='top', ha='right')
    return rtn

def Cross(x0=0.0, y0=0.0, clr='black', ls='dashed', lw=1, zorder=0):
    """
    Draw cross through zero
    =======================
    """
    axvline(x0, color=clr, linestyle=ls, linewidth=lw, zorder=zorder)
    axhline(y0, color=clr, linestyle=ls, linewidth=lw, zorder=zorder)

def video(data, Ti=None, Tf=None):
    fig, ax = py.subplots()
    py.xlim([-0.5, 9.5])
    py.ylim([-0.5, 9.5])
    py.axvline(xc[1] - 0.5, color="black", linestyle="--")
    py.axhline(yc[2] - 0.5, color="black", linestyle="--")

    if Tf == None:
        Tf = data.shape[2] - 1
    if Ti == None:
        Ti = 0

    lines = []
    for i in range(num_walkers):
        x = data[i, 0, Ti : Ti + 1]
        y = data[i, 1, Tf : Tf + 1]
        line, = py.plot(x, y, "o")
        lines.append(line)

    def animate(i):
        py.title(i)
        for j, line in enumerate(lines):
            line.set_xdata(data[j, 0, i : i + 1])
            line.set_ydata(data[j, 1, i : i + 1])
        return lines

    ani = animation.FuncAnimation(fig, animate, np.arange(Ti, Tf), interval=300, blit=False)

    py.show()

def test():
    if 0:
        from pandas import DataFrame
        X = np.linspace(0.01, 1.0, 10)
        Y = np.log(X)
        Y -= Y.min()
        Y /= Y.max()
        Y *= 0.95

        df = DataFrame({'X': X, 'Y': Y})
        P = Pareto(df, 'X', 'Y')

        data = []
        for val in np.linspace(0,1,15):
            data.append(dict(val=val, x=P.lookup_x(val), y=P.lookup_y(val)))
            pl.axvline(val, alpha=.5)
            pl.axhline(val, alpha=.5)
        dd = DataFrame(data)
        pl.scatter(dd.y, dd.val, lw=0, c='r')
        pl.scatter(dd.val, dd.x, lw=0, c='g')
        print dd

        #P.scatter(c='r', lw=0)
        P.show_frontier(c='r', lw=4)
        pl.show()

    X,Y = np.random.normal(0,1,size=(2, 30))

    for maxX in [0,1]:
        for maxY in [0,1]:
            pl.figure()
            pl.title('max x: %s, max y: %s' % (maxX, maxY))
            pl.scatter(X,Y,lw=0)
            show_frontier(X, Y, maxX=maxX, maxY=maxY)
            pl.show()

 def plot_complex (content,title):
   """Plots x vs y""";
   # plot errors bars, if available
   pylab.axhline(0,color='lightgrey')
   pylab.axvline(1,color='lightgrey')
   pylab.errorbar(
     [x.real for l1,l2,(x,xe),(y,ye) in content],[x.imag for l1,l2,(x,xe),(y,ye) in content],
     [xe for l1,l2,(x,xe),(y,ye) in content],[xe for l1,l2,(x,xe),(y,ye) in content],
     fmt=None,ecolor="lightgrey"
   );
   pylab.errorbar(
     [y.real for l1,l2,(x,xe),(y,ye) in content],[y.imag for l1,l2,(x,xe),(y,ye) in content],
     [ye for l1,l2,(x,xe),(y,ye) in content],[ye for l1,l2,(x,xe),(y,ye) in content],
     fmt=None,ecolor="lightgrey"
   );
   # max plot amplitude -- max point plus 1/4th of the error bar
   maxa = max([ max(abs(x),abs(y)) for l1,l2,(x,xe),(y,ye) in content ]);
   # plotlim = max([ abs(numpy.array([ 
   #                  getattr(v,attr)+sign*e/4 for v,e in (x,xe),(y,ye) for attr in 'real','imag' for sign in 1,-1 
   #                ])).max() 
   #   for l1,l2,(x,xe),(y,ye) in content ]);
   minre, maxre, minim, maxim = 2, -2, 2, -2
   for l1,l2,(x,xe),(y,ye) in content:
       offs = numpy.array([ getattr(v,attr)+sign*e/4 for v,e in (x,xe),(y,ye) 
                 for attr in 'real','imag' for sign in 1,-1 ])
       minre, maxre = min(x.real-xe/4, y.real-ye/4, minre), max(x.real+xe/4, y.real+ye/4, maxre)

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 OffsetPlot():
    import pylab as P
    import scipy.stats as S
    
    offsp = S.spearmanr(data[:,dict['medianOffset']], telfocusCorrected)
    offreg = S.linregress(data[:,dict['medianOffset']], telfocusCorrected)
    offreg2 = S.linregress(data[:,dict['medianOffset']], telfocusOld)
    min = -50.
    max = 50.
    
    print '\nOffset Spearman rank-order:', offsp
    print 'Offset fit:', offreg
    print 'and For unCorrected data:', offreg2
    
    P.plot(data[:,dict['medianOffset']], telfocusCorrected, 'bo', label = 'Data')
    P.plot([min,max], [min*offreg[0] + offreg[1], max*offreg[0] + offreg[1]], 
           'r-', label ='Linear Fit (Corrected)', lw = 2.0)
    P.plot([min,max], [min*offreg2[0] + offreg2[1], max*offreg2[0] + offreg2[1]], 
           'g--', label ='Linear Fit (UnCorrected)', lw = 1.5)
    P.axhline(medianNew, color ='b')
    P.xlim(min, max)
    P.xlabel('Median Offset (telescope units)')
    P.ylabel('Temperature Corrected Telescope Focus + Median Offset')
    P.legend(shadow=True)
    P.savefig('offsetCorrelation.png')
    P.close()