def __init__(self, svd_data, sv_list, timebase, min_dphase = -np.pi, phase_pairs = None):
     # NOTE I'd prefer not to duplicate info here which is in svd_data - should be able to refer to that, once sqlalchemy is hooked in
     #self.topo_channels = svd_data.topo_channels
     self.channels = svd_data.channels
     #self.svs = sv_list
     if len(sv_list) == 1:
         self.a12 = 0
     else:
         # this ratio was inverted accidentally at first
         self.a12 = svd_data.svs[sv_list[1]]/svd_data.svs[sv_list[0]]
     self._binary_svs = list2bin(sv_list)
     # peak frequency for fluctuation structure
     self.timebase = timebase
     self.freq, self.freq_elmt = peak_freq(svd_data.chronos[sv_list[0]], self.timebase)
     self.phase_pairs = phase_pairs
     self.t0 = timebase[0]
     # singular value filtered signals
     self.signal = np.dot(np.transpose(svd_data.topos[sv_list,:]),
                        np.dot(np.diag(svd_data.svs.take(sv_list)), svd_data.chronos[sv_list,:]))
     # phase differences between nearest neighbour channels
     self.dphase, self.fourier_values = self._get_dphase(min_dphase=min_dphase, get_fourier_value = 1)
     self.p = np.sum(svd_data.svs.take(sv_list)**2)/svd_data.E
     self.H = svd_data.H
     self.E = svd_data.E
     debug_(pyfusion.DEBUG, 4, key='FlucStruc')
     super(FlucStruc, self).__init__()

def regenerate_dim(x):
    """ assume x in ns since epoch from the current time """
    msg = None  # msg allows us to see which shot/diag was at fault
    diffs = np.diff(x)
    # bincount needs a positive input and needs an array with N elts where N is the largest number input
    small = (diffs > 0) & (diffs < 1000000)
    sorted_diffs = np.sort(diffs[np.where(small)[0]])
    counts = np.bincount(sorted_diffs)
    bigcounts, bigvals = myhist(diffs[np.where(~small)[0]])

    if pyfusion.VERBOSE>0:
        print('[[diff, count],....]')
        print('small:', [[argc, counts[argc]] for argc in np.argsort(counts)[::-1][0:5]])
        print('big or negative:', [[bigvals[argc], bigcounts[argc]] for argc in np.argsort(bigcounts)[::-1][0:10]])

    dtns = 1 + np.argmax(counts[1:])  # skip the first position - it is 0
    # wgt0 = np.where(sorted_diffs > 0)[0]  # we are in ns, so no worry about rounding
    histo = plt.hist if pyfusion.DBG() > 1 else np.histogram
    cnts, vals = histo(x, bins=200)[0:2]
    # ignore the two end bins - hopefully there will be very few there
    wmin = np.where(cnts[1:-1] < np.max(cnts[1:-1]))[0]
    if len(wmin)>0:
        print('**********\n*********** Gap in data > {p:.2f}%'.format(p=100*len(wmin)/float(len(cnts))))
    x01111 = np.ones(len(x))  # x01111 will be all 1s except for the first elt.
    x01111[0] = 0
    errcnt = np.sum(bigcounts) + np.sum(np.sort(counts)[::-1][1:])
    if errcnt>0 or (pyfusion.VERBOSE > 0): 
        msg = str('** repaired length of {l:,}, dtns={dtns:,}, {e} erroneous utcs'
              .format(l=len(x01111), dtns=dtns, e=errcnt))

    fixedx = np.cumsum(x01111)*dtns
    wbad = np.where((x - fixedx)>1e8)[0]
    fixedx[wbad] = np.nan
    debug_(pyfusion.DEBUG, 3, key="repair", msg="repair of W7-X scrambled Langmuir timebase") 
    return(fixedx, msg)

