netcomplexity.h
Go to the documentation of this file.
1 /*
2  @copyright Russell Standish 2000-2013
3  @author Russell Standish
4  This file is part of EcoLab
5 
6  Open source licensed under the MIT license. See LICENSE for details.
7 */
8 
12 #ifndef NETCOMPLEXITY_H
13 #define NETCOMPLEXITY_H
14 
15 #ifdef MPI_SUPPORT
16 #include <classdescMP.h>
17 using classdesc::MPIbuf;
18 #endif
19 #include <pack_base.h>
20 
21 #include <ref.h>
22 #include <graph.h>
23 
24 #include <ref.h>
25 #include <graph.h>
26 
27 #include <vector>
28 #include <map>
29 #include <set>
30 #include <utility>
31 #include <algorithm>
32 #include <iostream>
33 #include <stdio.h>
34 
35 using std::vector;
36 using std::map;
37 using std::pair;
38 using std::insert_iterator;
39 using std::max_element;
40 using std::set_difference;
41 using std::ostream;
42 using std::cout;
43 using std::endl;
44 using std::max;
45 
46 #include <math.h>
47 #include <assert.h>
48 
49 /*
50  bits needed from nauty.h, macros replaced with inline C++ functions
51 */
52 #include "nauty_sizes.h"
53 
54 
55 namespace ecolab
56 {
57 #if SIZEOF_LONG>4
58 #define WORDSIZE 64
59 #else
60 #define WORDSIZE 32
61 #endif
62 
63 #if WORDSIZE==32
64 #if SIZEOF_INT>=4
65  typedef unsigned int setword;
66 #define SETWORD_INT
67 #else
68  typedef unsigned long setword;
69 #define SETWORD_LONG
70 #endif
71 #endif
72 
73 #if WORDSIZE==64
74 #if SIZEOF_INT>=8
75  typedef unsigned int setword;
76 #define SETWORD_INT
77 #else
78 #if SIZEOF_LONG>=8
79  typedef unsigned long setword;
80 #define SETWORD_LONG
81 #else
82  typedef unsigned long long setword;
83 #define SETWORD_LONGLONG
84 #endif
85 #endif
86 #endif
87 
88 #if (WORDSIZE==64)
89  const setword MSK64=0xFFFFFFFFFFFFFFFFUL;
90  const setword MSK63C=0x7FFFFFFFFFFFFFFFUL;
91  const setword ALLBITS=MSK64;
92  inline unsigned SETWD(unsigned pos) {return pos>>6;}
93  inline unsigned SETBT(unsigned pos) {return pos&0x3F;}
94  inline setword BITMASK(unsigned x) {return MSK63C >> x;}
95 
96 #ifdef SETWORD_LONGLONG
97  const setword bitt[] = {01000000000000000000000LL,0400000000000000000000LL,
98  0200000000000000000000LL,0100000000000000000000LL,
99  040000000000000000000LL,020000000000000000000LL,
100  010000000000000000000LL,04000000000000000000LL,
101  02000000000000000000LL,01000000000000000000LL,
102  0400000000000000000LL,0200000000000000000LL,
103  0100000000000000000LL,040000000000000000LL,
104  020000000000000000LL,010000000000000000LL,
105  04000000000000000LL,02000000000000000LL,
106  01000000000000000LL,0400000000000000LL,0200000000000000LL,
107  0100000000000000LL,040000000000000LL,020000000000000LL,
108  010000000000000LL,04000000000000LL,02000000000000LL,
109  01000000000000LL,0400000000000LL,0200000000000LL,
110  0100000000000LL,040000000000LL,020000000000LL,010000000000LL,
111  04000000000LL,02000000000LL,01000000000LL,0400000000LL,
112  0200000000LL,0100000000LL,040000000LL,020000000LL,
113  010000000LL,04000000LL,02000000LL,01000000LL,0400000LL,
114  0200000LL,0100000LL,040000LL,020000LL,010000LL,04000LL,
115  02000LL,01000LL,0400LL,0200LL,0100LL,040LL,020LL,010LL,
116  04LL,02LL,01LL};
117 #else
118  const setword bitt[] = {01000000000000000000000,0400000000000000000000,
119  0200000000000000000000,0100000000000000000000,
120  040000000000000000000,020000000000000000000,
121  010000000000000000000,04000000000000000000,
122  02000000000000000000,01000000000000000000,
123  0400000000000000000,0200000000000000000,
124  0100000000000000000,040000000000000000,020000000000000000,
125  010000000000000000,04000000000000000,02000000000000000,
126  01000000000000000,0400000000000000,0200000000000000,
127  0100000000000000,040000000000000,020000000000000,
128  010000000000000,04000000000000,02000000000000,
129  01000000000000,0400000000000,0200000000000,0100000000000,
130  040000000000,020000000000,010000000000,04000000000,
131  02000000000,01000000000,0400000000,0200000000,0100000000,
132  040000000,020000000,010000000,04000000,02000000,01000000,
133  0400000,0200000,0100000,040000,020000,010000,04000,
134  02000,01000,0400,0200,0100,040,020,010,04,02,01};
135 #endif
136 #else
137  const setword MSK32=0xFFFFFFFFUL;
138  const setword MSK31C=0x7FFFFFFFUL;
139  const setword ALLBITS=MSK32;
140  inline unsigned SETWD(unsigned pos) {return pos>>5;}
141  inline unsigned SETBT(unsigned pos) {return pos&0x1F;}
142  inline setword BITMASK(unsigned x) {return MSK31C >> x;}
143 
144  const setword bitt[] = {020000000000,010000000000,04000000000,02000000000,
145  01000000000,0400000000,0200000000,0100000000,040000000,
146  020000000,010000000,04000000,02000000,01000000,0400000,
147  0200000,0100000,040000,020000,010000,04000,02000,01000,
148  0400,0200,0100,040,020,010,04,02,01};
149 #endif
150 
151 
152  inline void ADDELEMENT(setword* setadd,unsigned pos)
153  {setadd[SETWD(pos)] |= bitt[SETBT(pos)];}
154 
155  inline void DELELEMENT(setword *setadd, unsigned pos)
156  {setadd[SETWD(pos)] &= ~bitt[SETBT(pos)];}
157 
158  inline void FLIPELEMENT(setword *setadd, unsigned pos)
159  {setadd[SETWD(pos)] ^= bitt[SETBT(pos)];}
160 
161  inline bool ISELEMENT(const setword *setadd, unsigned pos)
162  {return (setadd[SETWD(pos)] & bitt[SETBT(pos)]) != 0;}
163 
164  inline void EMPTYSET(setword *setadd, unsigned m)
165  {
166  setword *es;
167  for (es = setadd+m; --es >= setadd;) *es=0;
168  }
169 
170 #include <pack_base.h>
171 #include <pack_stl.h>
172 #include <classdesc_access.h>
173 #ifdef MPI_SUPPORT
174 #include <classdescMP.h>
175 #endif
176 
177  /* iota seems to be missing? */
178  template <class T>
179  void iota(typename T::iterator p, const typename T::iterator& end, typename T::value_type val)
180  {for (; p!=end; p++, val++) *p=val;}
181 
182  inline void iota(vector<unsigned>::iterator p, const vector<unsigned>::iterator& end, unsigned val)
183  {for (; p!=end; p++, val++) *p=val;}
184 
185  inline void iota(unsigned* p, const unsigned * end, unsigned val)
186  {for (; p!=end; p++, val++) *p=val;}
187 
188 
189  /* wrap nauty's bitset class with bitvector syntax */
190  class bitref
191  {
192  setword *addr;
193  unsigned idx;
194  CLASSDESC_ACCESS(bitref);
195  public:
196  bitref(setword *a, unsigned i): addr(a), idx(i) {}
197  operator bool () const {return ISELEMENT(addr,idx);}
198  bitref& operator=(bool x) {x? ADDELEMENT(addr,idx): DELELEMENT(addr,idx); return *this;}
199  bitref& operator=(const bitref& x) {return (*this)=(bool)x;}
200  bool operator==(const bitref x) const {return (bool)(*this)==(bool)x;}
201  bool operator!=(const bitref x) const {return (bool)(*this)==(bool)x;}
202  bool operator<(const bitref x) const {return (bool)(*this)<(bool)x;}
203  };
204 
205  class const_bitref
206  {
207  const setword *addr;
208  unsigned idx;
209  CLASSDESC_ACCESS(const_bitref);
210  public:
211  const_bitref(const setword *a, unsigned i): addr(a), idx(i) {}
212  operator bool () const {return ISELEMENT(addr,idx);}
213  bool operator==(const bitref x) const {return (bool)(*this)==(bool)x;}
214  bool operator!=(const bitref x) const {return (bool)(*this)==(bool)x;}
215  bool operator<(const bitref x) const {return (bool)(*this)<(bool)x;}
216  };
217 
218  class NautyRep;
219 
220  class bitvect
221  {
222  unsigned sz;
223  std::vector<setword> data;
224 
225  unsigned nwords() const {return sz? div(sz-1)+1: 0;}
226  unsigned nbytes() const {return nwords()*(WORDSIZE>>3);}
227  unsigned div(unsigned i) const {return SETWD(i);} /* i/WORDSIZE; */
228  unsigned mod(unsigned i) const {return SETBT(i);} /* i%WORDSIZE; */
229  unsigned mod(int i) const {return SETBT(i);} /* i%WORDSIZE; */
230  void null_last_bits() {if (sz) data[nwords()-1]&=~BITMASK(mod(sz-1));}
231  bool last_bits_null() const {
232  return sz==0 || !(data[nwords()-1]&BITMASK(mod(sz-1)));
233  }
234 
235  //needs to reach in to the data member
236  friend void ecolab_nauty(const NautyRep&,NautyRep&,double&,bool,int);
237  CLASSDESC_ACCESS(bitvect);
238  public:
239 
240  bitvect(): sz(0), data(0) {}
241  bitvect(const bitvect& x): sz(x.sz), data(x.data) {assert(last_bits_null());}
242  bitvect(unsigned s, bool x=false): sz(s), data(nwords(),x?ALLBITS:0) {
243  null_last_bits();}
244  template <class T>
245  bitvect(const vector<T>& x): sz(x.sz), data(nwords()) {
246  for (unsigned i=0; i<sz; i++) (*this)[i]=x[i];
247  assert(last_bits_null());
248  }
249 
250  bitvect& operator|=(const bitvect& x) {
251  assert(x.data.size()==data.size());
252  for (unsigned i=0; i<x.data.size(); ++i) data[i]|=x.data[i];
253  assert(last_bits_null());
254  return *this;
255  }
256  bitvect operator|(bitvect x) const {return x|=*this;}
257 
258  bitvect& operator&=(const bitvect& x) {
259  assert(x.data.size()==data.size());
260  for (unsigned i=0; i<x.data.size(); ++i) data[i]&=x.data[i];
261  assert(last_bits_null());
262  return *this;
263  }
264 
265  bitvect operator&(bitvect x) const {return x&=*this;}
266 
267  bitvect& operator^=(const bitvect& x) {
268  assert(x.data.size()==data.size());
269  for (unsigned i=0; i<x.data.size(); ++i) data[i]^=x.data[i];
270  null_last_bits();
271  return *this;
272  }
273  bitvect operator^(bitvect x) const {return x^=*this;}
274 
275  bitvect operator~() const {
276  bitvect r(*this);
277  return r.flip();
278  }
279 
280  bitvect& flip() {
281  for (unsigned i=0; i<data.size(); ++i) data[i]=~data[i];
282  null_last_bits();
283  return *this;
284  }
285 
286  bitvect operator>>(unsigned offs) const {
287  bitvect r(sz+offs);
288  for (unsigned i=0; i<nwords()-1; ++i) {
289  r.data[div(i*WORDSIZE+offs)+1] = data[i]>>mod(offs);
290  r.data[div(i*WORDSIZE+offs)] = data[i]<<mod(-offs);
291  }
292  r.data[div((nwords()-1)*WORDSIZE+offs)] = data[nwords()-1]<<mod(-offs);
293  return r;
294  }
295  bitvect num(unsigned start,unsigned nbits) const;
296  bool equal(const bitvect& y, unsigned len, unsigned offs) const;
297  void setrange(unsigned start,unsigned end,bool value);
298  bitref operator[] (unsigned i) {assert(i<sz); return bitref(&data[0],i);}
299  const_bitref operator[] (unsigned i) const {assert(i<sz); return const_bitref(&data[0],i);}
300  unsigned size() const {return sz;}
302  void operator++(int);
305  bool next_perm();
306  unsigned popcnt(unsigned start=0, int end=-1 /* -1 == sz */) const;
307  bool operator==(const bitvect& x) const
308  {
309  bool r=x.sz==sz;
310  for (unsigned i=0; r&&i<nwords(); i++) r&=x.data[i]==data[i];
311  return r;
312  }
313  bool operator!=(const bitvect& x) const {return !((*this)==x);}
314  bool operator<(const bitvect& x) const
315  {
316  assert(last_bits_null());
317  bool r=x.sz==sz;
318  unsigned i;
319  for (i=0; i<nwords(); i++) if (!(r&=x.data[i]==data[i])) break;
320  return !r && (sz<x.sz || (sz==x.sz && data[i]<x.data[i]));
321  }
322  void rand() { /* fill vector with random bits */
323  for (unsigned i=0; i<nwords(); data[i++]=::rand());
324  null_last_bits();
325  }
326  };
327 
328 
329 #ifdef _CLASSDESC
330 #pragma omit pack ecolab::bitref
331 #pragma omit unpack ecolab::bitref
332 #pragma omit pack ecolab::const_bitref
333 #pragma omit unpack ecolab::const_bitref
334 #endif
335 
336  /* return the number corresponding the nbit field pointed to by i */
337  bitvect num(const bitvect& x, unsigned i,unsigned nbits);
338 
339  /*
340  Common assignment between all of the *Rep types
341  */
342  template <class T, class U>
343  void asgRep(T& x, const U& y)
344  {
345  if (x.nodes()!=y.nodes()) x=T(y.nodes());
346  for (unsigned i=0; i<y.nodes(); ++i)
347  for (unsigned j=0; j<y.nodes(); ++j)
348  if (i!=j) x(i,j)=y(i,j);
349  }
350 
352  template <class T, class Graph>
353  void asgRepGraph(T& x, const Graph& y)
354  {
355  x=T(y.nodes());
356  for (typename Graph::const_iterator i=y.begin(); i!=y.end(); ++i)
357  x(i->source(),i->target())=true;
358  }
359 
361  template <class Rep>
362  class Rep_const_iterator
363  {
364  Edge e;
365  unsigned &i, &j;
366  const Rep* rep;
367  public:
368  Rep_const_iterator(const Rep& rp, bool end):
369  e(0, 1), i(e.source()), j(e.target()), rep(&rp)
370  {
371  if (end||rep->nodes()<2) {i=rep->nodes(); j=0;}
372  else if (!(*rep)(i,j)) operator++();//initialise ref to point to first arc
373  }
374  Rep_const_iterator(const Rep_const_iterator& x):
375  i(e.source()), j(e.target()) {*this=x;}
376  const Rep_const_iterator& operator=(const Rep_const_iterator& x)
377  {e=x.e; rep=x.rep; return *this;}
378 
379  Edge operator*() const {return e;}
380  const Edge* operator->() const {return &e;}
381  void operator++() { //skip to next edge
382  do {
383  ++j;
384  if (j==i) ++j;
385  if (j>=rep->nodes())
386  {
387  j=0;
388  ++i;
389  }
390  } while (i<rep->nodes() && !(*rep)(i,j));
391  }
392  bool operator==(const Rep_const_iterator& x) const {
393  return e==x.e && rep==x.rep;
394  }
395  bool operator!=(const Rep_const_iterator& x) const {
396  return !operator==(x);
397  }
398  };
399 
400  class BiDirectionalBitRep;
401  class NautyRep;
402 
406  class BitRep
407  {
408  bitvect linklist;
409  unsigned nodecnt;
410  CLASSDESC_ACCESS(BitRep);
411  unsigned idx(unsigned i, unsigned j) const {
412  assert(i!=j);
413  return (j<i)? (nodecnt-1)*i+j: (nodecnt-1)*i+j-1;
414  }
415  public:
416  typedef Rep_const_iterator<BitRep> const_iterator;
417  const_iterator begin() const {return const_iterator(*this,false);}
418  const_iterator end() const {return const_iterator(*this,true);}
419 
422  BitRep(unsigned nodes=0, unsigned links=0):
423  linklist( nodes * (nodes-1), false), nodecnt(nodes) {
424  assert(links<=linklist.size());
425  linklist.setrange(0,links,true);
426  }
427  BitRep(const Graph& x): nodecnt(0) {asgRepGraph(*this,x);}
428  BitRep(const DiGraph& x): nodecnt(0) {asgRepGraph(*this,x);}
429  BitRep(const BiDirectionalGraph& x): nodecnt(0) {asgRepGraph(*this,x);}
430  BitRep(const NautyRep& x): nodecnt(0) {asgRep(*this,x);}
431  BitRep(const BiDirectionalBitRep& x): nodecnt(0) {asgRep(*this,x);}
432  const BitRep& operator=(const DiGraph& x) {asgRepGraph(*this,x); return *this;}
433  const BitRep& operator=(const BiDirectionalGraph& x) {
434  asgRepGraph(*this,x); return *this;}
435  const BitRep& operator=(const NautyRep& x) {asgRep(*this,x); return *this;}
436  const BitRep& operator=(const BiDirectionalBitRep& x)
437  {asgRep(*this,x); return *this;}
438 
439  const BitRep& operator|=(const BitRep& x)
440  {linklist|=x.linklist; return *this;}
441  BitRep operator|(BitRep x) const {return x|=*this;}
442  const BitRep& operator&=(const BitRep& x)
443  {linklist&=x.linklist; return *this;}
444  BitRep operator&(BitRep x) const {return x&=*this;}
445  const BitRep& operator^=(const BitRep& x)
446  {linklist^=x.linklist; return *this;}
447  BitRep operator^(BitRep x) const {return x^=*this;}
448  BitRep operator~() const
449  {BitRep r; r.nodecnt=nodecnt; r.linklist=~linklist; return r;}
450 
451  BitRep operator+=(const BitRep& x) {
452  linklist|=x.linklist>>nodecnt;
453  nodecnt+=x.nodecnt;
454  return *this;
455  }
456  BiDirectionalBitRep symmetrise();
457 
458  bitref operator()(unsigned i, unsigned j) {
459  return linklist[idx(i,j)];}
460  const_bitref operator()(unsigned i, unsigned j) const {
461  return linklist[idx(i,j)];}
463  unsigned nodes() const {return nodecnt;} /* number of nodes */
465  unsigned links() const {return linklist.popcnt();}
467  bool next_perm() {return linklist.next_perm();}
468  bool operator==(const BitRep& x) const {
469  return nodecnt==x.nodecnt && linklist==x.linklist;
470  }
471  bool operator!=(const BitRep& x) const {return !operator==(x);}
472  bool operator<(const BitRep& x) const {
473  return nodecnt<x.nodecnt || (nodecnt==x.nodecnt && linklist<x.linklist);
474  }
475  bool contains(const Edge& e) const {return operator()(e.source(),e.target());}
476  void push_back(const Edge& e) {operator()(e.source(),e.target())=true;}
477  void clear (unsigned nodes=0) {
478  nodecnt=nodes; linklist.setrange( 0, nodes * (nodes-1), false);}
479  };
480 
484  class NautyRep
485  {
486  unsigned nodecnt;
487  CLASSDESC_ACCESS(NautyRep);
488  friend void ecolab_nauty(const NautyRep&,NautyRep&,double&,bool,int);
489  public:
490  typedef Rep_const_iterator<NautyRep> const_iterator;
491  const_iterator begin() const {return const_iterator(*this,false);}
492  const_iterator end() const {return const_iterator(*this,true);}
493 
494  unsigned m() const {return SETWD(nodecnt-1)+1;}
495  unsigned mw() const {return WORDSIZE*m();}
496  bitvect linklist;
497 
498  NautyRep(): nodecnt(0) {}
499  NautyRep(unsigned nodes):
500  nodecnt(nodes),linklist( nodes * mw(), false) {
501  }
502  NautyRep(const Graph& x): nodecnt(0) {asgRepGraph(*this,x);}
503  NautyRep(const DiGraph& x): nodecnt(0) {asgRepGraph(*this,x);}
504  NautyRep(const BiDirectionalGraph& x): nodecnt(0) {asgRepGraph(*this,x);}
505  NautyRep(const BitRep& x): nodecnt(0) {asgRep(*this,x);}
506  NautyRep(const BiDirectionalBitRep& x): nodecnt(0) {asgRep(*this,x);}
507  const NautyRep& operator=(const Graph& x) {
508  asgRepGraph(*this,x); return *this;}
509  const NautyRep& operator=(const DiGraph& x) {
510  asgRepGraph(*this,x); return *this;}
511  const NautyRep& operator=(const BiDirectionalGraph& x) {
512  asgRepGraph(*this,x); return *this;}
513  const NautyRep& operator=(const BitRep& x) {asgRep(*this,x); return *this;}
514  const NautyRep& operator=(const BiDirectionalBitRep& x) {
515  asgRep(*this,x); return *this;}
516 
517  // BiDirectionalBitRep symmetrise();
518  bitref operator()(unsigned i, unsigned j) {
519  return linklist[mw()*i+j];}
520  const_bitref operator()(unsigned i, unsigned j) const {
521  return linklist[mw()*i+j];}
523  unsigned nodes() const {return nodecnt;} /* number of nodes */
525  unsigned links() const {return linklist.popcnt();}
527  bool operator==(const NautyRep& x) const {
528  return nodecnt==x.nodecnt && linklist==x.linklist;
529  }
530  bool operator!=(const NautyRep& x) const {return !operator==(x);}
531  bool operator<(const NautyRep& x) const {
532  return nodecnt<x.nodecnt || (nodecnt==x.nodecnt && linklist<x.linklist);
533  }
535  double lnomega() const;
538  NautyRep canonicalise() const;
539  bool contains(const Edge& e) const {return operator()(e.source(),e.target());}
540  void push_back(const Edge& e) {operator()(e.source(),e.target())=true;}
541  void clear (unsigned nodes=0) {
542  nodecnt=nodes; linklist.setrange( 0, nodes * (nodes-1), false);}
543  };
544 
548  class BiDirectionalBitRep
549  {
550  bitvect linklist;
551  unsigned nodecnt;
552  unsigned idx(unsigned i, unsigned j) const {
553  return (j<i)? ((i*(i-1))>>1) + j: ((j*(j-1))>>1)+i;
554  }
555  CLASSDESC_ACCESS(BiDirectionalBitRep);
556  public:
557  typedef Rep_const_iterator<BiDirectionalBitRep> const_iterator;
558  const_iterator begin() const {return const_iterator(*this,false);}
559  const_iterator end() const {return const_iterator(*this,true);}
560 
561  operator BitRep() const {BitRep r; asgRep(r,*this); return r;}
564  BiDirectionalBitRep(unsigned nodes=0, unsigned links=0):
565  linklist( nodes*(nodes-1)/2, false), nodecnt(nodes) {
566  assert(links<=linklist.size());
567  linklist.setrange(0,links,true);
568  }
569  BiDirectionalBitRep(const Graph& x): nodecnt(0) {asgRepGraph(*this,x);}
570  BiDirectionalBitRep(const DiGraph& x): nodecnt(0) {asgRepGraph(*this,x);}
571  BiDirectionalBitRep(const BiDirectionalGraph& x): nodecnt(0) {
572  asgRepGraph(*this,x);}
573  BiDirectionalBitRep(const NautyRep& x): nodecnt(0) {asgRep(*this,x);}
574  BiDirectionalBitRep(const BitRep& x): nodecnt(0) {asgRep(*this,x);}
575  const BiDirectionalBitRep& operator=(const DiGraph& x)
576  {asgRepGraph(*this,x); return *this;}
577  const BiDirectionalBitRep& operator=(const BiDirectionalGraph& x)
578  {asgRepGraph(*this,x); return *this;}
579  const BiDirectionalBitRep& operator=(const NautyRep& x)
580  {asgRep(*this,x); return *this;}
581  const BiDirectionalBitRep& operator=(const BitRep& x)
582  {asgRep(*this,x); return *this;}
583 
584  bitref operator()(unsigned i, unsigned j) {
585  return linklist[idx(i,j)];}
586  const_bitref operator()(unsigned i, unsigned j) const {
587  return linklist[idx(i,j)];}
589  unsigned nodes() const {return nodecnt;} /* number of nodes */
591  unsigned links() const {return 2*linklist.popcnt();}
593  bool next_perm() {return linklist.next_perm();}
594  bool operator==(const BiDirectionalBitRep& x) const {
595  return nodecnt==x.nodecnt && linklist==x.linklist;
596  }
597  bool operator!=(const BiDirectionalBitRep& x) const {return !operator==(x);}
598  bool operator<(const BiDirectionalBitRep& x) const {
599  return nodecnt<x.nodecnt || (nodecnt==x.nodecnt && linklist<x.linklist);
600  }
601  bool contains(const Edge& e) const {return operator()(e.source(),e.target());}
602  void push_back(const Edge& e) {operator()(e.source(),e.target())=true;}
603  void clear (unsigned nodes=0) {
604  nodecnt=nodes; linklist.setrange( 0, nodes * (nodes-1), false);}
605  };
606 
607  template <> inline
608  bool GraphAdaptor<BiDirectionalBitRep>::directed() const {return false;}
609 
610  BitRep permute(const BitRep& x, const vector<unsigned>& perm);
611 
619  void ecolab_nauty(const NautyRep& g, NautyRep& canonical, double& lnomega,
620  bool do_canonical, int invMethod=0);
621 
622  // call saucy for lnomega (if available)
623  double saucy_lnomega(const DiGraph&);
624  // compute lnomega using bliss \a method
625  double igraph_lnomega(const DiGraph& g, int method);
626 
627  inline NautyRep canonicalise(const NautyRep& x) {return x.canonicalise();}
628 
629  /* the above operators allow TCL_obj to work */
630 #ifdef _CLASSDESC
631 #pragma omit TCL_obj BitRep
632 #pragma omit TCL_obj BiDirectionalBitRep
633 #endif
634 
635  BitRep start_perm(unsigned nodes, unsigned links);
636  double lnfactorial(unsigned x);
637  double lnCombine(unsigned x, unsigned y);
638 
639  inline double baseComplexity(unsigned n,unsigned l, bool bidirectional=false)
640  {
641  double ilog2=1/log(2.0);
642  unsigned N=n*(n-1);
643  //cout << "lnCombine("<<N<<","<<l<<")="<<lnCombine(N,l)<<", lnomega="<<x.lnomega()<<endl;
644  return 2*ceil(log(n+1)*ilog2) + 1 + ceil(log(N)*ilog2) +
645  ilog2 * (bidirectional? lnCombine(N/2,l/2): lnCombine(N,l));
646  }
647 
648  inline double complexity(NautyRep x, bool bidirectional=false)
649  {
650  double ilog2=1/log(2.0);
651  // cout << "NautyRep omega="<<exp(x.lnomega())<<endl;
652  return baseComplexity(x.nodes(),x.links(),bidirectional) - ilog2 * x.lnomega();
653  }
654 
662  double canonical(BitRep& canon, const Graph& g);
663 
664  template <class G> double canonical(BitRep& canon, const GraphAdaptor<G>& g)
665  {return canonical(canon, static_cast<const Graph&>(g));}
666 
667  template <class G> double canonical(BitRep& canon, const ConcreteGraph<G>& g)
668  {return canonical(canon, static_cast<const Graph&>(g));}
669 
670  template <class G> double canonical(BitRep& canon, const G& g)
671  {return canonical(canon, GraphAdaptor<const G>(g));}
672 
674 
675  double complexity(const Graph& g);
676  template <class G> double complexity(const GraphAdaptor<G>& g)
677  {return complexity(static_cast<const Graph&>(g));}
678 
679  double complexity(const Graph& g);
680  template <class G> double complexity(const ConcreteGraph<G>& g)
681  {return complexity(static_cast<const Graph&>(g));}
682 
683  template <class G> double complexity(const G& g)
684  {return complexity(GraphAdaptor<const G>(g));}
686 
687 #if WORDSIZE==64
688  // shift right bug (~0UL)>>64 gives the wrong answer
689  inline setword bitrange(unsigned i) {return i==64? 0: ALLBITS>>i;}
690 #else
691  // shift right bug (~0U)>>32 gives the wrong answer
692  inline setword bitrange(unsigned i) {return i==32? 0: ALLBITS>>i;}
693 #endif
694 
695  inline bitvect num(const bitvect& x, unsigned start,unsigned nbits)
696  {
697  return x.num(start,nbits);
698  }
699 
700  inline bitvect bitvect::num(unsigned start,unsigned nbits) const
701  {
702  bitvect r(nbits);
703  if (start+nbits<WORDSIZE)
704  r.data[0]= (data[0]<<start) & ~bitrange(nbits);
705  else
706  for (unsigned i=0; i<nbits; i++)
707  r[i]=(*this)[start+i];
708  return r;
709  }
710 
711  template <class T>
712  vector<T> num(const vector<T>& x, unsigned start,unsigned nbits)
713  {return vector<T>(x.begin()+start,x.begin()+start+nbits);}
714 
715  inline bool equal(const bitvect& x, const bitvect& y, unsigned len, unsigned offs)
716  {
717  return x.equal(y,len,offs);
718  }
719 
720  inline bool bitvect::equal(const bitvect& y, unsigned len, unsigned offs) const
721  {
722  bool r=true;
723  if (offs+len>WORDSIZE)
724  {
725  unsigned i;
726  for (i=0; r && i<div(len); i++)
727  r&= data[i]==(y.data[div(offs)+i]<<mod(offs)|
728  y.data[div(offs)+i+1]>>(WORDSIZE-mod(offs)));
729  unsigned lastword= (div(offs)+i+1<y.nwords())?
730  y.data[div(offs)+i+1]>>(WORDSIZE-mod(offs)): 0;
731  r&= !(~bitrange(mod(len)) &
732  (data[i] ^ (y.data[div(offs)+i]<<mod(offs)| lastword)));
733  }
734  else
735  r&= !(~bitrange(len) & (data[0] ^ (y.data[0]<<offs)));
736  return r;
737  }
738 
739  template <class T>
740  bool equal(const T& x, const T& y, unsigned len, unsigned offs)
741  {
742  // return std::equal(x.begin(),x.begin()+len,y.begin()+offs);
743  bool r=true;
744  for (unsigned i=0; r&&i<len; i++)
745  r&=x[i]==y[offs+i];
746  return r;
747  }
748 
749  template <class T>
750  void setvec(vector<T>& x,const bitvect& y)
751  {
752  x.resize(y.size());
753  for (unsigned i=0; i<y.size(); i++)
754  x[i]=y[i];
755  }
756 
758  double ODcomplexity(const DiGraph& x);
759 
761  double MA(const DiGraph& x);
762 
763  inline BiDirectionalBitRep BitRep::symmetrise()
764  {BiDirectionalBitRep r; asgRep(r,*this); return r;}
765 
766 } //namespace ecolab
767 
768 ostream& operator<<(ostream& s, const ecolab::bitvect& x);
769 inline ostream& operator<<(ostream& s, const ecolab::BitRep& x)
770 {return s<<ecolab::DiGraph(x);}
771 inline ostream& operator<<(ostream& s, const ecolab::BiDirectionalBitRep& x)
772 {return s<<ecolab::BiDirectionalGraph(x);}
773 inline ostream& operator<<(ostream& s, const ecolab::NautyRep& x)
774 {return s<<ecolab::DiGraph(x);}
775 
776 #include "netcomplexity.cd"
777 #endif
classdesc_access.h
descriptor access to a class&#39;s privates
ecolab::Graph::begin
virtual const_iterator begin() const =0
iterator to first edge
ecolab::bitref
Definition: netcomplexity.h:190
ecolab::Rep_const_iterator
An iterator suitable for "duck-typing" with Graph.
Definition: netcomplexity.h:362
ecolab::NautyRep::nodes
unsigned nodes() const
number of nodes this graph has
Definition: netcomplexity.h:523
ecolab::BiDirectionalBitRep::BiDirectionalBitRep
BiDirectionalBitRep(unsigned nodes=0, unsigned links=0)
Definition: netcomplexity.h:564
ecolab::Graph::const_iterator
iterator over edges
Definition: graph.h:111
ecolab::NautyRep::canonicalise
NautyRep canonicalise() const
ecolab::BitRep::next_perm
bool next_perm()
next permutation of links - returns false if no further permutations
Definition: netcomplexity.h:467
ecolab::ecolab_nauty
void ecolab_nauty(const NautyRep &g, NautyRep &canonical, double &lnomega, bool do_canonical, int invMethod=0)
EcoLab interface to Nauty.
ecolab::Graph::nodes
virtual unsigned nodes() const =0
number of nodes
ecolab::DiGraph
Definition: graph.h:272
ecolab::NautyRep::operator==
bool operator==(const NautyRep &x) const
next permutation of links - returns false if no further permutations
Definition: netcomplexity.h:527
classdesc::MPIbuf
buffer object providing MPI functionality
Definition: classdescMP.h:75
ref.h
Reference counted smart pointer classes.
ecolab::MA
double MA(const DiGraph &x)
Medium articulation.
ecolab::BiDirectionalGraph
A graph in which each link is bidirectional.
Definition: graph.h:375
pack_base.h
serialisation descriptor
ecolab::BiDirectionalBitRep::nodes
unsigned nodes() const
number of nodes this graph has
Definition: netcomplexity.h:589
ecolab::asgRepGraph
void asgRepGraph(T &x, const Graph &y)
convert between an EcoLab Graph type and a BitRep type
Definition: netcomplexity.h:353
ecolab::ConcreteGraph
Definition: graph.h:257
ecolab::GraphAdaptor::directed
bool directed() const
true if a bidirectional graph (all edges go both ways)
Definition: graph.h:269
ecolab::BiDirectionalBitRep
Definition: netcomplexity.h:548
ecolab::BitRep::BitRep
BitRep(unsigned nodes=0, unsigned links=0)
Definition: netcomplexity.h:422
ecolab::GraphAdaptor
Definition: graph.h:211
ecolab::Graph
Definition: graph.h:94
ecolab::bitvect
Definition: netcomplexity.h:220
ecolab::idx
Definition: TCL_obj_stl.h:185
ecolab::BiDirectionalBitRep::next_perm
bool next_perm()
next permutation of links - returns false if no further permutations
Definition: netcomplexity.h:593
CLASSDESC_ACCESS
#define CLASSDESC_ACCESS(type)
add friend statements for each accessor function
Definition: classdesc_access.h:36
ecolab::BitRep::nodes
unsigned nodes() const
number of nodes this graph has
Definition: netcomplexity.h:463
ecolab::NautyRep::lnomega
double lnomega() const
log of the number of bitstrings equivalent to this graph
ecolab::canonical
double canonical(BitRep &canon, const Graph &g)
SuperNOVA automorphism algorithm:
graph.h
EcoLab graph library.
ecolab::NautyRep
Definition: netcomplexity.h:484
ecolab::BitRep::links
unsigned links() const
number of links this graph has
Definition: netcomplexity.h:465
ecolab::NautyRep::links
unsigned links() const
number of links this graph has
Definition: netcomplexity.h:525
pack_stl.h
serialisation for standard containers
classdescMP.h
MPI parallel processing library.
ecolab::const_bitref
Definition: netcomplexity.h:205
ecolab
_OPENMP
Definition: accessor.h:16
ecolab::Graph::end
virtual const_iterator end() const =0
iterator beyond last edge
ecolab::ODcomplexity
double ODcomplexity(const DiGraph &x)
Offdiagonal complexity.
ecolab::BitRep
Definition: netcomplexity.h:406
ecolab::Edge
Definition: graph.h:45
ecolab::BiDirectionalBitRep::links
unsigned links() const
number of links this graph has
Definition: netcomplexity.h:591