Barnes-Hut树用于快速N体力计算

感谢您查看此包。这是Zhifeng实现加速N体力计算的Barnes-Hut树的编码实践。

该包由两部分组成：[BarnesHutTree]结构体和包含一些辅助函数的[utils]模块。[BarnesHutTree]是一种数据结构，通过近似处理相对较远的“物体”组为超级节点，无需遍历每个单个物体来加速N体力计算。

为了使力计算更加可定制，[BarnesHutTree]接受一个闭包（匿名函数）来决定当前节点是否可以被视为“足够远”，以及一个闭包来决定如何计算力或其他目标节点与超级节点之间的关系。[utils]模块包含在关于图绘制（论文）中提到的排斥力计算的几个闭包工厂函数。

使用示例

use zhifeng_impl_barnes_hut_tree as zbht;

use zbht::BarnesHutTree as BHTree;

let k = 1.0; // Custom Parameters For Calculation
let c = 0.2;

let mut bht: BHTree<2> = BHTree::new();
bht.push(&[-1.0,1.0]);
bht.push(&[1.0,1.0]);

let mut ans_displacement = [0.0; 2];

let calc_fn = zbht::utils::factory_of_repulsive_displacement_calc_fn::<2>(k, c);

// Since the closure implies every super node is not "far" enough,
// the line calculates the exact displacement acting on value 0.
let is_super_fn = |_: &[f64;2],_:&[f64;2],_:f64| -> bool {false}; // ignoring super nodes
let calc_fn = zbht::utils::factory_of_repulsive_displacement_calc_fn::<2>(k, c);

bht.calc_force_on_value(0, &is_super_fn, &calc_fn, &mut ans_displacement);
assert_eq!(ans_displacement, [(-2.0 * k * k * c) / (2.0 * 2.0), 0.0]);

let theta = 1.2; // Custom Parameter For Super Node Determination

let is_super_fn = zbht::utils::factory_of_is_super_node_fn::<2>(theta);
// Calculating the approximated displacement acting on value 0.
bht.calc_force_on_value(0, &is_super_fn, &calc_fn, &mut ans_displacement);

设计和命名

为了阐明构造函数中的参数，以下是我在实现中使用的关键术语。

• 在“n体”计算的上下文中，“物体”被称为value。（目前，树将每个物体视为相等。）

• Barnes Hut树使用超立方体树快速找到相对接近的节点组，并在计算过程中预先计算一组value的平均位置以加速计算过程。

• 节点包含的value位置的平均值被称为值中心，简称vc。

• 树节点的边界框中心，即超立方体的中心，简称bc。

• 树节点的边界框半径，即超立方体的半宽，简称br。

为了使设计更安全，我在内部使用vec来存储节点：内部节点和叶节点。当一个节点需要被删除时，如果它不是最后一个节点，则vec中的最后一个节点将替换其位置并更新基于索引的虚拟“指针”。

总体思路

首先，[BarnesHutTree] 以默认或指定的边界超立方体开始。对于每个维度，超立方体的中心将维度分成两部分，从而将超空间分成维度的平方。例如，对于二维情况，一个节点可以将所有值分成四个组。[BarnesHutTree] 将尝试递归地划分所有值，直到所有值都分开，最终没有两个值直接属于同一个节点。有了这种结构和关于目标节点与超节点值中心位置之间相对位置的某些标准，我们可以高效地找到超节点并计算力。

但是，存在许多潜在问题：如何处理位于边界上的值、定义维度、处理超出初始边界超立方体的值的添加，以及处理添加两个过于接近或甚至相同的值，这会导致试图将它们分开的无限循环。

处理边界值

超立方的下限是包含的，上限是排除的。例如，对于二维空间，位于超立方体中心的价值属于节点的第一个象限子节点。

处理维度

该树使用模板参数来定义值（物体）位置维度。内部，一个内部节点使用一个 vec 存储有关其子节点的信息。这种方法使得访问子节点变得快速，但所需空间是维数的平方。在使用大量维度时，我们需要小心。

处理超出根边界范围的值

该树将尝试通过创建一个具有更大半径的新内部节点并朝向要包含的值方向（同样适用于上面没有或只有叶节点的最小半径限制）来加倍其宽度。当要包含的值太远时，需要小心，因为这可能导致诸如“无穷大”或“NaN”之类的数值问题。

处理过于接近或相同的值

由于树通常试图将每个值放入分开的叶子中，以计算超节点，如果没有额外的检查和处理，尝试插入两个相同的值将创建一个无限循环，试图将两个值放入不同的叶子节点。

该树目前通过允许我们预先设置一个界限，当超立方体的半径（半宽）小于或等于该界限时，树节点将不会继续划分，来处理两个过于接近的值。在计算“力”时，树将循环遍历同一叶子节点中的每个值，除了目标值本身。

在设置此限制值时，需要格外小心。如果限制过大，单个叶子节点中将会保留过多的值，导致效率降低。如果所有值都在单个叶子节点中，树的行为和效率将与遍历所有节点相同：N体计算的默认实现。目前，即使 f64 在一些划分之后最终达到零，但限制应该大于零且为有限值。目前，默认限制是 1e-8。

特性

序列化

此功能使用 serde 和 serde_json 来帮助序列化树以进行测试、调试和制作可视化。

未检查

此功能使用 get_unchecked、get_unchecked_mut 等，以实现“虚拟”的“内部-vec”节点的快速访问。基于带有 --features unchecked 的 cargo bench 的结果，"unchecked" 功能比默认功能快约 6%。

性能

该包使用 criterion 进行基准测试和 rand 生成随机测试值。为了模拟 [BarnesHutTree] 的常用用例，我使用了一轮遍历所有值，计算它们相应的位移并更新它们的位置作为基准测试标准。

总体恐慌

如果任何 64 位浮点数在构建和值更新过程中变为 "NaN" 或 "Infinity"，则树会崩溃。在将值添加到树之前，我们应该小心处理值的数值。

参考

Hu, Y. (2005). 高效、高质量的力导向图绘制。 Mathematica journal, 10(1), 37-71。

我从这篇论文的第四章“Barnes-Hut 力计算”中获得了 Barnes-Hut 树的灵感。我在我的实现（可能存在错误）中加入了些自己的想法。例如，在论文中，作者没有关注 Barnes-Hut 树的“值移除”部分，因为重建八叉结构既高效又准确，所以我独立地考虑了这个问题。

许可证：AGPL-3.0-only

加速 N 体计算的 Zhifeng 对 Barnes-Hut 树的实现版权所有 (C) 2024 Zhifeng Wang

此程序是免费软件：您可以在自由软件基金会发布的 GNU Affero 通用公共许可证的条款下重新分发和/或修改它，许可证版本 3，或版本 3。

此程序是在希望它有用的前提下分发的，但没有任何保证；甚至没有关于其适销性或适用于特定目的的隐含保证。有关详细信息，请参阅 GNU Affero 通用公共许可证。

您应该已收到随此程序一起提供的 GNU Affero 通用公共许可证副本。如果没有，请参阅 https://www.gnu.org/licenses/。