def process_chirps(dd, shots=None, Ns=None, Ms=None, maxd = 4, dfdtunit=1000, minlen=2,
                   plot=0, hold=0, clear_dfdt=False, debug=0, verbose=0):
    """ M and N default to all values present
    >>> 
    """
    # add islead, lead, dfdt entries
    if plot>0 and hold==0: pl.clf()  # plot(hold=0) runs slower here
    num = len(dd['shot'])
    if not dd.has_key('indx'): raise LookupError('dd needs and indx!')
    indx = dd['indx']
    if clear_dfdt or not dd.has_key('dfdt'):
        print('******* initialising dfdt data ************')
        dd.update({'islead': np.zeros(num, dtype=int)})
        dd.update({'lead': np.zeros(num, dtype=int)})
        dd['lead'][:]= -1
        dd.update({'dfdt': np.zeros(num, dtype=np.float32)})
        dd['dfdt'][:]=np.nan
    # extract the shot, M and N

    warn("need to choose a 'nan' equivalent in process_chirps")
    for k in ['shot', 'M', 'N']: 

        w = np.where(dd[k]!= -4)[0]  # right now it is -4 - but effect on colorscale?
        exec("all_{0}=np.array(dd['{0}'])[w]".format(k))

    if Ms == None: 
        Ms = np.unique(all_M) 
        warn('defaulting Ms to {0}'.format(Ms))
    if Ns == None: 
        Ns = np.unique(all_N) 
        warn('defaulting Ns to {0}'.format(Ns))
    debug_(debug,2,key='start_chirps')
    for shot_ in shots:    
        if verbose>0: print("\n{0}:".format(shot_)),
        for N_ in Ns:
            if verbose>0: print(" {0}:".format(N_)),
            for M_ in Ms:
                inds = np.where((shot_ == all_shot) & 
                             (M_ == all_M) & (N_ == all_N))[0]
    # extract, but order in time
                tord=np.argsort(dd['t_mid'][inds])
                if verbose > 4:
                    print("shot {0}, len tord={1}".format(shot_,len(tord)))
                for k in dd.keys(): 
                    if not hasattr(dd[k],'keys'):
                        exec("{0}=np.array(dd['{0}'])[inds[tord]]".format(k))
                (islead, lead, dfdt) = \
                    process_chirps_m_shot(t_mid, freq, indx, maxd=maxd, 
                                          dfdtunit=dfdtunit, minlen=minlen,
                                          plot=plot, debug=debug)
    # write them back into the dict dd
                if len(islead) != 0:            
                    if verbose > 2: print("indx={0}".format(indx))
                    dd['islead'][indx] = islead
                    dd['lead'][indx] = lead
                    dd['dfdt'][indx] = dfdt
                debug_(debug, 3, key='write_back')
    print("minlen = {m}, maxd = {d}, dfdtunit={df}".
          format(m=minlen, d=maxd, df=dfdtunit))
    if plot>0: pl.show()

