def set_irf( self, irf=None, wraptime=None, dispersion=None ):
        """
        The detector response isn't a delta-function, meaning that what
        you measure isn't the true time-dependence of the system you are
        measuring. It's the time dependence of the system convolved with
        the response of the detector. This method sets the measured
        trace of the detector response that will be used to convolve
        the mult-exponential model before fitting (thereby taking this
        convolution into account without doing nearly-impossible
        numerical deconvolution).
        """
        if isinstance( irf, Trace ):
            self.irf = irf
        elif type(irf) == str:
            self.irf = Trace( irf )
            if wraptime is not None:
                self.irf.wrapcurves( wraptime )
            elif self.wraptime is not None:
                self.irf.wrapcurves( self.wraptime )
        
        self.irf_dispersion = dispersion
        if dispersion is not None:
            # this is meant to address chromatic dispersion within the setup
            # (e.g. optical fiber)
            # don't bother with normalization b/c it gets normalized to unit area below anyway.
            #original = self.irf.curves[0].copy()
            #self.irf.curves[0][dispersion:] += original[:-dispersion]
            #self.irf.curves[0][:-dispersion] += original[dispersion:]
            len1 = len(self.irf.curves[0])
            chain_of_three = pylab.zeros( 3*len1 ) # stack three curves end-to-end so cspline1d_eval doesn't have to extrapolate beyond data
            chain_of_three[:len1] = self.irf.curves[0][:]
            chain_of_three[len1:2*len1] = self.irf.curves[0][:]
            chain_of_three[-len1:] = self.irf.curves[0][:]
            g = cspline1d(chain_of_three)
            smoothed = pylab.zeros( len1 )
            std_dev = dispersion/1000.0
            for t0 in pylab.linspace(-2*std_dev, 2*std_dev, 50):
                weight = pylab.exp( -t0**2/2.0/std_dev**2 )
                smoothed += weight * cspline1d_eval( g, self.irf.t[0]-t0, dx=self.irf.t[0][1], x0=-self.irf.t[0][-1] )
            self.irf.curves[0] = smoothed
            
        normalized = self.irf.curves[0].astype(numpy.float)/float(sum(self.irf.curves[0])) # normalize integral to 1, just like delta-function!!!
        self.irf.curves[0] = normalized.copy()

        self.irf_generator = cspline1d(self.irf.curves[0])
        self.irf_dt = self.irf.t[0][1]-self.irf.t[0][0]
        self.irf_t0 = self.irf.t[0][0]
        
        if False:
            """not sure this matters if we do interpolation
            """
            # difference in degree of binning (e.g. 8ps vs. 4ps is bin difference of 2)
            bin_difference = pylab.np.int( self.resolution / self.irf.resolution )
            if bin_difference != 1:
                raise ValueError("Have not yet tested deconvolution with different resolution than detector trace!!!")
                d = self.irf.curves[0]
                detector_binned = pylab.zeros( len(d)/bin_difference )
                for i in range( len(detector_binned ) ):
                    detector_binned[i] = sum( d[i*bin_difference : i*bin_difference+bin_difference] )

def TransformSpectrum(args):
    flux, invvar, c0, c1, newwave = args
    smoother = 3.0

    fitc = cspline1d(flux, lamb=smoother)
    fitc_iv = cspline1d(invvar, lamb=smoother)
    newf = cspline1d_eval(fitc, newwave, dx=c1, x0=c0)
    newiv   = cspline1d_eval(fitc_iv, newwave, dx=c1, x0=c0)

    return (newf, newiv)

def splineIntLinExt(y,x,xNew,splinecoeff = 0.):
    """
    Use the scipy spline interpolation, but linearly extrapolate at the edges,
    since scipy.signal.cspline1d assumes periodic boundary conditions
    """    
    if len(x) < 4:
        return interpolateLin(y,x,xNew)

    if isinstance(xNew,float):
        wasfloat = 1
        xNew = M.array([xNew])
    else:
        wasfloat = 0

    whereSpline = N.where((xNew > x[0]) * (xNew < x[-1]))[0]
    whereLin = N.where((xNew <= x[0]) + (xNew >= x[-1]))[0]
    
    ans = xNew * 0.
    if len(whereSpline) > 0:
        if isinstance(splinecoeff,float): # not pre-calculated.
            splinecoeff = SS.cspline1d(y)
        ans[whereSpline] = SS.cspline1d_eval(splinecoeff, xNew[whereSpline], dx=x[1]-x[0], x0 = x[0])

    if len(whereLin) > 0:
        ans[whereLin] = interpolateLin(y,x,xNew[whereLin])

    if wasfloat:
        return ans[0]
    else:
        return ans

