int igraph_local_scan_0_them(const igraph_t *us, const igraph_t *them,
			     igraph_vector_t *res,
			     const igraph_vector_t *weights_them,
			     igraph_neimode_t mode) {

  igraph_t is;

  if (igraph_vcount(us) != igraph_vcount(them)) {
    IGRAPH_ERROR("Number of vertices don't match in scan-0", IGRAPH_EINVAL);
  }
  if (igraph_is_directed(us) != igraph_is_directed(them)) {
    IGRAPH_ERROR("Directedness don't match in scan-0", IGRAPH_EINVAL);
  }

  if (weights_them) {
    return igraph_i_local_scan_0_them_w(us, them, res, weights_them, mode);
  }

  igraph_intersection(&is, us, them, /*edgemap1=*/ 0, /*edgemap2=*/ 0);
  IGRAPH_FINALLY(igraph_destroy, &is);

  igraph_degree(&is, res, igraph_vss_all(), mode, IGRAPH_LOOPS);

  igraph_destroy(&is);
  IGRAPH_FINALLY_CLEAN(1);

  return 0;
}

/**
 * \function igraph_maximum_bipartite_matching
 * Calculates a maximum matching in a bipartite graph.
 *
 * A matching in a bipartite graph is a partial assignment of vertices
 * of the first kind to vertices of the second kind such that each vertex of
 * the first kind is matched to at most one vertex of the second kind and
 * vice versa, and matched vertices must be connected by an edge in the graph.
 * The size (or cardinality) of a matching is the number of edges.
 * A matching is a maximum matching if there exists no other matching with
 * larger cardinality. For weighted graphs, a maximum matching is a matching
 * whose edges have the largest possible total weight among all possible
 * matchings.
 *
 * </para><para>
 * Maximum matchings in bipartite graphs are found by the push-relabel algorithm
 * with greedy initialization and a global relabeling after every n/2 steps where
 * n is the number of vertices in the graph.
 *
 * </para><para>
 * References: Cherkassky BV, Goldberg AV, Martin P, Setubal JC and Stolfi J:
 * Augment or push: A computational study of bipartite matching and
 * unit-capacity flow algorithms. ACM Journal of Experimental Algorithmics 3,
 * 1998.
 *
 * </para><para>
 * Kaya K, Langguth J, Manne F and Ucar B: Experiments on push-relabel-based
 * maximum cardinality matching algorithms for bipartite graphs. Technical
 * Report TR/PA/11/33 of the Centre Europeen de Recherche et de Formation
 * Avancee en Calcul Scientifique, 2011.
 *
 * \param graph The input graph. It can be directed but the edge directions
 *              will be ignored.
 * \param types Boolean vector giving the vertex types of the graph.
 * \param matching_size The size of the matching (i.e. the number of matched
 *                      vertex pairs will be returned here). It may be \c NULL
 *                      if you don't need this.
 * \param matching_weight The weight of the matching if the edges are weighted,
 *                        or the size of the matching again if the edges are
 *                        unweighted. It may be \c NULL if you don't need this.
 * \param matching The matching itself. It must be a vector where element i
 *                 contains the ID of the vertex that vertex i is matched to,
 *                 or -1 if vertex i is unmatched.
 * \param weights A null pointer (=no edge weights), or a vector giving the
 *                weights of the edges. Note that the algorithm is stable
 *                only for integer weights.
 * \param eps A small real number used in equality tests in the weighted
 *            bipartite matching algorithm. Two real numbers are considered
 *            equal in the algorithm if their difference is smaller than
 *            \c eps. This is required to avoid the accumulation of numerical
 *            errors. It is advised to pass a value derived from the
 *            \c DBL_EPSILON constant in \c float.h here. If you are
 *            running the algorithm with no \c weights vector, this argument
 *            is ignored.
 * \return Error code.
 *
 * Time complexity: O(sqrt(|V|) |E|) for unweighted graphs (according to the
 * technical report referenced above), O(|V||E|) for weighted graphs.
 * 
 * \example examples/simple/igraph_maximum_bipartite_matching.c
 */