def plot_shot(DA, sh, ax=None, diags = None, debug=0, fontsize=None):
    """ more flexible - no need to check for errors
    """
    if fontsize != None:     
        pl.rcParams['legend.fontsize']=fontsize

    if diags == None:
        diags = 'i_p,w_p,flat_level,NBI,ech,p_frac'

    inds=np.where(sh==DA.da['shot'])[0]
    pl.rcParams['legend.fontsize']='small'
    if ax == None: ax = pl.gca()
    #(t_mid,w_p,dw_pdt,dw_pdt2,b_0,ech,NBI,p_frac,flat_level)=9*[None]
    t = DA.da['t_mid'][inds]
    b_0 = DA.da['b_0'][inds]
    if (len(np.shape(diags)) == 0): diags = diags.split(',')
    for diag in diags:
        if DA.da.has_key(diag):
            dat = DA.da[diag][inds] ; lab = diag; linestyle='-'
            if diag in 'p_frac': 
                dat=dat*100
                lab+="*100"
            elif diag in 'ech': 
                dat=dat*10
                lab+="*10"
            elif 'flat_level' in diag: 
                dat = 30+200*dat
                linestyle='--'
                
            if diag == 'p_frac': linestyle = ':'    
            ax.plot(t,dat, linestyle=linestyle, label=lab)
    pl.legend()
    debug_(debug,1,key='plot_shot')

    pl.title("{s} {b}T".format(s=sh,b=b_0[0]))

    def do_fetch(self):
        channel_length = int(self.length)
        outdata=np.zeros(1024*2*256+1)
        ##  !! really should put a wrapper around gethjdata to do common stuff          
        #  outfile is only needed if the direct passing of binary won't work
        #  with tempfile.NamedTemporaryFile(prefix="pyfusion_") as outfile:
        ierror, getrets = gethjdata.gethjdata(self.shot,channel_length,self.path,
                                              verbose=VERBOSE, opt=1, ierror=2,
                                              outdata=outdata, outname='')
        if ierror != 0:
            raise LookupError('hj Okada style data not found for {s}:{c}'.format(s=self.shot, c=self.path))

        ch = Channel(self.path,
                     Coords('dummy', (0,0,0)))

        # the intent statement causes the out var to be returned in the result lsit
        # looks like the time,data is interleaved in a 1x256 array
        # it is fed in as real*64, but returns as real*32! (as per fortran decl)
        debug_(pyfusion.DEBUG, 4, key='Heliotron_fetch', msg='after call to getdata')
        # timebase in secs (ms in raw data) - could add a preferred unit?
        # this is partly allowed for in savez_compressed, newload, and
        # for plotting, in the config file.
        # important that the 1e-3 be inside the Timebase()
        output_data = TimeseriesData(timebase=Timebase(1e-3 * getrets[1::2]),
                                 signal=Signal(getrets[2::2]), channels=ch)
        output_data.meta.update({'shot':self.shot})
        if pyfusion.VERBOSE>0: print('HJ config name',self.config_name)
        output_data.config_name = self.config_name         
        stprms = get_static_params(shot=self.shot,signal=self.path)
        if len(list(stprms)) == 0:  # maybe this should be ignored - how do we use it?
            raise LookupError(' failure to get params for {shot}:{path}'
                              .format(shot=self.shot, path=self.path))
        output_data.params = stprms
        return output_data

    def do_fetch(self):
        # Allow for movement of Mirnov signals from A14 to PXI crate
        if pyfusion.VERBOSE > 1:
            print("LHDfetch - timeseries")
        chnl = int(self.channel_number)
        dggn = self.diag_name
        # the clever "-" thing should only be used in members of multi signal diagnostics.
        # so I moved it to base.py
        # dggn = (self.diag_name.split('-'))[-1]  # remove -
        debug_(pyfusion.DEBUG, level=4, key="LHD fetch debug")

        if dggn == "FMD":
            if self.shot < 72380:
                dggn = "SX8O"
                if chnl != 0:
                    chnl = chnl + 33
                    if self.shot < 26208:
                        chnl = chnl + 1

        filename_dict = {"diag_name": dggn, "channel_number": chnl, "shot": self.shot}
        self.basename = path.join(self.filepath, data_fileformat % filename_dict)

        files_exist = path.exists(self.basename + ".dat") and path.exists(self.basename + ".prm")
        if not files_exist:
            if pyfusion.VERBOSE > 3:
                print("RETR: retrieving %d chn %d to %s" % (self.shot, int(chnl), self.filepath))
            tmp = retrieve_to_file(diagg_name=dggn, shot=self.shot, subshot=1, channel=int(chnl), outdir=self.filepath)
            if not path.exists(self.basename + ".dat") and path.exists(self.basename + ".prm"):
                raise Exception, "something is buggered."

        return fetch_data_from_file(self)

def find_clipped(sigs, fact):
    """ look for digitizer or amplifier saturation in all raw 
    signals, and return offending indices.  Digitizer saturation
    is easy, amplifier saturation is softer - need to be careful 

    sigs can be a signal, array or a list of those
    Note that restore_sin() can be used on the sweep, rather than omitting
    the points.
    """
    if isinstance(sigs,(list, tuple)):
        goodbits = np.ones(len(sigs[0]),dtype=np.bool)
        for sig in sigs:
            goodbits = goodbits & find_clipped(sig, fact)
            debug_(debug, 5, key='find_clipped')

        return(goodbits)
    
    cnts, vals = np.histogram(sigs, 50)   # about 100ms for 250k points
    
    if (np.sum(cnts[0]) > fact * np.sum(cnts[1])):
        wunder = (sigs < vals[1])
    else:
        wunder = np.zeros(len(sigs),dtype=np.bool)

    if (np.sum(cnts[-1]) > fact * np.sum(cnts[-2])):
        wover =  (sigs > vals[-2])
    else:
        wover = np.zeros(len(sigs),dtype=np.bool)

    return(~wunder & ~wover)

