class BaseData(object):
    """Base class for handling processed data.

    In general, specialised subclasses of BaseData will be used
    to handle processed data rather than BaseData itself.

    Usage: ..........
    """
    __metaclass__ = MetaMethods

    def __init__(self):
        self.meta = PfMetaData()
        self.history = "%s > New %s" %(datetime.now(), self.__class__.__name__)
        if not hasattr(self, 'channels'):
            self.channels = ChannelList()
        
    def save(self):
        if pyfusion.orm_manager.IS_ACTIVE:
            # this may be inefficient: get it working, then get it fast
            self.channels.save()
            session = pyfusion.orm_manager.Session()
            session.add(self)
            session.commit()
            session.close() class BaseOrderedDataSet(object):
    __metaclass__ = MetaMethods

    def __init__(self, label=''):
        #SH added this to make it more compatible with DataSet
        self.meta = PfMetaData()
        self.created = datetime.now()
        self.label = label
        self.history = "%s > New %s" %(self.created, self.__class__.__name__)
        if label == '':
            label = unique_id()
        if not pyfusion.orm_manager.IS_ACTIVE:
            self.data_items = []
        
    def save(self):
        if pyfusion.orm_manager.IS_ACTIVE:
            session = pyfusion.orm_manager.Session()
            session.add(self)
            session.commit()
            session.close()

    def append(self, item):
        if pyfusion.orm_manager.IS_ACTIVE:
            self.data_items[len(self)] = OrderedDataSetItem(item, len(self))
        else:
            self.data_items.append(item)
    def __len__(self):
        #if pyfusion.orm_manager.IS_ACTIVE:
        #    return self.data_items.count()
        #else:
        return self.data_items.__len__()

    def __getitem__(self, key):
        if pyfusion.orm_manager.IS_ACTIVE:
            return self.data_items[key].item
        else:
            return self.data_items.__getitem__(key) class FloatDelta(BaseData):
#  Note: channel_1, 2 are channels, not channel names
    def __init__(self, channel_1, channel_2, delta, **kwargs):
        self.channel_1 = channel_1
        self.channel_2 = channel_2
        self.delta = delta
        super(FloatDelta, self).__init__(**kwargs) def cps(a,b):
    return fft.fft(a)*conjugate(fft.fft(b)) def peak_freq(signal,timebase,minfreq=0,maxfreq=1.e18):
    """
    TODO: old code: needs review
    this function only has a basic unittest to make sure it returns
    the correct freq in a simple case.
    """
    timebase = array(timebase)
    sig_fft = fft.fft(signal)
    sample_time = float(mean(timebase[1:]-timebase[:-1]))

    #SRH modification, frequencies seemed a little bit off because of the -1 in the denominator
    #Here we are trusting numpy....
    #fft_freqs = (1./sample_time)*arange(len(sig_fft)).astype(float)/(len(sig_fft)-1)
    fft_freqs = fft.fftfreq(len(sig_fft),d=sample_time)
    # only show up to nyquist freq
    new_len = len(sig_fft)/2
    sig_fft = sig_fft[:new_len]
    fft_freqs = fft_freqs[:new_len]
    [minfreq_elmt,maxfreq_elmt] = searchsorted(fft_freqs,[minfreq,maxfreq])
    sig_fft = sig_fft[minfreq_elmt:maxfreq_elmt]
    fft_freqs = fft_freqs[minfreq_elmt:maxfreq_elmt]
    
    peak_elmt = (argsort(abs(sig_fft)))[-1]
    return [fft_freqs[peak_elmt], peak_elmt] def remap_periodic(input_array, min_val, period = 2*pi):
    while len(input_array[input_array<min_val]) > 0:
        input_array[input_array<min_val] += period
    while len(input_array[input_array>=min_val+period]) > 0:
        input_array[input_array>=min_val+period] -= period
    return input_array def list2bin(input_list):
    # we explicitly cast to int(), as numpy's integer type clashes with sqlalchemy
    return int(sum(2**array(input_list))) def bin2list(input_value):
    output_list = []
    bin_index_str = _bin(input_value)[2:][::-1]
    for ind,i in enumerate(bin_index_str):
        if i == '1':
            output_list.append(ind)
    return output_list class PfMetaData(dict):
    pass def orm_register():
    def reg_item(orm_func):
        pyfusion.orm_manager.add_reg_func(orm_func)
        return orm_func
    return reg_item
