midpip.go

// A simple experiment to check how fast point in polygon tests with
// a first level spatial index filter are.
//
// This is Free Software covered by the MIT license.
// See https://opensource.org/licenses/MIT for details.
// (c) 2017 by Sascha L. Teichmann (sascha.teichmann@intevation.de)
//
package main

import (
	"flag"
	"fmt"
	"log"
	"math/rand"
	"os"
	"runtime"
	"runtime/pprof"
	"sync"
	"time"

	"github.com/tidwall/rtree"

	shp "github.com/jonas-p/go-shp"

	geo "github.com/kellydunn/golang-geo"
)

type item struct {
	min  []float64
	max  []float64
	poly *geo.Polygon
}

func (it *item) Rect(interface{}) ([]float64, []float64) {
	return it.min, it.max
}

func (it *item) mid() (float64, float64) {
	return 0.5 * (it.min[0] + it.max[0]), 0.5 * (it.min[1] + it.max[1])
}

type point struct{ *geo.Point }

func (p point) Rect(interface{}) ([]float64, []float64) {
	const eps = 0.01
	x, y := p.Lat(), p.Lng()
	return []float64{x - eps, y - eps}, []float64{x + eps, y + eps}
}

type worker struct {
	tree *rtree.RTree
	rnd  *rand.Rand
	hits int
}

func newWorker(tree *rtree.RTree) *worker {
	return &worker{
		tree: tree,
		rnd:  rand.New(rand.NewSource(time.Now().Unix())),
	}
}

func (w *worker) run(n int, items []*item, wg *sync.WaitGroup) {
	defer wg.Done()

	for i := 0; i < n; i++ {
		it := items[w.rnd.Intn(len(items))]
		// The center of the bbox is not always inside the polygon.
		x, y := it.mid()
		gp := geo.NewPoint(x, y)
		w.tree.Search(point{gp}, func(it rtree.Item) bool {
			i := it.(*item)
			if i.poly.Contains(gp) {
				w.hits++
			}
			return true
		})
	}
}

func convert(p *shp.Polygon) *geo.Polygon {
	pts := make([]*geo.Point, len(p.Points))
	ps := p.Points
	for i := range ps {
		pts[i] = geo.NewPoint(ps[i].X, ps[i].Y)
	}
	return geo.NewPolygon(pts)
}

func loadShapefile(fname string) (*rtree.RTree, []*item, error) {
	reader, err := shp.Open(fname)
	if err != nil {
		return nil, nil, err
	}
	defer reader.Close()

	rt := rtree.New(nil)
	var items []*item

	for reader.Next() {
		_, s := reader.Shape()
		b := s.BBox()
		switch g := s.(type) {
		case *shp.Polygon:
			p := convert(g)
			it := &item{
				min:  []float64{b.MinX, b.MinY},
				max:  []float64{b.MaxX, b.MaxY},
				poly: p,
			}
			rt.Insert(it)
			items = append(items, it)
		default:
			return nil, nil, fmt.Errorf("Unknown %T", s)
		}
	}

	return rt, items, nil
}

func main() {

	var n int
	var w int
	var cpuProfile string

	flag.IntVar(&n, "n", 1000000, "Number of lookups")
	flag.IntVar(&w, "w", runtime.NumCPU(), "Number of workers")
	flag.StringVar(&cpuProfile, "cpu", "", "Dump pprof profile to file.")

	flag.Parse()

	if flag.NArg() < 1 {
		log.Fatalln("ERROR: missing shapefile")
	}

	if cpuProfile != "" {
		f, err := os.Create(cpuProfile)
		if err != nil {
			log.Fatal(err)
		}
		pprof.StartCPUProfile(f)
		defer pprof.StopCPUProfile()
	}

	start := time.Now()

	rt, items, err := loadShapefile(flag.Arg(0))
	if err != nil {
		log.Fatalf("ERROR: %v\n", err)
	}

	fmt.Printf("Loading took: %v\n", time.Since(start))
	fmt.Printf("item count: %d\n", rt.Count())

	start = time.Now()
	ws := make([]*worker, w)

	var wg sync.WaitGroup
	wg.Add(w)

	for i, m, rest := 0, n/w, n%w; i < w; i++ {
		k := m
		if rest > 0 {
			rest--
			k++
		}
		wrk := newWorker(rt)
		ws[i] = wrk
		go wrk.run(k, items, &wg)
	}

	wg.Wait()

	hits := 0
	for _, wrk := range ws {
		hits += wrk.hits
	}

	stop := time.Now()
	took := stop.Sub(start)

	// n / took = s / time.Second
	s := (float64(n) / float64(took)) * float64(time.Second)
	fmt.Printf("Took %v\n", took)
	fmt.Printf("Lookups per second: %.2f\n", s)
	fmt.Printf("Number of hits: %d\n", hits)
}

pip.go

