// This is Free Software under GNU Affero General Public License v >= 3.0
// without warranty, see README.md and license for details.
//
// SPDX-License-Identifier: AGPL-3.0-or-later
// License-Filename: LICENSES/AGPL-3.0.txt
//
// Copyright (C) 2018 by via donau
//   – Österreichische Wasserstraßen-Gesellschaft mbH
// Software engineering by Intevation GmbH
//
// Author(s):
//  * Sascha L. Teichmann <sascha.teichmann@intevation.de>

package main

import (
	"context"
	"database/sql"
	"errors"
	"flag"
	"fmt"
	"log"
	"runtime"
	"sync"
	"time"

	"gemma.intevation.de/gemma/pkg/common"
	"gemma.intevation.de/gemma/pkg/octree"
)

var (
	bottleneck = flag.String("bottleneck", "", "name of the bottleneck")
	first      = flag.String("first", "", "date of the first sounding result")
	second     = flag.String("second", "", "date of the second sounding result")
)

const contourTolerance = 0.1

const (
	loadOctreeSQL = `
SELECT sr.octree_index
FROM waterway.sounding_results sr JOIN waterway.bottlenecks bn
  ON sr.bottleneck_id = bn.id
WHERE bn.bottleneck_id = $1 AND sr.date_info = $2::date
  AND sr.octree_index IS NOT NULL`

	insertContourSQL = `
INSERT INTO diff_contour_lines (
  height,
  lines
)
SELECT
  $4,
  ST_Transform(
    ST_Multi(
	  ST_CollectionExtract(
		ST_SimplifyPreserveTopology(
		  ST_Multi(ST_Collectionextract(
			ST_MakeValid(ST_GeomFromWKB($1, $2::integer)), 2)),
		  $3
		),
		2
	  )
	),
    4326
  )
`
	clippingPolygonSQL = `
WITH joined AS (
  SELECT
    sr.area      AS area,
    sr.date_info AS date_info
  FROM waterway.sounding_results sr JOIN
     waterway.bottlenecks bn ON sr.bottleneck_id = bn.id
  WHERE bn.bottleneck_id = $1
)
SELECT ST_AsBinary(
  ST_intersection(
    (SELECT ST_Transform(area::geometry, $2::int) FROM joined WHERE date_info = $3::date),
    (SELECT ST_Transform(area::geometry, $2::int) FROM joined WHERE date_info = $4::date)
  )
  ) AS area
`
	insertContourSQLClipped = `
WITH joined AS (
  SELECT
    sr.area    AS area,
    sr.date_info AS date_info
  FROM waterway.sounding_results sr JOIN
     waterway.bottlenecks bn ON sr.bottleneck_id = bn.id
  WHERE bn.bottleneck_id = $1
),
inter AS (
  SELECT
  ST_Buffer(
    ST_intersection(
      (SELECT ST_Transform(area::geometry, $2::int) FROM joined WHERE date_info = $3::date),
      (SELECT ST_Transform(area::geometry, $2::int) FROM joined WHERE date_info = $4::date)
    ),
    0.001) AS area
)
INSERT INTO diff_contour_lines (
  height,
  lines
)
SELECT
  $5,
  ST_Transform(
    ST_Multi(
	  ST_intersection(
		  ST_CollectionExtract(
			ST_SimplifyPreserveTopology(
			  ST_Multi(ST_Collectionextract(
				ST_MakeValid(ST_GeomFromWKB($6, $2::integer)), 2)),
			  $7
			),
			2
		  ),
		area
	  )
	),
    4326
  )
  FROM inter`
)

func check(err error) {
	if err != nil {
		log.Fatalf("error: %v\n", err)
	}
}

type point struct {
	x float64
	y float64
}

type pointMap map[point]float64

func (pm pointMap) triangulate() (*octree.Triangulation, error) {
	vertices := make([]octree.Vertex, len(pm))
	var i int
	for p, z := range pm {
		vertices[i] = octree.Vertex{X: p.x, Y: p.y, Z: z}
		i++
	}
	return octree.Triangulate(vertices)
}

