← Back to branch summary

~ubuntu-branches/ubuntu/lucid/igraph/lucid

« back to all changes in this revision

Viewing changes to src/walktrap_communities.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Mathieu Malaterre
  • Date: 2009-11-16 18:12:42 UTC
  • Revision ID: james.westby@ubuntu.com-20091116181242-mzv9p5fz9uj57xd1
Tags: upstream-0.5.3
Import upstream version 0.5.3

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/walktrap_communities.cpp
 
1
/* -*- mode: C -*-  */
 
2
/* 
 
3
   IGraph library.
 
4
   Copyright (C) 2007  Gabor Csardi <csardi@rmki.kfki.hu>
 
5
   MTA RMKI, Konkoly-Thege Miklos st. 29-33, Budapest 1121, Hungary
 
6
   
 
7
   This program is free software; you can redistribute it and/or modify
 
8
   it under the terms of the GNU General Public License as published by
 
9
   the Free Software Foundation; either version 2 of the License, or
 
10
   (at your option) any later version.
 
11
   
 
12
   This program is distributed in the hope that it will be useful,
 
13
   but WITHOUT ANY WARRANTY; without even the implied warranty of
 
14
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
15
   GNU General Public License for more details.
 
16
   
 
17
   You should have received a copy of the GNU General Public License
 
18
   along with this program; if not, write to the Free Software
 
19
   Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 
 
20
   02110-1301 USA
 
21
 
 
22
*/
 
23
 
 
24
/* The original version of this file was written by Pascal Pons
 
25
   The original copyright notice follows here. The FSF address was
 
26
   fixed by Tamas Nepusz */
 
27
 
 
28
// File: communities.cpp
 
29
//-----------------------------------------------------------------------------
 
30
// Walktrap v0.2 -- Finds community structure of networks using random walks
 
31
// Copyright (C) 2004-2005 Pascal Pons
 
32
//
 
33
// This program is free software; you can redistribute it and/or modify
 
34
// it under the terms of the GNU General Public License as published by
 
35
// the Free Software Foundation; either version 2 of the License, or
 
36
// (at your option) any later version.
 
37
//
 
38
// This program is distributed in the hope that it will be useful,
 
39
// but WITHOUT ANY WARRANTY; without even the implied warranty of
 
40
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
41
// GNU General Public License for more details.
 
42
//
 
43
// You should have received a copy of the GNU General Public License
 
44
// along with this program; if not, write to the Free Software
 
45
// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 
 
46
// 02110-1301 USA
 
47
//-----------------------------------------------------------------------------
 
48
// Author   : Pascal Pons 
 
49
// Email    : pons@liafa.jussieu.fr
 
50
// Web page : http://www.liafa.jussieu.fr/~pons/
 
51
// Location : Paris, France
 
52
// Time     : June 2005
 
53
//-----------------------------------------------------------------------------
 
54
// see readme.txt for more details
 
55
 
 
56
#include "walktrap_communities.h"
 
57
#include <cstdlib>
 
58
#include <iostream>
 
59
#include <cmath>
 
60
#include <algorithm>
 
61
 
 
62
int Probabilities::length = 0;
 
63
Communities* Probabilities::C = 0;
 
64
float* Probabilities::tmp_vector1 = 0;
 
65
float* Probabilities::tmp_vector2 = 0;
 
66
int* Probabilities::id = 0;
 
67
int* Probabilities::vertices1 = 0;
 
68
int* Probabilities::vertices2 = 0;
 
69
int Probabilities::current_id = 0;
 
70
 
 
71
 
 
72
Neighbor::Neighbor() {
 
73
  next_community1 = 0;
 
74
  previous_community1 = 0;
 
75
  next_community2 = 0;       
 
76
  previous_community2 = 0;
 
77
  heap_index = -1;
 
78
}
 
79
 
 
80
Probabilities::~Probabilities() {
 
81
  C->memory_used -= memory();
 
82
  if(P) delete[] P;
 
83
  if(vertices) delete[] vertices;
 
84
}
 