def get_coords_for_channel(channel_name=None, **kwargs):
    config_dict = kwargs.copy()
    if channel_name:
        config_dict.update(get_config_as_dict("Diagnostic", channel_name))
    coord_name = "dummy"
    coord_values = (0.0, 0.0, 0.0)
    transforms = []
    for k in config_dict.keys():
        if k.startswith("coords_"):
            coord_name = k[7:]
            coord_values = tuple(map(float, config_dict[k].split(",")))
    coords_instance = Coords(coord_name, coord_values)
    if "coord_transform" in config_dict:
        transform_list = pyfusion.config.pf_options("CoordTransform", config_dict["coord_transform"])
        for transform_name in transform_list:
            # this seems to return all the globals too
            transform_class_str = pyfusion.config.pf_get(
                "CoordTransform", config_dict["coord_transform"], transform_name
            )
            # so tyr to exclude the globals
            if pyfusion.config.has_option("global", transform_name):
                continue
            transform_class = import_from_str(transform_class_str)
            # if not hasattr(transform_class, 'output_coord'):
            #    raise Exception('??')
            coords_instance.load_transform(transform_class)
    debug_(pyfusion.DEBUG, 1, key=["coord", "device_name"])
    return coords_instance

    def _get_dphase(self, min_dphase = -np.pi):
        """
        remap to [min_dphase, min_dphase+2pi]
        """
        phases = np.array([self._get_single_channel_phase(i) for i in range(self.signal.shape[0])])
        debug_(pyfusion.DEBUG, 2, key='_get_dphase')
        if self.phase_pairs != None:
            tmp = []
            for a,b in self.phase_pairs:
                ind_a = self.channels.get_channel_index(a)
                ind_b = self.channels.get_channel_index(b)
                tmp.append(phases[ind_b]-phases[ind_a])
            d_phases=remap_periodic(np.array(tmp), min_val=min_dphase)
                
        else:
            d_phases = remap_periodic(phases[1:]-phases[:-1], min_val = min_dphase)
        #d_phase_dataset = OrderedDataSet(ordered_by="channel_1.name")
        d_phase_dataset = BaseOrderedDataSet('d_phase_%s' %datetime.now())
        ## append then sort should be faster than ordereddataset.add() [ fewer sorts()]
        if self.phase_pairs != None:
            for i,phase_pair in enumerate(self.phase_pairs):
                ind_a = self.channels.get_channel_index(phase_pair[0])
                ind_b = self.channels.get_channel_index(phase_pair[1])
                d_phase_dataset.append(FloatDelta(self.channels[ind_a],self.channels[ind_b],d_phases[i]))
        else:
            for i, d_ph in enumerate(d_phases):
                d_phase_dataset.append(FloatDelta(self.channels[i], self.channels[i+1], d_ph))

        #d_phase_dataset.sort()
        return d_phase_dataset

def svd(input_data):
    from timeseries import SVDData
    svddata = SVDData(input_data.timebase, input_data.channels, linalg.svd(input_data.signal, 0))
    svddata.history = input_data.history
    svddata.scales = input_data.scales # need to pass it on to caller
    if pyfusion.DEBUG>4: print("input_data.scales",input_data.scales)
    debug_(pyfusion.DEBUG, key='svd',msg='about to return from svd')
    return svddata

 def reg_item(plot_method):
     for cl_name in class_names:
         if cl_name not in plot_reg:
             plot_reg[cl_name] = [plot_method]
         else:
             plot_reg[cl_name].append(plot_method)
     debug_(pyfusion.DEBUG, 2, key='register plot')
     return plot_method

    def to_sqlalchemy(self,db = 'sqlite:///:memory:',n_recs=1000, chunk=1000):
        """ Write to an sqlachemy database 
            chunk = 500: 2000 34 element recs/sec to (big) sqllite file, 
                    1600/sec to mysql
            'mysql://[email protected]/junk'  (mysql need snans cvtd to nul
        """        
        import sqlalchemy as SA
        def cvt(val):
            if np.isnan(val): return(None)
            elif np.issubdtype(val, np.float):
                return(float(val))
            elif np.issubdtype(val, np.int):
                return(int(val))
            else:
                if self.debug>0: print('unrecognised type {t} in cvt'
                                       .format(t=type(val)))
                return(None)

        # define the table
        self.engine = SA.create_engine(db, echo=self.debug>2)
        self.metadata = SA.MetaData()
        self.fs_table = SA.Table('fs_table', self.metadata)
        (dbkeys,dbtypes)=([],[])
        for k in self.da.keys():
            arr = self.da[k]
            typ = None
            print(k)
            if hasattr(arr,'dtype'):
                if np.issubdtype(arr.dtype, np.int): 
                    typ = SA.Integer
                elif np.issubdtype(arr.dtype, np.float): 
                    typ = SA.Float
                else:
                    print('unknown dtype {d}'.format(d=arr.dtype))
                    
                if typ != None: # if it gets here, it is recognised
                    dbkeys.append(k)
                    dbtypes.append(typ)
                    self.metadata.tables['fs_table'].append_column(SA.Column(k, typ))
                    debug_(self.debug, 1)

            if self.debug>0: print(self.metadata.tables)

        if len(dbkeys)==0: return('nothing to create')
        self.metadata.create_all(self.engine)
        conn=self.engine.connect()
        if self.len > n_recs: print('Warning - only storing n_rec = {n} records'
                                    .format(n=n_recs))
        for c in range(0,min(n_recs,len(self.da[dbkeys[0]])),chunk):
            print(c, min(c+chunk, self.len))
            lst = []
            for i in range(c,min(c+chunk, min(self.len,n_recs))):
                dct = {}
                for (k,key) in enumerate(dbkeys): 
                    dct.update({key: cvt(self.da[key][i])})
                lst.append(dct)    
            if self.debug>0: print(lst)
            conn.execute(self.fs_table.insert(),lst)