int igraph_maximum_bipartite_matching(const igraph_t* graph,
    const igraph_vector_bool_t* types, igraph_integer_t* matching_size,
    igraph_real_t* matching_weight, igraph_vector_long_t* matching,
    const igraph_vector_t* weights, igraph_real_t eps) {

  /* Sanity checks */
  if (igraph_vector_bool_size(types) < igraph_vcount(graph)) {
    IGRAPH_ERROR("types vector too short", IGRAPH_EINVAL);
  }
  if (weights && igraph_vector_size(weights) < igraph_ecount(graph)) {
    IGRAPH_ERROR("weights vector too short", IGRAPH_EINVAL);
  }

  if (weights == 0) {
    IGRAPH_CHECK(igraph_i_maximum_bipartite_matching_unweighted(graph, types,
        matching_size, matching));
    if (matching_weight != 0) {
      *matching_weight = *matching_size;
    }
    return IGRAPH_SUCCESS;
  } else {
    return igraph_i_maximum_bipartite_matching_weighted(graph, types,
        matching_size, matching_weight, matching, weights, eps);
  }
}

int igraph_adjlist_init(const igraph_t *graph, igraph_adjlist_t *al, 
			  igraph_neimode_t mode) {
  long int i;

  if (mode != IGRAPH_IN && mode != IGRAPH_OUT && mode != IGRAPH_ALL) {
    IGRAPH_ERROR("Cannot create adjlist view", IGRAPH_EINVMODE);
  }

  if (!igraph_is_directed(graph)) { mode=IGRAPH_ALL; }

  al->length=igraph_vcount(graph);
  al->adjs=igraph_Calloc(al->length, igraph_vector_t);
  if (al->adjs == 0) {
    IGRAPH_ERROR("Cannot create adjlist view", IGRAPH_ENOMEM);
  }

  IGRAPH_FINALLY(igraph_adjlist_destroy, al);
  for (i=0; i<al->length; i++) {
    IGRAPH_ALLOW_INTERRUPTION();
    IGRAPH_CHECK(igraph_vector_init(&al->adjs[i], 0));
    IGRAPH_CHECK(igraph_neighbors(graph, &al->adjs[i], i, mode));
  }

  IGRAPH_FINALLY_CLEAN(1);
  return 0;
}

int igraph_sample_dirichlet(igraph_integer_t n, const igraph_vector_t *alpha,
			    igraph_matrix_t *res) {

  igraph_integer_t len=igraph_vector_size(alpha);
  igraph_integer_t i;
  igraph_vector_t vec;

  if (n < 0) {
    IGRAPH_ERROR("Number of samples should be non-negative",
		 IGRAPH_EINVAL);
  }
  if (len < 2) {
    IGRAPH_ERROR("Dirichlet parameter vector too short, must "
		 "have at least two entries", IGRAPH_EINVAL);
  }
  if (igraph_vector_min(alpha) <= 0) {
    IGRAPH_ERROR("Dirichlet concentration parameters must be positive",
		 IGRAPH_EINVAL);
  }

  IGRAPH_CHECK(igraph_matrix_resize(res, len, n));

  RNG_BEGIN();

  for (i = 0; i < n; i++) {
    igraph_vector_view(&vec, &MATRIX(*res, 0, i), len);
    igraph_rng_get_dirichlet(igraph_rng_default(), alpha, &vec);
  }

  RNG_END();

  return 0;
}

int igraph_get_eid(const igraph_t *graph, igraph_integer_t *eid,
		   igraph_integer_t pfrom, igraph_integer_t pto,
		   igraph_bool_t directed) {

  long int from=pfrom, to=pto;
  long int nov=igraph_vcount(graph);

  if (from < 0 || to < 0 || from > nov-1 || to > nov-1) {
    IGRAPH_ERROR("cannot get edge id", IGRAPH_EINVVID);
  }

  *eid=-1;
  if (igraph_is_directed(graph)) {

    /* Directed graph */
    FIND_DIRECTED_EDGE(graph,from,to,eid);
    if (!directed && *eid < 0) {
      FIND_DIRECTED_EDGE(graph,to,from,eid);
    }
    
  } else {

    /* Undirected graph, they only have one mode */
    FIND_UNDIRECTED_EDGE(graph,from,to,eid);

  }

  if (*eid < 0) {
    IGRAPH_ERROR("Cannot get edge id, no such edge", IGRAPH_EINVAL);
  }
  
  return IGRAPH_SUCCESS;  
}

