/* -*- mode: C -*-  */
/*
   IGraph library.
   Copyright (C) 2013  Gabor Csardi <csardi.gabor@gmail.com>
   334 Harvard street, Cambridge, MA 02139 USA

   This program is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 2 of the License, or
   (at your option) any later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software
   Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
   02110-1301 USA

*/

#include "igraph_embedding.h"
#include "igraph_interface.h"
#include "igraph_adjlist.h"
#include "igraph_random.h"
#include "igraph_centrality.h"
#include "igraph_blas.h"

typedef struct {
    const igraph_t *graph;
    const igraph_vector_t *cvec;
    const igraph_vector_t *cvec2;
    igraph_adjlist_t *outlist, *inlist;
    igraph_inclist_t *eoutlist, *einlist;
    igraph_vector_t *tmp;
    const igraph_vector_t *weights;
} igraph_i_asembedding_data_t;

/* Adjacency matrix, unweighted, undirected.
   Eigendecomposition is used */
int igraph_i_asembeddingu(igraph_real_t *to, const igraph_real_t *from,
                          int n, void *extra) {
    igraph_i_asembedding_data_t *data = extra;
    igraph_adjlist_t *outlist = data->outlist;
    const igraph_vector_t *cvec = data->cvec;
    igraph_vector_int_t *neis;
    int i, j, nlen;

    /* to = (A+cD) from */
    for (i = 0; i < n; i++) {
        neis = igraph_adjlist_get(outlist, i);
        nlen = igraph_vector_int_size(neis);
        to[i] = 0.0;
        for (j = 0; j < nlen; j++) {
            long int nei = (long int) VECTOR(*neis)[j];
            to[i] += from[nei];
        }
        to[i] += VECTOR(*cvec)[i] * from[i];
    }

    return 0;
}

/* Adjacency matrix, weighted, undirected.
   Eigendecomposition is used. */
int igraph_i_asembeddinguw(igraph_real_t *to, const igraph_real_t *from,
                           int n, void *extra) {
    igraph_i_asembedding_data_t *data = extra;
    igraph_inclist_t *outlist = data->eoutlist;
    const igraph_vector_t *cvec = data->cvec;
    const igraph_vector_t *weights = data->weights;
    const igraph_t *graph = data->graph;
    igraph_vector_int_t *incs;
    int i, j, nlen;

    /* to = (A+cD) from */
    for (i = 0; i < n; i++) {
        incs = igraph_inclist_get(outlist, i);
        nlen = igraph_vector_int_size(incs);
        to[i] = 0.0;
        for (j = 0; j < nlen; j++) {
            long int edge = VECTOR(*incs)[j];
            long int nei = IGRAPH_OTHER(graph, edge, i);
            igraph_real_t w = VECTOR(*weights)[edge];
            to[i] += w * from[nei];
        }
        to[i] += VECTOR(*cvec)[i] * from[i];
    }

    return 0;
}

/* Adjacency matrix, unweighted, directed. SVD. */
int igraph_i_asembedding(igraph_real_t *to, const igraph_real_t *from,
                         int n, void *extra) {
    igraph_i_asembedding_data_t *data = extra;
    igraph_adjlist_t *outlist = data->outlist;
    igraph_adjlist_t *inlist = data->inlist;
    const igraph_vector_t *cvec = data->cvec;
    igraph_vector_t *tmp = data->tmp;
    igraph_vector_int_t *neis;
    int i, j, nlen;

    /* tmp = (A+cD)' from */
    for (i = 0; i < n; i++) {
        neis = igraph_adjlist_get(inlist, i);
        nlen = igraph_vector_int_size(neis);
        VECTOR(*tmp)[i] = 0.0;
        for (j = 0; j < nlen; j++) {
            long int nei = (long int) VECTOR(*neis)[j];
            VECTOR(*tmp)[i] += from[nei];
        }
        VECTOR(*tmp)[i] += VECTOR(*cvec)[i] * from[i];
    }

    /* to = (A+cD) tmp */
    for (i = 0; i < n; i++) {
        neis = igraph_adjlist_get(outlist, i);
        nlen = igraph_vector_int_size(neis);
        to[i] = 0.0;
        for (j = 0; j < nlen; j++) {
            long int nei = (long int) VECTOR(*neis)[j];
            to[i] += VECTOR(*tmp)[nei];
        }
        to[i] += VECTOR(*cvec)[i] * VECTOR(*tmp)[i];
    }

    return 0;
}

/* Adjacency matrix, unweighted, directed. SVD, right eigenvectors */
int igraph_i_asembedding_right(igraph_real_t *to, const igraph_real_t *from,
                               int n, void *extra) {
    igraph_i_asembedding_data_t *data = extra;
    igraph_adjlist_t *inlist = data->inlist;
    const igraph_vector_t *cvec = data->cvec;
    igraph_vector_int_t *neis;
    int i, j, nlen;

    /* to = (A+cD)' from */
    for (i = 0; i < n; i++) {
        neis = igraph_adjlist_get(inlist, i);
        nlen = igraph_vector_int_size(neis);
        to[i] = 0.0;
        for (j = 0; j < nlen; j++) {
            long int nei = (long int) VECTOR(*neis)[j];
            to[i] += from[nei];
        }
        to[i] += VECTOR(*cvec)[i] * from[i];
    }

    return 0;
}

