概要
学校で、lifegameのプログラムに独自の改良を施すという課題が出たため、Rust言語にてhashlifeアルゴリズムの実装を行った。lifegameは、二次元平面のセルの集合が各セルごとに隣接するセルの情報を元に時間発展する数理モデルだ。そして、Hashlifeアルゴリズムは、1984年にBill Gosperによって作成された、lifegameを実行するためのアルゴリズムだ。これを用いると、1.広大な空間を少ないメモリで表すことができ、2.時間をスキップしたシミュレーションを行うことができる。メモ化再帰と木構造を用いてこれを実現する。
詳細
lifegameをナイーブに実装すると、すべてのセルに対して演算をしないといけなかったために、範囲の一辺の長さを $n $とすると計算量が $O(n^2) $かかってしまう。アルゴリズムを改良することで、時間的、空間的に計算をスキップしながら計算可能となり、高速なシミュレーションが可能となる。
Hashlife構造体は、cache,top,generationの3つの要素を持つ。topには、作成する木の頂点が格納され、generationには今の世代、キャッシュには過去に利用したセルの進化の履歴が蓄積されていく。ここから、hashlifeアルゴリズムの解説をしていく。
四分木構造による二次元空間の表現
四分木のそれぞれの頂点、ノードは自分の深さの情報であるlevelと、その範囲に存在している生きているセルの数を表すpopulation, 自分自身よりひとつ深いノード群を表すchildren、そして子供の状態から計算されるhashの4つから構成されている。このノードは、深さに応じて大きさが変化する。具体的には、levelが0のノードは縦横それぞれ $2^0 $の長さで、levelが $n $のノードは、縦横がそれぞれ $2^n $の長さの範囲を表している。すなわち、levelが $n $であるノードの面積は $2^n \times 2^n = 2^{2n} $となり、そこのノード範囲の最大生存セル数は $2^{2n} $となる。
そして、childrenには、levelが一つ低いノードが4つ格納される。level0のノードは子供を持たない。レベルが1以上のノードは、縦と横の長さがそれぞれ $2^n $なので、levelが $n-1 $の、縦横の長さが $2^{n-1} $になった子供を4つずつもつ。これが四分木である所以だ。
ここで、lifegameは、隣り合うセルの情報のみしか用いることが出来ないため、情報が伝わる速度は最大でも、一世代に1セル分ずつということが分かる。この事実を用いて、ある正方形領域、ノードについて考察してみる。ここにLevelが $n $であるノードがあるとする。すると、このノードの外側にある領域の情報が今注目しているノードの内側へとどのような速度で侵食していくだろうか。注目しているノードを縦横4分割、合計16分割し、その中心部分の4つのエリアを新しいノードとして持ってくるとする。すると、そのエリアに、元ノードの外側の情報がたどり着くのは、 $2^n/4 $世代先になる。これが何を表しているかというと、levelがnのノードがあったとき、その内側のlevelが $n-1 $のノードの $2^n/4 $世代後までの状態は、既にそのノードのみで決定されているということだ。
以上の考えより、あるlevelが $n $のノードの内側のレベル $n-1 $のノードの $2^n/4 $世代先のノードの状態は既にそのノードの状態から確定しているということは納得していただけたと思う。ここから、その状態を如何にして高速に計算するかを説明する。ここがhashlifeアルゴリズムの肝となる部分だ。
ここで、ノナント(nonants)というものが登場する。これは、Levelが2以上のノードについて一対一で対応させることができるノード群のようなものだ。ノナントは、9つのノードの集合体である。実際には、Levelが3以上のノードについて、これを用いて考える。具体例として、8x8のレベル3のノードについて考えていく。ノナントは以下の画像のようにノードを16分割し、それを2x2の正方形領域9個を取り出すことで、レベル2のノード９個を得ることができる。
今、これらのノナントを構成するノード群はそれぞれ隣り合うノードと、 $2^n/4 $ずつ重なって構成されている。これより、これら各ノードに対して前述した、levelがnのノードがあったとき、その内側のlevelがn-1のノードの $2^n/4 $世代後までの状態は、既にそのノードのみで決定されているを適用すると、ノナントを構成するそれぞれのレベル2のノードに対して、Level1のノードの1世代後の状態が既に決定されていることとなる。
以上より、この6x6の領域に関して、1世代後の領域を得ることができたこととなる。続いて、この6x6の領域を更に9分割し、4x4のノードを4つ用意する。そして、同様にlevelがnのノードがあったとき、その内側のlevelがn-1のノードの $2^n/4 $世代後までの状態は、既にそのノードのみで決定されているを適用すると、新たに作成したレベル2のノードに対して、level1の内側ノードの更に1世代後の状態が既に決定されているということとなる。
以上の議論を一般化すると、Level nのノードがあった時に、Level n-1のノナントによるノード群に変換、 $n-2 $の $(\frac{3}{4}n \times \frac{3}{4}n) $の大きさの $2^{n-3} $世代後の状態を取得。更にそこから4つのLevel n-1のノードに変換し、それぞれに対して世代を進めることで、更に $2^{n-3} $世代進めることができる。それらを合成し、結果として、levelが $n-1 $の世代が $2^{n-2} $進んだノードを得ることができる。ということである。
さて、以上で仕組みは説明できた。それではどのようにしてそれを高速に実行するかということを考えていく。以上のものを純粋に実装すると、計算量がnに応じて指数関数的に上昇してしまう。それではどうするかと言うと、キャッシュを用いるのである。今回の私の実装では、キャッシュを3つ用意した。joinキャッシュと、そのノードから可能な最大数のステップをした時に利用できるstepキャッシュ、そして、指定した回数の世代分だけ進むときに利用できるstep_partialキャッシュである。それぞれのキャッシュについて、一つずつ説明していく。
キャッシュしてメモ化再帰
joinキャッシュは、４つのノードから計算されるハッシュ値をキーとしてノードへの参照を返すハッシュマップとなっている。これを利用するのは、新しいレイアウトのノードが現れた時である。キャッシュには、ノードがそれぞれ独立して格納されている。キャッシュ内の全てのノードは、それらの子孫全ての情報を保持している。このキャッシュを用いることで、同じ状態のノードが複数生成されないようになっている。すなわち、四分木は、木に見えて、同じ子供を共通して保持した木なのである。
stepキャッシュは、ノードLevel nをキーとして、それを $2^{n-2} $世代進んだノードLevel n-1への参照を返すハッシュマップとなっている。これによって、一度計算したノードに対するステップは、キャッシュから取り出すだけで計算が終了することとなる。
step_partialキャッシュは、stepキャッシュと同様に、ノードを受け取ってノードを返すが、キーとして、kを追加して、 $2^{k-2} $世代進めたノードへの参照を返却する。このキャッシュによって、大きなLifegameのフィールドに対しても刻んだ世代先の情報を得ることができるようになる。
以上がHashlifeの動作の概要である。詳細なプログラムは、コードhashlife.rsを参照されたい。
参考文献
コード
hashlife.rs
use std::collections::HashMap;
use std::collections::hash_map::DefaultHasher;
use std::hash::{Hash, Hasher};
use std::rc::Rc;
//Node->木の頂点。状態は、ノードとその子どもの４つの状態木によって表現される。
//ノードには深さがある。level。0は葉で、それ以降は分岐となる。
//各ノードには密度がある。どれくらいのAutomataが生きているかを表している。ノードが葉の時、１か０をとる。
//ハッシュライフのアルゴリズムにあるノードは四分木として表される。Level０はNoneとなる。
#[derive(Debug)]
pub struct Node {
level: usize,
population: usize,
hash: u64,
children: Option<Children>,
}
impl PartialEq for Node {
fn eq(&self, other: &Node) -> bool {
self.level == other.level && self.population == other.population && self.hash == other.hash
}
}
impl Eq for Node {}
#[derive(Clone)]
pub struct Viewport {
pub center_x: f64,
pub center_y: f64,
/// ズームレベル
/// zoom_k > 0: ズームイン（1セルが2^(zoom_k+1)ピクセル）
/// zoom_k = 0: 1セル = 2ピクセル
/// zoom_k < 0: ズームアウト（2ピクセルに2^(-zoom_k)セル）
pub zoom_k: i32,
pub width_px: usize,
pub height_px: usize,
}
impl Viewport {
pub fn new(width_px: usize, height_px: usize, zoom_k: i32) -> Self {
Self {
center_x: 0.0,
center_y: 0.0,
zoom_k,
width_px,
height_px,
}
}
pub fn tile_size(&self) -> usize {
if self.zoom_k >= 0 {
2 << self.zoom_k as usize // 2, 4, 8, 16, ...
} else {
2 // ズームアウト時は2px固定
}
}
/// 1タイルがカバーするセル数（1方向）
pub fn cells_per_tile(&self) -> usize {
if self.zoom_k >= 0 {
1
} else {
1 << (-self.zoom_k) as usize // 2, 4, 8, ...
}
}
/// 描画に使用するNodeのレベル
pub fn target_level(&self) -> usize {
if self.zoom_k >= 0 {
0
} else {
(-self.zoom_k) as usize
}
}
/// タイル数（X方向）
pub fn tiles_x(&self) -> usize {
self.width_px / self.tile_size()
}
/// タイル数（Y方向）
pub fn tiles_y(&self) -> usize {
self.height_px / self.tile_size()
}
/// タイル座標からワールド座標へ変換
pub fn tile_to_world(&self, tile_x: usize, tile_y: usize) -> (isize, isize) {
let cells_per_tile = self.cells_per_tile() as isize;
let half_tiles_x = (self.tiles_x() / 2) as isize;
let half_tiles_y = (self.tiles_y() / 2) as isize;
let world_x = self.center_x as isize + (tile_x as isize - half_tiles_x) * cells_per_tile;
let world_y = self.center_y as isize + (half_tiles_y - tile_y as isize) * cells_per_tile;
(world_x, world_y)
}
/// ズームイン
pub fn zoom_in(&mut self) {
if self.zoom_k < 6 {
// 最大64pxまで
self.zoom_k += 1;
}
}
/// ズームアウト
pub fn zoom_out(&mut self) {
if self.zoom_k > -10 {
// 最小2^10=1024セル/タイルまで
self.zoom_k -= 1;
}
}
/// パン（移動）
pub fn pan(&mut self, dx: f64, dy: f64) {
let scale = self.cells_per_tile() as f64;
self.center_x += dx * scale;
self.center_y += dy * scale;
}
/// スクリーン座標（ピクセル）からワールド座標に変換
/// tile_to_worldの逆変換
pub fn screen_to_world(&self, screen_x: f32, screen_y: f32) -> (isize, isize) {
let tile_size = self.tile_size() as f32;
let cells_per_tile = self.cells_per_tile() as f32;
let half_tiles_x = (self.tiles_x() / 2) as f32;
let half_tiles_y = (self.tiles_y() / 2) as f32;
// スクリーン座標からタイル座標を計算
let tile_x = screen_x / tile_size;
let tile_y = screen_y / tile_size;
// タイル座標からワールド座標を計算（tile_to_worldの逆）
let world_x = self.center_x as f32 + (tile_x - half_tiles_x) * cells_per_tile;
let world_y = self.center_y as f32 + (half_tiles_y - tile_y) * cells_per_tile;
(world_x as isize, world_y as isize)
}
/// ワールド座標からスクリーン座標（ピクセル）に変換
/// tile_to_worldの逆で、タイルの左上端の座標を返す
pub fn world_to_screen(&self, world_x: isize, world_y: isize) -> (f32, f32) {
let tile_size = self.tile_size() as f32;
let cells_per_tile = self.cells_per_tile() as f32;
let half_tiles_x = (self.tiles_x() / 2) as f32;
let half_tiles_y = (self.tiles_y() / 2) as f32;
// ワールド座標からタイル座標を計算（tile_to_worldの逆）
let tile_x = (world_x as f32 - self.center_x as f32) / cells_per_tile + half_tiles_x;
let tile_y = half_tiles_y - (world_y as f32 - self.center_y as f32) / cells_per_tile;
// タイル座標からスクリーン座標を計算
let screen_x = tile_x * tile_size;
let screen_y = tile_y * tile_size;
(screen_x, screen_y)
}
/// 表示範囲のワールド座標を取得 (min_x, max_x, min_y, max_y)
pub fn visible_world_range(&self) -> (isize, isize, isize, isize) {
let (min_x, max_y) = self.screen_to_world(0.0, 0.0);
let (max_x, min_y) = self.screen_to_world(self.width_px as f32, self.height_px as f32);
(min_x, max_x, min_y, max_y)
}
}
#[derive(Hash, Debug)]
struct Children {
nw: Rc<Node>,
ne: Rc<Node>,
sw: Rc<Node>,
se: Rc<Node>,
}
struct GrandChildren {
nwnw: Rc<Node>,
nwne: Rc<Node>,
nwsw: Rc<Node>,
nwse: Rc<Node>,
nenw: Rc<Node>,
nene: Rc<Node>,
nesw: Rc<Node>,
nese: Rc<Node>,
swnw: Rc<Node>,
swne: Rc<Node>,
swsw: Rc<Node>,
swse: Rc<Node>,
senw: Rc<Node>,
sene: Rc<Node>,
sesw: Rc<Node>,
sese: Rc<Node>,
}
struct GrandAutomata {
nwnw: Automata,
nwne: Automata,
nwsw: Automata,
nwse: Automata,
nenw: Automata,
nene: Automata,
nesw: Automata,
nese: Automata,
swnw: Automata,
swne: Automata,
swsw: Automata,
swse: Automata,
senw: Automata,
sene: Automata,
sesw: Automata,
sese: Automata,
}
struct Nonants {
nw: Rc<Node>,
ne: Rc<Node>,
sw: Rc<Node>,
se: Rc<Node>,
n_: Rc<Node>,
s_: Rc<Node>,
w_: Rc<Node>,
e_: Rc<Node>,
c_: Rc<Node>,
}
pub struct Hashlife {
cache: Cache,
top: Option<Rc<Node>>,
generation: usize,
}
struct Cache {
join: HashMap<u64, Rc<Node>>,
step: HashMap<Rc<Node>, Rc<Node>>,
step_partial: HashMap<(u64, usize), Rc<Node>>,
dead: Option<Rc<Node>>,
alive: Option<Rc<Node>>,
}
impl Cache {
fn new() -> Self {
Self {
join: HashMap::new(),
step: HashMap::new(),
step_partial: HashMap::new(),
dead: None,
alive: None,
}
}
}
#[repr(u8)]
#[derive(Copy, Clone, Debug, Eq, PartialEq, Hash)]
pub enum Automata {
Dead = 0,
Alive = 1,
}
impl Automata {
pub fn from(number: usize) -> Self {
match number {
0 => Self::Dead,
1 => Self::Alive,
_ => panic!("not an automata"),
}
}
}
impl Default for Hashlife {
fn default() -> Self {
Self::new()
}
}
impl Hashlife {
pub fn new() -> Self {
Self {
cache: Cache::new(),
top: None,
generation: 0,
}
}
//４つのノードを受け取って、その上の一つ深い階層のノード参照を返す。ex.4x4のノード4つを受け取って、8x8のノードを返す。
fn join(&mut self, nw: Rc<Node>, ne: Rc<Node>, sw: Rc<Node>, se: Rc<Node>) -> Rc<Node> {
let children = Children::from(&nw, &ne, &sw, &se);
let level = nw.level + 1;
assert_eq!(nw.level, ne.level);
assert_eq!(nw.level, sw.level);
assert_eq!(nw.level, se.level);
let population = nw.population + ne.population + sw.population + se.population;
let hash = calculate_hash(&children);
if let Some(ref_to_node) = self.cache.join.get(&hash) {
return Rc::clone(ref_to_node);
}
let children = Some(children);
let node = Node {
level,
population,
hash,
children,
};
let node = Rc::new(node);
self.cache.join.insert(hash, Rc::clone(&node));
node
}
//ノード中心部分をノードとして取得
fn center(&mut self, node: Rc<Node>) -> Rc<Node> {
let c = node.get_children();
self.join(
Rc::clone(&c.nw.get_children().se),
Rc::clone(&c.ne.get_children().sw),
Rc::clone(&c.sw.get_children().ne),
Rc::clone(&c.se.get_children().nw),
)
}
//ノードを受け取って、そのノナントを返す。ノナントとは、次の世代を計算するための9個あるノード軍。
// 一つのノードからは、そのノードの1/4の範囲しか次のノードを計算できない。
// 9個のノードを組み合わせることで、9/16の範囲の次世代、2^hoge世代先の状態を計算できるようになる。
// ex.8x8のノードを受け取って、その左上4x4,上4x4,右上4x4、。。。と9このノードを組み合わせたノナントを返す。
fn make_nonants(&mut self, node: Rc<Node>) -> Nonants {
match &node.level {
0 => panic!("attempted to bread node into 9x9 at level 0"),
1 => panic!("attempted to bread node into 9x9 at level 1"),
2 => panic!("attempted to bread node into 9x9 at level 2"),
_ => (),
};
let c = node.get_children();
let g = node.get_grand_children();
Nonants {
nw: Rc::clone(&c.nw),
ne: Rc::clone(&c.ne),
sw: Rc::clone(&c.sw),
se: Rc::clone(&c.se),
n_: self.join(g.nwne, g.nenw, Rc::clone(&g.nwse), Rc::clone(&g.nesw)),
e_: self.join(Rc::clone(&g.nesw), g.nese, Rc::clone(&g.senw), g.sene),
s_: self.join(Rc::clone(&g.swne), Rc::clone(&g.senw), g.swse, g.sesw),
w_: self.join(g.nwsw, Rc::clone(&g.nwse), g.swnw, Rc::clone(&g.swne)),
c_: self.join(g.nwse, g.nesw, g.swne, g.senw),
}
}
//引数k: 2^k世代進める。
//ノードを一つ受け取って、ステップを進める。一つレベルの低い中心部分ノードを返す。一回の再帰で2ステップ進めることで、2^kステップを進めるようになっている。
//levelkのノードを受け取ると、2^(k-2)世代進めたlevelk-1の中心部分のノードを返すことができる。
//ex1.Level2のノードを受け取った時、内側の2x2ののノードを1ステップ進めて返す。
//ex2.level3のノード(8x8)を受け取った時、4x4のノナント9個を生成し、それぞれのノードを再帰でステップさせて2x2のノード9個、6x6の1ステップ進んだ領域を得る。
//ここからさらに2x2のノードを左上右上左下右下で4x4のノード4つにしてstepする。すると、初期状態から2世代進んだ中心部分4x4ノードが得られる。
//ex3.level4のノード(16x16)を受け取った時、8x8のノナント9個を生成し、それぞれのノードを再帰でステップさせて4x4のノード9個、12x12の2ステップ進んだ領域を得る。
//ここから更に4x4のノードを左上右上左下右下で8x8のノード4つにしてstepする。すると、初期状態から4世代進んだ中心部分4x4ノードが得られる。
fn step(&mut self, node: Rc<Node>) -> Rc<Node> {
if let Some(ref_to_node) = self.cache.step.get(&node) {
return Rc::clone(ref_to_node);
}
let step = match &node.level {
0 => panic!("atttempted to step a node with level 0"),
1 => panic!("atttempted to step a node with level 1"),
2 => {
let g = node.get_grand_automata();
let nw = simb3s23(
g.nwse, g.nwnw, g.nwne, g.nenw, g.nesw, g.senw, g.swne, g.swnw, g.nwsw,
);
let ne = simb3s23(
g.nesw, g.nwne, g.nenw, g.nene, g.nese, g.sene, g.senw, g.swne, g.nwse,
);
let sw = simb3s23(
g.swne, g.nwsw, g.nwse, g.nesw, g.senw, g.sesw, g.swse, g.swsw, g.swnw,
);
let se = simb3s23(
g.senw, g.nwse, g.nesw, g.nese, g.sene, g.sese, g.sesw, g.swse, g.swne,
);
let nw = self.make_automata(nw);
let ne = self.make_automata(ne);
let sw = self.make_automata(sw);
let se = self.make_automata(se);
self.join(nw, ne, sw, se)
}
_ => {
let mut g9x9 = self.make_nonants(Rc::clone(&node));
g9x9.nw = self.step(g9x9.nw);
g9x9.ne = self.step(g9x9.ne);
g9x9.sw = self.step(g9x9.sw);
g9x9.se = self.step(g9x9.se);
g9x9.n_ = self.step(g9x9.n_);
g9x9.e_ = self.step(g9x9.e_);
g9x9.s_ = self.step(g9x9.s_);
g9x9.w_ = self.step(g9x9.w_);
g9x9.c_ = self.step(g9x9.c_);
let mut res_nw = self.join(
Rc::clone(&g9x9.nw),
Rc::clone(&g9x9.n_),
Rc::clone(&g9x9.w_),
Rc::clone(&g9x9.c_),
);
let mut res_ne = self.join(
Rc::clone(&g9x9.n_),
Rc::clone(&g9x9.ne),
Rc::clone(&g9x9.c_),
Rc::clone(&g9x9.e_),
);
let mut res_sw = self.join(
Rc::clone(&g9x9.w_),
Rc::clone(&g9x9.c_),
Rc::clone(&g9x9.sw),
Rc::clone(&g9x9.s_),
);
let mut res_se = self.join(
Rc::clone(&g9x9.c_),
Rc::clone(&g9x9.e_),
Rc::clone(&g9x9.s_),
Rc::clone(&g9x9.se),
);
res_nw = self.step(res_nw);
res_ne = self.step(res_ne);
res_sw = self.step(res_sw);
res_se = self.step(res_se);
self.join(res_nw, res_ne, res_sw, res_se)
//self.join_nonants(g9x9)
}
};
self.cache.step.insert(node, Rc::clone(&step));
step
}
// 2^(k-2)ステップ進める関数。与えられるノードのレベルはk+2より大きいことを保証してもらう
// ex.level4,k=2なら1世代進める。
fn step_partial(&mut self, node: Rc<Node>, k: usize) -> Rc<Node> {
if node.level == k {
return self.step(node);
}
let cache_key = (node.hash, k);
if let Some(result) = self.cache.step_partial.get(&cache_key) {
return Rc::clone(result);
}
let mut g9x9 = self.make_nonants(Rc::clone(&node));
g9x9.nw = self.center(g9x9.nw);
g9x9.n_ = self.center(g9x9.n_);
g9x9.ne = self.center(g9x9.ne);
g9x9.w_ = self.center(g9x9.w_);
g9x9.c_ = self.center(g9x9.c_);
g9x9.e_ = self.center(g9x9.e_);
g9x9.sw = self.center(g9x9.sw);
g9x9.s_ = self.center(g9x9.s_);
g9x9.se = self.center(g9x9.se);
let mut res_nw = self.join(
Rc::clone(&g9x9.nw),
Rc::clone(&g9x9.n_),
Rc::clone(&g9x9.w_),
Rc::clone(&g9x9.c_),
);
let mut res_ne = self.join(
Rc::clone(&g9x9.n_),
Rc::clone(&g9x9.ne),
Rc::clone(&g9x9.c_),
Rc::clone(&g9x9.e_),
);
let mut res_sw = self.join(
Rc::clone(&g9x9.w_),
Rc::clone(&g9x9.c_),
Rc::clone(&g9x9.sw),
Rc::clone(&g9x9.s_),
);
let mut res_se = self.join(
Rc::clone(&g9x9.c_),
Rc::clone(&g9x9.e_),
Rc::clone(&g9x9.s_),
Rc::clone(&g9x9.se),
);
res_nw = self.step_partial(res_nw, k);
res_ne = self.step_partial(res_ne, k);
res_sw = self.step_partial(res_sw, k);
res_se = self.step_partial(res_se, k);
let result = self.join(res_nw, res_ne, res_sw, res_se);
self
.cache
.step_partial
.insert(cache_key, Rc::clone(&result));
result
}
//自分の領域を中心部分として、外側に空の領域を作る。
fn expand_empty_border(&mut self, node: Rc<Node>) -> Rc<Node> {
let c = node.get_children();
let e = self.empty(node.level - 1);
let e = || Rc::clone(&e);
let nw = self.join(e(), e(), e(), Rc::clone(&c.nw));
let ne = self.join(e(), e(), Rc::clone(&c.ne), e());
let sw = self.join(e(), Rc::clone(&c.sw), e(), e());
let se = self.join(Rc::clone(&c.se), e(), e(), e());
self.join(nw, ne, sw, se)
}
//2^k次の世代に進める。
//levelNのノードは、2なら1世代、3なら2世代、Nなら2^(N-2)世代進めることができる。
//topノードのレベルNがN<k+2の時、N==k+2になるまで拡張する必要がある。
fn next_2k_generation(&mut self, k: usize) {
let top = if let Some(top) = &self.top {
Rc::clone(top)
} else {
return;
};
let next = {
let mut expanded = self.expand_empty_border(Rc::clone(&top));
expanded = self.expand_empty_border(expanded);
while expanded.level < k + 2 {
expanded = self.expand_empty_border(expanded);
}
let step = self.step_partial(expanded, k + 2);
let g = step.get_grand_children();
let boarder_population = step.population
- g.nwse.population
- g.nesw.population
- g.swne.population
- g.senw.population;
if boarder_population == 0 {
self.join(g.nwse, g.nesw, g.swne, g.senw)
} else {
step
}
};
self.top = Some(next);
}
pub fn next_n_generation(&mut self, n: usize) {
let mut remaining = n;
let mut k = 0;
while remaining > 0 {
if remaining & 1 == 1 {
self.next_2k_generation(k);
}
remaining >>= 1;
k += 1;
}
self.generation += n;
}
//次の世代に進める
pub fn next_generation(&mut self) {
self.next_2k_generation(0);
self.generation += 1;
}
//オートマタのキャッシュを作るか、その生存状態のノード参照を返す
fn make_automata(&mut self, a: Automata) -> Rc<Node> {
match a {
Automata::Dead => {
if let Some(ref_to_node) = &self.cache.dead {
Rc::clone(ref_to_node)
} else {
let mut state = DefaultHasher::new();
a.hash(&mut state);
let node = Rc::new(Node {
level: 0,
population: a as usize,
children: None,
hash: state.finish(),
});
self.cache.dead = Some(Rc::clone(&node));
node
}
}
Automata::Alive => {
if let Some(ref_to_node) = &self.cache.alive {
Rc::clone(ref_to_node)
} else {
let mut state = DefaultHasher::new();
a.hash(&mut state);
let node = Rc::new(Node {
level: 0,
population: a as usize,
children: None,
hash: state.finish(),
});
self.cache.alive = Some(Rc::clone(&node));
node
}
}
}
}
//座標リストからHashlifeを構築する。座標リストを４分割し、空の座標リストがあった時は、
//そのリストはスキップすることで高速化している。
//計算量は一辺の長さをnとした時、最悪でO(4^log2(n))=O(n^2)となってしまう。
//キャッシュを用いて同じ区間パターンに同じノードの参照させるみたいな実装で、計算量を減らせるか。
//動作を説明する。
// ex1. pointsが[(0,0)]のときを見ていく。width,heightはどちらも1となる。レベルは1が設定される。
// root_x,root_yは0となり、construct_optimizedにはpointsと今から見たい範囲の左下の座標である(0,0)とレベルの1が渡される。
// レベル1のとき、Nodeの大きさは縦2横2である。再帰関数で、区間を上下左右に四分割し、1,1,1,1マスに対して処理し、結果を返す。
// ex2. pointsが[(0,0),(100,100)]のときを見ていく。min,maxは0,100となる、width,heightは101,101となる。levelはlog101をceilして7となる。
// レベルが7なので、Nodeの大きさは縦128横128である。これを、上下左右分割して再帰にかけていく
// 続きは下に記載する
pub fn from_points(points: &[(isize, isize)]) -> Self {
let mut hashlife = Hashlife::new();
if points.is_empty() {
hashlife.top = Some(hashlife.make_automata(Automata::Dead));
return hashlife;
}
let max_abs_x = points.iter().map(|(x, _)| x.abs()).max().unwrap();
let max_abs_y = points.iter().map(|(_, y)| y.abs()).max().unwrap();
let max_abs = std::cmp::max(max_abs_x, max_abs_y) as f64;
let level = if max_abs < 1.0 {
1
} else {
(max_abs + 1.0).log2().ceil() as usize + 1
};
let half_size = 1_isize << (level - 1);
let root_x = -half_size;
let root_y = -half_size;
let points_vec = points.to_vec();
let top_node = hashlife.construct(points_vec, root_x, root_y, level);
hashlife.top = Some(top_node);
hashlife
}
//ex2の続きを見ていく。[(0,0),(100,100)],0,0,7が入ってくる。
//half_sizeは1<<6で64,midx,midyは64となる。nw_pts,ne,sw,seはそれぞれ[],[(100,100)],[(0,0)],[]となる。
//これで、さらに四個の関数が呼ばれる。
// nw: [],0,64,6より、空なので、level6のemptyNodeへの参照が返される。
// ne: [(100,100)],64,64,6。(64,64)を起点としてlevel6(64x64)の4分割が行われる。また再帰される。
// sw: [(0,0)],0,0,6。(0,0)を起点としてlevel6の4分割再帰が行われる。
// se: [],64,0,6より、空なので、level6のemptyNodeへの参照が返される。
// これらの結果、空の部分は正方形で空のものが返される。
fn construct(
&mut self,
points: Vec<(isize, isize)>,
x: isize,
y: isize,
level: usize,
) -> Rc<Node> {
if points.is_empty() {
return self.empty(level);
}
if level == 0 {
return self.make_automata(Automata::Alive);
}
let half_size = 1_isize << (level - 1);
let mid_x = x + half_size;
let mid_y = y + half_size;
let mut nw_pts = Vec::new();
let mut ne_pts = Vec::new();
let mut sw_pts = Vec::new();
let mut se_pts = Vec::new();
for p in points {
let (px, py) = p;
if py >= mid_y {
if px < mid_x {
nw_pts.push(p);
} else {
ne_pts.push(p);
}
} else if px < mid_x {
sw_pts.push(p);
} else {
se_pts.push(p);
}
}
let nw = self.construct(nw_pts, x, mid_y, level - 1);
let ne = self.construct(ne_pts, mid_x, mid_y, level - 1);
let sw = self.construct(sw_pts, x, y, level - 1);
let se = self.construct(se_pts, mid_x, y, level - 1);
self.join(nw, ne, sw, se)
}
//任意深さの空ノードを返す関数
fn empty(&mut self, level: usize) -> Rc<Node> {
if level == 0 {
return self.make_automata(Automata::Dead);
}
let child = self.empty(level - 1);
let children = Children {
nw: Rc::clone(&child),
ne: Rc::clone(&child),
sw: Rc::clone(&child),
se: Rc::clone(&child),
};
let hash = calculate_hash(&children);
if let Some(ref_to_node) = self.cache.join.get(&hash) {
return Rc::clone(ref_to_node);
};
let empty = Rc::new(Node {
level,
population: 0,
children: Some(children),
hash,
});
self.cache.join.insert(hash, Rc::clone(&empty));
empty
}
//任意座標の生死状態を取得する関数
//ex. x=10,y=21の時を考える。topのレベルに応じてpositionsに値が格納される.level=5のときを考える。
//div_euclidの結果によって、positions=[(10,21),(5,10),(2,5),(1,2),(0,1)]となる。
//childrenには、level4,16x16が4こ入ったものが返ってくる。
//x,y>0より、0,0,positions,左上の子で補助関数が呼ばれる。
pub fn get(&self, x: isize, y: isize) -> Option<Automata> {
let top = if let Some(top) = self.top.as_ref() {
Rc::clone(top)
} else {
return None;
};
// ツリーの範囲をチェック
// 中心が(0,0)なので、範囲は -2^(L-1) から 2^(L-1) - 1
let half_size = 1_isize << (top.level - 1);
if x < -half_size || x >= half_size || y < -half_size || y >= half_size {
return None;
}
let mut positions = Vec::with_capacity(top.level);
let mut xx = x;
let mut yy = y;
for _ in 0..top.level {
positions.push((xx, yy));
xx = xx.div_euclid(2);
yy = yy.div_euclid(2);
}
let children = top.get_children();
if y < 0 {
if x < 0 {
Self::get_node_with(-1, -1, &positions, Rc::clone(&children.sw))
} else {
Self::get_node_with(0, -1, &positions, Rc::clone(&children.se))
}
} else if x < 0 {
Self::get_node_with(-1, 0, &positions, Rc::clone(&children.nw))
} else {
Self::get_node_with(0, 0, &positions, Rc::clone(&children.ne))
}
}
//getの補助関数
//上の続きを説明する。
// x,y>0より、0,0,positions,左上の子(level4)で補助関数が呼ばれた。
// [(10,21),(5,10),(2,5),(1,2),(0,1)]
//positionは(1,2)が入り、nw(0,1),ne(1,1),sw(0,0),se(1,0)となる。
//範囲外の場合はNoneを返す
fn get_node_with(
x: isize,
y: isize,
positions: &Vec<(isize, isize)>,
node: Rc<Node>,
) -> Option<Automata> {
if node.level == 0 {
return Some(node.as_automata());
}
let position = positions[node.level - 1];
let nw = (x * 2, y * 2 + 1);
let ne = (x * 2 + 1, y * 2 + 1);
let sw = (x * 2, y * 2);
let se = (x * 2 + 1, y * 2);
let children = node.get_children();
if position == nw {
Self::get_node_with(nw.0, nw.1, positions, Rc::clone(&children.nw))
} else if position == ne {
Self::get_node_with(ne.0, ne.1, positions, Rc::clone(&children.ne))
} else if position == sw {
Self::get_node_with(sw.0, sw.1, positions, Rc::clone(&children.sw))
} else if position == se {
Self::get_node_with(se.0, se.1, positions, Rc::clone(&children.se))
} else {
// 範囲外の座標の場合はNoneを返す（Dead扱い）
None
}
}
//現在の世代を返す関数
pub fn get_generation(&self) -> usize {
self.generation
}
/// topノードのレベルを返す
pub fn max_level(&self) -> usize {
self.top.as_ref().map(|t| t.level).unwrap_or(0)
}
/// 指定した座標とレベルでのNodeの密度を取得（0.0-1.0）
/// level=0の場合は0か1を返す
pub fn get_density_at(&self, x: isize, y: isize, level: usize) -> f32 {
let top = match self.top.as_ref() {
Some(t) => t,
None => return 0.0,
};
// levelがtopより大きい場合、topノードの密度を返す
if level >= top.level {
let max_cells = 1_usize << (2 * top.level);
return top.population as f32 / max_cells as f32;
}
if level == 0 {
// level 0の場合は既存のget関数を使う
match self.get(x, y) {
Some(Automata::Alive) => 1.0,
_ => 0.0,
}
} else {
let node = self.get_node_at(x, y, level);
match node {
Some(n) => {
let max_cells = 1_usize << (2 * level); // 4^level
n.population as f32 / max_cells as f32
}
None => 0.0,
}
}
}
/// 指定した座標とレベルでのNodeを取得
/// (x, y)を含むlevelサイズのノードを返す
fn get_node_at(&self, x: isize, y: isize, level: usize) -> Option<Rc<Node>> {
let top = self.top.as_ref()?;
if level > top.level {
return None;
}
// ツリーの範囲をチェック
// 中心が(0,0)なので、範囲は -2^(L-1) から 2^(L-1) - 1
let half_size = 1_isize << (top.level - 1);
if x < -half_size || x >= half_size || y < -half_size || y >= half_size {
return None;
}
// セルサイズ（1方向）
let cell_size = 1_isize << level;
// (x, y)が属するセルの左下座標
let cell_x = (x.div_euclid(cell_size)) * cell_size;
let cell_y = (y.div_euclid(cell_size)) * cell_size;
Self::get_node_at_recursive(Rc::clone(top), 0, 0, top.level, cell_x, cell_y, level)
}
#[allow(clippy::too_many_arguments)]
fn get_node_at_recursive(
node: Rc<Node>,
node_center_x: isize,
node_center_y: isize,
node_level: usize,
target_x: isize,
target_y: isize,
target_level: usize,
) -> Option<Rc<Node>> {
if node_level == target_level {
return Some(node);
}
// 子ノードの中心へのオフセット
let child_offset = 1_isize << (node_level - 2).max(0);
let children = node.get_children();
// ターゲットがどの象限にあるか判定（node_center_x/yが境界）
let (child, child_x, child_y) = if target_y >= node_center_y {
if target_x >= node_center_x {
// ne: 中心は (node_center + offset, node_center + offset)
(
Rc::clone(&children.ne),
node_center_x + child_offset,
node_center_y + child_offset,
)
} else {
// nw: 中心は (node_center - offset, node_center + offset)
(
Rc::clone(&children.nw),
node_center_x - child_offset,
node_center_y + child_offset,
)
}
} else if target_x >= node_center_x {
// se: 中心は (node_center + offset, node_center - offset)
(
Rc::clone(&children.se),
node_center_x + child_offset,
node_center_y - child_offset,
)
} else {
// sw: 中心は (node_center - offset, node_center - offset)
(
Rc::clone(&children.sw),
node_center_x - child_offset,
node_center_y - child_offset,
)
};
Self::get_node_at_recursive(
child,
child_x,
child_y,
node_level - 1,
target_x,
target_y,
target_level,
)
}
}
impl Node {
fn get_children(&self) -> &Children {
self.children.as_ref().unwrap()
}
fn get_grand_children(&self) -> GrandChildren {
let err1 = "unable to unwrap child (and expecting grand-children)";
let err2 = "unable to unwrap grand-children";
GrandChildren {
nwnw: Rc::clone(
&self
.children
.as_ref()
.expect(err1)
.nw
.children
.as_ref()
.expect(err2)
.nw,
),
nwne: Rc::clone(
&self
.children
.as_ref()
.expect(err1)
.nw
.children
.as_ref()
.expect(err2)
.ne,
),
nwsw: Rc::clone(
&self
.children
.as_ref()
.expect(err1)
.nw
.children
.as_ref()
.expect(err2)
.sw,
),
nwse: Rc::clone(
&self
.children
.as_ref()
.expect(err1)
.nw
.children
.as_ref()
.expect(err2)
.se,
),
nenw: Rc::clone(
&self
.children
.as_ref()
.expect(err1)
.ne
.children
.as_ref()
.expect(err2)
.nw,
),
nene: Rc::clone(
&self
.children
.as_ref()
.expect(err1)
.ne
.children
.as_ref()
.expect(err2)
.ne,
),
nesw: Rc::clone(
&self
.children
.as_ref()
.expect(err1)
.ne
.children
.as_ref()
.expect(err2)
.sw,
),
nese: Rc::clone(
&self
.children
.as_ref()
.expect(err1)
.ne
.children
.as_ref()
.expect(err2)
.se,
),
swnw: Rc::clone(
&self
.children
.as_ref()
.expect(err1)
.sw
.children
.as_ref()
.expect(err2)
.nw,
),
swne: Rc::clone(
&self
.children
.as_ref()
.expect(err1)
.sw
.children
.as_ref()
.expect(err2)
.ne,
),
swsw: Rc::clone(
&self
.children
.as_ref()
.expect(err1)
.sw
.children
.as_ref()
.expect(err2)
.sw,
),
swse: Rc::clone(
&self
.children
.as_ref()
.expect(err1)
.sw
.children
.as_ref()
.expect(err2)
.se,
),
senw: Rc::clone(
&self
.children
.as_ref()
.expect(err1)
.se
.children
.as_ref()
.expect(err2)
.nw,
),
sene: Rc::clone(
&self
.children
.as_ref()
.expect(err1)
.se
.children
.as_ref()
.expect(err2)
.ne,
),
sesw: Rc::clone(
&self
.children
.as_ref()
.expect(err1)
.se
.children
.as_ref()
.expect(err2)
.sw,
),
sese: Rc::clone(
&self
.children
.as_ref()
.expect(err1)
.se
.children
.as_ref()
.expect(err2)
.se,
),
}
}
fn get_grand_automata(&self) -> GrandAutomata {
if self.level != 2 {
panic!("node must be at level 2 to get automatas");
}
let grand_children = self.get_grand_children();
GrandAutomata {
nwnw: grand_children.nwnw.as_automata(),
nwne: grand_children.nwne.as_automata(),
nwsw: grand_children.nwsw.as_automata(),
nwse: grand_children.nwse.as_automata(),
nenw: grand_children.nenw.as_automata(),
nene: grand_children.nene.as_automata(),
nesw: grand_children.nesw.as_automata(),
nese: grand_children.nese.as_automata(),
swnw: grand_children.swnw.as_automata(),
swne: grand_children.swne.as_automata(),
swsw: grand_children.swsw.as_automata(),
swse: grand_children.swse.as_automata(),
senw: grand_children.senw.as_automata(),
sene: grand_children.sene.as_automata(),
sesw: grand_children.sesw.as_automata(),
sese: grand_children.sese.as_automata(),
}
}
fn as_automata(&self) -> Automata {
Automata::from(self.population)
}
}
impl Hash for Node {
fn hash<H: Hasher>(&self, state: &mut H) {
self.hash.hash(state);
}
}
impl Children {
fn from(nw: &Rc<Node>, ne: &Rc<Node>, sw: &Rc<Node>, se: &Rc<Node>) -> Self {
Self {
nw: Rc::clone(nw),
ne: Rc::clone(ne),
sw: Rc::clone(sw),
se: Rc::clone(se),
}
}
}
fn calculate_hash(children: &Children) -> u64 {
let mut state = DefaultHasher::new();
children.hash(&mut state);
state.finish()
}
#[allow(clippy::too_many_arguments)]
fn simb3s23(
center: Automata,
nw: Automata,
n: Automata,
ne: Automata,
w: Automata,
e: Automata,
sw: Automata,
s: Automata,
se: Automata,
) -> Automata {
let neighbors = (nw as u8)
+ (n as u8)
+ (ne as u8)
+ (w as u8)
+ (e as u8)
+ (sw as u8)
+ (s as u8)
+ (se as u8);
match (center, neighbors) {
(Automata::Alive, 2) | (Automata::Alive, 3) => Automata::Alive,
(Automata::Dead, 3) => Automata::Alive,
_ => Automata::Dead,
}
}