be_utils.py

# -*- coding: utf-8 -*-
"""
Created on Tue Jan 05 07:55:56 2016

@author: Suhas Somnath, Chris Smith, Rama K. Vasudevan
"""

from __future__ import division, print_function, absolute_import, unicode_literals

from os import path
from warnings import warn

import h5py
import matplotlib.pyplot as plt
import numpy as np
import xlrd as xlreader

from sidpy.hdf.hdf_utils import write_simple_attrs, link_h5_objects_as_attrs, \
    get_attr, get_auxiliary_datasets
from sidpy.proc.comp_utils import get_available_memory, parallel_compute

from pyUSID.io.hdf_utils import find_dataset, create_indexed_group, \
    write_main_dataset, get_unit_values
from pyUSID.io.write_utils import create_spec_inds_from_vals, Dimension

from ....processing.histogram import build_histogram
from ....analysis.utils.be_sho import SHOestimateGuess
from ....viz.be_viz_utils import plot_1d_spectrum, plot_2d_spectrogram, \
    plot_histograms
from ...hdf_writer import HDFwriter
from ...virtual_data import VirtualDataset, VirtualGroup

nf32 = np.dtype({'names': ['super_band', 'inter_bin_band', 'sub_band'],
                 'formats': [np.float32, np.float32, np.float32]})


