motifcluster.cpp File Reference

#include "stdafx.h"
#include "motifcluster.h"

Go to the source code of this file.

Macros
#define	F77_NAME(name)   name

Functions
void F77_NAME()	dsaupd (int *ido, char *bmat, int *n, char *which, int *nev, double *tol, double *resid, int *ncv, double *V, int *ldv, int *iparam, int *ipntr, double *workd, double *workl, int *lworkl, int *info)
void F77_NAME()	dseupd (int *rvec, char *HowMny, int *select, double *d, double *Z, int *ldz, double *sigma, char *bmat, int *n, char *which, int *nev, double *tol, double *resid, int *ncv, double *V, int *ldv, int *iparam, int *ipntr, double *workd, double *workl, int *lworkl, int *info)
static void	IncrementWeight (int i, int j, WeightVH &weights)
static PUNGraph	UnweightedGraphRepresentation (const WeightVH &weights)
static void	MapIdsToFirstN (const TIntV &ids, THash< TInt, TInt > &id_map, TIntV &rev_id_map)
static void	Sweep (const TSparseColMatrix &W, const TFltV &fvec, TFltV &conds, TIntV &order)
void	SymeigsSmallest (const TSparseColMatrix &A, int nev, TFltV &evals, TFullColMatrix &evecs, double tol, int maxiter)

Macro Definition Documentation

#define F77_NAME

(

name

)

name

Definition at line 7 of file motifcluster.cpp.

Function Documentation

void F77_NAME() dsaupd	(	int *	ido,
		char *	bmat,
		int *	n,
		char *	which,
		int *	nev,
		double *	tol,
		double *	resid,
		int *	ncv,
		double *	V,
		int *	ldv,
		int *	iparam,
		int *	ipntr,
		double *	workd,
		double *	workl,
		int *	lworkl,
		int *	info
	)

void F77_NAME() dseupd	(	int *	rvec,
		char *	HowMny,
		int *	select,
		double *	d,
		double *	Z,
		int *	ldz,
		double *	sigma,
		char *	bmat,
		int *	n,
		char *	which,
		int *	nev,
		double *	tol,
		double *	resid,
		int *	ncv,
		double *	V,
		int *	ldv,
		int *	iparam,
		int *	ipntr,
		double *	workd,
		double *	workl,
		int *	lworkl,
		int *	info
	)

static void IncrementWeight ( int  i, int  j, WeightVH &  weights  )

static

Definition at line 29 of file motifcluster.cpp.

                                                             {
  int minval = MIN(i, j);
  int maxval = MAX(i, j);
  THash<TInt, TInt>& edge_weights = weights[minval];
  edge_weights(maxval) += 1;
}

static void MapIdsToFirstN ( const TIntV &  ids, THash< TInt, TInt > &  id_map, TIntV &  rev_id_map  )

static

Definition at line 796 of file motifcluster.cpp.

                                              {
  id_map = THash<TInt, TInt>();
  rev_id_map = TIntV(ids.Len());
  for (int i = 0; i < ids.Len(); i++) {
    int id = ids[i];
    if (id_map.IsKey(id)) {
      TExcept::Throw("List of ids is not unique.");
    }
    id_map(id) = i;
    rev_id_map[i] = id;
  }
}

static void Sweep ( const TSparseColMatrix &  W, const TFltV &  fvec, TFltV &  conds, TIntV &  order  )

static

Definition at line 813 of file motifcluster.cpp.

                                {
  // Get ordering of nodes
  
  TVec< TKeyDat<TFlt, TInt> > fvec_inds(fvec.Len());
  for (int i = 0; i < fvec.Len(); i++) {
    fvec_inds[i] = TKeyDat<TFlt, TInt>(fvec[i], i);
  }
  fvec_inds.Sort();
  order = TIntV(fvec.Len());
  TIntV rank(fvec.Len());
  for (int i = 0; i < fvec.Len(); i++) {
    order[i] = fvec_inds[i].Dat;
    rank[order[i]] = i;
  }

  // Sweep by adjusting cut and volume
  conds = TFltV(order.Len() - 1);
  double cut = 0;
  double vol = 0;
  double total_vol = 0;
  for (int ind = 0; ind < order.Len(); ind++) {
    const TIntFltKdV& nbr_weights = W.ColSpVV[ind];
    for (TIntFltKdV::TIter it = nbr_weights.BegI(); it < nbr_weights.EndI();
         it++) {
      total_vol += it->Dat;
    }
  }
  double vol_comp = total_vol;
  
  for (int ind = 0; ind < order.Len() - 1; ind++) {
    int node = order[ind];
    const TIntFltKdV& nbr_weights = W.ColSpVV[node];
    for (TIntFltKdV::TIter it = nbr_weights.BegI(); it < nbr_weights.EndI();
         it++) {
      int nbr = it->Key;
      if (node == nbr) { continue; }
      double val = it->Dat;
      // Adjust the cut amount
      if (rank[nbr] > ind) {
        // nbr is on the other side: add to the cut
        cut += val;
      } else {
        // now on the same side as nbr: subtract from the cut
        cut -= val;
      }
      vol += val;
      vol_comp -= val;
    }

    double mvol = MIN(vol, vol_comp);
    if (mvol <= 0.0) {
      TExcept::Throw("Nonpositive set volume.");
    }
    conds[ind] = cut / mvol;
  }
}

