def plot_thresholds(rawdata, scan_values, plane='horizontal',
                    xlabel='turns', ylabel='intensity [particles]', zlabel='normalized emittance',
                    xlimits=((0.,8192)), ylimits=((0.,7.1e11)), zlimits=((0., 10.))):

    # Prepare input data.
    # x axis
    t = rawdata[0,:,:]
    turns = plt.ones(t.shape).T * plt.arange(len(t))
    turns = turns.T

    # z axis
    epsn_abs = {}
    epsn_abs['horizontal'] = plt.absolute(rawdata[11,:,:])
    epsn_abs['vertical']   = plt.absolute(rawdata[12,:,:])

    # Prepare plot environment.
    ax11, ax13 = _create_axes(xlabel, ylabel, zlabel, xlimits, ylimits, zlimits)
    cmap = plt.cm.get_cmap('jet', 2)
    ax11.patch.set_facecolor(cmap(range(2))[-1])
    cmap = plt.cm.get_cmap('jet')

    x, y = plt.meshgrid(turns[:,0], scan_values)
    z = epsn_abs[plane]

    threshold_plot = ax11.contourf(x, y, z.T, levels=plt.linspace(zlimits[0], zlimits[1], 201),
                                   vmin=zlimits[0], vmax=zlimits[1], cmap=cmap)
    cb = plt.colorbar(threshold_plot, ax13, orientation='vertical')
    cb.set_label(zlabel)

    plt.tight_layout()

def add_to_results(model, name):
    df = getattr(model, name)
    model.results['param'].append(name)
    model.results['bias'].append(df['abs_err'].mean())
    model.results['mae'].append((pl.median(pl.absolute(df['abs_err'].dropna()))))
    model.results['mare'].append(pl.median(pl.absolute(df['rel_err'].dropna())))
    model.results['pc'].append(df['covered?'].mean())