// A simple experiment to check how fast point in polygon tests with
// a first level spatial index filter are.
//
// This is Free Software covered by the MIT license.
// See https://opensource.org/licenses/MIT for details.
// (c) 2017 by Sascha L. Teichmann (sascha.teichmann@intevation.de)
//
package main

import (
	"flag"
	"fmt"
	"log"
	"math/rand"
	"os"
	"runtime"
	"runtime/pprof"
	"sync"
	"time"

	"github.com/tidwall/rtree"

	shp "github.com/jonas-p/go-shp"

	geo "github.com/kellydunn/golang-geo"
)

type item struct {
	min  []float64
	max  []float64
	poly *geo.Polygon
}

func (it *item) Rect(interface{}) ([]float64, []float64) {
	return it.min, it.max
}

func randRange(rnd *rand.Rand, min, max float64) func() float64 {
	w := max - min
	return func() float64 {
		return min + w*rnd.Float64()
	}
}

type point struct{ *geo.Point }

func (p point) Rect(interface{}) ([]float64, []float64) {
	const eps = 0.01
	x, y := p.Lat(), p.Lng()
	return []float64{x - eps, y - eps}, []float64{x + eps, y + eps}
}

type worker struct {
	tree *rtree.RTree
	xr   func() float64
	yr   func() float64
	hits int
}

func newWorker(tree *rtree.RTree, bbox shp.Box) *worker {
	rnd := rand.New(rand.NewSource(time.Now().Unix()))
	return &worker{
		tree: tree,
		xr:   randRange(rnd, bbox.MinX, bbox.MaxX),
		yr:   randRange(rnd, bbox.MinY, bbox.MaxY),
	}
}

func (w *worker) run(n int, wg *sync.WaitGroup) {
	defer wg.Done()

	for i := 0; i < n; i++ {
		x, y := w.xr(), w.yr()
		gp := geo.NewPoint(x, y)
		w.tree.Search(point{gp}, func(it rtree.Item) bool {
			i := it.(*item)
			if i.poly.Contains(gp) {
				w.hits++
			}
			return true
		})
	}
}

func convert(p *shp.Polygon) *geo.Polygon {
	pts := make([]*geo.Point, len(p.Points))
	ps := p.Points
	for i := range ps {
		pts[i] = geo.NewPoint(ps[i].X, ps[i].Y)
	}
	return geo.NewPolygon(pts)
}

func loadShapefile(fname string) (*rtree.RTree, shp.Box, error) {
	reader, err := shp.Open(fname)
	if err != nil {
		return nil, shp.Box{}, err
	}
	defer reader.Close()

	bbox := reader.BBox()

	rt := rtree.New(nil)
	for reader.Next() {
		_, s := reader.Shape()
		b := s.BBox()
		switch g := s.(type) {
		case *shp.Polygon:
			p := convert(g)
			rt.Insert(&item{
				min:  []float64{b.MinX, b.MinY},
				max:  []float64{b.MaxX, b.MaxY},
				poly: p,
			})
		default:
			return nil, shp.Box{}, fmt.Errorf("Unknown %T", s)
		}
	}

	return rt, bbox, nil
}

func main() {

	var n int
	var w int
	var cpuProfile string

	flag.IntVar(&n, "n", 1000000, "Number of lookups")
	flag.IntVar(&w, "w", runtime.NumCPU(), "Number of workers")
	flag.StringVar(&cpuProfile, "cpu", "", "Dump pprof profile to file.")

	flag.Parse()

	if flag.NArg() < 1 {
		log.Fatalln("ERROR: missing shapefile")
	}

	if cpuProfile != "" {
		f, err := os.Create(cpuProfile)
		if err != nil {
			log.Fatal(err)
		}
		pprof.StartCPUProfile(f)
		defer pprof.StopCPUProfile()
	}

	start := time.Now()

	rt, bbox, err := loadShapefile(flag.Arg(0))
	if err != nil {
		log.Fatalf("ERROR: %v\n", err)
	}

	fmt.Printf("Loading took: %v\n", time.Since(start))
	fmt.Printf("item count: %d\n", rt.Count())

	start = time.Now()
	ws := make([]*worker, w)

	var wg sync.WaitGroup
	wg.Add(w)

	for i, m, rest := 0, n/w, n%w; i < w; i++ {
		k := m
		if rest > 0 {
			rest--
			k++
		}
		wrk := newWorker(rt, bbox)
		ws[i] = wrk
		go wrk.run(k, &wg)
	}

	wg.Wait()

	hits := 0
	for _, wrk := range ws {
		hits += wrk.hits
	}

	stop := time.Now()
	took := stop.Sub(start)

	// n / took = s / time.Second
	s := (float64(n) / float64(took)) * float64(time.Second)
	fmt.Printf("Took %v\n", took)
	fmt.Printf("Lookups per second: %.2f\n", s)
	fmt.Printf("Number of hits: %d\n", hits)
}