func sliceWork(
	vs []octree.Vertex,
	dst pointMap,
	fn func([]octree.Vertex, func([]octree.Vertex) []octree.Vertex),
) {
	n := runtime.NumCPU()

	wg := new(sync.WaitGroup)

	slices := make(chan []octree.Vertex)
	out := make(chan []octree.Vertex)

	pool := make(chan []octree.Vertex, n)

	const pageSize = 2048

	turn := func(p []octree.Vertex) []octree.Vertex {
		if p != nil {
			out <- p
		}
		select {
		case p = <-pool:
		default:
			p = make([]octree.Vertex, 0, pageSize)
		}
		return p
	}

	for i := 0; i < n; i++ {
		wg.Add(1)
		go func() {
			defer wg.Done()
			for slice := range slices {
				fn(slice, turn)
			}
		}()
	}
	done := make(chan struct{})
	go func() {
		defer close(done)
		for s := range out {
			for i := range s {
				v := &s[i]
				key := point{x: v.X, y: v.Y}
				if z, found := dst[key]; found {
					dst[key] = (z + v.Z) * 0.5
				} else {
					dst[key] = v.Z
				}
			}
			select {
			case pool <- s[:0:pageSize]:
			default:
			}
		}
	}()

	size := len(vs)/n + 1
	for len(vs) > 0 {
		var l int
		if len(vs) < size {
			l = len(vs)
		} else {
			l = size
		}
		slices <- vs[:l]
		vs = vs[l:]
	}
	close(slices)
	wg.Wait()
	close(out)
	<-done
}