/* Adjacency matrix, weighted, directed. SVD. */
int igraph_i_asembeddingw(igraph_real_t *to, const igraph_real_t *from,
                          int n, void *extra) {
    igraph_i_asembedding_data_t *data = extra;
    igraph_inclist_t *outlist = data->eoutlist;
    igraph_inclist_t *inlist = data->einlist;
    const igraph_vector_t *cvec = data->cvec;
    const igraph_vector_t *weights = data->weights;
    const igraph_t *graph = data->graph;
    igraph_vector_t *tmp = data->tmp;
    igraph_vector_int_t *incs;
    int i, j, nlen;

    /* tmp = (A+cD)' from */
    for (i = 0; i < n; i++) {
        incs = igraph_inclist_get(inlist, i);
        nlen = igraph_vector_int_size(incs);
        VECTOR(*tmp)[i] = 0.0;
        for (j = 0; j < nlen; j++) {
            long int edge = VECTOR(*incs)[j];
            long int nei = IGRAPH_OTHER(graph, edge, i);
            igraph_real_t w = VECTOR(*weights)[edge];
            VECTOR(*tmp)[i] += w * from[nei];
        }
        VECTOR(*tmp)[i] += VECTOR(*cvec)[i] * from[i];
    }

    /* to = (A+cD) tmp */
    for (i = 0; i < n; i++) {
        incs = igraph_inclist_get(outlist, i);
        nlen = igraph_vector_int_size(incs);
        to[i] = 0.0;
        for (j = 0; j < nlen; j++) {
            long int edge = VECTOR(*incs)[j];
            long int nei = IGRAPH_OTHER(graph, edge, i);
            igraph_real_t w = VECTOR(*weights)[edge];
            to[i] += w * VECTOR(*tmp)[nei];
        }
        to[i] += VECTOR(*cvec)[i] * VECTOR(*tmp)[i];
    }

    return 0;
}

/* Adjacency matrix, weighted, directed. SVD, right eigenvectors. */
int igraph_i_asembeddingw_right(igraph_real_t *to, const igraph_real_t *from,
                                int n, void *extra) {
    igraph_i_asembedding_data_t *data = extra;
    igraph_inclist_t *inlist = data->einlist;
    const igraph_vector_t *cvec = data->cvec;
    const igraph_vector_t *weights = data->weights;
    const igraph_t *graph = data->graph;
    igraph_vector_int_t *incs;
    int i, j, nlen;

    /* to = (A+cD)' from */
    for (i = 0; i < n; i++) {
        incs = igraph_inclist_get(inlist, i);
        nlen = igraph_vector_int_size(incs);
        to[i] = 0.0;
        for (j = 0; j < nlen; j++) {
            long int edge = VECTOR(*incs)[j];
            long int nei = IGRAPH_OTHER(graph, edge, i);
            igraph_real_t w = VECTOR(*weights)[edge];
            to[i] += w * from[nei];
        }
        to[i] += VECTOR(*cvec)[i] * from[i];
    }

    return 0;
}

/* Laplacian D-A, unweighted, undirected. Eigendecomposition. */
int igraph_i_lsembedding_da(igraph_real_t *to, const igraph_real_t *from,
                            int n, void *extra) {
    igraph_i_asembedding_data_t *data = extra;
    igraph_adjlist_t *outlist = data->outlist;
    const igraph_vector_t *cvec = data->cvec;
    igraph_vector_int_t *neis;
    int i, j, nlen;

    /* to = (D-A) from */
    for (i = 0; i < n; i++) {
        neis = igraph_adjlist_get(outlist, i);
        nlen = igraph_vector_int_size(neis);
        to[i] = 0.0;
        for (j = 0; j < nlen; j++) {
            long int nei = (long int) VECTOR(*neis)[j];
            to[i] -= from[nei];
        }
        to[i] += VECTOR(*cvec)[i] * from[i];
    }

    return 0;
}

/* Laplacian D-A, weighted, undirected. Eigendecomposition. */
int igraph_i_lsembedding_daw(igraph_real_t *to, const igraph_real_t *from,
                             int n, void *extra) {
    igraph_i_asembedding_data_t *data = extra;
    igraph_inclist_t *outlist = data->eoutlist;
    const igraph_vector_t *cvec = data->cvec;
    const igraph_vector_t *weights = data->weights;
    const igraph_t *graph = data->graph;
    igraph_vector_int_t *incs;
    int i, j, nlen;

    /* to = (D-A) from */
    for (i = 0; i < n; i++) {
        incs = igraph_inclist_get(outlist, i);
        nlen = igraph_vector_int_size(incs);
        to[i] = 0.0;
        for (j = 0; j < nlen; j++) {
            long int edge = VECTOR(*incs)[j];
            long int nei = IGRAPH_OTHER(graph, edge, i);
            igraph_real_t w = VECTOR(*weights)[edge];
            to[i] -= w * from[nei];
        }
        to[i] += VECTOR(*cvec)[i] * from[i];
    }

    return 0;
}

/* Laplacian DAD, unweighted, undirected. Eigendecomposition. */
int igraph_i_lsembedding_dad(igraph_real_t *to, const igraph_real_t *from,
                             int n, void *extra) {

    igraph_i_asembedding_data_t *data = extra;
    igraph_adjlist_t *outlist = data->outlist;
    const igraph_vector_t *cvec = data->cvec;
    igraph_vector_t *tmp = data->tmp;
    igraph_vector_int_t *neis;
    int i, j, nlen;

    /* to = D^1/2 from */
    for (i = 0; i < n; i++) {
        to[i] = VECTOR(*cvec)[i] * from[i];
    }

    /* tmp = A to */
    for (i = 0; i < n; i++) {
        neis = igraph_adjlist_get(outlist, i);
        nlen = igraph_vector_int_size(neis);
        VECTOR(*tmp)[i] = 0.0;
        for (j = 0; j < nlen; j++) {
            long int nei = (long int) VECTOR(*neis)[j];
            VECTOR(*tmp)[i] += to[nei];
        }
    }

    /* to = D tmp */
    for (i = 0; i < n; i++) {
        to[i] = VECTOR(*cvec)[i] * VECTOR(*tmp)[i];
    }

    return 0;
}