def scalar_validation_statistics(results, groups):
    """ plot absolute z transformation of p_theta values,
    grouped by dictionary groups

    Parameters
    ----------
    results : pandas.DataFrame with row called 'z'
    groups : list of lists of columns of results
    """
    pl.figure()

    width = max(pl.absolute(results.ix['z']))
    for row, (g_name, g) in enumerate(reversed(groups)):
        z = pl.absolute(results.ix['z', g].__array__())
        pl.plot([pl.mean(z)], [row], 'o', color='k', mec='k', mew=1)
        pl.plot(z, [row]*len(z), 'o', color='none', mec='k', mew=1)

        msg = 'p: %s' % ', '.join(['%.3f'%p for p in sorted(results.ix['p', g]*len(g))])
        #msg += 'MAE: %s' % str(
        pl.text(1.1*width, row, msg, va='center', fontsize='small')

    pl.yticks(range(len(groups)), ['%s %d' % (g_name, len(g)) for (g_name, g) in reversed(groups)], fontsize='large')
    pl.axis([-.05*width, width*1.05, -.5, len(groups)-.5])
    pl.xlabel(r'Absolute $z$-score of $p_\theta$ values', fontsize='large')

    pl.subplots_adjust(right=.5)

    def rank_by_distance_bhatt(self, qkeys, ikeys, rkeys, dists):
        """
        ::

            Reduce timbre-channel distances to ranks list by ground-truth key indices
            Bhattacharyya distance on timbre-channel probabilities and Kullback distances
        """
        # timbre-channel search using pre-computed distances
        ranks_list = []
        t_keys, t_lens = self.get_adb_lists(0) 
        rdists=pylab.ones(len(t_keys))*float('inf')
        qk = self._get_probs_tc(qkeys)
        for i in range(len(ikeys[0])): # number of include keys
            ikey=[]
            dk = pylab.zeros(self.timbre_channels)
            for t_chan in range(self.timbre_channels): # timbre channels
                ikey.append(ikeys[t_chan][i])
                try: 
                    # find dist of key i for query
                    i_idx = rkeys[t_chan].index( ikey[t_chan] ) # dataset include-key match
                    # the reduced distance function in include_keys order
                    # distance is Bhattacharyya distance on probs and dists
                    dk[t_chan] = dists[t_chan][i_idx]
                except:
                    print "Key not found in result list: ", ikey, "for query:", qkeys[t_chan]
                    raise error.BregmanError()
            rk = self._get_probs_tc(ikey)
            a_idx = t_keys.index( ikey[0] ) # audiodb include-key index
            rdists[a_idx] = distance.bhatt(pylab.sqrt(pylab.absolute(dk)), pylab.sqrt(pylab.absolute(qk*rk)))
        #search for the index of the relevant keys
        rdists = pylab.absolute(rdists)
        sort_idx = pylab.argsort(rdists)   # Sort fields into database order
        for r in self.ground_truth: # relevant keys
            ranks_list.append(pylab.where(sort_idx==r)[0][0]) # Rank of the relevant key
        return ranks_list, rdists

  def update_design(self):
    ax = self.ax
    ax.cla()
    ax2 = self.ax2
    ax2.cla()

    wp = self.wp
    ws = self.ws
    gpass = self.gpass
    gstop = self.gstop
    b, a = ss.iirdesign(wp, ws, gpass, gstop, ftype=self.ftype, output='ba')
    self.a = a
    self.b = b
    #b = [1,2]; a = [1,2]
    #Print this on command line so we can use it in our programs
    print 'b = ', pylab.array_repr(b)
    print 'a = ', pylab.array_repr(a)

    my_w = pylab.logspace(pylab.log10(.1*self.ws[0]), 0.0, num=512)
    #import pdb;pdb.set_trace()
    w, h = freqz(b, a, worN=my_w*pylab.pi)
    gp = 10**(-gpass/20.)#Go from db to regular
    gs = 10**(-gstop/20.)
    self.design_line, = ax.plot([.1*self.ws[0], self.ws[0], wp[0], wp[1], ws[1], 1.0], [gs, gs, gp, gp, gs, gs], 'ko:', lw=2, picker=5)
    ax.semilogx(w/pylab.pi, pylab.absolute(h),lw=2)
    ax.text(.5,1.0, '{:d}/{:d}'.format(len(b), len(a)))
    pylab.setp(ax, 'xlim', [.1*self.ws[0], 1.2], 'ylim', [-.1, max(1.1,1.1*pylab.absolute(h).max())], 'xticklabels', [])

    ax2.semilogx(w/pylab.pi, pylab.unwrap(pylab.angle(h)),lw=2)
    pylab.setp(ax2, 'xlim', [.1*self.ws[0], 1.2])
    ax2.set_xlabel('Normalized frequency')

    pylab.draw()

def fwhm(x, y):
	hm = pl.amax(y/2.0);
	y_diff = pl.absolute(y-hm);
	y_diff_sorted = pl.sort(y_diff);
	i1 = pl.where(y_diff==y_diff_sorted[0]);
	i2 = pl.where(y_diff==y_diff_sorted[1]);
	fwhm = pl.absolute(x[i1]-x[i2]);
	return hm, fwhm

    def fresnelSingleTransformVW(self,d) :
        # compute new window
        x2 = self.nx*pl.absolute(d)*self.wl/(self.endx-self.startx)
        y2 = self.ny*pl.absolute(d)*self.wl/(self.endy-self.starty)

        # create new intensity object
        i2 = Intensity2D(self.nx,-x2/2,x2/2,
                         self.ny,-y2/2,y2/2,
                         self.wl)

        # compute intensity
        u1p   = self.i*pl.exp(-1j*pl.pi/(d*self.wl)*(self.xgrid**2+self.ygrid**2))
        ftu1p = pl.fftshift(pl.fft2(pl.fftshift(u1p)))
        i2.i  = ftu1p*1j/(d*i2.wl)*pl.exp(-1j*pl.pi/(d*i2.wl)*(i2.xgrid**2+i2.ygrid**2))
        return i2

def fwhm_2gauss(x, y, dx=0.001):
	'''
	Finds the FWHM for the profile y(x), with accuracy dx=0.001
	Uses a 2-Gauss 1D fit.
	'''
	popt, pcov = curve_fit(gauss2, x, y);
	xx = pl.arange(pl.amin(x), pl.amax(x)+dx, dx);
	ym = gauss2(xx, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])
	hm = pl.amax(ym/2.0);
	y_diff = pl.absolute(ym-hm);
	y_diff_sorted = pl.sort(y_diff);
	i1 = pl.where(y_diff==y_diff_sorted[0]);
	i2 = pl.where(y_diff==y_diff_sorted[1]);
	fwhm = pl.absolute(xx[i1]-xx[i2]);
	return hm, fwhm, xx, ym

