<commit_before>"""This class provides often needed analysis functions, for analysis that is done with python.
"""

import logging
import re
import os
import time
import collections
import numpy as np
import progressbar
import glob
import tables as tb
import numexpr as ne
from operator import itemgetter
from scipy.interpolate import interp1d
from scipy.interpolate import splrep, splev
from tables import dtype_from_descr
from stat import ST_CTIME

from pybar_fei4_interpreter import analysis_functions

from pybar.daq.fei4_record import FEI4Record
from pybar.analysis.plotting import plotting
from pybar.daq.readout_utils import is_fe_word, is_data_header


class AnalysisError(Exception):

    """Base class for exceptions in this module.
    """


class IncompleteInputError(AnalysisError):

    """Exception raised for errors in the input.
    """


class InvalidInputError(AnalysisError):

    """Exception raised for errors in the input.
    """


class NotSupportedError(AnalysisError):

    """Exception raised for not supported actions.
    """


def generate_threshold_mask(hist):
    '''Masking array elements when equal 0.0 or greater than 10 times the median

    Parameters
    ----------
    hist : array_like
        Input data.

    Returns
    -------
    masked array
        Returns copy of the array with masked elements.
    '''
    masked_array = np.ma.masked_values(hist, 0)
    masked_array = np.ma.masked_greater(masked_array, 10 * np.ma.median(hist))
    logging.info('Masking %d pixel(s)', np.ma.count_masked(masked_array))
    return np.ma.getmaskarray(masked_array)


def unique_row(array, use_columns=None, selected_columns_only=False):
    '''Takes a numpy array and returns the array reduced to unique rows. If columns are defined only these columns are taken to define a unique row.
    The returned array can have all columns of the original array or only the columns defined in use_columns.
    Parameters
    ----------
    array : numpy.ndarray
    use_columns : list
        Index of columns to be used to define a unique row
    selected_columns_only : bool
        If true only the columns defined in use_columns are returned

    Returns
    -------
    numpy.ndarray
    '''
    if array.dtype.names is None:  # normal array has no named dtype
        if use_columns is not None:
            a_cut = array[:, use_columns]
        else:
            a_cut = array
        if len(use_columns) > 1:
            b = np.ascontiguousarray(a_cut).view(np.dtype((np.void, a_cut.dtype.itemsize * a_cut.shape[1])))
        else:
            b = np.ascontiguousarray(a_cut)
        _, index = np.unique(b, return_index=True)
        if not selected_columns_only:
            return array[np.sort(index)]  # sort to preserve order
        else:
            return a_cut[np.sort(index)]  # sort to preserve order
    else:  # names for dtype founnd --> array is recarray
        names = list(array.dtype.names)
        if use_columns is not None:
            new_names = [names[i] for i in use_columns]
        else:
            new_names = names
        a_cut, index = np.unique(array[new_names], return_index=True)
        if not selected_columns_only:
            return array[np.sort(index)]  # sort to preserve order
        else:
            return array[np.sort(index)][new_names]  # sort to preserve order


def get_ranges_from_array(arr, append_last=True):
    '''Takes an array and calculates ranges [start, stop[. The last range end is none to keep the same length.

    Parameters
    ----------
    arr : array like
    append_last: bool
        If True, append item with a pair of last array item and None.

    Returns
    -------
    numpy.array
        The array formed by pairs of values by the given array.

    Example
    -------
    >>> a = np.array((1,2,3,4))
    >>> get_ranges_from_array(a, append_last=True)
    array([[1, 2],
           [2, 3],
           [3, 4],
           [4, None]])
    >>> get_ranges_from_array(a, append_last=False)
    array([[1, 2],
           [2, 3],
           [3, 4]])
    '''
    right = arr[1:]
    if append_last:
        left = arr[:]
        right = np.append(right, None)
    else:
        left = arr[:-1]
    return np.column_stack((left, right))


def get_mean_from_histogram(counts, bin_positions):
    return np.dot(counts, np.array(bin_positions)) / np.sum(counts).astype('f4')


def get_median_from_histogram(counts, bin_positions):
    values = []
    for index, one_bin in enumerate(counts):
        for _ in range(one_bin):
            values.append(bin_positions[index])
    return np.median(values)


def get_rms_from_histogram(counts, bin_positions):
    values = []
    for index, one_bin in enumerate(counts):
        for _ in range(one_bin):
            values.append(bin_positions[index])
    return np.std(values)


def in1d_sorted(ar1, ar2):
    """
    Does the same than np.in1d but uses the fact that ar1 and ar2 are sorted. Is therefore much faster.

    """
    if ar1.shape[0] == 0 or ar2.shape[0] == 0:  # check for empty arrays to avoid crash
        return []
    inds = ar2.searchsorted(ar1)
    inds[inds == len(ar2)] = 0
    return ar2[inds] == ar1


