// Copyright ©2016 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 distmv

import (
	"math"
	"sort"

	"golang.org/x/exp/rand"
	"golang.org/x/tools/container/intsets"

	"gonum.org/v1/gonum/floats"
	"gonum.org/v1/gonum/mat"
	"gonum.org/v1/gonum/stat"
	"gonum.org/v1/gonum/stat/distuv"
)

// StudentsT is a multivariate Student's T distribution. It is a distribution over
// ℝ^n with the probability density
//  p(y) = (Γ((ν+n)/2) / Γ(ν/2)) * (νπ)^(-n/2) * |Ʃ|^(-1/2) *
//             (1 + 1/ν * (y-μ)ᵀ * Ʃ^-1 * (y-μ))^(-(ν+n)/2)
// where ν is a scalar greater than 2, μ is a vector in ℝ^n, and Ʃ is an n×n
// symmetric positive definite matrix.
//
// In this distribution, ν sets the spread of the distribution, similar to
// the degrees of freedom in a univariate Student's T distribution. As ν → ∞,
// the distribution approaches a multi-variate normal distribution.
// μ is the mean of the distribution, and the covariance is  ν/(ν-2)*Ʃ.
//
// See https://en.wikipedia.org/wiki/Student%27s_t-distribution and
// http://users.isy.liu.se/en/rt/roth/student.pdf for more information.
type StudentsT struct {
	nu float64
	mu []float64
	// If src is altered, rnd must be updated.
	src rand.Source
	rnd *rand.Rand

	sigma mat.SymDense // only stored if needed

	chol       mat.Cholesky
	lower      mat.TriDense
	logSqrtDet float64
	dim        int
}

// NewStudentsT creates a new StudentsT with the given nu, mu, and sigma
// parameters.
//
// NewStudentsT panics if len(mu) == 0, or if len(mu) != sigma.Symmetric(). If
// the covariance matrix is not positive-definite, nil is returned and ok is false.
func NewStudentsT(mu []float64, sigma mat.Symmetric, nu float64, src rand.Source) (dist *StudentsT, ok bool) {
	if len(mu) == 0 {
		panic(badZeroDimension)
	}
	dim := sigma.Symmetric()
	if dim != len(mu) {
		panic(badSizeMismatch)
	}

	s := &StudentsT{
		nu:  nu,
		mu:  make([]float64, dim),
		dim: dim,
		src: src,
	}
	if src != nil {
		s.rnd = rand.New(src)
	}
	copy(s.mu, mu)

	ok = s.chol.Factorize(sigma)
	if !ok {
		return nil, false
	}
	s.sigma = *mat.NewSymDense(dim, nil)
	s.sigma.CopySym(sigma)
	s.chol.LTo(&s.lower)
	s.logSqrtDet = 0.5 * s.chol.LogDet()
	return s, true
}

// ConditionStudentsT returns the Student's T distribution that is the receiver
// conditioned on the input evidence, and the success of the operation.
// The returned Student's T has dimension
// n - len(observed), where n is the dimension of the original receiver.
// The dimension order is preserved during conditioning, so if the value
// of dimension 1 is observed, the returned normal represents dimensions {0, 2, ...}
// of the original Student's T distribution.
//
// ok indicates whether there was a failure during the update. If ok is false
// the operation failed and dist is not usable.
// Mathematically this is impossible, but can occur with finite precision arithmetic.
func (s *StudentsT) ConditionStudentsT(observed []int, values []float64, src rand.Source) (dist *StudentsT, ok bool) {
	if len(observed) == 0 {
		panic("studentst: no observed value")
	}
	if len(observed) != len(values) {
		panic(badInputLength)
	}

	for _, v := range observed {
		if v < 0 || v >= s.dim {
			panic("studentst: observed value out of bounds")
		}
	}

	newNu, newMean, newSigma := studentsTConditional(observed, values, s.nu, s.mu, &s.sigma)
	if newMean == nil {
		return nil, false
	}

	return NewStudentsT(newMean, newSigma, newNu, src)

}

