def get_data(self,):
        self.data = pf.getDevice('H1').acq.getdata(self.shot, self.array).reduce_time([self.start_time, self.end_time])
        self.data = self.data.subtract_mean(copy=False).normalise(method='v',separate=True,copy=False)
        self.timebase = self.data.timebase
        self.data_fft = self.data.generate_frequency_series(self.samples,self.samples/self.overlap)

        print self.other_arrays, self.other_array_labels
        if self.other_arrays == None: self.other_array_labels = []
        if self.other_arrays == None: self.other_arrays = []; 
        if self.meta_data == None : self.meta_data = []

        self.other_arrays_fft = []
        for i in self.other_arrays:
            tmp = pf.getDevice('H1').acq.getdata(self.shot, i).change_time_base(self.data.timebase)
            self.other_arrays_fft.append(tmp.generate_frequency_series(self.samples,self.samples/self.overlap))

        self.instance_array_list = []
        self.misc_data_dict = {}

        #How to deal with the static case?
        for i in self.other_array_labels: 
            if i[0]!=None:  self.misc_data_dict[i[0]] = []
            if i[1]!=None:  self.misc_data_dict[i[1]] = []

        #self.fs_values = ['p','a12','H','freq','E']
        for i in self.meta_data: self.misc_data_dict[i]=[]