int igraph_i_lsembedding_dadw(igraph_real_t *to, const igraph_real_t *from,
                              int n, void *extra) {

    igraph_i_asembedding_data_t *data = extra;
    igraph_inclist_t *outlist = data->eoutlist;
    const igraph_vector_t *cvec = data->cvec;
    const igraph_vector_t *weights = data->weights;
    const igraph_t *graph = data->graph;
    igraph_vector_t *tmp = data->tmp;
    igraph_vector_int_t *incs;
    int i, j, nlen;

    /* to = D^-1/2 from */
    for (i = 0; i < n; i++) {
        to[i] = VECTOR(*cvec)[i] * from[i];
    }

    /* tmp = A' to */
    for (i = 0; i < n; i++) {
        incs = igraph_inclist_get(outlist, i);
        nlen = igraph_vector_int_size(incs);
        VECTOR(*tmp)[i] = 0.0;
        for (j = 0; j < nlen; j++) {
            long int edge = VECTOR(*incs)[j];
            long int nei = IGRAPH_OTHER(graph, edge, i);
            igraph_real_t w = VECTOR(*weights)[edge];
            VECTOR(*tmp)[i] += w * to[nei];
        }
    }

    /* to = D tmp */
    for (i = 0; i < n; i++) {
        to[i] = VECTOR(*cvec)[i] * VECTOR(*cvec)[i] * VECTOR(*tmp)[i];
    }

    /* tmp = A to */
    for (i = 0; i < n; i++) {
        incs = igraph_inclist_get(outlist, i);
        nlen = igraph_vector_int_size(incs);
        VECTOR(*tmp)[i] = 0.0;
        for (j = 0; j < nlen; j++) {
            long int edge = VECTOR(*incs)[j];
            long int nei = IGRAPH_OTHER(graph, edge, i);
            igraph_real_t w = VECTOR(*weights)[edge];
            VECTOR(*tmp)[i] += w * to[nei];
        }
    }

    /* to = D^-1/2 tmp */
    for (i = 0; i < n; i++) {
        to[i] = VECTOR(*cvec)[i] * VECTOR(*tmp)[i];
    }

    return 0;
}

/* Laplacian I-DAD, unweighted, undirected. Eigendecomposition. */
int igraph_i_lsembedding_idad(igraph_real_t *to, const igraph_real_t *from,
                              int n, void *extra) {

    int i;

    igraph_i_lsembedding_dad(to, from, n, extra);
    for (i = 0; i < n; i++) {
        to[i] = from[i] - to[i];
    }

    return 0;
}

int igraph_i_lsembedding_idadw(igraph_real_t *to, const igraph_real_t *from,
                               int n, void *extra) {
    int i;

    igraph_i_lsembedding_dadw(to, from, n, extra);
    for (i = 0; i < n; i++) {
        to[i] = from[i] - to[i];
    }

    return 0;
}

/* Laplacian OAP, unweighted, directed. SVD. */
int igraph_i_lseembedding_oap(igraph_real_t *to, const igraph_real_t *from,
                              int n, void *extra) {

    igraph_i_asembedding_data_t *data = extra;
    igraph_adjlist_t *outlist = data->outlist;
    igraph_adjlist_t *inlist = data->inlist;
    const igraph_vector_t *deg_in = data->cvec;
    const igraph_vector_t *deg_out = data->cvec2;
    igraph_vector_t *tmp = data->tmp;
    igraph_vector_int_t *neis;
    int i, j, nlen;

    /* tmp = O' from */
    for (i = 0; i < n; i++) {
        VECTOR(*tmp)[i] = VECTOR(*deg_out)[i] * from[i];
    }

    /* to = A' tmp */
    for (i = 0; i < n; i++) {
        neis = igraph_adjlist_get(inlist, i);
        nlen = igraph_vector_int_size(neis);
        to[i] = 0.0;
        for (j = 0; j < nlen; j++) {
            int nei = VECTOR(*neis)[j];
            to[i] += VECTOR(*tmp)[nei];
        }
    }

    /* tmp = P' to */
    for (i = 0; i < n; i++) {
        VECTOR(*tmp)[i] = VECTOR(*deg_in)[i] * to[i];
    }

    /* to = P tmp */
    for (i = 0; i < n; i++) {
        to[i] = VECTOR(*deg_in)[i] * VECTOR(*tmp)[i];
    }

    /* tmp = A to */
    for (i = 0; i < n; i++) {
        neis = igraph_adjlist_get(outlist, i);
        nlen = igraph_vector_int_size(neis);
        VECTOR(*tmp)[i] = 0.0;
        for (j = 0; j < nlen; j++) {
            int nei = VECTOR(*neis)[j];
            VECTOR(*tmp)[i] += to[nei];
        }
    }

    /* to = O tmp */
    for (i = 0; i < n; i++) {
        to[i] = VECTOR(*deg_out)[i] * VECTOR(*tmp)[i];
    }

    return 0;
}

