code for som helix #49

code/SOM/SOM.jl

@@ -0,0 +1,50 @@
+function classify(x::Array{Float64,1}, prototype::Array{Float64, 2})
+    dist = sum((prototype .- x').^2, dims = 2)
+    return findmin(dist[:])
+end
+
+function calc_dist(x::Array{Int, 1}, m::Array{Array{Int, 1}, 1})
+    res = zeros(length(m))
+    for i = 1:length(m)
+        res[i] = sum((m[i] - x).^2)
+    end
+    return res
+end
+
+function SOM(data; q1::Int = 5, q2::Int = 5, R = 2, niter = 40)
+    K = q1 * q2
+    # numbe of observations
+    N = size(data, 1)
+    idx = sample(1:N, K, replace=false)
+    coor = data[idx, :]
+    coor_grid = [ [i...] for i in Iterators.product(1:q1, 1:q2)][:]
+
+    total_niter = niter * N
+    for iter = 1:niter
+        for i = 1:N
+            iter_i = N * (iter - 1) + i
+            α = -1 / total_niter * iter_i + 1
+            r = -R / total_niter * iter_i + R
+
+            xi = data[i, :]
+            mj_val, mj_idx = classify(xi, coor)
+            mj = coor_grid[mj_idx]
+            # distance in Q1xQ2
+            mk_idx = findall( calc_dist(mj, coor_grid) .<= r )
+            mk = coor_grid[mk_idx]
+            # distance in R^p
+            coor[mk_idx, :] += α * (xi .- coor[mk_idx, :]')'
+        end
+    end
+    return coor
+end
+
+# data from ../PCA/principal_curve.jl
+coor = SOM(hcat(X1, X2, X3))
+p1 = scatter(coor[:,1], coor[:, 2], coor[:,3], legend=false, title = "q1 = q2 = 5")
+
+coor = SOM(hcat(X1, X2, X3), q1 = 10, q2 = 10)
+p2 = scatter(coor[:,1], coor[:, 2], coor[:,3], legend=false, title = "q1 = q2 = 10")
+
+plot(p1, p2)
+savefig("som_helix.png")