def splineIntBig0LittleLog(y,x,xNew,splinecoeff = 0.):
    """
    Use the scipy spline interpolation, but linearly extrapolate at the edges,
    since scipy.signal.cspline1d assumes periodic boundary conditions
    """
    if len(x) < 4:
        return interpolateLin(y,x,xNew)

    if isinstance(xNew,float):
        wasfloat = 1
        xNew = N.array([xNew])
    else:
        wasfloat = 0

    whereSpline = N.where((xNew >= x[0]) * (xNew <= x[-1]))[0]
    whereLittle = N.where(xNew < x[0])[0]
    whereBig = N.where(xNew >= x[-1])[0]
    
    ans = xNew * 0.
    if len(whereSpline) > 0:
        if isinstance(splinecoeff,float): # not pre-calculated.
            splinecoeff = SS.cspline1d(y)
        ans[whereSpline] = SS.cspline1d_eval(splinecoeff, xNew[whereSpline], dx=x[1]-x[0], x0 = x[0])

    if len(whereLittle) > 0:
        xw = xNew[whereLittle]
        logx,logy = N.log(x[:2]),N.log(y[:2])
        ans[whereLittle] = N.exp(logy[0] + (N.log(xw)-logx[0])/(logx[1]-logx[0])*(logy[1]-logy[0]))


    if wasfloat:
        return ans[0]
    else:
        return ans

def interp_T(T,Ts,PPs):
    xvals=zeros(1)
    xvals[0]=T
    Ps=list()
    for i in PPs:
        cj=cspline1d(i)
        Ps.append(cspline1d_eval(cj,xvals,dx=500.0,x0=Ts[0]))
    return Ps