/* Laplacian OAP, unweighted, directed. SVD, right eigenvectors. */
int igraph_i_lseembedding_oap_right(igraph_real_t *to,
                                    const igraph_real_t *from,
                                    int n, void *extra) {
    igraph_i_asembedding_data_t *data = extra;
    igraph_adjlist_t *inlist = data->inlist;
    const igraph_vector_t *deg_in = data->cvec;
    const igraph_vector_t *deg_out = data->cvec2;
    igraph_vector_t *tmp = data->tmp;
    igraph_vector_int_t *neis;
    int i, j, nlen;

    /* to = O' from */
    for (i = 0; i < n; i++) {
        to[i] = VECTOR(*deg_out)[i] * from[i];
    }

    /* tmp = A' to */
    for (i = 0; i < n; i++) {
        neis = igraph_adjlist_get(inlist, i);
        nlen = igraph_vector_int_size(neis);
        VECTOR(*tmp)[i] = 0.0;
        for (j = 0; j < nlen; j++) {
            int nei = VECTOR(*neis)[j];
            VECTOR(*tmp)[i] += to[nei];
        }
    }

    /* to = P' tmp */
    for (i = 0; i < n; i++) {
        to[i] = VECTOR(*deg_in)[i] * VECTOR(*tmp)[i];
    }

    return 0;
}

/* Laplacian OAP, weighted, directed. SVD. */
int igraph_i_lseembedding_oapw(igraph_real_t *to, const igraph_real_t *from,
                               int n, void *extra) {

    igraph_i_asembedding_data_t *data = extra;
    igraph_inclist_t *outlist = data->eoutlist;
    igraph_inclist_t *inlist = data->einlist;
    const igraph_vector_t *deg_in = data->cvec;
    const igraph_vector_t *deg_out = data->cvec2;
    const igraph_vector_t *weights = data->weights;
    const igraph_t *graph = data->graph;
    igraph_vector_t *tmp = data->tmp;
    igraph_vector_int_t *neis;
    int i, j, nlen;

    /* tmp = O' from */
    for (i = 0; i < n; i++) {
        VECTOR(*tmp)[i] = VECTOR(*deg_out)[i] * from[i];
    }

    /* to = A' tmp */
    for (i = 0; i < n; i++) {
        neis = igraph_inclist_get(inlist, i);
        nlen = igraph_vector_int_size(neis);
        to[i] = 0.0;
        for (j = 0; j < nlen; j++) {
            int edge = VECTOR(*neis)[j];
            int nei = IGRAPH_OTHER(graph, edge, i);
            igraph_real_t w = VECTOR(*weights)[edge];
            to[i] += w * VECTOR(*tmp)[nei];
        }
    }

    /* tmp = P' to */
    for (i = 0; i < n; i++) {
        VECTOR(*tmp)[i] = VECTOR(*deg_in)[i] * to[i];
    }

    /* to = P tmp */
    for (i = 0; i < n; i++) {
        to[i] = VECTOR(*deg_in)[i] * VECTOR(*tmp)[i];
    }

    /* tmp = A to */
    for (i = 0; i < n; i++) {
        neis = igraph_inclist_get(outlist, i);
        nlen = igraph_vector_int_size(neis);
        VECTOR(*tmp)[i] = 0.0;
        for (j = 0; j < nlen; j++) {
            int edge = VECTOR(*neis)[j];
            int nei = IGRAPH_OTHER(graph, edge, i);
            igraph_real_t w = VECTOR(*weights)[edge];
            VECTOR(*tmp)[i] += w * to[nei];
        }
    }

    /* to = O tmp */
    for (i = 0; i < n; i++) {
        to[i] = VECTOR(*deg_out)[i] * VECTOR(*tmp)[i];
    }

    return 0;
}

/* Laplacian OAP, weighted, directed. SVD, right eigenvectors. */
int igraph_i_lseembedding_oapw_right(igraph_real_t *to,
                                     const igraph_real_t *from,
                                     int n, void *extra) {
    igraph_i_asembedding_data_t *data = extra;
    igraph_inclist_t *inlist = data->einlist;
    const igraph_vector_t *deg_in = data->cvec;
    const igraph_vector_t *deg_out = data->cvec2;
    const igraph_vector_t *weights = data->weights;
    const igraph_t *graph = data->graph;
    igraph_vector_t *tmp = data->tmp;
    igraph_vector_int_t *neis;
    int i, j, nlen;

    /* to = O' from */
    for (i = 0; i < n; i++) {
        to[i] = VECTOR(*deg_out)[i] * from[i];
    }

    /* tmp = A' to */
    for (i = 0; i < n; i++) {
        neis = igraph_inclist_get(inlist, i);
        nlen = igraph_vector_int_size(neis);
        VECTOR(*tmp)[i] = 0.0;
        for (j = 0; j < nlen; j++) {
            int edge = VECTOR(*neis)[j];
            int nei = IGRAPH_OTHER(graph, edge, i);
            igraph_real_t w = VECTOR(*weights)[edge];
            VECTOR(*tmp)[i] += w * to[nei];
        }
    }

    /* to = P' tmp */
    for (i = 0; i < n; i++) {
        to[i] = VECTOR(*deg_in)[i] * VECTOR(*tmp)[i];
    }

    return 0;
}

