From c9e4c13590f5452abc3de0ca0af66fed03925e78 Mon Sep 17 00:00:00 2001 From: "facat@lab.com" Date: Tue, 9 Dec 2014 16:37:20 +0800 Subject: [PATCH] =?UTF-8?q?=E7=BB=8F=E8=BF=87=E6=97=A0=E6=95=B0=E4=B8=AAbu?= =?UTF-8?q?g=E7=9A=84=E6=B6=88=E7=81=AD=EF=BC=8C=E6=94=B6=E6=95=9B?= =?UTF-8?q?=E4=BA=86=E3=80=82?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: facat@lab.com --- .gitignore | 1 + m.txt | 1 - pam.go | 179 ++++++++++++++++++++++++++++++++++++++++++++++++++--- 3 files changed, 173 insertions(+), 8 deletions(-) create mode 100644 .gitignore diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..adb36c8 --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +*.exe \ No newline at end of file diff --git a/m.txt b/m.txt index e780a34..4b44e9a 100644 --- a/m.txt +++ b/m.txt @@ -74,7 +74,6 @@ while 1 if J==cadSetS(I) continue end - matDistance=repmat(d,1,length(SetS)); minD=matDistance-data(:,SetS); minD=diag(minD'*minD).^.5; diff --git a/pam.go b/pam.go index 1645a80..16590d9 100644 --- a/pam.go +++ b/pam.go @@ -3,6 +3,7 @@ package main import ( "fmt" + "log" "math" "math/rand" "time" @@ -11,9 +12,9 @@ import ( //先设置一些参数 const ( - dataN = 10 + dataN = 50 dim = 3 - clusterN = 3 + clusterN = 5 ) type Ctx struct { @@ -35,13 +36,13 @@ func initData() [][]float64 { //初始数据集 return r } -func copyIntMap(src map[int]interface{}, dst map[int]interface{}) { +func copyIntMap(src map[int]interface{}, dst *map[int]interface{}) { for k, v := range src { - dst[k] = v + (*dst)[k] = v } } -func minDistance(vec []float64, mat [][]float64) float64 { +func minDistance2Mat(vec []float64, mat [][]float64) float64 { var r float64 r = 1e20 var t float64 @@ -51,24 +52,92 @@ func minDistance(vec []float64, mat [][]float64) float64 { t += math.Pow(vec[v]-mat[m][v], 2) } t = math.Sqrt(t) + //log.Println("当前距离", t) if t < r { r = t } } + //log.Println("最小距离", r) return r } -func subData(set map[int]interface{}, ctx *Ctx) [][]float64 { +func minDistance2Vec(vec []float64, mat []float64) float64 { + t := make([][]float64, 1) + t[0] = mat + return minDistance2Mat(vec, t) +} + +func min2minDistance(vec []float64, mat [][]float64) (min float64, min2 float64) { + var r float64 + r = 1e20 + var r2 float64 + r2 = r + var t float64 + for m := 0; m < len(mat); m++ { + t = 0 + for v := 0; v < len(vec); v++ { + t += math.Pow(vec[v]-mat[m][v], 2) + } + t = math.Sqrt(t) + //log.Println("当前距离", t) + if t < r { + r2 = r + r = t + } else if t < r2 { + r2 = t + } + + } + //log.Println("最小距离", r) + min = r + min2 = r2 + //log.Println("r r2 t", r, r2, t) + return +} + +func subDataMap(set map[int]interface{}, ctx *Ctx) [][]float64 { r := make([][]float64, len(set)) ind := 0 for k, _ := range set { r[ind] = ctx.data[k] ind++ } + return r +} + +func subDataInt(i int, ctx *Ctx) []float64 { + //r := make([][]float64, 1) + //r[0] = ctx.data[i] + return ctx.data[i] +} + +func ones(m int, n int, val float64) [][]float64 { + r := make([][]float64, n) + for i := 0; i < n; i++ { + r[i] = make([]float64, m) + for j := 0; j < m; j++ { + r[i][j] = val + } + } + return r +} + +func maxInd(vec []float64) int { + var r float64 + r = -1.1e20 + var ri int + for i, v := range vec { + if v > r { + r = v + ri = i + } + } + return ri } func selectInitCluster(ctx *Ctx) { for clusterI := 0; clusterI < clusterN-1; clusterI++ { + maxG := ones(dataN, 1, -100) //等价于maxG=-100*ones(dataN,1); for cluster := 0; cluster < dataN; cluster++ { if _, has := ctx.SetS[cluster]; has { continue @@ -76,7 +145,7 @@ func selectInitCluster(ctx *Ctx) { //选一个候选数据 cddtI := cluster cadSetS := make(map[int]interface{}) - copyIntMap(cadSetS, ctx.SetS) + copyIntMap(ctx.SetS, &cadSetS) //src,dst cadSetS[cddtI] = nil Cij := 0.0 for J := 0; J < dataN; J++ { @@ -85,9 +154,102 @@ func selectInitCluster(ctx *Ctx) { } d := ctx.data[J] //寻找最短距离 + minD := minDistance2Mat(d, subDataMap(ctx.SetS, ctx)) + distanceIJ := minDistance2Vec(d, subDataInt(cddtI, ctx)) + if minD-distanceIJ > 0 { + Cij += minD - distanceIJ + } + } + maxG[0][cluster] = Cij + } + maxGInd := maxInd(maxG[0]) //等价于maxGInd=find(maxG==max(maxG)); + ctx.SetS[maxGInd] = nil + } +} + +func getUSetS(ctx *Ctx) map[int]interface{} { + r := make(map[int]interface{}) + for i := 0; i < dataN; i++ { + if _, has := ctx.SetS[i]; !has { + r[i] = nil + } + } + return r +} + +func swap(ctx *Ctx) { + USetS := getUSetS(ctx) + var minKjih float64 + + var minSetS map[int]interface{} + var minUSetS map[int]interface{} + for { + minKjih = 1e20 + for I, _ := range ctx.SetS { + for H, _ := range USetS { + cadSetS := make(map[int]interface{}) + copyIntMap(ctx.SetS, &cadSetS) //src,dst + cadUSetS := make(map[int]interface{}) //src,dst + copyIntMap(USetS, &cadUSetS) + cadSetS[H] = nil + cadUSetS[I] = nil + delete(cadSetS, I) + delete(cadUSetS, H) + sumKjih := 0.0 + for D, _ := range cadUSetS { + J := D + if J == I { + + continue + } + //log.Println("for ", D, "in cadUSetS") + d := subDataInt(J, ctx) + minD, min2D := min2minDistance(d, subDataMap(ctx.SetS, ctx)) + distanceIJ := minDistance2Vec(d, subDataInt(I, ctx)) + distanceHJ := minDistance2Vec(d, subDataInt(H, ctx)) + Kjih := 0.0 + //log.Println("distanceIJ minD", distanceIJ, minD) + //log.Println("distanceHJ min2D", distanceHJ, min2D) + if distanceIJ > minD { + if distanceHJ-minD < 0 { + Kjih = distanceHJ - minD + } else { + Kjih = 0 + } + } + if math.Abs(distanceIJ-minD) < 1e-5 { + Kjih = -minD + if distanceHJ < min2D { + Kjih += distanceHJ + } else { + Kjih += min2D + } + } else if distanceIJ < minD { + log.Fatalln("Input must be a string") + } + sumKjih += Kjih + //log.Println(Kjih) + + } + if sumKjih < minKjih { + minKjih = sumKjih + minSetS = cadSetS + minUSetS = cadUSetS + } } } + if minKjih < 0 { + log.Println("minKjih", minKjih) + ctx.SetS = minSetS + USetS = minUSetS + //SetS + } else { + log.Println("minKjih", minKjih) + //SetS + log.Println("clustering is done.") + break + } } } @@ -96,5 +258,8 @@ func main() { ctx.data = initData() ctx.SetS = make(map[int]interface{}) ctx.SetS[1] = nil + ctx.SetS[2] = nil + selectInitCluster(&ctx) + swap(&ctx) fmt.Println("Hello World!") }