def wiggle(x, origin=0, posFill='black', negFill=None, lineColor='black', 
        resampleRatio=10, rescale=False, ymin=0, ymax=None, ax=None):
    """Plots a "wiggle" trace
    Input:
        x: input data (1D numpy array)
        origin: (default, 0) value to fill above or below (float)
        posFill: (default, black) color to fill positive wiggles with (string 
            or None)
        negFill: (default, None) color to fill negative wiggles with (string 
            or None)
        lineColor: (default, black) color of wiggle trace (string or None)
        resampleRatio: (default, 10) factor to resample traces by before 
            plotting (1 = raw data) (float)
        rescale: (default, False) If True, rescale "x" to be between -1 and 1
        ymin: (default, 0) The minimum y to use for plotting
        ymax: (default, len(x)) The maximum y to use for plotting
        ax: (default, current axis) The matplotlib axis to plot onto
    Output:
        a matplotlib plot on the current axes
    """
    from matplotlib import pyplot as plt
    from scipy.signal import cspline1d, cspline1d_eval

    if ymax is None:
        ymax = x.size

    # Rescale so that x ranges from -1 to 1
    if rescale:
        x = x.astype(np.float)
        x -= x.min()
        x /= x.ptp()
        x *= 2
        x -= 1

    # Interpolate at resampleRatio x the previous density
    y = np.linspace(0, x.size, x.size)
    interp_y = np.linspace(0, x.size, x.size * resampleRatio)
    cj = cspline1d(x)
    interpX = cspline1d_eval(cj,interp_y) #,dx=1,x0=0
    newy = np.linspace(ymax, ymin, interp_y.size)
    if origin == None: 
        origin = interpX.mean()

    # Plot
    if ax is None:
        ax = plt.gca()
        plt.hold(True)
    if posFill is not None: 
        ax.fill_betweenx(newy, interpX, origin,
                where=interpX > origin,
                facecolor=posFill)
    if negFill is not None:
        ax.fill_betweenx(newy, interpX, origin,
                where=interpX < origin,
                facecolor=negFill)
    if lineColor is not None:
        ax.plot(interpX, newy, color=lineColor)

    def test_basic(self):
        y=array([1,2,3,4,3,2,1,2,3.0])
        x=arange(len(y))
        dx=x[1]-x[0]
        cj = signal.cspline1d(y)

        x2=arange(len(y)*10.0)/10.0
        y2=signal.cspline1d_eval(cj, x2, dx=dx,x0=x[0])

        # make sure interpolated values are on knot points
        assert_array_almost_equal(y2[::10], y, decimal=5)

    def run(self):
        """
        create an inifile from the parameters and run camb on it
        and store the results in k,pk
        """
        self.printIniFile()
        os.system(self.cambPath + "/camb " + self.iniName)
        self.k, self.pk = utils.readColumns(self.cp.output_root + "_matterpower.dat")

        self.logk, self.logpk = M.log(self.k), M.log(self.pk)
        self.pkSplineCoeff = SS.cspline1d(self.logpk)

    def runEHuShiftBAO(self, sig8, shiftfact=1.):
        #Output EHu file
        f = file('ehu.in','w')

        f.write((str(self.cp.omega_baryon + self.cp.omega_cdm))+', '+str(self.cp.omega_lambda)+', '+\
                str(self.cp.omega_neutrino)+', '+str(self.cp.omega_baryon)+'\n')
        f.write(str(self.cp.hubble/100.)+', '+str(self.cp.temp_cmb)+', '+str(self.cp.massless_neutrinos)+'\n')
        f.write(str(self.cp.transfer_redshift[0])+'\n')
        f.write(str(self.cp.transfer_kmax)+', '+str(self.cp.transfer_k_per_logint)+'\n')
        f.write('1\n')
        tilt = self.cp.scalar_spectral_index[0]
        f.write(str(tilt)+'\n')
        f.write('0\n')

        f.close()

        # run EHu code
        os.system('~/ehu/power < ehu.in')

        # read into c.k, c.pk
        eh = N.loadtxt('trans.dat')
        self.k = eh[:,0]*1.
        self.logk = M.log(self.k)

        if (getnowig):
            self.trans = eh[:,3]*1.
        else:
            self.trans = eh[:,1]*1.

        #M.loglog(self.k,self.trans)
        #M.show()
        if tilt == 1.:
            delH =  1.94e-5*(self.cp.omega_cdm + self.cp.omega_baryon)**(-0.785)
            delta = delH**2*(3000.0*self.k/(self.cp.hubble/100.))**4*self.trans**2
        else:
            delH =  1.94e-5*self.cp.omega_cdm**(-0.785 - 0.05*M.log(tilt))\
                   * M.exp(-0.95*(tilt - 1.) - 0.169*(tilt - 1)**2)
            delta = delH**2*(3000.0*self.k/(self.cp.hubble/100.))**(3 + tilt)*self.trans**2

        # Just an approximate normalization; really need sig8.
            
        self.pk = (2.*M.pi**2 * delta/self.k**3)*(self.cp.hubble/100.)**3

        if self.cp.transfer_redshift[0] > 0.:
            ps = PowerSpectrum(self.cp)
            sig8use = sig8*ps.d1(self.cp.transfer_redshift[0])/ps.d1(0.)
        else:
            sig8use = sig8
        normalizePk(self,sig8use) # sets c.logpk, too
        self.pkSplineCoeff = SS.cspline1d(self.logpk)

        return