int igraph_i_spectral_embedding(const igraph_t *graph,
                                igraph_integer_t no,
                                const igraph_vector_t *weights,
                                igraph_eigen_which_position_t which,
                                igraph_bool_t scaled,
                                igraph_matrix_t *X,
                                igraph_matrix_t *Y,
                                igraph_vector_t *D,
                                const igraph_vector_t *cvec,
                                const igraph_vector_t *cvec2,
                                igraph_arpack_options_t *options,
                                igraph_arpack_function_t *callback,
                                igraph_arpack_function_t *callback_right,
                                igraph_bool_t symmetric,
                                igraph_bool_t eigen,
                                igraph_bool_t zapsmall) {

    igraph_integer_t vc = igraph_vcount(graph);
    igraph_vector_t tmp;
    igraph_adjlist_t outlist, inlist;
    igraph_inclist_t eoutlist, einlist;
    int i, j, cveclen = igraph_vector_size(cvec);
    igraph_i_asembedding_data_t data = { graph, cvec, cvec2, &outlist, &inlist,
                                         &eoutlist, &einlist, &tmp, weights
                                       };
    igraph_vector_t tmpD;

    if (weights && igraph_vector_size(weights) != igraph_ecount(graph)) {
        IGRAPH_ERROR("Invalid weight vector length", IGRAPH_EINVAL);
    }

    if (which != IGRAPH_EIGEN_LM &&
        which != IGRAPH_EIGEN_LA &&
        which != IGRAPH_EIGEN_SA) {
        IGRAPH_ERROR("Invalid eigenvalue chosen, must be one of "
                     "`largest magnitude', `largest algebraic' or "
                     "`smallest algebraic'", IGRAPH_EINVAL);
    }

    if (no > vc) {
        IGRAPH_ERROR("Too many singular values requested", IGRAPH_EINVAL);
    }
    if (no <= 0) {
        IGRAPH_ERROR("No singular values requested", IGRAPH_EINVAL);
    }

    if (cveclen != 1 && cveclen != vc) {
        IGRAPH_ERROR("Augmentation vector size is invalid, it should be "
                     "the number of vertices or scalar", IGRAPH_EINVAL);
    }

    IGRAPH_CHECK(igraph_matrix_resize(X, vc, no));
    if (Y) {
        IGRAPH_CHECK(igraph_matrix_resize(Y, vc, no));
    }

    /* empty graph */
    if (igraph_ecount(graph) == 0) {
        igraph_matrix_null(X);
        if (Y) {
            igraph_matrix_null(Y);
        }
        return 0;
    }

    igraph_vector_init(&tmp, vc);
    IGRAPH_FINALLY(igraph_vector_destroy, &tmp);
    if (!weights) {
        IGRAPH_CHECK(igraph_adjlist_init(graph, &outlist, IGRAPH_OUT));
        IGRAPH_FINALLY(igraph_adjlist_destroy, &outlist);
        if (!symmetric) {
            IGRAPH_CHECK(igraph_adjlist_init(graph, &inlist, IGRAPH_IN));
            IGRAPH_FINALLY(igraph_adjlist_destroy, &inlist);
        }
    } else {
        IGRAPH_CHECK(igraph_inclist_init(graph, &eoutlist, IGRAPH_OUT));
        IGRAPH_FINALLY(igraph_inclist_destroy, &eoutlist);
        if (!symmetric) {
            IGRAPH_CHECK(igraph_inclist_init(graph, &einlist, IGRAPH_IN));
            IGRAPH_FINALLY(igraph_inclist_destroy, &einlist);
        }
    }
    IGRAPH_VECTOR_INIT_FINALLY(&tmpD, no);

    options->n = vc;
    options->start = 0;   /* random start vector */
    options->nev = no;
    switch (which) {
    case IGRAPH_EIGEN_LM:
        options->which[0] = 'L'; options->which[1] = 'M';
        break;
    case IGRAPH_EIGEN_LA:
        options->which[0] = 'L'; options->which[1] = 'A';
        break;
    case IGRAPH_EIGEN_SA:
        options->which[0] = 'S'; options->which[1] = 'A';
        break;
    default:
        break;
    }
    options->ncv = no + 3;
    if (options->ncv > vc) {
        options->ncv = vc;
    }

    IGRAPH_CHECK(igraph_arpack_rssolve(callback, &data, options, 0, &tmpD, X));

    if (!symmetric) {
        /* calculate left eigenvalues */
        IGRAPH_CHECK(igraph_matrix_resize(Y, vc, no));
        for (i = 0; i < no; i++) {
            igraph_real_t norm;
            igraph_vector_t v;
            callback_right(&MATRIX(*Y, 0, i), &MATRIX(*X, 0, i), vc, &data);
            igraph_vector_view(&v, &MATRIX(*Y, 0, i), vc);
            norm = 1.0 / igraph_blas_dnrm2(&v);
            igraph_vector_scale(&v, norm);
        }
    } else if (Y) {
        IGRAPH_CHECK(igraph_matrix_update(Y, X));
    }

    if (zapsmall) {
        igraph_vector_zapsmall(&tmpD, 0);
        igraph_matrix_zapsmall(X, 0);
        if (Y) {
            igraph_matrix_zapsmall(Y, 0);
        }
    }

    if (D) {
        igraph_vector_update(D, &tmpD);
        if (!eigen) {
            for (i = 0; i < no; i++) {
                VECTOR(*D)[i] = sqrt(VECTOR(*D)[i]);
            }
        }
    }

    if (scaled) {
        if (eigen) {
            /* eigenvalues were calculated */
            for (i = 0; i < no; i++) {
                VECTOR(tmpD)[i] = sqrt(fabs(VECTOR(tmpD)[i]));
            }
        } else {
            /* singular values were calculated */
            for (i = 0; i < no; i++) {
                VECTOR(tmpD)[i] = sqrt(sqrt(VECTOR(tmpD)[i]));
            }
        }

        for (j = 0; j < vc; j++) {
            for (i = 0; i < no; i++) {
                MATRIX(*X, j, i) *= VECTOR(tmpD)[i];
            }
        }

        if (Y) {
            for (j = 0; j < vc; j++) {
                for (i = 0; i < no; i++) {
                    MATRIX(*Y, j, i) *= VECTOR(tmpD)[i];
                }
            }
        }
    }

    igraph_vector_destroy(&tmpD);
    if (!weights) {
        if (!symmetric) {
            igraph_adjlist_destroy(&inlist);
            IGRAPH_FINALLY_CLEAN(1);
        }
        igraph_adjlist_destroy(&outlist);
    } else {
        if (!symmetric) {
            igraph_inclist_destroy(&einlist);
            IGRAPH_FINALLY_CLEAN(1);
        }
        igraph_inclist_destroy(&eoutlist);
    }
    igraph_vector_destroy(&tmp);
    IGRAPH_FINALLY_CLEAN(3);

    return 0;
}