85
 
 
86
Probabilities::Probabilities(int community) {
 
87
  Graph* G = C->G;
 
88
  int nb_vertices1 = 0;
 
89
  int nb_vertices2 = 0;
 
90
 
 
91
  float initial_proba = 1./float(C->communities[community].size);
 
92
  int last =  C->members[C->communities[community].last_member];  
 
93
  for(int m = C->communities[community].first_member; m != last; m = C->members[m]) {
 
94
    tmp_vector1[m] = initial_proba;
 
95
    vertices1[nb_vertices1++] = m;
 
96
  }
 
97
  
 
98
  for(int t = 0; t < length; t++) {
 
99
    current_id++;    
 
100
    if(nb_vertices1 > (G->nb_vertices/2)) {
 
101
      nb_vertices2 = G->nb_vertices;
 
102
      for(int i = 0; i < G->nb_vertices; i++)
 
103
        tmp_vector2[i] = 0.;
 
104
      if(nb_vertices1 == G->nb_vertices) {
 
105
        for(int i = 0; i < G->nb_vertices; i++) {
 
106
          float proba = tmp_vector1[i]/G->vertices[i].total_weight;
 
107
          for(int j = 0; j < G->vertices[i].degree; j++)
 
108
            tmp_vector2[G->vertices[i].edges[j].neighbor] += proba*G->vertices[i].edges[j].weight;
 
109
        }
 
110
      }
 
111
      else {
 
112
        for(int i = 0; i < nb_vertices1; i++) {
 
113
          int v1 = vertices1[i];
 
114
          float proba = tmp_vector1[v1]/G->vertices[v1].total_weight;
 
115
          for(int j = 0; j < G->vertices[v1].degree; j++)
 
116
            tmp_vector2[G->vertices[v1].edges[j].neighbor] += proba*G->vertices[v1].edges[j].weight;
 
117
        }
 
118
      }
 
119
    }
 
120
    else {
 
121
      nb_vertices2 = 0;
 
122
      for(int i = 0; i < nb_vertices1; i++) {
 
123
        int v1 = vertices1[i];
 
124
        float proba = tmp_vector1[v1]/G->vertices[v1].total_weight;
 
125
        for(int j = 0; j < G->vertices[v1].degree; j++) {
 
126
          int v2 = G->vertices[v1].edges[j].neighbor;
 
127
          if(id[v2] == current_id)
 
128
            tmp_vector2[v2] += proba*G->vertices[v1].edges[j].weight;
 
129
          else {
 
130
            tmp_vector2[v2] = proba*G->vertices[v1].edges[j].weight;
 
131
            id[v2] = current_id;
 
132
            vertices2[nb_vertices2++] = v2;
 
133
          }
 
134
        }
 
135
      }
 
136
    }
 
137
    float* tmp = tmp_vector2;
 
138
    tmp_vector2 = tmp_vector1;
 
139
    tmp_vector1 = tmp;
 
140
 
 
141
    int* tmp2 = vertices2;
 
142
    vertices2 = vertices1;
 
143
    vertices1 = tmp2;
 
144
 
 
145
    nb_vertices1 = nb_vertices2;
 
146
  }
 
147
 
 
148
  if(nb_vertices1 > (G->nb_vertices/2)) {
 
149
    P = new float[G->nb_vertices];
 
150
    size = G->nb_vertices;
 
151
    vertices = 0;
 
152
    if(nb_vertices1 == G->nb_vertices) {
 
153
      for(int i = 0; i < G->nb_vertices; i++)
 
154
        P[i] = tmp_vector1[i]/sqrt(G->vertices[i].total_weight);
 
155
    }
 
156
    else {
 
157
      for(int i = 0; i < G->nb_vertices; i++)
 
158
        P[i] = 0.;
 
159
      for(int i = 0; i < nb_vertices1; i++)
 
160
        P[vertices1[i]] = tmp_vector1[vertices1[i]]/sqrt(G->vertices[vertices1[i]].total_weight);
 
161
    }
 
162
  }
 
163
  else {
 
164
    P = new float[nb_vertices1];
 
165
    size = nb_vertices1;
 
166
    vertices = new int[nb_vertices1];
 
167
    int j = 0;
 
168
    for(int i = 0; i < G->nb_vertices; i++) {
 
169
      if(id[i] == current_id) {
 
170
        P[j] = tmp_vector1[i]/sqrt(G->vertices[i].total_weight);
 
171
        vertices[j] = i;
 
172
        j++;
 
173
      }
 
174
    }
 
175
  }
 
176
  C->memory_used += memory();
 
177
}
 
