view pkg/common/nashsutcliffe.go @ 3097:e6ba32b060df

Nash Sutcliffe: Treat last prediction date as valid afterwards. Small optimizsations.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Tue, 23 Apr 2019 18:10:19 +0200
parents cb3360653652
children 8f2ac24b0cb3
line wrap: on
line source

// 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) 2019 by via donau
//   – Österreichische Wasserstraßen-Gesellschaft mbH
// Software engineering by Intevation GmbH
//
// Author(s):
//  * Sascha L. Teichmann <sascha.teichmann@intevation.de>

package common

import (
	"fmt"
	"sort"
	"time"
)

type TimedValue struct {
	When  time.Time
	Value float64
}

type TimedValues []TimedValue

func (mvs TimedValues) Sort() {
	sort.Slice(mvs, func(i, j int) bool {
		return mvs[i].When.Before(mvs[j].When)
	})
}

var utc0 = time.Unix(0, 0).UTC()

func (mvs TimedValues) Interpolate(when time.Time) (float64, bool) {

	for lo, hi := 0, len(mvs)-2; lo <= hi; {
		mid := lo + (hi-lo)/2
		m1 := &mvs[mid]
		m2 := &mvs[mid+1]
		if m2.When.Before(when) {
			hi = mid - 1
			continue
		}
		if m1.When.After(when) {
			lo = mid + 1
			continue
		}
		if m1.When.Equal(when) {
			//log.Printf("matches first: %f\n", m1.Value)
			return m1.Value, true
		}
		if m2.When.Equal(when) {
			//log.Printf("matches second: %f\n", m1.Value)
			return m2.Value, true
		}

		// f(m1.When) = m1.Value
		// f(m2.When) = m2.Value
		// m1.Value = m1.When*a + b <=> b = m1.Value - m1.When*a
		// m2.Value = m2.When*a + b
		// m1.Value - m2.Value = a*(m1.When - m2-When)
		// a = (m1.Value - m2.Value)/(m1.When - m2.When) for m1.When != m2.When

		if m1.When.Equal(m2.When) {
			return (m1.Value + m2.Value) * 0.5, true
		}

		a := (m1.Value - m2.Value) / m1.When.Sub(m2.When).Seconds()
		b := m1.Value - m1.When.Sub(utc0).Seconds()*a
		m := when.Sub(utc0).Seconds()*a + b

		//log.Printf("%f %f %f\n", m1.Value, m, m2.Value)
		return m, true
	}

	if l := len(mvs); l > 0 {
		if when.Equal(mvs[0].When) {
			//log.Printf("at start\n")
			return mvs[0].Value, true
		}
		if !when.Before(mvs[l-1].When) {
			//log.Printf("after end\n")
		}
		return mvs[l-1].Value, true
	}

	//if len(mvs) > 0 {
	//	log.Printf("does not match %v %v %v\n",
	//		mvs[0].When, mvs[len(mvs)-1].When.Sub(mvs[0].When), when)
	//}

	return 0, false
}

func NashSutcliffe(predicted, observed []float64) float64 {

	if len(predicted) != len(observed) {
		panic(fmt.Sprintf(
			"NashSutcliffe: predicted and observed len differ: %d != %d",
			len(predicted),
			len(observed)))
	}

	if len(observed) == 0 {
		return 0
	}

	var mo float64
	for _, v := range observed {
		mo += v
	}

	mo /= float64(len(observed))

	var num, denom float64
	for i, o := range observed {
		d1 := predicted[i] - o
		num += d1 * d1
		d2 := o - mo
		denom += d2 * d2
	}

	return 1 - num/denom
}