package/gonum
1
0
Fork 0
mirror of https://github.com/gonum/gonum.git synced 2025-10-31 18:42:45 +08:00
Code Issues Actions Packages Projects Releases Wiki Activity
Files
085e2bc63c50abe5c561181678700829f358b200
gonum/spatial/barneshut/barneshut2_test.go

547 lines
15 KiB
Go
Raw Blame History

// Copyright ©2019 The Gonum Authors. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.
package barneshut
import (
"fmt"
"math"
"reflect"
"testing"
"golang.org/x/exp/rand"
"gonum.org/v1/gonum/floats/scalar"
"gonum.org/v1/gonum/spatial/r2"
)
type particle2 struct {
x, y, m float64
name string
}
func (p particle2) Coord2() r2.Vec { return r2.Vec{X: p.x, Y: p.y} }
func (p particle2) Mass() float64 { return p.m }
var planeTests = []struct {
name string
particles []particle2
want *Plane
}{
{
name: "nil",
particles: nil,
want: &Plane{},
},
{
name: "empty",
particles: []particle2{},
want: &Plane{Particles: []Particle2{}},
},
{
name: "one",
particles: []particle2{{m: 1}}, // Must have a mass to avoid vacuum decay.
want: &Plane{
root: tile{
particle: particle2{x: 0, y: 0, m: 1},
bounds: r2.Box{Min: r2.Vec{X: 0, Y: 0}, Max: r2.Vec{X: 0, Y: 0}},
center: r2.Vec{X: 0, Y: 0},
mass: 1,
},
Particles: []Particle2{
particle2{m: 1},
},
},
},
{
name: "3 corners",
particles: []particle2{
{x: 1, y: 1, m: 1},
{x: -1, y: 1, m: 1},
{x: -1, y: -1, m: 1},
},
want: &Plane{
root: tile{
bounds: r2.Box{Min: r2.Vec{X: -1, Y: -1}, Max: r2.Vec{X: 1, Y: 1}},
nodes: [4]*tile{
se: {
particle: particle2{x: 1, y: 1, m: 1},
bounds: r2.Box{Min: r2.Vec{X: 0, Y: 0}, Max: r2.Vec{X: 1, Y: 1}},
center: r2.Vec{X: 1, Y: 1}, mass: 1,
},
sw: {
particle: particle2{x: -1, y: 1, m: 1},
bounds: r2.Box{Min: r2.Vec{X: -1, Y: 0}, Max: r2.Vec{X: 0, Y: 1}},
center: r2.Vec{X: -1, Y: 1}, mass: 1,
},
nw: {
particle: particle2{x: -1, y: -1, m: 1},
bounds: r2.Box{Min: r2.Vec{X: -1, Y: -1}, Max: r2.Vec{X: 0, Y: 0}},
center: r2.Vec{X: -1, Y: -1}, mass: 1,
},
},
center: r2.Vec{X: -0.3333333333333333, Y: 0.3333333333333333},
mass: 3,
},
Particles: []Particle2{
particle2{x: 1, y: 1, m: 1},
particle2{x: -1, y: 1, m: 1},
particle2{x: -1, y: -1, m: 1},
},
},
},
{
name: "4 corners",
particles: []particle2{
{x: 1, y: 1, m: 1},
{x: -1, y: 1, m: 1},
{x: 1, y: -1, m: 1},
{x: -1, y: -1, m: 1},
},
want: &Plane{
root: tile{
bounds: r2.Box{Min: r2.Vec{X: -1, Y: -1}, Max: r2.Vec{X: 1, Y: 1}},
nodes: [4]*tile{
{
particle: particle2{x: 1, y: -1, m: 1},
bounds: r2.Box{Min: r2.Vec{X: 0, Y: -1}, Max: r2.Vec{X: 1, Y: 0}},
center: r2.Vec{X: 1, Y: -1},
mass: 1,
},
{
particle: particle2{x: 1, y: 1, m: 1},
bounds: r2.Box{Min: r2.Vec{X: 0, Y: 0}, Max: r2.Vec{X: 1, Y: 1}},
center: r2.Vec{X: 1, Y: 1},
mass: 1,
},
{
particle: particle2{x: -1, y: 1, m: 1},
bounds: r2.Box{Min: r2.Vec{X: -1, Y: 0}, Max: r2.Vec{X: 0, Y: 1}},
center: r2.Vec{X: -1, Y: 1},
mass: 1,
},
{
particle: particle2{x: -1, y: -1, m: 1},
bounds: r2.Box{Min: r2.Vec{X: -1, Y: -1}, Max: r2.Vec{X: 0, Y: 0}},
center: r2.Vec{X: -1, Y: -1},
mass: 1,
},
},
center: r2.Vec{X: 0, Y: 0},
mass: 4,
},
Particles: []Particle2{
particle2{x: 1, y: 1, m: 1},
particle2{x: -1, y: 1, m: 1},
particle2{x: 1, y: -1, m: 1},
particle2{x: -1, y: -1, m: 1},
},
},
},
{
name: "5 corners",
particles: []particle2{
{x: 1, y: 1, m: 1},
{x: -1, y: 1, m: 1},
{x: 1, y: -1, m: 1},
{x: -1, y: -1, m: 1},
{x: -1.1, y: -1, m: 1},
},
want: &Plane{
root: tile{
bounds: r2.Box{Min: r2.Vec{X: -1.1, Y: -1}, Max: r2.Vec{X: 1, Y: 1}},
nodes: [4]*tile{
{
particle: particle2{x: 1, y: -1, m: 1},
bounds: r2.Box{Min: r2.Vec{X: -0.050000000000000044, Y: -1}, Max: r2.Vec{X: 1, Y: 0}},
center: r2.Vec{X: 1, Y: -1},
mass: 1,
},
{
particle: particle2{x: 1, y: 1, m: 1},
bounds: r2.Box{Min: r2.Vec{X: -0.050000000000000044, Y: 0}, Max: r2.Vec{X: 1, Y: 1}},
center: r2.Vec{X: 1, Y: 1},
mass: 1,
},
{
particle: particle2{x: -1, y: 1, m: 1},
bounds: r2.Box{Min: r2.Vec{X: -1.1, Y: 0}, Max: r2.Vec{X: -0.050000000000000044, Y: 1}},
center: r2.Vec{X: -1, Y: 1},
mass: 1,
},
{
bounds: r2.Box{Min: r2.Vec{X: -1.1, Y: -1}, Max: r2.Vec{X: -0.050000000000000044, Y: 0}},
nodes: [4]*tile{
nw: {
bounds: r2.Box{Min: r2.Vec{X: -1.1, Y: -1}, Max: r2.Vec{X: -0.5750000000000001, Y: -0.5}},
nodes: [4]*tile{
nw: {
bounds: r2.Box{Min: r2.Vec{X: -1.1, Y: -1}, Max: r2.Vec{X: -0.8375000000000001, Y: -0.75}},
nodes: [4]*tile{
nw: {
bounds: r2.Box{Min: r2.Vec{X: -1.1, Y: -1}, Max: r2.Vec{X: -0.9687500000000001, Y: -0.875}},
nodes: [4]*tile{
ne: {
particle: particle2{x: -1, y: -1, m: 1},
bounds: r2.Box{Min: r2.Vec{X: -1.034375, Y: -1}, Max: r2.Vec{X: -0.9687500000000001, Y: -0.9375}},
center: r2.Vec{X: -1, Y: -1},
mass: 1,
},
nw: {
particle: particle2{x: -1.1, y: -1, m: 1},
bounds: r2.Box{Min: r2.Vec{X: -1.1, Y: -1}, Max: r2.Vec{X: -1.034375, Y: -0.9375}},
center: r2.Vec{X: -1.1, Y: -1},
mass: 1,
},
},
center: r2.Vec{X: -1.05, Y: -1},
mass: 2,
},
},
center: r2.Vec{X: -1.05, Y: -1},
mass: 2,
},
},
center: r2.Vec{X: -1.05, Y: -1},
mass: 2,
},
},
center: r2.Vec{X: -1.05, Y: -1},
mass: 2,
},
},
center: r2.Vec{X: -0.22000000000000003, Y: -0.2},
mass: 5,
},
Particles: []Particle2{
particle2{x: 1, y: 1, m: 1},
particle2{x: -1, y: 1, m: 1},
particle2{x: 1, y: -1, m: 1},
particle2{x: -1, y: -1, m: 1},
particle2{x: -1.1, y: -1, m: 1},
},
},
},
{
// Note that the code here subdivides the space differently to
// how it is split in the example, since Plane makes a minimum
// bounding box based on the data, while the example does not.
name: "http://arborjs.org/docs/barnes-hut example",
particles: []particle2{
{x: 64.5, y: 81.5, m: 1, name: "A"},
{x: 242, y: 34, m: 1, name: "B"},
{x: 199, y: 69, m: 1, name: "C"},
{x: 285, y: 106.5, m: 1, name: "D"},
{x: 170, y: 194.5, m: 1, name: "E"},
{x: 42.5, y: 334.5, m: 1, name: "F"},
{x: 147, y: 309, m: 1, name: "G"},
{x: 236.5, y: 324, m: 1, name: "H"},
},
want: &Plane{
root: tile{
bounds: r2.Box{Min: r2.Vec{X: 42.5, Y: 34}, Max: r2.Vec{X: 285, Y: 334.5}},
nodes: [4]*tile{
{
bounds: r2.Box{Min: r2.Vec{X: 163.75, Y: 34}, Max: r2.Vec{X: 285, Y: 184.25}},
nodes: [4]*tile{
ne: {
bounds: r2.Box{Min: r2.Vec{X: 224.375, Y: 34}, Max: r2.Vec{X: 285, Y: 109.125}},
nodes: [4]*tile{
se: {
particle: particle2{x: 285, y: 106.5, m: 1, name: "D"},
bounds: r2.Box{Min: r2.Vec{X: 254.6875, Y: 71.5625}, Max: r2.Vec{X: 285, Y: 109.125}},
center: r2.Vec{X: 285, Y: 106.5},
mass: 1,
},
nw: {
particle: particle2{x: 242, y: 34, m: 1, name: "B"},
bounds: r2.Box{Min: r2.Vec{X: 224.375, Y: 34}, Max: r2.Vec{X: 254.6875, Y: 71.5625}},
center: r2.Vec{X: 242, Y: 34},
mass: 1,
},
},
center: r2.Vec{X: 263.5, Y: 70.25},
mass: 2,
},
nw: {
particle: particle2{x: 199, y: 69, m: 1, name: "C"},
bounds: r2.Box{Min: r2.Vec{X: 163.75, Y: 34}, Max: r2.Vec{X: 224.375, Y: 109.125}},
center: r2.Vec{X: 199, Y: 69},
mass: 1,
},
},
center: r2.Vec{X: 242, Y: 69.83333333333333},
mass: 3,
},
{
bounds: r2.Box{Min: r2.Vec{X: 163.75, Y: 184.25}, Max: r2.Vec{X: 285, Y: 334.5}},
nodes: [4]*tile{
se: {
particle: particle2{x: 236.5, y: 324, m: 1, name: "H"},
bounds: r2.Box{Min: r2.Vec{X: 224.375, Y: 259.375}, Max: r2.Vec{X: 285, Y: 334.5}},
center: r2.Vec{X: 236.5, Y: 324},
mass: 1,
},
nw: {
particle: particle2{x: 170, y: 194.5, m: 1, name: "E"},
bounds: r2.Box{Min: r2.Vec{X: 163.75, Y: 184.25}, Max: r2.Vec{X: 224.375, Y: 259.375}},
center: r2.Vec{X: 170, Y: 194.5},
mass: 1,
},
},
center: r2.Vec{X: 203.25, Y: 259.25},
mass: 2,
},
{
bounds: r2.Box{Min: r2.Vec{X: 42.5, Y: 184.25}, Max: r2.Vec{X: 163.75, Y: 334.5}},
nodes: [4]*tile{
se: {
particle: particle2{x: 147, y: 309, m: 1, name: "G"},
bounds: r2.Box{Min: r2.Vec{X: 103.125, Y: 259.375}, Max: r2.Vec{X: 163.75, Y: 334.5}},
center: r2.Vec{X: 147, Y: 309},
mass: 1,
},
sw: {
particle: particle2{x: 42.5, y: 334.5, m: 1, name: "F"},
bounds: r2.Box{Min: r2.Vec{X: 42.5, Y: 259.375}, Max: r2.Vec{X: 103.125, Y: 334.5}},
center: r2.Vec{X: 42.5, Y: 334.5},
mass: 1,
},
},
center: r2.Vec{X: 94.75, Y: 321.75},
mass: 2,
},
{
particle: particle2{x: 64.5, y: 81.5, m: 1, name: "A"},
bounds: r2.Box{Min: r2.Vec{X: 42.5, Y: 34}, Max: r2.Vec{X: 163.75, Y: 184.25}},
center: r2.Vec{X: 64.5, Y: 81.5},
mass: 1,
},
},
center: r2.Vec{X: 173.3125, Y: 181.625},
mass: 8,
},
Particles: []Particle2{
particle2{x: 64.5, y: 81.5, m: 1, name: "A"},
particle2{x: 242, y: 34, m: 1, name: "B"},
particle2{x: 199, y: 69, m: 1, name: "C"},
particle2{x: 285, y: 106.5, m: 1, name: "D"},
particle2{x: 170, y: 194.5, m: 1, name: "E"},
particle2{x: 42.5, y: 334.5, m: 1, name: "F"},
particle2{x: 147, y: 309, m: 1, name: "G"},
particle2{x: 236.5, y: 324, m: 1, name: "H"},
},
},
},
}
func TestPlane(t *testing.T) {
t.Parallel()
const tol = 1e-15
for _, test := range planeTests {
var particles []Particle2
if test.particles != nil {
particles = make([]Particle2, len(test.particles))
}
for i, p := range test.particles {
particles[i] = p
}
got, err := NewPlane(particles)
if err != nil {
t.Errorf("unexpected error: %v", err)
continue
}
if test.want != nil && !reflect.DeepEqual(got, test.want) {
t.Errorf("unexpected result for %q: got:%v want:%v", test.name, got, test.want)
}
// Recursively check all internal centers of mass.
walkPlane(&got.root, func(tl *tile) {
var sub []Particle2
walkPlane(tl, func(tl *tile) {
if tl.particle != nil {
sub = append(sub, tl.particle)
}
})
center, mass := centerOfMass2(sub)
if !scalar.EqualWithinAbsOrRel(center.X, tl.center.X, tol, tol) || !scalar.EqualWithinAbsOrRel(center.Y, tl.center.Y, tol, tol) {
t.Errorf("unexpected result for %q for center of mass: got:%f want:%f", test.name, tl.center, center)
}
if !scalar.EqualWithinAbsOrRel(mass, tl.mass, tol, tol) {
t.Errorf("unexpected result for %q for total mass: got:%f want:%f", test.name, tl.mass, mass)
}
})
}
}
func centerOfMass2(particles []Particle2) (center r2.Vec, mass float64) {
for _, p := range particles {
m := p.Mass()
mass += m
c := p.Coord2()
center.X += c.X * m
center.Y += c.Y * m
}
if mass != 0 {
center.X /= mass
center.Y /= mass
}
return center, mass
}
func walkPlane(t *tile, fn func(*tile)) {
if t == nil {
return
}
fn(t)
for _, q := range t.nodes {
walkPlane(q, fn)
}
}
func TestPlaneForceOn(t *testing.T) {
t.Parallel()
const (
size = 1000
tol = 0.07
)
for _, n := range []int{3e3, 1e4, 3e4} {
rnd := rand.New(rand.NewSource(1))
particles := make([]Particle2, n)
for i := range particles {
particles[i] = particle2{x: size * rnd.Float64(), y: size * rnd.Float64(), m: 1}
}
moved := make([]r2.Vec, n)
for i, p := range particles {
var v r2.Vec
m := p.Mass()
pv := p.Coord2()
for _, e := range particles {
v = v.Add(Gravity2(p, e, m, e.Mass(), e.Coord2().Sub(pv)))
}
moved[i] = p.Coord2().Add(v)
}
plane, err := NewPlane(particles)
if err != nil {
t.Errorf("unexpected error: %v", err)
continue
}
for _, theta := range []float64{0, 0.3, 0.6, 0.9} {
theta := theta
t.Run(fmt.Sprintf("%d-body/theta=%v", len(particles), theta), func(t *testing.T) {
t.Parallel()
var ssd, sd float64
var calls int
for i, p := range particles {
v := plane.ForceOn(p, theta, func(p1, p2 Particle2, m1, m2 float64, v r2.Vec) r2.Vec {
calls++
return Gravity2(p1, p2, m1, m2, v)
})
pos := p.Coord2().Add(v)
d := moved[i].Sub(pos)
ssd += d.X*d.X + d.Y*d.Y
sd += math.Hypot(d.X, d.Y)
}
rmsd := math.Sqrt(ssd / float64(len(particles)))
if rmsd > tol {
t.Error("RMSD for approximation too high")
}
t.Logf("rmsd=%.4v md=%.4v calls/particle=%.5v",
rmsd, sd/float64(len(particles)), float64(calls)/float64(len(particles)))
})
}
}
}
var (
fv2sink r2.Vec
planeSink *Plane
)
func BenchmarkNewPlane(b *testing.B) {
for _, n := range []int{1e3, 1e4, 1e5, 1e6} {
rnd := rand.New(rand.NewSource(1))
particles := make([]Particle2, n)
for i := range particles {
particles[i] = particle2{x: rnd.Float64(), y: rnd.Float64(), m: 1}
}
b.ResetTimer()
var err error
b.Run(fmt.Sprintf("%d-body", len(particles)), func(b *testing.B) {
for i := 0; i < b.N; i++ {
planeSink, err = NewPlane(particles)
if err != nil {
b.Fatalf("unexpected error: %v", err)
}
}
})
}
}
func BenchmarkPlaneForceOn(b *testing.B) {
for _, n := range []int{1e3, 1e4, 1e5} {
for _, theta := range []float64{0, 0.1, 0.5, 1, 1.5} {
if n > 1e4 && theta < 0.5 {
// Don't run unreasonably long benchmarks.
continue
}
rnd := rand.New(rand.NewSource(1))
particles := make([]Particle2, n)
for i := range particles {
particles[i] = particle2{x: rnd.Float64(), y: rnd.Float64(), m: 1}
}
plane, err := NewPlane(particles)
if err != nil {
b.Fatalf("unexpected error: %v", err)
}
b.ResetTimer()
b.Run(fmt.Sprintf("%d-body/theta=%v", len(particles), theta), func(b *testing.B) {
for i := 0; i < b.N; i++ {
for _, p := range particles {
fv2sink = plane.ForceOn(p, theta, Gravity2)
}
}
})
}
}
}
func BenchmarkPlaneFull(b *testing.B) {
for _, n := range []int{1e3, 1e4, 1e5} {
for _, theta := range []float64{0, 0.1, 0.5, 1, 1.5} {
if n > 1e4 && theta < 0.5 {
// Don't run unreasonably long benchmarks.
continue
}
rnd := rand.New(rand.NewSource(1))
particles := make([]Particle2, n)
for i := range particles {
particles[i] = particle2{x: rnd.Float64(), y: rnd.Float64(), m: 1}
}
b.ResetTimer()
b.Run(fmt.Sprintf("%d-body/theta=%v", len(particles), theta), func(b *testing.B) {
for i := 0; i < b.N; i++ {
plane, err := NewPlane(particles)
if err != nil {
b.Fatalf("unexpected error: %v", err)
}
for _, p := range particles {
fv2sink = plane.ForceOn(p, theta, Gravity2)
}
}
})
}
}
}
Reference in New Issue View Git Blame Copy Permalink