178
 
 
179
Probabilities::Probabilities(int community1, int community2) {
 
180
  // The two following probability vectors must exist.
 
181
  // Do not call this function if it is not the case.
 
182
  Probabilities* P1 = C->communities[community1].P;
 
183
  Probabilities* P2 = C->communities[community2].P;
 
184
  
 
185
  float w1 = float(C->communities[community1].size)/float(C->communities[community1].size + C->communities[community2].size);
 
186
  float w2 = float(C->communities[community2].size)/float(C->communities[community1].size + C->communities[community2].size);
 
187
 
 
188
 
 
189
  if(P1->size == C->G->nb_vertices) {
 
190
    P = new float[C->G->nb_vertices];
 
191
    size = C->G->nb_vertices;
 
192
    vertices = 0;
 
193
    
 
194
    if(P2->size == C->G->nb_vertices) { // two full vectors
 
195
      for(int i = 0; i < C->G->nb_vertices; i++)
 
196
        P[i] = P1->P[i]*w1 + P2->P[i]*w2;
 
197
    }
 
198
    else {  // P1 full vector, P2 partial vector
 
199
      int j = 0;
 
200
      for(int i = 0; i < P2->size; i++) {
 
201
        for(; j < P2->vertices[i]; j++)
 
202
          P[j] = P1->P[j]*w1;
 
203
        P[j] = P1->P[j]*w1 + P2->P[i]*w2;
 
204
        j++;
 
205
      }
 
206
      for(; j < C->G->nb_vertices; j++)
 
207
        P[j] = P1->P[j]*w1;
 
208
    }
 
209
  }
 
210
  else {
 
211
    if(P2->size == C->G->nb_vertices) { // P1 partial vector, P2 full vector
 
212
      P = new float[C->G->nb_vertices];
 
213
      size = C->G->nb_vertices;
 
214
      vertices = 0;
 
215
 
 
216
      int j = 0;
 
217
      for(int i = 0; i < P1->size; i++) {
 
218
        for(; j < P1->vertices[i]; j++)
 
219
          P[j] = P2->P[j]*w2;
 
220
        P[j] = P1->P[i]*w1 + P2->P[j]*w2;
 
221
        j++;
 
222
      }
 
223
      for(; j < C->G->nb_vertices; j++)
 
224
        P[j] = P2->P[j]*w2;
 
225
    }
 
226
    else {  // two partial vectors
 
227
      int i = 0;
 
228
      int j = 0;
 
229
      int nb_vertices1 = 0;
 
230
      while((i < P1->size) && (j < P2->size)) {
 
231
        if(P1->vertices[i] < P2->vertices[j]) {
 
232
          tmp_vector1[P1->vertices[i]] = P1->P[i]*w1;
 
233
          vertices1[nb_vertices1++] = P1->vertices[i];
 
234
          i++;
 
235
          continue;
 
236
        }
 
237
        if(P1->vertices[i] > P2->vertices[j]) {
 
238
          tmp_vector1[P2->vertices[j]] = P2->P[j]*w2;
 
239
          vertices1[nb_vertices1++] = P2->vertices[j];
 
240
          j++;
 
241
          continue;
 
242
        }
 
243
        tmp_vector1[P1->vertices[i]] = P1->P[i]*w1 + P2->P[j]*w2;
 
244
        vertices1[nb_vertices1++] = P1->vertices[i];
 
245
        i++;
 
246
        j++;
 
247
      }
 
248
      if(i == P1->size) {
 
249
        for(; j < P2->size; j++) {
 
250
          tmp_vector1[P2->vertices[j]] = P2->P[j]*w2;
 
251
          vertices1[nb_vertices1++] = P2->vertices[j];
 
252
        }
 
253
      }
 
254
      else {
 
255
        for(; i < P1->size; i++) {
 
256
          tmp_vector1[P1->vertices[i]] = P1->P[i]*w1;
 
257
          vertices1[nb_vertices1++] = P1->vertices[i];
 
258
        }
 
259
      }
 
260
 
 
261
      if(nb_vertices1 > (C->G->nb_vertices/2)) {
 
262
        P = new float[C->G->nb_vertices];
 
263
        size = C->G->nb_vertices;
 
264
        vertices = 0;
 
265
        for(int i = 0; i < C->G->nb_vertices; i++)
 
266
          P[i] = 0.;
 
267
        for(int i = 0; i < nb_vertices1; i++)
 
268
          P[vertices1[i]] = tmp_vector1[vertices1[i]];
 
269
      }
 
270
      else {
 
271
        P = new float[nb_vertices1];
 
272
        size = nb_vertices1;
 
273
        vertices = new int[nb_vertices1];
 
274
        for(int i = 0; i < nb_vertices1; i++) {
 
275
          vertices[i] = vertices1[i];
 
276
          P[i] = tmp_vector1[vertices1[i]];
 
277
        }
 
278
      }
 
279
    }
 
280
  }
 
281
 
 
282
  C->memory_used += memory();
 
283
}
 
284
 
 
285
double Probabilities::compute_distance(const Probabilities* P2) const {
 
286
  double r = 0.;
 
287
  if(vertices) {
 
288
    if(P2->vertices) {  // two partial vectors
 
289
      int i = 0;
 
290
      int j = 0;
 
291
      while((i < size) && (j < P2->size)) {
 
292
        if(vertices[i] < P2->vertices[j]) {
 
293
          r += P[i]*P[i];
 
294
          i++;
 
295
          continue;
 
296
        }
 
297
        if(vertices[i] > P2->vertices[j]) {
 
298
          r += P2->P[j]*P2->P[j];
 
299
          j++;
 
300
          continue;
 
301
        }
 
302
        r += (P[i] - P2->P[j])*(P[i] - P2->P[j]);
 
303
        i++;
 
304
        j++;
 
305
      }
 
306
      if(i == size) {
 
307
        for(; j < P2->size; j++)
 
308
          r += P2->P[j]*P2->P[j];
 
309
      }
 
310
      else {
 
311
        for(; i < size; i++)
 
312
          r += P[i]*P[i];
 
313
      }
 
314
    }
 
315
    else {  // P1 partial vector, P2 full vector 
 
316
 
 
317
      int i = 0;
 
318
      for(int j = 0; j < size; j++) {
 
319
        for(; i < vertices[j]; i++)
 
320
          r += P2->P[i]*P2->P[i];
 
321
        r += (P[j] - P2->P[i])*(P[j] - P2->P[i]);
 
322
        i++;
 
323
      }
 
324
      for(; i < P2->size; i++)
 
325
        r += P2->P[i]*P2->P[i];      
 
326
    }
 
327
  }
 
328
  else {
 
329
    if(P2->vertices) {  // P1 full vector, P2 partial vector
 
330
      int i = 0;
 
331
      for(int j = 0; j < P2->size; j++) {
 
332
        for(; i < P2->vertices[j]; i++)
 
333
          r += P[i]*P[i];
 
334
        r += (P[i] - P2->P[j])*(P[i] - P2->P[j]);
 
335
        i++;
 
336
      }
 
337
      for(; i < size; i++)
 
338
        r += P[i]*P[i];
 
339
    }
 
340
    else {  // two full vectors
 
341
      for(int i = 0; i < size; i++)
 
342
        r += (P[i] - P2->P[i])*(P[i] - P2->P[i]);
 
343
    }
 
344
  }
 
345
  return r;
 
346
}
 