/**
 * \function igraph_adjacency_spectral_embedding
 * Adjacency spectral embedding
 *
 * Spectral decomposition of the adjacency matrices of graphs.
 * This function computes a \code{no}-dimensional Euclidean
 * representation of the graph based on its adjacency
 * matrix, A. This representation is computed via the singular value
 * decomposition of the adjacency matrix, A=UDV^T. In the case,
 * where the graph is a random dot product graph generated using latent
 * position vectors in R^no for each vertex, the embedding will
 * provide an estimate of these latent vectors.
 *
 * </para><para>
 * For undirected graphs the latent positions are calculated as
 * X=U^no D^(1/2) where U^no equals to the first no columns of U, and
 * D^(1/2) is a diagonal matrix containing the square root of the selected
 * singular values on the diagonal.
 *
 * </para><para>
 * For directed graphs the embedding is defined as the pair
 * X=U^no D^(1/2), Y=V^no D^(1/2). (For undirected graphs U=V,
 * so it is enough to keep one of them.)
 *
 * \param graph The input graph, can be directed or undirected.
 * \param no An integer scalar. This value is the embedding dimension of
 *        the spectral embedding. Should be smaller than the number of
 *        vertices. The largest no-dimensional non-zero
 *        singular values are used for the spectral embedding.
 * \param weights Optional edge weights. Supply a null pointer for
 *        unweighted graphs.
 * \param which Which eigenvalues (or singular values, for directed
 *        graphs) to use, possible values:
 *        \clist
 *          \cli IGRAPH_EIGEN_LM
 *          the ones with the largest magnitude
 *          \cli IGRAPH_EIGEN_LA
 *          the (algebraic) largest ones
 *          \cli IGRAPH_EIGEN_SA
 *          the (algebraic) smallest ones.
 *        \endclist
 *        For directed graphs, <code>IGRAPH_EIGEN_LM</code> and
 *        <code>IGRAPH_EIGEN_LA</code> are the same because singular
 *        values are used for the ordering instead of eigenvalues.
 * \param scaled Whether to return X and Y (if scaled is non-zero), or
 *        U and V.
 * \param X Initialized matrix, the estimated latent positions are
 *        stored here.
 * \param Y Initialized matrix or a null pointer. If not a null
 *        pointer, then the second half of the latent positions are
 *        stored here. (For undirected graphs, this always equals X.)
 * \param D Initialized vector or a null pointer. If not a null
 *        pointer, then the eigenvalues (for undirected graphs) or the
 *        singular values (for directed graphs) are stored here.
 * \param cvec A numeric vector, its length is the number vertices in the
 *        graph. This vector is added to the diagonal of the adjacency
 *        matrix, before performing the SVD.
 * \param options Options to ARPACK. See \ref igraph_arpack_options_t
 *        for details. Note that the function overwrites the
 *        <code>n</code> (number of vertices), <code>nev</code> and
 *        <code>which</code> parameters and it always starts the
 *        calculation from a random start vector.
 * \return Error code.
 *
 */

int igraph_adjacency_spectral_embedding(const igraph_t *graph,
                                        igraph_integer_t no,
                                        const igraph_vector_t *weights,
                                        igraph_eigen_which_position_t which,
                                        igraph_bool_t scaled,
                                        igraph_matrix_t *X,
                                        igraph_matrix_t *Y,
                                        igraph_vector_t *D,
                                        const igraph_vector_t *cvec,
                                        igraph_arpack_options_t *options) {

    igraph_arpack_function_t *callback, *callback_right;
    igraph_bool_t directed = igraph_is_directed(graph);

    if (directed) {
        callback = weights ? igraph_i_asembeddingw : igraph_i_asembedding;
        callback_right = (weights ? igraph_i_asembeddingw_right :
                          igraph_i_asembedding_right);
    } else {
        callback = weights ? igraph_i_asembeddinguw : igraph_i_asembeddingu;
        callback_right = 0;
    }

    return igraph_i_spectral_embedding(graph, no, weights, which, scaled,
                                       X, Y, D, cvec, /* deg2=*/ 0,
                                       options, callback, callback_right,
                                       /*symmetric=*/ !directed,
                                       /*eigen=*/ !directed, /*zapsmall=*/ 1);
}

