Skip to content

Commit

Permalink
fix: matrix Inverse err, doc: note url
Browse files Browse the repository at this point in the history
lvisei committed Dec 16, 2020
1 parent b501d03 commit be75f1e
Showing 3 changed files with 21 additions and 33 deletions.
2 changes: 1 addition & 1 deletion README.md
Original file line number Diff line number Diff line change
@@ -40,7 +40,7 @@ func main() {

## Variogram and Probability Model

The various variogram models can be interpreted as kernel functions for 2-dimensional coordinates a, b and parameters nugget, range, sill and A. Reparameterized as a linear function, with w = [nugget, (sill-nugget)/range], this becomes:
According to [sakitam-gis](https://sakitam-gis.github.io/kriging.js/examples/world.html), the various variogram models can be interpreted as kernel functions for 2-dimensional coordinates a, b and parameters nugget, range, sill and A. Reparameterized as a linear function, with w = [nugget, (sill-nugget)/range], this becomes:

- Gaussian: k(a,b) = w[0] + w[1] * ( 1 - exp{ -( ||a-b|| / range )2 / A } )
- Exponential: k(a,b) = w[0] + w[1] * ( 1 - exp{ -( ||a-b|| / range ) / A } )
50 changes: 19 additions & 31 deletions ordinarykriging/matrix-inverse.go
Original file line number Diff line number Diff line change
@@ -102,33 +102,22 @@ func gaussJordanInversion(X []float64, n int) bool {
func matrixInverseByCol(a [][]float64) ([][]float64, bool) {
// eMat 返回 n 阶单位矩阵
var eMat = func(n int) ([][]float64, bool) {
/*
返回n阶单位矩阵
输入 :
n 阶数
输出 :
sol 解值
err 解出标志:false-未解出或达到步数上限;
true-全部解出
*/
sol := make([][]float64, n)
for i := 0; i < n; i++ {
sol[i] = make([]float64, n)
}
var err bool = false

//判断阶数
// 判断阶数
if n < 1 {
return nil, err
return nil, false
}

//分配元素
for i := 0; i < n; i++ {
sol[i][i] = 1.0
}

err = true
return sol, err
return sol, true
}
// maxAbs 向量第一个绝对值最大值及其位置
var maxAbs = func(a []float64) (float64, int, bool) {
@@ -150,18 +139,19 @@ func matrixInverseByCol(a [][]float64) ([][]float64, bool) {
return sol, ii, err
}

var err bool = false
n := len(a)
temp0, _ := eMat(n)
b := temp0
sol := b
temp1 := make([]float64, n)

//判断是否方阵
if len(a) != len(a[0]) {
return sol, err
return nil, false
}

n := len(a)
sol, isOk := eMat(n)
if !isOk {
return nil, false
}

temp1 := make([]float64, n)

//主元消去
for i := 0; i < n; i++ {
//求第i列的主元素并调整行顺序
@@ -174,33 +164,31 @@ func matrixInverseByCol(a [][]float64) ([][]float64, bool) {
temp1 = a[ii+i]
a[ii+i] = a[i]
a[i] = temp1
temp1 = b[ii+i]
b[ii+i] = b[i]
b[i] = temp1
temp1 = sol[ii+i]
sol[ii+i] = sol[i]
sol[i] = temp1
}

//列消去
//本行主元置一
mul := a[i][i]
for j := 0; j < n; j++ {
a[i][j] = a[i][j] / mul
b[i][j] = b[i][j] / mul
sol[i][j] = sol[i][j] / mul
}
//其它列置零
for j := 0; j < n; j++ {
if j != i {
mul = a[j][i] / a[i][i]
for k := 0; k < n; k++ {
a[j][k] = a[j][k] - a[i][k]*mul
b[j][k] = b[j][k] - b[i][k]*mul
sol[j][k] = sol[j][k] - sol[i][k]*mul
}
}
}
}

sol = b
err = true
return sol, err
return sol, true
}

func matrixInverseByCol_(x []float64, n int) ([]float64, bool) {
@@ -228,7 +216,7 @@ func matrixInverse(x []float64, n int) ([]float64, bool) {
// Take the inverse of a and place the result in ia.
err := ia.Inverse(a)
if err != nil {
return x, false
return ia.RawMatrix().Data, false
}

return ia.RawMatrix().Data, true
2 changes: 1 addition & 1 deletion ordinarykriging/ordinarykriging.go
Original file line number Diff line number Diff line change
@@ -227,7 +227,7 @@ func (variogram *Variogram) Train(model ModelType, sigma2 float64, alpha float64
}

// Copy unprojected inverted matrix as K 复制未投影的逆矩阵为K
K = C
copy(K, C)
var M = matrixMultiply(C, variogram.t, n, n, 1)
variogram.K = K
variogram.M = M

0 comments on commit be75f1e

Please sign in to comment.