def calcAUC(data, y0, lag, mgr, asym, time):
    """
    Calculate the area under the curve of the logistic function
    using its integrated formula
    [ A( [A-y0] log[ exp( [4m(l-t)/A]+2 )+1 ]) / 4m ] + At
    """

    # First check that max growth rate is not zero
    # If so, calculate using the data instead of the equation
    if mgr == 0:
        auc = calcAUCData(data, time)
    else:
        timeS = time[0]
        timeE = time[-1]
        t1 = asym - y0
        #try:
        t2_s = py.log(py.exp((4 * mgr * (lag - timeS) / asym) + 2) + 1)
        t2_e = py.log(py.exp((4 * mgr * (lag - timeE) / asym) + 2) + 1)
        #except RuntimeWarning as rw:
            # Exponent is too large, setting to 10^3
        #    newexp = 1000
        #    t2_s = py.log(newexp + 1)
        #    t2_e = py.log(newexp + 1)
        t3 = 4 * mgr
        t4_s = asym * timeS
        t4_e = asym * timeE

        start = (asym * (t1 * t2_s) / t3) + t4_s
        end = (asym * (t1 * t2_e) / t3) + t4_e
        auc = end - start

    if py.absolute(auc) == float('Inf'):
        x = py.diff(time)
        auc = py.sum(x * data[1:])
    return auc

def merge_data_csvs(id):

    df = pandas.DataFrame()

    dir = dismod3.settings.JOB_WORKING_DIR % id
    #print dir
    for f in sorted(glob.glob('%s/posterior/data-*.csv'%dir)):
        #print 'merging %s' % f
        df2 = pandas.read_csv(f, index_col=None)
        df2.index = df2['index']
        df = df.drop(set(df.index)&set(df2.index)).append(df2)

    df['residual'] = df['value'] - df['mu_pred']
    df['scaled_residual'] = df['residual'] / pl.sqrt(df['value'] * (1 - df['value']) / df['effective_sample_size'])
    #df['scaled_residual'] = df['residual'] * pl.sqrt(df['effective_sample_size'])  # including 
    df['abs_scaled_residual'] = pl.absolute(df['scaled_residual'])

    d = .005 # TODO: save delta in these files, use negative binomial to calc logp
    df['logp'] = [mc.negative_binomial_like(x*n, (p+1e-3)*n, d*(p+1e-3)*n) for x,p,n in zip(df['value'], df['mu_pred'], df['effective_sample_size'])]
    df['logp'][df['data_type'] == 'rr'] = df['scaled_residual'][df['data_type'] == 'rr']

    df = df.sort('logp')

    #print df.filter('data_type area age_start age_end year_start sex effective_sample_size value residual logp'.split())[:25]
    return df

    def rank_by_distance_avg(self, qkeys, ikeys, rkeys, dists):
        """
        ::

            Reduce timbre-channel distances to ranks list by ground-truth key indices
            Kullback distances
        """
        # timbre-channel search using pre-computed distances
        ranks_list = []
        t_keys, t_lens = self.get_adb_lists(0) 
        rdists=pylab.ones(len(t_keys))*float('inf')
        for t_chan in range(self.timbre_channels): # timbre channels
            t_keys, t_lens = self.get_adb_lists(t_chan) 
            for i, ikey in enumerate(ikeys[t_chan]): # include keys, results
                try: 
                    # find dist of key i for query
                    i_idx = rkeys[t_chan].index( ikey ) # lower_bounded include-key index
                    a_idx = t_keys.index( ikey ) # audiodb include-key index
                    # the reduced distance function in include_keys order
                    # distance is the sum for now
                    if t_chan:
                        rdists[a_idx] += dists[t_chan][i_idx]
                    else:
                        rdists[a_idx] = dists[t_chan][i_idx]
                except:
                    print "Key not found in result list: ", ikey, "for query:", qkeys[t_chan]
                    raise error.BregmanError()
        #search for the index of the relevant keys
        rdists = pylab.absolute(rdists)
        sort_idx = pylab.argsort(rdists)   # Sort fields into database order
        for r in self.ground_truth: # relevant keys
            ranks_list.append(pylab.where(sort_idx==r)[0][0]) # Rank of the relevant key
        return ranks_list, rdists

 def unimodal_rate(f=rate, age_indices=age_indices, tau=1.e5):
     df = pl.diff(f[age_indices])
     sign_changes = pl.find((df[:-1] > NEARLY_ZERO) & (df[1:] < -NEARLY_ZERO))
     sign = pl.ones(len(age_indices)-2)
     if len(sign_changes) > 0:
         change_age = sign_changes[len(sign_changes)/2]
         sign[change_age:] = -1.
     return -tau*pl.dot(pl.absolute(df[:-1]), (sign * df[:-1] < 0))

