package main

import (
	"bufio"
	"compress/gzip"
	"encoding/binary"
	"fmt"
	"io"
	"math"
	"os"
	"strconv"
	"time"

	xmlparser "github.com/tamerh/xml-stream-parser"
	"github.com/tkrajina/gpxgo/gpx"
)

const srtmResolution = 3600 // udává velikost jednoho čtverce (máme čtverce 3601x3601 bodů)
const minElevation = 1      // kolik metrů na 100m délky musí cesta klesnout, aby byla "klesající"

// Jeden bod z OSM
type node struct {
	id  int
	lat float64
	lon float64
}

// Struktury pro DAG (acyklický graf)
type edge struct {
	wayID    int
	wayNodes []*node
	from     *vertex
	to       *vertex
	length   float64
}
type vertex struct {
	node         *node
	altitude     float64
	outgoingWays []*edge
	// Proměnné pro počítání a vypsání nejdelší cesty v DAGu
	inDegree         int
	longestWayLength float64
	longestWayFrom   *edge
}

// Funkce na počítání vzdálenosti
func haversineDistance(point1, point2 *node) float64 {
	R := float64(6378137) // Poloměr Země v metrech (na rovníku)

	// Převedeme na radiány (vynásobením pi/180)
	phi1 := point1.lat * math.Pi / 180
	phi2 := point2.lat * math.Pi / 180
	dphi := (point2.lat - point1.lat) * math.Pi / 180
	dlambda := (point2.lon - point1.lon) * math.Pi / 180

	a := math.Pow(math.Sin(dphi/2), 2) + math.Cos(phi1)*math.Cos(phi2)*math.Pow(math.Sin(dlambda/2), 2)

	return 2 * R * math.Atan2(math.Sqrt(a), math.Sqrt(1-a))
}

// Funkce na spočítání délky cesty pomocí haversineDistance
func computeLength(way []*node) float64 {
	var lastNode *node
	length := 0.0
	for _, nd := range way {
		if lastNode != nil {
			length += haversineDistance(lastNode, nd)
		}
		lastNode = nd
	}
	return length
}

////////////////////////////////////////////////////////////////////////////////

// SRTM věci
type srtmArea struct {
	latStart, lonStart float64
	stepsize           float64
	points             [][]int16 // výšky indexované lat,lon
}

// Funkce na načtení SRTM (z několika čtverců) do objektu srtmArea
func loadSRTM(minLat, maxLat, minLon, maxLon float64) (*srtmArea, error) {
	latStart, latStop := int(minLat), int(maxLat)
	lonStart, lonStop := int(minLon), int(maxLon)
	// Spočítáme celkový počet řádek - každý blok je 3601x3601, ale překrývají se,
	// takže třeba dva bloky vedle sebe mají šířku 7201
	rows := (latStop-latStart+1)*srtmResolution + 1
	columns := (lonStop-lonStart+1)*srtmResolution + 1
	area := &srtmArea{
		latStart: float64(latStart),
		lonStart: float64(lonStart),
		points:   make([][]int16, rows),
		stepsize: 1.0 / srtmResolution,
	}
	for i := range area.points {
		area.points[i] = make([]int16, columns)
	}

	for latSquare := latStart; latSquare <= latStop; latSquare++ {
		for lonSquare := lonStart; lonSquare <= lonStop; lonSquare++ {
			latOffset := (latSquare - latStart) * srtmResolution
			lonOffset := (lonSquare - lonStart) * srtmResolution

			// Otevřeme soubor s blokem
			filename := fmt.Sprintf("N%02dE%03d.hgt", latSquare, lonSquare)
			fmt.Fprintf(os.Stderr, "Loading SRTM file %s\n", filename)
			file, err := os.Open(filename)
			if err != nil {
				return nil, fmt.Errorf("Cannot open file '%s': %v", filename, err)
			}

			// Data ze souboru čteme od severu k jihu, ale nám se více hodí, aby
			// y=0 odpovídalo spodku čtverce, y tedy jedeme odzadu
			for y := 3600; y >= 0; y-- {
				// Načteme 3601 hodnot, zapisujeme do velké tabulky na správné místo
				// (podle offsetu tohoto načítaného čtverce)
				err := binary.Read(file, binary.BigEndian, area.points[latOffset+y][lonOffset:lonOffset+srtmResolution+1])
				if err != nil {
					return nil, err
				}
			}
		}
	}
	return area, nil
}