347
 
 
348
long Probabilities::memory() {
 
349
  if(vertices)
 
350
    return (sizeof(Probabilities) + long(size)*(sizeof(float) + sizeof(int)));
 
351
  else
 
352
    return (sizeof(Probabilities) + long(size)*sizeof(float));
 
353
}
 
354
 
 
355
Community::Community() {
 
356
  P = 0;
 
357
  first_neighbor = 0;
 
358
  last_neighbor = 0;
 
359
  sub_community_of = -1;
 
360
  sub_communities[0] = -1;
 
361
  sub_communities[1] = -1;
 
362
  sigma = 0.;
 
363
  internal_weight = 0.;
 
364
  total_weight = 0.;
 
365
}
 
366
 
 
367
Community::~Community() {
 
368
  if(P) delete P;
 
369
}
 
370
 
 
371
 
 
372
Communities::Communities(Graph* graph, int random_walks_length, 
 
373
                         long m, igraph_matrix_t *pmerges,
 
374
                         igraph_vector_t *pmodularity) {
 
375
  max_memory = m;
 
376
  memory_used = 0;
 
377
  G = graph;
 
378
  merges=pmerges;
 
379
  mergeidx=0;
 
380
  modularity=pmodularity;
 
381
  
 
382
  Probabilities::C = this;
 
383
  Probabilities::length = random_walks_length;
 
384
  Probabilities::tmp_vector1 = new float[G->nb_vertices];
 
385
  Probabilities::tmp_vector2 = new float[G->nb_vertices];
 
386
  Probabilities::id = new int[G->nb_vertices];
 
387
  for(int i = 0; i < G->nb_vertices; i++) Probabilities::id[i] = 0;
 
388
  Probabilities::vertices1 = new int[G->nb_vertices];
 
389
  Probabilities::vertices2 = new int[G->nb_vertices];
 
390
  Probabilities::current_id = 0;
 
391
 
 
392
  
 
393
  members = new int[G->nb_vertices];  
 
394
  for(int i = 0; i < G->nb_vertices; i++)
 
395
    members[i] = -1;
 
396
 
 
397
  H = new Neighbor_heap(G->nb_edges);
 
398
  communities = new Community[2*G->nb_vertices];
 
399
 
 
400
// init the n single vertex communities
 
401
 
 
402
  if(max_memory != -1)
 
403
    min_delta_sigma = new Min_delta_sigma_heap(G->nb_vertices*2);
 
404
  else min_delta_sigma = 0;
 
405
  
 
406
  for(int i = 0; i < G->nb_vertices; i++) {
 
407
    communities[i].this_community = i;
 
408
    communities[i].first_member = i;
 
409
    communities[i].last_member = i;
 
410
    communities[i].size = 1;
 
411
    communities[i].sub_community_of = 0;
 
412
  }
 
413
 
 
414
  nb_communities = G->nb_vertices;
 
415
  nb_active_communities = G->nb_vertices;
 
416
 
 
417
  for(int i = 0; i < G->nb_vertices; i++)
 
418
    for(int j = 0; j < G->vertices[i].degree; j++)
 
419
      if (i < G->vertices[i].edges[j].neighbor) {
 
420
        communities[i].total_weight += G->vertices[i].edges[j].weight/2.;
 
421
        communities[G->vertices[i].edges[j].neighbor].total_weight += G->vertices[i].edges[j].weight/2.;
 
422
        Neighbor* N = new Neighbor;
 
423
        N->community1 = i;
 
424
        N->community2 = G->vertices[i].edges[j].neighbor;
 
425
        N->delta_sigma = -1./double(min(G->vertices[i].degree,  G->vertices[G->vertices[i].edges[j].neighbor].degree));
 
426
        N->weight = G->vertices[i].edges[j].weight;
 
427
        N->exact = false;
 
428
        add_neighbor(N);
 
429
      }
 
430
 
 
431
  if(max_memory != -1) {
 
432
    memory_used += min_delta_sigma->memory();
 
433
    memory_used += 2*long(G->nb_vertices)*sizeof(Community);
 
434
    memory_used += long(G->nb_vertices)*(2*sizeof(float) + 3*sizeof(int)); // the static data of Probabilities class
 
435
    memory_used += H->memory() + long(G->nb_edges)*sizeof(Neighbor);
 
436
    memory_used += G->memory();    
 
437
  }
 
438
 
 
439
/*   int c = 0; */
 
440
  Neighbor* N = H->get_first();  
 
441
  while(!N->exact) {
 
442
    update_neighbor(N, compute_delta_sigma(N->community1, N->community2));
 
443
    N->exact = true;
 
444
    N = H->get_first();
 
445
    if(max_memory != -1) manage_memory();
 
446
    /* TODO: this could use igraph_progress */
 
447
/*     if(!silent) { */
 
448
/*       c++; */
 
449
/*       for(int k = (500*(c-1))/G->nb_edges + 1; k <= (500*c)/G->nb_edges; k++) { */
 
450
/*      if(k % 50 == 1) {cerr.width(2); cerr << endl << k/ 5 << "% ";} */
 
451
/*      cerr << "."; */
 
452
/*       } */
 
453
/*     } */
 
454
  }
 
455
  
 
456
}
 