void SymeigsSmallest	(	const TSparseColMatrix &	A,
		int	nev,
		TFltV &	evals,
		TFullColMatrix &	evecs,
		double	tol,
		int	maxiter
	)

Definition at line 955 of file motifcluster.cpp.

                                                                     {
  // type of problem
  int mode = 1;
  // communication variable
  int ido = 0;
  // standard eigenvalue problem
  char bmat = 'I';
  // dimension of problem
  int n = A.GetCols();
  // type of eigenvector (smallest algebraic)  
  char which[] = "SA";
  // Space for residual
  double *resid = new double[n];
  // Number of lanczos vectors  
  int ncv = MIN(MAX(2 * nev, 20), n - 1);
  // Space for Lanczos basis vectors
  double *V = new double[ncv * n];
  
  // dimension of basis vectors
  int ldv = n;
  // Various input / output data
  int *iparam = new int[11];
  // exact shifts  
  iparam[0] = 1;
  // maximum number of iterations
  iparam[2] = maxiter;
  // mode
  iparam[6] = mode;
  // Output data
  int *ipntr = new int[11];
  // work array
  double *workd = new double[3 * n];
  // output workspace
  int lworkl = ncv * (ncv + 8);
  double *workl = new double[lworkl];
  // Communication information
  int info = 0;

  TFltV Ax(n);
  // Communication loop.  We keep applying A * x until ARPACK tells us to stop.
  while (true) {
    F77_NAME(dsaupd)(&ido, &bmat, &n, &which[0], &nev, &tol, &resid[0], &ncv,
                     &V[0], &ldv, &iparam[0], &ipntr[0], &workd[0], &workl[0],
                     &lworkl, &info);
    TFltV result(n);
    result.PutAll(0.0);
    double *load = &workd[ipntr[0] - 1];
    for (int i = 0; i < n; i++) {
      result[i] = load[i];
    }

    if (ido == 1) {
      // Need another matrix-vector product.
      A.Multiply(result, Ax);
      double *store = &workd[ipntr[1] - 1];
      for (int i = 0; i < n; i++) {
        store[i] = Ax[i];
      }
    } else if (ido == 99) {
      // We have finished.
      break;
    } else {
      // There was an error.
      TExcept::Throw("ARPACK error");
    }
  }
  
  // Extract the eigenvectors and eigenvalues
  
  // compute Ritz vectors
  int rvec(true);
  // Form of Ritz vectors
  char howmny = 'A';
  // select doesn't matter with howmny set to A
  int *select = new int[ncv];
  // Store Ritz values
  double *d = new double[nev];
  // Shift
  double sigmar = 0;

  F77_NAME(dseupd)(&rvec, &howmny, select, &d[0], &V[0], &ldv, &sigmar, &bmat,
                   &n, &which[0], &nev, &tol, &resid[0], &ncv, &V[0], &ldv,
                   &iparam[0], &ipntr[0], &workd[0], &workl[0], &lworkl, &info);
  if (info != 0) {
    TExcept::Throw("ARPACK error");
  }

  // Get eigenvalues and eigenvalues in sorted order
  TVec< TKeyDat<TFlt, TInt> > d_inds(nev);
  for (int i = 0; i < nev; i++) {
    d_inds[i] = TKeyDat<TFlt, TInt>(d[i], i);
  }
  d_inds.Sort();

  evals = TFltV(nev);
  evecs.ColV = TVec<TFltV>(nev);
  for (int i = 0; i < nev; i++) {
    double eval = d_inds[i].Key;
    int ind = d_inds[i].Dat;
    evals[ind] = eval;
    TFltV& col = evecs.ColV[ind];
    for (int j = ind * n; j < ind * n + n; j++) {
      col.Add(V[j]);
    }
  }
  evecs.RowN = n;
  evecs.ColN = nev;

  // Clean up
  delete[] resid;
  delete[] V;
  delete[] iparam;
  delete[] ipntr;
  delete[] workd;
  delete[] workl;
  delete[] select;
  delete[] d;
}

static PUNGraph UnweightedGraphRepresentation ( const WeightVH &  weights )

static

Definition at line 772 of file motifcluster.cpp.

                                                                       {
  int num_edges = 0;
  for (int i = 0; i < weights.Len(); i++) {
    num_edges += weights[i].Len();
  }
  PUNGraph graph = TUNGraph::New(weights.Len(), num_edges);
  for (int i = 0; i < weights.Len(); i++) {
    graph->AddNode(i);
  }
  for (int i = 0; i < weights.Len(); i++) {
    const THash<TInt, TInt>& edge_list = weights[i];
    for (THash<TInt, TInt>::TIter it = edge_list.BegI(); it < edge_list.EndI();
         it++) {    
      graph->AddEdge(i, it->Key);
    }
  }

  return graph;
}