def _get_relative_locations(shapes, graph, level_str, verbose):
    r"""
    returns numpy.array of size 2 x n_images x n_edges
    """
    # convert given shapes to point graphs
    if isinstance(graph, Tree):
        point_graphs = [PointTree(shape.points, graph.adjacency_array,
                                  graph.root_vertex) for shape in shapes]
    else:
        point_graphs = [PointDirectedGraph(shape.points, graph.adjacency_array)
                        for shape in shapes]

    # initialize an output numpy array
    rel_loc_array = np.empty((2, graph.n_edges, len(point_graphs)))

    # get relative locations
    for c, pt in enumerate(point_graphs):
        # print progress
        if verbose:
            print_dynamic('{}Computing relative locations from '
                          'shapes - {}'.format(
                          level_str,
                          progress_bar_str(float(c + 1) / len(point_graphs),
                                           show_bar=False)))

        # get relative locations from this shape
        rl = pt.relative_locations()

        # store
        rel_loc_array[..., c] = rl.T

    # rollaxis and return
    return np.rollaxis(rel_loc_array, 2, 1)

    def _regression_data(self, images, gt_shapes, perturbed_shapes,
                         verbose=False):
        r"""
        Method that generates the regression data : features and delta_ps.

        Parameters
        ----------
        images : list of :map:`MaskedImage`
            The set of landmarked images.

        gt_shapes : :map:`PointCloud` list
            List of the ground truth shapes that correspond to the images.

        perturbed_shapes : :map:`PointCloud` list
            List of the perturbed shapes in order to regress.

        verbose : `boolean`, optional
            If ``True``, the progress is printed.
        """
        if verbose:
            print_dynamic('- Generating regression data')

        n_images = len(images)
        features = []
        delta_ps = []
        for j, (i, s, p_shape) in enumerate(zip(images, gt_shapes,
                                                perturbed_shapes)):
            if verbose:
                print_dynamic('- Generating regression data - {}'.format(
                    progress_bar_str((j + 1.) / n_images, show_bar=False)))
            for ps in p_shape:
                features.append(self.features(i, ps))
                delta_ps.append(self.delta_ps(s, ps))
        return np.asarray(features), np.asarray(delta_ps)

def build_shape_model(shapes, max_components=None, prefix='', verbose=False):
    r"""
    Builds a shape model given a set of shapes.

    Parameters
    ----------
    shapes: list of :map:`PointCloud`
        The set of shapes from which to build the model.
    max_components: None or int or float
        Specifies the number of components of the trained shape model.
        If int, it specifies the exact number of components to be retained.
        If float, it specifies the percentage of variance to be retained.
        If None, all the available components are kept (100% of variance).

    Returns
    -------
    shape_model: :class:`menpo.model.pca`
        The PCA shape model.
    """
    if verbose:
        print_dynamic('{}Building shape model'.format(prefix))
    # compute aligned shapes
    aligned_shapes = align_shapes(shapes)
    # build shape model
    shape_model = PCAModel(aligned_shapes)
    if max_components is not None:
        # trim shape model if required
        shape_model.trim_components(max_components)
    return shape_model

def _import_glob_generator(
    pattern,
    extension_map,
    max_assets=None,
    has_landmarks=False,
    landmark_resolver=None,
    importer_kwargs=None,
    verbose=False,
):
    filepaths = list(glob_with_suffix(pattern, extension_map))
    if max_assets:
        filepaths = filepaths[:max_assets]
    n_files = len(filepaths)
    if n_files == 0:
        raise ValueError("The glob {} yields no assets".format(pattern))
    for i, asset in enumerate(
        _multi_import_generator(
            filepaths,
            extension_map,
            has_landmarks=has_landmarks,
            landmark_resolver=landmark_resolver,
            importer_kwargs=importer_kwargs,
        )
    ):
        if verbose:
            print_dynamic(
                "- Loading {} assets: {}".format(n_files, progress_bar_str(float(i + 1) / n_files, show_bar=True))
            )
        yield asset