457
 
 
458
Communities::~Communities() {
 
459
  delete[] members;
 
460
  delete[] communities;
 
461
  delete H;
 
462
  if(min_delta_sigma) delete min_delta_sigma;
 
463
  
 
464
  delete[] Probabilities::tmp_vector1;
 
465
  delete[] Probabilities::tmp_vector2;
 
466
  delete[] Probabilities::id;
 
467
  delete[] Probabilities::vertices1;
 
468
  delete[] Probabilities::vertices2;
 
469
}
 
470
 
 
471
float Community::min_delta_sigma() {
 
472
  float r = 1.;
 
473
  for(Neighbor* N = first_neighbor; N != 0;) {
 
474
    if(N->delta_sigma < r) r = N->delta_sigma;
 
475
    if(N->community1 == this_community)
 
476
      N = N->next_community1;
 
477
    else
 
478
      N = N->next_community2;
 
479
  }
 
480
  return r;
 
481
}
 
482
 
 
483
 
 
484
void Community::add_neighbor(Neighbor* N) { // add a new neighbor at the end of the list
 
485
  if (last_neighbor) {
 
486
    if(last_neighbor->community1 == this_community)
 
487
      last_neighbor->next_community1 = N;
 
488
    else
 
489
      last_neighbor->next_community2 = N;
 
490
    
 
491
    if(N->community1 == this_community)
 
492
      N->previous_community1 = last_neighbor;
 
493
    else
 
494
      N->previous_community2 = last_neighbor;
 
495
  }
 
496
  else {
 
497
    first_neighbor = N;
 
498
    if(N->community1 == this_community)
 
499
      N->previous_community1 = 0;
 
500
    else
 
501
      N->previous_community2 = 0;
 
502
  }
 
503
  last_neighbor = N;
 
504
}
 
505
 
 
506
void Community::remove_neighbor(Neighbor* N) {  // remove a neighbor from the list
 
507
  if (N->community1 == this_community) {
 
508
    if(N->next_community1) {
 
509
//      if (N->next_community1->community1 == this_community)
 
510
        N->next_community1->previous_community1 = N->previous_community1;
 
511
//      else 
 
512
//      N->next_community1->previous_community2 = N->previous_community1;
 
513
    }
 
514
    else last_neighbor = N->previous_community1;
 
515
    if(N->previous_community1) {
 
516
      if (N->previous_community1->community1 == this_community)
 
517
        N->previous_community1->next_community1 = N->next_community1;
 
518
      else 
 
519
        N->previous_community1->next_community2 = N->next_community1;
 
520
    }
 
521
    else first_neighbor = N->next_community1;
 
522
  }
 
523
  else {
 
524
    if(N->next_community2) {
 
525
      if (N->next_community2->community1 == this_community)
 
526
        N->next_community2->previous_community1 = N->previous_community2;
 
527
      else 
 
528
        N->next_community2->previous_community2 = N->previous_community2;
 
529
    }
 
530
    else last_neighbor = N->previous_community2;
 
531
    if(N->previous_community2) {
 
532
//      if (N->previous_community2->community1 == this_community)
 
533
//      N->previous_community2->next_community1 = N->next_community2;
 
534
//      else 
 
535
        N->previous_community2->next_community2 = N->next_community2;
 
536
    }
 
537
    else first_neighbor = N->next_community2;
 
538
  }
 
539
}
 
540
 
 
541
void Communities::remove_neighbor(Neighbor* N) {
 
542
  communities[N->community1].remove_neighbor(N);
 
543
  communities[N->community2].remove_neighbor(N);
 
544
  H->remove(N);
 
545
 
 
546
  if(max_memory !=-1) {
 
547
    if(N->delta_sigma == min_delta_sigma->delta_sigma[N->community1]) {
 
548
      min_delta_sigma->delta_sigma[N->community1] = communities[N->community1].min_delta_sigma();
 
549
      if(communities[N->community1].P) min_delta_sigma->update(N->community1);
 
550
    }
 
551
 
 
552
    if(N->delta_sigma == min_delta_sigma->delta_sigma[N->community2]) {
 
553
      min_delta_sigma->delta_sigma[N->community2] = communities[N->community2].min_delta_sigma();
 
554
      if(communities[N->community2].P) min_delta_sigma->update(N->community2);
 
555
    }
 
556
  }
 
557
}
 