func process(bottleneck string, firstDate, secondDate time.Time) error {
	start := time.Now()
	defer func() { log.Printf("processing took %v\n", time.Since(start)) }()

	ctx := context.Background()

	return run(func(db *sql.DB) error {

		conn, err := db.Conn(ctx)
		if err != nil {
			return err
		}
		defer conn.Close()

		tx, err := conn.BeginTx(ctx, nil)
		if err != nil {
			return err
		}
		defer tx.Rollback()

		type load struct {
			date time.Time
			data []byte
			err  *error
			dst  **octree.Tree
		}

		out := make(chan *load)
		wg := new(sync.WaitGroup)

		n := runtime.NumCPU()
		if n > 2 {
			n = 2
		}

		for i := 0; i < n; i++ {
			wg.Add(1)
			go func() {
				defer wg.Done()
				for l := range out {
					if *l.err == nil {
						*l.dst, *l.err = octree.Deserialize(l.data)
					}
				}
			}()
		}

		var firstErr, secondErr error
		var first, second *octree.Tree

		for _, l := range []*load{
			{date: firstDate, dst: &first, err: &firstErr},
			{date: secondDate, dst: &second, err: &secondErr},
		} {
			var data []byte
			if err := tx.QueryRowContext(
				ctx,
				loadOctreeSQL,
				bottleneck,
				l.date,
			).Scan(&data); err != nil {
				*l.err = err
			} else {
				l.data = data
			}
			out <- l
		}
		close(out)

		wg.Wait()

		if firstErr != nil || secondErr != nil {
			if firstErr != nil && secondErr != nil {
				return fmt.Errorf("%v, %v", firstErr, secondErr)
			}
			if firstErr != nil {
				return firstErr
			}
			return secondErr
		}

		if first.EPSG != second.EPSG {
			return errors.New("EPSG codes mismatch. Needs transformation slow pass.")
		}

		now := time.Now()
		log.Printf("loading took %v\n", now.Sub(start))
		last := now

		firstVs, secondVs := first.Vertices(), second.Vertices()

		result := make(pointMap, len(firstVs)+len(secondVs))

		sliceWork(
			firstVs,
			result,
			func(
				slice []octree.Vertex,
				turn func([]octree.Vertex) []octree.Vertex,
			) {
				p := turn(nil)
				for i := range slice {
					v := &slice[i]
					if z, found := second.Value(v.X, v.Y); found {
						p = append(p, octree.Vertex{v.X, v.Y, v.Z - z})
						if len(p) == cap(p) {
							p = turn(p)
						}
					}
				}
				if len(p) > 0 {
					turn(p)
				}
			})

		sliceWork(
			secondVs,
			result,
			func(
				slice []octree.Vertex,
				turn func([]octree.Vertex) []octree.Vertex,
			) {
				p := turn(nil)
				for i := range slice {
					v := &slice[i]
					if z, found := first.Value(v.X, v.Y); found {
						p = append(p, octree.Vertex{v.X, v.Y, z - v.Z})
						if len(p) == cap(p) {
							p = turn(p)
						}
					}
				}
				if len(p) > 0 {
					turn(p)
				}
			})

		now = time.Now()
		log.Printf("setting in took %v\n", now.Sub(last))
		last = now
		log.Printf("num points: %d\n", len(result))

		var clip []byte

		if err := tx.QueryRowContext(
			ctx, clippingPolygonSQL,
			bottleneck,
			first.EPSG,
			firstDate, secondDate,
		).Scan(&clip); err != nil {
			return err
		}

		var clippingPolygon octree.Polygon
		if err := clippingPolygon.FromWKB(clip); err != nil {
			return err
		}
		clippingPolygon.Indexify()

		now = time.Now()
		log.Printf("loading clipping polygon took %v\n", now.Sub(last))
		last = now

		tri, err := result.triangulate()
		if err != nil {
			return err
		}

		now = time.Now()
		log.Printf("triangulation took %v\n", now.Sub(last))
		last = now

		builder := octree.NewBuilder(tri.Tin())
		builder.Build()

		now = time.Now()
		log.Printf("building octree took %v\n", now.Sub(last))
		last = now

		tree := builder.Tree()
		tree.Clip(&clippingPolygon)
		//tree.Dump()

		now = time.Now()
		log.Printf("clipping octree took %v\n", now.Sub(last))
		last = now

		log.Printf("min/max: %f %f\n", tree.Min.Z, tree.Max.Z)

		var heights []float64

		switch {
		case tree.Min.Z >= 0: // All values positive.
			for v := 0.0; v <= tree.Max.Z; v += 0.1 {
				if v >= tree.Min.Z {
					heights = append(heights, v)
				}
			}
		case tree.Max.Z <= 0: // All values negative.
			for v := 0.0; v >= tree.Min.Z; v -= 0.1 {
				if v <= tree.Max.Z {
					heights = append(heights, v)
				}
			}
		default: // Positive and negative.
			for v := 0.1; v <= tree.Max.Z; v += 0.1 {
				heights = append(heights, v)
			}
			for i, j := 0, len(heights)-1; i < j; i, j = i+1, j-1 {
				heights[i], heights[j] = heights[j], heights[i]
			}
			for v := 0.0; v >= tree.Min.Z; v -= 0.1 {
				heights = append(heights, v)
			}
		}

		log.Printf("num heights: %d\n", len(heights))

		var dataSize int

		stmt, err := tx.PrepareContext(ctx, insertContourSQL)
		if err != nil {
			return err
		}
		defer stmt.Close()

		log.Println("do contours")
		octree.DoContours(tree, heights, func(res *octree.ContourResult) {
			if err == nil && len(res.Lines) > 0 {
				log.Printf("%.2f: lines: %d\n", res.Height, len(res.Lines))
				wkb := res.Lines.AsWKB2D()
				dataSize += len(wkb)
				_, err = stmt.ExecContext(
					ctx,
					wkb,
					first.EPSG,
					contourTolerance,
					res.Height,
				)
			}
		})

		if err != nil {
			return err
		}

		now = time.Now()
		log.Printf("Number of iso lines: %d\n", len(heights))
		log.Printf("Total WKB size: %.2fMB\n", float64(dataSize)/(1024*1024))
		log.Printf("generating iso lines took %v\n", now.Sub(last))
		last = now

		return tx.Commit()
	})
}

func main() {

	flag.Parse()

	firstDate, err := time.Parse(common.DateFormat, *first)
	check(err)
	secondDate, err := time.Parse(common.DateFormat, *second)
	check(err)

	if *bottleneck == "" {
		log.Fatalln("Missing bottleneck name")
	}

	check(process(*bottleneck, firstDate, secondDate))
}