def test_model(model, test_images, num_init):
    face_detector = menpodetect.load_dlib_frontal_face_detector()
    test_gt_shapes = util.get_gt_shapes(test_images)
    test_boxes = util.get_bounding_boxes(test_images, test_gt_shapes, face_detector)

    initial_errors = []
    final_errors = []

    initial_shapes = []
    final_shapes = []

    for k, (im, gt_shape, box) in enumerate(zip(test_images, test_gt_shapes, test_boxes)):
        init_shapes, fin_shapes = model.apply(im, ([box], num_init, None))

        init_shape = util.get_median_shape(init_shapes)
        final_shape = fin_shapes[0]

        initial_shapes.append(init_shape)
        final_shapes.append(final_shape)

        initial_errors.append(compute_error(init_shape, gt_shape))
        final_errors.append(compute_error(final_shape, gt_shape))

        print_dynamic('{}/{}'.format(k + 1, len(test_images)))

    return initial_errors, final_errors, initial_shapes, final_shapes

def compute_reference_shape(shapes, diagonal, verbose=False):
    r"""
    Function that computes the reference shape as the mean shape of the provided
    shapes.

    Parameters
    ----------
    shapes : `list` of `menpo.shape.PointCloud`
        The set of shapes from which to build the reference shape.
    diagonal : `int` or ``None``
        If `int`, it ensures that the mean shape is scaled so that the diagonal
        of the bounding box containing it matches the provided value.
        If ``None``, then the mean shape is not rescaled.
    verbose : `bool`, optional
        If ``True``, then progress information is printed.

    Returns
    -------
    reference_shape : `menpo.shape.PointCloud`
        The reference shape.
    """
    # the reference_shape is the mean shape of the images' landmarks
    if verbose:
        print_dynamic('- Computing reference shape')
    reference_shape = mean_pointcloud(shapes)

    # fix the reference_shape's diagonal length if asked
    if diagonal:
        x, y = reference_shape.range()
        scale = diagonal / np.sqrt(x ** 2 + y ** 2)
        reference_shape = Scale(scale, reference_shape.n_dims).apply(
            reference_shape)

    return reference_shape

def apply_pyramid_on_images(generators, n_levels, verbose=False):
    r"""
    Exhausts the pyramid generators verbosely
    """
    all_images = []
    for j in range(n_levels):

        if verbose:
            level_str = '- Apply pyramid: '
            if n_levels > 1:
                level_str = '- Apply pyramid: [Level {} - '.format(j + 1)

        level_images = []
        for c, g in enumerate(generators):
            if verbose:
                print_dynamic(
                    '{}Computing feature space/rescaling - {}'.format(
                        level_str,
                        progress_bar_str((c + 1.) / len(generators),
                                         show_bar=False)))
            level_images.append(next(g))
        all_images.append(level_images)
    if verbose:
        print_dynamic('- Apply pyramid: Done\n')
    return all_images

def _build_appearance_model_sparse(all_patches_array, graph, patch_shape,
                                   n_channels, n_appearance_parameters,
                                   level_str, verbose):
    # build appearance model
    if verbose:
        print_dynamic('{}Training appearance distribution per '
                      'edge'.format(level_str))

    # compute mean appearance vector
    app_mean = np.mean(all_patches_array, axis=1)

    # appearance vector and patch vector lengths
    patch_len = np.prod(patch_shape) * n_channels

    # initialize block sparse covariance matrix
    all_cov = lil_matrix((graph.n_vertices * patch_len,
                          graph.n_vertices * patch_len))

    # compute covariance matrix for each edge
    for e in range(graph.n_edges):
        # print progress
        if verbose:
            print_dynamic('{}Training appearance distribution '
                          'per edge - {}'.format(
                          level_str,
                          progress_bar_str(float(e + 1) / graph.n_edges,
                                           show_bar=False)))

        # edge vertices
        v1 = np.min(graph.adjacency_array[e, :])
        v2 = np.max(graph.adjacency_array[e, :])

        # find indices in target covariance matrix
        v1_from = v1 * patch_len
        v1_to = (v1 + 1) * patch_len
        v2_from = v2 * patch_len
        v2_to = (v2 + 1) * patch_len

        # extract data
        edge_data = np.concatenate((all_patches_array[v1_from:v1_to, :],
                                    all_patches_array[v2_from:v2_to, :]))

        # compute covariance inverse
        icov = _covariance_matrix_inverse(np.cov(edge_data),
                                          n_appearance_parameters)

        # v1, v2
        all_cov[v1_from:v1_to, v2_from:v2_to] += icov[:patch_len, patch_len::]

        # v2, v1
        all_cov[v2_from:v2_to, v1_from:v1_to] += icov[patch_len::, :patch_len]

        # v1, v1
        all_cov[v1_from:v1_to, v1_from:v1_to] += icov[:patch_len, :patch_len]

        # v2, v2
        all_cov[v2_from:v2_to, v2_from:v2_to] += icov[patch_len::, patch_len::]

    return app_mean, all_cov.tocsr()