def central_difference(x, y):
    '''Returns the dy/dx(x) via central difference method

    Parameters
    ----------
    x : array like
    y : array like

    Returns
    -------
    dy/dx : array like
    '''
    if (len(x) != len(y)):
        raise ValueError("x, y must have the same length")
    z1 = np.hstack((y[0], y[:-1]))
    z2 = np.hstack((y[1:], y[-1]))
    dx1 = np.hstack((0, np.diff(x)))
    dx2 = np.hstack((np.diff(x), 0))
    return (z2 - z1) / (dx2 + dx1)


def get_profile_histogram(x, y, n_bins=100):
    '''Takes 2D point data (x,y) and creates a profile histogram similar to the TProfile in ROOT. It calculates
    the y mean for every bin at the bin center and gives the y mean error as error bars.

    Parameters
    ----------
    x : array like
        data x positions
    y : array like
        data y positions
    n_bins : int
        the number of bins used to create the histogram
    '''
    if len(x) != len(y):
        raise ValueError('x and y dimensions have to be the same')
    y = y.astype(np.float32)
    n, bin_edges = np.histogram(x, bins=n_bins)  # needed to calculate the number of points per bin
    sy = np.histogram(x, bins=n_bins, weights=y)[0]  # the sum of the bin values
    sy2 = np.histogram(x, bins=n_bins, weights=y * y)[0]  # the quadratic sum of the bin values
    bin_centers = (bin_edges[1:] + bin_edges[:-1]) / 2  # calculate the bin center for all bins
    mean = sy / n  # calculate the mean of all bins
    std = np.sqrt((sy2 / n - mean * mean))  # TODO: not understood, need check if this is really the standard deviation
    std_mean = std / np.sqrt((n - 1))
    mean[np.isnan(mean)] = 0.
    std_mean[np.isnan(std_mean)] = 0.
    return bin_centers, mean, std_mean


def get_rate_normalization(hit_file, parameter, reference='event', cluster_file=None, plot=False, chunk_size=500000):
    ''' Takes different hit files (hit_files), extracts the number of events or the scan time (reference) per scan parameter (parameter)
    and returns an array with a normalization factor. This normalization factor has the length of the number of different parameters.
    If a cluster_file is specified also the number of cluster per event are used to create the normalization factor.

    Parameters
    ----------
    hit_files : string
    parameter : string
    reference : string
    plot : bool

    Returns
    -------
    numpy.ndarray
    '''

    logging.info('Calculate the rate normalization')
    with tb.openFile(hit_file, mode="r+") as in_hit_file_h5:  # open the hit file
        meta_data = in_hit_file_h5.root.meta_data[:]
        scan_parameter = get_scan_parameter(meta_data)[parameter]
        event_numbers = get_meta_data_at_scan_parameter(meta_data, parameter)['event_number']  # get the event numbers in meta_data where the scan parameter changes
        event_range = get_ranges_from_array(event_numbers)
        normalization_rate = []
        normalization_multiplicity = []
        try:
            event_range[-1, 1] = in_hit_file_h5.root.Hits[-1]['event_number'] + 1
        except tb.NoSuchNodeError:
            logging.error('Cannot find hits table')
            return

        # calculate rate normalization from the event rate for triggered data / measurement time for self triggered data for each scan parameter
        if reference == 'event':
            n_events = event_range[:, 1] - event_range[:, 0]  # number of events for every parameter setting
            normalization_rate.extend(n_events)
        elif reference == 'time':
            time_start = get_meta_data_at_scan_parameter(meta_data, parameter)['timestamp_start']
            time_spend = np.diff(time_start)
            time_spend = np.append(time_spend, meta_data[-1]['timestamp_stop'] - time_start[-1])  # TODO: needs check, add last missing entry
            normalization_rate.extend(time_spend)
        else:
            raise NotImplementedError('The normalization reference ' + reference + ' is not implemented')

        if cluster_file:  # calculate the rate normalization from the mean number of hits per event per scan parameter, needed for beam data since a beam since the multiplicity is rarely constant
            cluster_table = in_hit_file_h5.root.Cluster
            index_event_number(cluster_table)
            index = 0  # index where to start the read out, 0 at the beginning, increased during looping, variable for read speed up
            best_chunk_size = chunk_size  # variable for read speed up
            total_cluster = 0
            progress_bar = progressbar.ProgressBar(widgets=['', progressbar.Percentage(), ' ', progressbar.Bar(marker='*', left='|', right='|'), ' ', ETA(smoothing=0.8)], maxval=cluster_table.shape[0], term_width=80)
            progress_bar.start()
            for start_event, stop_event in event_range:  # loop over the selected events
                readout_cluster_len = 0  # variable to calculate a optimal chunk size value from the number of hits for speed up
                n_cluster_per_event = None
                for clusters, index in data_aligned_at_events(cluster_table, start_event_number=start_event, stop_event_number=stop_event, start=index, chunk_size=best_chunk_size):
                    if n_cluster_per_event is None:
                        n_cluster_per_event = get_n_cluster_in_events(clusters['event_number'])[:, 1]  # array with the number of cluster per event, cluster per event are at least 1
                    else:
                        n_cluster_per_event = np.append(n_cluster_per_event, get_n_cluster_in_events(clusters['event_number'])[:, 1])
                    readout_cluster_len += clusters.shape[0]
                    total_cluster += clusters.shape[0]
                    progress_bar.update(index)
                best_chunk_size = int(1.5 * readout_cluster_len) if int(1.05 * readout_cluster_len) < chunk_size else chunk_size  # to increase the readout speed, estimated the number of hits for one read instruction
                normalization_multiplicity.append(np.mean(n_cluster_per_event))
            progress_bar.finish()
            if total_cluster != cluster_table.shape[0]:
                logging.warning('Analysis shows inconsistent number of cluster (%d != %d). Check needed!', total_cluster, cluster_table.shape[0])

    if plot:
        x = scan_parameter
        if reference == 'event':
            plotting.plot_scatter(x, normalization_rate, title='Events per ' + parameter + ' setting', x_label=parameter, y_label='# events', log_x=True, filename=hit_file[:-3] + '_n_event_normalization.pdf')
        elif reference == 'time':
            plotting.plot_scatter(x, normalization_rate, title='Measuring time per GDAC setting', x_label=parameter, y_label='time [s]', log_x=True, filename=hit_file[:-3] + '_time_normalization.pdf')
        if cluster_file:
            plotting.plot_scatter(x, normalization_multiplicity, title='Mean number of particles per event', x_label=parameter, y_label='number of hits per event', log_x=True, filename=hit_file[:-3] + '_n_particles_normalization.pdf')
    print len(normalization_rate), len(normalization_multiplicity)
    if cluster_file:
        normalization_rate = np.array(normalization_rate)
        normalization_multiplicity = np.array(normalization_multiplicity)
        return np.amax(normalization_rate * normalization_multiplicity).astype('f16') / (normalization_rate * normalization_multiplicity)
    return np.amax(np.array(normalization_rate)).astype('f16') / np.array(normalization_rate)


