68static double get_mq(
SparseMatrix A,
int *assignment,
int *ncluster0,
double *mq_in0,
double *mq_out0,
double **dout0){
78 bool test_pattern_symmetry_only =
false;
79 int *counts, *ia =
A->ia, *ja =
A->ja, k, i, j, jj;
80 double mq_in = 0, mq_out = 0, *a =
NULL, Vi, Vj;
86 (void)test_pattern_symmetry_only;
92 for (i = 0; i < n; i++){
93 assert(assignment[i] >= 0 && assignment[i] < n);
94 if (counts[assignment[i]] == 0) ncluster++;
95 counts[assignment[i]]++;
98 assert(ncluster <= n);
100 for (i = 0; i < n; i++){
101 assert(assignment[i] < ncluster);
104 for (j = ia[i] ; j < ia[i+1]; j++){
107 if (jj >= i)
continue;
108 assert(assignment[jj] < ncluster);
109 Vj = counts[assignment[jj]];
110 if (assignment[jj] == c){
112 mq_in += a[j]/(Vi*Vi);
118 mq_out += a[j]/(Vi*Vj);
120 mq_out += 1./(Vi*Vj);
129 for (i = 0; i < n; i++){
130 for (j = ia[i]; j < ia[i+1]; j++){
132 if (jj == i)
continue;
134 dout[i] += a[j]/(double) counts[assignment[jj]];
136 dout[i] += 1./(double) counts[assignment[jj]];
148 return 2*(mq_in/k - mq_out/(k*(k-1)));
228 int n =
grid->
n, level =
grid->level, nc = 0, nclusters = n;
229 double mq = 0, mq_in = 0, mq_out = 0, mq_new, mq_in_new, mq_out_new, mq_max = 0, mq_in_max = 0, mq_out_max = 0;
230 int *ia =
A->ia, *ja =
A->ja;
232 double *deg_intra =
grid->deg_intra, *wgt =
grid->wgt;
233 int i, j, k, jj, jc, jmax;
234 double gain = 0, *dout =
grid->dout, deg_in_i, deg_in_j, wgt_i, wgt_j, a_ij, dout_i, dout_j, dout_max = 0, wgt_jmax = 0;
236 double total_gain = 0;
242 mq_out =
grid->mq_out;
244 double *deg_intra_new =
gv_calloc(n,
sizeof(
double));
245 double *wgt_new =
gv_calloc(n,
sizeof(
double));
246 double *deg_inter =
gv_calloc(n,
sizeof(
double));
248 double *dout_new =
gv_calloc(n,
sizeof(
double));
249 for (i = 0; i < n; i++) mask[i] = -1;
252 for (i = 0; i < n; i++) matching[i] =
UNMATCHED;
279 for (i = 0; i < n; i++){
282 for (j = ia[i]; j < ia[i+1]; j++){
284 if (jj == i)
continue;
288 deg_inter[jc] = a[j];
290 deg_inter[jc] += a[j];
294 deg_in_i = deg_intra[i];
300 for (j = ia[i]; j < ia[i+1]; j++){
302 if (jj == i)
continue;
307 deg_in_j = deg_intra[jj];
309 }
else if (deg_inter[jc] < 0){
312 a_ij = deg_inter[jc];
315 deg_in_j = deg_intra_new[jc];
316 dout_j = dout_new[jc];
319 mq_in_new = mq_in - deg_in_i/pow(wgt_i, 2) - deg_in_j/pow(wgt_j,2)
320 + (deg_in_i + deg_in_j + a_ij)/pow(wgt_i + wgt_j,2);
322 mq_out_new = mq_out - dout_i/wgt_i - dout_j/wgt_j + (dout_i + dout_j)/(wgt_i + wgt_j);
325 mq_new = 2*(mq_in_new/(nclusters - 1) - mq_out_new/((nclusters - 1)*(nclusters - 2)));
327 mq_new = 2*mq_in_new/(nclusters - 1);
332 double mq2, mq_in2, mq_out2, *dout2;
334 int *matching2 =
gv_calloc(
A->m,
sizeof(
int));
335 memcpy(matching2, matching,
sizeof(
double)*
A->m);
343 for (k = 0; k < n; k++)
if (matching2[k] ==
UNMATCHED) matching2[k] =nc2++;
344 mq2 =
get_mq(
A, matching2, &ncluster, &mq_in2, &mq_out2, &dout2);
345 fprintf(stderr,
" {dout_i, dout_j}={%f,%f}, {predicted, calculated}: mq = {%f, %f}, mq_in ={%f,%f}, mq_out = {%f,%f}\n",dout_i, dout_j, mq_new, mq2, mq_in_new, mq_in2, mq_out_new, mq_out2);
353 if (
Verbose) fprintf(stderr,
"gain in merging node %d with node %d = %f-%f = %f\n", i, jj, mq, mq_new, gain);
354 if (j == ia[i] || gain > maxgain){
361 mq_in_max = mq_in_new;
362 mq_out_max = mq_out_new;
368 if (maxgain > 0 || (nc >= 1 && nc > maxcluster)){
369 total_gain += maxgain;
372 fprintf(stderr,
"maxgain=%f, merge %d, %d\n",maxgain, i, jmax);
375 dout_new[nc] = dout_i + dout_max;
376 matching[i] = matching[jmax] = nc;
377 wgt_new[nc] = wgt[i] + wgt[jmax];
378 deg_intra_new[nc] = deg_intra[i] + deg_intra[jmax] + amax;
381 fprintf(stderr,
"maxgain=%f, merge with existing cluster %d, %d\n",maxgain, i, jc);
383 dout_new[jc] = dout_i + dout_max;
384 wgt_new[jc] += wgt[i];
386 deg_intra_new[jc] += deg_intra[i] + amax;
393 fprintf(stderr,
"gain: %f -- no gain, skip merging node %d\n", maxgain, i);
394 assert(maxgain <= 0);
397 deg_intra_new[nc] = deg_intra[i];
398 wgt_new[nc] = wgt[i];
406 mask[
LIST_GET(&neighbors[jc], l)] = n + i;
411 for (j = ia[k]; j < ia[k+1]; j++){
413 if (mask[jj] == n+i)
continue;
416 dout[jj] += -a[j]/wgt_i + a[j]/(wgt_i + wgt_jmax);
418 dout[jj] += -a[j]/wgt_jmax + a[j]/(wgt_i + wgt_jmax);
422 dout_new[jc] += -a[j]/wgt_i + a[j]/(wgt_i + wgt_jmax);
424 dout_new[jc] += -a[j]/wgt_jmax + a[j]/(wgt_i + wgt_jmax);
432 fprintf(stderr,
"verbose=%d\n",
Verbose);
433 if (
Verbose) fprintf(stderr,
"mq = %f new mq = %f level = %d, n = %d, nc = %d, gain = %g, mq_in = %f, mq_out = %f\n", mq, mq + total_gain,
434 level, n, nc, total_gain, mq_in, mq_out);
439 mq =
get_mq(
A, matching, &ncluster, &mq_in, &mq_out, &dout);
440 fprintf(stderr,
" mq = %f\n",mq);
445 if (nc >= 1 && (total_gain > 0 || nc < n)){
452 for (i = 0; i < n; i++){
478 deg_intra_new =
gv_recalloc(deg_intra_new, n, nc,
sizeof(
double));
479 wgt_new =
gv_recalloc(wgt_new, n, nc,
sizeof(
double));
481 cgrid->
mq =
grid->mq + total_gain;
482 cgrid->
wgt = wgt_new;
483 dout_new =
gv_recalloc(dout_new, n, nc,
sizeof(
double));
484 cgrid->
dout = dout_new;
492 for (i = 0; i < n; i++) matching[i] = i;
500 for (i = 0; i < n; i++)
LIST_FREE(&neighbors[i]);