Skip to content

Commit

Permalink
feat: optimization matrix inverse algorithm and release pre version 🚀
Browse files Browse the repository at this point in the history
  • Loading branch information
lvisei committed Dec 13, 2020
1 parent 6e0a806 commit bd43c64
Show file tree
Hide file tree
Showing 24 changed files with 4,323 additions and 456 deletions.
4 changes: 2 additions & 2 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -16,8 +16,8 @@

.idea

temp/
tempdata/
.temp/
testdata/
data/

internal/service
64 changes: 62 additions & 2 deletions README.md
Original file line number Diff line number Diff line change
@@ -1,5 +1,65 @@
# go-kriging

[![GoDoc](https://godoc.org/github.com/liuvigongzuoshi/go-kriging?status.svg)](https://godoc.org/github.com/liuvigongzuoshi/go-kriging)
[![GoDoc](https://godoc.org/github.com/liuvigongzuoshi/go-kriging?status.svg)](https://pkg.go.dev/github.com/liuvigongzuoshi/go-kriging)

Golang Multi-Goroutine spatial interpolation algorithm library for geospatial prediction and mapping via ordinary kriging
Golang Multi-Goroutine spatial interpolation algorithm library for geospatial prediction and mapping via ordinary kriging.

Based on [oeo4b/kriging.js](https://github.com/oeo4b/kriging.js) refactoring and optimized the algorithm and added some new features.

## Fitting a Model

The train method with the new ordinaryKriging fits your input to whatever variogram model you specify - gaussian, exponential or spherical - and returns a variogram variable.


```go
import "github.com/liuvigongzuoshi/go-kriging/ordinarykriging"

func main() {
sigma2 := 0
alpha := 100
ordinaryKriging := ordinarykriging.NewOrdinary(values, x, y)
variogram = ordinaryKriging.Train(ordinarykriging.Spherical, sigma2, alpha)
}
```

## Predicting New Values

Values can be predicted for new coordinate pairs by using the predict method with the new ordinaryKriging.

```go
import "github.com/liuvigongzuoshi/go-kriging/ordinarykriging"

func main() {
// ...
// Pair of new coordinates to predict
xnew := 0.5481
ynew := 0.4455
tpredicted := ordinaryKriging.predict(xnew, ynew)
}
```

## 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:

- 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 } )
- Spherical: k(a,b) = w[0] + w[1] * ( 1.5 * ( ||a-b|| / range ) - 0.5 * ( ||a-b|| / range )3 )

The variance parameter α of the prior distribution for w should be manually set, according to:

- w ~ N(w|0, αI)

Using the fitted kernel function hyperparameters and setting K as the Gram matrix, the prior and likelihood for the gaussian process become:

- y ~ N(y|0, K)
- t|y ~ N(t|y, σ2I)

The variance parameter σ2 of the likelihood reflects the error in the gaussian process and should be manually set.


## Other

[kriging-wasm example](https://github.com/liuvigongzuoshi/kriging-wasm) - Test example used by wasm compiled with go-kriging algorithm code.

[go-kriging-service](https://github.com/liuvigongzuoshi/go-kriging-service) - Call the REST service written by the go-kriging algorithm package, which supports concurrent calls by multiple users, and has a simple logging and fault-tolerant recovery mechanism.
56 changes: 30 additions & 26 deletions examples/csv/main.go
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@ package main
import (
"encoding/csv"
"fmt"

"image/color"
"image/png"
"io/ioutil"
Expand All @@ -15,18 +16,21 @@ import (
"runtime/pprof"

"github.com/liuvigongzuoshi/go-kriging/canvas"
"github.com/liuvigongzuoshi/go-kriging/ordinary"
"github.com/liuvigongzuoshi/go-kriging/ordinarykriging"
"github.com/liuvigongzuoshi/go-kriging/pkg/json"
)

const dirPath = "testdata"
const testDataDirPath = "testdata"
const tempDataDirPath = "tempdata"
const cpuProfileFilePath = tempDataDirPath + "/cpu_profile"
const memProfileFilePath = tempDataDirPath + "/mem_profile"

func main() {
cpuProfile, _ := os.Create("./testdata/cpu_profile")
cpuProfile, _ := os.Create(cpuProfileFilePath)
if err := pprof.StartCPUProfile(cpuProfile); err != nil {
log.Fatal(err)
}
//memProfile, _ := os.Create("./testdata/mem_profile")
//memProfile, _ := os.Create(memProfileFilePath)
//if err := pprof.WriteHeapProfile(memProfile); err != nil {
// log.Fatal(err)
//}
Expand All @@ -36,19 +40,20 @@ func main() {
//memProfile.Close()
}()

data, err := readCsvFile("examples/csv/data/2045.csv")
data, err := readCsvFile("examples/csv/testdata/2045.csv")
if err != nil {
log.Fatal(err)
}
polygon, err := readGeoJsonFile("examples/csv/data/yn.json")
polygon, err := readGeoJsonFile("examples/csv/testdata/yn.json")
if err != nil {
log.Fatal(err)
}
defer timeCost()("训练模型加插值总耗时")
defer timeCost()("训练模型与插值生成网格图片总耗时")

ordinaryKriging := ordinary.NewOrdinary(data["values"], data["x"], data["y"])
_ = ordinaryKriging.Train(ordinary.Spherical, 0, 100)
ordinaryKriging := ordinarykriging.NewOrdinary(data["values"], data["x"], data["y"])
_ = ordinaryKriging.Train(ordinarykriging.Spherical, 0, 100)

_ = polygon
gridPlot(ordinaryKriging, polygon)

//var wg sync.WaitGroup
Expand All @@ -61,18 +66,17 @@ func main() {
// defer wg.Done()
// contourRectanglePng(ordinaryKriging)
//}()

//wg.Wait()
}

func gridPlot(ordinaryKriging *ordinary.Variogram, polygon ordinary.PolygonCoordinates) {
func gridPlot(ordinaryKriging *ordinarykriging.Variogram, polygon ordinarykriging.PolygonCoordinates) {
defer timeCost()("插值生成网格图片耗时")
gridMatrices := ordinaryKriging.Grid(polygon, 0.01)
ctx := ordinaryKriging.Plot(gridMatrices, 500, 500, gridMatrices.Xlim, gridMatrices.Ylim, ordinary.DefaultGridLevelColor)
ctx := ordinaryKriging.Plot(gridMatrices, 500, 500, gridMatrices.Xlim, gridMatrices.Ylim, ordinarykriging.DefaultGridLevelColor)

subTitle := &canvas.TextConfig{
Text: "球面半变异函数模型",
FontName: "testdata/fonts/source-han-sans-sc/regular.ttf",
FontName: testDataDirPath + "/fonts/source-han-sans-sc/regular.ttf",
FontSize: 28,
Color: color.RGBA{R: 0, G: 0, B: 0, A: 255},
OffsetX: 252,
Expand All @@ -87,18 +91,18 @@ func gridPlot(ordinaryKriging *ordinary.Variogram, polygon ordinary.PolygonCoord
if err != nil {
log.Fatal(err)
} else {
saveBuffer("grid.png", buffer)
saveBufferFile("grid.png", buffer)
}

//writeFile("gridMatrices.json", gridMatrices)
}

func contourRectanglePng(ordinaryKriging *ordinary.Variogram) {
func contourRectanglePng(ordinaryKriging *ordinarykriging.Variogram) {
defer timeCost()("插值生成矩形图片耗时")
xWidth, yWidth := 800, 800
contourRectangle := ordinaryKriging.Contour(xWidth, yWidth)
pngPath := fmt.Sprintf("%v/%v.png", dirPath, time.Now().Format("2006-01-02 15:04:05"))
ctx := ordinaryKriging.PlotRectangleContour(contourRectangle, 500, 500, contourRectangle.Xlim, contourRectangle.Ylim, ordinary.DefaultLegendColor)
pngPath := fmt.Sprintf("%v/%v.png", tempDataDirPath, time.Now().Format("2006-01-02 15:04:05"))
ctx := ordinaryKriging.PlotRectangleContour(contourRectangle, 500, 500, contourRectangle.Xlim, contourRectangle.Ylim, ordinarykriging.DefaultLegendColor)
img := ordinaryKriging.PlotPng(contourRectangle)

err := os.MkdirAll(filepath.Dir(pngPath), os.ModePerm)
Expand All @@ -116,7 +120,7 @@ func contourRectanglePng(ordinaryKriging *ordinary.Variogram) {
if err != nil {
log.Fatal(err)
} else {
saveBuffer("rectangle.png", buffer)
saveBufferFile("rectangle.png", buffer)
}
}

Expand All @@ -134,7 +138,7 @@ func readCsvFile(filePath string) (map[string][]float64, error) {

records, err := csv.NewReader(f).ReadAll()
if err != nil {
log.Fatal("Unable to parse file as CSV for "+filePath, err)
log.Fatal("Unable to parse file as CSV for\n "+filePath, err)
return nil, err
}

Expand All @@ -159,18 +163,18 @@ func readCsvFile(filePath string) (map[string][]float64, error) {
data := map[string][]float64{"values": values, "x": lons, "y": lats}

//fmt.Printf("values %#v\n lons %#v\n lats %#v\n", values, lons, lats)
// writeFile("data.json", data)
//writeFile("tempdata.json", tempdata)

return data, nil
}

func readGeoJsonFile(filePath string) (ordinary.PolygonCoordinates, error) {
func readGeoJsonFile(filePath string) (ordinarykriging.PolygonCoordinates, error) {
content, err := ioutil.ReadFile(filePath)
if err != nil {
log.Fatal("Unable to read input file "+filePath, err)
log.Fatal("Unable to read input file \n"+filePath, err)
return nil, err
}
var polygonGeometry ordinary.PolygonGeometry
var polygonGeometry ordinarykriging.PolygonGeometry
err = json.Unmarshal(content, &polygonGeometry)
if err != nil {
log.Fatalf("invalid json: %v", err)
Expand All @@ -189,7 +193,7 @@ func timeCost() func(name string) {
}

func writeFile(fileName string, v interface{}) {
filePath := fmt.Sprintf("%v/%v %v", dirPath, time.Now().Format("2006-01-02 15-04-05"), fileName)
filePath := fmt.Sprintf("%v/%v %v", tempDataDirPath, time.Now().Format("2006-01-02 15-04-05"), fileName)
fmt.Printf("%#v\n", filePath)
// fmt.Printf("%#v\n", v)
content, err := json.Marshal(v)
Expand All @@ -199,8 +203,8 @@ func writeFile(fileName string, v interface{}) {
ioutil.WriteFile(filePath, content, os.ModePerm)
}

func saveBuffer(fileName string, content []byte) {
filePath := fmt.Sprintf("%v/%v %v", dirPath, time.Now().Format("2006-01-02 15-04-05"), fileName)
func saveBufferFile(fileName string, content []byte) {
filePath := fmt.Sprintf("%v/%v %v", tempDataDirPath, time.Now().Format("2006-01-02 15-04-05"), fileName)
fmt.Printf("%#v\n", filePath)
ioutil.WriteFile(filePath, content, os.ModePerm)
}
File renamed without changes.
File renamed without changes.
File renamed without changes.
Loading

0 comments on commit bd43c64

Please sign in to comment.