# Copyright (c) Microsoft Corporation and contributors.
# Licensed under the MIT License.

import numpy as np
from sklearn.utils.validation import check_is_fitted

from ..utils import import_graph, is_almost_symmetric
from .base import BaseEmbedMulti
from .svd import select_dimension, selectSVD


class MultipleASE(BaseEmbedMulti):
    r"""
    Multiple Adjacency Spectral Embedding (MASE) embeds arbitrary number of input
    graphs with matched vertex sets.

    For a population of undirected graphs, MASE assumes that the population of graphs
    is sampled from :math:`VR^{(i)}V^T` where :math:`V \in \mathbb{R}^{n\times d}` and
    :math:`R^{(i)} \in \mathbb{R}^{d\times d}`. Score matrices, :math:`R^{(i)}`, are
    allowed to vary for each graph, but are symmetric. All graphs share a common a
    latent position matrix :math:`V`.

    For a population of directed graphs, MASE assumes that the population is sampled
    from :math:`UR^{(i)}V^T` where :math:`U \in \mathbb{R}^{n\times d_1}`,
    :math:`V \in \mathbb{R}^{n\times d_2}`, and
    :math:`R^{(i)} \in \mathbb{R}^{d_1\times d_2}`. In this case, score matrices
    :math:`R^{(i)}` can be assymetric and non-square, but all graphs still share a
    common latent position matrices :math:`U` and :math:`V`.

    Parameters
    ----------
    n_components : int or None, default = None
        Desired dimensionality of output data. If "full",
        ``n_components`` must be ``<= min(X.shape)``. Otherwise, ``n_components`` must be
        ``< min(X.shape)``. If None, then optimal dimensions will be chosen by
        :func:`~graspologic.embed.select_dimension` using ``n_elbows`` argument.

    n_elbows : int, optional, default: 2
        If ``n_components`` is None, then compute the optimal embedding dimension using
        :func:`~graspologic.embed.select_dimension`. Otherwise, ignored.

    algorithm : {'randomized' (default), 'full', 'truncated'}, optional
        SVD solver to use:

        - 'randomized'
            Computes randomized svd using
            :func:`sklearn.utils.extmath.randomized_svd`
        - 'full'
            Computes full svd using :func:`scipy.linalg.svd`
        - 'truncated'
            Computes truncated svd using :func:`scipy.sparse.linalg.svds`

    n_iter : int, optional (default = 5)
        Number of iterations for randomized SVD solver. Not used by 'full' or
        'truncated'. The default is larger than the default in randomized_svd
        to handle sparse matrices that may have large slowly decaying spectrum.

    scaled : bool, optional (default=True)
        Whether to scale individual eigenvectors with eigenvalues in first embedding
        stage.

    diag_aug : bool, optional (default = True)
        Whether to replace the main diagonal of each adjacency matrices with
        a vector corresponding to the degree (or sum of edge weights for a
        weighted network) before embedding.

    concat : bool, optional (default False)
        If graph(s) are directed, whether to concatenate each graph's left and right (out and in) latent positions
        along axis 1.


    Attributes
    ----------
    n_graphs_ : int
        Number of graphs

    n_vertices_ : int
        Number of vertices in each graph

    latent_left_ : array, shape (n_samples, n_components)
        Estimated left latent positions of the graph.

    latent_right_ : array, shape (n_samples, n_components), or None
        Estimated right latent positions of the graph. Only computed when the an input
        graph is directed, or adjacency matrix is assymetric. Otherwise, None.

    scores_ : array, shape (n_samples, n_components, n_components)
        Estimated :math:`\hat{R}` matrices for each input graph.

    singular_values_ : array, shape (n_components) OR length 2 tuple of arrays
        If input graph is undirected, equal to the singular values of the concatenated
        adjacency spectral embeddings. If input graph is directed, :attr:`singular_values_`
        is a tuple of length 2, where :attr:`singular_values_[0]` corresponds to
        the singular values of the concatenated left adjacency spectral embeddings,
        and :attr:`singular_values_[1]` corresponds to
        the singular values of the concatenated right adjacency spectral embeddings.

    Notes
    -----
    When an input graph is directed, ``n_components`` of :attr:`latent_left_` may not be equal
    to ``n_components`` of :attr:`latent_right_`.
    """

    def __init__(
        self,
        n_components=None,
        n_elbows=2,
        algorithm="randomized",
        n_iter=5,
        scaled=True,
        diag_aug=True,
        concat=False,
    ):
        if not isinstance(scaled, bool):
            msg = "scaled must be a boolean, not {}".format(scaled)
            raise TypeError(msg)

        super().__init__(
            n_components=n_components,
            n_elbows=n_elbows,
            algorithm=algorithm,
            n_iter=n_iter,
            diag_aug=diag_aug,
            concat=concat,
        )
        self.scaled = scaled

    def _reduce_dim(self, graphs):
        # first embed into log2(n_vertices) for each graph
        n_components = int(np.ceil(np.log2(np.min(self.n_vertices_))))

        # embed individual graphs
        embeddings = [
            selectSVD(
                graph,
                n_components=n_components,
                algorithm=self.algorithm,
                n_iter=self.n_iter,
            )
            for graph in graphs
        ]
        Us, Ds, Vs = zip(*embeddings)

        # Choose the best embedding dimension for each graphs
        if self.n_components is None:
            embedding_dimensions = []
            for D in Ds:
                elbows, _ = select_dimension(D, n_elbows=self.n_elbows)
                embedding_dimensions.append(elbows[-1])

            # Choose the max of all of best embedding dimension of all graphs
            best_dimension = int(np.ceil(np.max(embedding_dimensions)))
        else:
            best_dimension = self.n_components

        if not self.scaled:
            Us = np.hstack([U[:, :best_dimension] for U in Us])
            Vs = np.hstack([V.T[:, :best_dimension] for V in Vs])
        else:
            # Equivalent to ASE
            Us = np.hstack(
                [
                    U[:, :best_dimension] @ np.diag(np.sqrt(D[:best_dimension]))
                    for U, D in zip(Us, Ds)
                ]
            )
            Vs = np.hstack(
                [
                    V.T[:, :best_dimension] @ np.diag(np.sqrt(D[:best_dimension]))
                    for V, D in zip(Vs, Ds)
                ]
            )

        # Second SVD for vertices
        # The notation is slightly different than the paper
        Uhat, sing_vals_left, _ = selectSVD(
            Us,
            n_components=self.n_components,
            n_elbows=self.n_elbows,
            algorithm=self.algorithm,
            n_iter=self.n_iter,
        )

        Vhat, sing_vals_right, _ = selectSVD(
            Vs,
            n_components=self.n_components,
            n_elbows=self.n_elbows,
            algorithm=self.algorithm,
            n_iter=self.n_iter,
        )
        return Uhat, Vhat, sing_vals_left, sing_vals_right

    def fit(self, graphs, y=None):
        """
        Fit the model with graphs.

        Parameters
        ----------
        graphs : list of nx.Graph or ndarray, or ndarray
            If list of nx.Graph, each Graph must contain same number of nodes.
            If list of ndarray, each array must have shape (n_vertices, n_vertices).
            If ndarray, then array must have shape (n_graphs, n_vertices, n_vertices).

        Returns
        -------
        self : object
            Returns an instance of self.
        """
        graphs = self._check_input_graphs(graphs)

        # Check if undirected
        undirected = all(is_almost_symmetric(g) for g in graphs)

        # Diag augment
        if self.diag_aug:
            graphs = self._diag_aug(graphs)

        # embed
        Uhat, Vhat, sing_vals_left, sing_vals_right = self._reduce_dim(graphs)
        self.latent_left_ = Uhat
        if not undirected:
            self.latent_right_ = Vhat
            self.scores_ = Uhat.T @ graphs @ Vhat
            self.singular_values_ = (sing_vals_left, sing_vals_right)
        else:
            self.latent_right_ = None
            self.scores_ = Uhat.T @ graphs @ Uhat
            self.singular_values_ = sing_vals_left

        return self

    def fit_transform(self, graphs, y=None):
        """
        Fit the model with graphs and apply the embedding on graphs.
        n_components is either automatically determined or based on user input.

        Parameters
        ----------
        graphs : list of nx.Graph or ndarray, or ndarray
            If list of nx.Graph, each Graph must contain same number of nodes.
            If list of ndarray, each array must have shape (n_vertices, n_vertices).
            If ndarray, then array must have shape (n_graphs, n_vertices, n_vertices).

        Returns
        -------
        out : np.ndarray or length 2 tuple of np.ndarray.
            If input graphs were symmetric shape (n_vertices, n_components).
            If graphs were directed and ``concat`` is False, returns tuple of two arrays (same shape as above).
            The first corresponds to the left latent positions, and the second to the right latent positions.
            When ``concat`` is True left and right (out and in) latent positions are concatenated along axis 1.
        """
        return self._fit_transform(graphs)