int igraph_inclist_init(const igraph_t *graph, 
			      igraph_inclist_t *il, 
			      igraph_neimode_t mode) {
  long int i;

  if (mode != IGRAPH_IN && mode != IGRAPH_OUT && mode != IGRAPH_ALL) {
    IGRAPH_ERROR("Cannot create incidence list view", IGRAPH_EINVMODE);
  }

  if (!igraph_is_directed(graph)) { mode=IGRAPH_ALL; }

  il->length=igraph_vcount(graph);
  il->incs=igraph_Calloc(il->length, igraph_vector_t);
  if (il->incs == 0) {
    IGRAPH_ERROR("Cannot create incidence list view", IGRAPH_ENOMEM);
  }

  IGRAPH_FINALLY(igraph_inclist_destroy, il);  
  for (i=0; i<il->length; i++) {
    IGRAPH_ALLOW_INTERRUPTION();
    IGRAPH_CHECK(igraph_vector_init(&il->incs[i], 0));
    IGRAPH_CHECK(igraph_incident(graph, &il->incs[i], i, mode));
  }
  
  IGRAPH_FINALLY_CLEAN(1);
  return 0;
}

/** 
 * \ingroup interface
 * \function igraph_empty_attrs
 * \brief Creates an empty graph with some vertices, no edges and some graph attributes.
 *
 * </para><para>
 * Use this instead of \ref igraph_empty() if you wish to add some graph
 * attributes right after initialization. This function is currently
 * not very interesting for the ordinary user, just supply 0 here or 
 * use \ref igraph_empty().
 * \param graph Pointer to a not-yet initialized graph object.
 * \param n The number of vertices in the graph, a non-negative
 *          integer number is expected.
 * \param directed Whether the graph is directed or not.
 * \param attr The attributes. 
 * \return Error code:
 *         \c IGRAPH_EINVAL: invalid number of vertices.
 * 
 * Time complexity: O(|V|) for a graph with
 * |V| vertices (and no edges).
 */
int igraph_empty_attrs(igraph_t *graph, igraph_integer_t n, igraph_bool_t directed, void* attr) {

  if (n<0) {
    IGRAPH_ERROR("cannot create empty graph with negative number of vertices",
		  IGRAPH_EINVAL);
  }
  
  if (!IGRAPH_FINITE(n)) {
    IGRAPH_ERROR("number of vertices is not finite (NA, NaN or Inf)", IGRAPH_EINVAL);
  }

  graph->n=0;
  graph->directed=directed;
  IGRAPH_VECTOR_INIT_FINALLY(&graph->from, 0);
  IGRAPH_VECTOR_INIT_FINALLY(&graph->to, 0);
  IGRAPH_VECTOR_INIT_FINALLY(&graph->oi, 0);
  IGRAPH_VECTOR_INIT_FINALLY(&graph->ii, 0);
  IGRAPH_VECTOR_INIT_FINALLY(&graph->os, 1);
  IGRAPH_VECTOR_INIT_FINALLY(&graph->is, 1);

  VECTOR(graph->os)[0]=0;
  VECTOR(graph->is)[0]=0;

  /* init attributes */
  graph->attr=0;
  IGRAPH_CHECK(igraph_i_attribute_init(graph, attr));

  /* add the vertices */
  IGRAPH_CHECK(igraph_add_vertices(graph, n, 0));
  
  IGRAPH_FINALLY_CLEAN(6);
  return 0;
}

// Choose vertices in the order of their IDs.
static int igraph_i_kleitman_wang_index(const igraph_vector_t *outdeg, const igraph_vector_t *indeg, igraph_vector_t *edges) {
    long n = igraph_vector_size(indeg); // number of vertices

    long ec = 0; // number of edges added so far

    typedef std::list<vbd_pair> vlist;
    vlist vertices;
    for (int i=0; i < n; ++i)
        vertices.push_back(vbd_pair(i, bidegree(VECTOR(*indeg)[i], VECTOR(*outdeg)[i])));

    std::vector<vlist::iterator> pointers;
    pointers.reserve(n);
    for (vlist::iterator it = vertices.begin(); it != vertices.end(); ++it)
        pointers.push_back(it);

    for (std::vector<vlist::iterator>::iterator pt = pointers.begin(); pt != pointers.end(); ++pt) {
        // sort vertices by (in, out) degree pairs in decreasing order
        // note: std::list::sort does a stable sort
        vertices.sort(degree_greater<vbd_pair>);

        // choose a vertex the out-stubs of which will be connected
        vbd_pair &vd = **pt;

        if (vd.degree.second == 0)
            continue;

        if (vd.degree.first < 0 || vd.degree.second < 0)
            IGRAPH_ERROR("Vertex degrees must be positive", IGRAPH_EINVAL);

        int k = 0;
        vlist::iterator it;
        for (it = vertices.begin();
             k != vd.degree.second && it != vertices.end();
             ++it)
        {
            if (it->vertex == vd.vertex)
                continue;

            if (--(it->degree.first) < 0)
                goto fail;

            VECTOR(*edges)[2*(ec+k)] = vd.vertex;
            VECTOR(*edges)[2*(ec+k)+1] = it->vertex;

            ++k;
        }
        if (it == vertices.end() && k < vd.degree.second)
            goto fail;

        ec += vd.degree.second;
        vd.degree.second = 0;
    }

    return IGRAPH_SUCCESS;

fail:
    IGRAPH_ERROR("The given directed degree sequence is not realizable", IGRAPH_EINVAL);
}