// Funkce na spočítání výšky bodu aproximací okolních SRTM dat
func (a *srtmArea) getAltitude(point *node) float64 {
	latIndex := (point.lat - a.latStart) / a.stepsize
	lonIndex := (point.lon - a.lonStart) / a.stepsize

	// Spočítáme si okolní mřížkové body
	// (nevadí nám, když některé budou stejné)
	latA := int(latIndex)
	latB := int(math.Ceil(latIndex))
	lonA := int(lonIndex)
	lonB := int(math.Ceil(lonIndex))

	topL := float64(a.points[latA][lonA])
	topR := float64(a.points[latA][lonB])
	bottomL := float64(a.points[latB][lonA])
	bottomR := float64(a.points[latB][lonB])

	// Spočítáme si průměr left-right podle lon
	// (tím čtverec splácneme do úsečky top-bottom)
	pomer := lonIndex - float64(lonA)
	top := topL + (topR-topL)*pomer
	bottom := bottomL + (bottomR-bottomL)*pomer

	// Spočítáme si průměr top-bottom podle lat
	return top + (bottom-top)*(latIndex-float64(latA))
}

////////////////////////////////////////////////////////////////////////////////

// Funkce na kontrolu jestli má element daný tag
func hasTag(el *xmlparser.XMLElement, name string) bool {
	for _, tag := range el.Childs["tag"] {
		if tag.Attrs["k"] == name {
			return true
		}
	}
	return false
}

// Funkce na kontrolu jestli má tag elementu danou hodnotu
func hasTagValue(el *xmlparser.XMLElement, name string, value string) bool {
	for _, tag := range el.Childs["tag"] {
		if tag.Attrs["k"] == name {
			return tag.Attrs["v"] == value
		}
	}
	return false
}


func die(exitcode int, format string, a ...interface{}) {
	fmt.Fprintf(os.Stderr, format, a...)
	os.Exit(1)
}