graspologic/embed/mase.py
# Copyright (c) Microsoft Corporation and contributors.
# Licensed under the MIT License.

import numpy as np
from sklearn.utils.validation import check_is_fitted

from ..utils import import_graph, is_almost_symmetric
from .base import BaseEmbedMulti
from .svd import select_dimension, selectSVD


class MultipleASE(BaseEmbedMulti):
    r"""
    Multiple Adjacency Spectral Embedding (MASE) embeds arbitrary number of input
    graphs with matched vertex sets.

    For a population of undirected graphs, MASE assumes that the population of graphs
    is sampled from :math:`VR^{(i)}V^T` where :math:`V \in \mathbb{R}^{n\times d}` and
    :math:`R^{(i)} \in \mathbb{R}^{d\times d}`. Score matrices, :math:`R^{(i)}`, are
    allowed to vary for each graph, but are symmetric. All graphs share a common a
    latent position matrix :math:`V`.

    For a population of directed graphs, MASE assumes that the population is sampled
    from :math:`UR^{(i)}V^T` where :math:`U \in \mathbb{R}^{n\times d_1}`,
    :math:`V \in \mathbb{R}^{n\times d_2}`, and
    :math:`R^{(i)} \in \mathbb{R}^{d_1\times d_2}`. In this case, score matrices
    :math:`R^{(i)}` can be assymetric and non-square, but all graphs still share a
    common latent position matrices :math:`U` and :math:`V`.

    Read more in the :ref:`tutorials <embed_tutorials>`

    Parameters
    ----------
    n_components : int or None, default = None
        Desired dimensionality of output data. If "full",
        ``n_components`` must be ``<= min(X.shape)``. Otherwise, ``n_components`` must be
        ``< min(X.shape)``. If None, then optimal dimensions will be chosen by
        :func:`~graspologic.embed.select_dimension` using ``n_elbows`` argument.

    n_elbows : int, optional, default: 2
        If ``n_components`` is None, then compute the optimal embedding dimension using
        :func:`~graspologic.embed.select_dimension`. Otherwise, ignored.

    algorithm : {'randomized' (default), 'full', 'truncated'}, optional
        SVD solver to use:

        - 'randomized'
            Computes randomized svd using
            :func:`sklearn.utils.extmath.randomized_svd`
        - 'full'
            Computes full svd using :func:`scipy.linalg.svd`
        - 'truncated'
            Computes truncated svd using :func:`scipy.sparse.linalg.svds`

    n_iter : int, optional (default = 5)
        Number of iterations for randomized SVD solver. Not used by 'full' or
        'truncated'. The default is larger than the default in randomized_svd
        to handle sparse matrices that may have large slowly decaying spectrum.

    scaled : bool, optional (default=True)
        Whether to scale individual eigenvectors with eigenvalues in first embedding
        stage.

    diag_aug : bool, optional (default = True)
        Whether to replace the main diagonal of each adjacency matrices with
        a vector corresponding to the degree (or sum of edge weights for a
        weighted network) before embedding.

    concat : bool, optional (default False)
        If graph(s) are directed, whether to concatenate each graph's left and right (out and in) latent positions
        along axis 1.


    Attributes
    ----------
    n_graphs_ : int
        Number of graphs

    n_vertices_ : int
        Number of vertices in each graph

    latent_left_ : array, shape (n_samples, n_components)
        Estimated left latent positions of the graph.

    latent_right_ : array, shape (n_samples, n_components), or None
        Estimated right latent positions of the graph. Only computed when the an input
        graph is directed, or adjacency matrix is assymetric. Otherwise, None.

    scores_ : array, shape (n_samples, n_components, n_components)
        Estimated :math:`\hat{R}` matrices for each input graph.

    singular_values_ : array, shape (n_components) OR length 2 tuple of arrays
        If input graph is undirected, equal to the singular values of the concatenated
        adjacency spectral embeddings. If input graph is directed, :attr:`singular_values_`
        is a tuple of length 2, where :attr:`singular_values_[0]` corresponds to
        the singular values of the concatenated left adjacency spectral embeddings,
        and :attr:`singular_values_[1]` corresponds to
        the singular values of the concatenated right adjacency spectral embeddings.

    Notes
    -----
    When an input graph is directed, ``n_components`` of :attr:`latent_left_` may not be equal
    to ``n_components`` of :attr:`latent_right_`.
    """

    def __init__(
        self,
        n_components=None,
        n_elbows=2,
        algorithm="randomized",
        n_iter=5,
        scaled=True,
        diag_aug=True,
        concat=False,
    ):
        if not isinstance(scaled, bool):
            msg = "scaled must be a boolean, not {}".format(scaled)
            raise TypeError(msg)

        super().__init__(
            n_components=n_components,
            n_elbows=n_elbows,
            algorithm=algorithm,
            n_iter=n_iter,
            diag_aug=diag_aug,
            concat=concat,
        )
        self.scaled = scaled

    def _reduce_dim(self, graphs):
        # first embed into log2(n_vertices) for each graph
        n_components = int(np.ceil(np.log2(np.min(self.n_vertices_))))

        # embed individual graphs
        embeddings = [
            selectSVD(
                graph,
                n_components=n_components,
                algorithm=self.algorithm,
                n_iter=self.n_iter,
            )
            for graph in graphs
        ]
        Us, Ds, Vs = zip(*embeddings)

        # Choose the best embedding dimension for each graphs
        if self.n_components is None:
            embedding_dimensions = []
            for D in Ds:
                elbows, _ = select_dimension(D, n_elbows=self.n_elbows)
                embedding_dimensions.append(elbows[-1])

            # Choose the max of all of best embedding dimension of all graphs
            best_dimension = int(np.ceil(np.max(embedding_dimensions)))
        else:
            best_dimension = self.n_components

        if not self.scaled:
            Us = np.hstack([U[:, :best_dimension] for U in Us])
            Vs = np.hstack([V.T[:, :best_dimension] for V in Vs])
        else:
            # Equivalent to ASE
            Us = np.hstack(
                [
                    U[:, :best_dimension] @ np.diag(np.sqrt(D[:best_dimension]))
                    for U, D in zip(Us, Ds)
                ]
            )
            Vs = np.hstack(
                [
                    V.T[:, :best_dimension] @ np.diag(np.sqrt(D[:best_dimension]))
                    for V, D in zip(Vs, Ds)
                ]
            )

        # Second SVD for vertices
        # The notation is slightly different than the paper
        Uhat, sing_vals_left, _ = selectSVD(
            Us,
            n_components=self.n_components,
            n_elbows=self.n_elbows,
            algorithm=self.algorithm,
            n_iter=self.n_iter,
        )

        Vhat, sing_vals_right, _ = selectSVD(
            Vs,
            n_components=self.n_components,
            n_elbows=self.n_elbows,
            algorithm=self.algorithm,
            n_iter=self.n_iter,
        )
        return Uhat, Vhat, sing_vals_left, sing_vals_right

    def fit(self, graphs, y=None):
        """
        Fit the model with graphs.

        Parameters
        ----------
        graphs : list of nx.Graph or ndarray, or ndarray
            If list of nx.Graph, each Graph must contain same number of nodes.
            If list of ndarray, each array must have shape (n_vertices, n_vertices).
            If ndarray, then array must have shape (n_graphs, n_vertices, n_vertices).

        Returns
        -------
        self : object
            Returns an instance of self.
        """
        graphs = self._check_input_graphs(graphs)

        # Check if undirected
        undirected = all(is_almost_symmetric(g) for g in graphs)

        # Diag augment
        if self.diag_aug:
            graphs = self._diag_aug(graphs)

        # embed
        Uhat, Vhat, sing_vals_left, sing_vals_right = self._reduce_dim(graphs)
        self.latent_left_ = Uhat
        if not undirected:
            self.latent_right_ = Vhat
            self.scores_ = Uhat.T @ graphs @ Vhat
            self.singular_values_ = (sing_vals_left, sing_vals_right)
        else:
            self.latent_right_ = None
            self.scores_ = Uhat.T @ graphs @ Uhat
            self.singular_values_ = sing_vals_left

        return self

    def fit_transform(self, graphs, y=None):
        """
        Fit the model with graphs and apply the embedding on graphs.
        n_components is either automatically determined or based on user input.

        Parameters
        ----------
        graphs : list of nx.Graph or ndarray, or ndarray
            If list of nx.Graph, each Graph must contain same number of nodes.
            If list of ndarray, each array must have shape (n_vertices, n_vertices).
            If ndarray, then array must have shape (n_graphs, n_vertices, n_vertices).

        Returns
        -------
        out : np.ndarray or length 2 tuple of np.ndarray.
            If input graphs were symmetric shape (n_vertices, n_components).
            If graphs were directed and ``concat`` is False, returns tuple of two arrays (same shape as above).
            The first corresponds to the left latent positions, and the second to the right latent positions.
            When ``concat`` is True left and right (out and in) latent positions are concatenated along axis 1.
        """
        return self._fit_transform(graphs)