def train_aps(experiments_path, fast, group, training_images_options,
              training_options, save_model, verbose):
    # update training_images_options
    training_images_options['save_path'] = os.path.join(experiments_path,
                                                        'Databases')
    training_images_options['fast'] = fast
    training_images_options['group'] = group
    training_images_options['verbose'] = verbose

    # parse training options
    adj, rv = parse_deformation_graph(training_options['graph_deformation'])
    training_options['adjacency_array_deformation'] = adj
    training_options['root_vertex_deformation'] = rv
    adj, fl = parse_appearance_graph(training_options['graph_appearance'])
    training_options['adjacency_array_appearance'] = adj
    training_options['gaussian_per_patch'] = fl
    training_options['features'] = parse_features(training_options['features'],
                                                  fast)
    graph_deformation_str = training_options['graph_deformation']
    graph_appearance_str = training_options['graph_appearance']
    del training_options['graph_deformation']
    del training_options['graph_appearance']

    # Load training images
    training_images = load_database(**training_images_options)

    # make model filename
    filename = model_filename(training_images_options, training_options, group,
                              fast, graph_deformation_str, graph_appearance_str)
    save_path = os.path.join(experiments_path, 'Models', filename)

    # train model
    if file_exists(save_path):
        if verbose:
            print_dynamic('Loading model...')
        aps = pickle_load(save_path)
        if verbose:
            print_dynamic('Model loaded.')
    else:
        training_options['max_shape_components'] = None

        # Train model
        if fast:
            from antonakoscvpr2015.menpofast.builder import APSBuilder
        else:
            from antonakoscvpr2015.menpo.builder import APSBuilder
        if group is not None:
            aps = APSBuilder(**training_options).build(training_images,
                                                       group=group.__name__,
                                                       verbose=verbose)
        else:
            aps = APSBuilder(**training_options).build(training_images,
                                                       verbose=verbose)

        # save model
        if save_model:
            pickle_dump(aps, save_path)

    return aps, filename, training_images

