kmeans

import%20marimo%0A%0A__generated_with%20%3D%20%220.15.2%22%0Aapp%20%3D%20marimo.App()%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22!%5BMOSEK%20ApS%5D(https%3A%2F%2Fwww.mosek.com%2Fstatic%2Fimages%2Fbranding%2Fwebgraphmoseklogocolor.png%20)%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20%23%20Clustering%20using%20Disjunctive%20Constraints%0A%0A%20%20%20%20%5BK-Means%20clustering%5D(https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FK-means_clustering)%20is%20one%20of%20the%20most%20used%20clustering%20problems%20in%20unsupervised%20learning.%20Typically%2C%20a%20heuristic%20algorithm%20is%20used%20to%20solve%20the%20K-Means%20clustering%20problem.%20Such%20algorithms%20however%20do%20not%20guarantee%20a%20global%20optimum.%20In%20this%20notebook%2C%20we%20show%20how%20K-Means%20can%20be%20expressed%20as%20a%20Generalized%20Disjunctive%20Program%20(GDP)%20with%20a%20Quadratic%20Rotated%20Cone%2C%20and%20we%20demonstrate%20how%20to%20implement%20it%20in%20MOSEK%20using%20Disjunctive%20Constraints%20(DJC).%20Such%20problem%20was%20for%20example%20studied%20by%20%5BPapageorgiou%2C%20Trespalacios%5D(https%3A%2F%2Flink.springer.com%2Farticle%2F10.1007%2Fs13675-017-0088-0)%20and%20%5BKronqvist%2C%20Misener%2C%20Tsay%5D(https%3A%2F%2Flink.springer.com%2Fchapter%2F10.1007%2F978-3-030-78230-6_19).%20We%20further%20show%20a%20modification%20called%20Euclidean%20clustering%20by%20changing%20only%20a%20few%20lines%20of%20code.%20%0A%0A%20%20%20%20We%20assume%20a%20set%20of%20points%20%24%5Ctextbf%7Bp%7D_1%2C%20%5Cldots%2C%20%5Ctextbf%7Bp%7Dn%20%5Cin%20%5Cmathbb%7BR%7D%5E%5Cmathcal%7BD%7D%24%20and%20a%20natural%20number%20%24%5Cmathcal%7BK%7D%20%5Cin%20%5C%7B1%2C%202%2C%20...%2C%20n%5C%7D%24%20which%20specifies%20number%20of%20centroids.%20We%20want%20to%20find%20positions%20of%20the%20centroids%20such%20that%20the%20overall%20squared%20Euclidean%20distance%20from%20each%20point%20to%20the%20closest%20centroid%20is%20minimized.%20The%20formulation%20using%20disjunctions%20can%20look%20as%20follows%0A%0A%20%20%20%20%24%24%5Cbegin%7Barray%7D%7Brll%7D%0A%20%20%20%20%5Ctext%7Bminimize%7D%20%26%20%5Csum_%7Bi%3D1%7D%5En%20d_i%20%26%20%5C%5C%0A%20%20%20%20%5Ctext%7Bsubject%20to%7D%20%26%20%5Cbigvee_%7Bj%20%5Cin%20%5C%7B1%2C%20..%2C%20%5Cmathcal%7BK%7D%5C%7D%7D%20%5CBigl%5B%20d_i%20%5Cgeq%20%7C%7C%20%5Ctextbf%7Bc%7D_j%20-%20%5Ctextbf%7Bp%7D_i%20%7C%7C_2%5E2%20%20%5Cwedge%20Y_i%20%3D%20j%20%5CBigr%5D%2C%20%26%20%5Cforall%20i%20%5Cin%20%5C%7B1%2C%20...%2C%20n%5C%7D%2C%5C%5C%0A%20%20%20%20%26%20c_%7Bj-1%2C%201%7D%20%5Cleq%20c_%7Bj%2C%201%7D%2C%20%26%20%5Cforall%20j%20%5Cin%20%5C%7B2%2C%20..%2C%20%5Cmathcal%7BK%7D%5C%7D%2C%5C%5C%0A%20%20%20%20%5C%5C%0A%20%20%20%20%26%20d_1%2C%20...%2C%20d_n%20%5Cin%20%5Cmathbb%7BR%7D%2C%20%26%20%5C%5C%0A%20%20%20%20%26%20%5Ctextbf%7Bc%7D_1%2C%20...%2C%20%5Ctextbf%7Bc%7D_%7B%5Cmathcal%7BK%7D%7D%20%5Cin%20%5Cmathbb%7BR%7D%5E%5Cmathcal%7BD%7D%2C%20%5C%5C%0A%20%20%20%20%26%20Y_1%2C%20...%2C%20Y_n%20%5Cin%20%5C%7B1%2C%20..%2C%20%5Cmathcal%7BK%7D%5C%7D%2C%20%26%0A%20%20%20%20%5Cend%7Barray%7D%24%24%0A%0A%20%20%20%20where%20%24%5Ctextbf%7Bc%7D_1%2C%20...%2C%20%5Ctextbf%7Bc%7D_%7B%5Cmathcal%7BK%7D%7D%24%20are%20the%20positions%20of%20the%20centroids%2C%20%24d_1%2C%20...%2C%20d_n%24%20are%20auxiliary%20variables%20representing%20the%20shortest%20squared%20distance%20to%20the%20nearest%20centroid%20for%20each%20point%20and%20%24Y_1%2C%20...%2C%20Y_n%24%20are%20classification%20labels%20for%20each%20point%2C%20indicating%20the%20index%20of%20the%20nearest%20centroid.%20The%20first%20constraint%20is%20a%20disjunctive%20constraint%20representing%20the%20choice%20of%20a%20centroid%20for%20each%20point.%20Exactly%20one%20of%20the%20clauses%20in%20each%20of%20%24n%24%20disjunctions%20is%20%22active%22%20and%20determines%20the%20index%20of%20the%20nearest%20centroid%20and%20the%20distance%20to%20is.%20The%20second%20constraint%20in%20the%20formulation%20is%20a%20symmetry-breaking%20constraint.%20%0A%0A%20%20%20%20**MOSEK**%20only%20supports%20Disjunctive%20Normal%20Form%20(DNF)%20of%20affine%20constraints.%20Formally%2C%20this%20means%20that%20each%20Disjunctive%20Constraint%20(DJC)%20is%20of%20the%20form%20%24%5Cbigvee_i%20%5Cbigwedge_j%20T_%7Bi%2C%20j%7D%24%2C%20where%20%24T_%7Bi%2Cj%7D%24%20is%20an%20affine%20constraint.%20Such%20constraint%20is%20satisfied%20if%20and%20only%20if%20there%20exists%20at%20least%20one%20term%20%24i%24%2C%20such%20that%20all%20affine%20constraints%20%24T_%7Bi%2Cj%7D%24%20are%20satisfied.%20We%20therefore%20need%20to%20move%20the%20non-linearity%20out%20of%20the%20disjunction.%20This%20can%20be%20tackled%20by%20a%20using%20new%20auxiliary%20variables%20%24dAux_%7Bi%2C%20j%7D%24%20and%20constraining%20them%20outside%20of%20the%20dicjunctions.%20The%20program%20then%20looks%20in%20the%20following%20way%3A%20%0A%0A%20%20%20%20%24%24%5Cbegin%7Barray%7D%7Brll%7D%0A%20%20%20%20%5Ctext%7Bminimize%7D%20%26%20%5Csum_%7Bi%3D1%7D%5En%20d_i%20%26%20%5C%5C%0A%20%20%20%20%5Ctext%7Bsubject%20to%7D%20%26%20%5Cbigvee_%7Bj%20%5Cin%20%5C%7B1%2C%20..%2C%20%5Cmathcal%7BK%7D%5C%7D%7D%20%5CBigl%5B%20d_i%20%5Cgeq%20dAux_%7Bi%2C%20j%7D%20%5Cwedge%20%5Chspace%7B0.2cm%7D%20Y_i%20%3D%20j%20%5CBigr%5D%2C%20%26%20%5Cforall%20i%20%5Cin%20%5C%7B1%2C%20...%2C%20n%5C%7D%2C%5C%5C%0A%20%20%20%20%26%20dAux_%7Bi%2C%20j%7D%20%5Cgeq%20%7C%7C%20%5Ctextbf%7Bc%7D_j%20-%20%5Ctextbf%7Bp%7D_i%20%7C%7C_2%5E2%2C%20%26%20%5Cforall%20j%20%5Cin%20%5C%7B1%2C%20..%2C%20%5Cmathcal%7BK%7D%5C%7D%2C%20%5Cforall%20i%20%5Cin%20%5C%7B1%2C%20...%2C%20n%5C%7D%20%5C%5C%0A%20%20%20%20%26%20c_%7Bj-1%2C%201%7D%20%5Cleq%20c_%7Bj%2C%201%7D%2C%20%26%20%5Cforall%20j%20%5Cin%20%5C%7B2%2C%20..%2C%20%5Cmathcal%7BK%7D%5C%7D%2C%20%5C%5C%0A%20%20%20%20%5C%5C%0A%20%20%20%20%26%20dAux%20%5Cin%20%5Cmathbb%7BR%7D%5E%7Bn%20%5Ctimes%20%5Cmathcal%7BK%7D%7D%2C%20%26%20%20%5C%5C%0A%20%20%20%20%26%20d_1%2C%20...%2C%20d_n%20%5Cin%20%5Cmathbb%7BR%7D%2C%20%26%20%5C%5C%0A%20%20%20%20%26%20%5Ctextbf%7Bc%7D_1%2C%20...%2C%20%5Ctextbf%7Bc%7D_%7B%5Cmathcal%7BK%7D%7D%20%5Cin%20%5Cmathbb%7BR%7D%5E%5Cmathcal%7BD%7D%2C%20%26%20%5C%5C%0A%20%20%20%20%26%20Y_1%2C%20...%2C%20Y_n%20%5Cin%20%5C%7B1%2C%20..%2C%20%5Cmathcal%7BK%7D%5C%7D.%20%26%0A%20%20%20%20%5Cend%7Barray%7D%24%24%0A%0A%20%20%20%20%23%23%23%20Preparing%20synthetic%20data%0A%0A%20%20%20%20To%20prepare%20the%20synthetic%20data%2C%20we%20generate%203%20clusters%20with%20the%20same%20number%20of%20points.%20These%20clusters%20are%20generated%20randomly%20according%20to%20the%20Multivariate%20normal%20distribution%20each%20with%20an%20appropriate%20mean%20vector%20and%20a%20covariance%20matrix.%0A%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20import%20numpy%20as%20np%0A%20%20%20%20from%20mosek.fusion%20import%20Model%2C%20Domain%2C%20Expr%2C%20ObjectiveSense%2C%20Matrix%2C%20Var%2C%20DJC%0A%20%20%20%20import%20mosek.fusion.pythonic%0A%20%20%20%20import%20matplotlib.pyplot%20as%20plt%0A%20%20%20%20import%20sys%0A%0A%20%20%20%20%23%20make%20the%20randomness%20deterministic%0A%20%20%20%20np.random.seed(0)%0A%20%20%20%20return%20DJC%2C%20Domain%2C%20Expr%2C%20Model%2C%20ObjectiveSense%2C%20np%2C%20plt%2C%20sys%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20plt)%3A%0A%20%20%20%20%23%20specify%20the%20number%20of%20points%0A%20%20%20%20nData%20%3D%2021%0A%0A%20%20%20%20%23%20generate%20data%0A%20%20%20%20numberOfClusters%20%3D%203%0A%20%20%20%20size%20%3D%20nData%20%2F%2F%20numberOfClusters%0A%0A%20%20%20%20pointsClass1%20%3D%20np.random.multivariate_normal(np.array(%5B1%2C%201%5D)*2%2C%20np.array(%5B%5B1%2C%200%5D%2C%20%5B0%2C%201%5D%5D)*0.7%2C%20size%3Dsize)%0A%20%20%20%20pointsClass2%20%3D%20np.random.multivariate_normal(np.array(%5B-1%2C%20-1%5D)*2%2C%20np.array(%5B%5B1%2C%200%5D%2C%20%5B0%2C%201%5D%5D)*0.7%2C%20size%3Dsize)%0A%20%20%20%20pointsClass3%20%3D%20np.random.multivariate_normal(np.array(%5B-1%2C%203%5D)*2%2C%20np.array(%5B%5B1%2C%200%5D%2C%20%5B0%2C%201%5D%5D)*0.7%2C%20size%3Dsize)%0A%0A%20%20%20%20points%20%3D%20np.vstack((pointsClass1%2C%20pointsClass2%2C%20pointsClass3))%0A%0A%20%20%20%20%23%20plot%20the%20generated%20data%0A%0A%20%20%20%20labels%20%3D%20np.zeros(3*size)%0A%20%20%20%20labels%5Bsize%3A2*size%5D%20%3D%201%0A%20%20%20%20labels%5B2*size%3A%5D%20%3D%202%0A%0A%20%20%20%20plt.title(%22Generated%20Data%22)%0A%20%20%20%20plt.scatter(pointsClass1%5B%3A%2C%200%5D%2C%20pointsClass1%5B%3A%2C%201%5D)%0A%20%20%20%20plt.scatter(pointsClass2%5B%3A%2C%200%5D%2C%20pointsClass2%5B%3A%2C%201%5D)%0A%20%20%20%20plt.scatter(pointsClass3%5B%3A%2C%200%5D%2C%20pointsClass3%5B%3A%2C%201%5D)%0A%20%20%20%20plt.show()%0A%20%20%20%20return%20numberOfClusters%2C%20points%2C%20pointsClass1%2C%20pointsClass2%2C%20pointsClass3%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20%23%23%23%20Fusion%20Model%0A%20%20%20%20In%20the%20following%20block%2C%20we%20show%20the%20K-Means%20model%20in%20the%20**Mosek%20Fusion%20for%20Python**%20in%20a%20vectorized%20fashion.%0A%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20DJC%2C%0A%20%20%20%20Domain%2C%0A%20%20%20%20Expr%2C%0A%20%20%20%20Model%2C%0A%20%20%20%20ObjectiveSense%2C%0A%20%20%20%20np%2C%0A%20%20%20%20numberOfClusters%2C%0A%20%20%20%20points%2C%0A%20%20%20%20sys%2C%0A)%3A%0A%20%20%20%20n%2C%20d%20%3D%20points.shape%0A%20%20%20%20xs%20%3D%20np.repeat(points%2C%20numberOfClusters%2C%20axis%3D0)%0A%20%20%20%20with%20Model()%20as%20M%3A%0A%20%20%20%20%20%20%20%20centroids%20%3D%20M.variable(%5BnumberOfClusters%2C%20d%5D%2C%20Domain.unbounded())%0A%20%20%20%20%20%20%20%20distances%20%3D%20M.variable(n%2C%20Domain.greaterThan(0))%0A%20%20%20%20%20%20%20%20distanceAux%20%3D%20M.variable(%5Bn%2C%20numberOfClusters%5D%2C%20Domain.greaterThan(0))%0A%20%20%20%20%20%20%20%20Y%20%3D%20M.variable(n%2C%20Domain.integral(Domain.inRange(0%2C%20numberOfClusters%20-%201)))%0A%20%20%20%20%20%20%20%20a1%20%3D%20Expr.flatten(distanceAux)%0A%20%20%20%20%20%20%20%20a2%20%3D%20Expr.constTerm(%5Bn%20*%20numberOfClusters%2C%201%5D%2C%200.5)%0A%20%20%20%20%20%20%20%20a3%20%3D%20Expr.repeat(centroids%2C%20n%2C%200)%20-%20xs%0A%20%20%20%20%20%20%20%20hstack%20%3D%20Expr.hstack(%5Ba1%2C%20a2%2C%20a3%5D)%0A%20%20%20%20%20%20%20%20M.constraint(hstack%2C%20Domain.inRotatedQCone())%0A%20%20%20%20%20%20%20%20for%20pointInd%20in%20range(0%2C%20n)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20label%20%3D%20Y.index(pointInd)%0A%20%20%20%20%20%20%20%20%20%20%20%20di%20%3D%20distances.index(pointInd)%0A%20%20%20%20%20%20%20%20%20%20%20%20ANDs%20%3D%20%7B%7D%0A%20%20%20%20%20%20%20%20%20%20%20%20for%20clusterInd%20in%20range(0%2C%20numberOfClusters)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20dijAux%20%3D%20distanceAux.index(%5BpointInd%2C%20clusterInd%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20ANDs%5BclusterInd%5D%20%3D%20DJC.AND(DJC.term(di%20%3E%3D%20dijAux)%2C%20DJC.term(label%20%3D%3D%20clusterInd))%0A%20%20%20%20%20%20%20%20%20%20%20%20M.disjunction(%5BANDs%5Bi%5D%20for%20i%20in%20range(0%2C%20numberOfClusters)%5D)%0A%20%20%20%20%20%20%20%20M.constraint(centroids%5B%3A-1%2C%200%5D%20%3C%3D%20centroids%5B1%3A%2C%200%5D)%0A%20%20%20%20%20%20%20%20M.objective(ObjectiveSense.Minimize%2C%20Expr.sum(distances))%0A%20%20%20%20%20%20%20%20M.setLogHandler(sys.stdout)%0A%20%20%20%20%20%20%20%20M.solve()%0A%20%20%20%20%20%20%20%20labels_KMEANS%20%3D%20Y.level()%0A%20%20%20%20%20%20%20%20centroids_KMEANS%20%3D%20np.array(centroids.level()).reshape(numberOfClusters%2C%20d)%0A%20%20%20%20return%20centroids_KMEANS%2C%20labels_KMEANS%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%23%23%23%20Plotting%20the%20results%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(centroids_KMEANS%2C%20labels_KMEANS%2C%20np%2C%20numberOfClusters%2C%20plt%2C%20points)%3A%0A%20%20%20%20def%20plotTheResults(centroids%2C%20points%2C%20labels%2C%20kmeans%3DTrue)%3A%0A%20%20%20%20%20%20%20%20plt.scatter(centroids%5B%3A%2C%200%5D%2C%20centroids%5B%3A%2C%201%5D%2C%20marker%3D'x'%2C%20color%3D'r')%0A%20%20%20%20%20%20%20%20plt.scatter(points%5B%3A%2C%200%5D%2C%20points%5B%3A%2C%201%5D%2C%20marker%3D'o')%0A%20%20%20%20%20%20%20%20for%20i%20in%20range(0%2C%20numberOfClusters)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20centroid%20%3D%20centroids%5Bi%2C%20%3A%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20cluster%20%3D%20points%5Blabels%20%3D%3D%20i%2C%20%3A%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20nPointsInCluster%20%3D%20cluster.shape%5B0%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20centroidCopied%20%3D%20np.tile(centroid%2C%20(nPointsInCluster%2C%201))%0A%20%20%20%20%20%20%20%20%20%20%20%20xx%20%3D%20np.vstack(%5BcentroidCopied%5B%3A%2C%200%5D%2C%20cluster%5B%3A%2C%200%5D%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20yy%20%3D%20np.vstack(%5BcentroidCopied%5B%3A%2C%201%5D%2C%20cluster%5B%3A%2C%201%5D%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20plt.plot(xx%2C%20yy%2C%20'--'%2C%20color%3D'green'%2C%20alpha%3D0.4)%0A%20%20%20%20%20%20%20%20if%20kmeans%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20plt.title('Solution%20to%20the%20K-Means%20clustering%20problem')%0A%20%20%20%20%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20plt.title('Solution%20to%20the%20Euclidean%20clustering%20problem')%0A%20%20%20%20%20%20%20%20plt.legend(%5B'centroids'%2C%20'data'%5D)%0A%20%20%20%20%20%20%20%20plt.show()%0A%20%20%20%20plotTheResults(centroids_KMEANS%2C%20points%2C%20labels_KMEANS%2C%20kmeans%3DTrue)%0A%20%20%20%20return%20(plotTheResults%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20%23%23%20Euclidean%20Clustering%0A%0A%20%20%20%20In%20order%20to%20change%20the%20proposed%20model%20to%20Euclidean%20Clustering%2C%20we%20need%20to%20only%20change%20the%20squared%20norm%20of%20the%20distances%20to%20standard%20Euclidean%20norm.%20This%20approach%20is%20more%20robust%20to%20outliers%20compared%20to%20the%20K-Means%20algorithm.%20Euclidean%20clustering%20problem%20has%20the%20following%20form%3A%0A%0A%20%20%20%20%24%24%5Cbegin%7Barray%7D%7Brll%7D%0A%20%20%20%20%5Ctext%7Bminimize%7D%20%26%20%5Csum_%7Bi%3D1%7D%5En%20d_i%20%26%20%5C%5C%0A%20%20%20%20%5Ctext%7Bsubject%20to%7D%20%26%20%5Cbigvee_%7Bj%20%5Cin%20%5C%7B1%2C%20..%2C%20%5Cmathcal%7BK%7D%5C%7D%7D%20%5CBigl%5B%20%5Chspace%7B0.2cm%7D%20d_i%20%5Cgeq%20dAux_%7Bi%2C%20j%7D%20%5Cwedge%20Y_i%20%3D%20j%20%5CBigr%5D%2C%20%26%20%5Cforall%20i%20%5Cin%20%5C%7B1%2C%20...%2C%20n%5C%7D%2C%5C%5C%0A%20%20%20%20%26%20dAux_%7Bi%2C%20j%7D%20%5Cgeq%20%7C%7C%20%5Ctextbf%7Bc%7D_j%20-%20%5Ctextbf%7Bp%7D_i%20%7C%7C_2%2C%20%26%20%5Cforall%20j%20%5Cin%20%5C%7B1%2C%20..%2C%20%5Cmathcal%7BK%7D%5C%7D%2C%20%26%20%5Cforall%20i%20%5Cin%20%5C%7B1%2C%20...%2C%20n%5C%7D%2C%20%5C%5C%0A%20%20%20%20%26%20c_%7Bj-1%2C%201%7D%20%5Cleq%20c_%7Bj%2C%201%7D%20%2C%20%26%20%5Cforall%20j%20%5Cin%20%5C%7B2%2C%20..%2C%20%5Cmathcal%7BK%7D%5C%7D%2C%20%5C%5C%0A%20%20%20%20%5C%5C%0A%20%20%20%20%26%20dAux%20%5Cin%20%5Cmathbb%7BR%7D%5E%7Bn%20%5Ctimes%20%5Cmathcal%7BK%7D%7D_%2B%20%20%26%20%5C%5C%0A%20%20%20%20%26%20d_1%2C%20...%2C%20d_n%20%5Cin%20%5Cmathbb%7BR%7D_%7B%2B%7D%2C%20%26%20%5C%5C%0A%20%20%20%20%26%20%5Ctextbf%7Bc%7D_1%2C%20...%2C%20%5Ctextbf%7Bc%7D_%7B%5Cmathcal%7BK%7D%7D%20%5Cin%20%5Cmathbb%7BR%7D%5E%5Cmathcal%7BD%7D%2C%20%26%20%5C%5C%0A%20%20%20%20%26%20Y_1%2C%20...%2C%20Y_n%20%5Cin%20%5C%7B1%2C%20..%2C%20%5Cmathcal%7BK%7D%5C%7D.%20%26%0A%20%20%20%20%5Cend%7Barray%7D%24%24%0A%0A%20%20%20%20Heuristic%20algorithms%20usually%20solve%20such%20a%20problem%20by%20finding%20a%20geometric%20median.%20In%20the%20**Fusion%20API**%20model%2C%20we%20need%20to%20only%20swap%20the%20Rotated%20Quadratic%20Cone%20for%20Quadratic%20Cone.%0A%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%23%23%20Data%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(plt%2C%20pointsClass1%2C%20pointsClass2%2C%20pointsClass3)%3A%0A%20%20%20%20%23%20use%20the%20same%20data%0A%0A%20%20%20%20plt.title(%22Generated%20Data%22)%0A%20%20%20%20plt.scatter(pointsClass1%5B%3A%2C%200%5D%2C%20pointsClass1%5B%3A%2C%201%5D)%0A%20%20%20%20plt.scatter(pointsClass2%5B%3A%2C%200%5D%2C%20pointsClass2%5B%3A%2C%201%5D)%0A%20%20%20%20plt.scatter(pointsClass3%5B%3A%2C%200%5D%2C%20pointsClass3%5B%3A%2C%201%5D)%0A%20%20%20%20plt.show()%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%23%23%20Fusion%20Model%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20DJC%2C%0A%20%20%20%20Domain%2C%0A%20%20%20%20Expr%2C%0A%20%20%20%20Model%2C%0A%20%20%20%20ObjectiveSense%2C%0A%20%20%20%20np%2C%0A%20%20%20%20numberOfClusters%2C%0A%20%20%20%20points%2C%0A%20%20%20%20sys%2C%0A)%3A%0A%20%20%20%20def%20_()%3A%0A%20%20%20%20%20%20%20%20n%2C%20d%20%3D%20points.shape%0A%20%20%20%20%20%20%20%20xs%20%3D%20np.repeat(points%2C%20numberOfClusters%2C%20axis%3D0)%0A%20%20%20%20%20%20%20%20with%20Model()%20as%20M%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20centroids%20%3D%20M.variable(%5BnumberOfClusters%2C%20d%5D%2C%20Domain.unbounded())%0A%20%20%20%20%20%20%20%20%20%20%20%20distances%20%3D%20M.variable(n%2C%20Domain.greaterThan(0))%0A%20%20%20%20%20%20%20%20%20%20%20%20distanceAux%20%3D%20M.variable(%5Bn%2C%20numberOfClusters%5D%2C%20Domain.greaterThan(0))%0A%20%20%20%20%20%20%20%20%20%20%20%20Y%20%3D%20M.variable(n%2C%20Domain.integral(Domain.inRange(0%2C%20numberOfClusters%20-%201)))%0A%20%20%20%20%20%20%20%20%20%20%20%20a1%20%3D%20Expr.flatten(distanceAux)%0A%20%20%20%20%20%20%20%20%20%20%20%20a2%20%3D%20Expr.repeat(centroids%2C%20n%2C%200)%20-%20xs%0A%20%20%20%20%20%20%20%20%20%20%20%20hstack%20%3D%20Expr.hstack(%5Ba1%2C%20a2%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20M.constraint(hstack%2C%20Domain.inQCone())%0A%20%20%20%20%20%20%20%20%20%20%20%20for%20pointInd%20in%20range(0%2C%20n)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20label%20%3D%20Y.index(pointInd)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20di%20%3D%20distances.index(pointInd)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20ANDs%20%3D%20%7B%7D%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20for%20clusterInd%20in%20range(0%2C%20numberOfClusters)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20dijAux%20%3D%20distanceAux.index(%5BpointInd%2C%20clusterInd%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20ANDs%5BclusterInd%5D%20%3D%20DJC.AND(DJC.term(di%20%3E%3D%20dijAux)%2C%20DJC.term(label%20%3D%3D%20clusterInd))%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20M.disjunction(%5BANDs%5Bi%5D%20for%20i%20in%20range(0%2C%20numberOfClusters)%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20M.constraint(centroids%5B%3A-1%2C%200%5D%20%3C%3D%20centroids%5B1%3A%2C%200%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20M.objective(ObjectiveSense.Minimize%2C%20Expr.sum(distances))%0A%20%20%20%20%20%20%20%20%20%20%20%20M.setLogHandler(sys.stdout)%0A%20%20%20%20%20%20%20%20%20%20%20%20M.solve()%0A%20%20%20%20%20%20%20%20%20%20%20%20labels_EUCL%20%3D%20Y.level()%0A%20%20%20%20%20%20%20%20%20%20%20%20centroids_EUCL%20%3D%20np.array(centroids.level()).reshape(numberOfClusters%2C%20d)%0A%20%20%20%20%20%20%20%20return%20centroids_EUCL%2C%20labels_EUCL%0A%20%20%20%20centroids_EUCL%2Clabels_EUCL%20%3D%20_()%0A%20%20%20%20return%20centroids_EUCL%2C%20labels_EUCL%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%23%23%23%20Plot%20the%20data%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20centroids_EUCL%2C%0A%20%20%20%20centroids_KMEANS%2C%0A%20%20%20%20labels_EUCL%2C%0A%20%20%20%20labels_KMEANS%2C%0A%20%20%20%20plotTheResults%2C%0A%20%20%20%20points%2C%0A)%3A%0A%20%20%20%20%23%20call%20the%20plotting%20functions%0A%20%20%20%20plotTheResults(centroids_KMEANS%2C%20points%2C%20labels_KMEANS%2C%20kmeans%3DTrue)%0A%20%20%20%20plotTheResults(centroids_EUCL%2C%20points%2C%20labels_EUCL%2C%20kmeans%3DFalse)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%3Ca%20rel%3D%22license%22%20href%3D%22http%3A%2F%2Fcreativecommons.org%2Flicenses%2Fby%2F4.0%2F%22%3E%3Cimg%20alt%3D%22Creative%20Commons%20License%22%20style%3D%22border-width%3A0%22%20src%3D%22https%3A%2F%2Fi.creativecommons.org%2Fl%2Fby%2F4.0%2F80x15.png%22%20%2F%3E%3C%2Fa%3E%3Cbr%20%2F%3EThis%20work%20is%20licensed%20under%20a%20%3Ca%20rel%3D%22license%22%20href%3D%22http%3A%2F%2Fcreativecommons.org%2Flicenses%2Fby%2F4.0%2F%22%3ECreative%20Commons%20Attribution%204.0%20International%20License%3C%2Fa%3E.%20The%20**MOSEK**%20logo%20and%20name%20are%20trademarks%20of%20%3Ca%20href%3D%22http%3A%2F%2Fmosek.com%22%3EMosek%20ApS%3C%2Fa%3E.%20The%20code%20is%20provided%20as-is.%20Compatibility%20with%20future%20release%20of%20**MOSEK**%20or%20the%20%60Fusion%20API%60%20are%20not%20guaranteed.%20For%20more%20information%20contact%20our%20%5Bsupport%5D(mailto%3Asupport%40mosek.com).%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_()%3A%0A%20%20%20%20import%20marimo%20as%20mo%0A%20%20%20%20return%20(mo%2C)%0A%0A%0Aif%20__name__%20%3D%3D%20%22__main__%22%3A%0A%20%20%20%20app.run()%0A