def get_transform_fxn(data_beads,
                      mef_values,
                      mef_channels,
                      clustering_fxn=clustering_gmm,
                      clustering_params={},
                      clustering_channels=None,
                      statistic_fxn=FlowCal.stats.median,
                      statistic_params={},
                      selection_fxn=selection_std,
                      selection_params={},
                      fitting_fxn=fit_beads_autofluorescence,
                      fitting_params={},
                      verbose=False,
                      plot=False,
                      plot_dir=None,
                      plot_filename=None,
                      full_output=False):
    """
    Get a transformation function to convert flow cytometry data to MEF.

    Parameters
    ----------
    data_beads : FCSData object
        Flow cytometry data describing calibration beads.
    mef_values : sequence of sequences
        Known MEF values for the calibration bead subpopulations, for each
        channel specified in `mef_channels`. The innermost sequences must
        have the same length (the same number of bead subpopulations must
        exist for each channel).
    mef_channels : int, or str, or list of int, or list of str
        Channels for which to generate transformation functions.
    verbose : bool, optional
        Flag specifying whether to print information about step completion
        and warnings.
    plot : bool, optional
        Flag specifying whether to produce diagnostic plots.
    plot_dir : str, optional
        Directory where to save diagnostics plots. Ignored if `plot` is
        False. If ``plot==True`` and ``plot_dir is None``, plot without
        saving.
    plot_filename : str, optional
        Name to use for plot files. If None, use ``str(data_beads)``.
    full_output : bool, optional
        Flag specifying whether to include intermediate results in the
        output. If `full_output` is True, the function returns a
        `MEFOutput` ``namedtuple`` with fields as described below. If
        `full_output` is False, the function only returns the calculated
        transformation function.

    Returns
    -------
    transform_fxn : function
        Transformation function to convert flow cytometry data from RFI
        units to MEF. This function has the following signature::

            data_mef = transform_fxn(data_rfi, channels)

    mef_channels : int, or str, or list, only if ``full_output==True``
        Channels on which the transformation function has been generated.
        Directly copied from the `mef_channels` argument.

    clustering : dict, only if ``full_output==True``
        Results of the clustering step. The structure of this dictionary
        is::

            clustering = {"labels": np.array}

        A description of each ``"key": value`` is given below.

        "labels" : array
            Array of length ``N``, where ``N`` is the number of events in
            `data_beads`. This array contains labels indicating which
            subpopulation each event has been assigned to by the clustering
            algorithm. Labels range from ``0`` to ``M - 1``, where ``M`` is
            the number of MEF values specified, and therefore the number of
            subpopulations identified by the clustering algorithm.

    statistic : dict, only if ``full_output==True``
        Results of the calculation of bead subpopulations' fluorescence.
        The structure of this dictionary is::

            statistic = {"values": [np.array, ...]}

        A description of each ``"key": value`` is given below.

        "values" : list of arrays
            Each array contains the representative fluorescence values of
            all subpopulations, for a specific fluorescence channel from
            `mef_channels`. Therefore, each array has a length equal to the
            number of subpopulations, and the outer list has as many arrays
            as the number of channels in `mef_channels`.

    selection : dict, only if ``full_output==True``
        Results of the subpopulation selection step. The structure of this
        dictionary is::

            selection = {"rfi": [np.array, ...],
                         "mef": [np.array, ...]}

        A description of each ``"key": value`` is given below.

        "rfi" : list of arrays
            Each array contains the fluorescence values of each selected
            subpopulation in RFI units, for a specific fluorescence channel
            from `mef_channels`. The outer list has as many arrays as the
            number of channels in `mef_channels`. Because the selection
            step may discard subpopulations, each array has a length less
            than or equal to the total number of subpopulations.
            Furthermore, different arrays in this list may not have the
            same length. However, the length of each array is consistent
            with the corresponding array in ``selection["mef"]`` (see
            below).
        "mef" : list of arrays
            Each array contains the fluorescence values of each selected
            subpopulation in MEF units, for a specific fluorescence channel
            from `mef_channels`. The outer list has as many arrays as the
            number of channels in `mef_channels`. Because the selection
            step may discard subpopulations, each array has a length less
            than or equal to the total number of subpopulations.
            Furthermore, different arrays in this list may not have the
            same length. However, the length of each array is consistent
            with the corresponding array in ``selection["rfi"]`` (see
            above).

    fitting : dict, only if ``full_output==True``
        Results of the model fitting step. The structure of this dictionary
        is::

            selection = {"std_crv": [func, ...],
                         "beads_model": [func, ...],
                         "beads_params": [np.array, ...],
                         "beads_model_str": [str, ...],
                         "beads_params_names": [[], ...]}

        A description of each ``"key": value`` is given below.

        "std_crv" : list of functions
            Functions encoding the fitted standard curves, for each channel
            in `mef_channels`. Each element of this list is the ``std_crv``
            output of the fitting function (see required signature of the
            ``fitting_fxn`` optional parameter), after applying it to the
            MEF and RFI fluorescence values of a specific channel from
            `mef_channels` .
        "beads_model" : list of functions
            Functions encoding the fluorescence model of the calibration
            beads, for each channel in `mef_channels`. Each element of this
            list is the ``beads_model`` output of the fitting function (see
            required signature of the ``fitting_fxn`` optional parameter),
            after applying it to the MEF and RFI fluorescence values of a
            specific channel from `mef_channels` .
        "beads_params" : list of arrays
            Fitted parameter values of the bead fluorescence model, for
            each channel in `mef_chanels`. Each element of this list is the
            ``beads_params`` output of the fitting function (see required
            signature of the ``fitting_fxn`` optional parameter), after
            applying it to the MEF and RFI fluorescence values of a
            specific channel from `mef_channels`.
        "beads_model_str" : list of str
            String representation of the bead models used, for each channel
            in `mef_channels`. Each element of this list is the
            ``beads_model_str`` output of the fitting function (see
            required signature of the ``fitting_fxn`` optional parameter),
            after applying it to the MEF and RFI fluorescence values of a
            specific channel from `mef_channels` .
        "beads_params_names" : list of list
            Names of the parameters given in `beads_params`, for each
            channel in `mef_channels`. Each element of this list is the
            ``beads_params_names`` output of the fitting function (see
            required signature of the ``fitting_fxn`` optional parameter),
            after applying it to the MEF and RFI fluorescence values of a
            specific channel from `mef_channels` .

    Other parameters
    ----------------
    clustering_fxn : function, optional
        Function used for clustering, or identification of subpopulations.
        Must have the following signature::

            labels = clustering_fxn(data, n_clusters, **clustering_params)

        where `data` is a NxD FCSData object or numpy array, `n_clusters`
        is the expected number of bead subpopulations, and `labels` is a 1D
        numpy array of length N, assigning each event in `data` to one
        subpopulation.

    clustering_params : dict, optional
        Additional keyword parameters to pass to `clustering_fxn`.

    clustering_channels : list, optional
        Channels used for clustering. If not specified, use `mef_channels`.
        If more than three channels are specified and `plot` is True, only
        a 3D scatter plot will be produced using the first three channels.

    statistic_fxn : function, optional
        Function used to calculate the representative fluorescence of each
        subpopulation. Must have the following signature::

            s = statistic_fxn(data, **statistic_params)

        where `data` is a 1D FCSData object or numpy array, and `s` is a
        float. Statistical functions from numpy, scipy, or FlowCal.stats
        are valid options.

    statistic_params : dict, optional
        Additional keyword parameters to pass to `statistic_fxn`.

    selection_fxn : function, optional
        Function to use for bead population selection. Must have the
        following signature::

            selected_mask = selection_fxn(data_list, **selection_params)

        where `data_list` is a list of FCSData objects, each one containing
        the events of one population, and `selected_mask` is a boolean
        array indicating whether the population has been selected (True) or
        discarded (False). If None, don't use a population selection
        procedure.

    selection_params : dict, optional
        Additional keyword parameters to pass to `selection_fxn`.

    fitting_fxn : function, optional
        Function used to fit the beads fluorescence model and obtain a
        standard curve. Must have the following signature::

            std_crv, beads_model, beads_params, \\
            beads_model_str, beads_params_names = fitting_fxn(
                fl_rfi, fl_mef, **fitting_params)

        where `std_crv` is a function implementing the standard curve,
        `beads_model` is a function implementing the beads fluorescence
        model, `beads_params` is an array containing the fitted parameters
        of the beads model, `beads_model_str` is a string representation
        of the beads model used, `beads_params_names` is a list with the
        parameter names in the same order as they are given in
        `beads_params`, and `fl_rfi` and `fl_mef` are the fluorescence
        values of the beads in RFI units and MEF units, respectively.
        Note that the standard curve and the fitted beads model are not
        necessarily the same.

    fitting_params : dict, optional
        Additional keyword parameters to pass to `fitting_fxn`.

    Notes
    -----
    The steps involved in generating the MEF transformation function are:

    1. The individual subpopulations of beads are first identified using a
       clustering method of choice. Clustering is performed in all
       specified channels simultaneously.
    2. The fluorescence of each subpopulation is calculated, for each
       channel in `mef_channels`.
    3. Some subpopulations are then discarded if they are close to either
       the minimum or the maximum channel range limits. In addition, if the
       MEF value of some subpopulation is unknown (represented as a
       ``np.nan`` in `mef_values`), the whole subpopulation is also
       discarded.
    4. The measured fluorescence of each subpopulation is compared with
       the known MEF values in `mef_values`, and a standard curve function
       is generated using the appropriate MEF model.

    At the end, a transformation function is generated using the calculated
    standard curves, `mef_channels`, and ``FlowCal.transform.to_mef()``.

    Note that applying the resulting transformation function to other
    flow cytometry samples only yields correct results if they have been
    taken at the same settings as the calibration beads, for all channels
    in `mef_channels`.

    Examples
    --------
    Here is a simple application of this function:

    >>> transform_fxn = FlowCal.mef.get_transform_fxn(
    ...     beads_data,
    ...     mef_channels=['FL1', 'FL3'],
    ...     mef_values=[np.array([    0,   646,   1704,   4827,
    ...                           15991, 47609, 135896, 273006],
    ...                 np.array([    0,  1614,   4035,   12025,
    ...                           31896, 95682, 353225, 1077421]],
    ...     )
    >>> sample_mef = transform_fxn(data=sample_rfi,
    ...                            channels=['FL1', 'FL3'])

    Here, we first generate ``transform_fxn`` from flow cytometry data
    contained in ``FCSData`` object ``beads_data``, for channels FL1 and
    FL3, using provided MEF values for each one of these channels. In the
    next line, we use the resulting transformation function to transform
    cell sample data in RFI to MEF.

    More data about intermediate steps can be obtained with the option
    ``full_output=True``:

    >>> get_transform_output = FlowCal.mef.get_transform_fxn(
    ...     beads_data,
    ...     mef_channels=['FL1', 'FL3'],
    ...     mef_values=[np.array([    0,   646,   1704,   4827,
    ...                           15991, 47609, 135896, 273006],
    ...                 np.array([    0,  1614,   4035,   12025,
    ...                           31896, 95682, 353225, 1077421]],
    ...     full_output=True)

    In this case, the output ``get_transform_output`` will be a
    `MEFOutput` ``namedtuple`` similar to the following::

        FlowCal.mef.MEFOutput(
            transform_fxn=<functools.partial object>,
            mef_channels=['FL1', 'FL3'],
            clustering={
                'labels' : [7, 2, 2, ... 4, 3, 5]
            },
            statistic={
                'values' : [np.array([ 101,  150,  231,  433,
                                      1241, 3106, 7774, 9306]),
                            np.array([   3,   30,   71,  204,
                                       704, 2054, 6732, 9912])]
            },
            selection={
                'rfi' : [np.array([  101,    150,    231,    433,
                                    1241,   3106,   7774]),
                         np.array([  30,      71,    204,    704,
                                   2054,    6732])]
                'mef' : [np.array([    0,    646,   1704,   4827,
                                   15991,  47609, 135896]),
                         np.array([ 1614,   4035,  12025,  31896,
                                   95682, 353225])]
            },
            fitting={
                'std_crv' : [<function <lambda>>,
                             <function <lambda>>]
                'beads_model' : [<function <lambda>>,
                                 <function <lambda>>]
                'beads_params' : [np.array([ 1.09e0, 2.02e0, 1.15e3]),
                                  np.array([9.66e-1, 4.17e0, 6.63e1])]
                'beads_model_str' : ['m*log(fl_rfi) + b =\
 log(fl_mef_auto + fl_mef)',
                                     'm*log(fl_rfi) + b =\
 log(fl_mef_auto + fl_mef)']
                'beads_params_names' : [['m', 'b', 'fl_mef_auto],
                                        ['m', 'b', 'fl_mef_auto]]
            },
        )

    """
    if verbose:
        prev_precision = np.get_printoptions()['precision']
        np.set_printoptions(precision=2)

    # Create directory if plot is True
    if plot and plot_dir is not None:
        if not os.path.exists(plot_dir):
            os.makedirs(plot_dir)

    # Default plot filename
    if plot_filename is None:
        plot_filename = str(data_beads)

    # mef_channels and mef_values should be iterables.
    if hasattr(mef_channels, '__iter__') \
            and not isinstance(mef_channels, six.string_types):
        mef_channels = list(mef_channels)
    else:
        mef_channels = [mef_channels]
        mef_values   = [mef_values]

    # Ensure matching number of `mef_values` for all channels (this implies
    # that the calibration beads have the same number of subpopulations for
    # all channels).
    if not np.all([len(mef_values_channel)==len(mef_values[0])
                   for mef_values_channel in mef_values]):
        msg  = "innermost sequences of mef_values must have the same length"
        msg += " (same number of bead subpopulations must exist for each"
        msg += " channel)"
        raise ValueError(msg)

    ###
    # 1. Clustering
    ###

    # If clustering channels not specified, use channels in mef_channels
    if clustering_channels is None:
        clustering_channels = mef_channels

    # Get number of clusters from number of specified MEF values
    n_clusters = len(mef_values[0])

    # Run clustering function
    labels = clustering_fxn(data_beads[:, clustering_channels],
                            n_clusters,
                            **clustering_params)

    # Separate events corresponding to each cluster
    unique_labels = np.array(list(set(labels)))
    populations = [data_beads[labels == i] for i in unique_labels]

    # Sort populations based on distance to the origin
    population_dist = [np.sum((np.mean(population[:,clustering_channels],
                                       axis=0))**2)
                       for population in populations]
    population_sorted_idx = np.argsort(population_dist)
    populations = [populations[i] for i in population_sorted_idx]

    # Print information
    if verbose:
        # Calculate and display percentage of events on each population
        population_count = np.array([population.shape[0]
                                     for population in populations])
        population_perc = population_count * 100.0 / population_count.sum()

        # Print information
        print("Step 1: Clustering")
        print("  Number of populations to find: {}".format(n_clusters))
        print("  Percentage of events in each population:")
        print("    " + str(population_perc))

    # Plot
    if plot:
        if plot_dir is not None:
            savefig = '{}/clustering_{}.png'.format(plot_dir, plot_filename)
        else:
            savefig = None

        # If used one channel for clustering, make histogram
        if len(clustering_channels) == 1:
            plt.figure(figsize=(8,4))
            FlowCal.plot.hist1d(
                populations,
                channel=clustering_channels[0],
                xscale='logicle',
                bins=256,
                alpha=0.75,
                savefig=savefig)

        # If used two channels for clustering, make 2D scatter plot
        elif len(clustering_channels) == 2:
            plt.figure(figsize=(6,4))
            FlowCal.plot.scatter2d(
                populations,
                channels=clustering_channels,
                xscale='logicle',
                yscale='logicle',
                savefig=savefig)

        # If used three channels or more for clustering, make 3D scatter plot
        # with the first three.
        elif len(clustering_channels) >= 3:
            plt.figure(figsize=(8,6))
            FlowCal.plot.scatter3d_and_projections(
                populations,
                channels=clustering_channels[:3],
                xscale='logicle',
                yscale='logicle',
                zscale='logicle',
                savefig=savefig)

        if plot_dir is not None:
            plt.close()

    # Initialize lists to acumulate results
    std_crv_res = []
    if full_output:
        stats_values_res = []
        selected_rfi_res = []
        selected_mef_res = []
        beads_model_res = []
        beads_params_res =[]
        beads_model_str_res =[]
        beads_params_names_res =[]

    # Iterate through each mef channel
    for mef_channel, mef_values_channel in zip(mef_channels, mef_values):

        populations_channel = [population[:, mef_channel]
                               for population in populations]

        ###
        # 2. Calculate statistics in each subpopulation.
        ###

        # Calculate statistics
        stats_values = [statistic_fxn(population, **statistic_params)
                        for population in populations_channel]
        stats_values = np.array(stats_values)

        # Accumulate results
        if full_output:
            stats_values_res.append(stats_values)

        # Print information
        if verbose:
            print("({}) Step 2: Population Statistic".format(mef_channel))
            print("  Fluorescence of each population (RFI):")
            print("    " + str(stats_values))

        ###
        # 3. Select populations to be used for fitting
        ###

        # Select populations based on selection_fxn
        if selection_fxn is not None:
            selected_mask = selection_fxn(
                [population for population in populations_channel],
                **selection_params)
        else:
            selected_mask = np.ones(n_clusters, dtype=bool)

        # Discard values specified as nan in mef_values_channel
        selected_mask = np.logical_and(selected_mask,
                                       ~np.isnan(mef_values_channel))

        # Get selected rfi and mef values
        selected_rfi = stats_values[selected_mask]
        selected_mef = mef_values_channel[selected_mask]

        # Accumulate results
        if full_output:
            selected_rfi_res.append(selected_rfi)
            selected_mef_res.append(selected_mef)

        # Print information
        if verbose:
            print("({}) Step 3: Population Selection".format(mef_channel))
            print("  {} populations selected.".format(len(selected_rfi)))
            print("  Fluorescence of selected populations (RFI):")
            print("    " + str(selected_rfi))
            print("  Fluorescence of selected populations (MEF):")
            print("    " + str(selected_mef))

        # Plot
        if plot:
            # Get colors for each population. Colors are taken from the default
            # colormap in FlowCal.plot, if the population has been selected.
            # Otherwise, the population is displayed in gray.
            color_levels = np.linspace(0, 1, n_clusters)
            colors = [FlowCal.plot.cmap_default(level)
                          if selected else (0.6, 0.6, 0.6)
                      for selected, level in zip(selected_mask, color_levels)]

            # Plot histograms
            plt.figure(figsize=(8,4))
            FlowCal.plot.hist1d(populations,
                                channel=mef_channel,
                                xscale='logicle',
                                bins=256,
                                alpha=0.75,
                                facecolor=colors)

            # Plot a vertical line for each population, with an x coordinate
            # corresponding to their statistic value.
            ylim = plt.ylim()
            for val, color in zip(stats_values, colors):
                plt.plot([val, val], [ylim[0], ylim[1]], color=color)
            plt.ylim(ylim)

            # Save and close
            if plot_dir is not None:
                plt.tight_layout()
                plt.savefig('{}/populations_{}_{}.png'.format(plot_dir,
                                                              mef_channel,
                                                              plot_filename),
                            dpi=FlowCal.plot.savefig_dpi)
                plt.close()

        ###
        # 4. Get standard curve
        ###

        # Fit
        fitting_output = fitting_fxn(selected_rfi,
                                     selected_mef,
                                     **fitting_params)
        std_crv = fitting_output[0]
        beads_model = fitting_output[1]
        beads_params = fitting_output[2]
        beads_model_str = fitting_output[3]
        beads_params_names = fitting_output[4]
        # Accumulate results
        std_crv_res.append(std_crv)
        if full_output:
            beads_model_res.append(beads_model)
            beads_params_res.append(beads_params)
            beads_model_str_res.append(beads_model_str)
            beads_params_names_res.append(beads_params_names)

        # Print information
        if verbose:
            print("({}) Step 4: Standard Curve Fitting".format(mef_channel))
            print("  Parameters of bead fluorescence model:")
            print("    " + str(beads_params))

        # Plot
        if plot:
            # Get channel range
            xlim = populations[0].range(channels=mef_channel)
            # The plot will be made in log scale. If the lower limit of the
            # range is zero or less, replace by one or some lower value, such
            # that the range covers at least five decades.
            if xlim[0] <= 0:
                xlim[0] = min(1., xlim[1]/1e5)

            # Plot standard curve
            plt.figure(figsize=(6,4))
            plot_standard_curve(selected_rfi,
                                selected_mef,
                                beads_model,
                                std_crv,
                                xscale='log',
                                yscale='log',
                                xlim=xlim)
            plt.xlabel('{} (a.u.)'.format(mef_channel))
            plt.ylabel('{} (MEF)'.format(mef_channel))

            # Save if required
            if plot_dir is not None:
                plt.tight_layout()
                plt.savefig('{}/std_crv_{}_{}.png'.format(plot_dir,
                                                          mef_channel,
                                                          plot_filename),
                            dpi=FlowCal.plot.savefig_dpi)
                plt.close()

    # Make output transformation function
    transform_fxn = functools.partial(FlowCal.transform.to_mef,
                                      sc_list=std_crv_res,
                                      sc_channels=mef_channels)

    if verbose:
        np.set_printoptions(precision=prev_precision)

    if full_output:
        # Clustering results
        clustering_res = {}
        clustering_res['labels'] = labels
        # Population stats results
        statistic_res = {}
        statistic_res['values'] = stats_values_res
        # Population selection results
        selection_res = {}
        selection_res['rfi'] = selected_rfi_res
        selection_res['mef'] = selected_mef_res
        # Fitting results
        fitting_res = {}
        fitting_res['std_crv'] = std_crv_res
        fitting_res['beads_model'] = beads_model_res
        fitting_res['beads_params'] = beads_params_res
        fitting_res['beads_model_str'] = beads_model_str_res
        fitting_res['beads_params_names'] = beads_params_names_res

        # Make namedtuple
        fields = ['mef_channels',
                  'transform_fxn',
                  'clustering',
                  'statistic',
                  'selection',
                  'fitting']
        MEFOutput = collections.namedtuple('MEFOutput', fields)
        out = MEFOutput(mef_channels=mef_channels,
                        transform_fxn=transform_fxn,
                        clustering=clustering_res,
                        statistic=statistic_res,
                        selection=selection_res,
                        fitting=fitting_res)
        return out
    else:
        return transform_fxn
