import%20marimo%0A%0A__generated_with%20%3D%20%220.24.0%22%0Aapp%20%3D%20marimo.App(width%3D%22medium%22)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%20Constrained%20%24k%24-means%20clustering%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20import%20warnings%0A%20%20%20%20from%20pathlib%20import%20Path%0A%0A%20%20%20%20warnings.filterwarnings(%22ignore%22)%0A%0A%20%20%20%20import%20cvxpy%20as%20cp%0A%20%20%20%20import%20marimo%20as%20mo%0A%20%20%20%20import%20matplotlib.pyplot%20as%20plt%0A%20%20%20%20import%20numpy%20as%20np%0A%0A%20%20%20%20from%20dbcp%20import%20BiconvexProblem%0A%0A%20%20%20%20plt.style.use(Path(__file__).resolve().parent%20%2F%20%22zhlatex.mplstyle%22)%0A%20%20%20%20figure_directory%20%3D%20Path(__file__).resolve().parent%20%2F%20%22figures%22%0A%20%20%20%20figure_directory.mkdir(parents%3DTrue%2C%20exist_ok%3DTrue)%0A%0A%20%20%20%20np.random.seed(10015)%0A%20%20%20%20return%20BiconvexProblem%2C%20cp%2C%20figure_directory%2C%20mo%2C%20np%2C%20plt%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Introduction%0A%0A%20%20%20%20We%20consider%20the%20problem%20of%20%24k%24-means%20clustering%20with%20constraints%3A%0A%0A%20%20%20%20%5C%5B%0A%20%20%20%20%20%20%20%20%5Cbegin%7Barray%7D%7Bll%7D%0A%20%20%20%20%20%20%20%20%20%20%20%20%5Ctext%7Bminimize%7D%20%26%20%5Csum_%7Bi%20%3D%201%7D%5E%7Bm%7D%20z_i%5ET%20(%7B%5C%7C%5Cbar%7Bx%7D_1%20-%20x_i%5C%7C%7D_2%5E2%2C%20%5Cldots%2C%20%7B%5C%7C%5Cbar%7Bx%7D_k%20-%20x_i%5C%7C%7D_2%5E2)%5C%5C%0A%20%20%20%20%20%20%20%20%20%20%20%20%5Ctext%7Bsubject%20to%7D%20%26%200%20%5Cpreceq%20z_i%20%5Cpreceq%20%5Cmathbf%7B1%7D%2C%5Cquad%20%5Cmathbf%7B1%7D%5ET%20z_i%20%3D%201%2C%5Cquad%20i%20%3D%201%2C%20%5Cldots%2C%20m%5C%5C%0A%20%20%20%20%20%20%20%20%20%20%20%20%26%20%7B%5C%7C%5Cbar%7Bx%7D_i%20-%20%5Cmu_i%5C%7C%7D_2%20%5Cleq%20r_i%0A%20%20%20%20%20%20%20%20%5Cend%7Barray%7D%0A%20%20%20%20%5C%5D%0A%0A%20%20%20%20with%20variables%20%24%5Cbar%7Bx%7D_i%20%5Cin%20%5Cmathbf%7BR%7D%5En%24%2C%20%24i%20%3D%201%2C%20%5Cldots%2C%20k%24%2C%20and%20%24z_i%20%5Cin%20%5Cmathbf%7BR%7D%5Ek%24%2C%20%24i%20%3D%201%2C%20%5Cldots%2C%20m%24.%0A%20%20%20%20The%20constraints%20%24%7B%5C%7C%5Cbar%7Bx%7D_i%20-%20%5Cmu_i%5C%7C%7D_2%20%5Cleq%20r_i%24%20can%20be%20interpreted%20as%20limiting%20the%20cluster%20center%0A%20%20%20%20%24%5Cbar%7Bx%7D_i%24%20to%20be%20within%20the%20Euclidean%20ball%20with%20(given)%20center%20%24%5Cmu_i%24%20and%20radius%20%24r_i%24.%0A%20%20%20%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(r%22%22%22%0A%20%20%20%20%23%23%20Generate%20dataset%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np)%3A%0A%20%20%20%20n%20%3D%202%0A%20%20%20%20k%20%3D%203%0A%20%20%20%20_m_cluster%20%3D%20%5B180%2C%20180%2C%20180%5D%0A%20%20%20%20_m_outlier%20%3D%2024%0A%20%20%20%20m%20%3D%20np.sum(_m_cluster)%20%2B%20k%20*%20_m_outlier%0A%0A%20%20%20%20mu_true%20%3D%20np.array(%0A%20%20%20%20%20%20%20%20%5B%0A%20%20%20%20%20%20%20%20%20%20%20%20%5B0.25%2C%20-0.15%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%5B4.65%2C%200.20%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%5B2.60%2C%204.30%5D%2C%0A%20%20%20%20%20%20%20%20%5D%0A%20%20%20%20)%0A%20%20%20%20mus%20%3D%20mu_true%20%2B%20np.random.normal(0%2C%201%2C%20size%3Dmu_true.shape)%0A%20%20%20%20r%20%3D%200.5%0A%0A%20%20%20%20%23%20core%20data%20points%0A%20%20%20%20_covs%20%3D%20%5B%0A%20%20%20%20%20%20%20%20np.array(%5B%5B0.45%2C%200.10%5D%2C%20%5B0.10%2C%200.30%5D%5D)%2C%0A%20%20%20%20%20%20%20%20np.array(%5B%5B0.35%2C%20-0.08%5D%2C%20%5B-0.08%2C%200.50%5D%5D)%2C%0A%20%20%20%20%20%20%20%20np.array(%5B%5B0.40%2C%200.00%5D%2C%20%5B0.00%2C%200.40%5D%5D)%2C%0A%20%20%20%20%5D%0A%0A%20%20%20%20_xs%20%3D%20%5B%5D%0A%20%20%20%20for%20_p%20in%20range(k)%3A%0A%20%20%20%20%20%20%20%20_xs.append(np.random.multivariate_normal(mu_true%5B_p%5D%2C%20_covs%5B_p%5D%2C%20size%3D_m_cluster%5B_p%5D))%0A%20%20%20%20x_core%20%3D%20np.vstack(_xs)%0A%0A%20%20%20%20%23%20outlier%20points%0A%20%20%20%20_pull%20%3D%20np.array(%0A%20%20%20%20%20%20%20%20%5B%0A%20%20%20%20%20%20%20%20%20%20%20%20%5B3.0%2C%20-2.0%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%5B-2.5%2C%202.0%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%5B1.5%2C%20-3.0%5D%2C%0A%20%20%20%20%20%20%20%20%5D%0A%20%20%20%20)%0A%0A%20%20%20%20_xs%20%3D%20%5B%5D%0A%20%20%20%20for%20_p%20in%20range(k)%3A%0A%20%20%20%20%20%20%20%20_base%20%3D%20mu_true%5B_p%5D%20%2B%20_pull%5B_p%5D%0A%20%20%20%20%20%20%20%20_xs.append(_base%20%2B%20np.random.normal(0%2C%200.35%2C%20size%3D(_m_outlier%2C%20n)))%0A%0A%20%20%20%20x_out%20%3D%20np.vstack(_xs)%0A%0A%20%20%20%20%23%20combine%0A%20%20%20%20xs%20%3D%20np.vstack(%5Bx_core%2C%20x_out%5D)%0A%20%20%20%20return%20k%2C%20m%2C%20mus%2C%20n%2C%20r%2C%20xs%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Specify%20the%20problem%20and%20solve%0A%20%20%20%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(r%22%22%22%0A%20%20%20%20%23%23%23%20Standard%20%24k%24-means%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(BiconvexProblem%2C%20cp%2C%20k%2C%20m%2C%20n%2C%20xs)%3A%0A%20%20%20%20xbars%20%3D%20cp.Variable((k%2C%20n))%0A%20%20%20%20zs%20%3D%20cp.Variable((m%2C%20k)%2C%20nonneg%3DTrue)%0A%20%20%20%20_d%20%3D%20cp.sum_squares(xs%5B%3A%2C%20None%2C%20%3A%5D%20-%20xbars%5BNone%2C%20%3A%2C%20%3A%5D%2C%20axis%3D2)%0A%20%20%20%20_obj%20%3D%20cp.Minimize(cp.sum(cp.multiply(zs%2C%20_d)))%0A%20%20%20%20_constr%20%3D%20%5Bzs%20%3C%3D%201%2C%20cp.sum(zs%2C%20axis%3D1)%20%3D%3D%201%5D%0A%20%20%20%20_prob%20%3D%20BiconvexProblem(_obj%2C%20%5Bxbars%5D%2C%20%5Bzs%5D%2C%20_constr)%0A%20%20%20%20_prob.solve(cp.CLARABEL%2C%20lbd%3D2)%0A%20%20%20%20return%20(xbars%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%20%24k%24-means%20with%20constraints%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(BiconvexProblem%2C%20cp%2C%20k%2C%20m%2C%20mus%2C%20n%2C%20r%2C%20xs)%3A%0A%20%20%20%20xbars_constr%20%3D%20cp.Variable((k%2C%20n))%0A%20%20%20%20zs_constr%20%3D%20cp.Variable((m%2C%20k)%2C%20nonneg%3DTrue)%0A%20%20%20%20_d%20%3D%20cp.sum_squares(xs%5B%3A%2C%20None%2C%20%3A%5D%20-%20xbars_constr%5BNone%2C%20%3A%2C%20%3A%5D%2C%20axis%3D2)%0A%20%20%20%20_obj%20%3D%20cp.Minimize(cp.sum(cp.multiply(zs_constr%2C%20_d)))%0A%20%20%20%20_constr%20%3D%20%5Bzs_constr%20%3C%3D%201%2C%20cp.sum(zs_constr%2C%20axis%3D1)%20%3D%3D%201%5D%0A%20%20%20%20for%20_c%2C%20_mu%20in%20zip(xbars_constr%2C%20mus)%3A%0A%20%20%20%20%20%20%20%20_constr.append(cp.norm2(_c%20-%20_mu)%20%3C%3D%20r)%0A%20%20%20%20_prob%20%3D%20BiconvexProblem(_obj%2C%20%5Bxbars_constr%5D%2C%20%5Bzs_constr%5D%2C%20_constr)%0A%20%20%20%20_prob.solve(cp.CLARABEL%2C%20lbd%3D10%2C%20abs_tol%3D1e-1)%0A%20%20%20%20return%20(xbars_constr%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Plot%20the%20results%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(figure_directory%2C%20mus%2C%20plt%2C%20xbars%2C%20xbars_constr%2C%20xs)%3A%0A%20%20%20%20fig%2C%20axs%20%3D%20plt.subplots(1%2C%201%2C%20figsize%3D(5%2C%205))%0A%20%20%20%20axs.scatter(xs%5B%3A%2C%200%5D%2C%20xs%5B%3A%2C%201%5D%2C%20s%3D20%2C%20alpha%3D0.5)%0A%20%20%20%20axs.scatter(mus%5B%3A%2C%200%5D%2C%20mus%5B%3A%2C%201%5D%2C%20s%3D100%2C%20color%3D%22r%22%2C%20marker%3D%22%5E%22%2C%20label%3Dr%22%24%5Cmu_i%24%22)%0A%20%20%20%20axs.scatter(%0A%20%20%20%20%20%20%20%20xbars_constr.value%5B%3A%2C%200%5D%2C%0A%20%20%20%20%20%20%20%20xbars_constr.value%5B%3A%2C%201%5D%2C%0A%20%20%20%20%20%20%20%20s%3D100%2C%0A%20%20%20%20%20%20%20%20color%3D%22g%22%2C%0A%20%20%20%20%20%20%20%20marker%3D%22*%22%2C%0A%20%20%20%20%20%20%20%20edgecolors%3D%22k%22%2C%0A%20%20%20%20%20%20%20%20linewidth%3D0.5%2C%0A%20%20%20%20%20%20%20%20label%3Dr%22%24%5Cbar%7Bx%7D_i%5E%7B%5Crm%20con%7D%24%22%2C%0A%20%20%20%20)%0A%20%20%20%20axs.scatter(xbars.value%5B%3A%2C%200%5D%2C%20xbars.value%5B%3A%2C%201%5D%2C%20s%3D50%2C%20color%3D%22k%22%2C%20marker%3D%22x%22%2C%20label%3Dr%22%24%5Cbar%7Bx%7D_i%5E%7B%5Crm%20unc%7D%24%22)%0A%20%20%20%20for%20_x%2C%20_y%20in%20mus%3A%0A%20%20%20%20%20%20%20%20axs.add_patch(plt.Circle((_x%2C%20_y)%2C%200.5%2C%20fill%3DFalse%2C%20edgecolor%3D%22k%22%2C%20linewidth%3D1%2C%20linestyle%3D%22dashed%22))%0A%20%20%20%20axs.set_xlabel(%22%24x_1%24%22)%0A%20%20%20%20axs.set_ylabel(%22%24x_2%24%22)%0A%0A%20%20%20%20axs.legend(frameon%3DFalse%2C%20handlelength%3D0.5%2C%20loc%3D%22upper%20left%22)%0A%20%20%20%20fig.tight_layout()%0A%20%20%20%20fig.savefig(figure_directory%20%2F%20%22kmeans_constr.pdf%22%2C%20bbox_inches%3D%22tight%22)%0A%20%20%20%20plt.show()%0A%20%20%20%20return%0A%0A%0Aif%20__name__%20%3D%3D%20%22__main__%22%3A%0A%20%20%20%20app.run()%0A
22fce6d4cb3e6e4374c9a65ff7531d02