int igraph_i_lse_und(const igraph_t *graph,
                     igraph_integer_t no,
                     const igraph_vector_t *weights,
                     igraph_eigen_which_position_t which,
                     igraph_neimode_t degmode,
                     igraph_laplacian_spectral_embedding_type_t type,
                     igraph_bool_t scaled,
                     igraph_matrix_t *X,
                     igraph_matrix_t *Y,
                     igraph_vector_t *D,
                     igraph_arpack_options_t *options) {

    igraph_arpack_function_t *callback;
    igraph_vector_t deg;

    switch (type) {
    case IGRAPH_EMBEDDING_D_A:
        callback = weights ? igraph_i_lsembedding_daw : igraph_i_lsembedding_da;
        break;
    case IGRAPH_EMBEDDING_DAD:
        callback = weights ? igraph_i_lsembedding_dadw : igraph_i_lsembedding_dad;
        break;
    case IGRAPH_EMBEDDING_I_DAD:
        callback = weights ? igraph_i_lsembedding_idadw : igraph_i_lsembedding_idad;
        break;
    default:
        IGRAPH_ERROR("Invalid Laplacian spectral embedding type",
                     IGRAPH_EINVAL);
        break;
    }

    IGRAPH_VECTOR_INIT_FINALLY(&deg, 0);
    igraph_strength(graph, &deg, igraph_vss_all(), IGRAPH_ALL, /*loops=*/ 1,
                    weights);

    switch (type) {
    case IGRAPH_EMBEDDING_D_A:
        break;
    case IGRAPH_EMBEDDING_DAD:
    case IGRAPH_EMBEDDING_I_DAD: {
        int i, n = igraph_vector_size(&deg);
        for (i = 0; i < n; i++) {
            VECTOR(deg)[i] = 1.0 / sqrt(VECTOR(deg)[i]);
        }
    }
    break;
    default:
        break;
    }

    IGRAPH_CHECK(igraph_i_spectral_embedding(graph, no, weights, which,
                 scaled, X, Y, D, /*cvec=*/ &deg, /*deg2=*/ 0,
                 options, callback, 0, /*symmetric=*/ 1,
                 /*eigen=*/ 1, /*zapsmall=*/ 1));

    igraph_vector_destroy(&deg);
    IGRAPH_FINALLY_CLEAN(1);

    return 0;
}

int igraph_i_lse_dir(const igraph_t *graph,
                     igraph_integer_t no,
                     const igraph_vector_t *weights,
                     igraph_eigen_which_position_t which,
                     igraph_neimode_t degmode,
                     igraph_laplacian_spectral_embedding_type_t type,
                     igraph_bool_t scaled,
                     igraph_matrix_t *X,
                     igraph_matrix_t *Y,
                     igraph_vector_t *D,
                     igraph_arpack_options_t *options) {

    igraph_arpack_function_t *callback =
        weights ? igraph_i_lseembedding_oapw : igraph_i_lseembedding_oap;
    igraph_arpack_function_t *callback_right =
        weights ? igraph_i_lseembedding_oapw_right :
        igraph_i_lseembedding_oap_right;
    igraph_vector_t deg_in, deg_out;
    int i, n = igraph_vcount(graph);

    if (type != IGRAPH_EMBEDDING_OAP) {
        IGRAPH_ERROR("Invalid Laplacian spectral embedding type", IGRAPH_EINVAL);
    }

    IGRAPH_VECTOR_INIT_FINALLY(&deg_in, n);
    IGRAPH_VECTOR_INIT_FINALLY(&deg_out, n);
    igraph_strength(graph, &deg_in, igraph_vss_all(), IGRAPH_IN, /*loops=*/ 1,
                    weights);
    igraph_strength(graph, &deg_out, igraph_vss_all(), IGRAPH_OUT, /*loops=*/ 1,
                    weights);

    for (i = 0; i < n; i++) {
        VECTOR(deg_in)[i] = 1.0 / sqrt(VECTOR(deg_in)[i]);
        VECTOR(deg_out)[i] = 1.0 / sqrt(VECTOR(deg_out)[i]);
    }

    IGRAPH_CHECK(igraph_i_spectral_embedding(graph, no, weights, which,
                 scaled, X, Y, D, /*cvec=*/ &deg_in,
                 /*deg2=*/ &deg_out, options, callback,
                 callback_right, /*symmetric=*/ 0, /*eigen=*/ 0,
                 /*zapsmall=*/ 1));

    igraph_vector_destroy(&deg_in);
    igraph_vector_destroy(&deg_out);
    IGRAPH_FINALLY_CLEAN(2);

    return 0;
}

/**
 * \function igraph_laplacian_spectral_embedding
 * Spectral embedding of the Laplacian of a graph
 *
 * This function essentially does the same as
 * \ref igraph_adjacency_spectral_embedding, but works on the Laplacian
 * of the graph, instead of the adjacency matrix.
 * \param graph The input graph.
 * \param no The number of eigenvectors (or singular vectors if the graph
 *        is directed) to use for the embedding.
 * \param weights Optional edge weights. Supply a null pointer for
 *        unweighted graphs.
 * \param which Which eigenvalues (or singular values, for directed
 *        graphs) to use, possible values:
 *        \clist
 *          \cli IGRAPH_EIGEN_LM
 *          the ones with the largest magnitude
 *          \cli IGRAPH_EIGEN_LA
 *          the (algebraic) largest ones
 *          \cli IGRAPH_EIGEN_SA
 *          the (algebraic) smallest ones.
 *        \endclist
 *        For directed graphs, <code>IGRAPH_EIGEN_LM</code> and
 *        <code>IGRAPH_EIGEN_LA</code> are the same because singular
 *        values are used for the ordering instead of eigenvalues.
 * \param type The type of the Laplacian to use. Various definitions
 *        exist for the Laplacian of a graph, and one can choose
 *        between them with this argument. Possible values:
 *        \clist
 *          \cli IGRAPH_EMBEDDING_D_A
 *           means D - A where D is the
 *           degree matrix and A is the adjacency matrix
 *          \cli IGRAPH_EMBEDDING_DAD
 *           means Di times A times Di,
 *           where Di is the inverse of the square root of the degree matrix;
 *          \cli IGRAPH_EMBEDDING_I_DAD
 *          means I - Di A Di, where I
 *          is the identity matrix.
 *        \endclist
 * \param scaled Whether to return X and Y (if scaled is non-zero), or
 *        U and V.
 * \param X Initialized matrix, the estimated latent positions are
 *        stored here.
 * \param Y Initialized matrix or a null pointer. If not a null
 *        pointer, then the second half of the latent positions are
 *        stored here. (For undirected graphs, this always equals X.)
 * \param D Initialized vector or a null pointer. If not a null
 *        pointer, then the eigenvalues (for undirected graphs) or the
 *        singular values (for directed graphs) are stored here.
 * \param options Options to ARPACK. See \ref igraph_arpack_options_t
 *        for details. Note that the function overwrites the
 *        <code>n</code> (number of vertices), <code>nev</code> and
 *        <code>which</code> parameters and it always starts the
 *        calculation from a random start vector.
 * \return Error code.
 *
 * \sa \ref igraph_adjacency_spectral_embedding to embed the adjacency
 * matrix.
 */

