blob: 91a6694876e5bc7473304b35314a3173d0bdbfd0 [file] [log] [blame]
// Copyright ©2013 The Gonum Authors. All rights reserved.
// Use of this code is governed by a BSD-style
// license that can be found in the LICENSE file.
package cmplxs
import (
"errors"
"math"
"math/cmplx"
"gonum.org/v1/gonum/cmplxs/cscalar"
"gonum.org/v1/gonum/internal/asm/c128"
)
const (
zeroLength = "cmplxs: zero length slice"
shortSpan = "cmplxs: slice length less than 2"
badLength = "cmplxs: slice lengths do not match"
badDstLength = "cmplxs: destination slice length does not match input"
)
// Abs calculates the absolute values of the elements of s, and stores them in dst.
// It panics if the argument lengths do not match.
func Abs(dst []float64, s []complex128) {
if len(dst) != len(s) {
panic(badDstLength)
}
for i, v := range s {
dst[i] = cmplx.Abs(v)
}
}
// Add adds, element-wise, the elements of s and dst, and stores the result in dst.
// It panics if the argument lengths do not match.
func Add(dst, s []complex128) {
if len(dst) != len(s) {
panic(badLength)
}
c128.AxpyUnitaryTo(dst, 1, s, dst)
}
// AddTo adds, element-wise, the elements of s and t and
// stores the result in dst.
// It panics if the argument lengths do not match.
func AddTo(dst, s, t []complex128) []complex128 {
if len(s) != len(t) {
panic(badLength)
}
if len(dst) != len(s) {
panic(badDstLength)
}
c128.AxpyUnitaryTo(dst, 1, s, t)
return dst
}
// AddConst adds the scalar c to all of the values in dst.
func AddConst(c complex128, dst []complex128) {
c128.AddConst(c, dst)
}
// AddScaled performs dst = dst + alpha * s.
// It panics if the slice argument lengths do not match.
func AddScaled(dst []complex128, alpha complex128, s []complex128) {
if len(dst) != len(s) {
panic(badLength)
}
c128.AxpyUnitaryTo(dst, alpha, s, dst)
}
// AddScaledTo performs dst = y + alpha * s, where alpha is a scalar,
// and dst, y and s are all slices.
// It panics if the slice argument lengths do not match.
//
// At the return of the function, dst[i] = y[i] + alpha * s[i]
func AddScaledTo(dst, y []complex128, alpha complex128, s []complex128) []complex128 {
if len(s) != len(y) {
panic(badLength)
}
if len(dst) != len(y) {
panic(badDstLength)
}
c128.AxpyUnitaryTo(dst, alpha, s, y)
return dst
}
// Count applies the function f to every element of s and returns the number
// of times the function returned true.
func Count(f func(complex128) bool, s []complex128) int {
var n int
for _, val := range s {
if f(val) {
n++
}
}
return n
}
// Complex fills each of the elements of dst with the complex number
// constructed from the corresponding elements of real and imag.
// It panics if the argument lengths do not match.
func Complex(dst []complex128, real, imag []float64) []complex128 {
if len(real) != len(imag) {
panic(badLength)
}
if len(dst) != len(real) {
panic(badDstLength)
}
if len(dst) == 0 {
return dst
}
for i, r := range real {
dst[i] = complex(r, imag[i])
}
return dst
}
// CumProd finds the cumulative product of elements of s and store it in
// place into dst so that
// dst[i] = s[i] * s[i-1] * s[i-2] * ... * s[0]
// It panics if the argument lengths do not match.
func CumProd(dst, s []complex128) []complex128 {
if len(dst) != len(s) {
panic(badDstLength)
}
if len(dst) == 0 {
return dst
}
return c128.CumProd(dst, s)
}
// CumSum finds the cumulative sum of elements of s and stores it in place
// into dst so that
// dst[i] = s[i] + s[i-1] + s[i-2] + ... + s[0]
// It panics if the argument lengths do not match.
func CumSum(dst, s []complex128) []complex128 {
if len(dst) != len(s) {
panic(badDstLength)
}
if len(dst) == 0 {
return dst
}
return c128.CumSum(dst, s)
}
// Distance computes the L-norm of s - t. See Norm for special cases.
// It panics if the slice argument lengths do not match.
func Distance(s, t []complex128, L float64) float64 {
if len(s) != len(t) {
panic(badLength)
}
if len(s) == 0 {
return 0
}
var norm float64
switch {
case L == 2:
return c128.L2DistanceUnitary(s, t)
case L == 1:
for i, v := range s {
norm += cmplx.Abs(t[i] - v)
}
return norm
case math.IsInf(L, 1):
for i, v := range s {
absDiff := cmplx.Abs(t[i] - v)
if absDiff > norm {
norm = absDiff
}
}
return norm
default:
for i, v := range s {
norm += math.Pow(cmplx.Abs(t[i]-v), L)
}
return math.Pow(norm, 1/L)
}
}
// Div performs element-wise division dst / s
// and stores the result in dst.
// It panics if the argument lengths do not match.
func Div(dst, s []complex128) {
if len(dst) != len(s) {
panic(badLength)
}
c128.Div(dst, s)
}
// DivTo performs element-wise division s / t
// and stores the result in dst.
// It panics if the argument lengths do not match.
func DivTo(dst, s, t []complex128) []complex128 {
if len(s) != len(t) {
panic(badLength)
}
if len(dst) != len(s) {
panic(badDstLength)
}
return c128.DivTo(dst, s, t)
}
// Dot computes the dot product of s1 and s2, i.e.
// sum_{i = 1}^N conj(s1[i])*s2[i].
// It panics if the argument lengths do not match.
func Dot(s1, s2 []complex128) complex128 {
if len(s1) != len(s2) {
panic(badLength)
}
return c128.DotUnitary(s1, s2)
}
// Equal returns true when the slices have equal lengths and
// all elements are numerically identical.
func Equal(s1, s2 []complex128) bool {
if len(s1) != len(s2) {
return false
}
for i, val := range s1 {
if s2[i] != val {
return false
}
}
return true
}
// EqualApprox returns true when the slices have equal lengths and
// all element pairs have an absolute tolerance less than tol or a
// relative tolerance less than tol.
func EqualApprox(s1, s2 []complex128, tol float64) bool {
if len(s1) != len(s2) {
return false
}
for i, a := range s1 {
if !cscalar.EqualWithinAbsOrRel(a, s2[i], tol, tol) {
return false
}
}
return true
}
// EqualFunc returns true when the slices have the same lengths
// and the function returns true for all element pairs.
func EqualFunc(s1, s2 []complex128, f func(complex128, complex128) bool) bool {
if len(s1) != len(s2) {
return false
}
for i, val := range s1 {
if !f(val, s2[i]) {
return false
}
}
return true
}
// EqualLengths returns true when all of the slices have equal length,
// and false otherwise. It also eturns true when there are no input slices.
func EqualLengths(slices ...[]complex128) bool {
// This length check is needed: http://play.golang.org/p/sdty6YiLhM
if len(slices) == 0 {
return true
}
l := len(slices[0])
for i := 1; i < len(slices); i++ {
if len(slices[i]) != l {
return false
}
}
return true
}
// Find applies f to every element of s and returns the indices of the first
// k elements for which the f returns true, or all such elements
// if k < 0.
// Find will reslice inds to have 0 length, and will append
// found indices to inds.
// If k > 0 and there are fewer than k elements in s satisfying f,
// all of the found elements will be returned along with an error.
// At the return of the function, the input inds will be in an undetermined state.
func Find(inds []int, f func(complex128) bool, s []complex128, k int) ([]int, error) {
// inds is also returned to allow for calling with nil.
// Reslice inds to have zero length.
inds = inds[:0]
// If zero elements requested, can just return.
if k == 0 {
return inds, nil
}
// If k < 0, return all of the found indices.
if k < 0 {
for i, val := range s {
if f(val) {
inds = append(inds, i)
}
}
return inds, nil
}
// Otherwise, find the first k elements.
nFound := 0
for i, val := range s {
if f(val) {
inds = append(inds, i)
nFound++
if nFound == k {
return inds, nil
}
}
}
// Finished iterating over the loop, which means k elements were not found.
return inds, errors.New("cmplxs: insufficient elements found")
}
// HasNaN returns true when the slice s has any values that are NaN and false
// otherwise.
func HasNaN(s []complex128) bool {
for _, v := range s {
if cmplx.IsNaN(v) {
return true
}
}
return false
}
// Imag places the imaginary components of src into dst.
// It panics if the argument lengths do not match.
func Imag(dst []float64, src []complex128) []float64 {
if len(dst) != len(src) {
panic(badDstLength)
}
if len(dst) == 0 {
return dst
}
for i, z := range src {
dst[i] = imag(z)
}
return dst
}
// LogSpan returns a set of n equally spaced points in log space between,
// l and u where N is equal to len(dst). The first element of the
// resulting dst will be l and the final element of dst will be u.
// Panics if len(dst) < 2
// Note that this call will return NaNs if either l or u are negative, and
// will return all zeros if l or u is zero.
// Also returns the mutated slice dst, so that it can be used in range, like:
//
// for i, x := range LogSpan(dst, l, u) { ... }
func LogSpan(dst []complex128, l, u complex128) []complex128 {
Span(dst, cmplx.Log(l), cmplx.Log(u))
for i := range dst {
dst[i] = cmplx.Exp(dst[i])
}
return dst
}
// MaxAbs returns the maximum absolute value in the input slice.
// It panics if s is zero length.
func MaxAbs(s []complex128) complex128 {
return s[MaxAbsIdx(s)]
}
// MaxAbsIdx returns the index of the maximum absolute value in the input slice.
// If several entries have the maximum absolute value, the first such index is
// returned.
// It panics if s is zero length.
func MaxAbsIdx(s []complex128) int {
if len(s) == 0 {
panic(zeroLength)
}
max := math.NaN()
var ind int
for i, v := range s {
if cmplx.IsNaN(v) {
continue
}
if a := cmplx.Abs(v); a > max || math.IsNaN(max) {
max = a
ind = i
}
}
return ind
}
// MinAbs returns the minimum absolute value in the input slice.
// It panics if s is zero length.
func MinAbs(s []complex128) complex128 {
return s[MinAbsIdx(s)]
}
// MinAbsIdx returns the index of the minimum absolute value in the input slice. If several
// entries have the minimum absolute value, the first such index is returned.
// It panics if s is zero length.
func MinAbsIdx(s []complex128) int {
if len(s) == 0 {
panic(zeroLength)
}
min := math.NaN()
var ind int
for i, v := range s {
if cmplx.IsNaN(v) {
continue
}
if a := cmplx.Abs(v); a < min || math.IsNaN(min) {
min = a
ind = i
}
}
return ind
}
// Mul performs element-wise multiplication between dst
// and s and stores the result in dst.
// It panics if the argument lengths do not match.
func Mul(dst, s []complex128) {
if len(dst) != len(s) {
panic(badLength)
}
for i, val := range s {
dst[i] *= val
}
}
// MulTo performs element-wise multiplication between s
// and t and stores the result in dst.
// It panics if the argument lengths do not match.
func MulTo(dst, s, t []complex128) []complex128 {
if len(s) != len(t) {
panic(badLength)
}
if len(dst) != len(s) {
panic(badDstLength)
}
for i, val := range t {
dst[i] = val * s[i]
}
return dst
}
// NearestIdx returns the index of the element in s
// whose value is nearest to v. If several such
// elements exist, the lowest index is returned.
// It panics if s is zero length.
func NearestIdx(s []complex128, v complex128) int {
if len(s) == 0 {
panic(zeroLength)
}
switch {
case cmplx.IsNaN(v):
return 0
case cmplx.IsInf(v):
return MaxAbsIdx(s)
}
var ind int
dist := math.NaN()
for i, val := range s {
newDist := cmplx.Abs(v - val)
// A NaN distance will not be closer.
if math.IsNaN(newDist) {
continue
}
if newDist < dist || math.IsNaN(dist) {
dist = newDist
ind = i
}
}
return ind
}
// Norm returns the L-norm of the slice S, defined as
// (sum_{i=1}^N abs(s[i])^L)^{1/L}
// Special cases:
// L = math.Inf(1) gives the maximum absolute value.
// Does not correctly compute the zero norm (use Count).
func Norm(s []complex128, L float64) float64 {
// Should this complain if L is not positive?
// Should this be done in log space for better numerical stability?
// would be more cost
// maybe only if L is high?
if len(s) == 0 {
return 0
}
var norm float64
switch {
case L == 2:
return c128.L2NormUnitary(s)
case L == 1:
for _, v := range s {
norm += cmplx.Abs(v)
}
return norm
case math.IsInf(L, 1):
for _, v := range s {
norm = math.Max(norm, cmplx.Abs(v))
}
return norm
default:
for _, v := range s {
norm += math.Pow(cmplx.Abs(v), L)
}
return math.Pow(norm, 1/L)
}
}
// Prod returns the product of the elements of the slice.
// Returns 1 if len(s) = 0.
func Prod(s []complex128) complex128 {
prod := 1 + 0i
for _, val := range s {
prod *= val
}
return prod
}
// Real places the real components of src into dst.
// It panics if the argument lengths do not match.
func Real(dst []float64, src []complex128) []float64 {
if len(dst) != len(src) {
panic(badDstLength)
}
if len(dst) == 0 {
return dst
}
for i, z := range src {
dst[i] = real(z)
}
return dst
}
// Reverse reverses the order of elements in the slice.
func Reverse(s []complex128) {
for i, j := 0, len(s)-1; i < j; i, j = i+1, j-1 {
s[i], s[j] = s[j], s[i]
}
}
// Same returns true when the input slices have the same length and all
// elements have the same value with NaN treated as the same.
func Same(s, t []complex128) bool {
if len(s) != len(t) {
return false
}
for i, v := range s {
w := t[i]
if v != w && !(cmplx.IsNaN(v) && cmplx.IsNaN(w)) {
return false
}
}
return true
}
// Scale multiplies every element in dst by the scalar c.
func Scale(c complex128, dst []complex128) {
if len(dst) > 0 {
c128.ScalUnitary(c, dst)
}
}
// ScaleTo multiplies the elements in s by c and stores the result in dst.
// It panics if the slice argument lengths do not match.
func ScaleTo(dst []complex128, c complex128, s []complex128) []complex128 {
if len(dst) != len(s) {
panic(badDstLength)
}
if len(dst) > 0 {
c128.ScalUnitaryTo(dst, c, s)
}
return dst
}
// Span returns a set of N equally spaced points between l and u, where N
// is equal to the length of the destination. The first element of the destination
// is l, the final element of the destination is u.
// It panics if the length of dst is less than 2.
//
// Span also returns the mutated slice dst, so that it can be used in range expressions,
// like:
//
// for i, x := range Span(dst, l, u) { ... }
func Span(dst []complex128, l, u complex128) []complex128 {
n := len(dst)
if n < 2 {
panic(shortSpan)
}
// Special cases for Inf and NaN.
switch {
case cmplx.IsNaN(l):
for i := range dst[:len(dst)-1] {
dst[i] = cmplx.NaN()
}
dst[len(dst)-1] = u
return dst
case cmplx.IsNaN(u):
for i := range dst[1:] {
dst[i+1] = cmplx.NaN()
}
dst[0] = l
return dst
case cmplx.IsInf(l) && cmplx.IsInf(u):
for i := range dst {
dst[i] = cmplx.Inf()
}
return dst
case cmplx.IsInf(l):
for i := range dst[:len(dst)-1] {
dst[i] = l
}
dst[len(dst)-1] = u
return dst
case cmplx.IsInf(u):
for i := range dst[1:] {
dst[i+1] = u
}
dst[0] = l
return dst
}
step := (u - l) / complex(float64(n-1), 0)
for i := range dst {
dst[i] = l + step*complex(float64(i), 0)
}
return dst
}
// Sub subtracts, element-wise, the elements of s from dst.
// It panics if the argument lengths do not match.
func Sub(dst, s []complex128) {
if len(dst) != len(s) {
panic(badLength)
}
c128.AxpyUnitaryTo(dst, -1, s, dst)
}
// SubTo subtracts, element-wise, the elements of t from s and
// stores the result in dst.
// It panics if the argument lengths do not match.
func SubTo(dst, s, t []complex128) []complex128 {
if len(s) != len(t) {
panic(badLength)
}
if len(dst) != len(s) {
panic(badDstLength)
}
c128.AxpyUnitaryTo(dst, -1, t, s)
return dst
}
// Sum returns the sum of the elements of the slice.
func Sum(s []complex128) complex128 {
return c128.Sum(s)
}