// studentsTConditional updates a Student's T distribution based on the observed samples
// (see documentation for the public function). The Gaussian conditional update
// is treated as a special case when  nu == math.Inf(1).
func studentsTConditional(observed []int, values []float64, nu float64, mu []float64, sigma mat.Symmetric) (newNu float64, newMean []float64, newSigma *mat.SymDense) {
	dim := len(mu)
	ob := len(observed)

	unobserved := findUnob(observed, dim)

	unob := len(unobserved)
	if unob == 0 {
		panic("stat: all dimensions observed")
	}

	mu1 := make([]float64, unob)
	for i, v := range unobserved {
		mu1[i] = mu[v]
	}
	mu2 := make([]float64, ob) // really v - mu2
	for i, v := range observed {
		mu2[i] = values[i] - mu[v]
	}

	var sigma11, sigma22 mat.SymDense
	sigma11.SubsetSym(sigma, unobserved)
	sigma22.SubsetSym(sigma, observed)

	sigma21 := mat.NewDense(ob, unob, nil)
	for i, r := range observed {
		for j, c := range unobserved {
			v := sigma.At(r, c)
			sigma21.Set(i, j, v)
		}
	}

	var chol mat.Cholesky
	ok := chol.Factorize(&sigma22)
	if !ok {
		return math.NaN(), nil, nil
	}

	// Compute mu_1 + sigma_{2,1}ᵀ * sigma_{2,2}^-1 (v - mu_2).
	v := mat.NewVecDense(ob, mu2)
	var tmp, tmp2 mat.VecDense
	err := chol.SolveVecTo(&tmp, v)
	if err != nil {
		return math.NaN(), nil, nil
	}
	tmp2.MulVec(sigma21.T(), &tmp)

	for i := range mu1 {
		mu1[i] += tmp2.At(i, 0)
	}

	// Compute tmp4 = sigma_{2,1}ᵀ * sigma_{2,2}^-1 * sigma_{2,1}.
	// TODO(btracey): Should this be a method of SymDense?
	var tmp3, tmp4 mat.Dense
	err = chol.SolveTo(&tmp3, sigma21)
	if err != nil {
		return math.NaN(), nil, nil
	}
	tmp4.Mul(sigma21.T(), &tmp3)

	// Compute sigma_{1,1} - tmp4
	// TODO(btracey): If tmp4 can constructed with a method, then this can be
	// replaced with SubSym.
	for i := 0; i < len(unobserved); i++ {
		for j := i; j < len(unobserved); j++ {
			v := sigma11.At(i, j)
			sigma11.SetSym(i, j, v-tmp4.At(i, j))
		}
	}

	// The computed variables are accurate for a Normal.
	if math.IsInf(nu, 1) {
		return nu, mu1, &sigma11
	}

	// Compute beta = (v - mu_2)ᵀ * sigma_{2,2}^-1 * (v - mu_2)ᵀ
	beta := mat.Dot(v, &tmp)

	// Scale the covariance matrix
	sigma11.ScaleSym((nu+beta)/(nu+float64(ob)), &sigma11)

	return nu + float64(ob), mu1, &sigma11
}

// findUnob returns the unobserved variables (the complementary set to observed).
// findUnob panics if any value repeated in observed.
func findUnob(observed []int, dim int) (unobserved []int) {
	var setOb intsets.Sparse
	for _, v := range observed {
		setOb.Insert(v)
	}
	var setAll intsets.Sparse
	for i := 0; i < dim; i++ {
		setAll.Insert(i)
	}
	var setUnob intsets.Sparse
	setUnob.Difference(&setAll, &setOb)
	unobserved = setUnob.AppendTo(nil)
	sort.Ints(unobserved)
	return unobserved
}

// CovarianceMatrix calculates the covariance matrix of the distribution,
// storing the result in dst. Upon return, the value at element {i, j} of the
// covariance matrix is equal to the covariance of the i^th and j^th variables.
//  covariance(i, j) = E[(x_i - E[x_i])(x_j - E[x_j])]
// If the dst matrix is empty it will be resized to the correct dimensions,
// otherwise dst must match the dimension of the receiver or CovarianceMatrix
// will panic.
func (st *StudentsT) CovarianceMatrix(dst *mat.SymDense) {
	if dst.IsEmpty() {
		*dst = *(dst.GrowSym(st.dim).(*mat.SymDense))
	} else if dst.Symmetric() != st.dim {
		panic("studentst: input matrix size mismatch")
	}
	dst.CopySym(&st.sigma)
	dst.ScaleSym(st.nu/(st.nu-2), dst)
}