int igraph_laplacian_spectral_embedding(const igraph_t *graph,
                                        igraph_integer_t no,
                                        const igraph_vector_t *weights,
                                        igraph_eigen_which_position_t which,
                                        igraph_neimode_t degmode,
                                        igraph_laplacian_spectral_embedding_type_t type,
                                        igraph_bool_t scaled,
                                        igraph_matrix_t *X,
                                        igraph_matrix_t *Y,
                                        igraph_vector_t *D,
                                        igraph_arpack_options_t *options) {

    if (igraph_is_directed(graph)) {
        return igraph_i_lse_dir(graph, no, weights, which, degmode, type, scaled,
                                X, Y, D, options);
    } else {
        return igraph_i_lse_und(graph, no, weights, which, degmode, type, scaled,
                                X, Y, D, options);
    }
}

/**
 * \function igraph_dim_select
 * Dimensionality selection
 *
 * Dimensionality selection for singular values using
 * profile likelihood.
 *
 * </para><para>
 * The input of the function is a numeric vector which contains
 * the measure of "importance" for each dimension.
 *
 * </para><para>
 * For spectral embedding, these are the singular values of the adjacency
 * matrix. The singular values are assumed to be generated from a
 * Gaussian mixture distribution with two components that have different
 * means and same variance. The dimensionality d is chosen to
 * maximize the likelihood when the d largest singular values are
 * assigned to one component of the mixture and the rest of the singular
 * values assigned to the other component.
 *
 * </para><para>
 * This function can also be used for the general separation problem,
 * where we assume that the left and the right of the vector are coming
 * from two Normal distributions, with different means, and we want
 * to know their border.
 *
 * \param sv A numeric vector, the ordered singular values.
 * \param dim The result is stored here.
 * \return Error code.
 *
 * Time complexity: O(n), n is the number of values in sv.
 *
 * \sa \ref igraph_adjacency_spectral_embedding().
 */

int igraph_dim_select(const igraph_vector_t *sv, igraph_integer_t *dim) {

    int i, n = igraph_vector_size(sv);
    igraph_real_t x, x2, sum1 = 0.0, sum2 = igraph_vector_sum(sv);
    igraph_real_t sumsq1 = 0.0, sumsq2 = 0.0; /* to be set */
    igraph_real_t oldmean1, oldmean2, mean1 = 0.0, mean2 = sum2 / n;
    igraph_real_t varsq1 = 0.0, varsq2 = 0.0; /* to be set */
    igraph_real_t var1, var2, sd, profile, max = IGRAPH_NEGINFINITY;

    if (n == 0) {
        IGRAPH_ERROR("Need at least one singular value for dimensionality "
                     "selection", IGRAPH_EINVAL);
    }

    if (n == 1) {
        *dim = 1;
        return 0;
    }

    for (i = 0; i < n; i++) {
        x = VECTOR(*sv)[i];
        sumsq2 += x * x;
        varsq2 += (mean2 - x) * (mean2 - x);
    }

    for (i = 0; i < n - 1; i++) {
        int n1 = i + 1, n2 = n - i - 1, n1m1 = n1 - 1, n2m1 = n2 - 1;
        x = VECTOR(*sv)[i]; x2 = x * x;
        sum1 += x; sum2 -= x;
        sumsq1 += x2; sumsq2 -= x2;
        oldmean1 = mean1; oldmean2 = mean2;
        mean1 = sum1 / n1; mean2 = sum2 / n2;
        varsq1 += (x - oldmean1) * (x - mean1);
        varsq2 -= (x - oldmean2) * (x - mean2);
        var1 = i == 0 ? 0 : varsq1 / n1m1;
        var2 = i == n - 2 ? 0 : varsq2 / n2m1;
        sd = sqrt(( n1m1 * var1 + n2m1 * var2) / (n - 2));
        profile = /* - n * log(2.0*M_PI)/2.0 */ /* This is redundant */
            - n * log(sd) -
            ((sumsq1 - 2 * mean1 * sum1 + n1 * mean1 * mean1) +
             (sumsq2 - 2 * mean2 * sum2 + n2 * mean2 * mean2)) / 2.0 / sd / sd;
        if (profile > max) {
            max = profile;
            *dim = n1;
        }
    }

    /* Plus the last case, all elements in one group */
    x = VECTOR(*sv)[n - 1];
    sum1 += x;
    oldmean1 = mean1;
    mean1 = sum1 / n;
    sumsq1 += x * x;
    varsq1 += (x - oldmean1) * (x - mean1);
    var1 = varsq1 / (n - 1);
    sd = sqrt(var1);
    profile = /* - n * log(2.0*M_PI)/2.0 */ /* This is redundant */
        - n * log(sd) -
        (sumsq1 - 2 * mean1 * sum1 + n * mean1 * mean1) / 2.0 / sd / sd;
    if (profile > max) {
        max = profile;
        *dim = n;
    }

    return 0;
}