def plot_shot(da, sh=None, ax=None, diags = None, marker=None, extra_diags=None, debug=0, quiet=True, fontsize=None, hold=1, **kwargs):
    """ more flexible - no need to check for errors
    Also initial version tailored to LHD, but a fudge used to detect HJ
    extra_diags are simply added to the default list.  This way a default list
    can be used simultaneously with a one or more extra diags.
    """

    if fontsize is not None:     
        pl.rcParams['legend.fontsize']=fontsize

    if diags is None:
        if 'MICRO01' in da.keys():
            diags = 'MICRO01,DIA135'
        else:
            diags = 'i_p,w_p,flat_level,NBI,ech,p_frac'

    if extra_diags is not None:
        diags += "," + extra_diags

    if sh is None: sh = da['shot'][0]

    inds=np.where(sh==da['shot'])[0]
    pl.rcParams['legend.fontsize']='small'
    if ax is None: 
        if hold==0: pl.plot(hold=0)
        ax = pl.gca()
    #(t_mid,w_p,dw_pdt,dw_pdt2,b_0,ech,NBI,p_frac,flat_level)=9*[None]
    t = da['t_mid'][inds]
    b_0 = safe_get(da, 'b_0',inds)

    if (len(np.shape(diags)) == 0): diags = diags.split(',')
    for (i,diag) in enumerate(diags):
        lab = diag   # set default label, marker
        if marker is None: marker = '-' 

        (err, dat) = eval_diag(da=da, inds=inds, diag=diag, debug=debug)
        if quiet and err is not None:
            continue

        if diag in 'p_frac': 
            dat=dat*100
            lab+="*100"
        elif diag in 'ech': 
            dat=dat*10
            lab+="*10"
        elif 'flat_level' in diag: 
            dat = 30+200*dat
            marker='--'

        if diag == 'p_frac': 
            marker = ':'    

        ax.plot(t,dat, marker, label=lab, **kwargs)
        ## no hold keyword in ax.plot #, hold=(hold & (i==0)))
    pl.legend()
    debug_(debug,1,key='plot_shot')

    pl.title("{s} {b:.3f}T".format(s=sh,b=b_0[0]))

    def fetch(self):
        """Fetch each  channel and combine into  a multichannel instance
        of :py:class:`~pyfusion.data.timeseries.TimeseriesData`.

        :rtype: :py:class:`~pyfusion.data.timeseries.TimeseriesData`
        """
 
        ## initially, assume only single channel signals
        # this base debug breakpoint will apply to all flavours of acquisition
        debug_(pyfusion.DEBUG, level=2, key='base_multi_fetch')
        ordered_channel_names = self.ordered_channel_names()
        data_list = []
        channels = ChannelList()
        timebase = None
        meta_dict={}
        if hasattr(self, 't_min') and hasattr(self, 't_max'):
            t_range = [float(self.t_min), float(self.t_max)]
        else:
            t_range = []
        for chan in ordered_channel_names:
            sgn = 1
            if chan[0]=='-': sgn = -sgn
            bare_chan = (chan.split('-'))[-1]
            fetcher_class = import_setting('Diagnostic', bare_chan, 'data_fetcher')
            tmp_data = fetcher_class(self.acq, self.shot,
                                     config_name=bare_chan).fetch()

            if len(t_range) == 2:
                tmp_data = tmp_data.reduce_time(t_range)
            channels.append(tmp_data.channels)
            # two tricky things here - tmp.data.channels only gets one channel hhere
            #config_name for a channel is attached to the multi part -
            #we need to move it to the particular channel 
            # was  channels[-1].config_name = chan
            channels[-1].config_name = tmp_data.config_name
            meta_dict.update(tmp_data.meta)
            #print(tmp_data.signal[-1], sgn)
            tmp_data.signal = sgn * tmp_data.signal
            #print(tmp_data.signal[-1], sgn)
            if timebase == None:
                timebase = tmp_data.timebase
                data_list.append(tmp_data.signal)
            else:
                if hasattr(self, 'skip_timebase_check') and self.skip_timebase_check == 'true':
                    data_list.append(tmp_data.signal)
                else:
                    try:
                        assert_array_almost_equal(timebase, tmp_data.timebase)
                        data_list.append(tmp_data.signal)
                    except:
                        raise
        signal=Signal(data_list)
        output_data = TimeseriesData(signal=signal, timebase=timebase,
                                     channels=channels)
        #output_data.meta.update({'shot':self.shot})
        output_data.meta.update(meta_dict)
        return output_data

 def test_local_access(self):
     # This tests access to a small test tree included in this directory
     shot = -1
     debug_(pyfusion.DEBUG, level=3, key='local_access',msg='enter test_local_access')
     tree_path = os.path.join(TEST_DATA_PATH, 'test_tree')
     from pyfusion.acquisition.utils import get_acq_from_config
     test_acq_class = get_acq_from_config('test_local_tree')
     test_acq = test_acq_class('test_local_tree', test_tree_path=tree_path)
     #test_acq = pyfusion.getAcquisition('test_tree', 
     test_data = test_acq.getdata(shot, 'test_signal')

