be_odf.py

# -*- coding: utf-8 -*-
"""
Created on Tue Nov  3 15:24:12 2015

@author: Suhas Somnath, Stephen Jesse
"""

from __future__ import division, print_function, absolute_import, unicode_literals

from os import path, listdir, remove
import sys
import datetime
from warnings import warn
import h5py
import numpy as np
from scipy.io.matlab import loadmat  # To load parameters stored in Matlab .mat file

from .df_utils.be_utils import trimUDVS, getSpectroscopicParmLabel, parmsToDict, generatePlotGroups, \
    createSpecVals, requires_conjugate, nf32
from pyUSID.io.translator import Translator
from pyUSID.io.write_utils import INDICES_DTYPE, VALUES_DTYPE, Dimension, calc_chunks
from pyUSID.io.hdf_utils import write_ind_val_dsets, write_main_dataset, write_region_references, \
    create_indexed_group, write_simple_attrs, write_book_keeping_attrs, copy_attributes,\
    write_reduced_anc_dsets
from pyUSID.io.usi_data import USIDataset
from pyUSID.processing.comp_utils import get_available_memory

if sys.version_info.major == 3:
    unicode = str


class BEodfTranslator(Translator):
    """
    Translates either the Band Excitation (BE) scan or Band Excitation 
    Polarization Switching (BEPS) data format from the old data format(s) to .h5
    """

    def __init__(self, *args, **kwargs):
        super(BEodfTranslator, self).__init__(*args, **kwargs)
        self.h5_raw = None
        self.num_rand_spectra = kwargs.pop('num_rand_spectra', 1000)
        self._cores = kwargs.pop('cores', None)
        self.FFT_BE_wave = None
        self.signal_type = None
        self.expt_type = None

    @staticmethod
    def is_valid_file(data_path):
        """
        Checks whether the provided file can be read by this translator

        Parameters
        ----------
        data_path : str
            Path to raw data file

        Returns
        -------
        obj : str
            Path to file that will be accepted by the translate() function if
            this translator is indeed capable of translating the provided file.
            Otherwise, None will be returned
        """
        if not isinstance(data_path, (str, unicode)):
            raise TypeError('data_path must be a string')

        ndf = 'newdataformat'

        data_path = path.abspath(data_path)

        if path.isfile(data_path):
            ext = data_path.split('.')[-1]
            if ext.lower() not in ['jpg', 'png', 'jpeg', 'tiff', 'mat', 'txt',
                                   'dat', 'xls', 'xlsx']:
                return None
            # we only care about the folder names at this point...
            data_path, _ = path.split(data_path)

        # Check if the data is in the new or old format:
        # Check one level up:
        _, dir_name = path.split(data_path)
        if dir_name == ndf:
            # Though this translator could also read the files but the NDF Translator is more robust...
            return None
        # Check one level down:
        if ndf in listdir(data_path):
            # Though this translator could also read the files but the NDF Translator is more robust...
            return None

        file_path = path.join(data_path, listdir(path=data_path)[0])

        _, path_dict = BEodfTranslator._parse_file_path(file_path)

        if any([x.find('bigtime_0') > 0 and x.endswith('.dat') for x in path_dict.values()]):
            # This is a G-mode Line experiment:
            return None
        if any([x.find('bigtime_0') > 0 and x.endswith('.dat') for x in
                path_dict.values()]):
            # This is a G-mode Line experiment:
            return None

        parm_found = any([piece in path_dict.keys() for piece in
                          ['parm_txt', 'old_mat_parms']])
        real_found = any([piece in path_dict.keys() for piece in
                          ['read_real', 'write_real']])
        imag_found = any([piece in path_dict.keys() for piece in
                          ['read_imag', 'write_imag']])

        if parm_found and real_found and imag_found:
            if 'parm_txt' in path_dict.keys():
                return path_dict['parm_txt']
            else:
                return path_dict['old_mat_parms']
        else:
            return None

    def translate(self, file_path, show_plots=True, save_plots=True,
                  do_histogram=False, verbose=False):
        """
        Translates .dat data file(s) to a single .h5 file
        
        Parameters
        -------------
        file_path : String / Unicode
            Absolute file path for one of the data files. 
            It is assumed that this file is of the OLD data format.
        show_plots : (optional) Boolean
            Whether or not to show intermediate plots
        save_plots : (optional) Boolean
            Whether or not to save plots to disk
        do_histogram : (optional) Boolean
            Whether or not to construct histograms to visualize data quality. Note - this takes a fair amount of time
        verbose : (optional) Boolean
            Whether or not to print statements
            
        Returns
        ----------
        h5_path : String / Unicode
            Absolute path of the resultant .h5 file
        """
        file_path = path.abspath(file_path)
        (folder_path, basename) = path.split(file_path)
        (basename, path_dict) = self._parse_file_path(file_path)

        h5_path = path.join(folder_path, basename + '.h5')
        tot_bins_multiplier = 1
        udvs_denom = 2

        if 'parm_txt' in path_dict.keys():
            if verbose:
                print('\treading parameters from text file')
            (isBEPS, parm_dict) = parmsToDict(path_dict['parm_txt'])
        elif 'old_mat_parms' in path_dict.keys():
            if verbose:
                print('\treading parameters from old mat file')
            parm_dict = self.__get_parms_from_old_mat(path_dict['old_mat_parms'], verbose=verbose)
            if parm_dict['VS_steps_per_full_cycle'] == 0:
                isBEPS=False
            else:
                isBEPS=True
        else:
            raise IOError('No parameters file found! Cannot translate this dataset!')
        if verbose:
            print('\tisBEPS = {}'.format(isBEPS))

        ignored_plt_grps = []
        if isBEPS:
            parm_dict['data_type'] = 'BEPSData'

            field_mode = parm_dict['VS_measure_in_field_loops']
            std_expt = parm_dict['VS_mode'] != 'load user defined VS Wave from file'

            if not std_expt:
                raise ValueError('This translator does not handle user defined voltage spectroscopy')

            spec_label = getSpectroscopicParmLabel(parm_dict['VS_mode'])

            if parm_dict['VS_mode'] in ['DC modulation mode', 'current mode']:
                if field_mode == 'in and out-of-field':
                    tot_bins_multiplier = 2
                    udvs_denom = 1
                else:
                    if field_mode == 'out-of-field':
                        ignored_plt_grps = ['in-field']
                    else:
                        ignored_plt_grps = ['out-of-field']
            else:
                tot_bins_multiplier = 1
                udvs_denom = 1

        else:
            spec_label = 'None'
            parm_dict['data_type'] = 'BELineData'

        # Check file sizes:
        if verbose:
            print('\tChecking sizes of real and imaginary data files')

        if 'read_real' in path_dict.keys():
            real_size = path.getsize(path_dict['read_real'])
            imag_size = path.getsize(path_dict['read_imag'])
        else:
            real_size = path.getsize(path_dict['write_real'])
            imag_size = path.getsize(path_dict['write_imag'])

        if real_size != imag_size:
            raise ValueError("Real and imaginary file sizes DON'T match!. Ending")

        # Check here if a second channel for current is present
        # Look for the file containing the current data

        if verbose:
            print('\tLooking for secondary channels')
        file_names = listdir(folder_path)
        aux_files = []
        current_data_exists = False
        for fname in file_names:
            if 'AI2' in fname:
                if 'write' in fname:
                    current_file = path.join(folder_path, fname)
                    current_data_exists=True
                aux_files.append(path.join(folder_path, fname))

        add_pix = False
        num_rows = int(parm_dict['grid_num_rows'])
        num_cols = int(parm_dict['grid_num_cols'])
        if verbose:
            print('\tRows: {}, Cols: {}'.format(num_rows, num_cols))
        num_pix = num_rows * num_cols
        tot_bins = real_size / (num_pix * 4)
        # Check for case where only a single pixel is missing.
        if num_pix == 1:
            check_bins = real_size / (num_pix * 4)
        else:
            check_bins = real_size / ((num_pix - 1) * 4)

        if verbose:
            print('\tChecking bins: Total: {}, actual: {}'.format(tot_bins,
                                                                  check_bins))

        if tot_bins % 1 and check_bins % 1:
            raise ValueError('Aborting! Some parameter appears to have '
                             'changed in-between')
        elif not tot_bins % 1:
            # Everything's ok
            pass
        elif not check_bins % 1:
            tot_bins = check_bins
            warn('Warning:  A pixel seems to be missing from the data. '
                 'File will be padded with zeros.')
            add_pix = True

        tot_bins = int(tot_bins) * tot_bins_multiplier

        if 'parm_mat' in path_dict.keys():
            if verbose:
                print('\treading BE arrays from parameters text file')
            bin_inds, bin_freqs, bin_FFT, ex_wfm = self.__read_parms_mat(path_dict['parm_mat'], isBEPS)
        elif 'old_mat_parms' in path_dict.keys():
            if verbose:
                print('\treading BE arrays from old mat text file')
            bin_inds, bin_freqs, bin_FFT, ex_wfm, dc_amp_vec = self.__read_old_mat_be_vecs(path_dict['old_mat_parms'])
        else:
            if verbose:
                print('\tGenerating dummy BE arrays')
            band_width = parm_dict['BE_band_width_[Hz]'] * (0.5 - parm_dict['BE_band_edge_trim'])
            st_f = parm_dict['BE_center_frequency_[Hz]'] - band_width
            en_f = parm_dict['BE_center_frequency_[Hz]'] + band_width
            bin_freqs = np.linspace(st_f, en_f, tot_bins, dtype=np.float32)

            warn('No parms .mat file found.... Filling dummy values into ancillary datasets.')
            bin_inds = np.zeros(shape=tot_bins, dtype=np.int32)
            bin_FFT = np.zeros(shape=tot_bins, dtype=np.complex64)
            ex_wfm = np.zeros(shape=100, dtype=np.float32)

        # Forcing standardized datatypes:
        bin_inds = np.int32(bin_inds)
        bin_freqs = np.float32(bin_freqs)
        bin_FFT = np.complex64(bin_FFT)
        ex_wfm = np.float32(ex_wfm)

        self.FFT_BE_wave = bin_FFT

        if isBEPS:
            if verbose:
                print('\tBuilding UDVS table for BEPS')
            UDVS_labs, UDVS_units, UDVS_mat = self.__build_udvs_table(parm_dict, verbose=verbose)

            if verbose:
                print('\tTrimming UDVS table to remove unused plot group columns')
            UDVS_mat, UDVS_labs, UDVS_units = trimUDVS(UDVS_mat, UDVS_labs, UDVS_units, ignored_plt_grps)

            old_spec_inds = np.zeros(shape=(2, tot_bins), dtype=INDICES_DTYPE)

            # Will assume that all excitation waveforms have same num of bins
            num_actual_udvs_steps = UDVS_mat.shape[0] / udvs_denom
            bins_per_step = tot_bins / num_actual_udvs_steps
            if verbose:
                print('\t# UDVS steps: {}, # bins/step: {}'
                      ''.format(num_actual_udvs_steps, bins_per_step))

            if bins_per_step % 1:
                print('UDVS mat shape: {}, total bins: {}, bins per step: {}'.format(UDVS_mat.shape, tot_bins,
                                                                                     bins_per_step))
                raise ValueError('Non integer number of bins per step!')

            bins_per_step = int(bins_per_step)
            num_actual_udvs_steps = int(num_actual_udvs_steps)

            stind = 0
            for step_index in range(UDVS_mat.shape[0]):
                if UDVS_mat[step_index, 2] < 1E-3:  # invalid AC amplitude
                    continue
                # Bin step
                old_spec_inds[0, stind:stind + bins_per_step] = np.arange(bins_per_step, dtype=INDICES_DTYPE)
                # UDVS step
                old_spec_inds[1, stind:stind + bins_per_step] = step_index * np.ones(bins_per_step, dtype=INDICES_DTYPE)
                stind += bins_per_step
            del stind, step_index

        else:  # BE Line
            if verbose:
                print('\tPreparing supporting variables since BE-Line')
            self.signal_type = 1
            self.expt_type = 1  # Stephen has not used this index for some reason
            num_actual_udvs_steps = 1
            bins_per_step = tot_bins
            UDVS_labs = ['step_num', 'dc_offset', 'ac_amp', 'wave_type', 'wave_mod', 'be-line']
            UDVS_units = ['', 'V', 'A', '', '', '']
            UDVS_mat = np.array([1, 0, parm_dict['BE_amplitude_[V]'], 1, 1, 1],
                                dtype=np.float32).reshape(1, len(UDVS_labs))

            old_spec_inds = np.vstack((np.arange(tot_bins, dtype=INDICES_DTYPE),
                                       np.zeros(tot_bins, dtype=INDICES_DTYPE)))

        # Some very basic information that can help the processing / analysis crew
        parm_dict['num_bins'] = tot_bins
        parm_dict['num_pix'] = num_pix
        parm_dict['num_udvs_steps'] = num_actual_udvs_steps
        parm_dict['num_steps'] = num_actual_udvs_steps

        if verbose:
            print('\tPreparing UDVS slices for region references')
        udvs_slices = dict()
        for col_ind, col_name in enumerate(UDVS_labs):
            udvs_slices[col_name] = (slice(None), slice(col_ind, col_ind + 1))

        # Need to add the Bin Waveform type - infer from UDVS        
        exec_bin_vec = self.signal_type * np.ones(len(bin_inds), dtype=np.int32)

        if self.expt_type == 2:
            if verbose:
                print('\tExperiment type = 2. Doubling BE vectors')
            exec_bin_vec = np.hstack((exec_bin_vec, -1 * exec_bin_vec))
            bin_inds = np.hstack((bin_inds, bin_inds))
            bin_freqs = np.hstack((bin_freqs, bin_freqs))
            # This is wrong but I don't know what else to do
            bin_FFT = np.hstack((bin_FFT, bin_FFT))

        # Create Spectroscopic Values and Spectroscopic Values Labels datasets
        # This is an old and legacy way of doing things. Ideally, all we would need ot do is just get the unit values
        if verbose:
            print('\tCalculating spectroscopic values')
        spec_vals, spec_inds, spec_vals_labs, spec_vals_units, spec_vals_labs_names = createSpecVals(UDVS_mat,
                                                                                                     old_spec_inds,
                                                                                                     bin_freqs,
                                                                                                     exec_bin_vec,
                                                                                                     parm_dict,
                                                                                                     UDVS_labs,
                                                                                                     UDVS_units)

        if verbose:
            print('\t\tspec_vals_labs_names: {}'.format(spec_vals_labs_names))

        # Not sure what is happening here but this should work.
        spec_dim_dict = dict()
        for entry in spec_vals_labs_names:
            spec_dim_dict[entry[0] + '_parameters'] = entry[1]

        spec_vals_slices = dict()

        for row_ind, row_name in enumerate(spec_vals_labs):
            spec_vals_slices[row_name] = (slice(row_ind, row_ind + 1), slice(None))

        if path.exists(h5_path):
            if verbose:
                print('\tRemoving existing / old translated file: ' + h5_path)
            remove(h5_path)

        # First create the file
        h5_f = h5py.File(h5_path, mode='w')

        # Then write root level attributes
        global_parms = dict()
        global_parms['grid_size_x'] = parm_dict['grid_num_cols']
        global_parms['grid_size_y'] = parm_dict['grid_num_rows']
        try:
            global_parms['experiment_date'] = parm_dict['File_date_and_time']
        except KeyError:
            global_parms['experiment_date'] = '1:1:1'

        # assuming that the experiment was completed:
        global_parms['current_position_x'] = parm_dict['grid_num_cols'] - 1
        global_parms['current_position_y'] = parm_dict['grid_num_rows'] - 1
        global_parms['data_type'] = parm_dict['data_type']
        global_parms['translator'] = 'ODF'
        if verbose:
            print('\tWriting attributes to HDF5 file root')
        write_simple_attrs(h5_f, global_parms)
        write_book_keeping_attrs(h5_f)

        # Then create the measurement group
        h5_meas_group = create_indexed_group(h5_f, 'Measurement')

        # Write attributes at the measurement group level
        keys = list(parm_dict.keys())
        keys.sort()
        """
        for key in keys:
            print('{} : {}'.format(key, main_dsets[0].parent.parent.attrs[key]))
        """
        nest_parm_dict = dict()
        for key in ['FORC', 'VS', 'grid', 'BE', 'IO', 'File', 'Misc']:
            nest_parm_dict[key] = dict()
        for key in keys:
            parts = key.split('_')
            parent = 'Misc'
            rem_key = key
            if len(parts) > 1:
                if parts[0] in nest_parm_dict.keys():
                    parent = parts[0]
                    rem_key = '_'.join(parts[1:])
            nest_parm_dict[parent].update(
                {rem_key: parm_dict[key]})

        if verbose:
            print('\twriting attributes to Measurement group')
            keys = list(parm_dict.keys())
            keys.sort()
            for key in keys:
                print('\t\t{} : {}'.format(key, parm_dict[key]))
        write_simple_attrs(h5_meas_group, parm_dict)

        # Create the Channel group
        h5_chan_grp = create_indexed_group(h5_meas_group, 'Channel')

        # Write channel group attributes
        write_simple_attrs(h5_chan_grp, {'Channel_Input': 'IO_Analog_Input_1',
                                         'channel_type': 'BE'})

        # Now the datasets!
        if verbose:
            print('\tCreating ancillary datasets')
        h5_chan_grp.create_dataset('Excitation_Waveform', data=ex_wfm)

        h5_udvs = h5_chan_grp.create_dataset('UDVS', data=UDVS_mat)
        write_region_references(h5_udvs, udvs_slices, add_labels_attr=True, verbose=verbose)
        write_simple_attrs(h5_udvs, {'units': UDVS_units}, verbose=verbose)

        h5_chan_grp.create_dataset('UDVS_Indices', data=old_spec_inds[1])

        h5_chan_grp.create_dataset('Bin_Step', data=np.arange(bins_per_step, dtype=INDICES_DTYPE),
                                   dtype=INDICES_DTYPE)

        h5_chan_grp.create_dataset('Bin_Indices', data=bin_inds, dtype=INDICES_DTYPE)
        h5_chan_grp.create_dataset('Bin_Frequencies', data=bin_freqs)
        h5_chan_grp.create_dataset('Bin_FFT', data=bin_FFT)
        h5_chan_grp.create_dataset('Bin_Wfm_Type', data=exec_bin_vec)

        if verbose:
            print('\tWriting Position datasets')

        pos_dims = [Dimension('X', 'm', np.arange(num_cols)),
                    Dimension('Y', 'm', np.arange(num_rows))]
        h5_pos_ind, h5_pos_val = write_ind_val_dsets(h5_chan_grp, pos_dims, is_spectral=False, verbose=verbose)
        if verbose:
            print('\tPosition datasets of shape: {}'.format(h5_pos_ind.shape))

        if verbose:
            print('\tWriting Spectroscopic datasets of shape: {}'.format(spec_inds.shape))
        h5_spec_inds = h5_chan_grp.create_dataset('Spectroscopic_Indices', data=spec_inds, dtype=INDICES_DTYPE)        
        h5_spec_vals = h5_chan_grp.create_dataset('Spectroscopic_Values', data=np.array(spec_vals), dtype=VALUES_DTYPE)
        for dset in [h5_spec_inds, h5_spec_vals]:
            write_region_references(dset, spec_vals_slices, add_labels_attr=True, verbose=verbose)
            write_simple_attrs(dset, {'units': spec_vals_units}, verbose=verbose)
            write_simple_attrs(dset, spec_dim_dict)

        # Noise floor should be of shape: (udvs_steps x 3 x positions)
        if verbose:
            print('\tWriting noise floor dataset')
        h5_chan_grp.create_dataset('Noise_Floor', (num_pix, num_actual_udvs_steps), dtype=nf32,
                                   chunks=(1, num_actual_udvs_steps))

        """
        New Method for chunking the Main_Data dataset.  Chunking is now done in N-by-N squares
        of UDVS steps by pixels.  N is determined dynamically based on the dimensions of the
        dataset.  Currently it is set such that individual chunks are less than 10kB in size.

        Chris Smith -- csmith55@utk.edu
        """
        BEPS_chunks = calc_chunks([num_pix, tot_bins],
                                  np.complex64(0).itemsize,
                                  unit_chunks=(1, bins_per_step))
        if verbose:
            print('\tHDF5 dataset will have chunks of size: {}'.format(BEPS_chunks))
            print('\tCreating empty main dataset of shape: ({}, {})'.format(num_pix, tot_bins))
        self.h5_raw = write_main_dataset(h5_chan_grp, (num_pix, tot_bins), 'Raw_Data', 'Piezoresponse', 'V', None, None,
                                         dtype=np.complex64, chunks=BEPS_chunks, compression='gzip',
                                         h5_pos_inds=h5_pos_ind, h5_pos_vals=h5_pos_val, h5_spec_inds=h5_spec_inds,
                                         h5_spec_vals=h5_spec_vals, verbose=verbose)

        if verbose:
            print('\tReading data from binary data files into raw HDF5')
        self._read_data(UDVS_mat, parm_dict, path_dict, real_size, isBEPS,
                        add_pix, verbose=verbose)

        if verbose:
            print('\tGenerating plot groups')
        generatePlotGroups(self.h5_raw, self.mean_resp, folder_path, basename,
                           self.max_resp, self.min_resp, max_mem_mb=self.max_ram,
                           spec_label=spec_label, show_plots=show_plots, save_plots=save_plots,
                           do_histogram=do_histogram, debug=verbose)
        if verbose:
            print('\tUpgrading to USIDataset')
        self.h5_raw = USIDataset(self.h5_raw)

        # Go ahead and read the current data in the second (current) channel
        if current_data_exists:                     #If a .dat file matches
            if verbose:
                print('\tReading data in secondary channels (current)')
            self._read_secondary_channel(h5_meas_group, aux_files)

        if verbose:
            print('\tClosing HDF5 file')
        h5_f.close()

        return h5_path

    def _read_data(self, UDVS_mat, parm_dict, path_dict, real_size, isBEPS,
                   add_pix, verbose=False):
        """
        Checks if the data is BEPS or BELine and calls the correct function to read the data from
        file

        Parameters
        ----------
        UDVS_mat : numpy.ndarray of float
            UDVS table
        parm_dict : dict
            Experimental parameters
        path_dict : dict
            Dictionary of data files to be read
        real_size : dict
            Size of each data file
        isBEPS : boolean
            Is the data BEPS
        add_pix : boolean
            Does the reader need to add extra pixels to the end of the dataset
        verbose : bool, optional. Default = False
            Whether or not to print logs

        Returns
        -------
        None
        """
        # Now read the raw data files:
        if not isBEPS:
            # Do this for all BE-Line (always small enough to read in one shot)
            if verbose:
                print('\t\tReading all raw data for BE-Line in one shot')
            self.__quick_read_data(path_dict['read_real'], path_dict['read_imag'], parm_dict['num_udvs_steps'])
        elif real_size < self.max_ram and parm_dict['VS_measure_in_field_loops'] == 'out-of-field':
            # Do this for out-of-field BEPS ONLY that is also small (256 MB)
            if verbose:
                print('\t\tReading all raw BEPS (out-of-field) data in one shot')
            self.__quick_read_data(path_dict['read_real'], path_dict['read_imag'], parm_dict['num_udvs_steps'])
        elif real_size < self.max_ram and parm_dict['VS_measure_in_field_loops'] == 'in-field':
            # Do this for in-field only
            if verbose:
                print('\t\tReading all raw BEPS (in-field only) data in one shot')
            self.__quick_read_data(path_dict['write_real'], path_dict['write_imag'], parm_dict['num_udvs_steps'])
        else:
            # Large BEPS datasets OR those with in-and-out of field
            if verbose:
                print('\t\tReading all raw data for in and out of filed OR very large file')
            self.__read_beps_data(path_dict, UDVS_mat.shape[0], parm_dict['VS_measure_in_field_loops'], add_pix)
        self.h5_raw.file.flush()

    def __read_beps_data(self, path_dict, udvs_steps, mode, add_pixel=False):
        """
        Reads the imaginary and real data files pixelwise and writes to the H5 file 
        
        Parameters 
        --------------------
        path_dict : dictionary
            Dictionary containing the absolute paths of the real and imaginary data files
        udvs_steps : unsigned int
            Number of UDVS steps
        mode : String / Unicode
            'in-field', 'out-of-field', or 'in and out-of-field'
        add_pixel : boolean. (Optional; default is False)
            If an empty pixel worth of data should be written to the end             
        
        Returns 
        -------------------- 
        None
        """

        print('---- reading pixel-by-pixel ----------')

        bytes_per_pix = self.h5_raw.shape[1] * 4
        step_size = self.h5_raw.shape[1] / udvs_steps

        if mode == 'out-of-field':
            parsers = [BEodfParser(path_dict['read_real'], path_dict['read_imag'],
                                   self.h5_raw.shape[0], bytes_per_pix)]
        elif mode == 'in-field':
            parsers = [BEodfParser(path_dict['write_real'], path_dict['write_imag'],
                                   self.h5_raw.shape[0], bytes_per_pix)]
        elif mode == 'in and out-of-field':
            # each file will only have half the udvs steps:
            if 0.5 * udvs_steps % 1:
                raise ValueError('Odd number of UDVS')

            udvs_steps = int(0.5 * udvs_steps)
            # be careful - each pair contains only half the necessary bins - so read half
            parsers = [BEodfParser(path_dict['write_real'], path_dict['write_imag'],
                                   self.h5_raw.shape[0], int(bytes_per_pix / 2)),
                       BEodfParser(path_dict['read_real'], path_dict['read_imag'],
                                   self.h5_raw.shape[0], int(bytes_per_pix / 2))]

            if step_size % 1:
                raise ValueError('strange number of bins per UDVS step. Exiting')

            step_size = int(step_size)

        rand_spectra = self.__get_random_spectra(parsers, self.h5_raw.shape[0], udvs_steps, step_size,
                                                 num_spectra=self.num_rand_spectra)
        take_conjugate = requires_conjugate(rand_spectra, cores=self._cores)

        self.mean_resp = np.zeros(shape=(self.h5_raw.shape[1]), dtype=np.complex64)
        self.max_resp = np.zeros(shape=(self.h5_raw.shape[0]), dtype=np.float32)
        self.min_resp = np.zeros(shape=(self.h5_raw.shape[0]), dtype=np.float32)

        numpix = self.h5_raw.shape[0]
        """ 
        Don't try to do the last step if a pixel is missing.   
        This will be handled after the loop. 
        """
        if add_pixel:
            numpix -= 1

        for pix_indx in range(numpix):
            if self.h5_raw.shape[0] > 5:
                if pix_indx % int(round(self.h5_raw.shape[0] / 10)) == 0:
                    print('Reading... {} complete'.format(round(100 * pix_indx / self.h5_raw.shape[0])))

            # get the raw stream from each parser
            pxl_data = list()
            for prsr in parsers:
                pxl_data.append(prsr.read_pixel())

            # interleave if both in and out of field
            # we are ignoring user defined possibilities...
            if mode == 'in and out-of-field':
                in_fld = pxl_data[0]
                out_fld = pxl_data[1]

                in_fld_2 = in_fld.reshape(udvs_steps, step_size)
                out_fld_2 = out_fld.reshape(udvs_steps, step_size)
                raw_mat = np.empty((udvs_steps * 2, step_size), dtype=out_fld.dtype)
                raw_mat[0::2, :] = in_fld_2
                raw_mat[1::2, :] = out_fld_2
                raw_vec = raw_mat.reshape(in_fld.size + out_fld.size).transpose()
            else:
                raw_vec = pxl_data[0]  # only one parser
            self.max_resp[pix_indx] = np.max(np.abs(raw_vec))
            self.min_resp[pix_indx] = np.min(np.abs(raw_vec))
            self.mean_resp = (1 / (pix_indx + 1)) * (raw_vec + pix_indx * self.mean_resp)

            if take_conjugate:
                raw_vec = np.conjugate(raw_vec)
            self.h5_raw[pix_indx, :] = np.complex64(raw_vec[:])
            self.h5_raw.file.flush()

        # Add zeros to main_data for the missing pixel. 
        if add_pixel:
            self.h5_raw[-1, :] = 0 + 0j

        print('---- Finished reading files -----')

    def __quick_read_data(self, real_path, imag_path, udvs_steps):
        """
        Returns information about the excitation BE waveform present in the .mat file

        Parameters
        -----------
        real_path : String / Unicode
            Absolute file path of the real data file
        imag_path : String / Unicode
            Absolute file path of the real data file
        udvs_steps : unsigned int
            Number of UDVS steps
        """
        print('---- reading all data at once ----------')

        parser = BEodfParser(real_path, imag_path, self.h5_raw.shape[0], self.h5_raw.shape[1] * 4)

        step_size = self.h5_raw.shape[1] / udvs_steps
        rand_spectra = self.__get_random_spectra([parser], self.h5_raw.shape[0], udvs_steps, step_size,
                                                 num_spectra=self.num_rand_spectra)
        take_conjugate = requires_conjugate(rand_spectra, cores=self._cores)
        raw_vec = parser.read_all_data()
        if take_conjugate:
            print('Taking conjugate to ensure positive Quality factors')
            raw_vec = np.conjugate(raw_vec)

        if raw_vec.shape != np.prod(self.h5_raw.shape):
            percentage_padded = 100 * (np.prod(self.h5_raw.shape) - raw_vec.shape) / np.prod(self.h5_raw.shape)
            print('Warning! Raw data length {} is not matching placeholder length {}. '
                  'Padding zeros for {}% of the data!'.format(raw_vec.shape, np.prod(self.h5_raw.shape), percentage_padded))

            padded_raw_vec = np.zeros(np.prod(self.h5_raw.shape), dtype = np.complex64)

            padded_raw_vec[:raw_vec.shape[0]] = raw_vec
            raw_mat = padded_raw_vec.reshape(self.h5_raw.shape[0], self.h5_raw.shape[1])
        else:
            raw_mat = raw_vec.reshape(self.h5_raw.shape[0], self.h5_raw.shape[1])


        # Write to the h5 dataset:
        self.mean_resp = np.mean(raw_mat, axis=0)
        self.max_resp = np.amax(np.abs(raw_mat), axis=0)
        self.min_resp = np.amin(np.abs(raw_mat), axis=0)
        self.h5_raw[:, :] = np.complex64(raw_mat)
        self.h5_raw.file.flush()

        print('---- Finished reading files -----')

    @staticmethod
    def _parse_file_path(data_filepath):
        """
        Returns the basename and a dictionary containing the absolute file paths for the
        real and imaginary data files, text and mat parameter files in a dictionary
        
        Parameters 
        --------------------
        data_filepath: String / Unicode
            Absolute path of any file in the same directory as the .dat files
        
        Returns 
        --------------------
        basename : String / Unicode
            Basename of the dataset      
        path_dict : Dictionary
            Dictionary containing absolute paths of all necessary data and parameter files
        """
        (folder_path, basename) = path.split(data_filepath)
        (super_folder, basename) = path.split(folder_path)

        if basename.endswith('_d') or basename.endswith('_c'):
            # Old old data format where the folder ended with a _d or _c to denote a completed spectroscopic run
            basename = basename[:-2]
        """
        A single pair of real and imaginary files are / were generated for:
            BE-Line and BEPS (compiled version only generated out-of-field or 'read')
        Two pairs of real and imaginary files were generated for later BEPS datasets
            These have 'read' and 'write' prefixes to denote out or in field respectively
        """
        path_dict = dict()

        for file_name in listdir(folder_path):
            abs_path = path.join(folder_path, file_name)
            if file_name.endswith('.txt') and file_name.find('parm') > 0:
                path_dict['parm_txt'] = abs_path
            elif file_name.find('.mat') > 0:
                if file_name.find('more_parms') > 0:
                    path_dict['parm_mat'] = abs_path
                elif file_name == (basename + '.mat'):
                    path_dict['old_mat_parms'] = abs_path
            elif file_name.endswith('.dat'):
                # Need to account for the second AI channel here
                file_tag = 'read'
                if file_name.find('write') > 0:
                    file_tag = 'write'
                if file_name.find('real') > 0:
                    file_tag += '_real'
                elif file_name.find('imag') > 0:
                    file_tag += '_imag'
                path_dict[file_tag] = abs_path

        return basename, path_dict

    def _read_secondary_channel(self, h5_meas_group, aux_file_path):
        """
        Reads secondary channel stored in AI .mat file
        Currently works for in-field measurements only, but should be updated to
        include both in and out of field measurements

        Parameters
        -----------
        h5_meas_group : h5 group
            Reference to the Measurement group
        aux_file_path : String / Unicode
            Absolute file path of the secondary channel file.
        """
        print('---- Reading Secondary Channel  ----------')
        if len(aux_file_path)>1:
            print('Detected multiple files, assuming in and out of field')
            aux_file_paths = aux_file_path
        else:
            aux_file_paths = list(aux_file_path)

        freq_index = self.h5_raw.spec_dim_labels.index('Frequency')
        num_pix = self.h5_raw.shape[0]
        spectral_len = 1

        for i in range(len(self.h5_raw.spec_dim_sizes)):
            if i == freq_index:
                continue
            spectral_len = spectral_len * self.h5_raw.spec_dim_sizes[i]

        #num_forc_cycles = self.h5_raw.spec_dim_sizes[self.h5_raw.spec_dim_labels.index("FORC")]
        #num_dc_steps =  self.h5_raw.spec_dim_sizes[self.h5_raw.spec_dim_labels.index("DC_Offset")]

        # create a new channel
        h5_current_channel_group = create_indexed_group(h5_meas_group, 'Channel')

        # Copy attributes from the main channel
        copy_attributes(self.h5_raw.parent, h5_current_channel_group)

        # Modify attributes that are different
        write_simple_attrs(h5_current_channel_group, {'Channel_Input': 'IO_Analog_Input_2',
                                                      'channel_type': 'Current'}, verbose=True)

        #Get the reduced dimensions
        h5_current_spec_inds, h5_current_spec_values = write_reduced_anc_dsets(h5_current_channel_group,
                                                        self.h5_raw.h5_spec_inds,
                                                        self.h5_raw.h5_spec_vals, 'Frequency', is_spec=True)


        h5_current_main = write_main_dataset(h5_current_channel_group,  # parent HDF5 group
                                             (num_pix, spectral_len),  # shape of Main dataset
                                             'Raw_Data',  # Name of main dataset
                                             'Current',  # Physical quantity contained in Main dataset
                                             'nA',  # Units for the physical quantity
                                             None,  # Position dimensions
                                             None,  # Spectroscopic dimensions
                                             h5_pos_inds=self.h5_raw.h5_pos_inds,
                                             h5_pos_vals=self.h5_raw.h5_pos_vals,
                                             h5_spec_inds=h5_current_spec_inds,
                                             h5_spec_vals=h5_current_spec_values,
                                             dtype=np.float32,  # data type / precision
                                             main_dset_attrs={'IO_rate': 4E+6, 'Amplifier_Gain': 9})

        # Now calculate the number of positions that can be stored in memory in one go.
        b_per_position = np.float32(0).itemsize * spectral_len

        max_pos_per_read = int(np.floor((get_available_memory()) / b_per_position))

        # if self._verbose:
        print('Allowed to read {} pixels per chunk'.format(max_pos_per_read))

        #Open the read and write files and write them to the hdf5 file
        for aux_file in aux_file_paths:
            if 'write' in aux_file:
                infield = True
            else:
                infield=False

            cur_file = open(aux_file, "rb")

            start_pix = 0

            while start_pix < num_pix:
                end_pix = min(num_pix, start_pix + max_pos_per_read)

                # TODO: Fix for when it won't fit in memory.

                #if max_pos_per_read * b_per_position > num_pix * b_per_position:
                cur_data = np.frombuffer(cur_file.read(), dtype='f')
                #else:
                #cur_data = np.frombuffer(cur_file.read(max_pos_per_read * b_per_position), dtype='f')

                cur_data = cur_data.reshape(end_pix - start_pix, spectral_len//2)

                # Write to h5
                if infield:
                    h5_current_main[start_pix:end_pix, ::2] = cur_data
                else:
                    h5_current_main[start_pix:end_pix, 1::2] = cur_data
                start_pix = end_pix


    @staticmethod
    def __read_old_mat_be_vecs(file_path):
        """
        Returns information about the excitation BE waveform present in the 
        more parms.mat file
        
        Parameters 
        --------------------
        filepath : String or unicode
            Absolute filepath of the .mat parameter file
        
        Returns 
        --------------------
        bin_inds : 1D numpy unsigned int array
            Indices of the excited and measured frequency bins
        bin_w : 1D numpy float array
            Excitation bin Frequencies
        bin_FFT : 1D numpy complex array
            FFT of the BE waveform for the excited bins
        BE_wave : 1D numpy float array
            Band Excitation waveform
        dc_amp_vec_full : 1D numpy float array
            spectroscopic waveform. 
            This information will be necessary for fixing the UDVS for AC modulation for example
        """
        matread = loadmat(file_path, squeeze_me=True)
        BE_wave = matread['BE_wave']
        bin_inds = matread['bin_ind'] - 1  # Python base 0
        bin_w = matread['bin_w']
        dc_amp_vec_full = matread['dc_amp_vec_full']
        FFT_full = np.fft.fftshift(np.fft.fft(BE_wave))
        bin_FFT = np.conjugate(FFT_full[bin_inds])
        return bin_inds, bin_w, bin_FFT, BE_wave, dc_amp_vec_full

    @staticmethod
    def __infer_frac_phase(slopes):
        """
        Infers the VS cycle fraction and phase when parameters were
        stored in old mat files

        Parameters
        --------------------
        slopes : list / tuple
            Array of mean slopes of each fraction of a SINGLE cycle

        Returns
        --------------------
        tuple:
            fraction : float
                Fraction of VS cycle
            phase : float
                Phase offset for VS cycle
        """
        if all([_ > 0 for _ in slopes]):
            return 0.25, 0
        elif all([_ < 0 for _ in slopes]):
            return 0.25, 0.75
        elif all([_ > 0 for _ in slopes[:2]]) and all(
                [_ < 0 for _ in slopes[2:]]):
            return 0.5, 0
        elif all([_ < 0 for _ in slopes[:2]]) and all(
                [_ > 0 for _ in slopes[2:]]):
            return 0.5, 0.5
        elif all([_ > 0 for _ in slopes[:1]]) and all(
                [_ < 0 for _ in slopes[1:]]):
            return 0.75, 0
        elif all([_ > 0 for _ in slopes[:3]]) and all(
                [_ < 0 for _ in slopes[3:]]):
            return 0.75, 0.25
        elif all([_ < 0 for _ in slopes[:1]]) and all(
                [_ > 0 for _ in slopes[1:]]):
            return 0.75, 0.5
        elif all([_ < 0 for _ in slopes[:3]]) and all(
                [_ > 0 for _ in slopes[3:]]):
            return 0.75, 0.75
        elif slopes[0] > 0 and slopes[1] < 0 and slopes[2] < 0 and slopes[
            3] > 0:
            return 1, 0
        elif slopes[0] < 0 and slopes[1] > 0 and slopes[2] > 0 and slopes[
            3] < 0:
            return 1, 0.5
        else:
            return 0, 0

    @staticmethod
    def __get_parms_from_old_mat(file_path, verbose=False):
        """
        Formats parameters found in the old parameters .mat file into a dictionary
        as though the dataset had a parms.txt describing it

        Parameters
        --------------------
        file_path : Unicode / String
            absolute filepath of the .mat file containing the parameters
        verbose : bool, optional, default = False
            Whether or not to print statemetns for debugging purposes

        Returns
        --------------------
        parm_dict : dictionary
            Parameters describing experiment
        """
        parm_dict = dict()
        matread = loadmat(file_path, squeeze_me=True)

        parm_dict['IO_rate'] = str(int(matread['AO_rate'] / 1E+6)) + ' MHz'