int igraph_adjlist_init_complementer(const igraph_t *graph,
				       igraph_adjlist_t *al, 
				       igraph_neimode_t mode,
				       igraph_bool_t loops) {
  long int i, j, k, n;
  igraph_bool_t* seen;
  igraph_vector_t vec;

  if (mode != IGRAPH_IN && mode != IGRAPH_OUT && mode != IGRAPH_ALL) {
    IGRAPH_ERROR("Cannot create complementer adjlist view", IGRAPH_EINVMODE);
  }

  if (!igraph_is_directed(graph)) { mode=IGRAPH_ALL; }

  al->length=igraph_vcount(graph);
  al->adjs=igraph_Calloc(al->length, igraph_vector_t);
  if (al->adjs == 0) {
    IGRAPH_ERROR("Cannot create complementer adjlist view", IGRAPH_ENOMEM);
  }

  IGRAPH_FINALLY(igraph_adjlist_destroy, al);

  n=al->length;
  seen=igraph_Calloc(n, igraph_bool_t);
  if (seen==0) {
    IGRAPH_ERROR("Cannot create complementer adjlist view", IGRAPH_ENOMEM);
  }
  IGRAPH_FINALLY(igraph_free, seen);

  IGRAPH_VECTOR_INIT_FINALLY(&vec, 0);

  for (i=0; i<al->length; i++) {
    IGRAPH_ALLOW_INTERRUPTION();
    igraph_neighbors(graph, &vec, i, mode);
    memset(seen, 0, sizeof(igraph_bool_t)*al->length);
    n=al->length;
    if (!loops) { seen[i] = 1; n--; }
    for (j=0; j<igraph_vector_size(&vec); j++) {
      if (! seen [ (long int) VECTOR(vec)[j] ] ) {
	n--;
	seen[ (long int) VECTOR(vec)[j] ] = 1;
      }
    }
    IGRAPH_CHECK(igraph_vector_init(&al->adjs[i], n));
    for (j=0, k=0; k<n; j++) {
      if (!seen[j]) {
	VECTOR(al->adjs[i])[k++] = j;
      }
    }
  }

  igraph_Free(seen);
  igraph_vector_destroy(&vec);
  IGRAPH_FINALLY_CLEAN(3);
  return 0;
}

// Generate undirected realization as edge-list.
// If largest=true, always choose the vertex with the largest remaining degree to connect up next.
// Otherwise, always choose the one with the smallest remaining degree.
static int igraph_i_havel_hakimi(const igraph_vector_t *deg, igraph_vector_t *edges, bool largest) {
    long n = igraph_vector_size(deg);

    long ec = 0; // number of edges added so far

    std::vector<vd_pair> vertices;
    vertices.reserve(n);
    for (int i=0; i < n; ++i)
        vertices.push_back(vd_pair(i, VECTOR(*deg)[i]));

    while (! vertices.empty()) {
        if (largest)
            std::stable_sort(vertices.begin(), vertices.end(), degree_less<vd_pair>);
        else
            std::stable_sort(vertices.begin(), vertices.end(), degree_greater<vd_pair>);

        // take the next vertex to be connected up
        vd_pair vd = vertices.back();
        vertices.pop_back();

        if (vd.degree < 0)
            IGRAPH_ERROR("Vertex degrees must be positive", IGRAPH_EINVAL);

        if (vd.degree == 0)
            continue;

        if (vertices.size() < size_t(vd.degree))
            goto fail;

        if (largest) {
            for (int i=0; i < vd.degree; ++i) {
                if (--(vertices[vertices.size() - 1 - i].degree) < 0)
                    goto fail;

                VECTOR(*edges)[2*(ec+i)] = vd.vertex;
                VECTOR(*edges)[2*(ec+i)+1] = vertices[vertices.size() - 1 - i].vertex;
            }
        } else {
            // this loop can only be reached if all zero-degree nodes have already been removed
            // therefore decrementing remaining degrees is safe
            for (int i=0; i < vd.degree; ++i) {
                vertices[i].degree--;

                VECTOR(*edges)[2*(ec+i)] = vd.vertex;
                VECTOR(*edges)[2*(ec+i)+1] = vertices[i].vertex;
            }
        }

        ec += vd.degree;
    }

    return IGRAPH_SUCCESS;

fail:
    IGRAPH_ERROR("The given degree sequence is not realizable", IGRAPH_EINVAL);
}

