Mercurial > gemma
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 }