def doresample(orig_x, orig_y, new_x, method='cubic', padlen=0, antialias=False):
    """
    Resample data from one spacing to another.  By default, does not apply any antialiasing filter.

    Parameters
    ----------
    orig_x
    orig_y
    new_x
    method
    padlen

    Returns
    -------

    """
    pad_y = tide_filt.padvec(orig_y, padlen=padlen)
    tstep = orig_x[1] - orig_x[0]
    if padlen > 0:
        pad_x = np.concatenate((np.arange(orig_x[0] - padlen * tstep, orig_x[0], tstep),
                                orig_x,
                                np.arange(orig_x[-1] + tstep, orig_x[-1] + tstep * (padlen + 1), tstep)))
    else:
        pad_x = orig_x
    if padlen > 0:
        print('padlen=', padlen)
        print('tstep=', tstep)
        print(pad_x)

    # antialias and ringstop filter
    init_freq = len(pad_x) / (pad_x[-1] - pad_x[0])
    final_freq = len(new_x) / (new_x[-1] - new_x[0])
    if antialias and (init_freq > final_freq):
        aafilterfreq = final_freq / 2.0
        aafilter = tide_filt.noncausalfilter(filtertype='arb', usebutterworth=False)
        aafilter.setarb(0.0, 0.0, 0.95 * aafilterfreq, aafilterfreq)
        pad_y = aafilter.apply(init_freq, pad_y)

    if method == 'cubic':
        cj = signal.cspline1d(pad_y)
        return tide_filt.unpadvec(
            np.float64(signal.cspline1d_eval(cj, new_x, dx=(orig_x[1] - orig_x[0]), x0=orig_x[0])), padlen=padlen)
    elif method == 'quadratic':
        qj = signal.qspline1d(pad_y)
        return tide_filt.unpadvec(
            np.float64(signal.qspline1d_eval(qj, new_x, dx=(orig_x[1] - orig_x[0]), x0=orig_x[0])), padlen=padlen)
    elif method == 'univariate':
        interpolator = sp.interpolate.UnivariateSpline(pad_x, pad_y, k=3, s=0)  # s=0 interpolates
        return tide_filt.unpadvec(np.float64(interpolator(new_x)), padlen=padlen)
    else:
        print('invalid interpolation method')
        return None

    def pkInterp(self, knew, calcsplinecoeff = False):
        """
        P(knew) interpolated in log space for arbitrary k.
        Call with calcsplinecoeff = True if the power spectrum (self.pk) might have changed recently.
        """
        if (calcsplinecoeff):
            self.pkSplineCoeff = SS.cspline1d(self.logpk)

        w0 = N.where(knew > 0.)[0]
        newy = 0.*knew
        # will handle the case k = 0
        newy[w0] = M.exp(utils.splineIntLinExt(self.logpk, self.logk, M.log(knew[w0]),splinecoeff = self.pkSplineCoeff))
        return newy

 def run(self):
     _, fileExtension = os.path.splitext(self.fileName)
     if fileExtension == '.gmv':
         print('Geomagnetic variation')
         with open(self.fileName, 'rt') as csvdata:
             date = []
             value = []
             for row in csv.reader(csvdata):
                 if ('#' in row[0]):
                     self.header.append(row)
                 else:
                     date.append(row[0])
                     value.append(row[1])
         self.notifyProgress.emit(20)
     elif fileExtension == '.ske':
         print('Kp estimation')
         with open(self.fileName, 'rt') as csvdata:
             date = []
             value = []
             for row in csv.reader(csvdata, delimiter=' '):
                 if ('#' in row[0]):
                     self.header.append(row)
                 else:
                     print(row)
                     if int(row[7]) < 2:
                         date.append(
                             dt.datetime.strptime(
                                 ''.join((row[0], row[1], row[2],
                                         row[4])),
                                 '%Y%m%d%H%M')),
                         value.append(float(row[-1])-float(row[-14]))  # 4h
                         # value.append(float(row[-1])-float(row[19]))  # 1h
         self.notifyProgress.emit(20)
     signal_src = np.array((date, value), dtype=np.dtype('a25'))
     signal = signal_src[:, np.logical_not(
         np.isnan(signal_src[1, :].astype(np.float)))]
     self.notifyProgress.emit(60)
     if self.interpolate:
         self.time = signal_src[0, :].astype(np.datetime64).astype(
             dt.datetime)
         dx = dates.date2num(self.time[1])-dates.date2num(self.time[0])
         cj = cspline1d(signal[1, :].astype(float))
         self.value = cspline1d_eval(cj, dates.date2num(self.time),
                                     dx=dx,
                                     x0=dates.date2num(self.time[0]))
     else:
         self.time = dates.signal[0, :].astype(np.datetime64).astype(
             dt.datetime)
         self.value = signal[1, :].astype(np.float)
     self.notifyProgress.emit(80)
     self.loaded.emit()

