// Copyright ©2017 The gonum Authors. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.

package fd

import "gonum.org/v1/gonum/mat"

// Watson implements the Watson's function.
// Dimension of the problem should be 2 <= dim <= 31. For dim == 9, the problem
// of minimizing the function is very ill conditioned.
//
// This is copied from gonum.org/v1/optimize/functions for testing Hessian-like
// derivative methods.
//
// References:
//  - Kowalik, J.S., Osborne, M.R.: Methods for Unconstrained Optimization
//    Problems. Elsevier North-Holland, New York, 1968
//  - More, J., Garbow, B.S., Hillstrom, K.E.: Testing unconstrained
//    optimization software. ACM Trans Math Softw 7 (1981), 17-41
type Watson struct{}

func (Watson) Func(x []float64) (sum float64) {
	for i := 1; i <= 29; i++ {
		d1 := float64(i) / 29

		d2 := 1.0
		var s1 float64
		for j := 1; j < len(x); j++ {
			s1 += float64(j) * d2 * x[j]
			d2 *= d1
		}

		d2 = 1.0
		var s2 float64
		for _, v := range x {
			s2 += d2 * v
			d2 *= d1
		}

		t := s1 - s2*s2 - 1
		sum += t * t
	}
	t := x[1] - x[0]*x[0] - 1
	sum += x[0]*x[0] + t*t
	return sum
}

func (Watson) Grad(grad, x []float64) {
	if len(x) != len(grad) {
		panic("incorrect size of the gradient")
	}

	for i := range grad {
		grad[i] = 0
	}
	for i := 1; i <= 29; i++ {
		d1 := float64(i) / 29

		d2 := 1.0
		var s1 float64
		for j := 1; j < len(x); j++ {
			s1 += float64(j) * d2 * x[j]
			d2 *= d1
		}

		d2 = 1.0
		var s2 float64
		for _, v := range x {
			s2 += d2 * v
			d2 *= d1
		}

		t := s1 - s2*s2 - 1
		s3 := 2 * d1 * s2
		d2 = 2 / d1
		for j := range x {
			grad[j] += d2 * (float64(j) - s3) * t
			d2 *= d1
		}
	}
	t := x[1] - x[0]*x[0] - 1
	grad[0] += x[0] * (2 - 4*t)
	grad[1] += 2 * t
}

func (Watson) Hess(hess mat.MutableSymmetric, x []float64) {
	dim := len(x)
	if dim != hess.Symmetric() {
		panic("incorrect size of the Hessian")
	}

	for j := 0; j < dim; j++ {
		for k := j; k < dim; k++ {
			hess.SetSym(j, k, 0)
		}
	}
	for i := 1; i <= 29; i++ {
		d1 := float64(i) / 29
		d2 := 1.0
		var s1 float64
		for j := 1; j < dim; j++ {
			s1 += float64(j) * d2 * x[j]
			d2 *= d1
		}

		d2 = 1.0
		var s2 float64
		for _, v := range x {
			s2 += d2 * v
			d2 *= d1
		}

		t := s1 - s2*s2 - 1
		s3 := 2 * d1 * s2
		d2 = 2 / d1
		th := 2 * d1 * d1 * t
		for j := 0; j < dim; j++ {
			v := float64(j) - s3
			d3 := 1 / d1
			for k := 0; k <= j; k++ {
				hess.SetSym(k, j, hess.At(k, j)+d2*d3*(v*(float64(k)-s3)-th))
				d3 *= d1
			}
			d2 *= d1
		}
	}
	t1 := x[1] - x[0]*x[0] - 1
	hess.SetSym(0, 0, hess.At(0, 0)+8*x[0]*x[0]+2-4*t1)
	hess.SetSym(0, 1, hess.At(0, 1)-4*x[0])
	hess.SetSym(1, 1, hess.At(1, 1)+2)
}