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%20Smallest%20sphere%20enclosing%20a%20set%20of%20points.%0A%20%20%20%20%3D%3D%3D%3D%3D%3D%3D%3D%3D%3D%3D%3D%3D%3D%3D%3D%3D%3D%3D%3D%3D%3D%3D%3D%3D%3D%3D%0A%0A%20%20%20%20The%20aim%20of%20this%20tutorial%20is%20two-fold%0A%0A%20%20%20%201.%20Demostrate%20how%20to%20write%20a%20conic%20quadratic%20model%20in%20%60Fusion%60%20in%20a%20very%20simple%20and%20compact%20way.%0A%20%20%20%202.%20Show%20how%20and%20way%20the%20dual%20formulation%20may%20solved%20more%20efficiently.%0A%0A%0A%20%20%20%20Our%20problem%20is%20the%20defined%20as%3A%0A%0A%20%20%20%20**Find%20the%20smallest%20sphere%20that%20encloses%20a%20set%20of**%20%24k%24%20**points**%20%24p_i%20%5Cin%20%5Cmathbb%7BR%7D%5En%24.%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%20%23%20'%25matplotlib%20inline'%20command%20supported%20automatically%20in%20marimo%0A%20%20%20%20import%20matplotlib.pyplot%20as%20plt%0A%20%20%20%20import%20matplotlib.patches%20as%20mpatches%0A%20%20%20%20import%20random%0A%0A%20%20%20%20def%20plot_points(p%2C%20p0%3D%5B%5D%2C%20r0%3D0.)%3A%0A%20%20%20%20%20%20%20%20n%2Ck%3D%20len(p0)%2C%20len(p)%0A%0A%20%20%20%20%20%20%20%20plt.rc('savefig'%2Cdpi%3D120)%0A%0A%20%20%20%20%20%20%20%20fig%2C%20ax%20%3D%20plt.subplots()%0A%20%20%20%20%20%20%20%20ax.set_aspect('equal')%0A%20%20%20%20%20%20%20%20ax.plot(%5B%20p%5Bi%5D%5B0%5D%20for%20i%20in%20range(k)%5D%2C%20%5B%20p%5Bi%5D%5B1%5D%20for%20i%20in%20range(k)%5D%2C%20'b*')%0A%0A%20%20%20%20%20%20%20%20if%20len(p0)%3E0%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20ax.plot(%20%20p0%5B0%5D%2Cp0%5B1%5D%2C%20'r.')%0A%20%20%20%20%20%20%20%20%20%20%20%20ax.add_patch(%20mpatches.Circle(%20p0%2C%20%20r0%20%2C%20%20fc%3D%22w%22%2C%20ec%3D%22r%22%2C%20lw%3D1.5)%20)%0A%20%20%20%20%20%20%20%20plt.grid()%0A%20%20%20%20%20%20%20%20plt.show()%0A%0A%20%20%20%20n%20%3D%202%0A%20%20%20%20k%20%3D%20500%0A%0A%20%20%20%20p%3D%20%20%5B%20%5Brandom.gauss(0.%2C10.)%20for%20nn%20in%20range(n)%5D%20for%20kk%20in%20range(k)%5D%0A%0A%20%20%20%20plot_points(p)%0A%20%20%20%20return%20p%2C%20plot_points%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%20The%20problem%20boils%20down%20to%20determine%20the%20sphere%20center%20%24p_0%5Cin%20%5Cmathbb%7BR%7D%5En%24%20and%20its%20radius%20%24r_0%5Cgeq%200%24%2C%20i.e.%0A%0A%0A%20%20%20%20%5Cbegin%7Bequation%7D%0A%20%20%20%20%20%20%5Cbegin%7Baligned%7D%0A%20%20%20%20%5Cmin%20%5Cmax_%7Bi%3D1%2C%5Cdots%2Ck%7D%20%5C%7C%20p_0%20-%20p_i%5C%7C_2%20%5C%5C%0A%20%20%20%20%20%20%5Cend%7Baligned%7D%0A%20%20%20%20%5Cend%7Bequation%7D%0A%0A%20%20%20%20The%20maximum%20in%20the%20objective%20function%20can%20be%20easily%2C%20i.e.%0A%0A%20%20%20%20%5Cbegin%7Bequation%7D%0A%20%20%20%20%20%20%5Cbegin%7Baligned%7D%0A%20%20%20%20%5Cmin%20r_0%20%26%20%26%20%26%5C%5C%0A%20%20%20%20s.t.%26%20%26%20%26%5C%5C%0A%20%20%20%20%26%20r_0%20%5Cgeq%20%5C%7C%20p_0%20-%20p_i%5C%7C_2%20%2C%26%20%5Cquad%20%26%20i%3D1%2C%5Cldots%2Ck%5C%5C%0A%20%20%20%20%5Cend%7Baligned%7D%0A%20%20%20%20%5Cend%7Bequation%7D%0A%0A%20%20%20%20The%20SOCP%20formulation%20reads%0A%0A%20%20%20%20%5Cbegin%7Bequation%7D%0A%20%20%20%20%20%20%5Cbegin%7Baligned%7D%0A%20%20%20%20%5Cmin%20r_0%20%26%20%26%20%26%5C%5C%0A%20%20%20%20s.t.%26%20%26%20%26%5C%5C%0A%20%20%20%20%26%20%5Cleft%5Br_0%2Cp_0%20-%20p_i%5Cright%5D%20%5Cin%20Q%5E%7B(n%2B1)%7D%2C%20%26%20%5Cquad%20%26%20i%3D1%2C%5Cldots%2Ck.%0A%20%20%20%20%5Cend%7Baligned%7D%0A%20%20%20%20%5Cend%7Bequation%7D%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(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20Before%20defining%20the%20constraints%2C%20we%20note%20that%20we%20can%20write%0A%0A%0A%20%20%20%20%5Cbegin%7Bequation%7D%0A%20%20%20%20R_0%20%3D%20%5Cleft(%20%5Cbegin%7Barray%7D%7Bc%7D%20r_0%20%20%20%5C%5C%20%5Cvdots%20%5C%5C%20r_0%20%20%20%5Cend%7Barray%7D%20%5Cright)%20%5Cin%20%5Cmathbb%7BR%7D%5Ek%20%20%20%20%20%20%20%20%20%20%2C%20%5Cquad%0A%20%20%20%20P_0%20%3D%20%5Cleft(%20%5Cbegin%7Barray%7D%7Bc%7D%20p_0%5ET%20%5C%5C%20%5Cvdots%20%5C%5C%20p_0%5ET%20%5Cend%7Barray%7D%20%5Cright)%20%5Cin%20%5Cmathbb%7BR%7D%5E%7Bk%5Ctimes%20n%7D%2C%20%5Cquad%0A%20%20%20%20P%20%20%20%3D%20%5Cleft(%20%5Cbegin%7Barray%7D%7Bc%7D%20p_1%5ET%20%5C%5C%20%5Cvdots%20%5C%5C%20p_k%5ET%20%5Cend%7Barray%7D%20%5Cright)%20%5Cin%20%5Cmathbb%7BR%7D%5E%7Bk%5Ctimes%20n%7D.%0A%20%20%20%20%5Cend%7Bequation%7D%0A%0A%20%20%20%20so%20that%20%0A%0A%20%20%20%20%5Cbegin%7Bequation%7D%0A%20%20%20%20%5Cleft%5Br_0%2Cp_i%20-%20p_0%5Cright%5D%20%5Cin%20Q%5E%7B(n%2B1)%7D%2C%20%20%5Cquad%20%20i%3D1%2C%5Cldots%2Ck.%0A%20%20%20%20%5Cend%7Bequation%7D%0A%0A%20%20%20%20can%20be%20compactly%20expressed%20as%20%0A%0A%20%20%20%20%5Cbegin%7Bequation%7D%0A%20%20%20%20%5Cleft%5B%20R_0%2CP_0-P%5Cright%5D%20%5Cin%20%5CPi%20Q%5E%7B(n%2B1)%7D%2C%0A%20%20%20%20%5Cend%7Bequation%7D%0A%0A%20%20%20%20that%20means%2C%20with%20a%20little%20abuse%20of%20notation%2C%20that%20each%20rows%20belongs%20to%20a%20quadratic%20cone%20of%20dimension%20%24n%2B1%24.%0A%0A%0A%20%20%20%20Now%20we%20are%20ready%20for%20a%20compact%20implementation%20in%20%60Fusion%60%3A%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%20from%20mosek.fusion%20import%20Model%2C%20Domain%2C%20Expr%2C%20ObjectiveSense%2C%20Matrix%2C%20Var%2C%20SolutionStatus%0A%20%20%20%20import%20mosek.fusion.pythonic%0A%20%20%20%20import%20mosek%20as%20msk%0A%20%20%20%20import%20sys%0A%0A%0A%20%20%20%20def%20primal_problem(P)%3A%0A%0A%20%20%20%20%20%20%20%20print(Model.getVersion())%0A%0A%20%20%20%20%20%20%20%20k%3D%20len(P)%0A%20%20%20%20%20%20%20%20if%20k%3D%3D0%3A%20return%20-1%2C%5B%5D%0A%0A%20%20%20%20%20%20%20%20n%3D%20len(P%5B0%5D)%0A%0A%20%20%20%20%20%20%20%20with%20Model(%22minimal%20sphere%20enclosing%20a%20set%20of%20points%20-%20primal%22)%20as%20M%3A%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20r0%20%3D%20M.variable(1%20%20%20%20%2C%20Domain.greaterThan(0.))%0A%20%20%20%20%20%20%20%20%20%20%20%20p0%20%3D%20M.variable(%5B1%2Cn%5D%2C%20Domain.unbounded())%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20R0%20%3D%20Var.repeat(r0%2Ck)%0A%20%20%20%20%20%20%20%20%20%20%20%20P0%20%3D%20Var.repeat(p0%2Ck)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20M.constraint(%20Expr.hstack(%20R0%2C%20P0%20-%20P%20)%2C%20Domain.inQCone())%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20M.objective(ObjectiveSense.Minimize%2C%20r0)%0A%20%20%20%20%20%20%20%20%20%20%20%20M.setLogHandler(sys.stdout)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20M.solve()%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20r0.level()%5B0%5D%2C%20p0.level()%0A%20%20%20%20return%20(%0A%20%20%20%20%20%20%20%20Domain%2C%0A%20%20%20%20%20%20%20%20Expr%2C%0A%20%20%20%20%20%20%20%20Matrix%2C%0A%20%20%20%20%20%20%20%20Model%2C%0A%20%20%20%20%20%20%20%20ObjectiveSense%2C%0A%20%20%20%20%20%20%20%20Var%2C%0A%20%20%20%20%20%20%20%20primal_problem%2C%0A%20%20%20%20%20%20%20%20sys%2C%0A%20%20%20%20)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22We%20will%20also%20store%20the%20solver%20output%20in%20a%20file%20to%20use%20it%20later%20on.%20And%20then%20just%20solve%20the%20problem.%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(p%2C%20primal_problem)%3A%0A%20%20%20%20r0%2Cp0%20%3D%20primal_problem(p)%0A%0A%20%20%20%20print%20(%22r0%5E*%20%3D%20%22%2C%20r0)%0A%20%20%20%20print%20(%22p0%5E*%20%3D%20%22%2C%20p0)%0A%20%20%20%20return%20p0%2C%20r0%0A%0A%0A%40app.cell%0Adef%20_(p%2C%20p0%2C%20plot_points%2C%20r0)%3A%0A%20%20%20%20plot_points(p%2Cp0%2Cr0)%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%20Dual%20Formulation%20%0A%20%20%20%20-----------------------%0A%0A%20%20%20%20The%20dual%20problem%20can%20be%20determined%20in%20few%20steps%20following%20basic%20rules.%20Introducing%20dual%20variables%0A%0A%20%20%20%20%5Cbegin%7Bequation%7D%0A%20%20%20%20%20Y%20%3D%20%5Cleft(%20%5Cbegin%7Barray%7D%7Bc%7D%20y_1%5ET%5C%5C%20%5Cvdots%20%5C%5Cy_k%20%20%5Cend%7Barray%7D%5Cright)%2C%20%5Cquad%20z%20%3D%20%5Cleft(%20%5Cbegin%7Barray%7D%7Bc%7D%20z_1%5C%5C%20%5Cvdots%20%5C%5Cz_k%20%20%5Cend%7Barray%7D%5Cright)%2C%20%0A%20%20%20%20%5Cend%7Bequation%7D%0A%0A%20%20%20%20the%20dual%20is%3A%0A%0A%20%20%20%20%5Cbegin%7Baligned%7D%0A%20%20%20%20%20%20%20%20%5Cmax%20%26%20%5Cleft%5Clangle%20P%2C%20Y%20%5Cright%5Crangle%20%5C%5C%0A%20%20%20%20%20%20%20%20%26%20e_k%5ET%20z%20%3D%201%5C%5C%0A%20%20%20%20%20%20%20%20%26%20Y%5ET%20e_k%20%3D%20%5Cmathbf%7B0%7D_n%20%5C%5C%0A%20%20%20%20%20%20%20%20%26%20%5Cleft%5Bz_i%20%2C%20y_i%5Cright%5D%20%5Cin%20%5Cmathcal%7BQ%7D%5E%7Bn%2B1%7D%5C%5C%0A%20%20%20%20%20%20%20%20%26%20z_i%5Cin%20%5Cmathbb%7BR%7D%2C%20y_i%5Cin%20%5Cmathbb%7BR%7D%5En%2C%0A%20%20%20%20%5Cend%7Baligned%7D%0A%0A%20%20%20%20where%20%24e_k%5Cin%20%5Cmathbb%7BR%7D%5Ek%24%20is%20a%20vector%20of%20all%20ones.%0A%0A%20%20%20%20The%20%60%60Fusion%60%60%20code%20is%20the%20following%3A%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_(Domain%2C%20Expr%2C%20Matrix%2C%20Model%2C%20ObjectiveSense%2C%20Var%2C%20p%2C%20sys)%3A%0A%20%20%20%20def%20dual_problem(P)%3A%0A%0A%20%20%20%20%20%20%20%20k%3D%20len(P)%0A%20%20%20%20%20%20%20%20if%20k%3D%3D0%3A%20return%20-1%2C%5B%5D%0A%0A%20%20%20%20%20%20%20%20n%3D%20len(P%5B0%5D)%0A%0A%20%20%20%20%20%20%20%20with%20Model(%22minimal%20sphere%20enclosing%20a%20set%20of%20points%20-%20dual%22)%20as%20M%3A%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20Y%3D%20M.variable(%5Bk%2Cn%5D%2C%20Domain.unbounded())%0A%20%20%20%20%20%20%20%20%20%20%20%20z%3D%20M.variable(k%20%20%20%20%2C%20Domain.unbounded())%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20M.constraint(Expr.sum(z)%20%3D%3D%201.0)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20e%3D%20%5B1.0%20for%20i%20in%20range(k)%5D%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20M.constraint(Y.T%20%40%20Matrix.ones(k%2C1)%20%3D%3D%200)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20M.constraint(%20Var.hstack(z%2CY)%2C%20Domain.inQCone())%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20M.objective(%20ObjectiveSense.Maximize%2C%20Expr.dot(%20P%2C%20Y%20))%20%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20M.setLogHandler(sys.stdout)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20M.solve()%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20%0A%0A%20%20%20%20dual_problem(p)%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%20We%20can%20track%20and%20compare%20the%20progression%20of%20both%20the%20**primal**%20and%20**dual**%20problems%20from%20the%20log%20outputs.%20When%20we%20take%20a%20close%20look%20into%20the%20**primal%20objective%20(POBJ)**%20and%20**dual%20objective%20(DOBJ)**%20progression%2C%20we%20see%20that%3A%20%20%0A%0A%20%20%20%20**The%20solver%20has%20got%20the%20same%20result%20doing%20exactly%20the%20same%20steps!**%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(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20In%20fact%2C%20**MOSEK**%20solves%20the%20same%20problem%2C%20thanks%20to%20the%20automatic%20dualizer%20that%20decides%20whether%20it%20is%20more%20convenient%20to%20solve%20the%20**primal**%20or%20the%20**dual**%20of%20a%20given%20problem.%20%20%0A%0A%20%20%20%20The%20primal%20problem%20log%20output%20reports%20the%20solved%20problem%20as%20*dual*%3A%20%20%0A%20%20%20%20%60Optimizer%20%20-%20solved%20problem%20%20%20%20%20%20%20%20%20%3A%20the%20dual%60%0A%0A%20%20%20%20While%20the%20dual%20formulation%20reports%20solving%20the%20*primal*%3A%20%20%0A%20%20%20%20%60Optimizer%20%20-%20solved%20problem%20%20%20%20%20%20%20%20%20%3A%20the%20primal%60%0A%0A%20%20%20%20This%20proves%20the%20claim.%20Therefore%2C%20in%20both%20cases%2C%20we%20actually%20solve%20the%20**very%20same%20formulation**.%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%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
86ea595066ebab6e207a1bf80bba837a2641cc25486620cce961a067885d57f0