def SplineResample(y, n_x, x_range = None, smoother = 3):
    ''' Resample a signal y using spline.

    y: the signal must be sampled on uniform x
    n_x: number of points for the new signal
    x_range: tuple (start x, end x) of the signal
    smoother: spline smoothing strength
    '''

    if x_range is None: x_range = (0., 1.)
    spcoeffs = cspline1d(y, lamb = smoother)
    return cspline1d_eval(
        spcoeffs, linspace(x_range[0], x_range[1], n_x), 
        dx = (x_range[1]-x_range[0])/(len(y)-1.0), x0 = x_range[0])

def interpol(x1,y1,x_out,plot=False):
	from scipy.signal import cspline1d, cspline1d_eval
	#assumes that data points are evenly spaced

	dx=x1[1]-x1[0]
	cj=cspline1d(y1)
	y_out=cspline1d_eval(cj,x_out,dx=dx,x0=x1[0])
	if plot:
		from pylab import plot, show,legend
		plot(x_out,y_out,'ob',x1,y1,'xg',ms=2)
		legend(['interpolated','original'])
		show()

	return y_out

def interp_around(X_sc,s_fracpeak,s_before,s_after,kind='cubic'):
    n_c = X_sc.shape[1]
    n_s = s_before+s_after
    Out_sc = np.empty((n_s,n_c),dtype=np.float32)
    for i_c in xrange(n_c):
        if kind == 'cubic':
            coeffs = cspline1d(X_sc[:,i_c])
            Out_sc[:,i_c] = cspline1d_eval(coeffs,
                                       newx=np.arange(s_fracpeak - s_before,s_fracpeak+s_after,dtype=np.float32))
        elif kind == "linear":
            Out_sc[:,i_c] = interp1d(np.arange(X_sc.shape[0]),X_sc[:,i_c],
                                     bounds_error=True,kind=kind)(np.arange(s_fracpeak - s_before,s_fracpeak+s_after,dtype=np.float32))
        else: raise Exception("kind must be 'linear' or 'cubic'")
    return Out_sc

def evalBeams(arrBeams,guessShearLoad,spc,debugFlag = False):
    
    lastBeam = arrBeams[-1]
    arrStrainEnergy = np.zeros(np.size(arrBeams))
    arrSurfEnergy = np.zeros(np.size(arrBeams))
    lastIndex = np.size(arrBeams)
    maxSurfEnergy = -gamma * lastBeam.w*(lastBeam.Lt - lastBeam.L)
    arrResults = np.array([])
    
    for j,beam in enumerate(arrBeams):
        if debugFlag:
            print('Beam Length: %f of %f'%(beam.L*scale,beam.Lt*scale))
        
        dispSearch(beam=beam,initLoad = guessShearLoad,goal=spc/2/scale,tol=1e-7,right=0,left=0)
        guessShearLoad = beam.shearLoad
        
        if debugFlag:
            print('Solved Beam -- TipDisp: %s Goal: %s Force: %s' % (beam.yTipDisplacement()*scale,spc/2,beam.shearLoad))
        
        arrStrainEnergy[j] = beam.calculateStrainEnergy()
        arrSurfEnergy[j] = -gamma * beam.w *(beam.Lt - beam.L)
        if arrStrainEnergy[j] >= np.abs(maxSurfEnergy): 
            # since there is more bending energy than surface energy stop computing 
            print('Super stiff beam')
            lastIndex = j
            break
    
    if lastIndex > 0:   # This ensures that we have more than one data point before trying to interpolate
        interpLens = np.linspace(arrBeamLens[0],arrBeamLens[lastIndex-1],num=100,endpoint=True) # Generate x values for which to interpolate
        csFit = cspline1d((arrStrainEnergy[0:lastIndex]+arrSurfEnergy[0:lastIndex]))    # Generate cubic spline fit to the sub dataset
        interpTotalEnergy = cspline1d_eval(csFit,interpLens,dx=(arrBeamLens[1]-arrBeamLens[0]), x0 = arrBeamLens[0])    # Generate the interpolated values from the fit and x points
        finalLen = interpLens[interpTotalEnergy.argmin()]   # find the minimum of the energy balance and grab index to choose the appropriate length
        
        if debugFlag:
            print('beamLens shape: %s arrStrain: %s'%(arrBeamLens[0:lastIndex].shape,arrStrainEnergy[0:lastIndex].shape))
            mpl.figure()
            mpl.hold(True)
            mpl.plot(arrBeamLens[0:lastIndex]*scale,arrStrainEnergy[0:lastIndex]*scale)
            mpl.plot(arrBeamLens[0:lastIndex]*scale,arrSurfEnergy[0:lastIndex]*scale)
            mpl.plot(interpLens*scale,interpTotalEnergy*scale,arrBeamLens[0:lastIndex]*scale,(arrStrainEnergy+arrSurfEnergy)[0:lastIndex]*scale,'o')
        arrResults = np.array([arrBeamLens[0:lastIndex],arrStrainEnergy[0:lastIndex]])
    else:   # since there is only one datapoint then use that as the value
        finalLen = arrBeamLens[lastIndex]
        arrResults = np.array([arrBeamLens[lastIndex],arrStrainEnergy[lastIndex]])
    
    
    
    return (finalLen,arrResults)