def get_total_n_data_words(files_dict, precise=False):
    n_words = 0
    if precise:  # open all files and determine the total number of words precicely, can take some time
        if len(files_dict) > 10:
            progress_bar = progressbar.ProgressBar(widgets=['', progressbar.Percentage(), ' ', progressbar.Bar(marker='*', left='|', right='|'), ' ', ETA()], maxval=len(files_dict), term_width=80)
            progress_bar.start()
        for index, file_name in enumerate(files_dict.iterkeys()):
            with tb.openFile(file_name, mode="r") as in_file_h5:  # open the actual file
                n_words += in_file_h5.root.raw_data.shape[0]
            if len(files_dict) > 10:
                progress_bar.update(index)
        if len(files_dict) > 10:
            progress_bar.finish()
        return n_words
    else:  # open just first an last file and take the mean to estimate the total numbe rof words
        with tb.openFile(files_dict.keys()[0], mode="r") as in_file_h5:  # open the actual file
            n_words += in_file_h5.root.raw_data.shape[0]
        with tb.openFile(files_dict.keys()[-1], mode="r") as in_file_h5:  # open the actual file
            n_words += in_file_h5.root.raw_data.shape[0]
        return n_words * len(files_dict) / 2


def create_parameter_table(files_dict):
    if not check_parameter_similarity(files_dict):
        raise RuntimeError('Cannot create table from file with different scan parameters.')
    # create the parameter names / format for the parameter table
    try:
        names = ','.join([name for name in files_dict.itervalues().next().keys()])
        formats = ','.join(['int32' for name in files_dict.itervalues().next().keys()])
        arrayList = [l for l in files_dict.itervalues().next().values()]
    except AttributeError:  # no parameters given, return None
        return
    parameter_table = None
    # create a parameter table with an entry for every read out
    for file_name, parameters in files_dict.iteritems():
        with tb.openFile(file_name, mode="r") as in_file_h5:  # open the actual file
            n_parameter_settings = max([len(i) for i in files_dict[file_name].values()])  # determine the number of different parameter settings from the list length of parameter values of the first parameter
            if n_parameter_settings == 0:  # no parameter values, first raw data file has only config info and no other data (meta, raw data, parameter data)
                continue
            try:  # try to combine the scan parameter tables
                if parameter_table is None:  # final parameter_table does not exists, so create is
                    parameter_table = in_file_h5.root.scan_parameters[:]
                else:  # final parameter table already exist, so append to existing
                    parameter_table.resize(parameter_table.shape[0] + in_file_h5.root.scan_parameters[:].shape[0], refcheck=False)  # fastest way to append, http://stackoverflow.com/questions/1730080/append-rows-to-a-numpy-record-array
                    parameter_table[-in_file_h5.root.scan_parameters.shape[0]:] = in_file_h5.root.scan_parameters[:]  # set table
            except tb.NoSuchNodeError:  # there is no scan parameter table, so create one
                read_out = in_file_h5.root.meta_data.shape[0]
                if parameter_table is None:  # final parameter_table does not exists, so create is
                    parameter_table = np.rec.fromarrays(arrayList, names=names, formats=formats)  # create recarray
                    parameter_table.resize(read_out, refcheck=False)
                    parameter_table[-read_out:] = np.rec.fromarrays(arrayList, names=names, formats=formats)
                else:  # final parameter table already exist, so append to existing
                    parameter_table.resize(parameter_table.shape[0] + read_out)  # fastest way to append, http://stackoverflow.com/questions/1730080/append-rows-to-a-numpy-record-array
                    parameter_table[-read_out:] = np.rec.fromarrays([l for l in parameters.values()], names=names, formats=formats)

    return parameter_table


