From a211e88530d4ca4a6efc5cc05fc48c581ac06a2f Mon Sep 17 00:00:00 2001 From: Will Charczuk Date: Tue, 18 Apr 2017 20:20:29 -0700 Subject: [PATCH] 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. --- _examples/poly_regression/main.go | 41 +++ linear_regression_series.go | 15 +- linear_regression_series_test.go | 2 +- matrix/matrix.go | 592 ++++++++++++++++++++++++++++++ matrix/matrix_test.go | 396 ++++++++++++++++++++ matrix/regression.go | 45 +++ matrix/regression_test.go | 22 ++ matrix/util.go | 36 ++ matrix/vector.go | 17 + polynomial_regression_series.go | 135 +++++++ 10 files changed, 1292 insertions(+), 9 deletions(-) create mode 100644 _examples/poly_regression/main.go create mode 100644 matrix/matrix.go create mode 100644 matrix/matrix_test.go create mode 100644 matrix/regression.go create mode 100644 matrix/regression_test.go create mode 100644 matrix/util.go create mode 100644 matrix/vector.go create mode 100644 polynomial_regression_series.go diff --git a/_examples/poly_regression/main.go b/_examples/poly_regression/main.go new file mode 100644 index 0000000..4ed25bf --- /dev/null +++ b/_examples/poly_regression/main.go @@ -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) +} diff --git a/linear_regression_series.go b/linear_regression_series.go index 142ea55..f1e3bc5 100644 --- a/linear_regression_series.go +++ b/linear_regression_series.go @@ -9,7 +9,7 @@ type LinearRegressionSeries struct { Style Style YAxis YAxisType - Window int + Limit int Offset int InnerSeries ValueProvider @@ -36,18 +36,18 @@ func (lrs LinearRegressionSeries) GetYAxis() YAxisType { // Len returns the number of elements in the series. 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. -func (lrs LinearRegressionSeries) GetWindow() int { - if lrs.Window == 0 { +// GetLimit returns the window size. +func (lrs LinearRegressionSeries) GetLimit() int { + if lrs.Limit == 0 { 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 { offset := lrs.GetOffset() + lrs.Len() innerSeriesLastIndex := lrs.InnerSeries.Len() - 1 @@ -105,7 +105,6 @@ func (lrs *LinearRegressionSeries) computeCoefficients() { xvalues := NewRingBufferWithCapacity(lrs.Len()) for index := startIndex; index < endIndex; index++ { - x, _ := lrs.InnerSeries.GetValue(index) xvalues.Enqueue(x) } diff --git a/linear_regression_series_test.go b/linear_regression_series_test.go index 4a72669..b8c39e1 100644 --- a/linear_regression_series_test.go +++ b/linear_regression_series_test.go @@ -62,7 +62,7 @@ func TestLinearRegressionSeriesWindowAndOffset(t *testing.T) { linRegSeries := &LinearRegressionSeries{ InnerSeries: mainSeries, Offset: 10, - Window: 10, + Limit: 10, } assert.Equal(10, linRegSeries.Len()) diff --git a/matrix/matrix.go b/matrix/matrix.go new file mode 100644 index 0000000..5fa857b --- /dev/null +++ b/matrix/matrix.go @@ -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 +} diff --git a/matrix/matrix_test.go b/matrix/matrix_test.go new file mode 100644 index 0000000..bc896be --- /dev/null +++ b/matrix/matrix_test.go @@ -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)) +} diff --git a/matrix/regression.go b/matrix/regression.go new file mode 100644 index 0000000..7aaea7c --- /dev/null +++ b/matrix/regression.go @@ -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 +} diff --git a/matrix/regression_test.go b/matrix/regression_test.go new file mode 100644 index 0000000..c55a480 --- /dev/null +++ b/matrix/regression_test.go @@ -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) +} diff --git a/matrix/util.go b/matrix/util.go new file mode 100644 index 0000000..5f15779 --- /dev/null +++ b/matrix/util.go @@ -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) +} diff --git a/matrix/vector.go b/matrix/vector.go new file mode 100644 index 0000000..c8aa380 --- /dev/null +++ b/matrix/vector.go @@ -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 +} diff --git a/polynomial_regression_series.go b/polynomial_regression_series.go new file mode 100644 index 0000000..70498c4 --- /dev/null +++ b/polynomial_regression_series.go @@ -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) +}