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(r%22%22%22%23%20Regularized%20Wasserstein%20Barycenters%20using%20Mosek%20and%20the%20exponential%20cone%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%22In%20a%20%5Bprevious%20notebook%20related%20to%20Wasserstein%20distances%5D(https%3A%2F%2Fgithub.com%2FMOSEK%2FTutorials%2Ftree%2Fmaster%2Fwasserstein)%20we%20defined%20the%20linear%20optimization%20problem%20of%20computing%20the%20Wasserstein%20barycenter%20of%20a%20set%20of%20discrete%20measures.%20Here%20we%20solve%20an%20entropy-regularized%20variant%20of%20the%20same%20problem%20and%20to%20demonstrate%20the%20exponential%20cone%20capabilities%20of%20MOSEK.%20We%20also%20use%20this%20problem%20to%20compare%20Fusion%20and%20CVXPY.%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%20As%20a%20reminder%2C%20the%20%24p%24-th%20order%20Wasserstein%20distance%20%24W_p(%5Cmu%2C%5Cupsilon)%24%20between%20discrete%20probability%20distributions%20%24%5Cmu%2C%5Cupsilon%24%20is%20the%20objective%20value%20of%20the%20following%20problem%3A%0A%0A%0A%20%20%20%20%24%24%20%5Ctext%7Bminimize%7D%20%5Cquad%20%5Csum_%7Bi%3D1%7D%5Csum_%7Bj%3D1%7D%20D(X_i%2CY_j)%5Ep%5Cpi_%7Bij%7D%24%24%0A%0A%0A%20%20%20%20%24%24%20%5Ctext%7Bst.%7D%20%5Cquad%20%5Csum_%7Bj%3D1%7D%20%5Cpi_%7Bij%7D%20%3D%20%5Cmu_i%20%2C%20%5Cquad%20i%20%3D%201%2C2%2C..n%20%24%24%0A%0A%20%20%20%20%24%24%20%5Cquad%20%5Csum_%7Bi%3D1%7D%20%5Cpi_%7Bij%7D%20%3D%20%5Cupsilon_j%2C%20%5Cquad%20j%20%3D%201%2C2%2C..m%20%24%24%0A%0A%20%20%20%20%24%24%20%5Cpi_%7Bij%7D%20%5Cgeq%200%2C%20%5Cquad%20%5Cforall_%7Bi%2Cj%7D%24%24%0A%0A%20%20%20%20where%20%24D(X_i%2CY_j)%24%20is%20the%20distance%20function.%0A%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%20Wasserstein%20Barycenter%20with%20regularization%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%20The%20entropy%20regularized%20barycenter%20problem%20with%20%24p%3D2%24%20is%3A%0A%0A%20%20%20%20%24%24%20%5Ctext%7Bminimize%7D%20%5Cquad%20%5Cfrac1N%20%5Csum_%7Bi%2Cj%2Ck%7D%5E%7BN%7D%20D(X_i%2CY_j)%5E2%5Cpi_%7Bij%7D%5Ek%20%2B%20%5Cfrac1%5Clambda%5Csum_%7Bi%2Cj%2Ck%7D%20%5Cpi_%7Bij%7D%5Ek%5Clog(%5Cpi_%7Bij%7D%5Ek)%24%24%0A%0A%20%20%20%20%24%24%5Ctext%7Bst.%7D%20%20%5Cquad%20%5Csum_%7Bj%3D1%7D%20%5Cpi_%7Bij%7D%5E%7Bk%7D%20%3D%20%5Cmu_i%2C%20%5Cquad%20%5Cforall_%7Bk%2Ci%7D%20%5Cquad%20(1)%24%24%0A%0A%20%20%20%20%24%24%20%5Cquad%20%5Csum_%7Bi%3D1%7D%20%5Cpi_%7Bij%7D%5E%7Bk%7D%20%3D%20%5Cupsilon_j%5E%7Bk%7D%2C%20%5Cquad%20%5Cforall_%7Bk%2Cj%7D%20%5Cquad%20(2)%20%24%24%0A%0A%20%20%20%20%24%24%20%5Cpi_%7Bij%7D%5E%7Bk%7D%20%5Cgeq%200%20%5Cquad%20%5Cforall_%7Bk%2Ci%2Cj%7D%24%24%0A%0A%20%20%20%20where%20%24D(X_i%2CY_j)%24%20is%20the%20euclidian%20distance%20between%20pixels%2C%20%24%5Clambda%20%3D%20median(D(X_i%2CY_j))%24%20and%20%24N%24%20is%20the%20number%20of%20samples.%0A%0A%20%20%20%20Without%20the%20entropy%20term%20the%20problem%20is%20just%20the%20linear%20problem%20of%20computing%20a%20distribution%20%24%5Cmu%24%20minimizing%20the%20sum%20of%20distances%20to%20%24%5Cupsilon_i%24%2C%20as%20studied%20in%20our%20other%20notebook.%20Entropy%20regularization%20was%20suggested%20to%20us%20by%20Stefano%20Gualandi%20and%20appears%20for%20example%20in%20the%20paper%20by%20Cuturi%20and%20Doucet%20http%3A%2F%2Fproceedings.mlr.press%2Fv32%2Fcuturi14.pdf.%20This%20paper%20contains%20also%20more%20details%20about%20the%20choice%20of%20%24%5Clambda%24.%20Also%20more%20detailed%20information%20about%20LP%20aproach%20to%20Wasserstein%20metric%20can%20be%20found%20in%20%5BStefano%20Gualandi's%20blogpost%5D(http%3A%2F%2Fstegua.github.io%2Fblog%2F2018%2F12%2F31%2Fwasserstein-distances-an-operations-research-perspective%2F).%0A%0A%20%20%20%20In%20this%20problem%2C%20Wasserstein%20Barycenter%20of%20Three's%20are%20visualized%20using%20images%20with%20size%2028x28%20using%20%242%24%20handwriten%20'3'%20digits%20from%20MNIST%20database.%20Computations%20are%20carried%20out%20by%20Intel(R)%20Xeon(R)%20CPU%20E5-2687W%20v4%20%40%203.00GHz%20processor.%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%20struct%0A%20%20%20%20import%20numpy%20as%20np%0A%20%20%20%20import%20pandas%20as%20pd%0A%20%20%20%20%23%20'%25matplotlib%20inline'%20command%20supported%20automatically%20in%20marimo%0A%20%20%20%20import%20matplotlib.pyplot%20as%20plt%0A%0A%20%20%20%20%23Define%20the%20number%20of%20images%20for%20the%20barycenter%20calculation%0A%20%20%20%20n%3D2%0A%20%20%20%20number%20%3D%203%0A%0A%20%20%20%20%23Read%20the%20images%20from%20the%20file%0A%20%20%20%20def%20read_idx(filename)%3A%0A%20%20%20%20%20%20%20%20with%20open(filename%2C%20'rb')%20as%20f%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20zero%2C%20data_type%2C%20dims%20%3D%20struct.unpack('%3EHBB'%2C%20f.read(4))%0A%20%20%20%20%20%20%20%20%20%20%20%20shape%20%3D%20tuple(struct.unpack('%3EI'%2C%20f.read(4))%5B0%5D%20for%20d%20in%20range(dims))%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20np.frombuffer(f.read()%2C%20dtype%3Dnp.uint8).reshape(shape)%0A%0A%20%20%20%20data%20%3D%20read_idx('train-images.idx3-ubyte')%0A%20%20%20%20labels%20%3D%20read_idx('train-labels.idx1-ubyte')%0A%20%20%20%20%23Select%20the%20images%0A%20%20%20%20digits%20%3D%20data%5Blabels%20%3D%3D%20number%5D%0A%20%20%20%20train%20%3D%20digits%5B%3An%5D%0A%0A%20%20%20%20plt.figure(figsize%3D(20%2C10))%0A%20%20%20%20for%20i%20in%20range(10)%3A%0A%20%20%20%20%20%20%20%20plt.subplot(2%2C5%2Ci%2B1)%0A%20%20%20%20%20%20%20%20plt.imshow(digits%5Bnp.random.randint(0%2Cdigits.shape%5B0%5D)%5D)%0A%20%20%20%20plt.show()%0A%20%20%20%20return%20np%2C%20pd%2C%20plt%2C%20train%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%23%23%20Regularized%20Barycenters%20using%20Mosek%20Fusion%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np)%3A%0A%20%20%20%20from%20mosek.fusion%20import%20Model%2C%20Domain%2C%20Expr%2C%20ObjectiveSense%2C%20Matrix%2C%20Var%2C%20SolutionStatus%2C%20SolutionType%0A%20%20%20%20import%20time%0A%20%20%20%20import%20sys%0A%0A%20%20%20%20class%20Wasserstein_Fusion%3A%0A%0A%20%20%20%20%20%20%20%20def%20__init__(self)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20self.time%20%3D%200.0%0A%20%20%20%20%20%20%20%20%20%20%20%20self.M%20%3D%20Model('Wasserstein')%0A%20%20%20%20%20%20%20%20%20%20%20%20self.result%20%3D%20None%0A%0A%0A%20%20%20%20%20%20%20%20def%20single_pmf(self%2C%20data%20%3D%20None%2C%20img%3DFalse)%3A%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20'''%20Takes%20a%20image%20or%20array%20of%20images%20and%20extracts%20the%20probabilty%20mass%20function'''%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20not%20img%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20v%3D%5B%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20for%20image%20in%20data%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20arr%20%3D%20np.asarray(image).ravel(order%3D'K')%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20v.append(arr%2Fnp.sum(arr))%0A%20%20%20%20%20%20%20%20%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20v%20%3D%20np.asarray(data).ravel(order%3D'K')%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20v%20%3D%20v%2Fnp.sum(v)%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20v%0A%0A%20%20%20%20%20%20%20%20def%20ms_distance(self%2C%20m%20%2Cn%2C%20constant%3DFalse)%3A%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20'''%20Squared%20Euclidean%20distance%20calculation%20between%20the%20pixels%20'''%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20constant%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20d%20%3D%20np.ones((m%2Cm))%0A%20%20%20%20%20%20%20%20%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20d%20%3D%20np.empty((m%2Cm))%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20coor%20%3D%20%5B%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20for%20i%20in%20range(n)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20for%20j%20in%20range(n)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20coor.append(np.array(%5Bi%2Cj%5D))%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20for%20i%20in%20range(m)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20for%20j%20in%20range(m)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20d%5Bi%5D%5Bj%5D%20%3D%20np.linalg.norm(coor%5Bi%5D-coor%5Bj%5D)**2%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20d%0A%0A%20%20%20%20%20%20%20%20def%20Wasserstein_Distance(self%2C%20bc%20%2Cdata%2C%20img%20%3D%20False)%3A%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20'''%20Calculation%20of%20wasserstein%20distance%20between%20a%20barycenter%20and%20an%20image%20by%20solving%20the%20minimization%20problem%20'''%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20v%20%3D%20np.array(self.single_pmf(data%2C%20img))%0A%20%20%20%20%20%20%20%20%20%20%20%20n%20%3D%20v.shape%5B0%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20d%20%3D%20self.ms_distance(n%2Cdata.shape%5B1%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20with%20Model('Wasserstein')%20as%20M%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%23Add%20variable%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20pi%20%3D%20M.variable('pi'%2C%5Bn%2Cn%5D%2C%20Domain.greaterThan(0.0))%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%23Add%20constraints%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20M.constraint('c1'%20%2C%20Expr.sum(pi%2C0)%2C%20Domain.equalsTo(v))%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20M.constraint('c2'%20%2C%20Expr.sum(pi%2C1)%2C%20Domain.equalsTo(bc))%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20M.objective('Obj.'%20%2C%20ObjectiveSense.Minimize%2C%20Expr.dot(d%2C%20pi))%0A%0A%20%20%20%20%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%20%20%20%20%20objective%20%3D%20M.primalObjValue()%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20objective%0A%0A%20%20%20%20%20%20%20%20def%20Wasserstein_BaryCenter(self%2Cdata)%3A%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20M%20%3D%20self.M%0A%20%20%20%20%20%20%20%20%20%20%20%20start_time%20%3D%20time.time()%0A%20%20%20%20%20%20%20%20%20%20%20%20k%20%3D%20data.shape%5B0%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20v%20%3D%20np.array(self.single_pmf(data))%0A%20%20%20%20%20%20%20%20%20%20%20%20n%20%3D%20v.shape%5B1%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20d%20%3D%20self.ms_distance(n%2Cdata.shape%5B1%5D)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23Add%20variables%20%20%20%0A%20%20%20%20%20%20%20%20%20%20%20%20mu%20%3D%20M.variable('Mu'%2C%20n%2C%20Domain.greaterThan(0.0))%20%20%0A%20%20%20%20%20%20%20%20%20%20%20%20pi%20%3D%20(M.variable('Pi'%2C%20%5Bk%2Cn%2Cn%5D%20%2C%20Domain.greaterThan(0.0)))%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23Add%20constraints%20%20%20%20%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23Constraint%20(1)%0A%20%20%20%20%20%20%20%20%20%20%20%20M.constraint('B'%2C%20Expr.sub(Expr.sum(pi%2C1)%20%2C%20Var.repeat(mu%2C1%2Ck).transpose())%2C%20Domain.equalsTo(0.0))%0A%20%20%20%20%20%20%20%20%20%20%20%20%23Constraint%20(2)%0A%20%20%20%20%20%20%20%20%20%20%20%20M.constraint('C'%2C%20Expr.sum(pi%2C2)%2C%20Domain.equalsTo(v))%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20M.objective('Obj'%20%2C%20ObjectiveSense.Minimize%2C%20Expr.sum(Expr.mul(Expr.mul(Expr.reshape(pi.asExpr()%2C%20k%2C%20n*n)%20%2C%20d.ravel())%2C%201%2Fk)))%0A%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%20self.result%20%3D%20mu.level()%0A%20%20%20%20%20%20%20%20%20%20%20%20M.selectedSolution(SolutionType.Interior)%0A%20%20%20%20%20%20%20%20%20%20%20%20self.objective%20%3D%20M.primalObjValue()%0A%20%20%20%20%20%20%20%20%20%20%20%20self.time%20%3D%20time.time()%20-%20start_time%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20mu.level()%0A%0A%20%20%20%20%20%20%20%20def%20Wasserstein_regBaryCenter(self%2Cdata%2C%20_lambda%20%3D%20None%2C%20relgap%20%3D%20None)%3A%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20M%20%3D%20self.M%0A%20%20%20%20%20%20%20%20%20%20%20%20start_time%20%3D%20time.time()%0A%20%20%20%20%20%20%20%20%20%20%20%20k%20%3D%20data.shape%5B0%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20v%20%3D%20np.array(self.single_pmf(data))%0A%20%20%20%20%20%20%20%20%20%20%20%20n%20%3D%20v.shape%5B1%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20d%20%3D%20self.ms_distance(n%2Cdata.shape%5B1%5D)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20not%20_lambda%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_lambda%20%3D%2060%2Fnp.median(d.ravel())%0A%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23Add%20variables%20%20%20%0A%20%20%20%20%20%20%20%20%20%20%20%20mu%20%3D%20M.variable('Mu'%2C%20n%2C%20Domain.greaterThan(0.0))%20%20%0A%20%20%20%20%20%20%20%20%20%20%20%20pi%20%3D%20(M.variable('Pi'%2C%20%5Bk%2Cn%2Cn%5D%20%2C%20Domain.greaterThan(0.0)))%0A%20%20%20%20%20%20%20%20%20%20%20%20z%20%3D%20M.variable('z'%2C%20%5Bk%2Cn*n%5D)%20%23Artificial%20variable%20%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23Add%20constraints%0A%20%20%20%20%20%20%20%20%20%20%20%20%23Intermediate%20conic%20constraints%20in%20form%20z%20%3C%3D%20-pi%20log(pi)%20%0A%20%20%20%20%20%20%20%20%20%20%20%20for%20i%20in%20range(1%2Ck%2B1)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20M.constraint(Expr.hstack(Expr.constTerm(n*n%2C%201.0)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20Expr.reshape(pi.asExpr()%2C%20k*n*n).slice(0%2B(i-1)*n*n%2C%20n*n*i)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20Expr.reshape(z.asExpr()%2C%20k*n*n).slice(0%2B(i-1)*n*n%2C%20n*n*i))%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20Domain.inPExpCone())%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23Constraint%20(1)%0A%20%20%20%20%20%20%20%20%20%20%20%20M.constraint('B'%2C%20Expr.sub(Expr.sum(pi%2C1)%20%2C%20Var.repeat(mu%2C1%2Ck).transpose())%2C%20Domain.equalsTo(0.0))%0A%20%20%20%20%20%20%20%20%20%20%20%20%23Constraint%20(2)%0A%20%20%20%20%20%20%20%20%20%20%20%20M.constraint('C'%2C%20Expr.sum(pi%2C2)%2C%20Domain.equalsTo(v))%0A%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20M.objective('Obj'%20%2C%20ObjectiveSense.Minimize%2C%20Expr.sum(Expr.mul(Expr.add(Expr.mul(Expr.reshape(pi.asExpr()%2C%20k%2C%20n*n)%2C%20d.ravel())%2CExpr.mul(Expr.sum(z%2C1)%2C-1%2F_lambda))%2C%201%2Fk)))%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23relgap%20is%20set%20in%20case%20of%20approximation%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20relgap%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20M.setSolverParam(%22intpntCoTolRelGap%22%2C%20relgap)%0A%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%20self.result%20%3D%20mu.level()%0A%20%20%20%20%20%20%20%20%20%20%20%20self.objective%20%3D%20M.primalObjValue()%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20self.time%20%3D%20time.time()%20-%20start_time%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20self.result%0A%0A%0A%0A%20%20%20%20%20%20%20%20def%20reset(self)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20self.M%20%3D%20Model('Wasserstein')%0A%20%20%20%20return%20Wasserstein_Fusion%2C%20time%0A%0A%0A%40app.cell%0Adef%20_(Wasserstein_Fusion%2C%20train)%3A%0A%20%20%20%20fusion_model%20%3D%20Wasserstein_Fusion()%0A%20%20%20%20f_bc%20%3D%20fusion_model.Wasserstein_regBaryCenter(train)%0A%20%20%20%20print('%5CnTime%20Spent%20to%20solve%20problem%20with%20Fusion%3A%20%5Cn%20%7B0%7D'.format(fusion_model.time))%0A%20%20%20%20print('Time%20Spent%20in%20solver%3A%20%5Cn%20%7B0%7D'.format(fusion_model.M.getSolverDoubleInfo(%22optimizerTime%22)))%0A%20%20%20%20print('Objective%3A%20%5Cn%20%7B0%7D'.format(fusion_model.objective))%0A%20%20%20%20return%20f_bc%2C%20fusion_model%0A%0A%0A%40app.cell%0Adef%20_(Wasserstein_Fusion%2C%20train)%3A%0A%20%20%20%20nonReg_model%20%3D%20Wasserstein_Fusion()%0A%20%20%20%20nonReg%20%3D%20nonReg_model.Wasserstein_BaryCenter(train)%0A%20%20%20%20print('%5CnTime%20Spent%20to%20solve%20non-regularized%20problem%20with%20Fusion%3A%20%5Cn%20%7B0%7D'.format(nonReg_model.time))%0A%20%20%20%20print('Time%20Spent%20in%20solver%3A%20%5Cn%20%7B0%7D'.format(nonReg_model.M.getSolverDoubleInfo(%22optimizerTime%22)))%0A%20%20%20%20print('The%20average%20Wasserstein%20distance%20between%20digits%20and%20the%20barycenter%3A%20%5Cn%20%7B0%7D'.format(nonReg_model.objective))%0A%20%20%20%20return%20nonReg%2C%20nonReg_model%0A%0A%0A%40app.cell%0Adef%20_(f_bc%2C%20nonReg%2C%20nonReg_model%2C%20np%2C%20plt)%3A%0A%20%20%20%20plt.figure(figsize%3D(15%2C8))%0A%20%20%20%20plt.subplot(1%2C3%2C1)%0A%20%20%20%20plt.imshow(np.reshape(f_bc%2C(28%2C28)))%0A%20%20%20%20plt.title('Regularized%20Barycenter')%0A%20%20%20%20plt.subplot(1%2C3%2C2)%0A%20%20%20%20plt.imshow(np.reshape(nonReg%2C%20(28%2C28)))%0A%20%20%20%20plt.title('Non-Regularized%20Interior%20Point%20Barycenter')%0A%20%20%20%20plt.imshow(np.reshape(nonReg%2C%20(28%2C28)))%0A%20%20%20%20plt.subplot(1%2C3%2C3)%0A%20%20%20%20plt.title('Non-Regularized%20Basic%20Point%20Barycenter')%0A%20%20%20%20plt.imshow(np.reshape(nonReg_model.result%2C%20(28%2C28)))%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%24%5Cquad%24%20The%20interiror%20point%20solution%20is%20different%20than%20the%20basic%20solution%2C%20however%20this%20is%20just%20a%20coincidence.%20The%20interior%20point%20solution%20gives%20the%20convex%20combination%20of%20the%20extreme%20points%20if%20there%20is%20infinetly%20many%20optimal%20solutions%20but%20this%20is%20not%20always%20the%20case%20in%20this%20problem.%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%24%5Cquad%20%24%20Solving%20the%20problem%20even%20for%20just%202%20images%20takes%20a%20long%20time.%20However%2C%20when%20the%20output%20of%20the%20iterations%20are%20investigated%20the%20solver%20is%20spending%20most%20of%20the%20time%20in%20earning%20very%20little%20improvements.%20Since%20the%20regularization%20term%20is%20an%20artificial%20addition%20to%20the%20problem%2C%20it%20is%20sensible%20to%20test%20if%20the%20approximation%20with%20a%20little%20loss%20of%20accuracy%20effects%20the%20values%20and%20images%20significantly%20or%20not%2C%20by%20using%20%22intpntCoTolRelGap%22%20parameter.%20This%20parameter%20controls%20the%20relative%20gap%20termination%20tolerance%20of%20the%20conic%20optimizer.%0A%20%20%20%20In%20the%20next%20problem%20the%20parameter%20is%20increased%20from%201.0e-7%20to%201.0e-3%20in%20order%20to%20terminate%20Mosek%20with%20a%20little%20bit%20less%20accuracy%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%20Obtaining%20approximate%20values%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(Wasserstein_Fusion%2C%20train)%3A%0A%20%20%20%20fusion_model2%20%3D%20Wasserstein_Fusion()%0A%20%20%20%20f_bc2%20%3D%20fusion_model2.Wasserstein_regBaryCenter(train%2C%20relgap%20%3D%20%221.0e-3%22)%0A%20%20%20%20print('%5CnTime%20Spent%20to%20solve%20problem%20with%20Fusion%3A%20%5Cn%20%7B0%7D'.format(fusion_model2.time))%0A%20%20%20%20print('Time%20Spent%20in%20solver%3A%20%5Cn%20%7B0%7D'.format(fusion_model2.M.getSolverDoubleInfo(%22optimizerTime%22)))%0A%20%20%20%20print('Objective%3A%20%5Cn%20%7B0%7D'.format(fusion_model2.objective))%0A%20%20%20%20return%20f_bc2%2C%20fusion_model2%0A%0A%0A%40app.cell%0Adef%20_(f_bc%2C%20np%2C%20plt)%3A%0A%20%20%20%20plt.figure(figsize%3D(10%2C10))%0A%20%20%20%20plt.subplot(1%2C2%2C1)%0A%20%20%20%20fus_bc%20%3D%20np.reshape(f_bc%2C(28%2C28))%0A%20%20%20%20plt.imshow(fus_bc)%0A%20%20%20%20plt.title('Optimal%20Barycenter')%0A%20%20%20%20plt.subplot(1%2C2%2C2)%0A%20%20%20%20fus_bc2%20%3D%20np.reshape(f_bc%2C(28%2C28))%0A%20%20%20%20plt.imshow(fus_bc)%0A%20%20%20%20plt.title('Approximate%20Barycenter')%0A%20%20%20%20plt.show()%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(f_bc%2C%20f_bc2%2C%20pd%2C%20plt)%3A%0A%20%20%20%20error%20%3D%20pd.Series(f_bc%20-%20f_bc2)%0A%20%20%20%20plt.plot(error)%0A%20%20%20%20plt.title('Error%20of%20Aproximation')%0A%20%20%20%20plt.xlabel(%22Pixels%22)%0A%20%20%20%20plt.ylabel('Error%20Values')%0A%20%20%20%20pd.DataFrame(error.describe()%2C%20columns%3D%5B'Stats'%5D).transpose()%0A%20%20%20%20plt.show()%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(fusion_model%2C%20fusion_model2%2C%20plt)%3A%0A%20%20%20%20plt.figure(figsize%3D(12%2C7))%0A%0A%20%20%20%20total_t%20%3D%20%5Bfusion_model.time%2C%20fusion_model2.time%5D%0A%20%20%20%20solver_t%20%3D%20%5Bfusion_model.M.getSolverDoubleInfo(%22optimizerTime%22)%2C%20fusion_model2.M.getSolverDoubleInfo(%22optimizerTime%22)%5D%0A%20%20%20%20%23Total%20time%20plot%0A%20%20%20%20plt.subplot(1%2C2%2C1)%0A%20%20%20%20plt.bar(%5B'Optimal'%2C%20'Approximation'%5D%2C%20height%3D%20total_t%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20width%3D0.4%2C%20color%3D(0.3%2C%200.6%2C%200.2%2C%200.5))%0A%20%20%20%20plt.ylabel(%22Total%20Time%20(s)%22)%0A%20%20%20%20plt.title(%22Comparison%20of%20Total%20Time%22)%0A%0A%20%20%20%20%23Solver%20time%20plot%0A%20%20%20%20plt.subplot(1%2C2%2C2)%0A%20%20%20%20plt.bar(%5B'Optimal'%2C%20'Approximation'%5D%2C%20height%3Dsolver_t%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20width%3D0.4%2C%20color%3D(0.5%2C%200.6%2C%200.9%2C%200.8))%0A%20%20%20%20plt.ylabel(%22Solver%20Time%20(s)%22)%0A%20%20%20%20plt.title(%22Comparison%20of%20Solver%20Time%22)%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%24%5Cquad%24The%20two%20barycenter%20images%20that%20obtained%20from%20the%20optimal%20solution%20and%20approximate%20solution%20seem%20almost%20identical%20to%20each%20other.%20In%20addition%20statistical%20description%20of%20values%20of%20the%20error%20between%20them%20is%20presented%20above%20with%20a%20plot.%20The%20mean%20and%20the%20even%20extreme%20values%20are%20small%20and%20indicates%20that%20the%20reduction%20of%20time%20obtained%20by%20approximation%20totally%20compensates%20the%20error.%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%23%23%23%20Modeling%20the%20same%20problem%20with%20CVXPY%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20time)%3A%0A%20%20%20%20import%20cvxpy%20as%20cp%0A%20%20%20%20class%20Wasserstein_CVXPY%3A%0A%0A%20%20%20%20%20%20%20%20def%20__init__(self)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20self.time%20%3D%200.0%0A%20%20%20%20%20%20%20%20%20%20%20%20self.result%20%3D%20None%0A%20%20%20%20%20%20%20%20%20%20%20%20self.prob%20%3D%20None%0A%0A%20%20%20%20%20%20%20%20def%20single_pmf(self%2C%20data%20%3D%20None%2C%20img%3DFalse)%3A%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20'''%20Takes%20a%20image%20or%20array%20of%20images%20and%20extracts%20the%20probabilty%20mass%20function'''%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20not%20img%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20v%3D%5B%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20for%20image%20in%20data%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20arr%20%3D%20np.asarray(image).ravel(order%3D'K')%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20v.append(arr%2Fnp.sum(arr))%0A%20%20%20%20%20%20%20%20%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20v%20%3D%20np.asarray(data).ravel(order%3D'K')%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20v%20%3D%20v%2Fnp.sum(v)%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20v%0A%0A%20%20%20%20%20%20%20%20def%20ms_distance(self%2C%20m%20%2Cn%2C%20constant%3DFalse)%3A%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20'''%20Squared%20Euclidean%20distance%20calculation%20between%20the%20pixels%20'''%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20constant%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20d%20%3D%20np.ones((m%2Cm))%0A%20%20%20%20%20%20%20%20%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20d%20%3D%20np.empty((m%2Cm))%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20coor%20%3D%20%5B%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20for%20i%20in%20range(n)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20for%20j%20in%20range(n)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20coor.append(np.array(%5Bi%2Cj%5D))%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20for%20i%20in%20range(m)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20for%20j%20in%20range(m)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20d%5Bi%5D%5Bj%5D%20%3D%20np.linalg.norm(coor%5Bi%5D-coor%5Bj%5D)**2%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20d%0A%0A%20%20%20%20%20%20%20%20def%20Wasserstein_Distance(self%2C%20bc%20%2Cdata%2C%20img%20%3D%20False)%3A%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20'''%20Calculation%20of%20wasserstein%20distance%20between%20a%20barycenter%20and%20an%20image%20by%20solving%20%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20the%20minimization%20problem%20'''%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20v%20%3D%20np.array(self.single_pmf(data%2C%20img))%0A%20%20%20%20%20%20%20%20%20%20%20%20n%20%3D%20v.shape%5B0%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20d%20%3D%20self.ms_distance(n%2Cdata.shape%5B1%5D)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20pi%20%3D%20cp.Variable((n%2Cn)%2C%20nonneg%3DTrue)%0A%20%20%20%20%20%20%20%20%20%20%20%20obj%20%3D%20cp.Minimize((np.ones(n).T%20%40%20cp.multiply(d%2Cpi)%20%40%20np.ones(n)))%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20Cons%3D%5B%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20Cons.append((np.ones(n)%20%40%20pi).T%20%3D%3D%20bc)%0A%20%20%20%20%20%20%20%20%20%20%20%20Cons.append((pi%20%40%20np.ones(n))%20%3D%3D%20v)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20prob%20%3D%20cp.Problem(obj%2C%20constraints%3D%20Cons)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20prob.solve(solver%3Dcp.MOSEK%2C%20verbose%20%3D%20True)%0A%0A%20%20%20%20%20%20%20%20def%20Wasserstein_BaryCenter(self%2Cdata)%3A%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20'''%20Calculation%20of%20wasserstein%20barycenter%20of%20given%20images%20by%20solving%20the%20minimization%20problem%20'''%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20start_time%20%3D%20time.time()%0A%20%20%20%20%20%20%20%20%20%20%20%20k%20%3D%20data.shape%5B0%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20v%20%3D%20np.array(self.single_pmf(data))%0A%20%20%20%20%20%20%20%20%20%20%20%20n%20%3D%20v.shape%5B1%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20d%20%3D%20self.ms_distance(n%2Cdata.shape%5B1%5D)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23Add%20variables%0A%20%20%20%20%20%20%20%20%20%20%20%20pi%3D%20%5B%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20t%3D%20%5B%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20mu%20%3D%20cp.Variable(n%2C%20nonneg%20%3D%20True)%0A%20%20%20%20%20%20%20%20%20%20%20%20for%20i%20in%20range(k)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20pi.append(cp.Variable((n%2Cn)%2C%20nonneg%20%3D%20True))%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20t.append(cp.Variable(nonneg%20%3D%20True))%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20obj%20%3D%20cp.Minimize(np.sum(t)%2Fk)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23Add%20constraints%0A%20%20%20%20%20%20%20%20%20%20%20%20Cons%3D%5B%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20for%20i%20in%20range(k)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20Cons.append(%20t%5Bi%5D%20%3E%3D%20np.ones(n).T%20%40%20%20cp.multiply(d%2Cpi%5Bi%5D)%20%40%20np.ones(n)%20)%20%23Constraint%20(1)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20Cons.append(%20(np.ones(n)%20%40%20pi%5Bi%5D).T%20%3D%3D%20mu)%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%23Constraint%20(2)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20Cons.append(%20(pi%5Bi%5D%20%40%20np.ones(n))%20%3D%3D%20v%5Bi%5D)%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%23Constraint%20(3)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20self.prob%20%3D%20cp.Problem(obj%2C%20constraints%3D%20Cons)%0A%20%20%20%20%20%20%20%20%20%20%20%20self.result%20%3D%20self.prob.solve(solver%3Dcp.MOSEK%2Cverbose%20%3D%20True)%0A%20%20%20%20%20%20%20%20%20%20%20%20self.time%20%3D%20time.time()%20-%20start_time%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20mu.value%0A%0A%20%20%20%20%20%20%20%20def%20Wasserstein_regBaryCenter(self%2Cdata%2C%20_lambda%20%3D%20None%2C%20relgap%3D%221.0e-7%22)%3A%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20'''%20Calculation%20of%20wasserstein%20barycenter%20of%20given%20%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20images%20by%20solving%20a%20entropy%20regularized%20minimization%20problem%20'''%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20start_time%20%3D%20time.time()%0A%20%20%20%20%20%20%20%20%20%20%20%20k%20%3D%20data.shape%5B0%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20v%20%3D%20np.array(self.single_pmf(data))%0A%20%20%20%20%20%20%20%20%20%20%20%20n%20%3D%20v.shape%5B1%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20d%20%3D%20self.ms_distance(n%2Cdata.shape%5B1%5D)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20not%20_lambda%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_lambda%20%3D%2060%2Fnp.median(d.ravel())%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23Add%20variables%0A%20%20%20%20%20%20%20%20%20%20%20%20pi%3D%20%5B%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20t%3D%20%5B%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20mu%20%3D%20cp.Variable(n%2C%20nonneg%20%3D%20True)%0A%20%20%20%20%20%20%20%20%20%20%20%20for%20i%20in%20range(k)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20pi.append(cp.Variable((n%2Cn)%2C%20nonneg%20%3D%20True))%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20t.append(cp.Variable(nonneg%20%3D%20True))%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20obj%20%3D%20cp.Minimize((np.sum(t)%20-%20(1%2F_lambda)*np.sum(cp.sum(cp.entr(pi%5Bi%5D))%20for%20i%20in%20range(k)))%2Fk)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23Add%20constraints%0A%20%20%20%20%20%20%20%20%20%20%20%20Cons%3D%5B%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20for%20i%20in%20range(k)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20Cons.append(%20t%5Bi%5D%20%3E%3D%20(np.ones(n).T%20%40%20%20cp.multiply(d%2Cpi%5Bi%5D)%20%40%20np.ones(n)))%23Constraint%20(1)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20Cons.append(%20(np.ones(n)%20%40%20pi%5Bi%5D).T%20%3D%3D%20mu)%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%23Constraint%20(2)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20Cons.append(%20(pi%5Bi%5D%20%40%20np.ones(n))%20%3D%3D%20v%5Bi%5D)%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%23Constraint%20(3)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20self.prob%20%3D%20cp.Problem(obj%2C%20constraints%3D%20Cons)%0A%20%20%20%20%20%20%20%20%20%20%20%20self.result%20%3D%20self.prob.solve(solver%3Dcp.MOSEK%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20verbose%20%3D%20True%2C%20%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20mosek_params%20%3D%20%7B%22MSK_DPAR_INTPNT_CO_TOL_REL_GAP%22%20%3A%20relgap%20%7D)%0A%20%20%20%20%20%20%20%20%20%20%20%20self.time%20%3D%20time.time()%20-%20start_time%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20mu.value%0A%0A%20%20%20%20%20%20%20%20def%20reset(self)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20self.prob%20%3D%20None%0A%20%20%20%20%20%20%20%20%20%20%20%20self.result%20%3D%20None%0A%20%20%20%20return%20(Wasserstein_CVXPY%2C)%0A%0A%0A%40app.cell%0Adef%20_(Wasserstein_CVXPY%2C%20train)%3A%0A%20%20%20%20cvxpy_model%20%3D%20Wasserstein_CVXPY()%0A%20%20%20%20cvxpy_result%20%3D%20cvxpy_model.Wasserstein_regBaryCenter(train%2C%20relgap%20%3D%201.0e-3)%0A%20%20%20%20print('%5CnTime%20Spent%20to%20solve%20problem%20with%20CVXPY%3A%20%5Cn%20%7B0%7D'.format(cvxpy_model.time))%0A%20%20%20%20print('Time%20Spent%20in%20solver%3A%20%5Cn%20%7B0%7D'.format(cvxpy_model.prob.solver_stats.solve_time))%0A%20%20%20%20print('The%20average%20Wasserstein%20distance%20between%20digits%20and%20the%20barycenter%3A%20%5Cn%20%7B0%7D'.format(cvxpy_model.result))%0A%20%20%20%20return%20cvxpy_model%2C%20cvxpy_result%0A%0A%0A%40app.cell%0Adef%20_(cvxpy_result%2C%20np%2C%20plt)%3A%0A%20%20%20%20plt.imshow(np.reshape(cvxpy_result.squeeze()%2C%20(28%2C28)))%0A%20%20%20%20plt.title('Barycenter')%0A%20%20%20%20plt.show()%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(cvxpy_model%2C%20fusion_model2%2C%20plt)%3A%0A%20%20%20%20plt.figure(figsize%3D(12%2C6))%0A%0A%20%20%20%20total_t50%20%3D%20%5Bfusion_model2.time%2C%20cvxpy_model.time%5D%0A%20%20%20%20solver_t50%20%3D%20%5Bfusion_model2.M.getSolverDoubleInfo(%22optimizerTime%22)%2C%20cvxpy_model.prob.solver_stats.solve_time%5D%0A%0A%20%20%20%20%23Total%20time%20plot%0A%20%20%20%20plt.subplot(1%2C2%2C1)%0A%20%20%20%20plt.bar(%5B'Fusion'%2C%20'CVXPY'%5D%2C%20height%3D%20total_t50%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20width%3D0.4%2C%20color%3D(0.3%2C%200.6%2C%200.2%2C%200.5))%0A%20%20%20%20plt.ylabel(%22Total%20Time%20(s)%22)%0A%20%20%20%20plt.title(%22Comparison%20of%20Total%20Time%22)%0A%0A%20%20%20%20%23Solver%20time%20plot%0A%20%20%20%20plt.subplot(1%2C2%2C2)%0A%20%20%20%20plt.bar(%5B'Fusion'%2C'CVXPY'%5D%2C%20height%3Dsolver_t50%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20width%3D0.4%2C%20color%3D(0.5%2C%200.6%2C%200.9%2C%200.8))%0A%20%20%20%20plt.ylabel(%22Solver%20Time%20(s)%22)%0A%20%20%20%20plt.title(%22Comparison%20of%20Solver%20Time%22)%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%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
a850d0399eba50c4a71cbe40178c7468abc43bf35e6e9b1bcc92d2595f2c3f6b