blob: bca0d4726cd80f68bb32f3c818338cce91e238f7 [file] [log] [blame]
 // 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) } } }) } } }