1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
use crate::debug::debug;
use crate::traits::{Int, Scalar};
use crate::Matrix;
// Depth-first-search workspace.
pub struct DFS {
root_list: Vec<usize>,
// ptr_list: Vec<usize>,
ptr_list2: Vec<(usize, usize)>,
flag: Vec<bool>,
}
impl DFS {
pub fn new(n: usize) -> Self {
Self {
root_list: vec![0; n],
// ptr_list: vec![0; n],
ptr_list2: vec![(0, 0); n],
flag: vec![false; n],
}
}
pub fn ludfs<I: Int, S: Scalar>(
&mut self,
l_mat: &Matrix<I, S>,
b_rowidx: &[I],
rperm: &[Option<usize>],
) -> &[usize] {
debug!("b = {:?}", b_rowidx);
// let csgraph = matrix_to_csc(l_mat);
// debug!("L =\n{}", csgraph.to_table());
let n = l_mat.len();
// let n = csgraph.cols();
// let indices = csgraph.rowidx();
// let indptr = csgraph.colptr();
// let n = l_mat.cols();
// let indices = l_mat.rowidx();
// let indptr = l_mat.colptr();
// println!("\nindices = {:?}", indices);
// println!("indptr = {:?}", indptr);
let mut i_rl_start = n;
for &e0 in b_rowidx {
// The depth-first search must mark the vertices it
// has reached, to avoid repeating parts of the search.
if self.flag[e0.to_index()] {
continue;
}
// self.dfs(e0, indices, indptr, &mut i_rl_start, rperm);
self.dfs(
e0.to_index(),
/*indices, indptr,*/ l_mat,
&mut i_rl_start,
rperm,
);
}
let found = &self.root_list[i_rl_start..];
debug!(
"found = {:?} {} {:?}",
found.to_vec(),
i_rl_start,
self.root_list
);
// debug!("flag = {:?}", self.flag);
found.iter().for_each(|i| self.flag[*i] = false);
found
}
// Based on `depth_first_directed` from SciPy v1.11.
// Modified to support repeated calls with different start nodes
// and the same permuted input matrix.
//
// Author: Jake Vanderplas -- <vanderplas@astro.washington.edu>
// License: BSD, (C) 2012
pub(crate) fn dfs<I: Int, S: Scalar>(
&mut self,
head_node: usize,
// indices: &[usize],
// indptr: &[usize],
l_mat: &Matrix<I, S>,
i_rl_start: &mut usize,
rperm: &[Option<usize>],
) {
// let n = node_list.len();
// node_list[0] = head_node;
self.root_list[0] = head_node;
let mut i_root_opt: Option<usize> = Some(0);
// let mut i_nl_end = 1;
// flag[head_node] = true;
while let Some(mut i_root) = i_root_opt {
let pnode = self.root_list[i_root];
let pnode_p = rperm[pnode];
if !self.flag[pnode] {
self.flag[pnode] = true;
match pnode_p {
Some(pnode_p) => {
// self.ptr_list[i_root] = indptr[pnode_p];
self.ptr_list2[i_root] = (pnode_p, 0);
}
None => {
// self.ptr_list[i_root] = 0;
self.ptr_list2[i_root] = (0, 0);
}
}
}
let mut no_children = true;
// let indptr1 = self.ptr_list[i_root];
// let indptr2 = match pnode_p {
// Some(pnode_p) => indptr[pnode_p + 1],
// None => 0,
// };
// println!(
// "\nind[{}..{}] = {:?}",
// indptr1,
// indptr2,
// indices[indptr1..indptr2].to_vec()
// );
if pnode_p.is_some() {
let (lcolind, lrow_offset) = self.ptr_list2[i_root];
let lcol = &l_mat[lcolind];
// println!(
// "col = {} {} ({}) - {:?}",
// lcolind,
// lrow_offset,
// pnode_p.is_none(),
// lcol
// );
// debug!(
// "pnode = {}, pnode_p = {:?}, p1 = {}, p2 = {}",
// pnode, pnode_p, indptr1, indptr2
// );
let mut k = 0;
// for i in indptr1..indptr2 {
for j in lrow_offset..lcol.len() {
k += 1;
let (cnode, _) = &lcol[j];
// let cnode = indices[i];
// let cnode = l_mat[indptr1.0][j].0;
if self.flag[cnode.to_index()] {
continue;
} else {
// self.ptr_list[i_root] = i;
self.ptr_list2[i_root] = (lcolind, k - 1);
i_root += 1;
i_root_opt = Some(i_root);
self.root_list[i_root] = cnode.to_index();
// node_list[i_nl_end] = cnode;
// flag[cnode] = true;
// i_nl_end += 1;
// debug!("i_root = {}, cnode = {}", i_root, cnode);
no_children = false;
break;
}
}
}
if *i_rl_start == 0 {
break;
}
if no_children {
i_root_opt = if i_root > 0 { Some(i_root - 1) } else { None };
*i_rl_start -= 1;
self.root_list[*i_rl_start] = pnode;
// debug!("i_rl_start = {}, pnode = {}", *i_rl_start, pnode);
}
}
}
}