DMRITool  v0.1.1-139-g860d86b4
Diffusion MRI Tool
dag.h
Go to the documentation of this file.
1 
20 #ifndef DAG_H
21 #define DAG_H
22 
23 #include <linalg.h>
24 #include <list.h>
25 
26 namespace spams
27 {
28 
29 template <typename T>
30 int count_cc_graph(const SpMatrix<T>& G, Vector<T>& active) {
31  int nzmax=0;
32  for (int i = 0; i<active.n(); ++i)
33  if (active[i]) nzmax++;
34  list_int** cc = new list_int*[nzmax];
35  int* pr_list = new int[active.n()];
36  memset(pr_list,-1,active.n()*sizeof(int));
37  int count=0;
38  list_int list;
39  for (int i = 0; i<active.n(); ++i) {
40  if (active[i]) {
41  list.push_back(i);
42  pr_list[i]=count;
43  cc[count] = new list_int();
44  cc[count++]->push_back(i);
45  }
46  }
47  const int* pB = G.pB();
48  const int* r = G.r();
49  const T* v = G.v();
50 
51  while (!list.empty()) {
52  int node=list.front();
53  list.pop_front();
54  for (int j = pB[node]; j<pB[node+1]; ++j) {
55  int child=r[j];
56  if (active[child]) {
57  if (pr_list[node] != pr_list[child]) {
58  // fusion
59  const int pr = pr_list[child];
60  for (const_iterator_int it = cc[pr]->begin(); it != cc[pr]->end(); ++it) {
61  pr_list[*it]=pr_list[node];
62  }
63  cc[pr_list[node]]->fusion(*(cc[pr]));
64  delete(cc[pr]);
65  cc[pr]=NULL;
66  }
67  }
68  }
69  }
70 
71  int num_cc=0;
72  for (int i = 0; i<nzmax; ++i) if (cc[i]) num_cc++;
73  for (int i = 0; i<nzmax; ++i) delete(cc[i]);
74  delete[](cc);
75  delete[](pr_list);
76  return num_cc;
77 }
78 
79 template <typename T>
80 void remove_cycles(const SpMatrix<T>& G1, SpMatrix<T>& G2) {
81  G2.copy(G1);
82  int n = G1.n();
83  int* color = new int[n];
84  memset(color,0,n*sizeof(int));
85  int next=0;
86  list_int list;
87  int* pB = G2.pB();
88  int* r = G2.r();
89  T* v = G2.v();
90 
91  list_int current_path;
92  bool cycle_detected=true;
93  int cycles_removed=0;
94  while (cycle_detected) {
95  cycle_detected=false;
96  list.clear();
97  for (int i = 0; i<n; ++i) if (color[i] != 2) {
98  color[i]=0;
99  list.push_back(i);
100  }
101  current_path.clear();
102  while (!list.empty()) {
103  const int node=list.front();
104  if (color[node] == 0) {
105  current_path.push_front(node);
106  color[node]=1;
107  for (int i = pB[node]; i<pB[node+1]; ++i) {
108  if (v[i]) {
109  const int child = r[i];
110  if (color[child]==1) {
111  cycle_detected=true;
112  list_int reverse_path;
113  current_path.reverse(reverse_path);
114  while (true) {
115  const int current_node=reverse_path.front();
116  if (current_node == child) break;
117  reverse_path.pop_front();
118  }
119  // remove beginning of the path
120 
121  reverse_path.push_back(child);
122  T min_link = INFINITY;
123  int min_node= -1;
124  const_iterator_int it = reverse_path.begin();
125  // detect weakest link
126  while (true) {
127  const int current_node=*it;
128  ++it;
129  if (it == reverse_path.end()) break;
130  for (int j = pB[current_node]; j<pB[current_node+1]; ++j) {
131  if (r[j]==*it) {
132  if (min_link > v[j]) {
133  min_link=v[j];
134  min_node=j;
135  }
136  break;
137  }
138  }
139  }
140  v[min_node]=0;
141  list.clear();
142  cycles_removed++;
143  cerr << " \r" << cycles_removed;
144  break;
145  } else if (color[child]==0) {
146  list.push_front(child);
147  }
148  }
149  }
150  } else if (color[node] == 1) {
152  color[node]=2;
153  list.pop_front();
154  current_path.pop_front();
155  } else if (color[node] == 2) {
156  list.pop_front();
157  }
158  }
159  }
160  int current_remove=0;
161  for (int i = 0; i<n; ++i) {
162  int old_remove=current_remove;
163  for (int j = pB[i]; j < pB[i+1]; ++j) {
164  if (v[j]) {
165  r[j-current_remove]=r[j];
166  v[j-current_remove]=v[j];
167  } else {
168  current_remove++;
169  }
170  }
171  pB[i]-=old_remove;
172  }
173  pB[n]-=current_remove;
174 
175  delete[](color);
176 }
177 
178 template <typename T>
180  const int n = G.n();
181  T* num_paths = new T[n];
182  memset(num_paths,0,n*sizeof(T));
183  const int* pB = G.pB();
184  const int* r = G.r();
185  int* color = new int[n];
186  memset(color,0,n*sizeof(int));
187  list_int list;
188  for (int i = 0; i<n; ++i)
189  list.push_back(i);
190 
191  while (!list.empty()) {
192  const int node=list.front();
193  if (color[node]==0) {
194  for (int i = pB[node]; i<pB[node+1]; ++i) {
195  list.push_front(r[i]);
196  }
197  color[node]++;
198  } else if (color[node]==1) {
199  num_paths[node]=T(1.0);
200  for (int i = pB[node]; i<pB[node+1]; ++i) {
201  num_paths[node]+=num_paths[r[i]];
202  }
203  color[node]++;
204  } else if (color[node]==2) {
205  list.pop_front();
206  }
207  }
208 
209  T sum=cblas_asum<T>(n,num_paths,1);
210  delete[](num_paths);
211  delete[](color);
212  return sum;
213 }
214 
215 }
216 
217 #endif
218 
Sparse Matrix class.
Definition: linalg.h:63
List< int > list_int
Definition: list.h:146
void pop_front()
Definition: list.h:62
Definition: dag.h:26
T count_paths_dags(const SpMatrix< T > &G)
Definition: dag.h:179
int pB(const int i) const
returns pB[i]
Definition: linalg.h:779
void reverse(List< T > &list)
Definition: list.h:124
bool empty() const
Definition: list.h:58
void push_back(T elem)
Definition: list.h:69
ListIterator< T > & begin() const
Definition: list.h:117
int n() const
returns the size of the vector
Definition: linalg.h:591
void push_front(T elem)
Definition: list.h:80
int count_cc_graph(const SpMatrix< T > &G, Vector< T > &active)
Definition: dag.h:30
T v(const int i) const
returns v[i]
Definition: linalg.h:783
Dense Vector class.
Definition: linalg.h:65
void copy(const SpMatrix< T > &mat)
Definition: linalg.h:4102
int r(const int i) const
returns r[i]
Definition: linalg.h:781
void clear()
Definition: list.h:87
T front() const
Definition: list.h:59
void * end() const
Definition: list.h:118
int n() const
returns the number of rows
Definition: linalg.h:787
void fusion(const List< T > &list)
Definition: list.h:119
void remove_cycles(const SpMatrix< T > &G1, SpMatrix< T > &G2)
Definition: dag.h:80
#define INFINITY
Definition: utils.h:62