latools/filtering/signal_optimiser.py
def signal_optimiser(d, analytes, min_points=5,
                     threshold_mode='kde_first_max',
                     threshold_mult=1., x_bias=0,
                     weights=None, ind=None, mode='minimise'):
    """
    Optimise data selection based on specified analytes.

    Identifies the longest possible contiguous data region in
    the signal where the relative standard deviation (std) and 
    concentration of all analytes is minimised.

    Optimisation is performed via a grid search of all possible
    contiguous data regions. For each region, the mean std and
    mean scaled analyte concentration ('amplitude') are calculated. 
    
    The size and position of the optimal data region are identified 
    using threshold std and amplitude values. Thresholds are derived
    from all calculated stds and amplitudes using the method specified
    by `threshold_mode`. For example, using the 'kde_max' method, a
    probability density function (PDF) is calculated for std and
    amplitude values, and the threshold is set as the maximum of the
    PDF. These thresholds are then used to identify the size and position
    of the longest contiguous region where the std is below the threshold, 
    and the amplitude is either below the threshold.

    All possible regions of the data that have at least
    `min_points` are considered.

    For a graphical demonstration of the action of signal_optimiser, 
    use `optimisation_plot`. 

    Parameters
    ----------
    d : latools.D object
        An latools data object.
    analytes : str or array-like
        Which analytes to consider.
    min_points : int
        The minimum number of contiguous points to
        consider.
    threshold_mode : str
        The method used to calculate the optimisation
        thresholds. Can be 'mean', 'median', 'kde_max'
        or 'bayes_mvs', or a custom function. If a
        function, must take a 1D array, and return a
        single, real number.
    threshood_mult : float or tuple
        A multiplier applied to the calculated threshold
        before use. If a tuple, the first value is applied
        to the mean threshold, and the second is applied to
        the standard deviation threshold. Reduce this to make
        data selection more stringent.
    x_bias : float
        If non-zero, a bias is applied to the calculated statistics
        to prefer the beginning (if > 0) or end (if < 0) of the
        signal. Should be between zero and 1.
    weights : array-like of length len(analytes)
        An array of numbers specifying the importance of
        each analyte considered. Larger number makes the
        analyte have a greater effect on the optimisation.
        Default is None.
    ind : boolean array
        A boolean array the same length as the data. Where
        false, data will not be included.
    mode : str
        Whether to 'minimise' or 'maximise' the concentration
        of the elements.

    Returns
    -------
    dict, str : optimisation result, error message
    """
    errmsg = ''

    if isinstance(analytes, str):
        analytes = [analytes]
    
    if ind is None:
        ind = np.full(len(d.Time), True)

    # initial catch
    if not any(ind) or (np.diff(bool_2_indices(ind)).max() < min_points):
        errmsg = 'Optmisation failed. No contiguous data regions longer than {:.0f} points.'.format(min_points)
        return Bunch({'means': np.nan,
                      'stds': np.nan,
                      'mean_threshold': np.nan,
                      'std_threshold': np.nan,
                      'lims': np.nan,
                      'filt': ind,
                      'threshold_mode': threshold_mode,
                      'min_points': min_points,
                      'analytes': analytes,
                      'opt_centre': np.nan,
                      'opt_n_points': np.nan,
                      'weights': weights,
                      'optimisation_success': False,
                      'errmsg': errmsg}), errmsg

    msmeans, msstds = calculate_optimisation_stats(d, analytes, min_points, weights, ind, x_bias)
    
    # second catch
    if all(np.isnan(msmeans).flat) or all(np.isnan(msmeans).flat):
        errmsg = 'Optmisation failed. No contiguous data regions longer than {:.0f} points.'.format(min_points)
        return Bunch({'means': np.nan,
                      'stds': np.nan,
                      'mean_threshold': np.nan,
                      'std_threshold': np.nan,
                      'lims': np.nan,
                      'filt': ind,
                      'threshold_mode': threshold_mode,
                      'min_points': min_points,
                      'analytes': analytes,
                      'opt_centre': np.nan,
                      'opt_n_points': np.nan,
                      'weights': weights,
                      'optimisation_success': False,
                      'errmsg': errmsg}), errmsg

    # define thresholds
    valid = ['kde_first_max', 'kde_max', 'median', 'bayes_mvs', 'mean']
    n_under = 0
    i = np.argwhere(np.array(valid) == threshold_mode)[0, 0]
    o_threshold_mode = threshold_mode
    while (n_under <= 0) & (i < len(valid)):
        if threshold_mode == 'median':
            # median - OK, but best?
            std_threshold = np.nanmedian(msstds)
            mean_threshold = np.nanmedian(msmeans)
        elif threshold_mode == 'mean':
            # mean
            std_threshold = np.nanmean(msstds)
            mean_threshold = np.nanmean(msmeans)
        elif threshold_mode == 'kde_max':
            # maximum of gaussian kernel density estimator
            mkd = gaussian_kde(msmeans[~np.isnan(msmeans)].flat)
            xm = np.linspace(*np.percentile(msmeans.flatten()[~np.isnan(msmeans.flatten())], (1, 99)), 100)
            mdf = mkd.pdf(xm)
            mean_threshold = xm[np.argmax(mdf)]

            rkd = gaussian_kde(msstds[~np.isnan(msstds)])
            xr = np.linspace(*np.percentile(msstds.flatten()[~np.isnan(msstds.flatten())], (1, 99)), 100)
            rdf = rkd.pdf(xr)
            std_threshold = xr[np.argmax(rdf)]
        elif threshold_mode == 'kde_first_max':
            # first local maximum of gaussian kernel density estimator
            mkd = gaussian_kde(msmeans[~np.isnan(msmeans)].flat)
            xm = np.linspace(*np.percentile(msmeans.flatten()[~np.isnan(msmeans.flatten())], (1, 99)), 100)
            mdf = mkd.pdf(xm)
            inds = np.argwhere(np.r_[False, mdf[1:] > mdf[:-1]] & 
                            np.r_[mdf[:-1] > mdf[1:], False] & 
                            (mdf > 0.25 * mdf.max()))
            mean_threshold = xm[np.min(inds)]

            rkd = gaussian_kde(msstds[~np.isnan(msstds)])
            xr = np.linspace(*np.percentile(msstds.flatten()[~np.isnan(msstds.flatten())], (1, 99)), 100)
            rdf = rkd.pdf(xr)
            inds = np.argwhere(np.r_[False, rdf[1:] > rdf[:-1]] & 
                            np.r_[rdf[:-1] > rdf[1:], False] & 
                            (rdf > 0.25 * rdf.max()))
            std_threshold = xr[np.min(inds)]
        elif threshold_mode == 'bayes_mvs':
            # bayesian mvs.
            bm, _, bs = bayes_mvs(msstds[~np.isnan(msstds)])
            std_threshold = bm.statistic

            bm, _, bs = bayes_mvs(msmeans[~np.isnan(msmeans)])
            mean_threshold = bm.statistic
        elif callable(threshold_mode):
            std_threshold = threshold_mode(msstds[~np.isnan(msstds)].flatten())
            mean_threshold = threshold_mode(msmeans[~np.isnan(msmeans)].flatten())
        else:
            try:
                mean_threshold, std_threshold = threshold_mode
            except:
                raise ValueError('\nthreshold_mode must be one of:\n   ' + ', '.join(valid) + ',\na custom function, or a \n(mean_threshold, std_threshold) tuple.')

        # apply threshold_mult
        if isinstance(threshold_mult, (int, float)):
            std_threshold *= threshold_mult
            mean_threshold *= threshold_mult
        elif len(threshold_mult) == 2:
            mean_threshold *= threshold_mult[0]
            std_threshold *= threshold_mult[1]
        else:
            raise ValueError('\nthreshold_mult must be a float, int or tuple of length 2.')

        rind = (msstds < std_threshold)

        if mode == 'minimise':
            mind = (msmeans < mean_threshold)
        else:
            mind = (msmeans > mean_threshold)

        ind = rind & mind

        n_under = ind.sum()
        if n_under == 0:
            i += 1
            if i <= len(valid) - 1:
                threshold_mode = valid[i]
            else:
                errmsg = 'Optimisation failed. No of the threshold_mode would work. Try reducting min_points.'
                return Bunch({'means': np.nan,
                              'stds': np.nan,
                              'mean_threshold': np.nan,
                              'std_threshold': np.nan,
                              'lims': np.nan,
                              'filt': ind,
                              'threshold_mode': threshold_mode,
                              'min_points': min_points,
                              'analytes': analytes,
                              'opt_centre': np.nan,
                              'opt_n_points': np.nan,
                              'weights': weights,
                              'optimisation_success': False,
                              'errmsg': errmsg}), errmsg
    
    if i > 0:
        errmsg = "optimisation failed using threshold_mode='{:}', falling back to '{:}'".format(o_threshold_mode, threshold_mode)

    # identify max number of points within thresholds
    passing = np.argwhere(ind)
    opt_n_points = passing[:, 0].max()
    opt_centre = passing[passing[:, 0] == opt_n_points, 1].min()

    opt_n_points += min_points

    # centres, npoints = np.meshgrid(np.arange(msmeans.shape[1]),
    #                                np.arange(min_points, min_points + msmeans.shape[0]))
    # opt_n_points = npoints[ind].max()
    # plus/minus one point to allow some freedom to shift selection window.
    # cind = ind & (npoints == opt_n_points)
    # opt_centre = centres[cind].min()

    if opt_n_points % 2 == 0:
        lims = (opt_centre - opt_n_points // 2,
                opt_centre + opt_n_points // 2)
    else:
        lims = (opt_centre - opt_n_points // 2,
                opt_centre + opt_n_points // 2 + 1)

    filt = np.zeros(d.Time.shape, dtype=bool)
    filt[lims[0]:lims[1]] = True

    return Bunch({'means': msmeans,
                  'stds': msstds,
                  'mean_threshold': mean_threshold,
                  'std_threshold': std_threshold,
                  'lims': lims,
                  'filt': filt,
                  'threshold_mode': threshold_mode,
                  'min_points': min_points,
                  'analytes': analytes,
                  'opt_centre': opt_centre,
                  'opt_n_points': opt_n_points,
                  'weights': weights,
                  'optimisation_success': True,
                  'errmsg': errmsg}), errmsg
def signal_optimiser(d, analytes, min_points=5,
                     threshold_mode='kde_first_max',
                     threshold_mult=1., x_bias=0,
                     weights=None, ind=None, mode='minimise'):
    """
    Optimise data selection based on specified analytes.

    Identifies the longest possible contiguous data region in
    the signal where the relative standard deviation (std) and 
    concentration of all analytes is minimised.

    Optimisation is performed via a grid search of all possible
    contiguous data regions. For each region, the mean std and
    mean scaled analyte concentration ('amplitude') are calculated. 
    
    The size and position of the optimal data region are identified 
    using threshold std and amplitude values. Thresholds are derived
    from all calculated stds and amplitudes using the method specified
    by `threshold_mode`. For example, using the 'kde_max' method, a
    probability density function (PDF) is calculated for std and
    amplitude values, and the threshold is set as the maximum of the
    PDF. These thresholds are then used to identify the size and position
    of the longest contiguous region where the std is below the threshold, 
    and the amplitude is either below the threshold.

    All possible regions of the data that have at least
    `min_points` are considered.

    For a graphical demonstration of the action of signal_optimiser, 
    use `optimisation_plot`. 

    Parameters
    ----------
    d : latools.D object
        An latools data object.
    analytes : str or array-like
        Which analytes to consider.
    min_points : int
        The minimum number of contiguous points to
        consider.
    threshold_mode : str
        The method used to calculate the optimisation
        thresholds. Can be 'mean', 'median', 'kde_max'
        or 'bayes_mvs', or a custom function. If a
        function, must take a 1D array, and return a
        single, real number.
    threshood_mult : float or tuple
        A multiplier applied to the calculated threshold
        before use. If a tuple, the first value is applied
        to the mean threshold, and the second is applied to
        the standard deviation threshold. Reduce this to make
        data selection more stringent.
    x_bias : float
        If non-zero, a bias is applied to the calculated statistics
        to prefer the beginning (if > 0) or end (if < 0) of the
        signal. Should be between zero and 1.
    weights : array-like of length len(analytes)
        An array of numbers specifying the importance of
        each analyte considered. Larger number makes the
        analyte have a greater effect on the optimisation.
        Default is None.
    ind : boolean array
        A boolean array the same length as the data. Where
        false, data will not be included.
    mode : str
        Whether to 'minimise' or 'maximise' the concentration
        of the elements.

    Returns
    -------
    dict, str : optimisation result, error message
    """
    errmsg = ''

    if isinstance(analytes, str):
        analytes = [analytes]
    
    if ind is None:
        ind = np.full(len(d.Time), True)

    # initial catch
    if not any(ind) or (np.diff(bool_2_indices(ind)).max() < min_points):
        errmsg = 'Optmisation failed. No contiguous data regions longer than {:.0f} points.'.format(min_points)
        return Bunch({'means': np.nan,
                      'stds': np.nan,
                      'mean_threshold': np.nan,
                      'std_threshold': np.nan,
                      'lims': np.nan,
                      'filt': ind,
                      'threshold_mode': threshold_mode,
                      'min_points': min_points,
                      'analytes': analytes,
                      'opt_centre': np.nan,
                      'opt_n_points': np.nan,
                      'weights': weights,
                      'optimisation_success': False,
                      'errmsg': errmsg}), errmsg

    msmeans, msstds = calculate_optimisation_stats(d, analytes, min_points, weights, ind, x_bias)
    
    # second catch
    if all(np.isnan(msmeans).flat) or all(np.isnan(msmeans).flat):
        errmsg = 'Optmisation failed. No contiguous data regions longer than {:.0f} points.'.format(min_points)
        return Bunch({'means': np.nan,
                      'stds': np.nan,
                      'mean_threshold': np.nan,
                      'std_threshold': np.nan,
                      'lims': np.nan,
                      'filt': ind,
                      'threshold_mode': threshold_mode,
                      'min_points': min_points,
                      'analytes': analytes,
                      'opt_centre': np.nan,
                      'opt_n_points': np.nan,
                      'weights': weights,
                      'optimisation_success': False,
                      'errmsg': errmsg}), errmsg

    # define thresholds
    valid = ['kde_first_max', 'kde_max', 'median', 'bayes_mvs', 'mean']
    n_under = 0
    i = np.argwhere(np.array(valid) == threshold_mode)[0, 0]
    o_threshold_mode = threshold_mode
    while (n_under <= 0) & (i < len(valid)):
        if threshold_mode == 'median':
            # median - OK, but best?
            std_threshold = np.nanmedian(msstds)
            mean_threshold = np.nanmedian(msmeans)
        elif threshold_mode == 'mean':
            # mean
            std_threshold = np.nanmean(msstds)
            mean_threshold = np.nanmean(msmeans)
        elif threshold_mode == 'kde_max':
            # maximum of gaussian kernel density estimator
            mkd = gaussian_kde(msmeans[~np.isnan(msmeans)].flat)
            xm = np.linspace(*np.percentile(msmeans.flatten()[~np.isnan(msmeans.flatten())], (1, 99)), 100)
            mdf = mkd.pdf(xm)
            mean_threshold = xm[np.argmax(mdf)]

            rkd = gaussian_kde(msstds[~np.isnan(msstds)])
            xr = np.linspace(*np.percentile(msstds.flatten()[~np.isnan(msstds.flatten())], (1, 99)), 100)
            rdf = rkd.pdf(xr)
            std_threshold = xr[np.argmax(rdf)]
        elif threshold_mode == 'kde_first_max':
            # first local maximum of gaussian kernel density estimator
            mkd = gaussian_kde(msmeans[~np.isnan(msmeans)].flat)
            xm = np.linspace(*np.percentile(msmeans.flatten()[~np.isnan(msmeans.flatten())], (1, 99)), 100)
            mdf = mkd.pdf(xm)
            inds = np.argwhere(np.r_[False, mdf[1:] > mdf[:-1]] & 
                            np.r_[mdf[:-1] > mdf[1:], False] & 
                            (mdf > 0.25 * mdf.max()))
            mean_threshold = xm[np.min(inds)]

            rkd = gaussian_kde(msstds[~np.isnan(msstds)])
            xr = np.linspace(*np.percentile(msstds.flatten()[~np.isnan(msstds.flatten())], (1, 99)), 100)
            rdf = rkd.pdf(xr)
            inds = np.argwhere(np.r_[False, rdf[1:] > rdf[:-1]] & 
                            np.r_[rdf[:-1] > rdf[1:], False] & 
                            (rdf > 0.25 * rdf.max()))
            std_threshold = xr[np.min(inds)]
        elif threshold_mode == 'bayes_mvs':
            # bayesian mvs.
            bm, _, bs = bayes_mvs(msstds[~np.isnan(msstds)])
            std_threshold = bm.statistic

            bm, _, bs = bayes_mvs(msmeans[~np.isnan(msmeans)])
            mean_threshold = bm.statistic
        elif callable(threshold_mode):
            std_threshold = threshold_mode(msstds[~np.isnan(msstds)].flatten())
            mean_threshold = threshold_mode(msmeans[~np.isnan(msmeans)].flatten())
        else:
            try:
                mean_threshold, std_threshold = threshold_mode
            except:
                raise ValueError('\nthreshold_mode must be one of:\n   ' + ', '.join(valid) + ',\na custom function, or a \n(mean_threshold, std_threshold) tuple.')

        # apply threshold_mult
        if isinstance(threshold_mult, (int, float)):
            std_threshold *= threshold_mult
            mean_threshold *= threshold_mult
        elif len(threshold_mult) == 2:
            mean_threshold *= threshold_mult[0]
            std_threshold *= threshold_mult[1]
        else:
            raise ValueError('\nthreshold_mult must be a float, int or tuple of length 2.')

        rind = (msstds < std_threshold)

        if mode == 'minimise':
            mind = (msmeans < mean_threshold)
        else:
            mind = (msmeans > mean_threshold)

        ind = rind & mind

        n_under = ind.sum()
        if n_under == 0:
            i += 1
            if i <= len(valid) - 1:
                threshold_mode = valid[i]
            else:
                errmsg = 'Optimisation failed. No of the threshold_mode would work. Try reducting min_points.'
                return Bunch({'means': np.nan,
                              'stds': np.nan,
                              'mean_threshold': np.nan,
                              'std_threshold': np.nan,
                              'lims': np.nan,
                              'filt': ind,
                              'threshold_mode': threshold_mode,
                              'min_points': min_points,
                              'analytes': analytes,
                              'opt_centre': np.nan,
                              'opt_n_points': np.nan,
                              'weights': weights,
                              'optimisation_success': False,
                              'errmsg': errmsg}), errmsg
    
    if i > 0:
        errmsg = "optimisation failed using threshold_mode='{:}', falling back to '{:}'".format(o_threshold_mode, threshold_mode)

    # identify max number of points within thresholds
    passing = np.argwhere(ind)
    opt_n_points = passing[:, 0].max()
    opt_centre = passing[passing[:, 0] == opt_n_points, 1].min()

    opt_n_points += min_points

    # centres, npoints = np.meshgrid(np.arange(msmeans.shape[1]),
    #                                np.arange(min_points, min_points + msmeans.shape[0]))
    # opt_n_points = npoints[ind].max()
    # plus/minus one point to allow some freedom to shift selection window.
    # cind = ind & (npoints == opt_n_points)
    # opt_centre = centres[cind].min()

    if opt_n_points % 2 == 0:
        lims = (opt_centre - opt_n_points // 2,
                opt_centre + opt_n_points // 2)
    else:
        lims = (opt_centre - opt_n_points // 2,
                opt_centre + opt_n_points // 2 + 1)

    filt = np.zeros(d.Time.shape, dtype=bool)
    filt[lims[0]:lims[1]] = True

    return Bunch({'means': msmeans,
                  'stds': msstds,
                  'mean_threshold': mean_threshold,
                  'std_threshold': std_threshold,
                  'lims': lims,
                  'filt': filt,
                  'threshold_mode': threshold_mode,
                  'min_points': min_points,
                  'analytes': analytes,
                  'opt_centre': opt_centre,
                  'opt_n_points': opt_n_points,
                  'weights': weights,
                  'optimisation_success': True,
                  'errmsg': errmsg}), errmsg
Optimise data selection based on specified analytes.

    Identifies the longest possible contiguous data region in
    the signal where the relative standard deviation (std) and 
    concentration of all analytes is minimised.

    Optimisation is performed via a grid search of all possible
    contiguous data regions. For each region, the mean std and
    mean scaled analyte concentration ('amplitude') are calculated. 
    
    The size and position of the optimal data region are identified 
    using threshold std and amplitude values. Thresholds are derived
    from all calculated stds and amplitudes using the method specified
    by `threshold_mode`. For example, using the 'kde_max' method, a
    probability density function (PDF) is calculated for std and
    amplitude values, and the threshold is set as the maximum of the
    PDF. These thresholds are then used to identify the size and position
    of the longest contiguous region where the std is below the threshold, 
    and the amplitude is either below the threshold.

    All possible regions of the data that have at least
    `min_points` are considered.

    For a graphical demonstration of the action of signal_optimiser, 
    use `optimisation_plot`. 

    Parameters
    ----------
    d : latools.D object
        An latools data object.
    analytes : str or array-like
        Which analytes to consider.
    min_points : int
        The minimum number of contiguous points to
        consider.
    threshold_mode : str
        The method used to calculate the optimisation
        thresholds. Can be 'mean', 'median', 'kde_max'
        or 'bayes_mvs', or a custom function. If a
        function, must take a 1D array, and return a
        single, real number.
    threshood_mult : float or tuple
        A multiplier applied to the calculated threshold
        before use. If a tuple, the first value is applied
        to the mean threshold, and the second is applied to
        the standard deviation threshold. Reduce this to make
        data selection more stringent.
    x_bias : float
        If non-zero, a bias is applied to the calculated statistics
        to prefer the beginning (if > 0) or end (if < 0) of the
        signal. Should be between zero and 1.
    weights : array-like of length len(analytes)
        An array of numbers specifying the importance of
        each analyte considered. Larger number makes the
        analyte have a greater effect on the optimisation.
        Default is None.
    ind : boolean array
        A boolean array the same length as the data. Where
        false, data will not be included.
    mode : str
        Whether to 'minimise' or 'maximise' the concentration
        of the elements.

    Returns
    -------
    dict, str : optimisation result, error message