def _calculate_spectra_sussix(sx, sy, Q_x, Q_y, Q_s, n_lines):

    n_turns, n_files = sx.shape

    # Allocate memory for output.        
    oxx, axx = plt.zeros((n_lines, n_files)), plt.zeros((n_lines, n_files))
    oyy, ayy = plt.zeros((n_lines, n_files)), plt.zeros((n_lines, n_files))

    # Initialise Sussix object.
    SX = PySussix.Sussix()
    
    x, xp, y, yp = sx.real, sx.imag, sy.real, sy.imag
    for file_i in xrange(n_files):
        SX.sussix_inp(nt1=1, nt2=n_turns, idam=2, ir=0, tunex=Q_x[file_i] % 1, tuney=Q_y[file_i] % 1)
        SX.sussix(x[:,file_i], xp[:,file_i], y[:,file_i], yp[:,file_i], sx[:,file_i], sx[:,file_i])

        # Amplitude normalisation
        SX.ax /= plt.amax(SX.ax)
        SX.ay /= plt.amax(SX.ay)

        # Tunes
        SX.ox = plt.absolute(SX.ox)
        SX.oy = plt.absolute(SX.oy)
        if file_i==0:
            tunexsx = SX.ox[plt.argmax(SX.ax)]
            tuneysx = SX.oy[plt.argmax(SX.ay)]
            print "\n*** Tunes from Sussix"
            print "    tunex", tunexsx, ", tuney", tuneysx, "\n"

        # Tune normalisation
        SX.ox = (SX.ox - (Q_x[file_i] % 1)) / Q_s[file_i]
        SX.oy = (SX.oy - (Q_y[file_i] % 1)) / Q_s[file_i]
    
        # Sort
        CX = plt.rec.fromarrays([SX.ox, SX.ax], names='ox, ax')
        CX.sort(order='ax')
        CY = plt.rec.fromarrays([SX.oy, SX.ay], names='oy, ay')
        CY.sort(order='ay')
        ox, ax, oy, ay = CX.ox, CX.ax, CY.oy, CY.ay
        oxx[:,file_i], axx[:,file_i], oyy[:,file_i], ayy[:,file_i] = ox, ax, oy, ay

    spectra = {}
    spectra['horizontal'] = (oxx, axx)
    spectra['vertical']   = (oyy, ayy)
        
    return spectra