def generate_bipolar_triangular_waveform(cycle_pts, cycle_frac=1, phase=0,
                                         amplitude=1, cycles=1, offset=0):
    """
    Generates a bi-polar triangular waveform based on basic parameters

    Parameters
    ----------
    cycle_pts : int
        Number of parameters in a single cycle
    cycle_frac : float, optional. Default = 1
        Fraction of the cycle - ranges as (0, 1]
    phase : float, optional. Default = 0
        Phase offset for the waveform. Ranges as [0, 1)
    amplitude : float, optional. Default = 1
        Maximum amplitude for the waveform
    offset : float, optional. Default = 0
        DC offset for the waveform
    cycles : int, optional. Default = 1
        Number of repetitions or cycles

    Returns
    -------
    dc_vec
    """
    # Scale the number of points in the cycle by the (inverse of) the fraction
    actual_cycle_pts = int(cycle_pts // cycle_frac)

    # Build the default bi-polar triangular waveform
    dc_vec = np.hstack(
        (np.linspace(0, 1, actual_cycle_pts // 4, endpoint=False),
         np.linspace(1, 0, actual_cycle_pts // 4, endpoint=False),
         np.linspace(0, -1, actual_cycle_pts // 4, endpoint=False),
         np.linspace(-1, 0, actual_cycle_pts // 4, endpoint=False)))

    # Apply phase offset via a roll:
    dc_vec = np.roll(dc_vec, int(phase * actual_cycle_pts))

    # Next truncate by the fraction
    dc_vec = dc_vec[:cycle_pts]

    # Next, scale by the amplitude:
    dc_vec = amplitude * dc_vec

    # Next offset:
    dc_vec += offset

    # Finally, tile by the number of cycles
    dc_vec = np.tile(dc_vec, int(cycles))

    return dc_vec


def infer_bipolar_triangular_fraction_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


def flat_parm_dict_to_nested(parm_dict):
    """
    Creates a nested dictionary of parameters for Band EXcitation from the
    standard flat dictionary

    Parameters
    ----------
    parm_dict : dict
        Dictionary with parameters

    Returns
    -------
    nest_parm_dict : dict
        Nested dictionary of the same parameters
    """
    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]})
    return nest_parm_dict


def remove_non_exist_spec_dim_labs(h5_spec_inds, h5_spec_vals,
                                   h5_meas_grp, verbose=False):
    """
    Removes non-existent spectroscopic dimension name and units from
    attributes of spectroscopic datasets.

    Notes
    -----
    This was written mainly to clean up
    after the buggy Labview HDF5 acquisition software. This will be used in
    the LabviewHDF5Patcher Translator

    Parameters
    ----------
    h5_spec_inds : h5py.Dataset
        Dataset containing the spectroscopic indices
    h5_spec_vals : h5py.Dataset
        Dataset containing the spectroscopic values
    h5_meas_grp : h5py.Group
        Group containing all the parameters for the BE measurement
    verbose : bool, optional. Default = False
        Whether or not to print statements aiding in debugging

    Returns
    -------
    None
    """
    for obj, name, exp_type in zip([h5_spec_inds, h5_spec_vals, h5_meas_grp],
                                   ['h5_spec_inds', 'h5_spec_vals', 'h5_meas_grp'],
                                   [h5py.Dataset, h5py.Dataset, h5py.Group]):
        if not isinstance(obj, exp_type):
            raise TypeError('{} should be a {}. Provided object was: {}'
                            ''.format(name, exp_type, type(obj)))

    spec_dim_names = get_attr(h5_spec_inds, 'labels')
    spec_dim_units = get_attr(h5_spec_inds, 'units')
    if len(spec_dim_names) != len(spec_dim_units):
        raise ValueError('Unqeual lengths for the spec dim names and units')
    if len(spec_dim_names) <= h5_spec_inds.shape[0]:
        return

    if verbose:
        print(
            'extra dimensions in the list of names attributes: {} compared to '
            'rows in dataset: {}'
            ''.format(len(spec_dim_names), h5_spec_inds.shape[0]))

    if len(spec_dim_names) - h5_spec_inds.shape[0] > 1:
        raise NotImplementedError(
            'Cannot handle case when more than one dimensions are fake')

    # Gather basic parameters for each dimension
    # h5_meas_grp = h5_raw.parent.parent
    field_type = get_attr(h5_meas_grp, 'VS_measure_in_field_loops')
    num_freq_bins = int(get_attr(h5_meas_grp, 'num_bins') /
                        get_attr(h5_meas_grp, 'num_UDVS_steps'))
    num_fields = 1 + int(all([targ in field_type for targ in ['in', 'out']]))
    dc_off_steps = get_attr(h5_meas_grp, 'VS_steps_per_full_cycle')
    num_cycles = get_attr(h5_meas_grp, 'VS_number_of_cycles')
    try:
        num_forc_cycles = get_attr(h5_meas_grp, 'VS_num_of_FORC_cycles')
    except KeyError:
        # This is named differently in some h5 files
        num_forc_cycles = get_attr(h5_meas_grp, 'VS_num_of_Forc_cycles')

    size_dict = {'Frequency': num_freq_bins,
                 'DC_Offset': dc_off_steps,
                 'Field': num_fields,
                 'Cycle': num_cycles,
                 'FORC': num_forc_cycles,
                 }

    if verbose:
        for key, val in size_dict.items():
            print('{} : {}'.format(key, val))

    matched_dims = list()
    matched_units = list()

    row_ind = 0

    for dim_name, dim_units in zip(spec_dim_names, spec_dim_units):
        # Pass one row at a time:
        this_dict = get_unit_values(np.expand_dims(h5_spec_inds[row_ind],
                                                   axis=0),
                                    np.expand_dims(h5_spec_vals[row_ind],
                                                   axis=0),
                                    is_spec=True, all_dim_names=[dim_name])

        std_dim_name = dim_name
        if dim_name == 'FORC_Cycle':
            std_dim_name = 'FORC'

        if verbose:
            print(dim_name, len(this_dict[dim_name]), size_dict[std_dim_name])

        if len(this_dict[dim_name]) == size_dict[std_dim_name]:
            row_ind += 1
            matched_dims.append(dim_name)
            matched_units.append(dim_units)
        else:
            if verbose:
                print(
                    '\tDimension: "' + dim_name + '" did not match with what '
                                                  'was in dataset')

    if verbose:
        print('Writing new attributes to spectroscopic datasets')

    new_attrs = {'labels': matched_dims, 'units': matched_units}
    for h5_dset in [h5_spec_inds, h5_spec_vals]:
        write_simple_attrs(h5_dset, new_attrs)


def parmsToDict(filepath, parms_to_remove=[]):
    """
    Translates the parameters in the text file into a dictionary. 
    Also indentifies whether this is a BEPS or BELine dataset.
    
    Parameters
    -----------
    filepath : String / Unicode
        Absolute path of the parameters text or spreadsheet file.
    parms_to_remove : List of string (Optional)
        keys that this function should attempt to remove from the dictionary
    
    Returns
    ----------
    isBEPS : Boolean
        whether this dataset is BEPS or BE Line.
    parm_dict : Dictionary
        experimental parameters
    """
    verbose = False

    lines = list()

    if filepath.lower().endswith('.txt'):
        file_handle = open(filepath, 'r')
        raw_lines = file_handle.readlines()
        file_handle.close()
        for line in raw_lines:
            line = line.rstrip()
            if len(line) > 0:
                lines.append(line.split(" : "))

    elif filepath.lower().endswith('.xls') or filepath.lower().endswith('.xlsx'):
        workbook = xlreader.open_workbook(filepath)
        worksheet = workbook.sheet_by_index(0)
        for row in range(worksheet.nrows):
            temp = list()
            for col in range(worksheet.ncols):
                try:
                    val = str(worksheet.cell(row, col).value).strip()
                    if len(val) > 0:
                        temp.append(val)
                except ValueError:
                    pass
                except:
                    raise
            lines.append(temp)
    else:
        warn('Parameter file not of expected format: text or spreadsheet')
        return None

    if verbose:
        print("Finished reading the file")

    is_beps = False
    parm_dict = dict()

    prefix = 'File_'
    for fields in lines:
        # Ignore the parameters describing the GUI choices
        if prefix == 'Multi_':
            continue

        # Check if the line is a group header or parameter/value pair
        if len(fields) == 2:
            # Get the current name/value pair, and clean up the name
            name = fields[0].strip().replace('# ', 'num_').replace('#', 'num_').replace(' ', '_')
            value = fields[1]

            # Rename specific parameters
            if name == '1_mode':
                name = 'mode'
            if name == 'IO_rate':
                name = 'IO_rate_[Hz]'
                value = int(value.split()[0]) * 1E6
            if name == 'AO_range':
                name = 'AO_range_[V]'
                value = ' '.join(value.split()[:-1])
            if name == 'AO_amplifier':
                value = value.split()[0]
            if name == 'cycle_time[s]':
                name = 'cycle_time_[s]'
            if name == 'FORC_V_high1[V]':
                name = 'FORC_V_high1_[V]'
            if name == 'FORC_V_high2[V]':
                name = 'FORC_V_high2_[V]'
            if name == 'FORC_V_low1[V]':
                name = 'FORC_V_low1_[V]'
            if name == 'FORC_V_low2[V]':
                name = 'FORC_V_low2_[V]'
            if name == 'amplitude[V]':
                name = 'amplitude_[V]'
            if name == 'offset[V]':
                name = 'offset_[V]'
            if name == 'set_pulst_amplitude[V]':
                name = 'set_pulst_amplitude_[V]'
            if name == 'set_pulst_duration[s]':
                name = 'set_pulst_duration_[s]'
            if name == 'step_edge_smoothing[s]':
                name = 'step_edge_smoothing_[s]'

            # Append the prefix to the name
            name = prefix + name.lstrip(prefix)

            # Write parameter to parm_dict
            try:
                num = float(value)
                parm_dict[name] = num
                if num == int(num):
                    parm_dict[name] = int(num)
            except ValueError:
                parm_dict[name] = value
            except:
                raise
        elif len(fields) == 1:
            # Change the parameter prefix to the new one from the group header
            prefix = fields[0].strip('<').strip('>')
            prefix = prefix.split()[0] + '_'

            # Check if there are VS Parameters.  Set isBEPS to true if so.
            is_beps = is_beps or prefix == 'VS_'

    if is_beps:
        useless_parms = \
            ['Multi_1_BE_response_spectra',
             'Multi_2_BE_amplitude_spectra',
             'Multi_3_BE_phase_spectra',
             'Multi_4_BE_response_spectrogram',
             'Multi_5_BE_amplitude_spectrogram',
             'Multi_6_BE_phase_spectrogram',
             'Multi_7_VS_response_loops',
             'Multi_8_VS_amplitude_loops',
             'Multi_9_VS_phase_loops',
             'Multi_10_topography',
             'Multi_11_mean_channel_2']
    else:
        useless_parms = \
            ['Multi_1_BE_response_spectra',
             'Multi_2_BE_amplitude_spectra',
             'Multi_3_BE_phase_spectra',
             'Multi_4_BE_response_spectrogram',
             'Multi_5_BE_amplitude_spectrogram',
             'Multi_6_BE_phase_spectrogram',
             'Multi_7_amplitude_map',
             'Multi_8_resonance_map',
             'Multi_9_Q_map',
             'Multi_10_phase_map',
             'Multi_11_topography',
             'Multi_12_AI_time_domain',
             'Multi_13_AI_Fourier_amplitude']

    if verbose:
        print("Finished parsing the text pairs. isBEPS = {}".format(is_beps))

    useless_parms.extend(parms_to_remove)

    # Now remove the list of useless parameters:
    for uparm in useless_parms:
        try:
            del parm_dict[uparm]
        except KeyError:
            # warn('Parameter to be deleted does not exist')
            pass
        except:
            raise
    del uparm, useless_parms

    if verbose:
        print("Finished removing useless parameters")

    if is_beps:
        # fix the DC type in the parms:
        if parm_dict['VS_measure_in_field_loops'] == 'out-of-field only':
            parm_dict['VS_measure_in_field_loops'] = 'out-of-field'
        elif parm_dict['VS_measure_in_field_loops'] == 'in-field only':
            parm_dict['VS_measure_in_field_loops'] = 'in-field'

    return is_beps, parm_dict


###############################################################################


def requires_conjugate(chosen_spectra, default_q=10, cores=None):
    """
    Determines whether or not the conjugate of the data needs to be taken based on the quality factor

    Parameters
    ----------
    chosen_spectra : 2D complex numpy array
        N random spectra arranged as [instance, frequency]
    default_q : unsigned int, Optional
        Default value of Q factor that the SHO guess function results in for poor guesses
    cores : uint, Optional
        Number of CPU cores to use for the SHO guess.
        Default = None - ask pyUSID.processing.comp_utils.recommend_cpu_cores for automatic assignment

    Returns
    -------
    do_conjugate : Boolean
        Whether or not to take the conjugate of the data
    """
    # Do the SHO Guess for each of these
    sho_guess = parallel_compute(chosen_spectra, SHOestimateGuess, cores=cores,
                                 func_args=[np.arange(chosen_spectra.shape[1]),
                                            5])

    q_results = np.array(sho_guess)[:, 2]
    good_q = q_results[np.where(q_results != default_q)]

    if np.mean(good_q) < 0:
        return True
    else:
        return False


###############################################################################

def getSpectroscopicParmLabel(expt_type):
    """
    Returns the label for the spectroscopic parameter in the plot group. 
    
    Parameters
    ----------
    expt_type : str
        Type of the experiment - found in the parms.txt file
    
    Returns
    -------
    str
        label for the spectroscopic parameter axis in the plot
    """

    if expt_type in ['DC modulation mode', 'current mode']:
        return 'DC Bias'
    elif expt_type == 'AC modulation mode with time reversal':
        return 'AC amplitude'
    return 'User Defined'


def normalizeBEresponse(spectrogram_mat, FFT_BE_wave, harmonic):
    """
    This function normalizes the BE waveform to correct the phase by diving by the excitation
    
    Parameters
    ------------
    spectrogram_mat : 2D complex numpy array
        BE response arranged as [bins, steps]
    FFT_BE_wave : 1D complex numpy array
        FFT of the BE waveform at the appropriate bins. Number of bins must match with spectrogram_mat
    harmonic : unsigned int
        nth harmonic of the excitation waveform
        
    Returns
    ----------
    spectrogram_mat : 2D complex numpy array
        Normalized BE response spectrogram

    """
    BE_wave = np.fft.ifftshift(np.fft.ifft(FFT_BE_wave))
    scaling_factor = 1

    if harmonic == 2:
        scaling_factor = np.fft.fftshift(np.fft.fft(BE_wave ** 2)) / (2 * np.exp(1j * 3 * np.pi * 0.5))
    elif harmonic == 3:
        scaling_factor = np.fft.fftshift(np.fft.fft(BE_wave ** 3)) / (4 * np.exp(1j * np.pi))
    elif harmonic >= 4:
        print("Warning these high harmonics are not supported in translator.")

    # Generate transfer functions
    F_AO_spectrogram = np.transpose(np.tile(FFT_BE_wave / scaling_factor, [spectrogram_mat.shape[1], 1]))
    # Divide by transfer function
    spectrogram_mat = spectrogram_mat / F_AO_spectrogram

    return spectrogram_mat


def generatePlotGroups(h5_main, mean_resp, folder_path, basename, max_resp=[], min_resp=[],
                       max_mem_mb=1024, spec_label='None',
                       show_plots=True, save_plots=True, do_histogram=False,
                       debug=False):
    """
    Generates the spatially averaged datasets for the given raw dataset. 
    The averaged datasets are necessary for quick visualization of the quality of data. 
    
    Parameters
    ----------
    h5_main : H5 reference
        to the main dataset
    mean_resp : 1D numpy array
        spatially averaged amplitude
    folder_path : String
        Absolute path of the data folder
    basename : String
        base name of the dataset
    max_resp : 1D numpy array
        Maximum amplitude for all pixels
    min_resp : 1D numpy array
        Minimum amplitude for all pixels
    max_mem_mb : Unisigned integer
        Maximum memory that can be used for generating histograms
    spec_label : String
        Parameter that is varying
    show_plots : (optional) Boolean
        Whether or not to show plots
    save_plots : (optional) Boolean
        Whether or not to save generated plots
    do_histogram : Boolean (Optional. Default = False)
        Whether or not to generate hisograms. 
        Caution - Histograms can take a fair amount of time to compute.
    debug : Boolean, Optional
        If True, then extra debug statements are printed.
        Default False
    """
    # Too
    assert isinstance(h5_main, h5py.Dataset)
    h5_f = h5_main.file

    grp = h5_main.parent
    h5_freq = grp['Bin_Frequencies']
    UDVS = grp['UDVS']
    spec_inds = h5_main.h5_spec_inds
    UDVS_inds = grp['UDVS_Indices']
    spec_vals = h5_main.h5_spec_vals

    #     std_cols = ['wave_type','Frequency','DC_Offset','wave_mod','AC_Amplitude','dc_offset','ac_amplitude']

    col_names = UDVS.attrs['labels']

    if len(col_names) <= 5:
        """
        No plot groups are defined in the UDVS table.
        All plot group datasets will be written to the 
        Channel group
        """
        return

    # Removing the standard columns
    col_names = get_attr(UDVS, 'labels')[5:]

    #     col_names = [col for col in col_names if col not in std_cols + ignore_plot_groups]

    freq_inds = spec_inds[spec_inds.attrs['Frequency']].flatten()

    for col_name in col_names:

        # Make sure the column name is a regular Python string
        col_name = str(col_name)

        ref = UDVS.attrs[col_name]
        #         Make sure we're actually dealing with a reference of some type
        if not isinstance(ref, h5py.RegionReference):
            continue
        # 4. Access that column of the data through region reference
        steps = np.where(np.isfinite(UDVS[ref]))[0]
        step_inds = np.array([np.where(UDVS_inds[()] == step)[0] for step in steps]).flatten()
        """selected_UDVS_steps = UDVS[ref]
        selected_UDVS_steps = selected_UDVS_steps[np.isfinite(selected_UDVS_steps)]"""

        (step_averaged_vec, mean_spec) = reshape_mean_data(spec_inds, step_inds, mean_resp)

        """ 
        Need to account for cases with multiple excitation waveforms
        This will affect the frequency indices / values
        We are assuming that there is only one excitation waveform per plot group
        """
        freq_slice = np.unique(freq_inds[step_inds])
        freq_vec = h5_freq[()][freq_slice]

        num_bins = len(freq_slice)  # int(len(freq_inds)/len(UDVS[ref]))
        pg_data = np.repeat(UDVS[ref], num_bins)

        plot_grp = create_indexed_group(grp, 'Spatially_Averaged_Plot_Group')
        write_simple_attrs(plot_grp, {'Name': col_name})

        h5_mean_spec = plot_grp.create_dataset('Mean_Spectrogram',
                                               data=mean_spec,
                                               dtype=np.complex64)
        h5_step_avg = plot_grp.create_dataset('Step_Averaged_Response',
                                              data=step_averaged_vec,
                                              dtype=np.complex64)
        # cannot assume that this is DC offset, could be AC amplitude....
        h5_spec_parm = plot_grp.create_dataset('Spectroscopic_Parameter',
                                               data=np.squeeze(pg_data[step_inds]),
                                               dtype=np.uint32)
        write_simple_attrs(h5_spec_parm, {'name': spec_label})
        h5_freq_vec = plot_grp.create_dataset('Bin_Frequencies',
                                          data=freq_vec,
                                          dtype=h5_freq.dtype)

        # Linking the datasets with the frequency and the spectroscopic variable:
        link_h5_objects_as_attrs(h5_mean_spec, [h5_spec_parm, h5_freq_vec])
        link_h5_objects_as_attrs(h5_step_avg, [h5_freq_vec])

        """
        Create Region Reference for the plot group in the Raw_Data, Spectroscopic_Indices 
        and Spectroscopic_Values Datasets
        """
        raw_ref = h5_main.regionref[:, step_inds]
        spec_inds_ref = spec_inds.regionref[:, step_inds]
        spec_vals_ref = spec_vals.regionref[:, step_inds]

        ref_name = col_name.replace(' ', '_').replace('-', '_') + '_Plot_Group'
        h5_main.attrs[ref_name] = raw_ref
        spec_inds.attrs[ref_name] = spec_inds_ref
        spec_vals.attrs[ref_name] = spec_vals_ref

        h5_f.flush()

        if do_histogram:
            """
            Build the histograms for the current plot group
            """
            hist = BEHistogram()
            hist_mat, hist_labels, hist_indices, hist_indices_labels = \
                hist.buildPlotGroupHist(h5_main, step_inds, max_response=max_resp,
                                        min_response=min_resp, max_mem_mb=max_mem_mb, debug=debug)

            hist_grp = create_indexed_group(plot_grp, 'Histogram')

            hist_spec_dims = list()
            hist_units = ['V', '', 'V', 'V']
            for hist_ind, hist_dim in enumerate(hist_labels):
                hist_spec_dims.append(Dimension(hist_dim,
                                                hist_units[hist_ind],
                                                hist_indices[hist_ind]))

            h5_hist = write_main_dataset(hist_grp, hist_mat, 'Histograms',
                                         'Counts', 'a.u.',
                                         None, hist_spec_dims,
                                         h5_pos_inds=h5_main.h5_pos_inds, h5_pos_vals=h5_main.h5_pos_vals,
                                         dtype=np.int32,
                                         chunking=(1, hist_mat.shape[1]),
                                         compression='gzip')

        else:
            """
            Write the min and max response vectors so that histograms can be generated later.
            """
            h5_max_resp = plot_grp.create_dataset('Max_Response',
                                                  data=max_resp)
            h5_min_resp = plot_grp.create_dataset('Min_Response',
                                                  data=min_resp)

        if save_plots or show_plots:
            fig_title = '_'.join(grp.name[1:].split('/') + [col_name])

            fig_1d, axes_1d = plot_1d_spectrum(step_averaged_vec, freq_vec, fig_title)
            if save_plots:
                path_1d = path.join(folder_path, basename + '_Step_Avg_' + fig_title + '.png')
                path_2d = path.join(folder_path, basename + '_Mean_Spec_' + fig_title + '.png')
                path_hist = path.join(folder_path, basename + '_Histograms_' + fig_title + '.png')

                fig_1d.savefig(path_1d, format='png', dpi=300)
            if mean_spec.shape[0] > 1:
                fig_2d, axes_2d = plot_2d_spectrogram(mean_spec, freq_vec, title=fig_title)
                if save_plots:
                    fig_2d.savefig(path_2d, format='png', dpi=300)

            if do_histogram:
                plot_histograms(hist_mat, hist_indices, grp.name, figure_path=path_hist)

            if show_plots:
                plt.show()
            else:
                plt.close('all')

            # print('Generated spatially average data for group: %s' %(col_name))
    print('Completed generating spatially averaged plot groups')


###############################################################################


def reshape_mean_data(spec_inds, step_inds, mean_resp):
    """
    Takes in the mean data vector and rearranges that data according to 
    plot group as [step number,bins]
    
    Parameters
    -----------
    spec_inds : 2D numpy array
        UDVS_Indices as a 2D mat
    step_inds : 1D numpy array or list
        UDVS step indices corresponding to this plot group
    mean_resp : 1D numpy complex array
        position averaged BE data
            
    Returns
    ----------
    step_averaged_vec : 1D complex numpy array
        Mean (position averaged) spectrogram averaged over the UDVS steps as well
    mean_spectrogram : 2D complex numpy array
        Position averaged data arranged as [step number,bins]
    """
    num_bins = len(np.unique(spec_inds[0, step_inds]))
    # Stephen says that we can assume that the number of bins will NOT change in a plot group
    mean_spectrogram = mean_resp[step_inds].reshape(-1, num_bins)

    step_averaged_vec = np.mean(mean_spectrogram, axis=0)
    return step_averaged_vec, mean_spectrogram


###############################################################################


def visualize_plot_groups(h5_filepath):
    """
    Visualizes the plot groups present in the provided BE data file
    
    Parameters
    ----------
    h5_filepath : String / Uniciode
        Absolute path of the h5 file
    """
    with h5py.File(h5_filepath, mode='r') as h5f:
        expt_type = h5f.attrs.get('data_type')
        if expt_type not in ['BEPSData', 'BELineData']:
            warn('Invalid data format')
            return
        for grp_name in h5f.keys():
            grp = h5f[grp_name]['Channel_000']
            for plt_grp_name in grp.keys():
                if plt_grp_name.startswith('Spatially_Averaged_Plot_Group_'):
                    plt_grp = grp[plt_grp_name]
                    if expt_type == 'BEPSData':
                        spect_data = plt_grp['Mean_Spectrogram'].value
                        _ = plot_2d_spectrogram(spect_data, plt_grp['Bin_Frequencies'].value,
                                                title=plt_grp.attrs['Name'])
                    step_avg_data = plt_grp['Step_Averaged_Response']
                    _ = plot_1d_spectrum(step_avg_data, plt_grp['Bin_Frequencies'].value, plt_grp.attrs['Name'])
                    try:
                        hist_data = plt_grp['Histograms']
                        hist_bins = plt_grp['Histograms_Indicies']
                        plot_histograms(hist_data, hist_bins, plt_grp.attrs['Name'])
                    except:
                        pass

    plt.show()
    plt.close('all')


###############################################################################


def trimUDVS(udvs_mat, udvs_labs, udvs_units, target_col_names):
    """
    Removes unused (typically default) plot groups
    
    Parameters
    ----------
    udvs_mat : 2D numpy array
        UDVS table arranged as [steps, col]
    udvs_labs : list of strings
        Column names of the UDVS table
    udvs_units : list of strings
        Units for the columns of the UDVS table
    target_col_names : list of strings
        Column names that need to be removed
        
    Returns
    -----------
    udvs_mat : 2D numpy array
        Truncated UDVS table
    udvs_labs : list of strings 
        Truncated list of UDVS column names
    udvs_units : list of strings
        Truncated list of UDVS column units
    """

    if len(target_col_names) == 0:
        return udvs_mat, udvs_labs, udvs_units

    if len(udvs_labs) != udvs_mat.shape[1]:
        warn('Error: Incompatible UDVS matrix and labels. Not truncating!')
        return udvs_mat, udvs_labs, udvs_units

    # First figure out the column indices
    col_inds = []
    found_cols = []
    for target in target_col_names:
        for ind, udvs_col in enumerate(udvs_labs):
            if udvs_col == target:
                col_inds.append(ind)
                found_cols.append(udvs_col)
                break
    col_inds = np.sort(np.unique(col_inds))

    # Now remove from the labels and the matrix
    udvs_mat = np.delete(udvs_mat, col_inds, axis=1)
    # col_inds.sort(reverse=True)
    [udvs_units.pop(ind) for ind in range(len(col_inds), 0, -1)]
    udvs_labs = [col for col in udvs_labs if col not in found_cols]

    return udvs_mat, udvs_labs, udvs_units


def createSpecVals(udvs_mat, spec_inds, bin_freqs, bin_wfm_type, parm_dict,
                   udvs_labs, udvs_units, verbose=False):
    """
    This function will determine the proper Spectroscopic Value array for the 
    dataset
    
    Chris Smith -- csmith55@utk.edu
    
    Parameters
    ----------
    udvs_mat : numpy array 
        UDVS table from dataset
    spec_inds : numpy array 
        Spectroscopic Indices table from dataset
    bin_freqs : numpy array 
        Bin frequencies
    bin_wfm_type : numpy array 
        waveform type for each frequency index
    parm_dict : dictionary
        parameters for dataset
    udvs_labs : list of strings
        labels for the columns of the UDVS matrix
    udvs_units : list of strings
        units for the columns of the UDVS matrix
    verbose : bool, optional. Default = False
        Whether or not to print debugging print statements
    
    Returns
    -----------
    ds_spec_val_mat : numpy array 
        Spectroscopic Values table
    ds_spec_val_labs : list of strings
        names of the columns of the Spectroscopic Values table
    ds_spec_val_units : list of strings
        units of the columns of the Spectroscopic Values table
    ds_spec_val_labs_names : list of Strings
        labels with the names of their parameters
    """

    def __FindSpecValIndices(udvs_mat, spec_inds, usr_defined=False,
                             verbose=False):
        """
        This function finds the Spectroscopic Values associated with the dataset that
        have more than one unique value
        
        Parameters
        ----------
        udvs_mat : numpy array containing the UDVS table
        spec_inds : numpy array contain Spectroscopic indices table
    
        Returns
        -------
        iSpec_var : integer array holding column indices in UDVS that change
        ds_spec_val_mat : array holding all spectral values for columns in iSpec_var
        """
        udvs_cols = ['step_num', 'dc_offset', 'ac_ampli', 'wave_type',
                     'wave_mod', 'in-field', 'out-of-field']
        if verbose:
            print('\t' * 3 + '__FindSpecValIndices:')
            print('\t' * 4 + 'UDVS matrix of shape: {}'.format(udvs_mat.shape))
            print('\t' * 4 + 'spec_inds of shape: {}'.format(spec_inds.shape))
            if False:
                print('\t'.join(udvs_cols))
                for ud_row in udvs_mat:
                    print('\t\t'.join(['{:04.2f}'.format(item) for item in ud_row]))
            if False:  # Turn this off if necessary
                fig, axes = plt.subplots(nrows=2, figsize=(10, 5))
                for ind, axis in enumerate(axes.flat):
                    axis.plot(spec_inds[ind, :])
                    axis.set_title('spec_inds[{}]'.format(ind))
                fig.tight_layout()

        """
        icheck is an array containing all UDVS steps which should be checked.
        """
        icheck = np.unique(spec_inds[1])
        if verbose and False:
            print('\t' * 4 + 'UDVS steps that will be checked: {}'.format(icheck))
        if len(icheck) < 1:
            raise ValueError('No row in the spectroscopic indices varied.\n'
                             'Cannot build spectroscopic datasets further')

        # Copy even step values of DC_offset into odd steps
        UDVS = np.copy(udvs_mat)

        if not usr_defined:
            DC = UDVS[:, 1]
            for step in range(0, DC.size, 2):
                DC[step + 1] = DC[step]
            UDVS[:, 1] = DC

        """
        Keep only the UDVS values for steps which we care about and the 
        first 5 columns
        """
        UDVS = UDVS[icheck, :5]
        udvs_cols = udvs_cols[:5]
        #         UDVS = np.array([UDVS[i] for i in icheck])
        if verbose and False:
            print('\t' * 4 + 'UDVS matrix after down-selecting rows: {}'
                             ''.format(UDVS.shape))
            print('\t'.join(udvs_cols))
            for ud_row in UDVS:
                print('\t\t'.join(['{:04.2f}'.format(item) for item in ud_row]))

        """
        Transpose UDVS for ease of looping later on and store the number of steps
        as num_cols
        """
        num_cols = np.size(UDVS, 1)
        """
        Initialize the iSpec_var as an empty array.  It will store the index of the 
        UDVS label for any column which has more than one unique value
        """
        iSpec_var = []

        """
        Loop over all columns in udvs_mat
        """
        for col_ind, col_name in zip(range(1, num_cols), udvs_cols[1:]):
            """
            Find all unique values in the current column
            """
            toosmall = np.where(abs(UDVS[:, col_ind]) < 1E-5)[0]
            UDVS[toosmall, col_ind] = 0
            uvals = np.unique(UDVS[:, col_ind])
            """
            np.unique considers all NaNs to be unique values
            These two lines find the indices of all NaNs in the unique value array 
            and removes all but the first
            """
            nanvals = np.where(np.isnan(uvals))[0]
            if verbose:
                print('\t\t\tColumn {}: {} has {} unique values and {} NaNs'
                      ''.format(col_ind, col_name, len(uvals), len(nanvals)))
            uvals = np.delete(uvals, nanvals[1:])
            """
            Check if more that one unique value
            Append column number to iSpec_var if true
            """
            if verbose:
                print('\t\t\tColumn {}: {} had {} actually unique values'
                      ''.format(col_ind, col_name, len(uvals)))
            if uvals.size > 1:
                iSpec_var = np.append(iSpec_var, int(col_ind))

        if verbose:
            print('\t' * 4 + 'UDVS matrix of shape: {}'.format(UDVS.shape))
            print('\t' * 4 + 'Taking columns specified by iSpec_var: {}'
                             ''.format(iSpec_var))

        if len(iSpec_var) < 1:
            raise ValueError('Invalid UDVS inputs. No variables were varied')

        iSpec_var = np.asarray(iSpec_var, np.int)
        ds_spec_val_mat = UDVS[:, iSpec_var]

        return iSpec_var, ds_spec_val_mat

    def __BEPSVals(udvs_mat, spec_inds, bin_freqs, bin_wfm_type, parm_dict,
                   udvs_labs, udvs_units, verbose=False):
        """
        Returns the Spectroscopic Value array for a BEPS dataset
    
        Parameters
        ----------
        udvs_mat : hdf5 dataset reference to UDVS dataset
        spec_inds : numpy array containing Spectroscopic indices table 
        bin_freqs : 1D numpy array of frequencies
        bin_wfm_type : numpy array containing the waveform type for each
            frequency index
        parm_dict : parameter dictinary for dataset
        udvs_labs : list of labels for the columns of the UDVS matrix
        udvs_units : list of units for the columns of the UDVS matrix
        
        Returns
        -------
        ds_spec_val_mat : list holding final Spectroscopic Value table 
        """

        """
        Check the mode to determine if DC, AC, or something else
        """
        mode = parm_dict['VS_mode']

        if mode == 'DC modulation mode' or mode == 'current mode':
            """ 
            First we call the FindSpecVals function to get the columns in UDVS 
            of interest and return a first pass on the spectral value array 
            """
            if verbose:
                print('\t' * 3 + 'DC modulation mode / current mode')
            iSpecVals, inSpecVals = __FindSpecValIndices(udvs_mat, spec_inds,
                                                         verbose=verbose)
            if verbose:
                print('\t' * 3 + 'Generated spec vals. Calling __BEPSDC')

            return __BEPSDC(udvs_mat, inSpecVals, bin_freqs, bin_wfm_type,
                            parm_dict, verbose=verbose)

        elif mode == 'AC modulation mode with time reversal':
            """ 
            First we call the FindSpecVals function to get the columns in UDVS
            of interest and return a first pass on the spectral value array 
            """
            if verbose:
                print('\t' * 3 + 'AC modulation')
            iSpecVals, inSpecVals = __FindSpecValIndices(udvs_mat, spec_inds,
                                                         verbose=verbose)
            if verbose:
                print('\t' * 3 + 'Generated spec vals. Calling __BEPSAC')

            return __BEPSAC(udvs_mat, inSpecVals, bin_freqs, bin_wfm_type,
                            parm_dict, verbose=verbose)
        else:
            """ 
            First we call the FindSpecVals function to get the columns in UDVS
            of interest and return a first pass on the spectral value array 
            """
            if verbose:
                print('\t' * 3 + 'user defined voltage spectroscopy')
            iSpecVals, inSpecVals = __FindSpecValIndices(udvs_mat, spec_inds,
                                                         usr_defined=True,
                                                         verbose=verbose)
            if verbose:
                print('\t' * 3 + 'Generated spec vals. Calling __BEPSgen')

            return __BEPSgen(udvs_mat, inSpecVals, bin_freqs, bin_wfm_type,
                             parm_dict, udvs_labs, iSpecVals, udvs_units)

    def __BEPSDC(udvs_mat, inSpecVals, bin_freqs, bin_wfm_type, parm_dict,
                 verbose=False):
        """
        Calculates Spectroscopic Values for BEPS data in DC modulation mode
        
        Parameters
        ----------
        udvs_mat : hdf5 dataset reference to UDVS dataset
        inSpecVals : list holding initial guess at spectral values 
        bin_freqs : 1D numpy array of frequencies
        bin_wfm_type : numpy array containing the waveform type for each frequency index
        parm_dict : parameter dictinary for dataset
        verbose: bool, optional. Default = False
            Whether or not to print statements for debugging
                        
        Returns
        -------
        ds_spec_val_mat : list holding final Spectroscopic Value table 
        SpecValsLabels : list holding labels of column names of ds_spec_val_mat
        """
        hascycles = False
        hasFORCS = False

        # print('in shape:',np.shape(inSpecVals))
        """
        All DC datasets will need Spectroscopic Value fields for Bin, DC, and Field
        
        Bin     - Bin number
        DC      - dc offset from UDVS
        Field   - 0 if in-field, 1 if out-of-field 
        """
        nrow = 2
        ds_spec_val_labs = ['Frequency', 'DC_Offset']
        ds_spec_val_units = ['Hz', 'V']

        """
        Get the number of and steps in the current dataset
        """
        numsteps = udvs_mat.shape[0]

        """
        Get the wave form for each step from udvs_mat
        """
        wave_form = udvs_mat[:, 3]

        """
        Define list of attribute names needed from the group metadata
        """
        field_type = parm_dict['VS_measure_in_field_loops']
        numcycles = parm_dict['VS_number_of_cycles']
        numFORCs = parm_dict['FORC_num_of_FORC_cycles']
        numcyclesteps = parm_dict['VS_steps_per_full_cycle']
        cycle_fraction = parm_dict['VS_cycle_fraction']

        if field_type == 'in and out-of-field':
            ds_spec_val_labs.append('Field')
            ds_spec_val_units.append('')
            nrow += 1

        frac = {'full': 1.0, '1/2': 0.5, '1/4': 0.25, '3/4': 0.75}

        numcyclesteps = frac[cycle_fraction] * numcyclesteps

        """
        Check the number of cycles and FORCs
        Add to Spectroscopic Values Labels as needed 
        """
        if numcycles > 1:
            hascycles = True
            nrow += 1
            ds_spec_val_labs.append('Cycle')
            ds_spec_val_units.append('')

        if numFORCs > 1:
            hasFORCS = True
            """
        It's possible to have 1 cycle with multiple FORCs so force cycle tracking 
        if FROCs exist
            """
            nrow += 1
            ds_spec_val_labs.append('FORC')
            ds_spec_val_units.append('')
            numFORCsteps = numcycles * numcyclesteps * 2

        """
        Check the field type
        For in-field and out-of-field we know all values of field ahead of time
        For in and out-of-field we must check at each step
        If something else is in field_type, we default to in and out-of-field and print message
        """
        if field_type == 'out-of-field':
            field = 1
            numsteps = int(numsteps / 2)
            numcyclesteps = int(numcyclesteps / 2)
            swapfield = [1, 1]
            field_names = ['out-of-field']
        elif field_type == 'in-field':
            field = 0
            numsteps = int(numsteps / 2)
            numcyclesteps = int(numcyclesteps / 2)
            swapfield = [0, 0]
            field_names = ['in-field']
        elif field_type == 'in and out-of-field':
            field = 0
            swapfield = [1, 0]
            field_names = ['out-of-field', 'in-field']
        else:
            warn('{} is not a known field type'.format(field_type))
            field = 0