// Copyright (c) July 12, 2020 CareerMonk Publications and others.
// E-Mail           		: info@careermonk.com
// Creation Date    		: 2020-07-12 06:15:46
// Last modification		: 2020-07-12
//               by		: Narasimha Karumanchi
// Book Title			: Data Structures And Algorithmic Thinking With Go
// Warranty         		: This software is provided "as is" without any
// 				   warranty without even the implied warranty of
// 				    merchantability or fitness for a particular purpose.

package main

import (
	"fmt"
	"math"
	"math/rand"
	"time"
)

// number of points to search for closest pair
const n = 1e6

// size of bounding box for points.
// x and y will be random with uniform distribution in the range [0,scale).
const scale = 100.

// point struct
type xy struct {
	x, y float64 // coordinates
	key  int64   // an annotation used in the algorithm
}

func d(p1, p2 xy) float64 {
	return math.Hypot(p2.x-p1.x, p2.y-p1.y)
}

func closestPair(s []xy) (p1, p2 xy) {
	if len(s) < 2 {
		panic("2 points required")
	}
	var dxi float64
	// step 0
	for s1, i := s, 1; ; i++ {
		// step 1: compute min distance to a random point
		// (for the case of random data, it's enough to just try
		// to pick a different point)
		rp := i % len(s1)
		xi := s1[rp]
		dxi = 2 * scale
		for p, xn := range s1 {
			if p != rp {
				if dq := d(xi, xn); dq < dxi {
					dxi = dq
				}
			}
		}

		// step 2: filter
		invB := 3 / dxi             // b is size of a mesh cell
		mx := int64(scale*invB) + 1 // mx is number of cells along a side
		// construct map as a histogram:
		// key is index into mesh.  value is count of points in cell
		hm := map[int64]int{}
		for ip, p := range s1 {
			key := int64(p.x*invB)*mx + int64(p.y*invB)
			s1[ip].key = key
			hm[key]++
		}
		// construct s2 = s1 less the points without neighbors
		s2 := make([]xy, 0, len(s1))
		nx := []int64{-mx - 1, -mx, -mx + 1, -1, 0, 1, mx - 1, mx, mx + 1}
		for i, p := range s1 {
			nn := 0
			for _, ofs := range nx {
				nn += hm[p.key+ofs]
				if nn > 1 {
					s2 = append(s2, s1[i])
					break
				}
			}
		}

		// step 3: done?
		if len(s2) == 0 {
			break
		}
		s1 = s2
	}
	// step 4: compute answer from approximation
	invB := 1 / dxi
	mx := int64(scale*invB) + 1
	hm := map[int64][]int{}
	for i, p := range s {
		key := int64(p.x*invB)*mx + int64(p.y*invB)
		s[i].key = key
		hm[key] = append(hm[key], i)
	}
	nx := []int64{-mx - 1, -mx, -mx + 1, -1, 0, 1, mx - 1, mx, mx + 1}
	var min = scale * 2
	for ip, p := range s {
		for _, ofs := range nx {
			for _, iq := range hm[p.key+ofs] {
				if ip != iq {
					if d1 := d(p, s[iq]); d1 < min {
						min = d1
						p1, p2 = p, s[iq]
					}
				}
			}
		}
	}
	return p1, p2
}

func main() {
	rand.Seed(time.Now().Unix())
	points := make([]xy, n)
	for i := range points {
		points[i] = xy{rand.Float64() * scale, rand.Float64() * scale, 0}
	}
	p1, p2 := closestPair(points)
	fmt.Println(p1, p2)
	fmt.Println("distance:", d(p1, p2))
}
