经过无数个bug的消灭，收敛了。

Signed-off-by: facat@lab.com <facat@lab.com>

.gitignore

@@ -0,0 +1 @@
+*.exe

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;

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!")
 }