forked from jump-dev/JuMP.jl
-
Notifications
You must be signed in to change notification settings - Fork 0
/
cluster.jl
85 lines (71 loc) · 1.79 KB
/
cluster.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
##############################################################################
# JuMP
# An algebraic modeling langauge for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
# cluster.jl
#
# From "Approximating K-means-type clustering via semidefinite programming"
# By Jiming Peng and Yu Wei
#
# Given a set of points a_1, ..., a_m in R_n, allocate them to k clusters
#############################################################################
using JuMP
using SCS
solver = SCSSolver()
# Data points
n = 2
m = 6
a = Any[[2.0, 2.0], [2.5, 2.1], [7.0, 7.0],
[2.2, 2.3], [6.8, 7.0], [7.2, 7.5]]
k = 2
# Weight matrix
W = zeros(m,m)
for i = 1:m
for j = i+1:m
dist = exp(-norm(a[i] - a[j])/1.0)
W[i,j] = dist
W[j,i] = dist
end
end
mod = Model(solver=SCSSolver())
# Z >= 0, PSD
@variable(mod, Z[1:m,1:m], SDP)
@constraint(mod, Z .>= 0)
# min Tr(W(I-Z))
@objective(mod, Min, trace(W * (eye(m,m) - Z)))
# Z e = e
@constraint(mod, Z*ones(m) .== ones(m))
#for i = 1:m
# addConstraint(mod, sum([ Z[i,j] for j = 1:m]) == 1)
#end
# Tr(Z) = k
@constraint(mod, trace(Z) == k)
solve(mod)
Z_val = getvalue(Z)[:,:]
println("Raw solution")
println(round(Z_val,4))
# A simple rounding scheme
which_cluster = zeros(Int,m)
num_clusters = 0
for i = 1:m
Z_val[i,i] <= 1e-6 && continue
if which_cluster[i] == 0
num_clusters += 1
which_cluster[i] = num_clusters
for j = i+1:m
if norm(Z_val[i,j] - Z_val[i,i]) <= 1e-6
which_cluster[j] = num_clusters
end
end
end
end
# Print results
for cluster = 1:k
println("Cluster $cluster")
for i = 1:m
if which_cluster[i] == cluster
println(a[i])
end
end
end