Rust で階層的クラスタリング
こちらの続き。
階層的クラスタリングについてはこちらを。サンプルデータもリンク先のものを利用させていただく。
再帰的なデータ構造の作成
階層的クラスタリングでは もっとも距離の小さいクラスタ同士を順に併合していく。そのため併合元となるクラスタ (子) と 作成される新しいクラスタ (親) に親子関係をもたせたい。
こういったとき、親から子に参照を持たせればいいかな? と思うのだが、Rust でこれをやろうとすると "親のインスタンスが存在している間は子への参照が有効となっている" ということを lifetime として記述しなければならず、すこし面倒である。
そのため、参照ではなく所有権を移してしまうのが一般的なようだ。
参考 The idea behind a recursive struct - help - The Rust Programming Language Forum
階層的クラスタリング
場合分けをシンプルにするため、ほぼ同一のロジックである最近隣法、最遠隣法、群平均法を書いた。
extern crate csv; extern crate nalgebra; extern crate num; extern crate rand; use csv::Reader; use nalgebra::{DVec, DMat, RowSlice, Iterable}; use num::{Signed, Zero, Float}; use std::f64; use std::io; use std::ops::Index; use std::str; use std::str::FromStr; fn main() { let data = "名前,算数,理科,国語,英語,社会 田中,89,90,67,46,50 佐藤,57,70,80,85,90 鈴木,80,90,35,40,50 本田,40,60,50,45,55 川端,78,85,45,55,60 吉野,55,65,80,75,85 斉藤,90,85,88,92,95"; // http://burntsushi.net/rustdoc/csv/ let mut reader = csv::Reader::from_string(data).has_headers(true); // http://nalgebra.org/doc/nalgebra/struct.DMat.html let dx = read_csv_f64(&mut reader); println!("最近隣法"); let mut hclust = HClust::new(ClusterDistance::Single); hclust.fit(&dx); println!("最遠隣法"); let mut hclust = HClust::new(ClusterDistance::Complete); hclust.fit(&dx); println!("群平均法"); let mut hclust = HClust::new(ClusterDistance::Average); hclust.fit(&dx); } enum ClusterDistance { Single, // 最近隣法 Complete, // 最遠隣法 Average, // 群平均法 } struct HClust { dist_mat: DMat<f64>, // 距離行列 method: ClusterDistance, // クラスタ間距離の計算方法 clusters: Vec<Cluster> } impl HClust { fn new(method: ClusterDistance) -> HClust { HClust { // dummy dist_mat: DMat::from_elem(1, 1, 1.), method: method, clusters: vec![] } } fn fit(&mut self, data: &DMat<f64>) { self.dist_mat = self.get_dist_matrix(&data); // クラスタを初期化 for i in 0..data.nrows() { let c = Cluster::from_nodes(vec![i]); self.clusters.push(c); } while self.clusters.len() > 1 { self.fit_step(); } } /// もっとも距離の小さいクラスタを併合 fn fit_step(&mut self) { let mut tmp_i = 0; let mut tmp_j = 0; let mut current_dist = f64::MAX; for i in 0..self.clusters.len() { for j in (i + 1)..self.clusters.len() { let d = self.get_cluster_dist(&self.clusters[i], &self.clusters[j]); if d < current_dist { current_dist = d; tmp_i = i; tmp_j = j; } } } // R と比較するための出力を表示 println!("Mergeing clusters: {:?} and {:?}: Distance {:?}", &self.clusters[tmp_i].nodes, &self.clusters[tmp_j].nodes, (current_dist * 1e+6).round() / 1e+6); let mut new_clusters: Vec<Cluster> = vec![]; for (i, n) in self.clusters.iter().enumerate() { if i != tmp_i && i != tmp_j { let n2 = Cluster::from_nodes(n.nodes.clone()); new_clusters.push(n2); } } // Vec から要素を取り出し & 所有権を取得 // コンストラクタに所有権ごと渡す let new = Cluster::from_clusters(self.clusters.swap_remove(tmp_j), self.clusters.swap_remove(tmp_i)); new_clusters.push(new); self.clusters = new_clusters; } /// クラスタに含まれるノード間の距離を計算 fn get_cluster_dist(&self, c1: &Cluster, c2: &Cluster) -> f64 { let mut dists: Vec<f64> = vec![]; for i in &c1.nodes { for j in &c2.nodes { dists.push(self.get_node_dist(*i, *j)); } } return self.fold_dist_vec(dists); } /// クラスタに含まれるノード間の距離から、クラスタ間の距離を計算 fn fold_dist_vec(&self, dists: Vec<f64>) -> f64 { match self.method { ClusterDistance::Single => { return dists.iter().fold(f64::MAX, |a, b| a.min(*b)); }, ClusterDistance::Complete => { return dists.iter().fold(f64::MIN, |a, b| a.max(*b)); }, ClusterDistance::Average => { return dists.iter().fold(0., |a, b| a + b) / (dists.len() as f64); }, } } /// 距離行列を作成 fn get_dist_matrix(&self, data: &DMat<f64>) -> DMat<f64> { // 各列は 0番目 から ノード数-1 番目までのノードに対応 - // 各行は 1番目 から ノード数番目 までのノードに対応 return DMat::from_fn(data.nrows() - 1, data.nrows() - 1, |i, j| if i >= j { euc_dist(&data.row_slice(i + 1, 0, data.ncols()), &data.row_slice(j, 0, data.ncols())) } else { 0. }); } /// 距離行列を利用して ノード間の距離を求める fn get_node_dist(&self, i: usize, j: usize) -> f64 { match i > j { true => *self.dist_mat.index((i - 1, j)), false => *self.dist_mat.index((j - 1, i)) } } } struct Cluster { nodes: Vec<usize>, children: Vec<Cluster> } impl Cluster { fn from_nodes(nodes: Vec<usize>) -> Cluster { Cluster { nodes: nodes, children: vec![] } } /// 二つのクラスタを併合 fn from_clusters(left: Cluster, right: Cluster) -> Cluster { let mut id = vec![]; for i in &left.nodes { id.push(*i); } for j in &right.nodes { id.push(*j); } Cluster { nodes: id, children: vec![left, right] } } } // ユークリッド距離 fn euc_dist<T: Float + Signed>(vec1: &DVec<T>, vec2: &DVec<T>) -> T { let mut val: T = Zero::zero(); for (v1, v2) in vec1.iter().zip(vec2.iter()) { val = val + num::pow(num::abs(*v1 - *v2), 2); } return val.sqrt(); } /// csv::Reder から f64 に変換できるカラムのみ読み込み pub fn read_csv_f64<R: io::Read>(reader: &mut Reader<R>) -> DMat<f64> { let mut x:Vec<f64> = vec![]; let mut nrows: usize = 0; for record in reader.byte_records().map(|r| r.unwrap()) { // f64 に変換できる列のみ読み込み for item in record.iter().map(|i| str::from_utf8(i).unwrap()) { match f64::from_str(item) { Ok(v) => x.push(v), Err(_) => {} }; } nrows += 1; } let ncols = x.len() / nrows; // http://nalgebra.org/doc/nalgebra/struct.DMat.html return DMat::from_row_vec(nrows, ncols, &x); }
出力。
# 最近隣法 # Mergeing clusters: [1] and [5]: Distance 12.409674 # Mergeing clusters: [2] and [4]: Distance 21.307276 # Mergeing clusters: [0] and [4, 2]: Distance 28.478062 # Mergeing clusters: [6] and [5, 1]: Distance 38.105118 # Mergeing clusters: [3] and [4, 2, 0]: Distance 47.106263 # Mergeing clusters: [5, 1, 6] and [4, 2, 0, 3]: Distance 54.313902 # 最遠隣法 # Mergeing clusters: [1] and [5]: Distance 12.409674 # Mergeing clusters: [2] and [4]: Distance 21.307276 # Mergeing clusters: [0] and [4, 2]: Distance 33.778692 # Mergeing clusters: [6] and [5, 1]: Distance 45.585085 # Mergeing clusters: [3] and [4, 2, 0]: Distance 60.133186 # Mergeing clusters: [5, 1, 6] and [4, 2, 0, 3]: Distance 91.531415 # 群平均法 # Mergeing clusters: [1] and [5]: Distance 12.409674 # Mergeing clusters: [2] and [4]: Distance 21.307276 # Mergeing clusters: [0] and [4, 2]: Distance 31.128377 # Mergeing clusters: [6] and [5, 1]: Distance 41.845102 # Mergeing clusters: [3] and [4, 2, 0]: Distance 53.305906 # Mergeing clusters: [5, 1, 6] and [4, 2, 0, 3]: Distance 69.922956
R の結果。コメントは ノードが 0 から始まるとして (上の出力にあわせて ) 記載。作成されるクラスタはどの方法でも同じだが、クラスタ間の距離 (3 列目) は異なる。
# 最近隣法 hc <- hclust(seiseki.d, method = "single") cbind(hc$merge, hc$height) # [,1] [,2] [,3] # [1,] -2 -6 12.40967 # ノード 1 と ノード 5 を併合 # [2,] -3 -5 21.30728 # ノード 2 と ノード 4 を併合 # [3,] -1 2 28.47806 # ノード 0 と ノード [2, 4] からなるクラスタを併合 # [4,] -7 1 38.10512 # ノード 6 と ノード [1, 5] からなるクラスタを併合 # [5,] -4 3 47.10626 # ノード 3 と ノード [0, 2, 4] からなるクラスタを併合 # [6,] 4 5 54.31390 # ノード [1, 5, 6] と ノード [0, 2, 3, 4] からなるクラスタを併合 # 最遠隣法 hc <- hclust(seiseki.d, method = "complete") cbind(hc$merge, hc$height) # [,1] [,2] [,3] # [1,] -2 -6 12.40967 # ノード 1 と ノード 5 を併合 # [2,] -3 -5 21.30728 # ノード 2 と ノード 4 を併合 # [3,] -1 2 33.77869 # ノード 0 と ノード [2, 4] からなるクラスタを併合 # [4,] -7 1 45.58509 # ノード 6 と ノード [1, 5] からなるクラスタを併合 # [5,] -4 3 60.13319 # ノード 3 と ノード [0, 2, 4] からなるクラスタを併合 # [6,] 4 5 91.53142 # ノード [1, 5, 6] と ノード [0, 2, 3, 4] からなるクラスタを併合 # 群平均法 hc <- hclust(seiseki.d, method = "average") cbind(hc$merge, hc$height) # [,1] [,2] [,3] # [1,] -2 -6 12.40967 # ノード 1 と ノード 5 を併合 # [2,] -3 -5 21.30728 # ノード 2 と ノード 4 を併合 # [3,] -1 2 31.12838 # ノード 0 と ノード [2, 4] からなるクラスタを併合 # [4,] -7 1 41.84510 # ノード 6 と ノード [1, 5] からなるクラスタを併合 # [5,] -4 3 53.30591 # ノード 3 と ノード [0, 2, 4] からなるクラスタを併合 # [6,] 4 5 69.92296 # ノード [1, 5, 6] と ノード [0, 2, 3, 4] からなるクラスタを併合
2016/12/20 追記 上記以外の手法を追加 / Lance-Williams formulation で書き直した。