def get_parameter_value_from_file_names(files, parameters=None, unique=False, sort=True):
    """
    Takes a list of files, searches for the parameter name in the file name and returns a ordered dict with the file name
    in the first dimension and the corresponding parameter value in the second.
    The file names can be sorted by the parameter value, otherwise the order is kept. If unique is true every parameter is unique and
    mapped to the file name that occurred last in the files list.

    Parameters
    ----------
    files : list of strings
    parameter : string or list of strings
    unique : bool
    sort : bool

    Returns
    -------
    collections.OrderedDict

    """
#     unique=False
    logging.debug('Get the parameter: ' + str(parameters) + ' values from the file names of ' + str(len(files)) + ' files')
    files_dict = collections.OrderedDict()
    if parameters is None:  # special case, no parameter defined
        return files_dict
    if isinstance(parameters, basestring):
        parameters = (parameters, )
    search_string = '_'.join(parameters)
    for _ in parameters:
        search_string += r'_(-?\d+)'
    result = {}
    for one_file in files:
        parameter_values = re.findall(search_string, one_file)
        if parameter_values:
            if isinstance(parameter_values[0], tuple):
                parameter_values = list(reduce(lambda t1, t2: t1 + t2, parameter_values))
            parameter_values = [[int(i), ] for i in parameter_values]  # convert string value to list with int
            files_dict[one_file] = dict(zip(parameters, parameter_values))
            if unique:  # reduce to the files with different scan parameters
                for key, value in files_dict.items():
                    if value not in result.values():
                        result[key] = value
            else:
                result[one_file] = files_dict[one_file]
    return collections.OrderedDict(sorted(result.iteritems(), key=itemgetter(1)) if sort else files_dict)  # with PEP 265 solution of sorting a dict by value


def get_data_file_names_from_scan_base(scan_base, filter_str=['_analyzed.h5', '_interpreted.h5', '_cut.h5', '_result.h5', '_hists.h5'], sort_by_time=True, meta_data_v2=True):
    """
    Generate a list of .h5 files which have a similar file name.

    Parameters
    ----------
    scan_base : list, string
        List of string or string of the scan base names. The scan_base will be used to search for files containing the string. The .h5 file extension will be added automatically.
    filter : list, string
        List of string or string which are used to filter the returned filenames. File names containing filter_str in the file name will not be returned. Use None to disable filter.
    sort_by_time : bool
        If True, return file name list sorted from oldest to newest. The time from meta table will be used to sort the files.
    meta_data_v2 : bool
        True for new (v2) meta data format, False for the old (v1) format.

    Returns
    -------
    data_files : list
        List of file names matching the obove conditions.
    """
    data_files = []
    if scan_base is None:
        return data_files
    if isinstance(scan_base, basestring):
        scan_base = [scan_base]
    for scan_base_str in scan_base:
        if '.h5' == os.path.splitext(scan_base_str)[1]:
            data_files.append(scan_base_str)
        else:
            data_files.extend(glob.glob(scan_base_str + '*.h5'))

    if filter_str:
        if isinstance(filter_str, basestring):
            filter_str = [filter_str]
        data_files = filter(lambda data_file: not any([(True if x in data_file else False) for x in filter_str]), data_files)
    if sort_by_time and len(data_files) > 1:
        f_list = {}
        for data_file in data_files:
            with tb.open_file(data_file, mode="r") as h5_file:
                try:
                    meta_data = h5_file.root.meta_data
                except IndexError:
                    logging.info("File %s is missing meta_data" % h5_file.filename)
                else:
                    if meta_data_v2:
                        timestamp = meta_data[0]["timestamp_start"]
                    else:
                        timestamp = meta_data[0]["timestamp"]
                    f_list[data_file] = timestamp

        data_files = list(sorted(f_list, key=f_list.__getitem__, reverse=False))
    return data_files