def xyamb(xytab,qu,xyout=''):

    mytb=taskinit.tbtool()

    if not isinstance(qu,tuple):
        raise Exception,'qu must be a tuple: (Q,U)'

    if xyout=='':
        xyout=xytab
    if xyout!=xytab:
        os.system('cp -r '+xytab+' '+xyout)

    QUexp=complex(qu[0],qu[1])
    print 'Expected QU = ',qu   # , '  (',pl.angle(QUexp)*180/pi,')'

    mytb.open(xyout,nomodify=False)

    QU=mytb.getkeyword('QU')['QU']
    P=pl.sqrt(QU[0,:]**2+QU[1,:]**2)

    nspw=P.shape[0]
    for ispw in range(nspw):
        st=mytb.query('SPECTRAL_WINDOW_ID=='+str(ispw))
        if (st.nrows()>0):
            q=QU[0,ispw]
            u=QU[1,ispw]
            qufound=complex(q,u)
            c=st.getcol('CPARAM')
            fl=st.getcol('FLAG')
            xyph0=pl.angle(pl.mean(c[0,:,:][pl.logical_not(fl[0,:,:])]),True)
            print 'Spw = '+str(ispw)+': Found QU = '+str(QU[:,ispw])  # +'   ('+str(pl.angle(qufound)*180/pi)+')'
            #if ( (abs(q)>0.0 and abs(qu[0])>0.0 and (q/qu[0])<0.0) or
            #     (abs(u)>0.0 and abs(qu[1])>0.0 and (u/qu[1])<0.0) ):
            if ( pl.absolute(pl.angle(qufound/QUexp)*180/pi)>90.0 ):
                c[0,:,:]*=-1.0
                xyph1=pl.angle(pl.mean(c[0,:,:][pl.logical_not(fl[0,:,:])]),True)
                st.putcol('CPARAM',c)
                QU[:,ispw]*=-1
                print '   ...CONVERTING X-Y phase from '+str(xyph0)+' to '+str(xyph1)+' deg'
            else:
                print '      ...KEEPING X-Y phase '+str(xyph0)+' deg'
            st.close()
    QUr={}
    QUr['QU']=QU
    mytb.putkeyword('QU',QUr)
    mytb.close()
    QUm=pl.mean(QU[:,P>0],1)
    QUe=pl.std(QU[:,P>0],1)
    Pm=pl.sqrt(QUm[0]**2+QUm[1]**2)
    Xm=0.5*atan2(QUm[1],QUm[0])*180/pi

    print 'Ambiguity resolved (spw mean): Q=',QUm[0],'U=',QUm[1],'(rms=',QUe[0],QUe[1],')','P=',Pm,'X=',Xm

    stokes=[1.0,QUm[0],QUm[1],0.0]
    print 'Returning the following Stokes vector: '+str(stokes)
    
    return stokes

 def _power(self):
     if not self._have_stft:
         if not self._stft():
             return False
     fp = self._check_feature_params()
     self.POWER=(pylab.absolute(self.STFT)**2).sum(0)
     self._have_power=True
     if fp['verbosity']:
         print "Extracted POWER"
     return True

 def _power(self):
     if not self._stft():
         return False
     fp = self._check_feature_params()
     self.POWER=(P.absolute(self.STFT)**2).sum(0)
     self._have_power=True
     if self.verbosity:
         print "Extracted POWER"
     self.X=self.POWER
     return True

