haskell-igraph-0.8.0: igraph/src/embedding.c
/* -*- 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(°, 0);
igraph_strength(graph, °, 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(°);
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=*/ °, /*deg2=*/ 0,
options, callback, 0, /*symmetric=*/ 1,
/*eigen=*/ 1, /*zapsmall=*/ 1));
igraph_vector_destroy(°);
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(°_in, n);
IGRAPH_VECTOR_INIT_FINALLY(°_out, n);
igraph_strength(graph, °_in, igraph_vss_all(), IGRAPH_IN, /*loops=*/ 1,
weights);
igraph_strength(graph, °_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=*/ °_in,
/*deg2=*/ °_out, options, callback,
callback_right, /*symmetric=*/ 0, /*eigen=*/ 0,
/*zapsmall=*/ 1));
igraph_vector_destroy(°_in);
igraph_vector_destroy(°_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;
}