// Dim returns the dimension of the distribution.
func (s *StudentsT) Dim() int {
	return s.dim
}

// LogProb computes the log of the pdf of the point x.
func (s *StudentsT) LogProb(y []float64) float64 {
	if len(y) != s.dim {
		panic(badInputLength)
	}

	nu := s.nu
	n := float64(s.dim)
	lg1, _ := math.Lgamma((nu + n) / 2)
	lg2, _ := math.Lgamma(nu / 2)

	t1 := lg1 - lg2 - n/2*math.Log(nu*math.Pi) - s.logSqrtDet

	mahal := stat.Mahalanobis(mat.NewVecDense(len(y), y), mat.NewVecDense(len(s.mu), s.mu), &s.chol)
	mahal *= mahal
	return t1 - ((nu+n)/2)*math.Log(1+mahal/nu)
}

// MarginalStudentsT returns the marginal distribution of the given input variables,
// and the success of the operation.
// That is, MarginalStudentsT returns
//  p(x_i) = \int_{x_o} p(x_i | x_o) p(x_o) dx_o
// where x_i are the dimensions in the input, and x_o are the remaining dimensions.
// See https://en.wikipedia.org/wiki/Marginal_distribution for more information.
//
// The input src is passed to the created StudentsT.
//
// ok indicates whether there was a failure during the marginalization. If ok is false
// the operation failed and dist is not usable.
// Mathematically this is impossible, but can occur with finite precision arithmetic.
func (s *StudentsT) MarginalStudentsT(vars []int, src rand.Source) (dist *StudentsT, ok bool) {
	newMean := make([]float64, len(vars))
	for i, v := range vars {
		newMean[i] = s.mu[v]
	}
	var newSigma mat.SymDense
	newSigma.SubsetSym(&s.sigma, vars)
	return NewStudentsT(newMean, &newSigma, s.nu, src)
}

// MarginalStudentsTSingle returns the marginal distribution of the given input variable.
// That is, MarginalStudentsTSingle returns
//  p(x_i) = \int_{x_o} p(x_i | x_o) p(x_o) dx_o
// where i is the input index, and x_o are the remaining dimensions.
// See https://en.wikipedia.org/wiki/Marginal_distribution for more information.
//
// The input src is passed to the call to NewStudentsT.
func (s *StudentsT) MarginalStudentsTSingle(i int, src rand.Source) distuv.StudentsT {
	return distuv.StudentsT{
		Mu:    s.mu[i],
		Sigma: math.Sqrt(s.sigma.At(i, i)),
		Nu:    s.nu,
		Src:   src,
	}
}

// TODO(btracey): Implement marginal single. Need to modify univariate StudentsT
// to be three-parameter.

// Mean returns the mean of the probability distribution at x. If the
// input argument is nil, a new slice will be allocated, otherwise the result
// will be put in-place into the receiver.
func (s *StudentsT) Mean(x []float64) []float64 {
	x = reuseAs(x, s.dim)
	copy(x, s.mu)
	return x
}

// Nu returns the degrees of freedom parameter of the distribution.
func (s *StudentsT) Nu() float64 {
	return s.nu
}

// Prob computes the value of the probability density function at x.
func (s *StudentsT) Prob(y []float64) float64 {
	return math.Exp(s.LogProb(y))
}

// Rand generates a random number according to the distributon.
// If the input slice is nil, new memory is allocated, otherwise the result is stored
// in place.
func (s *StudentsT) Rand(x []float64) []float64 {
	// If Y is distributed according to N(0,Sigma), and U is chi^2 with
	// parameter ν, then
	//  X = mu + Y * sqrt(nu / U)
	// X is distributed according to this distribution.

	// Generate Y.
	x = reuseAs(x, s.dim)
	tmp := make([]float64, s.dim)
	if s.rnd == nil {
		for i := range x {
			tmp[i] = rand.NormFloat64()
		}
	} else {
		for i := range x {
			tmp[i] = s.rnd.NormFloat64()
		}
	}
	xVec := mat.NewVecDense(s.dim, x)
	tmpVec := mat.NewVecDense(s.dim, tmp)
	xVec.MulVec(&s.lower, tmpVec)

	u := distuv.ChiSquared{K: s.nu, Src: s.src}.Rand()
	floats.Scale(math.Sqrt(s.nu/u), x)

	floats.Add(x, s.mu)
	return x
}