graspologic/embed/mase.py
diff --git a/docs/tutorial.rst b/docs/tutorial.rst
index 104580b46..467f55518 100644
--- a/docs/tutorial.rst
+++ b/docs/tutorial.rst
@@ -41,6 +41,7 @@ Inference on random graphs depends on low-dimensional Euclidean representation o
    
    tutorials/embedding/AdjacencySpectralEmbed
    tutorials/embedding/Omnibus
+   tutorials/embedding/MASE
    
 .. _inference_tutorials: 
 
diff --git a/docs/tutorials/embedding/MASE.ipynb b/docs/tutorials/embedding/MASE.ipynb
new file mode 100644
index 000000000..149c26051
--- /dev/null
+++ b/docs/tutorials/embedding/MASE.ipynb
@@ -0,0 +1,380 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Multiple Adjacency Spectral Embedding (MASE)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Multiple Adjacency Spectral Embedding (MASE) is an extension of Adjacency Spectral Embedding (see ASE [tutorial](https://graspy.neurodata.io/tutorials/embedding/adjacencyspectralembed)) for an arbitrary number of graphs.  Once graphs are embedded, the low-dimensional Euclidean representation can be used to visualize the latent positions of vertices, perform inference, etc."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "%matplotlib inline"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Replicating Omnibus Embedding Tutorial"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "In order to demonstrate how to use the ``MultipleASE`` class, we will use a simple example: two 2-block stochastic block models (SBMs) with different block probabilities.  Indeed, we will use the same models as in the omnibus [tutorial](https://graspy.neurodata.io/tutorials/embedding/omnibus) to facilitate direct comparison. "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Simulate two different graphs using stochastic block models (SBM)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We sample two 2-block SBMs (undirected, no self-loops) with 50 vertices, each block containing 25 vertices ($n = [25, 25]$), and with the following block connectivity matrices:\n",
+    "\n",
+    "\\begin{align*}\n",
+    "B_1 = \n",
+    "\\begin{bmatrix}0.3 & 0.1\\\\\n",
+    "0.1 & 0.7\n",
+    "\\end{bmatrix},~\n",
+    "B_2 = \\begin{bmatrix}0.3 & 0.1\\\\\n",
+    "0.1 & 0.3\n",
+    "\\end{bmatrix}\n",
+    "\\end{align*}\n",
+    "\n",
+    "The only difference between the two is the within-block probability for the second block. We sample $G_1 \\sim \\text{SBM}(n, B_1)$ and $G_2 \\sim \\text{SBM}(n, B_2)$."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from graspologic.simulations import sbm\n",
+    "\n",
+    "n = [25, 25]\n",
+    "B1 = [[0.3, 0.1],\n",
+    "      [0.1, 0.7]]\n",
+    "B2 = [[0.3, 0.1],\n",
+    "      [0.1, 0.3]]\n",
+    "\n",
+    "np.random.seed(8)\n",
+    "G1 = sbm(n, B1)\n",
+    "G2 = sbm(n, B2)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Visualize the graphs using heatmap\n",
+    "\n",
+    "We visualize the sampled graphs using the ``heatmap()`` function. ``heatmap()`` will plot the adjacency matrix, where the colors represent the weight of each edge. In this case, we have binary graphs so the values will be either 0 or 1. \n",
+    "\n",
+    "We see that there is clear block structure to the graphs. Furthermore, the lower right quarter of $G_1$, representing the within-group connections for the second group, is more dense than that of $G_2$."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from graspologic.plot import heatmap\n",
+    "\n",
+    "heatmap(G1, figsize=(7,7), title=\"Visualization of Graph 1\", cbar=False)\n",
+    "_ = heatmap(G2, figsize=(7,7), title=\"Visualization of Graph 2\", cbar=False)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Embed the two graphs using MASE"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Just as ASE fits a single input graph to the Random Dot Product Graph (RDPG) model, MASE fits many input graphs to the common subspace independent-edge (COSIE) model, which is an extension of the RDPG model to multiple graphs. The COSIE model assumes that all of the graphs have important shared properties, but can also have meaningful differences.  As a result, COSIE decomposes the set of graphs into a set of shared latent positions $V$ that describe similarities, and a set of score matrices $R^{(i)}$ that describe how each individual graph is different.  Mathematically, for the adjacency matrix of an undirected graph $A^{(i)}$,\n",
+    "$$A^{(i)} = VR^{(i)}V^T$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "In general, 2 clusters can be faithfully represented in 2 dimensions, so we specify the number of components as 2. Then, the ``fit_transform()`` method returns the estimated latent positions $\\hat{V}$, which we can see has 50 entries corresponding to the 50 vertices embedded in 2 dimensions."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from graspologic.embed import MultipleASE as MASE\n",
+    "\n",
+    "embedder = MASE(n_components=2)\n",
+    "V_hat = embedder.fit_transform([G1, G2])\n",
+    "print(V_hat.shape)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Visualize Common Latent Positions"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Since the two graphs have clear block structures, with higher within-block probabilities than between-block, we see two distinct \"clusters\" when we visualize the shared latent positions."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "def vis_latent(V, title, predicted=False, pred_labels=None):\n",
+    "    fig, ax = plt.subplots(figsize=(10,10))\n",
+    "    if predicted:\n",
+    "        ax.scatter(V_hat[pred_labels==0, 0], V_hat[pred_labels==0, 1], c=\"blue\", label=\"Predicted Block 1\")\n",
+    "        ax.scatter(V_hat[pred_labels==1, 0], V_hat[pred_labels==1, 1], c=\"red\", label=\"Predicted Block 2\")\n",
+    "    else:\n",
+    "        ax.scatter(V[:25, 0], V[:25, 1], c=\"blue\", label=\"Block 1\")\n",
+    "        ax.scatter(V[25:, 0], V[25:, 1], c=\"red\", label=\"Block 2\")\n",
+    "        \n",
+    "    ax.legend(prop={'size':20})\n",
+    "    ax.set_xticklabels([])\n",
+    "    ax.set_yticklabels([])\n",
+    "    ax.set_xticks([])\n",
+    "    ax.set_yticks([])\n",
+    "    _ = ax.set_title(title, fontsize=20)\n",
+    "\n",
+    "vis_latent(V_hat, \"Common Latent Positions from MASE Embedding\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Recover Original Communities"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We can apply a simple clustering algorithm (kmeans) to recover the original communities."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#Cluster V_hat\n",
+    "from graspologic.cluster import KMeansCluster\n",
+    "clusterer = KMeansCluster()\n",
+    "pred_labels = clusterer.fit_predict(V_hat)\n",
+    "\n",
+    "#Remap labels\n",
+    "from graspologic.utils import remap_labels\n",
+    "true_labels = [0]*25 + [1]*25\n",
+    "remapped_labels = remap_labels(true_labels, pred_labels)\n",
+    "\n",
+    "#Visualize latent positions\n",
+    "vis_latent(V_hat, \"Predicted Communities from Common Latent Positions\", True, remapped_labels)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Even with only one graph sampled from each SBM, we can perfectly recover the two communities!"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Sampling More Graphs"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We can recover a more accurate representation of the two SBM populations with more sample graphs.  Let's try with 100 graphs sampled from each."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "G1 = [sbm(n, B1) for i in range(100)]\n",
+    "G2 = [sbm(n, B2) for i in range(100)]\n",
+    "\n",
+    "V_hat = embedder.fit_transform(G1 + G2)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Visualizing Common Latent Subspace"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now, the clustering is even more distinct, as we would expect with more samples."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "vis_latent(V_hat, \"Common Latent Subspace from MASE Embedding with Many Sampled Graphs\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Visualize Score Matrices"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now that we have more samples, we can apply the same process to visualize the individual score matrices.  Note, since the score matrices are symmetric, $2 \\times 2$ matrices, we only need a 3-dimensional vector to represent each.  In this case, we can clearly see two distinct clusters: one for the score matrices from the first SBM population another for the score matrices from the second."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from graspologic.plot import pairplot\n",
+    "\n",
+    "#Extract Score Matrices (in vector form)\n",
+    "R_1_hat = np.reshape(embedder.scores_[:100], (100, 4))\n",
+    "R_2_hat = np.reshape(embedder.scores_[100:], (100, 4))\n",
+    "\n",
+    "#Omitting redundant entries due to symmetry\n",
+    "R_1_hat = np.delete(R_1_hat, 1, axis=1)\n",
+    "R_2_hat = np.delete(R_2_hat, 1, axis=1)\n",
+    "\n",
+    "#Combining into a single array\n",
+    "R_hat = np.concatenate((R_1_hat, R_2_hat), axis=0)\n",
+    "\n",
+    "#Visualize\n",
+    "labels = [\"SBM 1\"]*100 + [\"SBM 2\"]*100\n",
+    "_ = pairplot(R_hat, labels, title=\"Score Matrices from MASE Embedding\",\n",
+    "             palette={\"green\", \"purple\"}, legend_name=\"Population\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Recover Original SBM Membership"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "As before, we can apply a simple clustering algorithm (kmeans) to recover from which SBM population each graph was generated."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#Cluster R_hat\n",
+    "pred_labels = clusterer.fit_predict(R_hat)\n",
+    "true_labels = [0]*100 + [1]*100\n",
+    "remapped_labels = remap_labels(true_labels, pred_labels)\n",
+    "pred_sbm = []\n",
+    "for label in remapped_labels:\n",
+    "    if label == 0:\n",
+    "        pred_sbm.append(\"SBM 1\")\n",
+    "    else:\n",
+    "        pred_sbm.append(\"SBM 2\")\n",
+    "\n",
+    "#Visualize Results\n",
+    "_ = pairplot(R_hat, pred_sbm, title=\"Predicted SBM Membership from Score Matrices\",\n",
+    "             palette={\"green\", \"purple\"}, legend_name=\"Population\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "With this method we can perfectly recover the original SBM population to which each graph belongs!"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.7.0"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/docs/tutorials/embedding/Omnibus.ipynb b/docs/tutorials/embedding/Omnibus.ipynb
index 0581b9e40..ba93fef5c 100644
--- a/docs/tutorials/embedding/Omnibus.ipynb
+++ b/docs/tutorials/embedding/Omnibus.ipynb
@@ -4,7 +4,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# Embedding Multiple Graphs\n",
+    "# Omnibus Embedding for Multiple Graphs\n",
     "\n",
     "This demo shows you how to simultaneously embed two graphs using omnibus embedding from two graphs sampled from different stochastic block models (SBM). We will also compare the results to that of adjacency spectral embedding, and show why it is useful to embed the graphs simultaneously."
    ]
diff --git a/graspologic/embed/mase.py b/graspologic/embed/mase.py
index d14cc4f9e..17e7f27e4 100644
--- a/graspologic/embed/mase.py
+++ b/graspologic/embed/mase.py
@@ -27,6 +27,8 @@ class MultipleASE(BaseEmbedMulti):
     :math:`R^{(i)}` can be assymetric and non-square, but all graphs still share a
     common latent position matrices :math:`U` and :math:`V`.
 
+    Read more in the :ref:`tutorials <embed_tutorials>`
+
     Parameters
     ----------
     n_components : int or None, default = None

Add tutorial for MASE

