Finds the outliers for a given population.

    :param mds: the ``mds`` information about each samples.
    :param centers: the centers of the three reference population clusters.
    :param center_info: the label of the three reference population clusters.
    :param ref_pop: the reference population for which we need the outliers
                    from.
    :param options: the options

    :type mds: numpy.recarray
    :type centers: numpy.array
    :type center_info: dict
    :type ref_pop: str
    :type options: argparse.Namespace

    :returns: a :py:class:`set` of outliers from the ``ref_pop`` population.

    Perform a ``KMeans`` classification using the three centers from the three
    reference population cluster.

    Samples are outliers of the required reference population (``ref_pop``) if:

    * the sample is part of another reference population cluster;
    * the sample is an outlier of the desired reference population
      (``ref_pop``).

    A sample is an outlier of a given cluster :math:`C_j` if the distance
    between this sample and the center of the cluster :math:`C_j` (:math:`O_j`)
    is bigger than a constant times the cluster's standard deviation
    :math:`\\sigma_j`.

    .. math::
        \\sigma_j = \\sqrt{\\frac{\\sum{d(s_i,O_j)^2}}{||C_j|| - 1}}

    where :math:`||C_j||` is the number of samples in the cluster :math:`C_j`,
    and :math:`d(s_i,O_j)` is the distance between the sample :math:`s_i` and
    the center :math:`O_j` of the cluster :math:`C_j`.

    .. math::
        d(s_i, O_j) = \\sqrt{(x_{O_j} - x_{s_i})^2 + (y_{O_j} - y_{s_i})^2}

    Using a constant equals to one ensure we remove 100% of the outliers from
    the cluster. Using a constant of 1.6 or 1.9 ensures we remove 99% and 95%
    of outliers, respectively (an error rate of 1% and 5%, respectively).