int igraph_eigen_matrix_symmetric(const igraph_matrix_t *A,
				  const igraph_sparsemat_t *sA,
				  igraph_arpack_function_t *fun, int n,
				  void *extra,
				  igraph_eigen_algorithm_t algorithm,
				  const igraph_eigen_which_t *which,
				  igraph_arpack_options_t *options,
				  igraph_arpack_storage_t *storage,
				  igraph_vector_t *values, 
				  igraph_matrix_t *vectors) {

  IGRAPH_CHECK(igraph_i_eigen_checks(A, sA, fun, n));
  
  if (which->pos != IGRAPH_EIGEN_LM && 
      which->pos != IGRAPH_EIGEN_SM && 
      which->pos != IGRAPH_EIGEN_LA && 
      which->pos != IGRAPH_EIGEN_SA && 
      which->pos != IGRAPH_EIGEN_BE && 
      which->pos != IGRAPH_EIGEN_ALL && 
      which->pos != IGRAPH_EIGEN_INTERVAL && 
      which->pos != IGRAPH_EIGEN_SELECT) {
    IGRAPH_ERROR("Invalid 'pos' position in 'which'", IGRAPH_EINVAL);
  }

  switch (algorithm) {
  case IGRAPH_EIGEN_AUTO:
    if (which->howmany==n || n < 100) {
      IGRAPH_CHECK(igraph_i_eigen_matrix_symmetric_lapack(A, sA, fun, n,
							  extra, which, 
							  values, vectors));
    } else {
      IGRAPH_CHECK(igraph_i_eigen_matrix_symmetric_arpack(A, sA, fun, n, 
							  extra, which, 
							  options, storage,
							  values, vectors));
    }
    break;
  case IGRAPH_EIGEN_LAPACK:
    IGRAPH_CHECK(igraph_i_eigen_matrix_symmetric_lapack(A, sA, fun, n ,extra,
							which, values, 
							vectors));
    break;
  case IGRAPH_EIGEN_ARPACK:
    IGRAPH_CHECK(igraph_i_eigen_matrix_symmetric_arpack(A, sA, fun, n, extra,
							which, options, 
							storage,
							values, vectors));
    break;
  default:
    IGRAPH_ERROR("Unknown 'algorithm'", IGRAPH_EINVAL);
  }
    
  return 0;
}

int igraph_i_maximal_or_largest_cliques_or_indsets(const igraph_t *graph,
                                        igraph_vector_ptr_t *res,
                                        igraph_integer_t *clique_number,
                                        igraph_bool_t keep_only_largest,
                                        igraph_bool_t complementer) {
  igraph_i_max_ind_vsets_data_t clqdata;
  long int no_of_nodes = igraph_vcount(graph), i;

  if (igraph_is_directed(graph))
    IGRAPH_WARNING("directionality of edges is ignored for directed graphs");

  clqdata.matrix_size=no_of_nodes;
  clqdata.keep_only_largest=keep_only_largest;

  if (complementer)
    IGRAPH_CHECK(igraph_adjlist_init_complementer(graph, &clqdata.adj_list, IGRAPH_ALL, 0));
  else
    IGRAPH_CHECK(igraph_adjlist_init(graph, &clqdata.adj_list, IGRAPH_ALL));
  IGRAPH_FINALLY(igraph_adjlist_destroy, &clqdata.adj_list);

  clqdata.IS = igraph_Calloc(no_of_nodes, igraph_integer_t);
  if (clqdata.IS == 0)
    IGRAPH_ERROR("igraph_i_maximal_or_largest_cliques_or_indsets failed", IGRAPH_ENOMEM);
  IGRAPH_FINALLY(igraph_free, clqdata.IS);

  IGRAPH_VECTOR_INIT_FINALLY(&clqdata.deg, no_of_nodes);
  for (i=0; i<no_of_nodes; i++)
    VECTOR(clqdata.deg)[i] = igraph_vector_size(igraph_adjlist_get(&clqdata.adj_list, i));

  clqdata.buckets = igraph_Calloc(no_of_nodes+1, igraph_set_t);
  if (clqdata.buckets == 0)
    IGRAPH_ERROR("igraph_maximal_or_largest_cliques_or_indsets failed", IGRAPH_ENOMEM);
  IGRAPH_FINALLY(igraph_i_free_set_array, clqdata.buckets);

  for (i=0; i<no_of_nodes; i++)
    IGRAPH_CHECK(igraph_set_init(&clqdata.buckets[i], 0));

  if (res) igraph_vector_ptr_clear(res);
  
  /* Do the show */
  clqdata.largest_set_size=0;
  IGRAPH_CHECK(igraph_i_maximal_independent_vertex_sets_backtrack(graph, res, &clqdata, 0));

  /* Cleanup */
  for (i=0; i<no_of_nodes; i++) igraph_set_destroy(&clqdata.buckets[i]);
  igraph_adjlist_destroy(&clqdata.adj_list);
  igraph_vector_destroy(&clqdata.deg);
  igraph_free(clqdata.IS);
  igraph_free(clqdata.buckets);
  IGRAPH_FINALLY_CLEAN(4);

  if (clique_number) *clique_number = clqdata.largest_set_size;
  return 0;
}