def plot_risetimes(a, b, **kwargs):

    # plt.ion()
    # if kwargs is not None:
    #     for key, value in kwargs.iteritems():
    #         if key == 'file_list':
    #             file_list = value
    #         if key == 'scan_line':
    #             scan_line = value
    # varray = plt.array(get_value_from_cfg(file_list, scan_line))

    n_files = a.shape[-1]
    cmap = plt.get_cmap('jet')
    c = [cmap(i) for i in plt.linspace(0, 1, n_files)]

    fig1, (ax1, ax2) = plt.subplots(1, 2, figsize=(20, 10))
    [ax.set_color_cycle(c) for ax in (ax1, ax2)]

    r = []
    for i in xrange(n_files):
        x, y = a[:,i], b[:,i]
        # xo, yo = x, y #, get_envelope(x, y)
        xo, yo = get_envelope(x, y)
        p = plt.polyfit(xo, np.log(yo), 1)

        # Right way to fit... a la Nicolas - the fit expert!
        l = ax1.plot(x, plt.log(plt.absolute(y)))
        lcolor = l[-1].get_color()
        ax1.plot(xo, plt.log(yo), color=lcolor, marker='o', mec=None)
        ax1.plot(x, p[1] + x * p[0], color=lcolor, ls='--', lw=3)

        l = ax2.plot(x, y)
        lcolor = l[-1].get_color()
        ax2.plot(xo, yo, 'o', color=lcolor)
        xi = plt.linspace(plt.amin(x), plt.amax(x))
        yi = plt.exp(p[1] + p[0] * xi)
        ax2.plot(xi, yi, color=lcolor, ls='--', lw=3)

        print p[1], p[0], 1 / p[0]
        # plt.draw()
        # ax1.cla()
        # ax2.cla()

        r.append(1/p[0])

    ax2.set_ylim(0, 1000)
    plt.figure(2)
    plt.plot(r, lw=3, c='purple')
    # plt.gca().set_ylim(0, 10000)

    # ax3 = plt.subplot(111)
    # ax3.semilogy(x, y)
    # ax3.semilogy(xo, yo)

    return r

 def _corrIntensity(self, laserNum, intensity, dist):
     #iScale = self.calib.maxIntensity - self.calib.minIntensity
     focalOffset = 256 * (1 - self.calib.caliParams[laserNum].focalDist/13100)**2
     corrIntensity = intensity + self.calib.caliParams[laserNum].focalSlope * pl.absolute(focalOffset - 256*(1-dist/65535)**2)
     
     if corrIntensity<self.calib.minIntensity: corrIntensity = self.calib.minIntensity
     if corrIntensity>self.calib.maxIntensity: corrIntensity = self.calib.maxIntensity
     
     corrIntensity = (corrIntensity-self.calib.minIntensity)/self.calib.maxIntensity
     
     return corrIntensity

    def _calculate_sussix_spectrum(self, turn, window_width):
    
        # Initialise Sussix object
        SX = PySussix.Sussix()
        SX.sussix_inp(nt1=1, nt2=window_width, idam=2, ir=0, tunex=self.q_x, tuney=self.q_y)

        tunes_x = plt.zeros(self.n_particles_to_analyse)
        tunes_y = plt.zeros(self.n_particles_to_analyse)

        n_particles_to_analyse_10th = int(self.n_particles_to_analyse/10)
        print 'Running SUSSIX analysis ...'
        for i in xrange(self.n_particles_to_analyse):
            if not i%n_particles_to_analyse_10th: print '  Particle', i
            SX.sussix(self.x[i,turn:turn+window_width], self.xp[i,turn:turn+window_width],
                      self.y[i,turn:turn+window_width], self.yp[i,turn:turn+window_width],
                      self.x[i,turn:turn+window_width], self.xp[i,turn:turn+window_width]) # this line is not used by sussix!
            tunes_x[i] = plt.absolute(SX.ox[plt.argmax(SX.ax)])
            tunes_y[i] = plt.absolute(SX.oy[plt.argmax(SX.ay)])
                
        return tunes_x, tunes_y

def _calculate_spectra_fft(sx, sy, Q_x, Q_y, Q_s, n_lines):

    n_turns, n_files = sx.shape
        
    # Allocate memory for output.
    oxx, axx = plt.zeros((n_lines, n_files)), plt.zeros((n_lines, n_files))
    oyy, ayy = plt.zeros((n_lines, n_files)), plt.zeros((n_lines, n_files))

    for file_i in xrange(n_files):
        t = plt.linspace(0, 1, n_turns)
        ax = plt.absolute(plt.fft(sx[:, file_i]))
        ay = plt.absolute(plt.fft(sy[:, file_i]))

        # Amplitude normalisation
        ax /= plt.amax(ax, axis=0)
        ay /= plt.amax(ay, axis=0)
    
        # Tunes
        if file_i==0:
            tunexfft = t[plt.argmax(ax[:n_turns/2], axis=0)]
            tuneyfft = t[plt.argmax(ay[:n_turns/2], axis=0)]
            print "\n*** Tunes from FFT"
            print "    tunex:", tunexfft, ", tuney:", tuneyfft, "\n"

        # Tune normalisation
        ox = (t - (Q_x[file_i] % 1)) / Q_s[file_i]
        oy = (t - (Q_y[file_i] % 1)) / Q_s[file_i]
    
        # Sort
        CX = plt.rec.fromarrays([ox, ax], names='ox, ax')
        CX.sort(order='ax')
        CY = plt.rec.fromarrays([oy, ay], names='oy, ay')
        CY.sort(order='ay')
        ox, ax, oy, ay = CX.ox[-n_lines:], CX.ax[-n_lines:], CY.oy[-n_lines:], CY.ay[-n_lines:]
        oxx[:,file_i], axx[:,file_i], oyy[:,file_i], ayy[:,file_i] = ox, ax, oy, ay

    spectra = {}
    spectra['horizontal'] = (oxx, axx)
    spectra['vertical']   = (oyy, ayy)
        
    return spectra