def extract_polarisation_data(current_shot):
    try:
        #MDSTree=MDS.Tree('mirnov',current_shot)
        coil_1x=pf.getDevice('H1').acq.getdata(current_shot,'H1ToroidalMirnov_1x')
        coil_1y=pf.getDevice('H1').acq.getdata(current_shot,'H1ToroidalMirnov_1y')
        coil_1z=pf.getDevice('H1').acq.getdata(current_shot,'H1ToroidalMirnov_1z')
        #Need to narrow the time down.... maybe extract the above data elsewhere
        #print 'successful extraction of pyrex coil'
        return coil_1x, coil_1y, coil_1z
    except:
        print 'Error getting polarisation data'
        return 0,0,0

    def __init__(self, shot, i_diag, v_diag, dev_name="W7X", debug=debug, plot=1, verbose=0, params=None):
        self.dev = pyfusion.getDevice(dev_name)
        self.shot = shot
        self.verbose = verbose
        self.i_diag = i_diag
        self.v_diag = v_diag
        self.debug = debug
        self.plot = plot
        self.select = None
        self.t_comp = (0.1,0.2)
        self.params = params
        self.figs = []
        self.suffix = ''  # this gets put at the end of the fig name (title bar)

        self.imeasfull = self.dev.acq.getdata(shot, i_diag)
        self.vmeasfull = self.dev.acq.getdata(shot, v_diag)
        comlen = min(len(self.vmeasfull.timebase), len(self.imeasfull.timebase))
        FFT_size = nice_FFT_size(comlen-2, -1)
        # the minus 2 is a fudge to hide small inconsistencies in reduce_time
        # e.g. 20160310 9 W7X_L5_LPALLI
        self.imeasfull = self.imeasfull.reduce_time([self.imeasfull.timebase[0], self.imeasfull.timebase[FFT_size]])
        self.vmeasfull = self.vmeasfull.reduce_time([self.vmeasfull.timebase[0], self.vmeasfull.timebase[FFT_size]])

        if self.params is not None:
            self.process_swept_Langmuir(**self.params)

 def test_single_mirnov_channel_kappah_from_metadata(self):
     h1test = pyfusion.getDevice("H1")
     # shot_kh = (58073, 0.74)
     shot_kh = (58123, 0.74)
     # TODO: why doesn't this work with thick client??
     data = h1test.acq.getdata(shot_kh[0], "H1_mirnov_array_1_coil_1")
     print(shot_kh)

 def test_plot_signals(self):
     # in this position, the local data/test.cfg is used, so to run separately,
     # you need to point your PYTHON_CONFIG_FILE at that file.
     # - doesn't work as of Mar 7   - 
     dev = pyfusion.getDevice('H1')
     print('\n'.join(pyfusion.conf.utils.dump(eol='')))
     dat = dev.acq.getdata(58123,'Test_H1_multi_small')
     dat.plot_signals()

    def test_single_mirnov_channel_kappah_as_argument(self):
        d = pyfusion.getDevice("H1")
        data = d.acq.getdata(58123, "H1_mirnov_array_1_coil_1")
        self.assertTrue(isinstance(data, TimeseriesData))
        from pyfusion.data.base import PfMetaData

        self.assertTrue(isinstance(data.meta, PfMetaData))
        """

    def test_device_getdata_single_shot(self):
        dev = pyfusion.getDevice("TestDevice")

        # what we want...
        expected_data = dev.acq.getdata(12345, "test_timeseries_shot_unique")
        # what we get....
        data = dev.getdata(12345, "test_timeseries_shot_unique")
        self.assertEqual(expected_data, data)

    def test_get_data(self):
        """Check that we end up with the correct data class starting from Device"""
        from pyfusion import getDevice

        test_device = getDevice(self.listed_device)
        test_data = test_device.acquisition.getdata(self.shot_number, timeseries_test_channel_1)
        from pyfusion.data.timeseries import TimeseriesData

        self.assertTrue(isinstance(test_data, TimeseriesData))

def get_array_data(current_shot, array_name, time_window=None,new_timebase=None):
    array_cutoff_locs = [0]
    data = pf.getDevice('H1').acq.getdata(current_shot, array_name)
    if new_timebase!=None:
        print('interpolating onto a new timebase....')
        data = data.change_time_base(new_timebase)
    if time_window!=None:
        data = data.reduce_time(time_window)
    return data

def extract_data(current_shot,array):
    tries,success=(0,0)
    while tries<10 and success==0:
        try:
            data=pf.getDevice('H1').acq.getdata(current_shot,array)
            success=1
            #print 'Data extracted on Shot : %d'%(current_shot)
        except (MDS.TdiException, MDS.TreeException) as e:
            print current_shot, e
            tries=tries+1
            time.sleep(0.5)
            data=None
    return data

    def test_shot_flucstrucs(self):
        """Just check that the number  of flucstrucs is the same whether
        we  use flucstruc  directly on  the shot  data with  the segment
        kwarg or whether we explicitly call the segment method.

        """
        n_samples = 90
        dev = pyfusion.getDevice("TestDevice")

        # version with explicit calls to segment:
        explicit_data = dev.getdata(12345,"test_multichannel_timeseries_large")
        explicit_dataset = pyfusion.data.base.DataSet()        
        for seg in explicit_data.segment(n_samples):
            explicit_dataset.update(seg.flucstruc())

        # version using flucstruc segment shortcut
        shortcut_flucstrucs = dev.getdata(12345,"test_multichannel_timeseries_large").flucstruc(segment=n_samples)
        self.assertEqual(len(explicit_dataset), len(shortcut_flucstrucs))

    def test_device_getdatat_multishot(self):
        dev = pyfusion.getDevice("TestDevice")

        diag = "test_timeseries_shot_unique"

        shot_list = [100, 200, 300]

        # expected to get a dataset with the data from each shot
        expected_dataset = pyfusion.data.base.DataSet()
        for shot in shot_list:
            expected_dataset.add(dev.acq.getdata(shot, diag))
        
        dataset = dev.getdata(shot_list, diag)

        # TODO: better checking if two datasets are same.
        ## using sets checks if the object is same, not using __eq__ ?
        for i in dataset:
            is_in_other = False
            for j in expected_dataset:
                if i == j:
                    is_in_other = True
            if not is_in_other:
                assert False

def single_shot_fluc_strucs(shot=None, array=None, other_arrays=None, other_array_labels=None, start_time=0.001, end_time = 0.08, samples=1024, power_cutoff = 0.1, n_svs = 2, overlap = 4, meta_data=None):
    '''This function will extract all the important information from a
    flucstruc and put it into the form that is useful for clustering
    using hte clustering module.

    SH: 8June2013 '''
    print(os.getpid(), shot)
    time_bounds = [start_time, end_time]

    #extract data for array, naked_coil and ne_array, then reduce_time, interpolate etc...
    data = pf.getDevice('H1').acq.getdata(shot, array).reduce_time([start_time, end_time])
    data = data.subtract_mean(copy=False).normalise(method='v',separate=True,copy=False)
    data_segmented = data.segment(samples,overlap=overlap, datalist = 1)
    print(other_arrays, other_array_labels)
    if other_arrays is None: other_array_labels = []
    if other_arrays is None: other_arrays = []; 
    if meta_data is None : meta_data = []

    #Get the naked coil and interferometer array if required
    #Need to include the standard interferometer channels somehow.
    other_arrays_segmented = []
    for i in other_arrays:
        tmp = pf.getDevice('H1').acq.getdata(shot, i).change_time_base(data.timebase)
        other_arrays_segmented.append(tmp.segment(samples, overlap = overlap, datalist = 1))
    instance_array_list = []
    misc_data_dict = {'RMS':[],'time':[], 'svs':[]}

    #How to deal with the static case?
    for i in other_array_labels: 
        if i[0]!=None:  misc_data_dict[i[0]] = []
        if i[1]!=None:  misc_data_dict[i[1]] = []
    #This should probably be hard coded in... 
    fs_values = ['p','a12','H','freq','E']
    for i in meta_data: misc_data_dict[i]=[]
    for i in fs_values: misc_data_dict[i]=[]

    #Cycle through the time segments looking for flucstrucs
    for seg_loc in range(len(data_segmented)):
        data_seg = data_segmented[seg_loc]
        time_seg_average_time = np.mean([data_seg.timebase[0],data_seg.timebase[-1]])
        fs_set = data_seg.flucstruc()
        #Need to check my usage of rfft.... seems different to scipy.fftpack.rfft approach
        other_arrays_data_fft = []
        for i in other_arrays_segmented:
            other_arrays_data_fft.append(np.fft.rfft(i[seg_loc].signal)/samples)
            if not np.allclose(i[seg_loc].timebase,data_seg.timebase): 
                print("WARNING possible timebase mismatch between other array data and data!!!")
        d = (data_seg.timebase[1] - data_seg.timebase[0])
        val = 1.0/(samples*d)
        N = samples//2 + 1
        frequency_base = np.round((np.arange(0, N, dtype=int)) * val,4)
        #get the valid flucstrucs
        valid_fs = []
        for fs in fs_set:
            if (fs.p > power_cutoff) and (len(fs.svs()) >= n_svs): valid_fs.append(fs)
        #extract the useful information from the valid flucstrucs
        for fs in valid_fs:
            for i in fs_values: misc_data_dict[i].append(getattr(fs,i))
            misc_data_dict['svs'].append(len(fs.svs()))
            #for i in meta_values: misc_data_dict[i].append(eval(i))
            for i in meta_data:
                try:
                    misc_data_dict[i].append(copy.deepcopy(data.meta[i]))
                except KeyError:
                    misc_data_dict[i].append(None)
            misc_data_dict['RMS'].append((np.mean(data.scales**2))**0.5)
            misc_data_dict['time'].append(time_seg_average_time)
            #other array data
            tmp_loc = np.argmin(np.abs(misc_data_dict['freq'][-1]-frequency_base))
            for i,tmp_label in zip(other_arrays_data_fft, other_array_labels):
                if tmp_label[0]!=None: misc_data_dict[tmp_label[0]].append(np.abs(i[:,0]))
                if tmp_label[1]!=None: misc_data_dict[tmp_label[1]].append(np.abs(i[:,tmp_loc]))
            phases = np.array([tmp_phase.delta for tmp_phase in fs.dphase])
            phases[np.abs(phases)<0.001]=0
            instance_array_list.append(phases)
    #convert lists to arrays....
    for i in misc_data_dict.keys():misc_data_dict[i]=np.array(misc_data_dict[i])
    return np.array(instance_array_list), misc_data_dict

 def test_getDevice_from_pf(self):
     device = pyfusion.getDevice("H1")

exec(pf.utils.process_cmd_line_args())
if help==1: 
    print(__doc__) 
    exit()

#dev_name='LHD'
if dev_name == 'LHD': 
    if diag_name == '': diag_name= 'MP2010'
    if shot_number is None: shot_number = 27233
    #shot_range = range(90090, 90110)
elif dev_name.find('H1')>=0: 
    if diag_name == '': diag_name = "H1DTacqAxial"
    if shot_number is None: shot_number = 69270


device = pf.getDevice(dev_name)

try:
    old_shot
except:
    old_shot=0


if old_shot>0: # we can expect the variables to be still around, run with -i
    if (old_diag != diag_name) or (old_shot != shot_number): old_shot=0

if old_shot == 0: 
    d = device.acq.getdata(shot_number, diag_name) # ~ 50MB for 6ch 1MS. (27233MP)
    old_shot = shot_number
    old_diag = diag_name

""" Example 1, JSPF tutorial: simple density profile scan 