from datetime import datetime
from pyfusion.data.base import BaseData, BaseOrderedDataSet, FloatDelta
from pyfusion.data.utils import cps, peak_freq, remap_periodic, list2bin, bin2list
from pyfusion.data.base import PfMetaData
from pyfusion.orm.utils import orm_register
    from sqlalchemy import Table, Column, Integer, ForeignKey, Float
    from sqlalchemy.orm import mapper, relation
import numpy as np
import pyfusion
import scipy
        #http://docs.scipy.org/doc/numpy-dev/reference/generated/numpy.fft.rfftfreq.html
        d = (self.timebase[1] - self.timebase[0])
        val = 1.0/(NFFT*d)
        N = NFFT//2 + 1
        frequency_base = np.round((np.arange(0, N, dtype=int)) * val,4)
        return FrequencyseriesData(frequency_base = frequency_base, timebase=timebase, 
                                   signal=signal,channels=self.channels,NFFT=NFFT, window='hamming',step=step)


@orm_register()
def orm_load_timeseries_data(man):
    man.tsd_table = Table('timeseriesdata', man.metadata,
                          Column('basedata_id', Integer, ForeignKey('basedata.basedata_id'), primary_key=True))
    #man.metadata.create_all()
    mapper(TimeseriesData, man.tsd_table, inherits=BaseData, polymorphic_identity='tsd')



class FrequencyseriesData(BaseData):
    def __init__(self, timebase = None, frequency_base = None, signal=None, channels=None, NFFT=None,window=None,step=None,**kwargs):
        self.timebase = timebase
        self.frequency_base = frequency_base
        self.channels = channels
        self.NFFT = NFFT
        self.window = window
        self.step = step
        self.signal = signal
        super(FrequencyseriesData, self).__init__(**kwargs)


@orm_register()
def orm_load_frequencyseries_data(man):
    man.tsd_table = Table('frequencyseriesdata', man.metadata,
                          Column('basedata_id', Integer, ForeignKey('basedata.basedata_id'), primary_key=True))
    #man.metadata.create_all()
    mapper(FrequencyseriesData, man.tsd_table, inherits=BaseData, polymorphic_identity='fsd')


class SVDData(BaseData):
    def __init__(self, chrono_labels, topo_channels, svd_input):
        """
        svd_input is a tuple as outputted by numpy.linalg.svd(data, 0)
        """
        self.channels = topo_channels
        self.chrono_labels = chrono_labels
        self.topos = np.transpose(svd_input[0])
        self.svs = svd_input[1]
        self.chronos = svd_input[2]
        self.E = sum(self.svs*self.svs)
        self.p = self.svs**2/self.E
        self.H = float((-1./np.log(len(self.svs)))*sum(self.p*np.log(self.p)))
        super(SVDData, self).__init__()

        
    def self_cps(self):
        try:
            return self._self_cps
        except AttributeError:
            self._self_cps = np.array([np.mean(cps(i,i)) for i in self.chronos])
            return self._self_cps


@orm_register()
def orm_load_svd_data(man):
    man.svd_table = Table('svddata', man.metadata,
                      Column('basedata_id', Integer, ForeignKey('basedata.basedata_id'), primary_key=True))
    #man.metadata.create_all()
    mapper(SVDData, man.svd_table, inherits=BaseData, polymorphic_identity='svddata')


class FlucStruc(BaseData):
    def __init__(self, svd_data, sv_list, timebase, min_dphase = -np.pi):
        # 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.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
        super(FlucStruc, self).__init__()

    def save(self):
        self.dphase.save()
        super(FlucStruc, self).save()

    def svs(self):
        return bin2list(self._binary_svs)

    def _get_dphase(self, min_dphase = -np.pi, get_fourier_value=0):
        """
        remap to [min_dphase, min_dphase+2pi]
        #SRH modified 24June2013 to return the complex value also
        """
        if get_fourier_value==0:
            phases = np.array([self._get_single_channel_phase(i) for i in range(self.signal.shape[0])])
        else:
            fourier_vals = []; phases = []
            for i in range(self.signal.shape[0]):
                phases_tmp, fourier_vals_tmp = self._get_single_channel_phase(i, get_fourier_value = 1)
                phases.append(phases_tmp)
                fourier_vals.append(fourier_vals_tmp)
            phases = np.array(phases)
            fourier_vals = np.array(fourier_vals)
        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())