func main() {
	if len(os.Args) != 3 {
		die(1, "Usage: %s <file.xml.gz> <output.gpx>\n", os.Args[0])
	}
	// 1. Otevřeme OSM soubor pro vstup a GPX soubor pro výstup
	file, err := os.Open(os.Args[1])
	if err != nil {
		die(1, "Cannot open file: %v\n", err)
	}
	defer file.Close()
	gpxFile, err := os.Create(os.Args[2])
	if err != nil {
		die(1, "Cannot open gpx file for write: %v\n", err)
	}
	defer gpxFile.Close()

	// 2. Připravíme si funkci na resetování readeru, protože budeme muset XML soubor projít vícekrát
	var archiveReader io.Reader
	var parser *xmlparser.XMLParser
	resetReader := func(elements ...string) {
		if _, err := file.Seek(0, 0); err != nil {
			die(1, "Cannot seek to beginning of the file")
		}
		// Pořídíme si gzip reader
		if archiveReader, err = gzip.NewReader(file); err != nil {
			die(1, "Cannot open the gzip reader: %v", err)
		}
		// Dopustíme se triku - odzipování dat spustíme v samostatném vlákně
		// a pošleme si výsledek skrz rouru (pipe), zrychlí nás to v případě
		// gzipu o zhruba 10-20%.
		r, w := io.Pipe()
		go func() {
			io.Copy(w, archiveReader)
			w.Close()
		}()
		parser = xmlparser.NewXMLParser(bufio.NewReaderSize(r, 1024*1024*256), elements...)
	}

	ways := map[int][]*node{}
	wayNodes := map[int]*node{}
	var foundNodesCounter int

	start := time.Now()
	phaseStart := start
	// Funkce na vypisování statistiky, spustíme ji tentokrát na pozadí
	printStat := func() {
		megabytes := float64(parser.TotalReadSize) / 1024 / 1024
		elapsed := time.Now().Sub(phaseStart)
		fmt.Fprintf(os.Stderr,
			"Elapsed: %-20s Read: %10.2fMB (%.2fMB/s)\tWays: %10d\tWay nodes: %10d / %-10d\r",
			elapsed, megabytes, megabytes/elapsed.Seconds(), len(ways), foundNodesCounter, len(wayNodes),
		)
	}
	statTicker := time.NewTicker(500 * time.Millisecond)
	go func() {
		for range statTicker.C {
			printStat()
		}
	}()

	////////////////////////////////////////////////////////////////////////
	// PRVNÍ FÁZE - nalezení všech cest ////////////////////////////////////

	fmt.Fprintf(os.Stderr, "[Phase 1] Getting all ways\n")

	// 1. Pořídíme si reader od začátku souboru a budeme chtít teď pouze cesty
	resetReader("way", "node")
	parser.SkipElements([]string{"relation"}).SkipOuterElements()

	lonMin, lonMax := math.Inf(1), math.Inf(-1)
	latMin, latMax := math.Inf(1), math.Inf(-1)

	for xml := range parser.Stream() {
		switch xml.Name {
		case "way":
			// Chceme highway, ale pro pěknější cesty vyloučíme dálnice a víceproudé silnice
			if hasTag(xml, "highway") && !hasTagValue(xml, "highway", "motorway") && !hasTagValue(xml, "highway", "trunk") {
				wayID, _ := strconv.Atoi(xml.Attrs["id"])
				ways[wayID] = []*node{}
				for _, nd := range xml.Childs["nd"] {
					nodeID, _ := strconv.Atoi(nd.Attrs["ref"])
					// Každý node si pamatujeme jenom jednou a z cesty na něj odkazujeme
					// (to nám poté pomůže při plnění vlastností nodů)
					if _, found := wayNodes[nodeID]; !found {
						wayNodes[nodeID] = &node{id: nodeID}
					}
					ways[wayID] = append(ways[wayID], wayNodes[nodeID])
				}
			}
		}
	}
	printStat()
	fmt.Fprintf(os.Stderr, "\nFound %d ways\n", len(ways))

	////////////////////////////////////////////////////////////////////////
	// DRUHÁ FÁZE - načtení bodů cest //////////////////////////////////////

	fmt.Fprintf(os.Stderr, "\n[Phase 2] Getting all nodes of the ways\n")

	// 1. Resetujeme se na začátek souboru a připravíme se na další průchod
	resetReader("node")
	parser.SkipElements([]string{"way", "relation"}).SkipOuterElements()

	// 2. Projdeme všechny nodes a zaznamenáme si ty, které chceme
	phaseStart = time.Now()
	for xml := range parser.Stream() {
		switch xml.Name {
		case "node":
			id, _ := strconv.Atoi(xml.Attrs["id"])
			if nd, found := wayNodes[id]; found {
				nd.lat, _ = strconv.ParseFloat(xml.Attrs["lat"], 64)
				nd.lon, _ = strconv.ParseFloat(xml.Attrs["lon"], 64)
				lonMin = math.Min(lonMin, nd.lon)
				lonMax = math.Max(lonMax, nd.lon)
				latMin = math.Min(latMin, nd.lat)
				latMax = math.Max(latMax, nd.lat)
				foundNodesCounter++
				if foundNodesCounter == len(wayNodes) {
					break
				}
			}
		}
	}

	printStat()
	statTicker.Stop()
	fmt.Fprintf(os.Stderr, "\nCoordinates of %d nodes found\n", foundNodesCounter)
	fmt.Fprintf(os.Stderr, "Min: [%f, %f]\tMax: [%f,%f]\n", latMin, lonMin, latMax, lonMax)

	////////////////////////////////////////////////////////////////////////
	// TŘETÍ FÁZE - načtení SRTM dat ///////////////////////////////////////

	fmt.Fprintf(os.Stderr, "\n[Phase 3] Loading SRTM\n")

	srtm, err := loadSRTM(latMin, latMax, lonMin, lonMax)
	if err != nil {
		die(1, "Cannot load SRTM data: %v", err)
	}

	////////////////////////////////////////////////////////////////////////
	// ČTVRTÁ FÁZE - filtrace cest a sestavení grafu (DAGu) ////////////////

	fmt.Fprintf(os.Stderr, "\n[Phase 4] Constructing DAG from %d candidate ways\n", len(ways))

	vertices := map[int]*vertex{}
	edgeCounter := 0

	// Profiltrujeme a vytvoříme z nich hrany našeho grafu
	for wayID, way := range ways {
		// Musíme zkontrolovat, že cesta buď klesá nebo stoupá a že konce cesty se liší alespoň o danou elevaci (viz konstanta nahoře)
		length := computeLength(way)
		diff := srtm.getAltitude(way[0]) - srtm.getAltitude(way[len(way)-1])
		if math.Abs(diff) < length/100*minElevation {
			continue // cesta "není z kopce" -> nechceme ji
		} else if diff < 0 {
			// cesta "stoupá", otočíme ji
			for left, right := 0, len(way)-1; left < right; left, right = left+1, right-1 {
				way[left], way[right] = way[right], way[left]
			}
		}

		// Předpokládejme, že cesta "klesá", potřebujeme to ověřit
		last := srtm.getAltitude(way[0])
		ok := true
		for _, node := range way {
			current := srtm.getAltitude(node)
			if last < current {
				ok = false
				break // lokální kopec
			}
			last = current
		}
		if !ok {
			continue
		}
		// Přidáme klesající hranu do grafu
		fromNode := way[0]
		toNode := way[len(way)-1]
		fromVertex, found := vertices[fromNode.id]
		if !found {
			fromVertex = &vertex{
				node:             fromNode,
				altitude:         srtm.getAltitude(fromNode),
				outgoingWays:     []*edge{},
			}
			vertices[fromNode.id] = fromVertex
		}
		toVertex, found := vertices[toNode.id]
		if !found {
			toVertex = &vertex{
				node:             toNode,
				altitude:         srtm.getAltitude(toNode),
				outgoingWays:     []*edge{},
			}
			vertices[toNode.id] = toVertex
		}
		fromVertex.outgoingWays = append(fromVertex.outgoingWays, &edge{
			wayID:    wayID,
			wayNodes: way,
			from:     fromVertex,
			to:       toVertex,
			length:   length,
		})
		toVertex.inDegree++
		edgeCounter++
	}

	fmt.Fprintf(os.Stderr, "\nConstructed DAG with %d vertices and %d edges\n", len(vertices), edgeCounter)

	////////////////////////////////////////////////////////////////////////
	// PÁTÁ FÁZE - hledání nejdelší cesty v DAGu ///////////////////////////

	fmt.Fprintf(os.Stderr, "\n[Phase 5] Finding longest way\n")

	// 1. Najdeme vrcholy, do kterých nic nevede
	verticesQueue := []*vertex{}
	for _, v := range vertices {
		if v.inDegree == 0 {
			verticesQueue = append(verticesQueue, v)
		}
	}

	maxVertex := verticesQueue[0] // budeme si pamatovat vrchol ve kterém končí nejdelší nalezená cesta

	// 2. Zpracujeme postupně všechny vrcholy
	for true {
		if len(verticesQueue) == 0 {
			break
		}
		v := verticesQueue[0]
		verticesQueue = verticesQueue[1:]

		// Check nejdelší zatím nalezené cesty
		if maxVertex.longestWayLength < v.longestWayLength {
			maxVertex = v
		}

		// Zpracujeme všechny výstupní hrany
		for _, e := range v.outgoingWays {
			if e.to.longestWayLength < v.longestWayLength+e.length {
				e.to.longestWayLength = v.longestWayLength + e.length
				e.to.longestWayFrom = e
			}
			e.to.inDegree-- // snížíme počet nezpracovaných vstupních hran
			if e.to.inDegree == 0 {
				verticesQueue = append(verticesQueue, e.to)
			}
		}
	}

	// 3. Vypíšeme nejdelší nalezenou cestu (+ uložíme gpx)
	fmt.Fprintf(os.Stderr, "\nLongest way length: %f m\n", maxVertex.longestWayLength)
	finalEdges := []*edge{}
	v := maxVertex
	for v.longestWayFrom != nil {
		finalEdges = append(finalEdges, v.longestWayFrom)
		v = v.longestWayFrom.from
	}

	// Samotný výpis
	fmt.Fprintf(os.Stderr, "Start: https://www.openstreetmap.org/node/%d\tAltitude: %f m\n", v.node.id, v.altitude)
	segments := []gpx.GPXTrackSegment{}
	for i := len(finalEdges) - 1; i >= 0; i-- {
		e := finalEdges[i]
		fmt.Fprintf(os.Stderr, "Way:  https://www.openstreetmap.org/way/%d\tLength: %f m\n", e.wayID, e.length)
		fmt.Fprintf(os.Stderr, "Node: https://www.openstreetmap.org/node/%d\tAltitude: %f m\n", e.to.node.id, e.to.altitude)

		points := []gpx.GPXPoint{}
		for _, node := range e.wayNodes {
			p := gpx.GPXPoint{
				Point: gpx.Point{
					Latitude:  node.lat,
					Longitude: node.lon,
				},
			}
			p.Elevation.SetValue(srtm.getAltitude(node)) // oklika na nullable typ
			points = append(points, p)
		}
		segments = append(segments, gpx.GPXTrackSegment{Points: points})
	}

	// 4. Uložíme GPX
	gpxStruct := gpx.GPX{
		Name:    "Longest way",
		Creator: "KSP",
		Tracks: []gpx.GPXTrack{gpx.GPXTrack{
			Name:     "Longest way",
			Segments: segments,
		}},
	}
	gpxBytes, err := gpxStruct.ToXml(gpx.ToXmlParams{Version: "1.1", Indent: true})
	if err != nil {
		fmt.Fprintf(os.Stderr, "Error during converting GPX: %v", err)
	}
	gpxFile.Write(gpxBytes)
}