spatial/barneshut: catch finitely too far points

BarrensZeppelin · Jun 28, 2019 · 6006dfb · 6006dfb
1 parent e34e6b9
commit 6006dfb
Show file tree

Hide file tree

Showing 4 changed files with 110 additions and 46 deletions.
diff --git a/spatial/barneshut/barneshut2.go b/spatial/barneshut/barneshut2.go
@@ -5,6 +5,7 @@
 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
-}
-
-func isInfR2(v r2.Vec) bool {
-	return math.IsInf(v.X, 0) || math.IsInf(v.Y, 0)
+	err := q.Reset()
+	if err != nil {
+		return nil, err
+	}
+	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,15 +94,28 @@ 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.
 	for _, e := range q.Particles[1:] {
 		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
 }
 

diff --git a/spatial/barneshut/barneshut2_test.go b/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)
 					}

diff --git a/spatial/barneshut/barneshut3.go b/spatial/barneshut/barneshut3.go
@@ -5,6 +5,7 @@
 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
-}
-
-func isInfR3(v r3.Vec) bool {
-	return math.IsInf(v.X, 0) || math.IsInf(v.Y, 0) || math.IsInf(v.Z, 0)
+	err := q.Reset()
+	if err != nil {
+		return nil, err
+	}
+	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,15 +100,28 @@ 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.
 	for _, e := range q.Particles[1:] {
 		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
 }
 

diff --git a/spatial/barneshut/barneshut3_test.go b/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)
 					}