int igraph_random_walk(const igraph_t *graph, igraph_vector_t *walk,
		       igraph_integer_t start, igraph_neimode_t mode,
		       igraph_integer_t steps,
		       igraph_random_walk_stuck_t stuck) {

  /* TODO:
     - multiple walks potentially from multiple start vertices
     - weights
  */

  igraph_lazy_adjlist_t adj;
  igraph_integer_t vc = igraph_vcount(graph);
  igraph_integer_t i;

  if (start < 0 || start >= vc) {
    IGRAPH_ERROR("Invalid start vertex", IGRAPH_EINVAL);
  }
  if (steps < 0) {
    IGRAPH_ERROR("Invalid number of steps", IGRAPH_EINVAL);
  }

  IGRAPH_CHECK(igraph_lazy_adjlist_init(graph, &adj, mode,
					IGRAPH_DONT_SIMPLIFY));
  IGRAPH_FINALLY(igraph_lazy_adjlist_destroy, &adj);

  IGRAPH_CHECK(igraph_vector_resize(walk, steps));

  RNG_BEGIN();

  VECTOR(*walk)[0] = start;
  for (i = 1; i < steps; i++) {
    igraph_vector_t *neis;
    igraph_integer_t nn;
    neis = igraph_lazy_adjlist_get(&adj, start);
    nn = igraph_vector_size(neis);

    if (IGRAPH_UNLIKELY(nn == 0)) {
      igraph_vector_resize(walk, i);
      if (stuck == IGRAPH_RANDOM_WALK_STUCK_RETURN) {
	break;
      } else {
	IGRAPH_ERROR("Random walk got stuck", IGRAPH_ERWSTUCK);
      }
    }
    start = VECTOR(*walk)[i] = VECTOR(*neis)[ RNG_INTEGER(0, nn - 1) ];
  }

  RNG_END();

  igraph_lazy_adjlist_destroy(&adj);
  IGRAPH_FINALLY_CLEAN(1);

  return 0;
}

// Choose vertices in the order of their IDs.
static int igraph_i_havel_hakimi_index(const igraph_vector_t *deg, igraph_vector_t *edges) {
    long n = igraph_vector_size(deg);

    long ec = 0; // number of edges added so far

    typedef std::list<vd_pair> vlist;
    vlist vertices;
    for (int i=0; i < n; ++i)
        vertices.push_back(vd_pair(i, VECTOR(*deg)[i]));

    std::vector<vlist::iterator> pointers;
    pointers.reserve(n);
    for (vlist::iterator it = vertices.begin(); it != vertices.end(); ++it)
        pointers.push_back(it);

    for (std::vector<vlist::iterator>::iterator pt = pointers.begin(); pt != pointers.end(); ++pt) {
        vertices.sort(degree_greater<vd_pair>);

        vd_pair vd = **pt;
        vertices.erase(*pt);

        if (vd.degree < 0)
            IGRAPH_ERROR("Vertex degrees must be positive", IGRAPH_EINVAL);

        if (vd.degree == 0)
            continue;

        int k;
        vlist::iterator it;
        for (it = vertices.begin(), k = 0;
             k != vd.degree && it != vertices.end();
             ++it, ++k)
        {
            if (--(it->degree) < 0)
                goto fail;

            VECTOR(*edges)[2*(ec+k)] = vd.vertex;
            VECTOR(*edges)[2*(ec+k)+1] = it->vertex;
        }
        if (it == vertices.end() && k < vd.degree)
            goto fail;

        ec += vd.degree;
    }

    return IGRAPH_SUCCESS;

fail:
    IGRAPH_ERROR("The given degree sequence is not realizable", IGRAPH_EINVAL);
}

