SNAP Library 6.0, User Reference  2020-12-09 16:24:20
SNAP, a general purpose, high performance system for analysis and manipulation of large networks
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.

29  {
30  int minval = MIN(i, j);
31  int maxval = MAX(i, j);
32  THash<TInt, TInt>& edge_weights = weights[minval];
33  edge_weights(maxval) += 1;
34 }
#define MIN(a, b)
Definition: bd.h:346
#define MAX(a, b)
Definition: bd.h:350
static void MapIdsToFirstN ( const TIntV ids,
THash< TInt, TInt > &  id_map,
TIntV rev_id_map 
)
static

Definition at line 796 of file motifcluster.cpp.

797  {
798  id_map = THash<TInt, TInt>();
799  rev_id_map = TIntV(ids.Len());
800  for (int i = 0; i < ids.Len(); i++) {
801  int id = ids[i];
802  if (id_map.IsKey(id)) {
803  TExcept::Throw("List of ids is not unique.");
804  }
805  id_map(id) = i;
806  rev_id_map[i] = id;
807  }
808 }
TSizeTy Len() const
Returns the number of elements in the vector.
Definition: ds.h:575
static void Throw(const TStr &MsgStr)
Definition: ut.h:187
TVec< TInt > TIntV
Definition: ds.h:1594
bool IsKey(const TKey &Key) const
Definition: hash.h:258
static void Sweep ( const TSparseColMatrix W,
const TFltV fvec,
TFltV conds,
TIntV order 
)
static

Definition at line 813 of file motifcluster.cpp.

814  {
815  // Get ordering of nodes
816 
817  TVec< TKeyDat<TFlt, TInt> > fvec_inds(fvec.Len());
818  for (int i = 0; i < fvec.Len(); i++) {
819  fvec_inds[i] = TKeyDat<TFlt, TInt>(fvec[i], i);
820  }
821  fvec_inds.Sort();
822  order = TIntV(fvec.Len());
823  TIntV rank(fvec.Len());
824  for (int i = 0; i < fvec.Len(); i++) {
825  order[i] = fvec_inds[i].Dat;
826  rank[order[i]] = i;
827  }
828 
829  // Sweep by adjusting cut and volume
830  conds = TFltV(order.Len() - 1);
831  double cut = 0;
832  double vol = 0;
833  double total_vol = 0;
834  for (int ind = 0; ind < order.Len(); ind++) {
835  const TIntFltKdV& nbr_weights = W.ColSpVV[ind];
836  for (TIntFltKdV::TIter it = nbr_weights.BegI(); it < nbr_weights.EndI();
837  it++) {
838  total_vol += it->Dat;
839  }
840  }
841  double vol_comp = total_vol;
842 
843  for (int ind = 0; ind < order.Len() - 1; ind++) {
844  int node = order[ind];
845  const TIntFltKdV& nbr_weights = W.ColSpVV[node];
846  for (TIntFltKdV::TIter it = nbr_weights.BegI(); it < nbr_weights.EndI();
847  it++) {
848  int nbr = it->Key;
849  if (node == nbr) { continue; }
850  double val = it->Dat;
851  // Adjust the cut amount
852  if (rank[nbr] > ind) {
853  // nbr is on the other side: add to the cut
854  cut += val;
855  } else {
856  // now on the same side as nbr: subtract from the cut
857  cut -= val;
858  }
859  vol += val;
860  vol_comp -= val;
861  }
862 
863  double mvol = MIN(vol, vol_comp);
864  if (mvol <= 0.0) {
865  TExcept::Throw("Nonpositive set volume.");
866  }
867  conds[ind] = cut / mvol;
868  }
869 }
Definition: ds.h:346
TIter EndI() const
Returns an iterator referring to the past-the-end element in the vector.
Definition: ds.h:595
TSizeTy Len() const
Returns the number of elements in the vector.
Definition: ds.h:575
static void Throw(const TStr &MsgStr)
Definition: ut.h:187
#define MIN(a, b)
Definition: bd.h:346
TVec< TIntFltKdV > ColSpVV
Definition: linalg.h:60
TVal * TIter
Random access iterator to TVal.
Definition: ds.h:432
TVec< TFlt > TFltV
Definition: ds.h:1596
TIter BegI() const
Returns an iterator pointing to the first element in the vector.
Definition: ds.h:593
TVec< TInt > TIntV
Definition: ds.h:1594
Vector is a sequence TVal objects representing an array that can change in size.
Definition: ds.h:430
void SymeigsSmallest ( const TSparseColMatrix A,
int  nev,
TFltV evals,
TFullColMatrix evecs,
double  tol,
int  maxiter 
)

Definition at line 955 of file motifcluster.cpp.