In this file, data was downsampled to save space in the download package.
# The following line was used to downsample the data:
run pyfusion/examples/save_to_local.py shot_list=range(86507,86517+1)  overwrite_local=1 dev_name='H1Local' diag_name='ElectronDensity15' downsample=100 local_dir='pyfusion/examples/JSPF_tutorial/local_data'
"""
import pyfusion as pf    # (we will assume these three import lines in all future examples)
import numpy as np       
import matplotlib.pyplot as plt
import os

plt.figure('Example 1')

dev = pf.getDevice('H1Local')  # open the device (choose the experiment - e.g H-1, LHD, Heliotron-J)
ne_profile, t_mid, shot = [ ], [ ], [ ]  # prepare empty lists for ne_profile, shot and time of measurement

# next line redirects pyfusion to find downsampled local data in ./local_data
pf.config.set('global','localdatapath','local_data') 

for shot_number in range(86507, 86517+1):  # the +1 ensures 86517 is the last shot
        d = dev.acq.getdata(shot_number, 'ElectronDensity15')	# a multichannel diagnostic
        sections = d.segment(n_samples=.001)   # break into time segments
	# work through each time segment, extracting the average density during that time
        for seg in sections:
            ne_profile.append(np.average(seg.signal,axis=1))   #  axis=1 -> average over time, not channel
            t_mid.append(np.average(seg.timebase))
            shot.append(shot_number)

# store the data in a DA (Dictionary of Arrays) object, which is like a DataFrame in R or python panda
myDA = pf.data.DA_datamining.DA(dict(shot=shot, ne_profile=ne_profile, t_mid=t_mid)) 

""" LHD version of example2.py
Takes several minutes.  Need to be on site and have access to LHD data.
Should be able to paste in the next two parts from example2.py after extracting this.
"""
import pyfusion as pf  # (we will assume these three import lines in all future examples)
import numpy as np
import matplotlib.pyplot as plt

plt.figure("Example 1 - LHD")

from pyfusion.utils import get_local_shot_numbers

shots = np.sort(get_local_shot_numbers("fircall", local_path="/data/datamining/cache/fircall/", number_posn=[-13, -8]))
dev = pf.getDevice("LHD")  # open the device (choose the experiment)
ne_profile, t_mid, shot = [], [], []  # prepare empty lists for ne_profile, shot and time of measurement
for shot_number in shots:  # the +1 ensures 86517 is the last shot
    d = dev.acq.getdata(shot_number, "LHD_n_e_array")  # a multichannel diagnostic
    sections = d.segment(n_samples=128)  # break into time segments
    # work through each time segment, extracting the average density during that time
    for seg in sections:
        ne_profile.append(np.average(seg.signal, axis=1))  # axis=1 -> avg over time, not channel
        t_mid.append(np.average(seg.timebase))
        shot.append(shot_number)

# store the data in a DA (Dictionary of Arrays) object, which is like a DataFrame in R or panda
myDA = pf.data.DA_datamining.DA(dict(shot=shot, ne_profile=ne_profile, t_mid=t_mid))
myDA.save("LHD_ne_profile.npz")
myDA.write_arff("ne_profile.arff", ["ne_profile"])

myDA.extract(locals())

def single_shot_fluc_strucs_new(shot, array, time_bounds, samples, fs_dictionary = None, serial_number=0, power_cutoff = 0.1, n_svs = 2, ignore_ne_fail = 0, ne_array = 0, naked_coil = 0, overlap = 4):
    #This is supposed to replace the old single_shot_fluc struc with the new
    #instance_array, misc_details_dict approach....
    data = pf.getDevice('H1').acq.getdata(shot, array).reduce_time(time_bounds)
    data = data.subtract_mean(copy=False).normalise(method='v',separate=True,copy=False)
    data_segmented = data.segment(samples,overlap=overlap,datalist = 1)

    #Get the naked coil and interferometer array if required
    if ne_array:
        ne_data = pf.getDevice('H1').acq.getdata(shot, "ElectronDensity").change_time_base(data.timebase)
        ne_data_segmented = ne_data.segment(samples,overlap=overlap,datalist=1)
    else: ne_data_segmented = [None for i in data_segmented]
    if naked_coil:
        naked_coil = pf.getDevice('H1').acq.getdata(shot, "H1ToroidalNakedCoil").change_time_base(data.timebase)
        naked_coil_segmented = naked_coil.segment(samples,overlap=overlap,datalist=1)
    else: naked_coil_segmented = [None for i in data_segmented]

    meta_data = ['kh','heating_freq','main_current','sec_current', 'shot']
    #kh, heating_freq, main_current, sec_current = data.meta['kh'], data.meta['heating_freq'],data.meta['main_current'], data.meta['sec_current'] 
    #data.meta['shot'] = shot
    #remove the mean, normalise and segment the data....

    #Segment the data....
    instance_array_list = []
    misc_data_dict = {'RMS':[],'time':[], 'svs':[]}
    if naked_coil: misc_data_dict['naked_coil'] = []
    if ne_array: misc_data_dict['ne_mode'] = []; misc_data_dict['ne_static'] = []

    meta_values = ['kh','main_current','heating_freq','shot','sec_current']
    fs_values = ['p','a12','H','freq','E']
    for i in meta_values: misc_data_dict[i]=[]
    for i in fs_values: misc_data_dict[i]=[]

    for data_seg, ne_seg, naked_coil_seg in zip(data_segmented, ne_data_segmented, naked_coil_segmented):
        time_seg_average_time=np.mean([data_seg.timebase[0],data_seg.timebase[-1]])
        fs_set = data_seg.flucstruc()
        if ne_seg!=None:
            ne_fft = np.fft.rfft(ne_seg.signal)/samples
            if not np.allclose(ne_seg.timebase,data_seg.timebase): print "WARNING possible timebase mismatch between ne_data and data!!!"
        if naked_coil_seg!=None:
            naked_fft = np.fft.rfft(naked_coil_seg.signal)/samples
            if not np.allclose(naked_coil_seg.timebase,data_seg.timebase): print "WARNING possible timebase mismatch between ne_data and naked coil!!!"
        d = (data_seg.timebase[1] - data_seg.timebase[0])
        val = 1.0/(samples*d)
        N = samples//2 + 1
        frequency_base = np.round((np.arange(0, N, dtype=int)) * val,4)
        #get the valid flucstrucs
        valid_fs = []
        #make a list of the valid flucstrucs
        for fs in fs_set:
            if (fs.p>power_cutoff) and (len(fs.svs())>=n_svs): valid_fs.append(fs)
        #extract the useful information from the valid flucstrucs
        for fs in valid_fs:
            for i in fs_values: misc_data_dict[i].append(getattr(fs,i))
            misc_data_dict['svs'].append(fs.svs())
            #for i in meta_values: misc_data_dict[i].append(eval(i))
            for i in meta_values: #misc_data_dict[i].append(eval(i))
                try:
                    misc_data_dict[i].append(copy.deepcopy(data.meta[i]))
                except KeyError:
                    misc_data_dict[i].append(None)
                    
            phases = np.array([tmp_phase.delta for tmp_phase in fs.dphase])
            phases[np.abs(phases)<0.001]=0
            instance_array_list.append(phases)
            misc_data_dict['RMS'].append((np.mean(data.scales**2))**0.5)
            misc_data_dict['time'].append(time_seg_average_time)

            #ne_data and naked_coil_data
            tmp_loc = np.argmin(np.abs(misc_data_dict['freq'][-1]-frequency_base))
            if ne_seg!=None:
                misc_data_dict['ne_static'].append(np.abs(ne_fft[:,0]))
                misc_data_dict['ne_mode'].append(ne_fft[:,tmp_loc])
            if naked_coil_seg!=None:
                misc_data_dict['naked_coil'].append(naked_fft[:,tmp_loc])
    #convert lists to arrays....
    for i in misc_data_dict.keys():misc_data_dict[i]=np.array(misc_data_dict[i])
    return np.array(instance_array_list), misc_data_dict

import pyfusion as pf
h1=pf.getDevice("LHD")
data=h1.acq.getdata(27233,'MP')
data.plot_signals()

import pyfusion as pf

pf.config.get('global','database')  #  'sqlite:///sqlite.txt'
pf.read_config('shaun_feb_2010.cfg')

from pyfusion.conf.utils import get_config_as_dict
get_config_as_dict('Device','H1')

get_config_as_dict('Diagnostic','H1PoloidalAll')


import pyfusion as pf
h1=pf.getDevice("H1")
data=h1.acq.getdata(70071,'H1ToroidalAxial')
data.meta.keys()
data.plot_signals()

# overlay fs on spectrum
run examples/plot_specgram.py dev_name='LHD' shot_number=27233 hold=0 time_range=[.35,.5] NFFT=256 noverlap=220
run examples/plot_text_pyfusion.py filename='PF2_120229_MP_27233_27233_1_256.dat' hold=1 min_e=0.8 freq_scale=1e3
colorbar();xlim(0.35,.5);ylim(0,150000)

# mode identification overlaid on spectrum - used in LHD report Feb 2012
run examples/plot_specgram.py dev_name='LHD' shot_number=27233 hold=0 time_range=[.35,1.5] NFFT=256 noverlap=220 
clim(-210,-40)
run examples/mode_identify_example_2012.py hold=1 fsfile='PF2_120229_MP_27233_27233_1_256.dat'
xlim(0.35,.5);ylim(0,150000)
xlabel('Time (s)'); ylabel('Frequency (kHz)')

# chirp following example
run examples/plot_text_pyfusion.py filename='PF2_120229_MP_27233_27233_1_256.dat' hold=1 min_e=0.8 freq_scale=1e3 plot=1 time_range=[0.35,.4]