def process_chirps_m_shot(t_mid, freq, indx, maxd = 4, dfdtunit=1000, minlen=2, plot=0, debug=0    ):
    """
    >>> indx = np.array([1,2,3,4]); freq = np.array([5,4,2,1])
    >>> t_mid = np.array([1,2,3,10])
    >>> process_chirps_m_shot( t_mid, freq, indx, maxd=4, dfdtunit=1, debug=0)
    3 connected of 4 chirping,
    (array([3, 0, 0, 0]), array([ 1,  1,  1, -1]), array([ -1.,  -2.,  nan,  nan], dtype=float32))
    """
    def ftdist(i_from, i_to, tm, fr, dfdtunit=dfdtunit):
        return(np.sqrt((fr[i_to]-fr[i_from])**2 + 
                       ((tm[i_to]-tm[i_from])*dfdtunit)**2))

    lines = [] # list of lists of indices for points in lines
    line = [0] # start at the first point
    if debug>3: print("processed {0},".format(len(freq))),
    for i in range(1,len(freq)):
        # if the current point is still coincident with the start, skip
        # else if it is also within maxd in ft space, append it.
        # this method has the flaw that two points at the same time 
        # terminate the line - but that should not happen.
        if t_mid[i]!=t_mid[line[-1]] and ftdist(line[-1],i,t_mid,freq) < maxd:
            line.append(i)
        else:
            if len(line) > 1: lines.append(line) # save it if len > 1
            line = [i]  # start a new one

    for line in lines:
        if plot: pl.plot(t_mid[line],freq[line],'-o')

        # if there are at least 3 points, and the df dt agrees with neighbours
        #  within 1kHz/ms, set df/dt - otherwise set to 0 (this line not impl?)
        # also should record the id of the line (-1 for no line) the points 
        # belong to, a good choice is the index (rel to the original 
        # dataset) of the first point.

    num = len(indx)
    lead = indx*0 - 1
    islead = np.zeros(num, dtype=int)
    dfdt = np.zeros(num, dtype=np.float32)
    dfdt[:]=np.nan
    chirping = 0
    for line in lines:
        if len(line) >= minlen:
            islead[line[0]] = len(line)  # marks this as the lead point
            chirping += len(line)
            # marks the points as belonging together
            lead[line] = indx[line[0]]   
            # only distinct dfdts are recorded - the last 
            # point is left without a dfdt
            dfdt[line[0:-1]] = ((freq[line[1:]] - freq[line[0:-1]])/
                                (t_mid[line[1:]] - t_mid[line[0:-1]])/
                                dfdtunit).astype(np.float32)
    print("{0} connected of {1} chirping,".format(chirping, len(freq))),
    debug_(debug, 4, key='detail')
    return(islead, lead, dfdt)