def get_scan_parameter_names(scan_parameters):
    ''' Returns the scan parameter names of the scan_paraemeter table.

    Parameters
    ----------
    scan_parameters : numpy.array
        Can be None

    Returns<commit_msg>REG: fix check for node<commit_after>"""This class provides often needed analysis functions, for analysis that is done with python.
"""

import logging
import re
import os
import time
import collections
import numpy as np
import progressbar
import glob
import tables as tb
import numexpr as ne
from operator import itemgetter
from scipy.interpolate import interp1d
from scipy.interpolate import splrep, splev
from tables import dtype_from_descr
from stat import ST_CTIME

from pybar_fei4_interpreter import analysis_functions

from pybar.daq.fei4_record import FEI4Record
from pybar.analysis.plotting import plotting
from pybar.daq.readout_utils import is_fe_word, is_data_header


class AnalysisError(Exception):

    """Base class for exceptions in this module.
    """


class IncompleteInputError(AnalysisError):

    """Exception raised for errors in the input.
    """


class InvalidInputError(AnalysisError):

    """Exception raised for errors in the input.
    """


class NotSupportedError(AnalysisError):

    """Exception raised for not supported actions.
    """


def generate_threshold_mask(hist):
    '''Masking array elements when equal 0.0 or greater than 10 times the median

    Parameters
    ----------
    hist : array_like
        Input data.

    Returns
    -------
    masked array
        Returns copy of the array with masked elements.
    '''
    masked_array = np.ma.masked_values(hist, 0)
    masked_array = np.ma.masked_greater(masked_array, 10 * np.ma.median(hist))
    logging.info('Masking %d pixel(s)', np.ma.count_masked(masked_array))
    return np.ma.getmaskarray(masked_array)


def unique_row(array, use_columns=None, selected_columns_only=False):
    '''Takes a numpy array and returns the array reduced to unique rows. If columns are defined only these columns are taken to define a unique row.
    The returned array can have all columns of the original array or only the columns defined in use_columns.
    Parameters
    ----------
    array : numpy.ndarray
    use_columns : list
        Index of columns to be used to define a unique row
    selected_columns_only : bool
        If true only the columns defined in use_columns are returned

    Returns
    -------
    numpy.ndarray
    '''
    if array.dtype.names is None:  # normal array has no named dtype
        if use_columns is not None:
            a_cut = array[:, use_columns]
        else:
            a_cut = array
        if len(use_columns) > 1:
            b = np.ascontiguousarray(a_cut).view(np.dtype((np.void, a_cut.dtype.itemsize * a_cut.shape[1])))
        else:
            b = np.ascontiguousarray(a_cut)
        _, index = np.unique(b, return_index=True)
        if not selected_columns_only:
            return array[np.sort(index)]  # sort to preserve order
        else:
            return a_cut[np.sort(index)]  # sort to preserve order
    else:  # names for dtype founnd --> array is recarray
        names = list(array.dtype.names)
        if use_columns is not None:
            new_names = [names[i] for i in use_columns]
        else:
            new_names = names
        a_cut, index = np.unique(array[new_names], return_index=True)
        if not selected_columns_only:
            return array[np.sort(index)]  # sort to preserve order
        else:
            return array[np.sort(index)][new_names]  # sort to preserve order


def get_ranges_from_array(arr, append_last=True):
    '''Takes an array and calculates ranges [start, stop[. The last range end is none to keep the same length.

    Parameters
    ----------
    arr : array like
    append_last: bool
        If True, append item with a pair of last array item and None.

    Returns
    -------
    numpy.array
        The array formed by pairs of values by the given array.

    Example
    -------
    >>> a = np.array((1,2,3,4))
    >>> get_ranges_from_array(a, append_last=True)
    array([[1, 2],
           [2, 3],
           [3, 4],
           [4, None]])
    >>> get_ranges_from_array(a, append_last=False)
    array([[1, 2],
           [2, 3],
           [3, 4]])
    '''
    right = arr[1:]
    if append_last:
        left = arr[:]
        right = np.append(right, None)
    else:
        left = arr[:-1]
    return np.column_stack((left, right))


def get_mean_from_histogram(counts, bin_positions):
    return np.dot(counts, np.array(bin_positions)) / np.sum(counts).astype('f4')


def get_median_from_histogram(counts, bin_positions):
    values = []
    for index, one_bin in enumerate(counts):
        for _ in range(one_bin):
            values.append(bin_positions[index])
    return np.median(values)


def get_rms_from_histogram(counts, bin_positions):
    values = []
    for index, one_bin in enumerate(counts):
        for _ in range(one_bin):
            values.append(bin_positions[index])
    return np.std(values)


def in1d_sorted(ar1, ar2):
    """
    Does the same than np.in1d but uses the fact that ar1 and ar2 are sorted. Is therefore much faster.

    """
    if ar1.shape[0] == 0 or ar2.shape[0] == 0:  # check for empty arrays to avoid crash
        return []
    inds = ar2.searchsorted(ar1)
    inds[inds == len(ar2)] = 0
    return ar2[inds] == ar1


def central_difference(x, y):
    '''Returns the dy/dx(x) via central difference method

    Parameters
    ----------
    x : array like
    y : array like

    Returns
    -------
    dy/dx : array like
    '''
    if (len(x) != len(y)):
        raise ValueError("x, y must have the same length")
    z1 = np.hstack((y[0], y[:-1]))
    z2 = np.hstack((y[1:], y[-1]))
    dx1 = np.hstack((0, np.diff(x)))
    dx2 = np.hstack((np.diff(x), 0))
    return (z2 - z1) / (dx2 + dx1)


def get_profile_histogram(x, y, n_bins=100):
    '''Takes 2D point data (x,y) and creates a profile histogram similar to the TProfile in ROOT. It calculates
    the y mean for every bin at the bin center and gives the y mean error as error bars.

    Parameters
    ----------
    x : array like
        data x positions
    y : array like
        data y positions
    n_bins : int
        the number of bins used to create the histogram
    '''
    if len(x) != len(y):
        raise ValueError('x and y dimensions have to be the same')
    y = y.astype(np.float32)
    n, bin_edges = np.histogram(x, bins=n_bins)  # needed to calculate the number of points per bin
    sy = np.histogram(x, bins=n_bins, weights=y)[0]  # the sum of the bin values
    sy2 = np.histogram(x, bins=n_bins, weights=y * y)[0]  # the quadratic sum of the bin values
    bin_centers = (bin_edges[1:] + bin_edges[:-1]) / 2  # calculate the bin center for all bins
    mean = sy / n  # calculate the mean of all bins
    std = np.sqrt((sy2 / n - mean * mean))  # TODO: not understood, need check if this is really the standard deviation
    std_mean = std / np.sqrt((n - 1))
    mean[np.isnan(mean)] = 0.
    std_mean[np.isnan(std_mean)] = 0.
    return bin_centers, mean, std_mean


def get_rate_normalization(hit_file, parameter, reference='event', cluster_file=None, plot=False, chunk_size=500000):
    ''' Takes different hit files (hit_files), extracts the number of events or the scan time (reference) per scan parameter (parameter)
    and returns an array with a normalization factor. This normalization factor has the length of the number of different parameters.
    If a cluster_file is specified also the number of cluster per event are used to create the normalization factor.

    Parameters
    ----------
    hit_files : string
    parameter : string
    reference : string
    plot : bool

    Returns
    -------
    numpy.ndarray
    '''

    logging.info('Calculate the rate normalization')
    with tb.openFile(hit_file, mode="r+") as in_hit_file_h5:  # open the hit file
        meta_data = in_hit_file_h5.root.meta_data[:]
        scan_parameter = get_scan_parameter(meta_data)[parameter]
        event_numbers = get_meta_data_at_scan_parameter(meta_data, parameter)['event_number']  # get the event numbers in meta_data where the scan parameter changes
        event_range = get_ranges_from_array(event_numbers)
        normalization_rate = []
        normalization_multiplicity = []
        try:
            event_range[-1, 1] = in_hit_file_h5.root.Hits[-1]['event_number'] + 1
        except tb.NoSuchNodeError:
            logging.error('Cannot find hits table')
            return

        # calculate rate normalization from the event rate for triggered data / measurement time for self triggered data for each scan parameter
        if reference == 'event':
            n_events = event_range[:, 1] - event_range[:, 0]  # number of events for every parameter setting
            normalization_rate.extend(n_events)
        elif reference == 'time':
            time_start = get_meta_data_at_scan_parameter(meta_data, parameter)['timestamp_start']
            time_spend = np.diff(time_start)
            time_spend = np.append(time_spend, meta_data[-1]['timestamp_stop'] - time_start[-1])  # TODO: needs check, add last missing entry
            normalization_rate.extend(time_spend)
        else:
            raise NotImplementedError('The normalization reference ' + reference + ' is not implemented')

        if cluster_file:  # calculate the rate normalization from the mean number of hits per event per scan parameter, needed for beam data since a beam since the multiplicity is rarely constant
            cluster_table = in_hit_file_h5.root.Cluster
            index_event_number(cluster_table)
            index = 0  # index where to start the read out, 0 at the beginning, increased during looping, variable for read speed up
            best_chunk_size = chunk_size  # variable for read speed up
            total_cluster = 0
            progress_bar = progressbar.ProgressBar(widgets=['', progressbar.Percentage(), ' ', progressbar.Bar(marker='*', left='|', right='|'), ' ', ETA(smoothing=0.8)], maxval=cluster_table.shape[0], term_width=80)
            progress_bar.start()
            for start_event, stop_event in event_range:  # loop over the selected events
                readout_cluster_len = 0  # variable to calculate a optimal chunk size value from the number of hits for speed up
                n_cluster_per_event = None
                for clusters, index in data_aligned_at_events(cluster_table, start_event_number=start_event, stop_event_number=stop_event, start=index, chunk_size=best_chunk_size):
                    if n_cluster_per_event is None:
                        n_cluster_per_event = get_n_cluster_in_events(clusters['event_number'])[:, 1]  # array with the number of cluster per event, cluster per event are at least 1
                    else:
                        n_cluster_per_event = np.append(n_cluster_per_event, get_n_cluster_in_events(clusters['event_number'])[:, 1])
                    readout_cluster_len += clusters.shape[0]
                    total_cluster += clusters.shape[0]
                    progress_bar.update(index)
                best_chunk_size = int(1.5 * readout_cluster_len) if int(1.05 * readout_cluster_len) < chunk_size else chunk_size  # to increase the readout speed, estimated the number of hits for one read instruction
                normalization_multiplicity.append(np.mean(n_cluster_per_event))
            progress_bar.finish()
            if total_cluster != cluster_table.shape[0]:
                logging.warning('Analysis shows inconsistent number of cluster (%d != %d). Check needed!', total_cluster, cluster_table.shape[0])

    if plot:
        x = scan_parameter
        if reference == 'event':
            plotting.plot_scatter(x, normalization_rate, title='Events per ' + parameter + ' setting', x_label=parameter, y_label='# events', log_x=True, filename=hit_file[:-3] + '_n_event_normalization.pdf')
        elif reference == 'time':
            plotting.plot_scatter(x, normalization_rate, title='Measuring time per GDAC setting', x_label=parameter, y_label='time [s]', log_x=True, filename=hit_file[:-3] + '_time_normalization.pdf')
        if cluster_file:
            plotting.plot_scatter(x, normalization_multiplicity, title='Mean number of particles per event', x_label=parameter, y_label='number of hits per event', log_x=True, filename=hit_file[:-3] + '_n_particles_normalization.pdf')
    print len(normalization_rate), len(normalization_multiplicity)
    if cluster_file:
        normalization_rate = np.array(normalization_rate)
        normalization_multiplicity = np.array(normalization_multiplicity)
        return np.amax(normalization_rate * normalization_multiplicity).astype('f16') / (normalization_rate * normalization_multiplicity)
    return np.amax(np.array(normalization_rate)).astype('f16') / np.array(normalization_rate)


def get_total_n_data_words(files_dict, precise=False):
    n_words = 0
    if precise:  # open all files and determine the total number of words precicely, can take some time
        if len(files_dict) > 10:
            progress_bar = progressbar.ProgressBar(widgets=['', progressbar.Percentage(), ' ', progressbar.Bar(marker='*', left='|', right='|'), ' ', ETA()], maxval=len(files_dict), term_width=80)
            progress_bar.start()
        for index, file_name in enumerate(files_dict.iterkeys()):
            with tb.openFile(file_name, mode="r") as in_file_h5:  # open the actual file
                n_words += in_file_h5.root.raw_data.shape[0]
            if len(files_dict) > 10:
                progress_bar.update(index)
        if len(files_dict) > 10:
            progress_bar.finish()
        return n_words
    else:  # open just first an last file and take the mean to estimate the total numbe rof words
        with tb.openFile(files_dict.keys()[0], mode="r") as in_file_h5:  # open the actual file
            n_words += in_file_h5.root.raw_data.shape[0]
        with tb.openFile(files_dict.keys()[-1], mode="r") as in_file_h5:  # open the actual file
            n_words += in_file_h5.root.raw_data.shape[0]
        return n_words * len(files_dict) / 2


def create_parameter_table(files_dict):
    if not check_parameter_similarity(files_dict):
        raise RuntimeError('Cannot create table from file with different scan parameters.')
    # create the parameter names / format for the parameter table
    try:
        names = ','.join([name for name in files_dict.itervalues().next().keys()])
        formats = ','.join(['int32' for name in files_dict.itervalues().next().keys()])
        arrayList = [l for l in files_dict.itervalues().next().values()]
    except AttributeError:  # no parameters given, return None
        return
    parameter_table = None
    # create a parameter table with an entry for every read out
    for file_name, parameters in files_dict.iteritems():
        with tb.openFile(file_name, mode="r") as in_file_h5:  # open the actual file
            n_parameter_settings = max([len(i) for i in files_dict[file_name].values()])  # determine the number of different parameter settings from the list length of parameter values of the first parameter
            if n_parameter_settings == 0:  # no parameter values, first raw data file has only config info and no other data (meta, raw data, parameter data)
                continue
            try:  # try to combine the scan parameter tables
                if parameter_table is None:  # final parameter_table does not exists, so create is
                    parameter_table = in_file_h5.root.scan_parameters[:]
                else:  # final parameter table already exist, so append to existing
                    parameter_table.resize(parameter_table.shape[0] + in_file_h5.root.scan_parameters[:].shape[0], refcheck=False)  # fastest way to append, http://stackoverflow.com/questions/1730080/append-rows-to-a-numpy-record-array
                    parameter_table[-in_file_h5.root.scan_parameters.shape[0]:] = in_file_h5.root.scan_parameters[:]  # set table
            except tb.NoSuchNodeError:  # there is no scan parameter table, so create one
                read_out = in_file_h5.root.meta_data.shape[0]
                if parameter_table is None:  # final parameter_table does not exists, so create is
                    parameter_table = np.rec.fromarrays(arrayList, names=names, formats=formats)  # create recarray
                    parameter_table.resize(read_out, refcheck=False)
                    parameter_table[-read_out:] = np.rec.fromarrays(arrayList, names=names, formats=formats)
                else:  # final parameter table already exist, so append to existing
                    parameter_table.resize(parameter_table.shape[0] + read_out)  # fastest way to append, http://stackoverflow.com/questions/1730080/append-rows-to-a-numpy-record-array
                    parameter_table[-read_out:] = np.rec.fromarrays([l for l in parameters.values()], names=names, formats=formats)

    return parameter_table


def get_parameter_value_from_file_names(files, parameters=None, unique=False, sort=True):
    """
    Takes a list of files, searches for the parameter name in the file name and returns a ordered dict with the file name
    in the first dimension and the corresponding parameter value in the second.
    The file names can be sorted by the parameter value, otherwise the order is kept. If unique is true every parameter is unique and
    mapped to the file name that occurred last in the files list.

    Parameters
    ----------
    files : list of strings
    parameter : string or list of strings
    unique : bool
    sort : bool

    Returns
    -------
    collections.OrderedDict

    """
#     unique=False
    logging.debug('Get the parameter: ' + str(parameters) + ' values from the file names of ' + str(len(files)) + ' files')
    files_dict = collections.OrderedDict()
    if parameters is None:  # special case, no parameter defined
        return files_dict
    if isinstance(parameters, basestring):
        parameters = (parameters, )
    search_string = '_'.join(parameters)
    for _ in parameters:
        search_string += r'_(-?\d+)'
    result = {}
    for one_file in files:
        parameter_values = re.findall(search_string, one_file)
        if parameter_values:
            if isinstance(parameter_values[0], tuple):
                parameter_values = list(reduce(lambda t1, t2: t1 + t2, parameter_values))
            parameter_values = [[int(i), ] for i in parameter_values]  # convert string value to list with int
            files_dict[one_file] = dict(zip(parameters, parameter_values))
            if unique:  # reduce to the files with different scan parameters
                for key, value in files_dict.items():
                    if value not in result.values():
                        result[key] = value
            else:
                result[one_file] = files_dict[one_file]
    return collections.OrderedDict(sorted(result.iteritems(), key=itemgetter(1)) if sort else files_dict)  # with PEP 265 solution of sorting a dict by value


def get_data_file_names_from_scan_base(scan_base, filter_str=['_analyzed.h5', '_interpreted.h5', '_cut.h5', '_result.h5', '_hists.h5'], sort_by_time=True, meta_data_v2=True):
    """
    Generate a list of .h5 files which have a similar file name.

    Parameters
    ----------
    scan_base : list, string
        List of string or string of the scan base names. The scan_base will be used to search for files containing the string. The .h5 file extension will be added automatically.
    filter : list, string
        List of string or string which are used to filter the returned filenames. File names containing filter_str in the file name will not be returned. Use None to disable filter.
    sort_by_time : bool
        If True, return file name list sorted from oldest to newest. The time from meta table will be used to sort the files.
    meta_data_v2 : bool
        True for new (v2) meta data format, False for the old (v1) format.

    Returns
    -------
    data_files : list
        List of file names matching the obove conditions.
    """
    data_files = []
    if scan_base is None:
        return data_files
    if isinstance(scan_base, basestring):
        scan_base = [scan_base]
    for scan_base_str in scan_base:
        if '.h5' == os.path.splitext(scan_base_str)[1]:
            data_files.append(scan_base_str)
        else:
            data_files.extend(glob.glob(scan_base_str + '*.h5'))

    if filter_str:
        if isinstance(filter_str, basestring):
            filter_str = [filter_str]
        data_files = filter(lambda data_file: not any([(True if x in data_file else False) for x in filter_str]), data_files)
    if sort_by_time and len(data_files) > 1:
        f_list = {}
        for data_file in data_files:
            with tb.open_file(data_file, mode="r") as h5_file:
                try:
                    meta_data = h5_file.root.meta_data
                except NoSuchNodeError:
                    logging.warning("File %s is missing meta_data" % h5_file.filename)
                else:
                    try:
                        if meta_data_v2:
                            timestamp = meta_data[0]["timestamp_start"]
                        else:
                            timestamp = meta_data[0]["timestamp"]
                    except IndexError:
                        logging.info("File %s has empty meta_data" % h5_file.filename)
                    else:
                        f_list[data_file] = timestamp

        data_files = list(sorted(f_list, key=f_list.__getitem__, reverse=False))
    return data_files


def get_scan_parameter_names(scan_parameters):
    ''' Returns the scan parameter names of the scan_paraemeter table.

    Parameters
    ----------
    scan_parameters : numpy.array
        Can be None

    Returns
{'fraction_non_alphanumeric': 0.09372906802171152, 'fraction_numerical': 0.0095392077607114, 'mean_word_length': 2.7240667469464994, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 33, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4801600', 'n_tokens_mistral': 12767, 'n_tokens_neox': 12342, 'n_words': 3911}