def _build_deformation_model(graph, relative_locations, level_str, verbose):
    # build deformation model
    if verbose:
        print_dynamic('{}Training deformation distribution per '
                      'graph edge'.format(level_str))
    def_len = 2 * graph.n_vertices
    def_cov = np.zeros((def_len, def_len))
    for e in range(graph.n_edges):
        # print progress
        if verbose:
            print_dynamic('{}Training deformation distribution '
                          'per edge - {}'.format(
                          level_str,
                          progress_bar_str(float(e + 1) / graph.n_edges,
                                           show_bar=False)))

        # get vertices adjacent to edge
        parent = graph.adjacency_array[e, 0]
        child = graph.adjacency_array[e, 1]

        # compute covariance matrix
        edge_cov = np.linalg.inv(np.cov(relative_locations[..., e]))

        # store its values
        s1 = edge_cov[0, 0]
        s2 = edge_cov[1, 1]
        s3 = 2 * edge_cov[0, 1]

        # Fill the covariance matrix matrix
        # get indices
        p1 = 2 * parent
        p2 = 2 * parent + 1
        c1 = 2 * child
        c2 = 2 * child + 1

        # up-left block
        def_cov[p1, p1] += s1
        def_cov[p2, p2] += s2
        def_cov[p2, p1] += s3

        # up-right block
        def_cov[p1, c1] = - s1
        def_cov[p2, c2] = - s2
        def_cov[p1, c2] = - s3 / 2
        def_cov[p2, c1] = - s3 / 2

        # down-left block
        def_cov[c1, p1] = - s1
        def_cov[c2, p2] = - s2
        def_cov[c1, p2] = - s3 / 2
        def_cov[c2, p1] = - s3 / 2

        # down-right block
        def_cov[c1, c1] += s1
        def_cov[c2, c2] += s2
        def_cov[c1, c2] += s3

    return def_cov

    def _train(self, images, gt_shapes, current_shapes, increment=False,
               prefix='', verbose=False):

        if not increment:
            # Reset the regressors
            self.regressors = []

        n_perturbations = len(current_shapes[0])
        template_shape = gt_shapes[0]

        # obtain delta_x and gt_x (parameters rather than shapes)
        delta_x, gt_x = obtain_parametric_delta_x(gt_shapes, current_shapes,
                                                  self.transform)

        # Cascaded Regression loop
        for k in range(self.n_iterations):
            # generate regression data
            features = self._generate_features(
                images, current_shapes,
                prefix='{}(Iteration {}) - '.format(prefix, k),
                verbose=verbose)

            if verbose:
                print_dynamic('{}(Iteration {}) - Performing regression'.format(
                    prefix, k))

            if not increment:
                r = self._regressor_cls()
                r.train(features, delta_x)
                self.regressors.append(r)
            else:
                self.regressors[k].increment(features, delta_x)

            # Estimate delta_points
            estimated_delta_x = self.regressors[k].predict(features)
            if verbose:
                self._print_regression_info(template_shape, gt_shapes,
                                            n_perturbations, delta_x,
                                            estimated_delta_x, k,
                                            prefix=prefix)

            j = 0
            for shapes in current_shapes:
                for s in shapes:
                    # Estimate parameters
                    edx = estimated_delta_x[j]
                    # Current parameters
                    cx = _weights_for_target(self.transform, s) + edx

                    # Uses less memory to find updated target shape
                    self.transform.from_vector_inplace(cx)
                    # Update current shape inplace
                    s.from_vector_inplace(self.transform.target.as_vector())

                    delta_x[j] = gt_x[j] - cx
                    j += 1

        return current_shapes

    def _create_pyramid(cls, images, n_levels, downscale, pyramid_on_features,
                        feature_type, verbose=False):
        r"""
        Function that creates a generator function for Gaussian pyramid. The
        pyramid can be created either on the feature space or the original
        (intensities) space.

        Parameters
        ----------
        images: list of :class:`menpo.image.Image`
            The set of landmarked images from which to build the AAM.
        n_levels: int
            The number of multi-resolution pyramidal levels to be used.
        downscale: float
            The downscale factor that will be used to create the different
            pyramidal levels.
        pyramid_on_features: boolean
            If True, the features are extracted at the highest level and the
            pyramid is created on the feature images.
            If False, the pyramid is created on the original (intensities)
            space.
        feature_type: list of size 1 with str or function/closure or None
            The feature type to be used in case pyramid_on_features is enabled.
        verbose: bool, Optional
            Flag that controls information and progress printing.

            Default: False

        Returns
        -------
        generator: function
            The generator function of the Gaussian pyramid.
        """
        if pyramid_on_features:
            # compute features at highest level
            feature_images = []
            for c, i in enumerate(images):
                if verbose:
                    print_dynamic('- Computing feature space: {}'.format(
                        progress_bar_str((c + 1.) / len(images),
                                         show_bar=False)))
                feature_images.append(compute_features(i, feature_type[0]))
            if verbose:
                print_dynamic('- Computing feature space: Done\n')

            # create pyramid on feature_images
            generator = [i.gaussian_pyramid(n_levels=n_levels,
                                            downscale=downscale)
                         for i in feature_images]
        else:
            # create pyramid on intensities images
            # features will be computed per level
            generator = [i.gaussian_pyramid(n_levels=n_levels,
                                            downscale=downscale)
                         for i in images]
        return generator

 def post_process(self, ferns):
     basis = None
     if self.compress:
         print("\nPerforming fern compression.\n")
         # Create a new basis by randomly sampling from all fern outputs.
         basis = self._random_basis(ferns, self.basis_size)
         for i, fern in enumerate(ferns):
             print_dynamic("Compressing fern {}/{}.".format(i, len(ferns)))
             fern.compress(basis, self.compression_maxnonzero)
     return basis

 def _scale_images(cls, images, s, level_str, verbose):
     scaled_images = []
     for c, i in enumerate(images):
         if verbose:
             print_dynamic(
                 '{}Scaling features: {}'.format(
                     level_str, progress_bar_str((c + 1.) / len(images),
                                                 show_bar=False)))
         scaled_images.append(i.rescale(s))
     return scaled_images

 def _compute_reference_shape(self, images, group, label, verbose):
     # the reference_shape is the mean shape of the images' landmarks
     if verbose:
         print_dynamic('- Computing reference shape')
     shapes = [i.landmarks[group][label] for i in images]
     ref_shape = mean_pointcloud(shapes)
     # fix the reference_shape's diagonal length if specified
     if self.diagonal:
         x, y = ref_shape.range()
         scale = self.diagonal / np.sqrt(x**2 + y**2)
         ref_shape = Scale(scale, ref_shape.n_dims).apply(ref_shape)
     return ref_shape

    def _train(self, images, gt_shapes, current_shapes, increment=False,
               prefix='', verbose=False):

        if not increment:
            # Reset the regressors
            self.regressors = []

        n_perturbations = len(current_shapes[0])
        template_shape = gt_shapes[0]

        # obtain delta_x and gt_x
        delta_x, gt_x = obtain_delta_x(gt_shapes, current_shapes)

        # Cascaded Regression loop
        for k in range(self.n_iterations):
            # generate regression data
            features = features_per_image(
                images, current_shapes, self.patch_shape, self.patch_features,
                prefix='{}(Iteration {}) - '.format(prefix, k),
                verbose=verbose)

            if verbose:
                print_dynamic('{}(Iteration {}) - Performing regression'.format(
                    prefix, k))

            if not increment:
                r = self._regressor_cls()
                r.train(features, delta_x)
                self.regressors.append(r)
            else:
                self.regressors[k].increment(features, delta_x)

            # Estimate delta_points
            estimated_delta_x = self.regressors[k].predict(features)
            if verbose:
                self._print_regression_info(template_shape, gt_shapes,
                                            n_perturbations, delta_x,
                                            estimated_delta_x, k,
                                            prefix=prefix)

            j = 0
            for shapes in current_shapes:
                for s in shapes:
                    # update current x
                    current_x = s.as_vector() + estimated_delta_x[j]
                    # update current shape inplace
                    s.from_vector_inplace(current_x)
                    # update delta_x
                    delta_x[j] = gt_x[j] - current_x
                    # increase index
                    j += 1

        return current_shapes