Get a transformation function to convert flow cytometry data to MEF.

    Parameters
    ----------
    data_beads : FCSData object
        Flow cytometry data describing calibration beads.
    mef_values : sequence of sequences
        Known MEF values for the calibration bead subpopulations, for each
        channel specified in `mef_channels`. The innermost sequences must
        have the same length (the same number of bead subpopulations must
        exist for each channel).
    mef_channels : int, or str, or list of int, or list of str
        Channels for which to generate transformation functions.
    verbose : bool, optional
        Flag specifying whether to print information about step completion
        and warnings.
    plot : bool, optional
        Flag specifying whether to produce diagnostic plots.
    plot_dir : str, optional
        Directory where to save diagnostics plots. Ignored if `plot` is
        False. If ``plot==True`` and ``plot_dir is None``, plot without
        saving.
    plot_filename : str, optional
        Name to use for plot files. If None, use ``str(data_beads)``.
    full_output : bool, optional
        Flag specifying whether to include intermediate results in the
        output. If `full_output` is True, the function returns a
        `MEFOutput` ``namedtuple`` with fields as described below. If
        `full_output` is False, the function only returns the calculated
        transformation function.

    Returns
    -------
    transform_fxn : function
        Transformation function to convert flow cytometry data from RFI
        units to MEF. This function has the following signature::

            data_mef = transform_fxn(data_rfi, channels)

    mef_channels : int, or str, or list, only if ``full_output==True``
        Channels on which the transformation function has been generated.
        Directly copied from the `mef_channels` argument.

    clustering : dict, only if ``full_output==True``
        Results of the clustering step. The structure of this dictionary
        is::

            clustering = {"labels": np.array}

        A description of each ``"key": value`` is given below.

        "labels" : array
            Array of length ``N``, where ``N`` is the number of events in
            `data_beads`. This array contains labels indicating which
            subpopulation each event has been assigned to by the clustering
            algorithm. Labels range from ``0`` to ``M - 1``, where ``M`` is
            the number of MEF values specified, and therefore the number of
            subpopulations identified by the clustering algorithm.

    statistic : dict, only if ``full_output==True``
        Results of the calculation of bead subpopulations' fluorescence.
        The structure of this dictionary is::

            statistic = {"values": [np.array, ...]}

        A description of each ``"key": value`` is given below.

        "values" : list of arrays
            Each array contains the representative fluorescence values of
            all subpopulations, for a specific fluorescence channel from
            `mef_channels`. Therefore, each array has a length equal to the
            number of subpopulations, and the outer list has as many arrays
            as the number of channels in `mef_channels`.

    selection : dict, only if ``full_output==True``
        Results of the subpopulation selection step. The structure of this
        dictionary is::

            selection = {"rfi": [np.array, ...],
                         "mef": [np.array, ...]}

        A description of each ``"key": value`` is given below.

        "rfi" : list of arrays
            Each array contains the fluorescence values of each selected
            subpopulation in RFI units, for a specific fluorescence channel
            from `mef_channels`. The outer list has as many arrays as the
            number of channels in `mef_channels`. Because the selection
            step may discard subpopulations, each array has a length less
            than or equal to the total number of subpopulations.
            Furthermore, different arrays in this list may not have the
            same length. However, the length of each array is consistent
            with the corresponding array in ``selection["mef"]`` (see
            below).
        "mef" : list of arrays
            Each array contains the fluorescence values of each selected
            subpopulation in MEF units, for a specific fluorescence channel
            from `mef_channels`. The outer list has as many arrays as the
            number of channels in `mef_channels`. Because the selection
            step may discard subpopulations, each array has a length less
            than or equal to the total number of subpopulations.
            Furthermore, different arrays in this list may not have the
            same length. However, the length of each array is consistent
            with the corresponding array in ``selection["rfi"]`` (see
            above).

    fitting : dict, only if ``full_output==True``
        Results of the model fitting step. The structure of this dictionary
        is::

            selection = {"std_crv": [func, ...],
                         "beads_model": [func, ...],
                         "beads_params": [np.array, ...],
                         "beads_model_str": [str, ...],
                         "beads_params_names": [[], ...]}

        A description of each ``"key": value`` is given below.

        "std_crv" : list of functions
            Functions encoding the fitted standard curves, for each channel
            in `mef_channels`. Each element of this list is the ``std_crv``
            output of the fitting function (see required signature of the
            ``fitting_fxn`` optional parameter), after applying it to the
            MEF and RFI fluorescence values of a specific channel from
            `mef_channels` .
        "beads_model" : list of functions
            Functions encoding the fluorescence model of the calibration
            beads, for each channel in `mef_channels`. Each element of this
            list is the ``beads_model`` output of the fitting function (see
            required signature of the ``fitting_fxn`` optional parameter),
            after applying it to the MEF and RFI fluorescence values of a
            specific channel from `mef_channels` .
        "beads_params" : list of arrays
            Fitted parameter values of the bead fluorescence model, for
            each channel in `mef_chanels`. Each element of this list is the
            ``beads_params`` output of the fitting function (see required
            signature of the ``fitting_fxn`` optional parameter), after
            applying it to the MEF and RFI fluorescence values of a
            specific channel from `mef_channels`.
        "beads_model_str" : list of str
            String representation of the bead models used, for each channel
            in `mef_channels`. Each element of this list is the
            ``beads_model_str`` output of the fitting function (see
            required signature of the ``fitting_fxn`` optional parameter),
            after applying it to the MEF and RFI fluorescence values of a
            specific channel from `mef_channels` .
        "beads_params_names" : list of list
            Names of the parameters given in `beads_params`, for each
            channel in `mef_channels`. Each element of this list is the
            ``beads_params_names`` output of the fitting function (see
            required signature of the ``fitting_fxn`` optional parameter),
            after applying it to the MEF and RFI fluorescence values of a
            specific channel from `mef_channels` .

    Other parameters
    ----------------
    clustering_fxn : function, optional
        Function used for clustering, or identification of subpopulations.
        Must have the following signature::

            labels = clustering_fxn(data, n_clusters, **clustering_params)

        where `data` is a NxD FCSData object or numpy array, `n_clusters`
        is the expected number of bead subpopulations, and `labels` is a 1D
        numpy array of length N, assigning each event in `data` to one
        subpopulation.

    clustering_params : dict, optional
        Additional keyword parameters to pass to `clustering_fxn`.

    clustering_channels : list, optional
        Channels used for clustering. If not specified, use `mef_channels`.
        If more than three channels are specified and `plot` is True, only
        a 3D scatter plot will be produced using the first three channels.

    statistic_fxn : function, optional
        Function used to calculate the representative fluorescence of each
        subpopulation. Must have the following signature::

            s = statistic_fxn(data, **statistic_params)

        where `data` is a 1D FCSData object or numpy array, and `s` is a
        float. Statistical functions from numpy, scipy, or FlowCal.stats
        are valid options.

    statistic_params : dict, optional
        Additional keyword parameters to pass to `statistic_fxn`.

    selection_fxn : function, optional
        Function to use for bead population selection. Must have the
        following signature::

            selected_mask = selection_fxn(data_list, **selection_params)

        where `data_list` is a list of FCSData objects, each one containing
        the events of one population, and `selected_mask` is a boolean
        array indicating whether the population has been selected (True) or
        discarded (False). If None, don't use a population selection
        procedure.

    selection_params : dict, optional
        Additional keyword parameters to pass to `selection_fxn`.

    fitting_fxn : function, optional
        Function used to fit the beads fluorescence model and obtain a
        standard curve. Must have the following signature::

            std_crv, beads_model, beads_params, \\
            beads_model_str, beads_params_names = fitting_fxn(
                fl_rfi, fl_mef, **fitting_params)

        where `std_crv` is a function implementing the standard curve,
        `beads_model` is a function implementing the beads fluorescence
        model, `beads_params` is an array containing the fitted parameters
        of the beads model, `beads_model_str` is a string representation
        of the beads model used, `beads_params_names` is a list with the
        parameter names in the same order as they are given in
        `beads_params`, and `fl_rfi` and `fl_mef` are the fluorescence
        values of the beads in RFI units and MEF units, respectively.
        Note that the standard curve and the fitted beads model are not
        necessarily the same.

    fitting_params : dict, optional
        Additional keyword parameters to pass to `fitting_fxn`.

    Notes
    -----
    The steps involved in generating the MEF transformation function are:

    1. The individual subpopulations of beads are first identified using a
       clustering method of choice. Clustering is performed in all
       specified channels simultaneously.
    2. The fluorescence of each subpopulation is calculated, for each
       channel in `mef_channels`.
    3. Some subpopulations are then discarded if they are close to either
       the minimum or the maximum channel range limits. In addition, if the
       MEF value of some subpopulation is unknown (represented as a
       ``np.nan`` in `mef_values`), the whole subpopulation is also
       discarded.
    4. The measured fluorescence of each subpopulation is compared with
       the known MEF values in `mef_values`, and a standard curve function
       is generated using the appropriate MEF model.

    At the end, a transformation function is generated using the calculated
    standard curves, `mef_channels`, and ``FlowCal.transform.to_mef()``.

    Note that applying the resulting transformation function to other
    flow cytometry samples only yields correct results if they have been
    taken at the same settings as the calibration beads, for all channels
    in `mef_channels`.

    Examples
    --------
    Here is a simple application of this function:

    >>> transform_fxn = FlowCal.mef.get_transform_fxn(
    ...     beads_data,
    ...     mef_channels=['FL1', 'FL3'],
    ...     mef_values=[np.array([    0,   646,   1704,   4827,
    ...                           15991, 47609, 135896, 273006],
    ...                 np.array([    0,  1614,   4035,   12025,
    ...                           31896, 95682, 353225, 1077421]],
    ...     )
    >>> sample_mef = transform_fxn(data=sample_rfi,
    ...                            channels=['FL1', 'FL3'])

    Here, we first generate ``transform_fxn`` from flow cytometry data
    contained in ``FCSData`` object ``beads_data``, for channels FL1 and
    FL3, using provided MEF values for each one of these channels. In the
    next line, we use the resulting transformation function to transform
    cell sample data in RFI to MEF.

    More data about intermediate steps can be obtained with the option
    ``full_output=True``:

    >>> get_transform_output = FlowCal.mef.get_transform_fxn(
    ...     beads_data,
    ...     mef_channels=['FL1', 'FL3'],
    ...     mef_values=[np.array([    0,   646,   1704,   4827,
    ...                           15991, 47609, 135896, 273006],
    ...                 np.array([    0,  1614,   4035,   12025,
    ...                           31896, 95682, 353225, 1077421]],
    ...     full_output=True)

    In this case, the output ``get_transform_output`` will be a
    `MEFOutput` ``namedtuple`` similar to the following::

        FlowCal.mef.MEFOutput(
            transform_fxn=<functools.partial object>,
            mef_channels=['FL1', 'FL3'],
            clustering={
                'labels' : [7, 2, 2, ... 4, 3, 5]
            },
            statistic={
                'values' : [np.array([ 101,  150,  231,  433,
                                      1241, 3106, 7774, 9306]),
                            np.array([   3,   30,   71,  204,
                                       704, 2054, 6732, 9912])]
            },
            selection={
                'rfi' : [np.array([  101,    150,    231,    433,
                                    1241,   3106,   7774]),
                         np.array([  30,      71,    204,    704,
                                   2054,    6732])]
                'mef' : [np.array([    0,    646,   1704,   4827,
                                   15991,  47609, 135896]),
                         np.array([ 1614,   4035,  12025,  31896,
                                   95682, 353225])]
            },
            fitting={
                'std_crv' : [<function <lambda>>,
                             <function <lambda>>]
                'beads_model' : [<function <lambda>>,
                                 <function <lambda>>]
                'beads_params' : [np.array([ 1.09e0, 2.02e0, 1.15e3]),
                                  np.array([9.66e-1, 4.17e0, 6.63e1])]
                'beads_model_str' : ['m*log(fl_rfi) + b =\
 log(fl_mef_auto + fl_mef)',
                                     'm*log(fl_rfi) + b =\
 log(fl_mef_auto + fl_mef)']
                'beads_params_names' : [['m', 'b', 'fl_mef_auto],
                                        ['m', 'b', 'fl_mef_auto]]
            },
        )