def find_outliers(mds, centers, center_info, ref_pop, options):
    """Finds the outliers for a given population.

    :param mds: the ``mds`` information about each samples.
    :param centers: the centers of the three reference population clusters.
    :param center_info: the label of the three reference population clusters.
    :param ref_pop: the reference population for which we need the outliers
                    from.
    :param options: the options

    :type mds: numpy.recarray
    :type centers: numpy.array
    :type center_info: dict
    :type ref_pop: str
    :type options: argparse.Namespace

    :returns: a :py:class:`set` of outliers from the ``ref_pop`` population.

    Perform a ``KMeans`` classification using the three centers from the three
    reference population cluster.

    Samples are outliers of the required reference population (``ref_pop``) if:

    * the sample is part of another reference population cluster;
    * the sample is an outlier of the desired reference population
      (``ref_pop``).

    A sample is an outlier of a given cluster :math:`C_j` if the distance
    between this sample and the center of the cluster :math:`C_j` (:math:`O_j`)
    is bigger than a constant times the cluster's standard deviation
    :math:`\\sigma_j`.

    .. math::
        \\sigma_j = \\sqrt{\\frac{\\sum{d(s_i,O_j)^2}}{||C_j|| - 1}}

    where :math:`||C_j||` is the number of samples in the cluster :math:`C_j`,
    and :math:`d(s_i,O_j)` is the distance between the sample :math:`s_i` and
    the center :math:`O_j` of the cluster :math:`C_j`.

    .. math::
        d(s_i, O_j) = \\sqrt{(x_{O_j} - x_{s_i})^2 + (y_{O_j} - y_{s_i})^2}

    Using a constant equals to one ensure we remove 100% of the outliers from
    the cluster. Using a constant of 1.6 or 1.9 ensures we remove 99% and 95%
    of outliers, respectively (an error rate of 1% and 5%, respectively).

    """
    # Importing matplotlib for plotting purposes
    import matplotlib as mpl
    if options.format != "X11" and mpl.get_backend() != "agg":
        mpl.use("Agg")
    import matplotlib.pyplot as plt
    if options.format != "X11":
        plt.ioff()
    import matplotlib.mlab as mlab

    from sklearn.cluster import KMeans
    from sklearn.metrics.pairwise import euclidean_distances

    # Formatting the data
    data = np.array(zip(mds["c1"], mds["c2"]))

    # Configuring and running the KMeans
    k_means = KMeans(init=centers, n_clusters=3, n_init=1)
    k_means.fit_predict(data)

    # Creating the figure and axes
    fig_before, axe_before = plt.subplots(1, 1)
    fig_after, axe_after = plt.subplots(1, 1)
    fig_outliers, axe_outliers = plt.subplots(1, 1)

    # Setting the axe
    axe_before.xaxis.set_ticks_position("bottom")
    axe_before.yaxis.set_ticks_position("left")
    axe_before.spines["top"].set_visible(False)
    axe_before.spines["right"].set_visible(False)
    axe_before.spines["bottom"].set_position(("outward", 9))
    axe_before.spines["left"].set_position(("outward", 9))
    axe_after.xaxis.set_ticks_position("bottom")
    axe_after.yaxis.set_ticks_position("left")
    axe_after.spines["top"].set_visible(False)
    axe_after.spines["right"].set_visible(False)
    axe_after.spines["bottom"].set_position(("outward", 9))
    axe_after.spines["left"].set_position(("outward", 9))
    axe_outliers.xaxis.set_ticks_position("bottom")
    axe_outliers.yaxis.set_ticks_position("left")
    axe_outliers.spines["top"].set_visible(False)
    axe_outliers.spines["right"].set_visible(False)
    axe_outliers.spines["bottom"].set_position(("outward", 9))
    axe_outliers.spines["left"].set_position(("outward", 9))

    # Setting the title and labels
    axe_before.set_title("Before finding outliers", weight="bold")
    axe_before.set_xlabel(options.xaxis)
    axe_before.set_ylabel(options.yaxis)
    axe_after.set_title("After finding outliers\n($> {} "
                        "\\sigma$)".format(options.multiplier), weight="bold")
    axe_after.set_xlabel(options.xaxis)
    axe_after.set_ylabel(options.yaxis)
    axe_outliers.set_title("Outliers", weight="bold")
    axe_outliers.set_xlabel(options.xaxis)
    axe_outliers.set_ylabel(options.yaxis)

    # The population name
    ref_pop_name = ["CEU", "YRI", "JPT-CHB"]

    # The colors
    colors = ["#CC0000", "#669900", "#0099CC"]
    outlier_colors = ["#FFCACA", "#E2F0B6", "#C5EAF8"]

    # Plotting each of the clusters with the initial center
    plots_before = []
    plots_after = []
    plots_outliers = []
    plot_outliers = None
    plot_not_outliers = None
    outliers_set = set()
    for label in xrange(3):
        # Subsetting the data
        subset_mds = mds[k_means.labels_ == label]
        subset_data = data[k_means.labels_ == label]

        # Plotting the cluster
        p, = axe_before.plot(subset_mds["c1"], subset_mds["c2"], "o",
                             mec=colors[label], mfc=colors[label], ms=2,
                             clip_on=False)
        plots_before.append(p)

        # Plotting the cluster center (the real one)
        axe_before.plot(centers[label][0], centers[label][1], "o",
                        mec="#000000", mfc="#FFBB33", ms=6, clip_on=False)

        # Computing the distances
        distances = euclidean_distances(subset_data, centers[label])

        # Finding the outliers (that are not in the reference populations
        sigma = np.sqrt(np.true_divide(np.sum(distances ** 2),
                                       len(distances) - 1))
        outliers = np.logical_and(
            (distances > options.multiplier * sigma).flatten(),
            subset_mds["pop"] != ref_pop_name[label],
        )
        logger.info("  - {} outliers for the {} cluster".format(
            np.sum(outliers),
            ref_pop_name[label],
        ))

        # Saving the outliers
        if ref_pop_name[label] != options.outliers_of:
            # This is not the population we want, hence everybody is an outlier
            # (we don't include the reference population).
            outlier_mds = subset_mds[subset_mds["pop"] != ref_pop_name[label]]
            outliers_set |= set([(i["fid"], i["iid"]) for i in outlier_mds])

            # Plotting all samples that are not part of the reference
            # populations
            axe_outliers.plot(
                subset_mds["c1"][subset_mds["pop"] != ref_pop_name[label]],
                subset_mds["c2"][subset_mds["pop"] != ref_pop_name[label]],
                "o",
                mec="#555555",
                mfc="#555555",
                ms=2,
                clip_on=False,
            )
        else:
            # This is the population we want, hence only the real outliers are
            # outliers (we don't include the reference population)
            outlier_mds = subset_mds[
                np.logical_and(subset_mds["pop"] != ref_pop_name[label],
                               outliers)
            ]

            # Plotting the outliers
            plot_outliers, = axe_outliers.plot(outlier_mds["c1"],
                                               outlier_mds["c2"], "o",
                                               mec="#555555", mfc="#555555",
                                               ms=2, clip_on=False)

            # Plotting the not outliers
            plot_not_outliers, = axe_outliers.plot(
                subset_mds["c1"][np.logical_and(
                        ~outliers,
                        subset_mds["pop"] != ref_pop_name[label]
                )],
                subset_mds["c2"][np.logical_and(
                    ~outliers,
                    subset_mds["pop"] != ref_pop_name[label]
                )],
                "o",
                mec="#FFBB33",
                mfc="#FFBB33",
                ms=2,
                clip_on=False,
            )
            outliers_set |= set([(i["fid"], i["iid"]) for i in outlier_mds])

        # Plotting the cluster (without outliers)
        p, = axe_after.plot(
            subset_mds[~outliers]["c1"],
            subset_mds[~outliers]["c2"],
            "o",
            mec=colors[label],
            mfc=colors[label],
            ms=2,
            clip_on=False,
        )
        plots_after.append(p)

        # Plotting the cluster (only outliers)
        axe_after.plot(subset_mds[outliers]["c1"], subset_mds[outliers]["c2"],
                       "o", mec=outlier_colors[label],
                       mfc=outlier_colors[label], ms=2, clip_on=False)

        # Plotting only the reference populations
        p, = axe_outliers.plot(
            subset_mds["c1"][subset_mds["pop"] == ref_pop_name[label]],
            subset_mds["c2"][subset_mds["pop"] == ref_pop_name[label]],
            "o",
            mec=colors[label],
            mfc=colors[label],
            ms=2,
            clip_on=False,
        )
        plots_outliers.append(p)

        # Plotting the cluster center (the real one)
        axe_after.plot(centers[label][0], centers[label][1], "o",
                       mec="#000000", mfc="#FFBB33", ms=6)

    # The legends
    axe_before.legend(plots_before, ref_pop_name, loc="best", numpoints=1,
                      fancybox=True, fontsize=8).get_frame().set_alpha(0.5)
    axe_after.legend(plots_after, ref_pop_name, loc="best", numpoints=1,
                     fancybox=True, fontsize=8).get_frame().set_alpha(0.5)
    axe_outliers.legend(plots_outliers + [plot_not_outliers, plot_outliers],
                        ref_pop_name + ["SOURCE", "OUTLIERS"], loc="best",
                        numpoints=1, fancybox=True,
                        fontsize=8).get_frame().set_alpha(0.5)

    # Saving the figure
    fig_before.savefig("{}.before.{}".format(options.out, options.format),
                       dpi=300)
    fig_after.savefig("{}.after.{}".format(options.out, options.format),
                      dpi=300)
    fig_outliers.savefig("{}.outliers.{}".format(options.out, options.format),
                         dpi=300)

    return outliers_set