def normalizePk(c, sig8):
    """
    renormalize the power spectrum to the given sigma_8

    Note: this doesn't take into account the redshift; it just blindly sets c.pk to have
    sigma_8 = sig8.
    """
    sig2now = sigma2fromPk(c, 8.0)
    # print 'sig2now=',sig2now
    c.pk *= sig8 ** 2 / sig2now
    c.logpk = M.log(c.pk)

    c.cp.scalar_amp[0] = c.cp.scalar_amp[0] * sig8 ** 2 / sig2now[0]  # inelegant tuple change

    # for scipy splines
    c.pkSplineCoeff = SS.cspline1d(c.logpk)
    return sig2now

def sr_interpol2(x,y,ytarget,doplot=0,factor=10):
    dx = x[1]-x[0]
    newx = linspace(min(x),max(x),factor*len(x))

    cj = cspline1d(y)
    newy = cspline1d_eval(cj, newx, dx=dx,x0=x[0])

    ysq = (ytarget-newy)**2
    index = where(ysq == min(ysq))

    if doplot:
        clf()
        plot(x,y,'o')
        plot(newx,newy)
        plot(newx[index],newy[index],'o')
        show()

    return newx[index[0][0]]

    def pkInterp(self, knew, calcsplinecoeff=False):
        """
        P(knew) interpolated in log space for arbitrary k.
        Call with calcsplinecoeff = True if the power spectrum (self.pk) might have changed recently.
        """
        if calcsplinecoeff:
            self.pkSplineCoeff = SS.cspline1d(self.logpk)

        # newy = M.exp(utils.splineIntLinExt(self.logpk, self.logk, M.log(knew),splinecoeff = self.pkSplineCoeff))
        w0 = N.where(knew <= 0.0)
        wn0 = N.where(knew > 0.0)
        newy = 0.0 * knew
        newy[wn0] = M.exp(
            utils.splineIntLinExt(self.logpk, self.logk, M.log(knew[wn0]), splinecoeff=self.pkSplineCoeff)
        )
        # do 0
        newy[w0] = 0.0
        return newy

    def fitSpline(self,degree=2):
        """
        **SUMMARY**

        A function to generate a spline curve fitting over the points in LineScan with
        order of precision given by the parameter degree

        **PARAMETERS**

        * *degree* - the precision of the generated spline 

        **RETURNS**

        The spline as a LineScan fitting over the initial values of LineScan

        **EXAMPLE**

        >>> import matplotlib.pyplot as plt
        >>> img = Image("lenna")
        >>> ls = img.getLineScan(pt1=(10,10)),pt2=(20,20)).normalize()
        >>> spline = ls.fitSpline()
        >>> plt.plot(ls)
        >>> plt.show()
        >>> plt.plot(spline)
        >>> plt.show()
        
        **NOTES**

        Implementation taken from http://www.scipy.org/Cookbook/Interpolation  

        """
        if degree > 4:
            degree = 4  # No significant improvement with respect to time usage
        if degree < 1:
            warnings.warn('LineScan.fitSpline - degree needs to be >= 1')
            return None
        retVal = None
        y = np.array(self)
        x = np.arange(0,len(y),1)
        dx = 1
        newx = np.arange(0,len(y)-1,pow(0.1,degree))
        cj = sps.cspline1d(y)
        retVal = sps.cspline1d_eval(cj,newx,dx=dx,x0=x[0])
        return retVal