558
 
 
559
void Communities::add_neighbor(Neighbor* N) {
 
560
  communities[N->community1].add_neighbor(N);
 
561
  communities[N->community2].add_neighbor(N);
 
562
  H->add(N);
 
563
 
 
564
  if(max_memory !=-1) {
 
565
    if(N->delta_sigma < min_delta_sigma->delta_sigma[N->community1]) {
 
566
      min_delta_sigma->delta_sigma[N->community1] = N->delta_sigma;
 
567
      if(communities[N->community1].P) min_delta_sigma->update(N->community1);
 
568
    }
 
569
 
 
570
    if(N->delta_sigma < min_delta_sigma->delta_sigma[N->community2]) {
 
571
      min_delta_sigma->delta_sigma[N->community2] = N->delta_sigma;
 
572
      if(communities[N->community2].P) min_delta_sigma->update(N->community2);
 
573
    }
 
574
  }
 
575
}
 
576
 
 
577
void Communities::update_neighbor(Neighbor* N, float new_delta_sigma) {
 
578
  if(max_memory !=-1) {
 
579
    if(new_delta_sigma < min_delta_sigma->delta_sigma[N->community1]) {
 
580
      min_delta_sigma->delta_sigma[N->community1] = new_delta_sigma;
 
581
      if(communities[N->community1].P) min_delta_sigma->update(N->community1);
 
582
    }
 
583
   
 
584
    if(new_delta_sigma < min_delta_sigma->delta_sigma[N->community2]) {
 
585
      min_delta_sigma->delta_sigma[N->community2] = new_delta_sigma;
 
586
      if(communities[N->community2].P) min_delta_sigma->update(N->community2);
 
587
    }
 
588
 
 
589
    float old_delta_sigma = N->delta_sigma;
 
590
    N->delta_sigma = new_delta_sigma;
 
591
    H->update(N);
 
592
 
 
593
    if(old_delta_sigma == min_delta_sigma->delta_sigma[N->community1]) {
 
594
      min_delta_sigma->delta_sigma[N->community1] = communities[N->community1].min_delta_sigma();
 
595
      if(communities[N->community1].P) min_delta_sigma->update(N->community1);
 
596
    }
 
597
 
 
598
    if(old_delta_sigma == min_delta_sigma->delta_sigma[N->community2]) {
 
599
      min_delta_sigma->delta_sigma[N->community2] = communities[N->community2].min_delta_sigma();
 
600
      if(communities[N->community2].P) min_delta_sigma->update(N->community2);
 
601
    }
 
602
  }
 
603
  else {
 
604
    N->delta_sigma = new_delta_sigma;
 
605
    H->update(N);
 
606
  }
 
607
}
 
608
 
 
609
void Communities::manage_memory() {
 
610
  while((memory_used > max_memory) && !min_delta_sigma->is_empty()) {
 
611
    int c = min_delta_sigma->get_max_community();
 
612
    delete communities[c].P;
 
613
    communities[c].P = 0;
 
614
    min_delta_sigma->remove_community(c);
 
615
  }  
 
616
}
 