def get_flat_top(shot=54196, times=None, smooth_dt = None, maxddw = None, hold=0, debug=0):
    """ debug=1 gives a plot
    """
    if times==None: times=np.linspace(0.02,8,8020) ;  
    from pyfusion.data.signal_processing import smooth

    bp=get_basic_diagnostics(shot=shot,diags=['w_p','dw_pdt','b_0'],times=times)
    # assume sign is OK - at the moment, the code to fix sign is in merge
    # but it is inactive.  Probably should be in get_basic_diag..
    # so far, it seems that w_p and i_p are corrected - not sure 
    # about other flux loops.
    w_p = bp['w_p']
    dw = bp['dw_pdt']
    w=np.where(w_p < 1e6)[0]  # I guess this is to exclude nans
    len(w)
    cent = np.sum(w_p[w]*times[w])/np.sum(w_p[w])
    icent = np.where(times > cent)[0][0]
    print("centroid = {0:.1f}".format(cent))
    if maxddw == None: maxddw = 100
    if smooth_dt==None: smooth_dt = 0.1 # smooth for 0.05 sec
    dt = (times[1]-times[0])
    ns = int(smooth_dt/dt)
    smootharr = [ns,ns,ns]
    offs = len(smootharr*ns)  # correction for smoothing offset
    dwsm = smooth(dw,n_smooth=smootharr)  # smooth dwdt
    ddw = np.diff(dwsm)/dt  #second deriv
    # work away from the centroid until 2nd deriv exceeds maxddw
    # assume 100kJ /sec is ramp, and a change of this over a second

    wb = int(0.5*offs) + np.nanargmax(dwsm)
    we = int(0.1*offs) + np.nanargmin(dwsm) # 
    wpmax = np.nanmax(w_p)
    # used to be maxddw - too restrictive now try dwsm
    wgtrev = np.where(np.abs(dwsm[icent-offs/2::-1])> maxddw*wpmax/100)[0]
    wgtfor = np.where(np.abs(dwsm[icent-offs/2:])> maxddw*wpmax/100)[0]
    if (len(wgtrev) < 10) or (len(wgtfor) < 10): 
        print('*** flat_top not found on shot {s}'.format(s=shot))
        return (0,0,(0,0,0,0,0))

    wbf = icent - wgtrev[0]
    wef = icent + wgtfor[0]
    if debug>0:
        pl.plot(w_p,label='w_p',hold=hold)
        pl.plot(ddw,label='ddw')
        pl.plot(dwsm,linewidth=3,label='sm(dw)')
        pl.plot(dw/10,label='dw_pdt/10')
        pl.scatter([wb, wbf, icent, wef, we],[0,200,300,250,275])
        pl.plot([wb,we],[0,0],label='b--e')
        pl.plot([wbf,wef],np.ones(2)*maxddw*wpmax/100,'o-',linewidth=2,label='bf-ef')
        pl.ylim(np.array([-1.1,1.1])*max(abs(dwsm)))
        pl.title(shot)
        pl.legend()
    debug_(max(pyfusion.DEBUG, debug),2, key='flat_top')
    #return(times[wb], times[we],(wb,we,wbf,wef,icent)) # used to ignore wbf,wef
    return(times[wbf], times[wef],(wb,we,wbf,wef,icent))

    def process_constant_Langmuir(self, t_range=None, t_comp=[0, 0.1], leakage=None, threshold=0.01, threshchan=12, filename=None, amu=1, plot=None, return_data=False, suffix=''):
        '''
        This function should analyse the shots with constant Voltage, usually measuring I_sat. It will try to use as much of process_swept_Langmuir as possible and the pyfusion environment.
        '''
        self.figs = []  # reset the count of figures used to stop too many plots
        self.amu = amu
        if not isinstance(self.imeasfull.channels, (list, tuple, np.ndarray)):
            self.imeasfull.channels = [self.imeasfull.channels]

        self.i_chans = [ch.config_name for ch in self.imeasfull.channels]
        for ch in self.i_chans:  # only take the current channels, not U
            if ch[-1] != 'I':
                raise ValueError("Warning - removal of V chans doesn't work!!!")
                # hopefully we can ignore _U channels eventually, but not yet
                self.i_chans.remove(ch)

        self.coords = [ch.coords.w7_x_koord for ch in self.imeasfull.channels]
        self.get_iprobe(leakage=leakage, t_comp=t_comp)
        tb = self.iprobefull.timebase
        # the 3000 below tries to avoid glitches from Hilbert at both ends
        # only look at electron current - bigger (shot 0309.52 LP53 has a positive spike at 2s)
        w_plasma = np.where((-self.iprobefull.signal[threshchan] > threshold) & (tb > tb[3000]) &(tb < tb[-3000]))[0]
        
        if t_range is None:
            t_range = [tb[w_plasma[0]], tb[w_plasma[-1]]]

        self.t_range = t_range
        if t_range is not None:
            self.imeas = self.imeasfull.reduce_time(t_range)
            self.iprobe = self.iprobefull.reduce_time(t_range)
        else:
            self.imeas = self.imeasfull
            self.iprobe = self.iprobefull
        if self.debug>0: print(' {n} segments'.format(n=len(self.segs)))
        debug_(self.debug, 3, key='process_loop')
        if filename is not None:            
            if  '*' in filename:
                fmt = 'LP{s0}_{s1}_'
                if 'L5' in self.i_diag:
                    fmt += 'L5' + self.i_diag.split('L5')[1][0]
                filename = filename.replace('*',fmt+'_')
                if filename.endswith('_'):  # remove ending _
                    filename = filename[0:-1]
            if '{' in filename:
                filename = filename.format(s0=self.shot[0], s1=self.shot[1], i_diag=self.i_diag)
            print('writing {fn}'.format(fn=filename))
            self.write_DA(filename)
            
        if return_data:
            return(self.iprobe)

 def getdata(self, shots, diag):
     # check if shots is an iterable (list or set of shots)
     try:
         iterator = iter(shots)
     except TypeError:
         # not iterable -- assume single shot
         return self.acq.getdata(shots, diag)
     else:
         # iterable - create dataset and loop through shots
         dataset = pyfusion.data.base.DataSet()
         for shot in shots:
             dataset.add(self.acq.getdata(shot, diag))
         dataset.parent=self.name
         debug_(pyfusion.DEBUG,5, key='Device:getdata')
         return dataset