int igraph_create_bipartite(igraph_t *graph, const igraph_vector_bool_t *types,
			    const igraph_vector_t *edges, 
			    igraph_bool_t directed) {

  igraph_integer_t no_of_nodes=
    (igraph_integer_t) igraph_vector_bool_size(types);
  long int no_of_edges=igraph_vector_size(edges);
  igraph_real_t min_edge=0, max_edge=0;
  igraph_bool_t min_type=0, max_type=0;
  long int i;

  if (no_of_edges % 2 != 0) {
    IGRAPH_ERROR("Invalid (odd) edges vector", IGRAPH_EINVEVECTOR);
  }
  no_of_edges /= 2;
  
  if (no_of_edges != 0) {
    igraph_vector_minmax(edges, &min_edge, &max_edge);
  }
  if (min_edge < 0 || max_edge >= no_of_nodes) {
    IGRAPH_ERROR("Invalid (negative) vertex id", IGRAPH_EINVVID);
  }

  /* Check types vector */
  if (no_of_nodes != 0) {
    igraph_vector_bool_minmax(types, &min_type, &max_type);
    if (min_type < 0 || max_type > 1) {
      IGRAPH_WARNING("Non-binary type vector when creating a bipartite graph");
    }
  }

  /* Check bipartiteness */
  for (i=0; i<no_of_edges*2; i+=2) {
    long int from=(long int) VECTOR(*edges)[i];
    long int to=(long int) VECTOR(*edges)[i+1];
    long int t1=VECTOR(*types)[from];
    long int t2=VECTOR(*types)[to];
    if ( (t1 && t2) || (!t1 && !t2) ) {
      IGRAPH_ERROR("Invalid edges, not a bipartite graph", IGRAPH_EINVAL);
    }
  }
  
  IGRAPH_CHECK(igraph_empty(graph, no_of_nodes, directed));
  IGRAPH_FINALLY(igraph_destroy, graph);
  IGRAPH_CHECK(igraph_add_edges(graph, edges, 0));
  
  IGRAPH_FINALLY_CLEAN(1);
  return 0;
}

int igraph_adjacent(const igraph_t *graph, igraph_vector_t *eids, 
		    igraph_integer_t pnode, igraph_neimode_t mode) {
  
  long int length=0, idx=0;   
  long int no_of_edges;
  long int i, j;

  long int node=pnode;

  if (node<0 || node>igraph_vcount(graph)-1) {
    IGRAPH_ERROR("cannot get neighbors", IGRAPH_EINVVID);
  }
  if (mode != IGRAPH_OUT && mode != IGRAPH_IN && 
      mode != IGRAPH_ALL) {
    IGRAPH_ERROR("cannot get neighbors", IGRAPH_EINVMODE);
  }

  no_of_edges=igraph_vector_size(&graph->from);
  if (! graph->directed) {
    mode=IGRAPH_ALL;
  }

  /* Calculate needed space first & allocate it*/

  if (mode & IGRAPH_OUT) {
    length += (VECTOR(graph->os)[node+1] - VECTOR(graph->os)[node]);
  }
  if (mode & IGRAPH_IN) {
    length += (VECTOR(graph->is)[node+1] - VECTOR(graph->is)[node]);
  }
  
  IGRAPH_CHECK(igraph_vector_resize(eids, length));
  
  if (mode & IGRAPH_OUT) {
    j=VECTOR(graph->os)[node+1];
    for (i=VECTOR(graph->os)[node]; i<j; i++) {
      VECTOR(*eids)[idx++] = VECTOR(graph->oi)[i];
    }
  }
  if (mode & IGRAPH_IN) {
    j=VECTOR(graph->is)[node+1];
    for (i=VECTOR(graph->is)[node]; i<j; i++) {
      VECTOR(*eids)[idx++] = VECTOR(graph->ii)[i];
    }
  }

  return 0;
}