def increment_shape_model(shape_model, shapes, forgetting_factor=None,
                          max_components=None, prefix='', verbose=False):
    r"""
    """
    if verbose:
        print_dynamic('{}Incrementing shape model'.format(prefix))
    # compute aligned shapes
    aligned_shapes = align_shapes(shapes)
    # increment shape model
    shape_model.increment(aligned_shapes, forgetting_factor=forgetting_factor)
    if max_components is not None:
        shape_model.trim_components(max_components)
    return shape_model

 def _normalize_images(self, images, group, label, ref_shape, verbose):
     # normalize the scaling of all images wrt the reference_shape size
     norm_images = []
     for c, i in enumerate(images):
         if verbose:
             print_dynamic('- Normalizing images size: {}'.format(
                 progress_bar_str((c + 1.) / len(images), show_bar=False)))
         i = rescale_to_reference_shape(i, ref_shape, group=group,
                                          label=label)
         if self.sigma:
             i.pixels = fsmooth(i.pixels, self.sigma)
         norm_images.append(i)
     return norm_images

    def _compute_features(self, images, level_str, verbose):
        feature_images = []
        for c, i in enumerate(images):
            if verbose:
                print_dynamic(
                    '{}Computing feature space: {}'.format(
                        level_str, progress_bar_str((c + 1.) / len(images),
                                                    show_bar=False)))
            if self.features:
                i = self.features(i)
            feature_images.append(i)

        return feature_images

def _compute_minimum_spanning_tree(shapes, root_vertex, level_str, verbose):
    # initialize edges and weights matrix
    n_vertices = shapes[0].n_points
    n_edges = nchoosek(n_vertices, 2)
    weights = np.zeros((n_vertices, n_vertices))
    edges = np.empty((n_edges, 2), dtype=np.int32)

    # fill edges and weights
    e = -1
    for i in range(n_vertices-1):
        for j in range(i+1, n_vertices, 1):
            # edge counter
            e += 1

            # print progress
            if verbose:
                print_dynamic('{}Computing complete graph`s weights - {}'.format(
                    level_str,
                    progress_bar_str(float(e + 1) / n_edges,
                                     show_bar=False)))

            # fill in edges
            edges[e, 0] = i
            edges[e, 1] = j

            # create data matrix of edge
            diffs_x = [s.points[i, 0] - s.points[j, 0] for s in shapes]
            diffs_y = [s.points[i, 1] - s.points[j, 1] for s in shapes]
            coords = np.array([diffs_x, diffs_y])

            # compute mean
            m = np.mean(coords, axis=1)

            # compute covariance
            c = np.cov(coords)

            # get weight
            for im in range(len(shapes)):
                weights[i, j] += -np.log(multivariate_normal.pdf(coords[:, im],
                                                                 mean=m, cov=c))
            weights[j, i] = weights[i, j]

    # create undirected graph
    complete_graph = UndirectedGraph(edges)

    if verbose:
        print_dynamic('{}Minimum spanning graph computed.\n'.format(level_str))

    # compute minimum spanning graph
    return complete_graph.minimum_spanning_tree(weights, root_vertex)