617
 
 
618
 
 
619
 
 
620
void Communities::merge_communities(Neighbor* merge_N) {
 
621
  int c1 = merge_N->community1;
 
622
  int c2 = merge_N->community2;
 
623
  
 
624
  communities[nb_communities].first_member = communities[c1].first_member;      // merge the 
 
625
  communities[nb_communities].last_member = communities[c2].last_member;        // two lists   
 
626
  members[communities[c1].last_member] = communities[c2].first_member;          // of members
 
627
 
 
628
  communities[nb_communities].size = communities[c1].size + communities[c2].size;
 
629
  communities[nb_communities].this_community = nb_communities;
 
630
  communities[nb_communities].sub_community_of = 0;
 
631
  communities[nb_communities].sub_communities[0] = c1;
 
632
  communities[nb_communities].sub_communities[1] = c2;
 
633
  communities[nb_communities].total_weight = communities[c1].total_weight + communities[c2].total_weight;
 
634
  communities[nb_communities].internal_weight = communities[c1].internal_weight + communities[c2].internal_weight + merge_N->weight;
 
635
  communities[nb_communities].sigma = communities[c1].sigma + communities[c2].sigma + merge_N->delta_sigma;
 
636
  
 
637
  communities[c1].sub_community_of = nb_communities;
 
638
  communities[c2].sub_community_of = nb_communities;
 
639
 
 
640
// update the new probability vector...
 
641
  
 
642
  if(communities[c1].P && communities[c2].P) communities[nb_communities].P = new Probabilities(c1, c2);
 
643
 
 
644
  if(communities[c1].P) {
 
645
    delete communities[c1].P; 
 
646
    communities[c1].P = 0;
 
647
    if(max_memory != -1) min_delta_sigma->remove_community(c1);
 
648
  }
 
649
  if(communities[c2].P) {
 
650
    delete communities[c2].P;
 
651
    communities[c2].P = 0;
 
652
    if(max_memory != -1) min_delta_sigma->remove_community(c2);
 
653
  }
 
654
 
 
655
  if(max_memory != -1) {
 
656
    min_delta_sigma->delta_sigma[c1] = -1.;                 // to avoid to update the min_delta_sigma for these communities
 
657
    min_delta_sigma->delta_sigma[c2] = -1.;                 // 
 
658
    min_delta_sigma->delta_sigma[nb_communities] = -1.;
 
659
  }
 
660
  
 
661
// update the new neighbors
 
662
// by enumerating all the neighbors of c1 and c2
 
663
 
 
664
  Neighbor* N1 = communities[c1].first_neighbor;
 
665
  Neighbor* N2 = communities[c2].first_neighbor;
 
666
 
 
667
  while(N1 && N2) { 
 
668
    int neighbor_community1;
 
669
    int neighbor_community2;
 
670
    
 
671
    if (N1->community1 == c1) neighbor_community1 = N1->community2;
 
672
    else neighbor_community1 = N1->community1;
 
673
    if (N2->community1 == c2) neighbor_community2 = N2->community2;
 
674
    else neighbor_community2 = N2->community1;
 
675
 
 
676
    if (neighbor_community1 < neighbor_community2) {
 
677
      Neighbor* tmp = N1;
 
678
      if (N1->community1 == c1) N1 = N1->next_community1;
 
679
      else N1 = N1->next_community2;
 
680
      remove_neighbor(tmp);
 
681
      Neighbor* N = new Neighbor;
 
682
      N->weight = tmp->weight;
 
683
      N->community1 = neighbor_community1;
 
684
      N->community2 = nb_communities;
 
685
      N->delta_sigma = (double(communities[c1].size+communities[neighbor_community1].size)*tmp->delta_sigma + double(communities[c2].size)*merge_N->delta_sigma)/(double(communities[c1].size+communities[c2].size+communities[neighbor_community1].size));//compute_delta_sigma(neighbor_community1, nb_communities);
 
686
      N->exact = false;
 
687
      delete tmp;
 
688
      add_neighbor(N); 
 
689
    }
 
690
    
 
691
    if (neighbor_community2 < neighbor_community1) {
 
692
      Neighbor* tmp = N2;
 
693
      if (N2->community1 == c2) N2 = N2->next_community1;
 
694
      else N2 = N2->next_community2;
 
695
      remove_neighbor(tmp);
 
696
      Neighbor* N = new Neighbor;
 
697
      N->weight = tmp->weight;
 
698
      N->community1 = neighbor_community2;
 
699
      N->community2 = nb_communities;
 
700
      N->delta_sigma = (double(communities[c1].size)*merge_N->delta_sigma + double(communities[c2].size+communities[neighbor_community2].size)*tmp->delta_sigma)/(double(communities[c1].size+communities[c2].size+communities[neighbor_community2].size));//compute_delta_sigma(neighbor_community2, nb_communities);
 
701
      N->exact = false;
 
702
      delete tmp;
 
703
      add_neighbor(N); 
 
704
    }
 
705
    
 
706
    if (neighbor_community1 == neighbor_community2) {
 
707
      Neighbor* tmp1 = N1;
 
708
      Neighbor* tmp2 = N2;
 
709
      bool exact = N1->exact && N2->exact;
 
710
      if (N1->community1 == c1) N1 = N1->next_community1;
 
711
      else N1 = N1->next_community2;
 
712
      if (N2->community1 == c2) N2 = N2->next_community1;
 
713
      else N2 = N2->next_community2;
 
714
      remove_neighbor(tmp1);
 
715
      remove_neighbor(tmp2);
 
716
      Neighbor* N = new Neighbor;
 
717
      N->weight = tmp1->weight + tmp2->weight;
 
718
      N->community1 = neighbor_community1;
 
719
      N->community2 = nb_communities;
 
720
      N->delta_sigma = (double(communities[c1].size+communities[neighbor_community1].size)*tmp1->delta_sigma + double(communities[c2].size+communities[neighbor_community1].size)*tmp2->delta_sigma - double(communities[neighbor_community1].size)*merge_N->delta_sigma)/(double(communities[c1].size+communities[c2].size+communities[neighbor_community1].size));
 
721
      N->exact = exact;
 
722
      delete tmp1;
 
723
      delete tmp2;
 
724
      add_neighbor(N);
 
725
    }
 
726
  }
 
727
 
 
728
  
 
729
  if(!N1) {
 
730
    while(N2) {
 
731
//      double delta_sigma2 = N2->delta_sigma;
 
732
      int neighbor_community;
 
733
      if (N2->community1 == c2) neighbor_community = N2->community2;
 
734
      else neighbor_community = N2->community1;
 
735
      Neighbor* tmp = N2;
 
736
      if (N2->community1 == c2) N2 = N2->next_community1;
 
737
      else N2 = N2->next_community2;
 
738
      remove_neighbor(tmp);
 
739
      Neighbor* N = new Neighbor;
 
740
      N->weight = tmp->weight;
 
741
      N->community1 = neighbor_community;
 
742
      N->community2 = nb_communities;
 
743
      N->delta_sigma = (double(communities[c1].size)*merge_N->delta_sigma + double(communities[c2].size+communities[neighbor_community].size)*tmp->delta_sigma)/(double(communities[c1].size+communities[c2].size+communities[neighbor_community].size));//compute_delta_sigma(neighbor_community, nb_communities);
 
744
      N->exact = false;
 
745
      delete tmp;
 
746
      add_neighbor(N);
 
747
    }
 
748
  }
 
749
  if(!N2) {
 
750
    while(N1) {
 
751
//      double delta_sigma1 = N1->delta_sigma;
 
752
      int neighbor_community;
 
753
      if (N1->community1 == c1) neighbor_community = N1->community2;
 
754
      else neighbor_community = N1->community1;
 
755
      Neighbor* tmp = N1;
 
756
      if (N1->community1 == c1) N1 = N1->next_community1;
 
757
      else N1 = N1->next_community2;
 
758
      remove_neighbor(tmp);
 
759
      Neighbor* N = new Neighbor;
 
760
      N->weight = tmp->weight;
 
761
      N->community1 = neighbor_community;
 
762
      N->community2 = nb_communities;
 
763
      N->delta_sigma = (double(communities[c1].size+communities[neighbor_community].size)*tmp->delta_sigma + double(communities[c2].size)*merge_N->delta_sigma)/(double(communities[c1].size+communities[c2].size+communities[neighbor_community].size));//compute_delta_sigma(neighbor_community, nb_communities);
 
764
      N->exact = false;
 
765
      delete tmp;
 
766
      add_neighbor(N);
 
767
    }
 
768
  }
 
769
 
 
770
  if(max_memory != -1) {
 
771
    min_delta_sigma->delta_sigma[nb_communities] = communities[nb_communities].min_delta_sigma();
 
772
    min_delta_sigma->update(nb_communities);
 
773
  } 
 
774
 
 
775
  nb_communities++;
 
776
  nb_active_communities--;
 
777
}
 
