spatial/barneshut: catch finitely too far points
MODIFY
spatial/barneshut/barneshut2.go
spatial/barneshut/barneshut2.go
package barneshut
import (
+ "errors"
"fmt"
"math"
@@ -46,32 +47,28 @@ type Plane struct {
Particles []Particle2
}
-// NewPlane returns a new Plane.
-//
-// Points in p must not be infinitely distant, otherwise NewPlane will panic.
-func NewPlane(p []Particle2) *Plane {
- for _, l := range p {
- if isInfR2(l.Coord2()) {
- panic("barneshut: point at infinity")
- }
- }
+// NewPlane returns a new Plane. If the plane is too large to allow
+// particle coordinates to be distinguished due to floating point
+// precision limits, NewPlane will return a non-nil error.
+func NewPlane(p []Particle2) (*Plane, error) {
q := Plane{Particles: p}
- q.Reset()
- return &q
+ err := q.Reset()
+ if err != nil {
+ return nil, err
}
-
-func isInfR2(v r2.Vec) bool {
- return math.IsInf(v.X, 0) || math.IsInf(v.Y, 0)
+ return &q, nil
}
// Reset reconstructs the Barnes-Hut tree. Reset must be called if the
// Particles field or elements of Particles have been altered, unless
// ForceOn is called with theta=0 or no data structures have been
-// previously built.
-func (q *Plane) Reset() {
+// previously built. If the plane is too large to allow particle
+// coordinates to be distinguished due to floating point precision
+// limits, Reset will return a non-nil error.
+func (q *Plane) Reset() (err error) {
if len(q.Particles) == 0 {
q.root = tile{}
- return
+ return nil
}
q.root = tile{
@@ -97,6 +94,16 @@ func (q *Plane) Reset() {
}
}
+ defer func() {
+ switch r := recover(); r {
+ case nil:
+ case volumeTooBig:
+ err = volumeTooBig
+ default:
+ panic(r)
+ }
+ }()
+
// TODO(kortschak): Partially parallelise this by
// choosing the direction and using one of four
// goroutines to work on each root quadrant.
@@ -104,8 +111,11 @@ func (q *Plane) Reset() {
q.root.insert(e)
}
q.root.summarize()
+ return nil
}
+var planeTooBig = errors.New("barneshut: plane too big")
+
// ForceOn returns a force vector on p given p's mass and the force function, f,
// using the Barnes-Hut theta approximation parameter.
//
@@ -205,6 +215,7 @@ func quadrantOf(b r2.Box, p Particle2) int {
// splitPlane returns a quadrant subdivision of b in the given direction.
func splitPlane(b r2.Box, dir int) r2.Box {
+ old := b
halfX := (b.Max.X - b.Min.X) / 2
halfY := (b.Max.Y - b.Min.Y) / 2
switch dir {
@@ -221,6 +232,9 @@ func splitPlane(b r2.Box, dir int) r2.Box {
b.Max.X -= halfX
b.Max.Y -= halfY
}
+ if b == old {
+ panic(planeTooBig)
+ }
return b
}

MODIFY
spatial/barneshut/barneshut2_test.go
spatial/barneshut/barneshut2_test.go
@@ -353,7 +353,11 @@ func TestPlane(t *testing.T) {
particles[i] = p
}
- got := NewPlane(particles)
+ 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)
@@ -426,7 +430,11 @@ func TestPlaneForceOn(t *testing.T) {
moved[i] = p.Coord2().Add(v)
}
- plane := NewPlane(particles)
+ 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} {
t.Run(fmt.Sprintf("%d-body/theta=%v", len(particles), theta), func(t *testing.T) {
var ssd, sd float64
@@ -465,9 +473,13 @@ func BenchmarkNewPlane(b *testing.B) {
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 = NewPlane(particles)
+ planeSink, err = NewPlane(particles)
+ if err != nil {
+ b.Fatalf("unexpected error: %v", err)
+ }
}
})
}
@@ -485,7 +497,10 @@ func BenchmarkPlaneForceOn(b *testing.B) {
for i := range particles {
particles[i] = particle2{x: rnd.Float64(), y: rnd.Float64(), m: 1}
}
- plane := NewPlane(particles)
+ 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++ {
@@ -513,7 +528,10 @@ func BenchmarkPlaneFull(b *testing.B) {
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 := NewPlane(particles)
+ plane, err := NewPlane(particles)
+ if err != nil {
+ b.Fatalf("unexpected error: %v", err)
+ }
for _, p := range particles {
fv2sink = plane.ForceOn(p, theta, Gravity2)
}

MODIFY
spatial/barneshut/barneshut3.go
spatial/barneshut/barneshut3.go
package barneshut
import (
+ "errors"
"fmt"
"math"
@@ -46,32 +47,28 @@ type Volume struct {
Particles []Particle3
}
-// NewVolume returns a new Volume.
-//
-// Points in p must not be infinitely distant, otherwise NewVolume will panic.
-func NewVolume(p []Particle3) *Volume {
- for _, l := range p {
- if isInfR3(l.Coord3()) {
- panic("barneshut: point at infinity")
- }
- }
+// NewVolume returns a new Volume. If the volume is too large to allow
+// particle coordinates to be distinguished due to floating point
+// precision limits, NewVolume will return a non-nil error.
+func NewVolume(p []Particle3) (*Volume, error) {
q := Volume{Particles: p}
- q.Reset()
- return &q
+ err := q.Reset()
+ if err != nil {
+ return nil, err
}
-
-func isInfR3(v r3.Vec) bool {
- return math.IsInf(v.X, 0) || math.IsInf(v.Y, 0) || math.IsInf(v.Z, 0)
+ return &q, nil
}
// Reset reconstructs the Barnes-Hut tree. Reset must be called if the
// Particles field or elements of Particles have been altered, unless
// ForceOn is called with theta=0 or no data structures have been
-// previously built.
-func (q *Volume) Reset() {
+// previously built. If the volume is too large to allow particle
+// coordinates to be distinguished due to floating point precision
+// limits, Reset will return a non-nil error.
+func (q *Volume) Reset() (err error) {
if len(q.Particles) == 0 {
q.root = bucket{}
- return
+ return nil
}
q.root = bucket{
@@ -103,6 +100,16 @@ func (q *Volume) Reset() {
}
}
+ defer func() {
+ switch r := recover(); r {
+ case nil:
+ case volumeTooBig:
+ err = volumeTooBig
+ default:
+ panic(r)
+ }
+ }()
+
// TODO(kortschak): Partially parallelise this by
// choosing the direction and using one of eight
// goroutines to work on each root octant.
@@ -110,8 +117,11 @@ func (q *Volume) Reset() {
q.root.insert(e)
}
q.root.summarize()
+ return nil
}
+var volumeTooBig = errors.New("barneshut: volume too big")
+
// ForceOn returns a force vector on p given p's mass and the force function, f,
// using the Barnes-Hut theta approximation parameter.
//
@@ -233,6 +243,7 @@ func octantOf(b r3.Box, p Particle3) int {
// splitVolume returns an octant subdivision of b in the given direction.
func splitVolume(b r3.Box, dir int) r3.Box {
+ old := b
halfX := (b.Max.X - b.Min.X) / 2
halfY := (b.Max.Y - b.Min.Y) / 2
halfZ := (b.Max.Z - b.Min.Z) / 2
@@ -270,6 +281,9 @@ func splitVolume(b r3.Box, dir int) r3.Box {
b.Max.Y -= halfY
b.Min.Z += halfZ
}
+ if b == old {
+ panic(volumeTooBig)
+ }
return b
}

MODIFY
spatial/barneshut/barneshut3_test.go
spatial/barneshut/barneshut3_test.go
@@ -348,7 +348,11 @@ func TestVolume(t *testing.T) {
particles[i] = p
}
- got := NewVolume(particles)
+ got, err := NewVolume(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)
@@ -423,7 +427,11 @@ func TestVolumeForceOn(t *testing.T) {
moved[i] = p.Coord3().Add(v)
}
- volume := NewVolume(particles)
+ volume, err := NewVolume(particles)
+ if err != nil {
+ t.Errorf("unexpected error: %v", err)
+ continue
+ }
for _, theta := range []float64{0, 0.3, 0.6, 0.9} {
t.Run(fmt.Sprintf("%d-body/theta=%v", len(particles), theta), func(t *testing.T) {
var ssd, sd float64
@@ -462,9 +470,13 @@ func BenchmarkNewVolume(b *testing.B) {
particles[i] = particle3{x: rnd.Float64(), y: rnd.Float64(), z: 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++ {
- volumeSink = NewVolume(particles)
+ volumeSink, err = NewVolume(particles)
+ if err != nil {
+ b.Fatalf("unexpected error: %v", err)
+ }
}
})
}
@@ -482,7 +494,10 @@ func BenchmarkVolumeForceOn(b *testing.B) {
for i := range particles {
particles[i] = particle3{x: rnd.Float64(), y: rnd.Float64(), z: rnd.Float64(), m: 1}
}
- volume := NewVolume(particles)
+ volume, err := NewVolume(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++ {
@@ -510,7 +525,10 @@ func BenchmarkVolumeFull(b *testing.B) {
b.ResetTimer()
b.Run(fmt.Sprintf("%d-body/theta=%v", len(particles), theta), func(b *testing.B) {
for i := 0; i < b.N; i++ {
- volume := NewVolume(particles)
+ volume, err := NewVolume(particles)
+ if err != nil {
+ b.Fatalf("unexpected error: %v", err)
+ }
for _, p := range particles {
fv3sink = volume.ForceOn(p, theta, Gravity3)
}

spatial/barneshut: catch finitely too far points