def normalise(input_data, method=None, separate=False):
    """ method=None -> default, method=0 -> DON'T normalise
    """
    from numpy import mean, sqrt, max, abs, var, atleast_2d
    from pyfusion.data.base import DataSet
    # this allows method='0'(or 0) to prevent normalisation for cleaner code
    # elsewhere
    if pyfusion.DEBUG>3: print('separate = %d' % (separate))
    if (method == 0) or (method == '0'): return(input_data)
    if (method == None) or (method.lower() == "none"): method='rms'
    if isinstance(input_data, DataSet):
        output_dataset = DataSet(input_data.label+"_normalise")
        for d in input_data:
            output_dataset.add(normalise(d, method=method, separate=separate))
        return output_dataset
    if method.lower() in ['rms', 'r']:
        if input_data.signal.ndim == 1:
            norm_value = sqrt(mean(input_data.signal**2))
        else:
            rms_vals = sqrt(mean(input_data.signal**2, axis=1))
            if separate == False:
                rms_vals = max(rms_vals)
            norm_value = atleast_2d(rms_vals).T            
    elif method.lower() in ['peak', 'p']:
        if input_data.signal.ndim == 1:
            norm_value = abs(input_data.signal).max(axis=0)
        else:
            max_vals = abs(input_data.signal).max(axis=1)
            if separate == False:
                max_vals = max(max_vals)
            norm_value = atleast_2d(max_vals).T
    elif method.lower() in ['var', 'variance', 'v']:
        if input_data.signal.ndim == 1:
            norm_value = var(input_data.signal)
        else:
            var_vals = var(input_data.signal, axis=1)
            if separate == False:
                var_vals = max(var_vals)
            norm_value = atleast_2d(var_vals).T            
    input_data.signal = input_data.signal / norm_value
    #print('norm_value = %s' % norm_value)
    norm_hist = ','.join(["{0:.2g}".format(v) for v in norm_value.flatten()])
    input_data.history += "\n:: norm_value =[{0}]".format(norm_hist)
    input_data.history += ", method={0}, separate={1}".format(method, separate)
    input_data.scales = norm_value

    debug_(pyfusion.DEBUG, key='normalise',msg='about to return from normalise')
    return input_data