778
 
 
779
double Communities::merge_nearest_communities() {
 
780
  Neighbor* N = H->get_first();  
 
781
  while(!N->exact) {
 
782
    update_neighbor(N, compute_delta_sigma(N->community1, N->community2));
 
783
    N->exact = true;
 
784
    N = H->get_first();
 
785
    if(max_memory != -1) manage_memory();
 
786
  }
 
787
 
 
788
  double d = N->delta_sigma;
 
789
  remove_neighbor(N);
 
790
 
 
791
  merge_communities(N);
 
792
  if(max_memory != -1) manage_memory();
 
793
  
 
794
  if (merges) {
 
795
    MATRIX(*merges, mergeidx, 0)=N->community1;
 
796
    MATRIX(*merges, mergeidx, 1)=N->community2;
 
797
    mergeidx++;
 
798
  }
 
799
 
 
800
  if (modularity) {
 
801
    float Q = 0.;
 
802
    for(int i = 0; i < nb_communities; i++) {
 
803
      if(communities[i].sub_community_of == 0) {
 
804
        Q += (communities[i].internal_weight - communities[i].total_weight*communities[i].total_weight/G->total_weight)/G->total_weight;
 
805
      }
 
806
    }
 
807
    VECTOR(*modularity)[mergeidx]=Q;
 
808
  }
 
809
 
 
810
  delete N;
 
811
 
 
812
  /* This could use igraph_progress */
 
813
/*   if(!silent) { */
 
814
/*     for(int k = (500*(G->nb_vertices - nb_active_communities - 1))/(G->nb_vertices-1) + 1; k <= (500*(G->nb_vertices - nb_active_communities))/(G->nb_vertices-1); k++) { */
 
815
/*       if(k % 50 == 1) {cerr.width(2); cerr << endl << k/ 5 << "% ";} */
 
816
/*       cerr << "."; */
 
817
/*     } */
 
818
/*   } */
 
819
  return d;
 
820
}
 
821
 
 
822
double Communities::compute_delta_sigma(int community1, int community2) {
 
823
  if(!communities[community1].P) {
 
824
    communities[community1].P = new Probabilities(community1);
 
825
    if(max_memory != -1) min_delta_sigma->update(community1);
 
826
  }
 
827
  if(!communities[community2].P) {
 
828
    communities[community2].P = new Probabilities(community2);
 
829
    if(max_memory != -1) min_delta_sigma->update(community2);
 
830
  }
 
831
  
 
832
  return communities[community1].P->compute_distance(communities[community2].P)*double(communities[community1].size)*double(communities[community2].size)/double(communities[community1].size + communities[community2].size);
 
833
}