int igraph_i_cliquer_callback(const igraph_t *graph,
                    igraph_integer_t min_size, igraph_integer_t max_size,
                    igraph_clique_handler_t *cliquehandler_fn, void *arg)
{
    graph_t *g;
    struct callback_data cd;
    igraph_integer_t vcount = igraph_vcount(graph);

    if (vcount == 0)
        return IGRAPH_SUCCESS;

    if (min_size <= 0) min_size = 1;
    if (max_size <= 0) max_size = 0;

    if (max_size > 0 && max_size < min_size)
        IGRAPH_ERROR("max_size must not be smaller than min_size", IGRAPH_EINVAL);

    igraph_to_cliquer(graph, &g);
    IGRAPH_FINALLY(graph_free, g);

    cd.handler = cliquehandler_fn;
    cd.arg = arg;
    igraph_cliquer_opt.user_data = &cd;
    igraph_cliquer_opt.user_function = &callback_callback;

    CLIQUER_INTERRUPTABLE(clique_unweighted_find_all(g, min_size, max_size, /* maximal= */ FALSE, &igraph_cliquer_opt));

    graph_free(g);
    IGRAPH_FINALLY_CLEAN(1);

    return IGRAPH_SUCCESS;
}

int igraph_i_cliquer_cliques(const igraph_t *graph, igraph_vector_ptr_t *res,
                    igraph_integer_t min_size, igraph_integer_t max_size)
{
    graph_t *g;
    igraph_integer_t vcount = igraph_vcount(graph);

    if (vcount == 0) {
        igraph_vector_ptr_clear(res);
        return IGRAPH_SUCCESS;
    }

    if (min_size <= 0) min_size = 1;
    if (max_size <= 0) max_size = 0;

    if (max_size > 0 && max_size < min_size)
        IGRAPH_ERROR("max_size must not be smaller than min_size", IGRAPH_EINVAL);

    igraph_to_cliquer(graph, &g);
    IGRAPH_FINALLY(graph_free, g);

    igraph_vector_ptr_clear(res);
    igraph_cliquer_opt.user_data = res;
    igraph_cliquer_opt.user_function = &collect_cliques_callback;

    IGRAPH_FINALLY(free_clique_list, res);
    CLIQUER_INTERRUPTABLE(clique_unweighted_find_all(g, min_size, max_size, /* maximal= */ FALSE, &igraph_cliquer_opt));
    IGRAPH_FINALLY_CLEAN(1);

    graph_free(g);
    IGRAPH_FINALLY_CLEAN(1);

    return IGRAPH_SUCCESS;
}

int igraph_i_largest_cliques_store(const igraph_vector_t* clique, void* data, igraph_bool_t* cont) {
  igraph_vector_ptr_t* result = (igraph_vector_ptr_t*)data;
  igraph_vector_t* vec;
  long int i, n;

  /* Is the current clique at least as large as the others that we have found? */
  if (!igraph_vector_ptr_empty(result)) {
    n = igraph_vector_size(clique);
    if (n < igraph_vector_size(VECTOR(*result)[0]))
      return IGRAPH_SUCCESS;

    if (n > igraph_vector_size(VECTOR(*result)[0])) {
      for (i = 0; i < igraph_vector_ptr_size(result); i++)
        igraph_vector_destroy(VECTOR(*result)[i]);
      igraph_vector_ptr_free_all(result);
      igraph_vector_ptr_resize(result, 0);
    }
  }

  vec = igraph_Calloc(1, igraph_vector_t);
  if (vec == 0)
    IGRAPH_ERROR("cannot allocate memory for storing next clique", IGRAPH_ENOMEM);

  IGRAPH_CHECK(igraph_vector_copy(vec, clique));
  IGRAPH_CHECK(igraph_vector_ptr_push_back(result, vec));

  return IGRAPH_SUCCESS;
}

int igraph_running_mean(const igraph_vector_t *data, igraph_vector_t *res, 
			igraph_integer_t binwidth) {

  double sum=0;
  long int i;

  /* Check */
  if (igraph_vector_size(data) < binwidth) {
    IGRAPH_ERROR("Vector too short for this binwidth", IGRAPH_EINVAL); 
  }

  /* Memory for result */

  IGRAPH_CHECK(igraph_vector_resize(res, (long int)(igraph_vector_size(data)-binwidth+1)));
  
  /* Initial bin */
  for (i=0; i<binwidth; i++) {
    sum += VECTOR(*data)[i];
  }
  
  VECTOR(*res)[0]=sum/binwidth;
  
  for (i=1; i<igraph_vector_size(data)-binwidth+1; i++) {
    IGRAPH_ALLOW_INTERRUPTION();
    sum -= VECTOR(*data)[i-1];
    sum += VECTOR(*data)[ (long int)(i+binwidth-1)];
    VECTOR(*res)[i] = sum/binwidth;
  }
  
  return 0;
}