956  {
957  // type of problem
958  int mode = 1;
959  // communication variable
960  int ido = 0;
961  // standard eigenvalue problem
962  char bmat = 'I';
963  // dimension of problem
964  int n = A.GetCols();
965  // type of eigenvector (smallest algebraic)
966  char which[] = "SA";
967  // Space for residual
968  double *resid = new double[n];
969  // Number of lanczos vectors
970  int ncv = MIN(MAX(2 * nev, 20), n - 1);
971  // Space for Lanczos basis vectors
972  double *V = new double[ncv * n];
973 
974  // dimension of basis vectors
975  int ldv = n;
976  // Various input / output data
977  int *iparam = new int[11];
978  // exact shifts
979  iparam[0] = 1;
980  // maximum number of iterations
981  iparam[2] = maxiter;
982  // mode
983  iparam[6] = mode;
984  // Output data
985  int *ipntr = new int[11];
986  // work array
987  double *workd = new double[3 * n];
988  // output workspace
989  int lworkl = ncv * (ncv + 8);
990  double *workl = new double[lworkl];
991  // Communication information
992  int info = 0;
993 
994  TFltV Ax(n);
995  // Communication loop. We keep applying A * x until ARPACK tells us to stop.
996  while (true) {
997  F77_NAME(dsaupd)(&ido, &bmat, &n, &which[0], &nev, &tol, &resid[0], &ncv,
998  &V[0], &ldv, &iparam[0], &ipntr[0], &workd[0], &workl[0],
999  &lworkl, &info);
1000  TFltV result(n);
1001  result.PutAll(0.0);
1002  double *load = &workd[ipntr[0] - 1];
1003  for (int i = 0; i < n; i++) {
1004  result[i] = load[i];
1005  }
1006 
1007  if (ido == 1) {
1008  // Need another matrix-vector product.
1009  A.Multiply(result, Ax);
1010  double *store = &workd[ipntr[1] - 1];
1011  for (int i = 0; i < n; i++) {
1012  store[i] = Ax[i];
1013  }
1014  } else if (ido == 99) {
1015  // We have finished.
1016  break;
1017  } else {
1018  // There was an error.
1019  TExcept::Throw("ARPACK error");
1020  }
1021  }
1022 
1023  // Extract the eigenvectors and eigenvalues
1024 
1025  // compute Ritz vectors
1026  int rvec(true);
1027  // Form of Ritz vectors
1028  char howmny = 'A';
1029  // select doesn't matter with howmny set to A
1030  int *select = new int[ncv];
1031  // Store Ritz values
1032  double *d = new double[nev];
1033  // Shift
1034  double sigmar = 0;
1035 
1036  F77_NAME(dseupd)(&rvec, &howmny, select, &d[0], &V[0], &ldv, &sigmar, &bmat,
1037  &n, &which[0], &nev, &tol, &resid[0], &ncv, &V[0], &ldv,
1038  &iparam[0], &ipntr[0], &workd[0], &workl[0], &lworkl, &info);
1039  if (info != 0) {
1040  TExcept::Throw("ARPACK error");
1041  }
1042 
1043  // Get eigenvalues and eigenvalues in sorted order
1044  TVec< TKeyDat<TFlt, TInt> > d_inds(nev);
1045  for (int i = 0; i < nev; i++) {
1046  d_inds[i] = TKeyDat<TFlt, TInt>(d[i], i);
1047  }
1048  d_inds.Sort();
1049 
1050  evals = TFltV(nev);
1051  evecs.ColV = TVec<TFltV>(nev);
1052  for (int i = 0; i < nev; i++) {
1053  double eval = d_inds[i].Key;
1054  int ind = d_inds[i].Dat;
1055  evals[ind] = eval;
1056  TFltV& col = evecs.ColV[ind];
1057  for (int j = ind * n; j < ind * n + n; j++) {
1058  col.Add(V[j]);
1059  }
1060  }
1061  evecs.RowN = n;
1062  evecs.ColN = nev;
1063 
1064  // Clean up
1065  delete[] resid;
1066  delete[] V;
1067  delete[] iparam;
1068  delete[] ipntr;
1069  delete[] workd;
1070  delete[] workl;
1071  delete[] select;
1072  delete[] d;
1073 }
Definition: ds.h:346
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)
static void Throw(const TStr &MsgStr)
Definition: ut.h:187
#define MIN(a, b)
Definition: bd.h:346
int GetCols() const
Definition: linalg.h:47
void Multiply(const TFltVV &B, int ColId, TFltV &Result) const
Definition: linalg.h:24
TVec< TFltV > ColV
Definition: linalg.h:131
TVec< TFlt > TFltV
Definition: ds.h:1596
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)
#define F77_NAME(name)
Definition: motifcluster.cpp:7
TSizeTy Add()
Adds a new element at the end of the vector, after its current last element.
Definition: ds.h:602
#define MAX(a, b)
Definition: bd.h:350
static PUNGraph UnweightedGraphRepresentation ( const WeightVH weights)
static

Definition at line 772 of file motifcluster.cpp.

772  {
773  int num_edges = 0;
774  for (int i = 0; i < weights.Len(); i++) {
775  num_edges += weights[i].Len();
776  }
777  PUNGraph graph = TUNGraph::New(weights.Len(), num_edges);
778  for (int i = 0; i < weights.Len(); i++) {
779  graph->AddNode(i);
780  }
781  for (int i = 0; i < weights.Len(); i++) {
782  const THash<TInt, TInt>& edge_list = weights[i];
783  for (THash<TInt, TInt>::TIter it = edge_list.BegI(); it < edge_list.EndI();
784  it++) {
785  graph->AddEdge(i, it->Key);
786  }
787  }
788 
789  return graph;
790 }
TIter BegI() const
Definition: hash.h:213
TSizeTy Len() const
Returns the number of elements in the vector.
Definition: ds.h:575
TIter EndI() const
Definition: hash.h:218
static PUNGraph New()
Static constructor that returns a pointer to the graph. Call: PUNGraph Graph = TUNGraph::New().
Definition: graph.h:172
Definition: bd.h:196