Adds matrix
sub package & adds polynomial regression series (#36)
* updates * updates * tests. * test coverage * fixing test * stride not rows + cols * lu decomp implementation. * poly regression! * poly regression works. * typo.
This commit is contained in:
parent
b3dc3fef3c
commit
a211e88530
10 changed files with 1292 additions and 9 deletions
41
_examples/poly_regression/main.go
Normal file
41
_examples/poly_regression/main.go
Normal file
|
@ -0,0 +1,41 @@
|
||||||
|
package main
|
||||||
|
|
||||||
|
import (
|
||||||
|
"net/http"
|
||||||
|
|
||||||
|
"github.com/wcharczuk/go-chart"
|
||||||
|
)
|
||||||
|
|
||||||
|
func drawChart(res http.ResponseWriter, req *http.Request) {
|
||||||
|
|
||||||
|
/*
|
||||||
|
In this example we add a new type of series, a `PolynomialRegressionSeries` that takes another series as a required argument.
|
||||||
|
InnerSeries only needs to implement `ValueProvider`, so really you could chain `PolynomialRegressionSeries` together if you wanted.
|
||||||
|
*/
|
||||||
|
|
||||||
|
mainSeries := chart.ContinuousSeries{
|
||||||
|
Name: "A test series",
|
||||||
|
XValues: chart.Sequence.Float64(1.0, 100.0), //generates a []float64 from 1.0 to 100.0 in 1.0 step increments, or 100 elements.
|
||||||
|
YValues: chart.Sequence.Random(100, 100), //generates a []float64 randomly from 0 to 100 with 100 elements.
|
||||||
|
}
|
||||||
|
|
||||||
|
polyRegSeries := &chart.PolynomialRegressionSeries{
|
||||||
|
Degree: 3,
|
||||||
|
InnerSeries: mainSeries,
|
||||||
|
}
|
||||||
|
|
||||||
|
graph := chart.Chart{
|
||||||
|
Series: []chart.Series{
|
||||||
|
mainSeries,
|
||||||
|
polyRegSeries,
|
||||||
|
},
|
||||||
|
}
|
||||||
|
|
||||||
|
res.Header().Set("Content-Type", "image/png")
|
||||||
|
graph.Render(chart.PNG, res)
|
||||||
|
}
|
||||||
|
|
||||||
|
func main() {
|
||||||
|
http.HandleFunc("/", drawChart)
|
||||||
|
http.ListenAndServe(":8080", nil)
|
||||||
|
}
|
|
@ -9,7 +9,7 @@ type LinearRegressionSeries struct {
|
||||||
Style Style
|
Style Style
|
||||||
YAxis YAxisType
|
YAxis YAxisType
|
||||||
|
|
||||||
Window int
|
Limit int
|
||||||
Offset int
|
Offset int
|
||||||
InnerSeries ValueProvider
|
InnerSeries ValueProvider
|
||||||
|
|
||||||
|
@ -36,18 +36,18 @@ func (lrs LinearRegressionSeries) GetYAxis() YAxisType {
|
||||||
|
|
||||||
// Len returns the number of elements in the series.
|
// Len returns the number of elements in the series.
|
||||||
func (lrs LinearRegressionSeries) Len() int {
|
func (lrs LinearRegressionSeries) Len() int {
|
||||||
return Math.MinInt(lrs.GetWindow(), lrs.InnerSeries.Len()-lrs.GetOffset())
|
return Math.MinInt(lrs.GetLimit(), lrs.InnerSeries.Len()-lrs.GetOffset())
|
||||||
}
|
}
|
||||||
|
|
||||||
// GetWindow returns the window size.
|
// GetLimit returns the window size.
|
||||||
func (lrs LinearRegressionSeries) GetWindow() int {
|
func (lrs LinearRegressionSeries) GetLimit() int {
|
||||||
if lrs.Window == 0 {
|
if lrs.Limit == 0 {
|
||||||
return lrs.InnerSeries.Len()
|
return lrs.InnerSeries.Len()
|
||||||
}
|
}
|
||||||
return lrs.Window
|
return lrs.Limit
|
||||||
}
|
}
|
||||||
|
|
||||||
// GetEndIndex returns the effective window end.
|
// GetEndIndex returns the effective limit end.
|
||||||
func (lrs LinearRegressionSeries) GetEndIndex() int {
|
func (lrs LinearRegressionSeries) GetEndIndex() int {
|
||||||
offset := lrs.GetOffset() + lrs.Len()
|
offset := lrs.GetOffset() + lrs.Len()
|
||||||
innerSeriesLastIndex := lrs.InnerSeries.Len() - 1
|
innerSeriesLastIndex := lrs.InnerSeries.Len() - 1
|
||||||
|
@ -105,7 +105,6 @@ func (lrs *LinearRegressionSeries) computeCoefficients() {
|
||||||
|
|
||||||
xvalues := NewRingBufferWithCapacity(lrs.Len())
|
xvalues := NewRingBufferWithCapacity(lrs.Len())
|
||||||
for index := startIndex; index < endIndex; index++ {
|
for index := startIndex; index < endIndex; index++ {
|
||||||
|
|
||||||
x, _ := lrs.InnerSeries.GetValue(index)
|
x, _ := lrs.InnerSeries.GetValue(index)
|
||||||
xvalues.Enqueue(x)
|
xvalues.Enqueue(x)
|
||||||
}
|
}
|
||||||
|
|
|
@ -62,7 +62,7 @@ func TestLinearRegressionSeriesWindowAndOffset(t *testing.T) {
|
||||||
linRegSeries := &LinearRegressionSeries{
|
linRegSeries := &LinearRegressionSeries{
|
||||||
InnerSeries: mainSeries,
|
InnerSeries: mainSeries,
|
||||||
Offset: 10,
|
Offset: 10,
|
||||||
Window: 10,
|
Limit: 10,
|
||||||
}
|
}
|
||||||
|
|
||||||
assert.Equal(10, linRegSeries.Len())
|
assert.Equal(10, linRegSeries.Len())
|
||||||
|
|
592
matrix/matrix.go
Normal file
592
matrix/matrix.go
Normal file
|
@ -0,0 +1,592 @@
|
||||||
|
package matrix
|
||||||
|
|
||||||
|
import (
|
||||||
|
"bytes"
|
||||||
|
"errors"
|
||||||
|
"fmt"
|
||||||
|
"math"
|
||||||
|
)
|
||||||
|
|
||||||
|
const (
|
||||||
|
// DefaultEpsilon represents the minimum precision for matrix math operations.
|
||||||
|
DefaultEpsilon = 0.000001
|
||||||
|
)
|
||||||
|
|
||||||
|
var (
|
||||||
|
// ErrDimensionMismatch is a typical error.
|
||||||
|
ErrDimensionMismatch = errors.New("dimension mismatch")
|
||||||
|
|
||||||
|
// ErrSingularValue is a typical error.
|
||||||
|
ErrSingularValue = errors.New("singular value")
|
||||||
|
)
|
||||||
|
|
||||||
|
// New returns a new matrix.
|
||||||
|
func New(rows, cols int, values ...float64) *Matrix {
|
||||||
|
if len(values) == 0 {
|
||||||
|
return &Matrix{
|
||||||
|
stride: cols,
|
||||||
|
epsilon: DefaultEpsilon,
|
||||||
|
elements: make([]float64, rows*cols),
|
||||||
|
}
|
||||||
|
}
|
||||||
|
elems := make([]float64, rows*cols)
|
||||||
|
copy(elems, values)
|
||||||
|
return &Matrix{
|
||||||
|
stride: cols,
|
||||||
|
epsilon: DefaultEpsilon,
|
||||||
|
elements: elems,
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// Identity returns the identity matrix of a given order.
|
||||||
|
func Identity(order int) *Matrix {
|
||||||
|
m := New(order, order)
|
||||||
|
for i := 0; i < order; i++ {
|
||||||
|
m.Set(i, i, 1)
|
||||||
|
}
|
||||||
|
return m
|
||||||
|
}
|
||||||
|
|
||||||
|
// Zero returns a matrix of a given size zeroed.
|
||||||
|
func Zero(rows, cols int) *Matrix {
|
||||||
|
return New(rows, cols)
|
||||||
|
}
|
||||||
|
|
||||||
|
// Ones returns an matrix of ones.
|
||||||
|
func Ones(rows, cols int) *Matrix {
|
||||||
|
ones := make([]float64, rows*cols)
|
||||||
|
for i := 0; i < (rows * cols); i++ {
|
||||||
|
ones[i] = 1
|
||||||
|
}
|
||||||
|
|
||||||
|
return &Matrix{
|
||||||
|
stride: cols,
|
||||||
|
epsilon: DefaultEpsilon,
|
||||||
|
elements: ones,
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// Eye returns the eye matrix.
|
||||||
|
func Eye(n int) *Matrix {
|
||||||
|
m := Zero(n, n)
|
||||||
|
for i := 0; i < len(m.elements); i += n + 1 {
|
||||||
|
m.elements[i] = 1
|
||||||
|
}
|
||||||
|
return m
|
||||||
|
}
|
||||||
|
|
||||||
|
// NewFromArrays creates a matrix from a jagged array set.
|
||||||
|
func NewFromArrays(a [][]float64) *Matrix {
|
||||||
|
rows := len(a)
|
||||||
|
if rows == 0 {
|
||||||
|
return nil
|
||||||
|
}
|
||||||
|
cols := len(a[0])
|
||||||
|
m := New(rows, cols)
|
||||||
|
for row := 0; row < rows; row++ {
|
||||||
|
for col := 0; col < cols; col++ {
|
||||||
|
m.Set(row, col, a[row][col])
|
||||||
|
}
|
||||||
|
}
|
||||||
|
return m
|
||||||
|
}
|
||||||
|
|
||||||
|
// Matrix represents a 2d dense array of floats.
|
||||||
|
type Matrix struct {
|
||||||
|
epsilon float64
|
||||||
|
elements []float64
|
||||||
|
stride int
|
||||||
|
}
|
||||||
|
|
||||||
|
// String returns a string representation of the matrix.
|
||||||
|
func (m *Matrix) String() string {
|
||||||
|
buffer := bytes.NewBuffer(nil)
|
||||||
|
rows, cols := m.Size()
|
||||||
|
|
||||||
|
for row := 0; row < rows; row++ {
|
||||||
|
for col := 0; col < cols; col++ {
|
||||||
|
buffer.WriteString(f64s(m.Get(row, col)))
|
||||||
|
buffer.WriteRune(' ')
|
||||||
|
}
|
||||||
|
buffer.WriteRune('\n')
|
||||||
|
}
|
||||||
|
return buffer.String()
|
||||||
|
}
|
||||||
|
|
||||||
|
// Epsilon returns the maximum precision for math operations.
|
||||||
|
func (m *Matrix) Epsilon() float64 {
|
||||||
|
return m.epsilon
|
||||||
|
}
|
||||||
|
|
||||||
|
// WithEpsilon sets the epsilon on the matrix and returns a reference to the matrix.
|
||||||
|
func (m *Matrix) WithEpsilon(epsilon float64) *Matrix {
|
||||||
|
m.epsilon = epsilon
|
||||||
|
return m
|
||||||
|
}
|
||||||
|
|
||||||
|
// Each applies the action to each element of the matrix in
|
||||||
|
// rows => cols order.
|
||||||
|
func (m *Matrix) Each(action func(row, cow int, value float64)) {
|
||||||
|
rows, cols := m.Size()
|
||||||
|
for row := 0; row < rows; row++ {
|
||||||
|
for col := 0; col < cols; col++ {
|
||||||
|
action(row, col, m.Get(row, col))
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// Round rounds all the values in a matrix to it epsilon,
|
||||||
|
// returning a reference to the original
|
||||||
|
func (m *Matrix) Round() *Matrix {
|
||||||
|
rows, cols := m.Size()
|
||||||
|
for row := 0; row < rows; row++ {
|
||||||
|
for col := 0; col < cols; col++ {
|
||||||
|
m.Set(row, col, roundToEpsilon(m.Get(row, col), m.epsilon))
|
||||||
|
}
|
||||||
|
}
|
||||||
|
return m
|
||||||
|
}
|
||||||
|
|
||||||
|
// Arrays returns the matrix as a two dimensional jagged array.
|
||||||
|
func (m *Matrix) Arrays() [][]float64 {
|
||||||
|
rows, cols := m.Size()
|
||||||
|
a := make([][]float64, rows)
|
||||||
|
|
||||||
|
for row := 0; row < rows; row++ {
|
||||||
|
a[row] = make([]float64, cols)
|
||||||
|
|
||||||
|
for col := 0; col < cols; col++ {
|
||||||
|
a[row][col] = m.Get(row, col)
|
||||||
|
}
|
||||||
|
}
|
||||||
|
return a
|
||||||
|
}
|
||||||
|
|
||||||
|
// Size returns the dimensions of the matrix.
|
||||||
|
func (m *Matrix) Size() (rows, cols int) {
|
||||||
|
rows = len(m.elements) / m.stride
|
||||||
|
cols = m.stride
|
||||||
|
return
|
||||||
|
}
|
||||||
|
|
||||||
|
// IsSquare returns if the row count is equal to the column count.
|
||||||
|
func (m *Matrix) IsSquare() bool {
|
||||||
|
return m.stride == (len(m.elements) / m.stride)
|
||||||
|
}
|
||||||
|
|
||||||
|
// IsSymmetric returns if the matrix is symmetric about its diagonal.
|
||||||
|
func (m *Matrix) IsSymmetric() bool {
|
||||||
|
rows, cols := m.Size()
|
||||||
|
|
||||||
|
if rows != cols {
|
||||||
|
return false
|
||||||
|
}
|
||||||
|
|
||||||
|
for i := 0; i < rows; i++ {
|
||||||
|
for j := 0; j < i; j++ {
|
||||||
|
if m.Get(i, j) != m.Get(j, i) {
|
||||||
|
return false
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
return true
|
||||||
|
}
|
||||||
|
|
||||||
|
// Get returns the element at the given row, col.
|
||||||
|
func (m *Matrix) Get(row, col int) float64 {
|
||||||
|
index := (m.stride * row) + col
|
||||||
|
return m.elements[index]
|
||||||
|
}
|
||||||
|
|
||||||
|
// Set sets a value.
|
||||||
|
func (m *Matrix) Set(row, col int, val float64) {
|
||||||
|
index := (m.stride * row) + col
|
||||||
|
m.elements[index] = val
|
||||||
|
}
|
||||||
|
|
||||||
|
// Col returns a column of the matrix as a vector.
|
||||||
|
func (m *Matrix) Col(col int) Vector {
|
||||||
|
rows, _ := m.Size()
|
||||||
|
values := make([]float64, rows)
|
||||||
|
for row := 0; row < rows; row++ {
|
||||||
|
values[row] = m.Get(row, col)
|
||||||
|
}
|
||||||
|
return Vector(values)
|
||||||
|
}
|
||||||
|
|
||||||
|
// Row returns a row of the matrix as a vector.
|
||||||
|
func (m *Matrix) Row(row int) Vector {
|
||||||
|
_, cols := m.Size()
|
||||||
|
values := make([]float64, cols)
|
||||||
|
for col := 0; col < cols; col++ {
|
||||||
|
values[col] = m.Get(row, col)
|
||||||
|
}
|
||||||
|
return Vector(values)
|
||||||
|
}
|
||||||
|
|
||||||
|
// SubMatrix returns a sub matrix from a given outer matrix.
|
||||||
|
func (m *Matrix) SubMatrix(i, j, rows, cols int) *Matrix {
|
||||||
|
return &Matrix{
|
||||||
|
elements: m.elements[i*m.stride+j : i*m.stride+j+(rows-1)*m.stride+cols],
|
||||||
|
stride: m.stride,
|
||||||
|
epsilon: m.epsilon,
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// ScaleRow applies a scale to an entire row.
|
||||||
|
func (m *Matrix) ScaleRow(row int, scale float64) {
|
||||||
|
startIndex := row * m.stride
|
||||||
|
for i := startIndex; i < m.stride; i++ {
|
||||||
|
m.elements[i] = m.elements[i] * scale
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
func (m *Matrix) scaleAddRow(rd int, rs int, f float64) {
|
||||||
|
indexd := rd * m.stride
|
||||||
|
indexs := rs * m.stride
|
||||||
|
for col := 0; col < m.stride; col++ {
|
||||||
|
m.elements[indexd] += f * m.elements[indexs]
|
||||||
|
indexd++
|
||||||
|
indexs++
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// SwapRows swaps a row in the matrix in place.
|
||||||
|
func (m *Matrix) SwapRows(i, j int) {
|
||||||
|
var vi, vj float64
|
||||||
|
for col := 0; col < m.stride; col++ {
|
||||||
|
vi = m.Get(i, col)
|
||||||
|
vj = m.Get(j, col)
|
||||||
|
m.Set(i, col, vj)
|
||||||
|
m.Set(j, col, vi)
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// Augment concatenates two matrices about the horizontal.
|
||||||
|
func (m *Matrix) Augment(m2 *Matrix) (*Matrix, error) {
|
||||||
|
mr, mc := m.Size()
|
||||||
|
m2r, m2c := m2.Size()
|
||||||
|
if mr != m2r {
|
||||||
|
return nil, ErrDimensionMismatch
|
||||||
|
}
|
||||||
|
|
||||||
|
m3 := Zero(mr, mc+m2c)
|
||||||
|
for row := 0; row < mr; row++ {
|
||||||
|
for col := 0; col < mc; col++ {
|
||||||
|
m3.Set(row, col, m.Get(row, col))
|
||||||
|
}
|
||||||
|
for col := 0; col < m2c; col++ {
|
||||||
|
m3.Set(row, mc+col, m2.Get(row, col))
|
||||||
|
}
|
||||||
|
}
|
||||||
|
return m3, nil
|
||||||
|
}
|
||||||
|
|
||||||
|
// Copy returns a duplicate of a given matrix.
|
||||||
|
func (m *Matrix) Copy() *Matrix {
|
||||||
|
m2 := &Matrix{stride: m.stride, epsilon: m.epsilon, elements: make([]float64, len(m.elements))}
|
||||||
|
copy(m2.elements, m.elements)
|
||||||
|
return m2
|
||||||
|
}
|
||||||
|
|
||||||
|
// DiagonalVector returns a vector from the diagonal of a matrix.
|
||||||
|
func (m *Matrix) DiagonalVector() Vector {
|
||||||
|
rows, cols := m.Size()
|
||||||
|
rank := minInt(rows, cols)
|
||||||
|
values := make([]float64, rank)
|
||||||
|
|
||||||
|
for index := 0; index < rank; index++ {
|
||||||
|
values[index] = m.Get(index, index)
|
||||||
|
}
|
||||||
|
return Vector(values)
|
||||||
|
}
|
||||||
|
|
||||||
|
// Diagonal returns a matrix from the diagonal of a matrix.
|
||||||
|
func (m *Matrix) Diagonal() *Matrix {
|
||||||
|
rows, cols := m.Size()
|
||||||
|
rank := minInt(rows, cols)
|
||||||
|
m2 := New(rank, rank)
|
||||||
|
|
||||||
|
for index := 0; index < rank; index++ {
|
||||||
|
m2.Set(index, index, m.Get(index, index))
|
||||||
|
}
|
||||||
|
return m2
|
||||||
|
}
|
||||||
|
|
||||||
|
// Equals returns if a matrix equals another matrix.
|
||||||
|
func (m *Matrix) Equals(other *Matrix) bool {
|
||||||
|
if other == nil && m != nil {
|
||||||
|
return false
|
||||||
|
} else if other == nil {
|
||||||
|
return true
|
||||||
|
}
|
||||||
|
|
||||||
|
if m.stride != other.stride {
|
||||||
|
return false
|
||||||
|
}
|
||||||
|
|
||||||
|
msize := len(m.elements)
|
||||||
|
m2size := len(other.elements)
|
||||||
|
|
||||||
|
if msize != m2size {
|
||||||
|
return false
|
||||||
|
}
|
||||||
|
|
||||||
|
for i := 0; i < msize; i++ {
|
||||||
|
if m.elements[i] != other.elements[i] {
|
||||||
|
return false
|
||||||
|
}
|
||||||
|
}
|
||||||
|
return true
|
||||||
|
}
|
||||||
|
|
||||||
|
// L returns the matrix with zeros below the diagonal.
|
||||||
|
func (m *Matrix) L() *Matrix {
|
||||||
|
rows, cols := m.Size()
|
||||||
|
m2 := New(rows, cols)
|
||||||
|
for row := 0; row < rows; row++ {
|
||||||
|
for col := row; col < cols; col++ {
|
||||||
|
m2.Set(row, col, m.Get(row, col))
|
||||||
|
}
|
||||||
|
}
|
||||||
|
return m2
|
||||||
|
}
|
||||||
|
|
||||||
|
// U returns the matrix with zeros above the diagonal.
|
||||||
|
// Does not include the diagonal.
|
||||||
|
func (m *Matrix) U() *Matrix {
|
||||||
|
rows, cols := m.Size()
|
||||||
|
m2 := New(rows, cols)
|
||||||
|
for row := 0; row < rows; row++ {
|
||||||
|
for col := 0; col < row && col < cols; col++ {
|
||||||
|
m2.Set(row, col, m.Get(row, col))
|
||||||
|
}
|
||||||
|
}
|
||||||
|
return m2
|
||||||
|
}
|
||||||
|
|
||||||
|
// math operations
|
||||||
|
|
||||||
|
// Multiply multiplies two matrices.
|
||||||
|
func (m *Matrix) Multiply(m2 *Matrix) (m3 *Matrix, err error) {
|
||||||
|
if m.stride*m2.stride != len(m2.elements) {
|
||||||
|
return nil, ErrDimensionMismatch
|
||||||
|
}
|
||||||
|
|
||||||
|
m3 = &Matrix{epsilon: m.epsilon, stride: m2.stride, elements: make([]float64, (len(m.elements)/m.stride)*m2.stride)}
|
||||||
|
for m1c0, m3x := 0, 0; m1c0 < len(m.elements); m1c0 += m.stride {
|
||||||
|
for m2r0 := 0; m2r0 < m2.stride; m2r0++ {
|
||||||
|
for m1x, m2x := m1c0, m2r0; m2x < len(m2.elements); m2x += m2.stride {
|
||||||
|
m3.elements[m3x] += m.elements[m1x] * m2.elements[m2x]
|
||||||
|
m1x++
|
||||||
|
}
|
||||||
|
m3x++
|
||||||
|
}
|
||||||
|
}
|
||||||
|
return
|
||||||
|
}
|
||||||
|
|
||||||
|
// Pivotize does something i'm not sure what.
|
||||||
|
func (m *Matrix) Pivotize() *Matrix {
|
||||||
|
pv := make([]int, m.stride)
|
||||||
|
|
||||||
|
for i := range pv {
|
||||||
|
pv[i] = i
|
||||||
|
}
|
||||||
|
|
||||||
|
for j, dx := 0, 0; j < m.stride; j++ {
|
||||||
|
row := j
|
||||||
|
max := m.elements[dx]
|
||||||
|
for i, ixcj := j, dx; i < m.stride; i++ {
|
||||||
|
if m.elements[ixcj] > max {
|
||||||
|
max = m.elements[ixcj]
|
||||||
|
row = i
|
||||||
|
}
|
||||||
|
ixcj += m.stride
|
||||||
|
}
|
||||||
|
if j != row {
|
||||||
|
pv[row], pv[j] = pv[j], pv[row]
|
||||||
|
}
|
||||||
|
dx += m.stride + 1
|
||||||
|
}
|
||||||
|
p := Zero(m.stride, m.stride)
|
||||||
|
for r, c := range pv {
|
||||||
|
p.elements[r*m.stride+c] = 1
|
||||||
|
}
|
||||||
|
return p
|
||||||
|
}
|
||||||
|
|
||||||
|
// Times returns the product of a matrix and another.
|
||||||
|
func (m *Matrix) Times(m2 *Matrix) (*Matrix, error) {
|
||||||
|
mr, mc := m.Size()
|
||||||
|
m2r, m2c := m2.Size()
|
||||||
|
|
||||||
|
if mc != m2r {
|
||||||
|
return nil, fmt.Errorf("cannot multiply (%dx%d) and (%dx%d)", mr, mc, m2r, m2c)
|
||||||
|
//return nil, ErrDimensionMismatch
|
||||||
|
}
|
||||||
|
|
||||||
|
c := Zero(mr, m2c)
|
||||||
|
|
||||||
|
for i := 0; i < mr; i++ {
|
||||||
|
sums := c.elements[i*c.stride : (i+1)*c.stride]
|
||||||
|
for k, a := range m.elements[i*m.stride : i*m.stride+m.stride] {
|
||||||
|
for j, b := range m2.elements[k*m2.stride : k*m2.stride+m2.stride] {
|
||||||
|
sums[j] += a * b
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
return c, nil
|
||||||
|
}
|
||||||
|
|
||||||
|
// Decompositions
|
||||||
|
|
||||||
|
// LU performs the LU decomposition.
|
||||||
|
func (m *Matrix) LU() (l, u, p *Matrix) {
|
||||||
|
l = Zero(m.stride, m.stride)
|
||||||
|
u = Zero(m.stride, m.stride)
|
||||||
|
p = m.Pivotize()
|
||||||
|
m, _ = p.Multiply(m)
|
||||||
|
for j, jxc0 := 0, 0; j < m.stride; j++ {
|
||||||
|
l.elements[jxc0+j] = 1
|
||||||
|
for i, ixc0 := 0, 0; ixc0 <= jxc0; i++ {
|
||||||
|
sum := 0.
|
||||||
|
for k, kxcj := 0, j; k < i; k++ {
|
||||||
|
sum += u.elements[kxcj] * l.elements[ixc0+k]
|
||||||
|
kxcj += m.stride
|
||||||
|
}
|
||||||
|
u.elements[ixc0+j] = m.elements[ixc0+j] - sum
|
||||||
|
ixc0 += m.stride
|
||||||
|
}
|
||||||
|
for ixc0 := jxc0; ixc0 < len(m.elements); ixc0 += m.stride {
|
||||||
|
sum := 0.
|
||||||
|
for k, kxcj := 0, j; k < j; k++ {
|
||||||
|
sum += u.elements[kxcj] * l.elements[ixc0+k]
|
||||||
|
kxcj += m.stride
|
||||||
|
}
|
||||||
|
l.elements[ixc0+j] = (m.elements[ixc0+j] - sum) / u.elements[jxc0+j]
|
||||||
|
}
|
||||||
|
jxc0 += m.stride
|
||||||
|
}
|
||||||
|
return
|
||||||
|
}
|
||||||
|
|
||||||
|
// QR performs the qr decomposition.
|
||||||
|
func (m *Matrix) QR() (q, r *Matrix) {
|
||||||
|
defer func() {
|
||||||
|
q = q.Round()
|
||||||
|
r = r.Round()
|
||||||
|
}()
|
||||||
|
|
||||||
|
rows, cols := m.Size()
|
||||||
|
qr := m.Copy()
|
||||||
|
q = New(rows, cols)
|
||||||
|
r = New(rows, cols)
|
||||||
|
|
||||||
|
var i, j, k int
|
||||||
|
var norm, s float64
|
||||||
|
|
||||||
|
for k = 0; k < cols; k++ {
|
||||||
|
norm = 0
|
||||||
|
for i = k; i < rows; i++ {
|
||||||
|
norm = math.Hypot(norm, qr.Get(i, k))
|
||||||
|
}
|
||||||
|
|
||||||
|
if norm != 0 {
|
||||||
|
if qr.Get(k, k) < 0 {
|
||||||
|
norm = -norm
|
||||||
|
}
|
||||||
|
|
||||||
|
for i = k; i < rows; i++ {
|
||||||
|
qr.Set(i, k, qr.Get(i, k)/norm)
|
||||||
|
}
|
||||||
|
qr.Set(k, k, qr.Get(k, k)+1.0)
|
||||||
|
|
||||||
|
for j = k + 1; j < cols; j++ {
|
||||||
|
s = 0
|
||||||
|
for i = k; i < rows; i++ {
|
||||||
|
s += qr.Get(i, k) * qr.Get(i, j)
|
||||||
|
}
|
||||||
|
s = -s / qr.Get(k, k)
|
||||||
|
for i = k; i < rows; i++ {
|
||||||
|
qr.Set(i, j, qr.Get(i, j)+s*qr.Get(i, k))
|
||||||
|
|
||||||
|
if i < j {
|
||||||
|
r.Set(i, j, qr.Get(i, j))
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
r.Set(k, k, -norm)
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
//Q Matrix:
|
||||||
|
i, j, k = 0, 0, 0
|
||||||
|
|
||||||
|
for k = cols - 1; k >= 0; k-- {
|
||||||
|
q.Set(k, k, 1.0)
|
||||||
|
for j = k; j < cols; j++ {
|
||||||
|
if qr.Get(k, k) != 0 {
|
||||||
|
s = 0
|
||||||
|
for i = k; i < rows; i++ {
|
||||||
|
s += qr.Get(i, k) * q.Get(i, j)
|
||||||
|
}
|
||||||
|
s = -s / qr.Get(k, k)
|
||||||
|
for i = k; i < rows; i++ {
|
||||||
|
q.Set(i, j, q.Get(i, j)+s*qr.Get(i, k))
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
return
|
||||||
|
}
|
||||||
|
|
||||||
|
// Transpose flips a matrix about its diagonal, returning a new copy.
|
||||||
|
func (m *Matrix) Transpose() *Matrix {
|
||||||
|
rows, cols := m.Size()
|
||||||
|
m2 := Zero(cols, rows)
|
||||||
|
for i := 0; i < rows; i++ {
|
||||||
|
for j := 0; j < cols; j++ {
|
||||||
|
m2.Set(j, i, m.Get(i, j))
|
||||||
|
}
|
||||||
|
}
|
||||||
|
return m2
|
||||||
|
}
|
||||||
|
|
||||||
|
// Inverse returns a matrix such that M*I==1.
|
||||||
|
func (m *Matrix) Inverse() (*Matrix, error) {
|
||||||
|
if !m.IsSymmetric() {
|
||||||
|
return nil, ErrDimensionMismatch
|
||||||
|
}
|
||||||
|
|
||||||
|
rows, cols := m.Size()
|
||||||
|
|
||||||
|
aug, _ := m.Augment(Eye(rows))
|
||||||
|
for i := 0; i < rows; i++ {
|
||||||
|
j := i
|
||||||
|
for k := i; k < rows; k++ {
|
||||||
|
if math.Abs(aug.Get(k, i)) > math.Abs(aug.Get(j, i)) {
|
||||||
|
j = k
|
||||||
|
}
|
||||||
|
}
|
||||||
|
if j != i {
|
||||||
|
aug.SwapRows(i, j)
|
||||||
|
}
|
||||||
|
if aug.Get(i, i) == 0 {
|
||||||
|
return nil, ErrSingularValue
|
||||||
|
}
|
||||||
|
aug.ScaleRow(i, 1.0/aug.Get(i, i))
|
||||||
|
for k := 0; k < rows; k++ {
|
||||||
|
if k == i {
|
||||||
|
continue
|
||||||
|
}
|
||||||
|
aug.scaleAddRow(k, i, -aug.Get(k, i))
|
||||||
|
}
|
||||||
|
}
|
||||||
|
return aug.SubMatrix(0, cols, rows, cols), nil
|
||||||
|
}
|
396
matrix/matrix_test.go
Normal file
396
matrix/matrix_test.go
Normal file
|
@ -0,0 +1,396 @@
|
||||||
|
package matrix
|
||||||
|
|
||||||
|
import (
|
||||||
|
"testing"
|
||||||
|
|
||||||
|
assert "github.com/blendlabs/go-assert"
|
||||||
|
)
|
||||||
|
|
||||||
|
func TestNew(t *testing.T) {
|
||||||
|
assert := assert.New(t)
|
||||||
|
|
||||||
|
m := New(10, 5)
|
||||||
|
rows, cols := m.Size()
|
||||||
|
assert.Equal(10, rows)
|
||||||
|
assert.Equal(5, cols)
|
||||||
|
assert.Zero(m.Get(0, 0))
|
||||||
|
assert.Zero(m.Get(9, 4))
|
||||||
|
}
|
||||||
|
|
||||||
|
func TestNewWithValues(t *testing.T) {
|
||||||
|
assert := assert.New(t)
|
||||||
|
|
||||||
|
m := New(5, 2, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
|
||||||
|
rows, cols := m.Size()
|
||||||
|
assert.Equal(5, rows)
|
||||||
|
assert.Equal(2, cols)
|
||||||
|
assert.Equal(1, m.Get(0, 0))
|
||||||
|
assert.Equal(10, m.Get(4, 1))
|
||||||
|
}
|
||||||
|
|
||||||
|
func TestIdentitiy(t *testing.T) {
|
||||||
|
assert := assert.New(t)
|
||||||
|
|
||||||
|
id := Identity(5)
|
||||||
|
rows, cols := id.Size()
|
||||||
|
assert.Equal(5, rows)
|
||||||
|
assert.Equal(5, cols)
|
||||||
|
assert.Equal(1, id.Get(0, 0))
|
||||||
|
assert.Equal(1, id.Get(1, 1))
|
||||||
|
assert.Equal(1, id.Get(2, 2))
|
||||||
|
assert.Equal(1, id.Get(3, 3))
|
||||||
|
assert.Equal(1, id.Get(4, 4))
|
||||||
|
assert.Equal(0, id.Get(0, 1))
|
||||||
|
assert.Equal(0, id.Get(1, 0))
|
||||||
|
assert.Equal(0, id.Get(4, 0))
|
||||||
|
assert.Equal(0, id.Get(0, 4))
|
||||||
|
}
|
||||||
|
|
||||||
|
func TestNewFromArrays(t *testing.T) {
|
||||||
|
assert := assert.New(t)
|
||||||
|
|
||||||
|
m := NewFromArrays([][]float64{
|
||||||
|
{1, 2, 3, 4},
|
||||||
|
{5, 6, 7, 8},
|
||||||
|
})
|
||||||
|
assert.NotNil(m)
|
||||||
|
|
||||||
|
rows, cols := m.Size()
|
||||||
|
assert.Equal(2, rows)
|
||||||
|
assert.Equal(4, cols)
|
||||||
|
}
|
||||||
|
|
||||||
|
func TestOnes(t *testing.T) {
|
||||||
|
assert := assert.New(t)
|
||||||
|
|
||||||
|
ones := Ones(5, 10)
|
||||||
|
rows, cols := ones.Size()
|
||||||
|
assert.Equal(5, rows)
|
||||||
|
assert.Equal(10, cols)
|
||||||
|
|
||||||
|
for row := 0; row < rows; row++ {
|
||||||
|
for col := 0; col < cols; col++ {
|
||||||
|
assert.Equal(1, ones.Get(row, col))
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
func TestMatrixEpsilon(t *testing.T) {
|
||||||
|
assert := assert.New(t)
|
||||||
|
|
||||||
|
ones := Ones(2, 2)
|
||||||
|
ones = ones.WithEpsilon(0.001)
|
||||||
|
assert.Equal(0.001, ones.Epsilon())
|
||||||
|
}
|
||||||
|
|
||||||
|
func TestMatrixArrays(t *testing.T) {
|
||||||
|
assert := assert.New(t)
|
||||||
|
|
||||||
|
m := NewFromArrays([][]float64{
|
||||||
|
{1, 2, 3},
|
||||||
|
{4, 5, 6},
|
||||||
|
})
|
||||||
|
|
||||||
|
assert.NotNil(m)
|
||||||
|
|
||||||
|
arrays := m.Arrays()
|
||||||
|
|
||||||
|
assert.Equal(arrays, [][]float64{
|
||||||
|
{1, 2, 3},
|
||||||
|
{4, 5, 6},
|
||||||
|
})
|
||||||
|
}
|
||||||
|
|
||||||
|
func TestMatrixIsSquare(t *testing.T) {
|
||||||
|
assert := assert.New(t)
|
||||||
|
|
||||||
|
assert.False(NewFromArrays([][]float64{
|
||||||
|
{1, 2, 3},
|
||||||
|
{4, 5, 6},
|
||||||
|
}).IsSquare())
|
||||||
|
|
||||||
|
assert.False(NewFromArrays([][]float64{
|
||||||
|
{1, 2},
|
||||||
|
{3, 4},
|
||||||
|
{5, 6},
|
||||||
|
}).IsSquare())
|
||||||
|
|
||||||
|
assert.True(NewFromArrays([][]float64{
|
||||||
|
{1, 2},
|
||||||
|
{3, 4},
|
||||||
|
}).IsSquare())
|
||||||
|
}
|
||||||
|
|
||||||
|
func TestMatrixIsSymmetric(t *testing.T) {
|
||||||
|
assert := assert.New(t)
|
||||||
|
|
||||||
|
assert.False(NewFromArrays([][]float64{
|
||||||
|
{1, 2, 3},
|
||||||
|
{2, 1, 2},
|
||||||
|
}).IsSymmetric())
|
||||||
|
|
||||||
|
assert.False(NewFromArrays([][]float64{
|
||||||
|
{1, 2, 3},
|
||||||
|
{4, 5, 6},
|
||||||
|
{7, 8, 9},
|
||||||
|
}).IsSymmetric())
|
||||||
|
|
||||||
|
assert.True(NewFromArrays([][]float64{
|
||||||
|
{1, 2, 3},
|
||||||
|
{2, 1, 2},
|
||||||
|
{3, 2, 1},
|
||||||
|
}).IsSymmetric())
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
func TestMatrixGet(t *testing.T) {
|
||||||
|
assert := assert.New(t)
|
||||||
|
|
||||||
|
m := NewFromArrays([][]float64{
|
||||||
|
{1, 2, 3},
|
||||||
|
{4, 5, 6},
|
||||||
|
{7, 8, 9},
|
||||||
|
})
|
||||||
|
|
||||||
|
assert.Equal(1, m.Get(0, 0))
|
||||||
|
assert.Equal(2, m.Get(0, 1))
|
||||||
|
assert.Equal(3, m.Get(0, 2))
|
||||||
|
assert.Equal(4, m.Get(1, 0))
|
||||||
|
assert.Equal(5, m.Get(1, 1))
|
||||||
|
assert.Equal(6, m.Get(1, 2))
|
||||||
|
assert.Equal(7, m.Get(2, 0))
|
||||||
|
assert.Equal(8, m.Get(2, 1))
|
||||||
|
assert.Equal(9, m.Get(2, 2))
|
||||||
|
}
|
||||||
|
|
||||||
|
func TestMatrixSet(t *testing.T) {
|
||||||
|
assert := assert.New(t)
|
||||||
|
|
||||||
|
m := NewFromArrays([][]float64{
|
||||||
|
{1, 2, 3},
|
||||||
|
{4, 5, 6},
|
||||||
|
{7, 8, 9},
|
||||||
|
})
|
||||||
|
|
||||||
|
m.Set(1, 1, 99)
|
||||||
|
assert.Equal(99, m.Get(1, 1))
|
||||||
|
}
|
||||||
|
|
||||||
|
func TestMatrixCol(t *testing.T) {
|
||||||
|
assert := assert.New(t)
|
||||||
|
|
||||||
|
m := NewFromArrays([][]float64{
|
||||||
|
{1, 2, 3},
|
||||||
|
{4, 5, 6},
|
||||||
|
{7, 8, 9},
|
||||||
|
})
|
||||||
|
|
||||||
|
assert.Equal([]float64{1, 4, 7}, m.Col(0))
|
||||||
|
assert.Equal([]float64{2, 5, 8}, m.Col(1))
|
||||||
|
assert.Equal([]float64{3, 6, 9}, m.Col(2))
|
||||||
|
}
|
||||||
|
|
||||||
|
func TestMatrixRow(t *testing.T) {
|
||||||
|
assert := assert.New(t)
|
||||||
|
|
||||||
|
m := NewFromArrays([][]float64{
|
||||||
|
{1, 2, 3},
|
||||||
|
{4, 5, 6},
|
||||||
|
{7, 8, 9},
|
||||||
|
})
|
||||||
|
|
||||||
|
assert.Equal([]float64{1, 2, 3}, m.Row(0))
|
||||||
|
assert.Equal([]float64{4, 5, 6}, m.Row(1))
|
||||||
|
assert.Equal([]float64{7, 8, 9}, m.Row(2))
|
||||||
|
}
|
||||||
|
|
||||||
|
func TestMatrixSwapRows(t *testing.T) {
|
||||||
|
assert := assert.New(t)
|
||||||
|
|
||||||
|
m := NewFromArrays([][]float64{
|
||||||
|
{1, 2, 3},
|
||||||
|
{4, 5, 6},
|
||||||
|
{7, 8, 9},
|
||||||
|
})
|
||||||
|
|
||||||
|
m.SwapRows(0, 1)
|
||||||
|
|
||||||
|
assert.Equal([]float64{4, 5, 6}, m.Row(0))
|
||||||
|
assert.Equal([]float64{1, 2, 3}, m.Row(1))
|
||||||
|
assert.Equal([]float64{7, 8, 9}, m.Row(2))
|
||||||
|
}
|
||||||
|
|
||||||
|
func TestMatrixCopy(t *testing.T) {
|
||||||
|
assert := assert.New(t)
|
||||||
|
|
||||||
|
m := NewFromArrays([][]float64{
|
||||||
|
{1, 2, 3},
|
||||||
|
{4, 5, 6},
|
||||||
|
{7, 8, 9},
|
||||||
|
})
|
||||||
|
|
||||||
|
m2 := m.Copy()
|
||||||
|
assert.False(m == m2)
|
||||||
|
assert.True(m.Equals(m2))
|
||||||
|
}
|
||||||
|
|
||||||
|
func TestMatrixDiagonalVector(t *testing.T) {
|
||||||
|
assert := assert.New(t)
|
||||||
|
|
||||||
|
m := NewFromArrays([][]float64{
|
||||||
|
{1, 4, 7},
|
||||||
|
{4, 2, 8},
|
||||||
|
{7, 8, 3},
|
||||||
|
})
|
||||||
|
|
||||||
|
diag := m.DiagonalVector()
|
||||||
|
assert.Equal([]float64{1, 2, 3}, diag)
|
||||||
|
}
|
||||||
|
|
||||||
|
func TestMatrixDiagonalVectorLandscape(t *testing.T) {
|
||||||
|
assert := assert.New(t)
|
||||||
|
|
||||||
|
m := NewFromArrays([][]float64{
|
||||||
|
{1, 4, 7, 99},
|
||||||
|
{4, 2, 8, 99},
|
||||||
|
})
|
||||||
|
|
||||||
|
diag := m.DiagonalVector()
|
||||||
|
assert.Equal([]float64{1, 2}, diag)
|
||||||
|
}
|
||||||
|
|
||||||
|
func TestMatrixDiagonalVectorPortrait(t *testing.T) {
|
||||||
|
assert := assert.New(t)
|
||||||
|
|
||||||
|
m := NewFromArrays([][]float64{
|
||||||
|
{1, 4},
|
||||||
|
{4, 2},
|
||||||
|
{99, 99},
|
||||||
|
})
|
||||||
|
|
||||||
|
diag := m.DiagonalVector()
|
||||||
|
assert.Equal([]float64{1, 2}, diag)
|
||||||
|
}
|
||||||
|
|
||||||
|
func TestMatrixDiagonal(t *testing.T) {
|
||||||
|
assert := assert.New(t)
|
||||||
|
|
||||||
|
m := NewFromArrays([][]float64{
|
||||||
|
{1, 4, 7},
|
||||||
|
{4, 2, 8},
|
||||||
|
{7, 8, 3},
|
||||||
|
})
|
||||||
|
|
||||||
|
m2 := NewFromArrays([][]float64{
|
||||||
|
{1, 0, 0},
|
||||||
|
{0, 2, 0},
|
||||||
|
{0, 0, 3},
|
||||||
|
})
|
||||||
|
|
||||||
|
assert.True(m.Diagonal().Equals(m2))
|
||||||
|
}
|
||||||
|
|
||||||
|
func TestMatrixEquals(t *testing.T) {
|
||||||
|
assert := assert.New(t)
|
||||||
|
|
||||||
|
m := NewFromArrays([][]float64{
|
||||||
|
{1, 4, 7},
|
||||||
|
{4, 2, 8},
|
||||||
|
{7, 8, 3},
|
||||||
|
})
|
||||||
|
|
||||||
|
assert.False(m.Equals(nil))
|
||||||
|
var nilMatrix *Matrix
|
||||||
|
assert.True(nilMatrix.Equals(nil))
|
||||||
|
assert.False(m.Equals(New(1, 1)))
|
||||||
|
assert.False(m.Equals(New(3, 3)))
|
||||||
|
assert.True(m.Equals(New(3, 3, 1, 4, 7, 4, 2, 8, 7, 8, 3)))
|
||||||
|
}
|
||||||
|
|
||||||
|
func TestMatrixL(t *testing.T) {
|
||||||
|
assert := assert.New(t)
|
||||||
|
|
||||||
|
m := NewFromArrays([][]float64{
|
||||||
|
{1, 2, 3},
|
||||||
|
{4, 5, 6},
|
||||||
|
{7, 8, 9},
|
||||||
|
})
|
||||||
|
|
||||||
|
l := m.L()
|
||||||
|
assert.True(l.Equals(New(3, 3, 1, 2, 3, 0, 5, 6, 0, 0, 9)))
|
||||||
|
}
|
||||||
|
|
||||||
|
func TestMatrixU(t *testing.T) {
|
||||||
|
assert := assert.New(t)
|
||||||
|
|
||||||
|
m := NewFromArrays([][]float64{
|
||||||
|
{1, 2, 3},
|
||||||
|
{4, 5, 6},
|
||||||
|
{7, 8, 9},
|
||||||
|
})
|
||||||
|
|
||||||
|
u := m.U()
|
||||||
|
assert.True(u.Equals(New(3, 3, 0, 0, 0, 4, 0, 0, 7, 8, 0)))
|
||||||
|
}
|
||||||
|
|
||||||
|
func TestMatrixString(t *testing.T) {
|
||||||
|
assert := assert.New(t)
|
||||||
|
|
||||||
|
m := NewFromArrays([][]float64{
|
||||||
|
{1, 2, 3},
|
||||||
|
{4, 5, 6},
|
||||||
|
{7, 8, 9},
|
||||||
|
})
|
||||||
|
|
||||||
|
assert.Equal("1 2 3 \n4 5 6 \n7 8 9 \n", m.String())
|
||||||
|
}
|
||||||
|
|
||||||
|
func TestMatrixLU(t *testing.T) {
|
||||||
|
assert := assert.New(t)
|
||||||
|
|
||||||
|
m := NewFromArrays([][]float64{
|
||||||
|
{1, 3, 5},
|
||||||
|
{2, 4, 7},
|
||||||
|
{1, 1, 0},
|
||||||
|
})
|
||||||
|
|
||||||
|
l, u, p := m.LU()
|
||||||
|
assert.NotNil(l)
|
||||||
|
assert.NotNil(u)
|
||||||
|
assert.NotNil(p)
|
||||||
|
}
|
||||||
|
|
||||||
|
func TestMatrixQR(t *testing.T) {
|
||||||
|
assert := assert.New(t)
|
||||||
|
|
||||||
|
m := NewFromArrays([][]float64{
|
||||||
|
{12, -51, 4},
|
||||||
|
{6, 167, -68},
|
||||||
|
{-4, 24, -41},
|
||||||
|
})
|
||||||
|
|
||||||
|
q, r := m.QR()
|
||||||
|
assert.NotNil(q)
|
||||||
|
assert.NotNil(r)
|
||||||
|
}
|
||||||
|
|
||||||
|
func TestMatrixTranspose(t *testing.T) {
|
||||||
|
assert := assert.New(t)
|
||||||
|
|
||||||
|
m := NewFromArrays([][]float64{
|
||||||
|
{1, 2, 3},
|
||||||
|
{4, 5, 6},
|
||||||
|
{7, 8, 9},
|
||||||
|
{10, 11, 12},
|
||||||
|
})
|
||||||
|
|
||||||
|
m2 := m.Transpose()
|
||||||
|
|
||||||
|
rows, cols := m2.Size()
|
||||||
|
assert.Equal(3, rows)
|
||||||
|
assert.Equal(4, cols)
|
||||||
|
|
||||||
|
assert.Equal(1, m2.Get(0, 0))
|
||||||
|
assert.Equal(10, m2.Get(0, 3))
|
||||||
|
assert.Equal(3, m2.Get(2, 0))
|
||||||
|
}
|
45
matrix/regression.go
Normal file
45
matrix/regression.go
Normal file
|
@ -0,0 +1,45 @@
|
||||||
|
package matrix
|
||||||
|
|
||||||
|
import "errors"
|
||||||
|
|
||||||
|
var (
|
||||||
|
// ErrPolyRegArraysSameLength is a common error.
|
||||||
|
ErrPolyRegArraysSameLength = errors.New("polynomial array inputs must be the same length")
|
||||||
|
)
|
||||||
|
|
||||||
|
// Poly returns the polynomial regress of a given degree over the given values.
|
||||||
|
func Poly(xvalues, yvalues []float64, degree int) ([]float64, error) {
|
||||||
|
if len(xvalues) != len(yvalues) {
|
||||||
|
return nil, ErrPolyRegArraysSameLength
|
||||||
|
}
|
||||||
|
|
||||||
|
m := len(yvalues)
|
||||||
|
n := degree + 1
|
||||||
|
y := New(m, 1, yvalues...)
|
||||||
|
x := Zero(m, n)
|
||||||
|
|
||||||
|
for i := 0; i < m; i++ {
|
||||||
|
ip := float64(1)
|
||||||
|
for j := 0; j < n; j++ {
|
||||||
|
x.Set(i, j, ip)
|
||||||
|
ip *= xvalues[i]
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
q, r := x.QR()
|
||||||
|
qty, err := q.Transpose().Times(y)
|
||||||
|
if err != nil {
|
||||||
|
return nil, err
|
||||||
|
}
|
||||||
|
|
||||||
|
c := make([]float64, n)
|
||||||
|
for i := n - 1; i >= 0; i-- {
|
||||||
|
c[i] = qty.Get(i, 0)
|
||||||
|
for j := i + 1; j < n; j++ {
|
||||||
|
c[i] -= c[j] * r.Get(i, j)
|
||||||
|
}
|
||||||
|
c[i] /= r.Get(i, i)
|
||||||
|
}
|
||||||
|
|
||||||
|
return c, nil
|
||||||
|
}
|
22
matrix/regression_test.go
Normal file
22
matrix/regression_test.go
Normal file
|
@ -0,0 +1,22 @@
|
||||||
|
package matrix
|
||||||
|
|
||||||
|
import (
|
||||||
|
"testing"
|
||||||
|
|
||||||
|
assert "github.com/blendlabs/go-assert"
|
||||||
|
)
|
||||||
|
|
||||||
|
func TestPoly(t *testing.T) {
|
||||||
|
assert := assert.New(t)
|
||||||
|
var xGiven = []float64{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}
|
||||||
|
var yGiven = []float64{1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321}
|
||||||
|
var degree = 2
|
||||||
|
|
||||||
|
c, err := Poly(xGiven, yGiven, degree)
|
||||||
|
assert.Nil(err)
|
||||||
|
assert.Len(c, 3)
|
||||||
|
|
||||||
|
assert.InDelta(c[0], 0.999999999, DefaultEpsilon)
|
||||||
|
assert.InDelta(c[1], 2, DefaultEpsilon)
|
||||||
|
assert.InDelta(c[2], 3, DefaultEpsilon)
|
||||||
|
}
|
36
matrix/util.go
Normal file
36
matrix/util.go
Normal file
|
@ -0,0 +1,36 @@
|
||||||
|
package matrix
|
||||||
|
|
||||||
|
import (
|
||||||
|
"math"
|
||||||
|
"strconv"
|
||||||
|
)
|
||||||
|
|
||||||
|
func minInt(values ...int) int {
|
||||||
|
min := math.MaxInt32
|
||||||
|
|
||||||
|
for x := 0; x < len(values); x++ {
|
||||||
|
if values[x] < min {
|
||||||
|
min = values[x]
|
||||||
|
}
|
||||||
|
}
|
||||||
|
return min
|
||||||
|
}
|
||||||
|
|
||||||
|
func maxInt(values ...int) int {
|
||||||
|
max := math.MinInt32
|
||||||
|
|
||||||
|
for x := 0; x < len(values); x++ {
|
||||||
|
if values[x] > max {
|
||||||
|
max = values[x]
|
||||||
|
}
|
||||||
|
}
|
||||||
|
return max
|
||||||
|
}
|
||||||
|
|
||||||
|
func f64s(v float64) string {
|
||||||
|
return strconv.FormatFloat(v, 'f', -1, 64)
|
||||||
|
}
|
||||||
|
|
||||||
|
func roundToEpsilon(value, epsilon float64) float64 {
|
||||||
|
return math.Nextafter(value, value)
|
||||||
|
}
|
17
matrix/vector.go
Normal file
17
matrix/vector.go
Normal file
|
@ -0,0 +1,17 @@
|
||||||
|
package matrix
|
||||||
|
|
||||||
|
// Vector is just an array of values.
|
||||||
|
type Vector []float64
|
||||||
|
|
||||||
|
// DotProduct returns the dot product of two vectors.
|
||||||
|
func (v Vector) DotProduct(v2 Vector) (result float64, err error) {
|
||||||
|
if len(v) != len(v2) {
|
||||||
|
err = ErrDimensionMismatch
|
||||||
|
return
|
||||||
|
}
|
||||||
|
|
||||||
|
for i := 0; i < len(v); i++ {
|
||||||
|
result = result + (v[i] * v2[i])
|
||||||
|
}
|
||||||
|
return
|
||||||
|
}
|
135
polynomial_regression_series.go
Normal file
135
polynomial_regression_series.go
Normal file
|
@ -0,0 +1,135 @@
|
||||||
|
package chart
|
||||||
|
|
||||||
|
import (
|
||||||
|
"fmt"
|
||||||
|
"math"
|
||||||
|
|
||||||
|
"github.com/wcharczuk/go-chart/matrix"
|
||||||
|
)
|
||||||
|
|
||||||
|
// PolynomialRegressionSeries implements a polynomial regression over a given
|
||||||
|
// inner series.
|
||||||
|
type PolynomialRegressionSeries struct {
|
||||||
|
Name string
|
||||||
|
Style Style
|
||||||
|
YAxis YAxisType
|
||||||
|
|
||||||
|
Limit int
|
||||||
|
Offset int
|
||||||
|
Degree int
|
||||||
|
InnerSeries ValueProvider
|
||||||
|
|
||||||
|
coeffs []float64
|
||||||
|
}
|
||||||
|
|
||||||
|
// GetName returns the name of the time series.
|
||||||
|
func (prs PolynomialRegressionSeries) GetName() string {
|
||||||
|
return prs.Name
|
||||||
|
}
|
||||||
|
|
||||||
|
// GetStyle returns the line style.
|
||||||
|
func (prs PolynomialRegressionSeries) GetStyle() Style {
|
||||||
|
return prs.Style
|
||||||
|
}
|
||||||
|
|
||||||
|
// GetYAxis returns which YAxis the series draws on.
|
||||||
|
func (prs PolynomialRegressionSeries) GetYAxis() YAxisType {
|
||||||
|
return prs.YAxis
|
||||||
|
}
|
||||||
|
|
||||||
|
// Len returns the number of elements in the series.
|
||||||
|
func (prs PolynomialRegressionSeries) Len() int {
|
||||||
|
return Math.MinInt(prs.GetLimit(), prs.InnerSeries.Len()-prs.GetOffset())
|
||||||
|
}
|
||||||
|
|
||||||
|
// GetLimit returns the window size.
|
||||||
|
func (prs PolynomialRegressionSeries) GetLimit() int {
|
||||||
|
if prs.Limit == 0 {
|
||||||
|
return prs.InnerSeries.Len()
|
||||||
|
}
|
||||||
|
return prs.Limit
|
||||||
|
}
|
||||||
|
|
||||||
|
// GetEndIndex returns the effective limit end.
|
||||||
|
func (prs PolynomialRegressionSeries) GetEndIndex() int {
|
||||||
|
offset := prs.GetOffset() + prs.Len()
|
||||||
|
innerSeriesLastIndex := prs.InnerSeries.Len() - 1
|
||||||
|
return Math.MinInt(offset, innerSeriesLastIndex)
|
||||||
|
}
|
||||||
|
|
||||||
|
// GetOffset returns the data offset.
|
||||||
|
func (prs PolynomialRegressionSeries) GetOffset() int {
|
||||||
|
if prs.Offset == 0 {
|
||||||
|
return 0
|
||||||
|
}
|
||||||
|
return prs.Offset
|
||||||
|
}
|
||||||
|
|
||||||
|
// Validate validates the series.
|
||||||
|
func (prs *PolynomialRegressionSeries) Validate() error {
|
||||||
|
if prs.InnerSeries == nil {
|
||||||
|
return fmt.Errorf("linear regression series requires InnerSeries to be set")
|
||||||
|
}
|
||||||
|
|
||||||
|
endIndex := prs.GetEndIndex()
|
||||||
|
if endIndex >= prs.InnerSeries.Len() {
|
||||||
|
return fmt.Errorf("invalid window; inner series has length %d but end index is %d", prs.InnerSeries.Len(), endIndex)
|
||||||
|
}
|
||||||
|
|
||||||
|
return nil
|
||||||
|
}
|
||||||
|
|
||||||
|
// GetValue returns the series value for a given index.
|
||||||
|
func (prs *PolynomialRegressionSeries) GetValue(index int) (x, y float64) {
|
||||||
|
if prs.InnerSeries == nil || prs.InnerSeries.Len() == 0 {
|
||||||
|
return
|
||||||
|
}
|
||||||
|
|
||||||
|
if prs.coeffs == nil {
|
||||||
|
coeffs, err := prs.computeCoefficients()
|
||||||
|
if err != nil {
|
||||||
|
panic(err)
|
||||||
|
}
|
||||||
|
prs.coeffs = coeffs
|
||||||
|
}
|
||||||
|
|
||||||
|
offset := prs.GetOffset()
|
||||||
|
effectiveIndex := Math.MinInt(index+offset, prs.InnerSeries.Len())
|
||||||
|
x, y = prs.InnerSeries.GetValue(effectiveIndex)
|
||||||
|
y = prs.apply(x)
|
||||||
|
return
|
||||||
|
}
|
||||||
|
|
||||||
|
func (prs *PolynomialRegressionSeries) apply(v float64) (out float64) {
|
||||||
|
for index, coeff := range prs.coeffs {
|
||||||
|
out = out + (coeff * math.Pow(v, float64(index)))
|
||||||
|
}
|
||||||
|
return
|
||||||
|
}
|
||||||
|
|
||||||
|
func (prs *PolynomialRegressionSeries) computeCoefficients() ([]float64, error) {
|
||||||
|
xvalues, yvalues := prs.values()
|
||||||
|
return matrix.Poly(xvalues, yvalues, prs.Degree)
|
||||||
|
}
|
||||||
|
|
||||||
|
func (prs *PolynomialRegressionSeries) values() (xvalues, yvalues []float64) {
|
||||||
|
startIndex := prs.GetOffset()
|
||||||
|
endIndex := prs.GetEndIndex()
|
||||||
|
|
||||||
|
xvalues = make([]float64, endIndex-startIndex)
|
||||||
|
yvalues = make([]float64, endIndex-startIndex)
|
||||||
|
|
||||||
|
for index := startIndex; index < endIndex; index++ {
|
||||||
|
x, y := prs.InnerSeries.GetValue(index)
|
||||||
|
xvalues[index] = x
|
||||||
|
yvalues[index] = y
|
||||||
|
}
|
||||||
|
|
||||||
|
return
|
||||||
|
}
|
||||||
|
|
||||||
|
// Render renders the series.
|
||||||
|
func (prs *PolynomialRegressionSeries) Render(r Renderer, canvasBox Box, xrange, yrange Range, defaults Style) {
|
||||||
|
style := prs.Style.InheritFrom(defaults)
|
||||||
|
Draw.LineSeries(r, canvasBox, xrange, yrange, style, prs)
|
||||||
|
}
|
Loading…
Reference in a new issue