func (k *Sum) Diag(X mat.Matrix) *mat.DiagDense {
	K1 := k.K1.Diag(X)
	K2 := k.K2.Diag(X)
	nx, _ := K1.Dims()
	K1r, K2r := K1.RawBand(), K2.RawBand()
	for i, k1i, k2i := 0, 0, 0; i < nx; i, k1i, k2i = i+1, k1i+K1r.Stride, k2i+K2r.Stride {
		K1r.Data[k1i] += K2r.Data[k2i]
	}
	return K1
}
Diag returns the diagonal of the kernel k(X, X)
package kernels
import (
	"fmt"
	"github.com/pa-m/sklearn/base"
	"math"
	"gonum.org/v1/gonum/floats"
	"gonum.org/v1/gonum/mat"
	t "github.com/pa-m/sklearn/gaussian_process/tensor"
)
type hyperparameter struct {
	Name    string
	PValue  *float64
	PBounds *[2]float64
}
func (param hyperparameter) IsFixed() bool {
	return param.PBounds == nil || (*param.PBounds)[1] == (*param.PBounds)[0]
}
type hyperparameters []hyperparameter
func (params hyperparameters) notFixed() hyperparameters {
	notFixed := hyperparameters{}
	for _, p := range params {
		if !p.IsFixed() {
			notFixed = append(notFixed, p)
		}
	}
	return notFixed
}
func (params hyperparameters) Theta() mat.Matrix {
	notFixed := params.notFixed()
	return matFromFunc{
		r: len(notFixed),
		c: 1,
		at: func(i, j int) float64 {
			return math.Log(*notFixed[i].PValue)
		},
		set: func(i, j int, v float64) {
			*notFixed[i].PValue = math.Exp(v)
		},
	}
}
func (params hyperparameters) Bounds() (t mat.Matrix) {
	notFixed := params.notFixed()
	return matFromFunc{
		r: len(notFixed),
		c: 2,
		at: func(i, j int) float64 {
			return math.Log((*notFixed[i].PBounds)[j])
		},
		set: func(i, j int, v float64) {
			(*notFixed[i].PBounds)[j] = math.Exp(v)
		},
	}
}
type Kernel interface {
	hyperparameters() hyperparameters
	Theta() mat.Matrix
	Bounds() mat.Matrix
	CloneWithTheta(theta mat.Matrix) Kernel
	Eval(X, Y mat.Matrix, evalGradient bool) (*mat.Dense, *t.Dense)
	Diag(X mat.Matrix) (K *mat.DiagDense)
	IsStationary() bool
	String() string
}
type StationaryKernelMixin struct{}
func (StationaryKernelMixin) IsStationary() bool { return true }
type NormalizedKernelMixin struct{}
func (k NormalizedKernelMixin) Diag(X mat.Matrix) (K *mat.DiagDense) {
	nx, _ := X.Dims()
	K = mat.NewDiagDense(nx, nil)
	for ix := 0; ix < nx; ix++ {
		K.SetDiag(ix, 1)
	}
	return
}
type KernelOperator struct {
	K1, K2 Kernel
}
func (k KernelOperator) hyperparameters() hyperparameters {
	return append(k.K1.hyperparameters(), k.K2.hyperparameters()...)
}
func (k KernelOperator) Theta() mat.Matrix {
	return k.hyperparameters().Theta()
}
func (k KernelOperator) Bounds() mat.Matrix {
	return k.hyperparameters().Bounds()
}
func (k KernelOperator) CloneWithTheta(theta mat.Matrix) Kernel {
	if theta == mat.Matrix(nil) {
		return KernelOperator{K1: k.K1.CloneWithTheta(nil), K2: k.K2.CloneWithTheta(nil)}
	}
	var td = base.ToDense(theta)
	n1, _ := k.K1.Theta().Dims()
	n2, _ := k.K2.Theta().Dims()
	return KernelOperator{
		K1: k.K1.CloneWithTheta(td.Slice(0, n1, 0, 1)),
		K2: k.K2.CloneWithTheta(td.Slice(n1, n1+n2, 0, 1)),
	}
}
func (k KernelOperator) Eval(X, Y mat.Matrix, evalGradient bool) (*mat.Dense, *t.Dense) {
	panic("Eval must be implemented by wrapper")
}
func (k KernelOperator) Diag(X mat.Matrix) *mat.DiagDense {
	panic("Diag must be implemented by wrapper")
}
func (k KernelOperator) String() string {
	panic("Diag must be implemented by wrapper")
}
func (k KernelOperator) IsStationary() bool {
	return k.K1.IsStationary() && k.K2.IsStationary()
}
type Sum struct {
	KernelOperator
}
func (k Sum) CloneWithTheta(theta mat.Matrix) Kernel {
	return &Sum{KernelOperator: k.KernelOperator.CloneWithTheta(theta).(KernelOperator)}
}
func (k *Sum) Eval(X, Y mat.Matrix, evalGradient bool) (*mat.Dense, *t.Dense) {
	K1, K1g := k.K1.Eval(X, Y, evalGradient)
	K2, K2g := k.K2.Eval(X, Y, evalGradient)
	K1.Add(K1, K2)
	var Kg *t.Dense
	if evalGradient {
		s1, s2 := K1g.Shape(), K2g.Shape()
		s := t.Shape{s1[0], s1[1], s1[2] + s2[2]}
		Kg = t.NewDense(K1g.Dtype(), s)
		K1gdata, K2gdata, Kgdata := K1g.Data().([]float64), K2g.Data().([]float64), Kg.Data().([]float64)
		for i := range Kgdata {
			i2 := i % s[2]
			i1 := ((i - i2) / s[2]) % s[1]
			i0 := ((i - i2) / s[2]) / s[1]
			if i2 < s1[2] {
				Kgdata[i] = K1gdata[i0*s1[1]*s1[2]+i1*s1[2]+i2]
			} else {
				Kgdata[i] = K2gdata[i0*s2[1]*s2[2]+i1*s2[2]+(i2-s1[2])]
			}
		}
	}
	return K1, Kg
}