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%20Wasserstein%20Barycenters%20using%20Mosek%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%22Wasserstein%20Distance%20is%20a%20way%20to%20measure%20the%20distance%20between%20two%20probabilty%20distributions.%20It%20allows%20to%20summarize%2C%20compare%2C%20match%20and%20reduce%20the%20dimensionality%20of%20the%20emprical%20probability%20measures%20to%20carry%20out%20some%20machine%20learning%20fundamentals.%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%20Wasserstein%20distance%20of%20order%20%24p%24%20between%20two%20probabilty%20measures%20%24%5Cmu%24%20and%20%20%20%20%20%20%24%5Cupsilon%24%20in%20%24P(%5COmega)%24%20is%20defined%20as%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20W_p(%5Cmu%2C%5Cupsilon)%20%5Coverset%7B%5Cunderset%7B%5Cmathrm%7Bdef%7D%7D%7B%7D%7D%7B%3D%7D%20%5Cbigg(%20%5Cunderset%7B%5Cpi%20%5Cin%20%5CPi%7B(%5Cmu%2C%20%5Cupsilon)%7D%7D%7B%5Ctext%7Binf%7D%7D%20%5Cint_%7B%5COmega%5E2%7D%20D(X_i%2CY_j)%5Ep%20d%5Cpi(x%2Cy)%5Cbigg)%5E%7B1%2Fp%7D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20where%20%24%5CPi(%5Cmu%2C%20%5Cupsilon)%24%20is%20the%20set%20of%20all%20probability%20measures%20on%20%24%5COmega%5E2%24%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%20If%20the%20distributions%20are%20discrete%20%24W_p(%5Cmu%2C%5Cupsilon)%24%20is%20equilavent%20to%20the%20objective%20%20%20%20%20of%20the%20following%20LP%20model%3A%20%20%20%0A%0A%20%20%20%20%24%24%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%20%20%20%20%24%24%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%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%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%20on%20%24%5COmega%24.%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%22There%20are%20more%20efficient%20ways%20to%20approximate%20this%20metric%20but%20LP%20approach%20will%20be%20applied%20in%20order%20to%20compare%20the%20performance%20and%20modeling%20structure%20of%20Fusion%2C%20Pyomo%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(r%22%22%22%23%23%20Wasserstein%20Barycenter%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%20Wasserstein%20barycenter%20problem%20involves%20all%20Wasserstein%20distances%20from%20one%20to%20many%20measures.%20Given%20measures%20%24%5Cupsilon_1%2C%5Cldots%2C%5Cupsilon_N%24%20we%20want%20to%20find%20the%20measure%20which%20minimizes%20the%20sum%20of%20distances%20to%20the%20given%20measures%2C%20that%20is%20solve%20the%20problem%3A%0A%20%20%20%20%24%24%5Ctext%7Bminimize%7D_%5Cmu%5Csum_%7Bi%3D1%7D%20%5Clambda_i%20W_p(%5Cmu%2C%5Cupsilon_i)%24%24%0A%20%20%20%20for%20some%20fixed%20system%20%24%5Clambda_i%24%20of%20weights%20of%20distances%20to%20specific%20distributions%20that%20satisfies%20%24%24%5Csum_%7Bi%3D1%7D%5Clambda_i%20%3D%201.%24%24%0A%20%20%20%20For%20simplicity%20uniform%20weights%20are%20used%20in%20this%20problem.%20Then%20the%20barycenters%20problem%20becomes%3A%0A%20%20%20%20%24%24%5Ctext%7Bminimize%7D_%5Cmu%20%5Cfrac1N%20%5Csum_%7Bi%3D1%7D%5E%7BN%7D%20W_p(%5Cmu%2C%5Cupsilon_i).%24%24%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%22In%20this%20problem%2C%20Wasserstein%20Barycenter%20of%20One's%20are%20visualized%20using%20images%20with%20size%20%2428x28%24%20using%20%2420%24%20handwriten%20'1'%20digits%20from%20MNIST%20database%20http%3A%2F%2Fyann.lecun.com%2Fexdb%2Fmnist%2F.%20Computations%20are%20carried%20out%20by%20Intel(R)%20Xeon(R)%20CPU%20E5-2687W%20v4%20%40%203.00GHz%20processor.%20Similar%20experiments%20are%20performed%20by%20Cuturi%20and%20Doucet%20in%20http%3A%2F%2Fproceedings.mlr.press%2Fv32%2Fcuturi14.pdf.%22%22%22)%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%20matplotlib.pyplot%20as%20plt%0A%20%20%20%20n%20%3D%2020%0A%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%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%20ones%20%3D%20data%5Blabels%20%3D%3D%201%5D%0A%20%20%20%20train_1%20%3D%20ones%5B%3An%5D%0A%20%20%20%20plt.figure(figsize%3D(20%2C%2010))%0A%20%20%20%20for%20i%20in%20range(10)%3A%0A%20%20%20%20%20%20%20%20plt.subplot(2%2C%205%2C%20i%20%2B%201)%0A%20%20%20%20%20%20%20%20plt.imshow(ones%5Bnp.random.randint(0%2C%20ones.shape%5B0%5D)%5D%2C%20cmap%3D'viridis')%0A%20%20%20%20%20%20%20%20plt.axis('off')%0A%20%20%20%20plt.show()%0A%20%20%20%20return%20np%2C%20ones%2C%20plt%2C%20train_1%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%23%23%20Barycenters%20using%20Mosek%20Fusion%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%20final%20barycenter%20problem%20is%20as%20follows.%20We%20choose%20%24p%3D2%24.%0A%0A%20%20%20%20%24%24%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%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%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%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%20Euclidean%20distance%20between%20pixels%20and%20%24N%24%20is%20the%20number%20of%20samples.%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_(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%20%20%20%20%20%20%20%20def%20single_pmf(self%2C%20data%3DNone%2C%20img%3DFalse)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%22%22%22%20Takes%20a%20image%20or%20array%20of%20images%20and%20extracts%20the%20probabilty%20mass%20function%22%22%22%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%20%3D%20%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%20%2F%20np.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%20%2F%20np.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%2C%20n%2C%20constant%3DFalse)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%22%22%22%20Squared%20Euclidean%20distance%20calculation%20between%20the%20pixels%20%22%22%22%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%2C%20m))%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%2C%20m))%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%2C%20j%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%20-%20coor%5Bj%5D)%20**%202%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%2C%20data%2C%20img%3DFalse)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%22%22%22%20Calculation%20of%20wasserstein%20distance%20between%20a%20barycenter%20and%20an%20image%20by%20solving%20the%20minimization%20problem%20%22%22%22%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%2C%20data.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%20pi%20%3D%20M.variable('pi'%2C%20%5Bn%2C%20n%5D%2C%20Domain.greaterThan(0.0))%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20M.constraint('c1'%2C%20Expr.sum(pi%2C%200)%2C%20Domain.equalsTo(v))%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20M.constraint('c2'%2C%20Expr.sum(pi%2C%201)%2C%20Domain.equalsTo(bc))%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20M.objective('Obj.'%2C%20ObjectiveSense.Minimize%2C%20Expr.dot(d%2C%20pi))%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%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%2C%20data)%3A%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%2C%20data.shape%5B1%5D)%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))%0A%20%20%20%20%20%20%20%20%20%20%20%20pi%20%3D%20M.variable('Pi'%2C%20%5Bk%2C%20n%2C%20n%5D%2C%20Domain.greaterThan(0.0))%0A%20%20%20%20%20%20%20%20%20%20%20%20M.constraint('B'%2C%20Expr.sub(Expr.sum(pi%2C%201)%2C%20Var.repeat(mu%2C%201%2C%20k).transpose())%2C%20Domain.equalsTo(0.0))%0A%20%20%20%20%20%20%20%20%20%20%20%20M.constraint('C'%2C%20Expr.sum(pi%2C%202)%2C%20Domain.equalsTo(v))%0A%20%20%20%20%20%20%20%20%20%20%20%20M.objective('Obj'%2C%20ObjectiveSense.Minimize%2C%20Expr.sum(Expr.mul(Expr.mul(Expr.reshape(pi.asExpr()%2C%20k%2C%20n%20*%20n)%2C%20d.ravel())%2C%201%20%2F%20k)))%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%20%20%20%20%20%20%20%20%20%20%20%20return%20mu.level()%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_1)%3A%0A%20%20%20%20fusion_model%20%3D%20Wasserstein_Fusion()%0A%20%20%20%20f_bc%20%3D%20fusion_model.Wasserstein_BaryCenter(train_1)%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('The%20average%20Wasserstein%20distance%20between%20digits%20and%20the%20barycenter%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_(f_bc%2C%20np%2C%20plt)%3A%0A%20%20%20%20fus_bc%20%3D%20np.reshape(f_bc%2C%20(28%2C%2028))%0A%20%20%20%20plt.figure(figsize%3D(3%2C3))%0A%20%20%20%20plt.imshow(fus_bc%2C%20cmap%3D'viridis')%0A%20%20%20%20plt.title('Barycenter')%0A%20%20%20%20plt.axis('off')%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%20Modeling%20the%20same%20problem%20with%20Pyomo%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%22The%20same%20problem%20is%20formulated%20using%20Pyomo%20and%20solved%20using%20Mosek.%20Unlike%20Fusion%20Pyomo%20requires%20rules%20and%20summations%20to%20formulate%20the%20problem.%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20time)%3A%0A%20%20%20%20import%20pyomo.environ%20as%20pyo%0A%0A%20%20%20%20class%20Wasserstein_Pyomo%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._solver%20%3D%20'mosek'%0A%20%20%20%20%20%20%20%20%20%20%20%20self.M%20%3D%20pyo.ConcreteModel()%0A%0A%20%20%20%20%20%20%20%20def%20single_pmf(self%2C%20data%3DNone%2C%20img%3DFalse)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%22%22%22%20Takes%20a%20image%20or%20array%20of%20images%20and%20extracts%20the%20probabilty%20mass%20function%22%22%22%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%20%3D%20%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%20%2F%20np.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%20%2F%20np.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%2C%20n%2C%20constant%3DFalse)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%22%22%22%20Squared%20Euclidean%20distance%20calculation%20between%20the%20pixels%20%22%22%22%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%2C%20m))%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%2C%20m))%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%2C%20j%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%20-%20coor%5Bj%5D)%20**%202%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_BaryCenter(self%2C%20data)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%22%22%22%20Calculation%20of%20wasserstein%20barycenter%20of%20given%20images%20by%20solving%20the%20minimization%20problem%20%22%22%22%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%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%2C%20data.shape%5B1%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20M.i%20%3D%20range(n)%0A%20%20%20%20%20%20%20%20%20%20%20%20M.j%20%3D%20range(n)%0A%20%20%20%20%20%20%20%20%20%20%20%20M.k%20%3D%20range(k)%0A%20%20%20%20%20%20%20%20%20%20%20%20M.pi%20%3D%20pyo.Var(M.k%2C%20M.i%2C%20M.j%2C%20domain%3Dpyo.NonNegativeReals)%0A%20%20%20%20%20%20%20%20%20%20%20%20M.mu%20%3D%20pyo.Var(M.i%2C%20domain%3Dpyo.NonNegativeReals)%0A%20%20%20%20%20%20%20%20%20%20%20%20M.t%20%3D%20pyo.Var(M.k%2C%20domain%3Dpyo.NonNegativeReals)%0A%20%20%20%20%20%20%20%20%20%20%20%20M.obj%20%3D%20pyo.Objective(expr%3Dsum((M.t%5Bk%5D%20for%20k%20in%20M.k))%20%2F%20k%2C%20sense%3Dpyo.minimize)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20def%20c3_rule(model%2C%20k%2C%20j)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20return%20sum((model.pi%5Bk%2C%20i%2C%20j%5D%20for%20i%20in%20model.i))%20%3D%3D%20v%5Bk%5D%5Bj%5D%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20def%20c2_rule(model%2C%20k%2C%20i)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20return%20sum((model.pi%5Bk%2C%20i%2C%20j%5D%20for%20j%20in%20model.j))%20%3D%3D%20model.mu%5Bi%5D%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20def%20c1_rule(model%2C%20k)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20return%20sum((d%5Bi%5D%5Bj%5D%20*%20model.pi%5Bk%2C%20i%2C%20j%5D%20for%20i%20in%20model.i%20for%20j%20in%20model.j))%20%3C%3D%20model.t%5Bk%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20M.c3%20%3D%20pyo.Constraint(M.k%2C%20M.j%2C%20rule%3Dc3_rule)%0A%20%20%20%20%20%20%20%20%20%20%20%20M.c2%20%3D%20pyo.Constraint(M.k%2C%20M.i%2C%20rule%3Dc2_rule)%0A%20%20%20%20%20%20%20%20%20%20%20%20M.c1%20%3D%20pyo.Constraint(M.k%2C%20rule%3Dc1_rule)%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20M%0A%0A%20%20%20%20%20%20%20%20def%20run(self%2C%20data)%3A%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%20model%20%3D%20self.Wasserstein_BaryCenter(data)%0A%20%20%20%20%20%20%20%20%20%20%20%20opt%20%3D%20pyo.SolverFactory(self._solver)%0A%20%20%20%20%20%20%20%20%20%20%20%20self.result%20%3D%20opt.solve(model%2C%20tee%3DTrue)%0A%20%20%20%20%20%20%20%20%20%20%20%20self.time%20%3D%20time.time()%20-%20start_time%0A%20%20%20%20%20%20%20%20%20%20%20%20bc%20%3D%20%5B%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20%5Bbc.append(model.mu%5Bi%5D())%20for%20i%20in%20range(data.shape%5B1%5D%20*%20data.shape%5B2%5D)%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20np.array(bc)%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%20pyo.ConcreteModel()%0A%20%20%20%20return%20(Wasserstein_Pyomo%2C)%0A%0A%0A%40app.cell%0Adef%20_(Wasserstein_Pyomo%2C%20train_1)%3A%0A%20%20%20%20pyomo_model%20%3D%20Wasserstein_Pyomo()%0A%20%20%20%20p_bc%20%3D%20pyomo_model.run(train_1)%0A%20%20%20%20print('%5CnTime%20spent%20to%20solve%20problem%20with%20Pyomo%3A%20%5Cn%20%7B0%7D'.format(pyomo_model.time))%0A%20%20%20%20print('Time%20spent%20in%20solver%3A%20%5Cn%20%7B0%7D'.format(pyomo_model.result.solver.wallclock_time))%0A%20%20%20%20print('The%20average%20Wasserstein%20distance%20between%20digits%20and%20the%20barycenter%3A%20%5Cn%20%7B0%7D'.format(pyomo_model.M.obj()))%0A%20%20%20%20return%20p_bc%2C%20pyomo_model%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20p_bc%2C%20plt)%3A%0A%20%20%20%20pyo_bc%20%3D%20np.reshape(p_bc%2C%20(28%2C28))%0A%20%20%20%20%23print('Visualization%20of%20the%20barycenter%3A')%0A%20%20%20%20plt.imshow(pyo_bc%2C%20cmap%3D'viridis')%20%20%23%20or%20'plasma'%2C%20'magma'%0A%20%20%20%20plt.title('Barycenter')%0A%20%20%20%20plt.axis('off')%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%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%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%3DNone%2C%20img%3DFalse)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%22%22%22%20Takes%20a%20image%20or%20array%20of%20images%20and%20extracts%20the%20probabilty%20mass%20function%22%22%22%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%20%3D%20%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%20%2F%20np.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%20%2F%20np.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%2C%20n%2C%20constant%3DFalse)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%22%22%22%20Squared%20Euclidean%20distance%20calculation%20between%20the%20pixels%20%22%22%22%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%2C%20m))%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%2C%20m))%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%2C%20j%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%20-%20coor%5Bj%5D)%20**%202%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%2C%20data%2C%20img%3DFalse)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%22%22%22%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%22%22%22%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%2C%20data.shape%5B1%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20pi%20%3D%20cp.Variable((n%2C%20n)%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%2C%20pi)%20%40%20np.ones(n))%0A%20%20%20%20%20%20%20%20%20%20%20%20Cons%20%3D%20%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%20%20%20%20%20%20%20%20%20%20%20%20prob%20%3D%20cp.Problem(obj%2C%20constraints%3DCons)%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20prob.solve(solver%3Dcp.MOSEK%2C%20verbose%3DTrue)%0A%0A%20%20%20%20%20%20%20%20def%20Wasserstein_BaryCenter(self%2C%20data)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%22%22%22%20Calculation%20of%20wasserstein%20barycenter%20of%20given%20images%20by%20solving%20the%20minimization%20problem%20%22%22%22%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%2C%20data.shape%5B1%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20pi%20%3D%20%5B%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20t%20%3D%20%5B%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20mu%20%3D%20cp.Variable(n%2C%20nonneg%3DTrue)%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%2C%20n)%2C%20nonneg%3DTrue))%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20t.append(cp.Variable(nonneg%3DTrue))%0A%20%20%20%20%20%20%20%20%20%20%20%20obj%20%3D%20cp.Minimize(np.sum(t)%20%2F%20k)%0A%20%20%20%20%20%20%20%20%20%20%20%20Cons%20%3D%20%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(t%5Bi%5D%20%3E%3D%20np.ones(n).T%20%40%20cp.multiply(d%2C%20pi%5Bi%5D)%20%40%20np.ones(n))%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20Cons.append((np.ones(n)%20%40%20pi%5Bi%5D).T%20%3D%3D%20mu)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20Cons.append(pi%5Bi%5D%20%40%20np.ones(n)%20%3D%3D%20v%5Bi%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20self.prob%20%3D%20cp.Problem(obj%2C%20constraints%3DCons)%0A%20%20%20%20%20%20%20%20%20%20%20%20self.result%20%3D%20self.prob.solve(solver%3Dcp.MOSEK%2C%20verbose%3DTrue)%0A%20%20%20%20%20%20%20%20%20%20%20%20self.time%20%3D%20time.time()%20-%20start_time%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_1)%3A%0A%20%20%20%20cvxpy_model%20%3D%20Wasserstein_CVXPY()%0A%20%20%20%20result%20%3D%20cvxpy_model.Wasserstein_BaryCenter(train_1)%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%20result%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20plt%2C%20result)%3A%0A%20%20%20%20plt.imshow(np.reshape(result.squeeze()%2C%20(28%2C%2028))%2C%20cmap%3D'viridis')%0A%20%20%20%20plt.title('Barycenter')%0A%20%20%20%20plt.axis('off')%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%22Also%2C%20when%20the%20log%20data%20from%20Mosek%20is%20investigated%20the%20number%20of%20variables%20are%20same%20for%20each%20modeling%20language.%20However%2C%20while%20Pyomo%20and%20Fusion%20having%20same%20number%20of%20constraints%20CVXPY%20has%20a%20lot%20of%20extra%20constraints%20that%20resulted%20from%20nonnegativity%20of%20the%20variables.%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%20Comparison%20of%20Modeling%20Languages%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%20Time%20to%20solve%20the%20problem%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(cvxpy_model%2C%20fusion_model%2C%20pyomo_model)%3A%0A%20%20%20%20print('Total%20Time%3A')%0A%20%20%20%20print('Fusion%3A%20%7B0%7D'.format(fusion_model.time))%0A%20%20%20%20print('Pyomo%20%3A%20%7B0%7D'.format(pyomo_model.time))%0A%20%20%20%20print('CVXPY%20%3A%20%7B0%7D'.format(cvxpy_model.time))%0A%20%20%20%20print('%5CnTime%20spent%20in%20the%20solver%3A')%0A%20%20%20%20print('Fusion%3A%20%7B0%7D'.format(fusion_model.M.getSolverDoubleInfo(%22optimizerTime%22)))%0A%20%20%20%20print('Pyomo%20%3A%20%7B0%7D'.format(pyomo_model.result.solver.wallclock_time))%0A%20%20%20%20print('CVXPY%20%3A%20%7B0%7D'.format(cvxpy_model.prob.solver_stats.solve_time))%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(cvxpy_model%2C%20fusion_model%2C%20plt%2C%20pyomo_model)%3A%0A%20%20%20%20plt.style.use('seaborn-v0_8')%20%0A%20%20%20%20plt.figure(figsize%3D(10%2C4))%0A%0A%20%20%20%20total_t%20%3D%20%5Bfusion_model.time%2C%20pyomo_model.time%2C%20cvxpy_model.time%5D%0A%20%20%20%20solver_t%20%3D%20%5Bfusion_model.M.getSolverDoubleInfo(%22optimizerTime%22)%2C%20pyomo_model.result.solver.wallclock_time%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'Pyomo'%2C%20'CVXPY'%5D%2C%20height%3D%20total_t%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20width%3D0.5%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%20'Pyomo'%2C%20'CVXPY'%5D%2C%20height%3Dsolver_t%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20width%3D0.5%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(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20%23%23%20Discussion%0A%0A%20%20%20%20Apparently%2C%20**Fusion%20passes%20the%20model%20data%20to%20Mosek%20faster**%20than%20the%20other%20ones.%20CVXPY%20is%20close%20to%20Fusion%20but%20its%20solver%20time%20is%20longer%2C%20mainly%20due%20to%20presolve%20(for%20example%20CVXPY%20enters%20variable%20bounds%20as%20constraints).%20In%20terms%20of%20total%20time%20Pyomo%20is%20behind%20the%20others%20because%20all%20the%20transformations%20are%20made%20in%20Python%2C%20as%20opposed%20to%20Fusion%20and%20CVXPY%20which%20call%20a%20C%20library.%20However%2C%20this%20is%20a%20huge%20model%20with%2031%20thousand%20constraints%20and%2012%20million%20variables%20and%20the%20difference%20will%20not%20be%20that%20big%20for%20normal-sized%20models.%0A%0A%20%20%20%20On%20the%20other%20hand%2C%20Fusion%20and%20CVXPY%20allow%20you%20to%20express%20model%20in%20vectorized%20form%20(using%20**matrix%2C%20vector**%20notation).%20Pyomo%20is%20mainly%20based%20on%20**sum%20expressions**%20that%20can%20be%20defined%20by%20**rule%20functions**%20which%20makes%20modelling%20much%20easier.%20Therefore%20the%20**time%20and%20effort%20spent%20on%20constructing%20the%20model%20in%20Pyomo%20was%20much%20smaller**%20than%20with%20%20the%20other%20languages.%0A%0A%20%20%20%20The%20decision%20of%20which%20one%20to%20use%20depends%20on%20the%20problem%20and%20the%20preferences%20of%20the%20user.%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%20CVXPY%20and%20Fusion%20on%20huge%20models%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%22The%20plots%20above%20show%20the%20performance%20of%20the%20modeling%20languages%20for%2020%20images.%20It%20is%20also%20reasonable%20to%20test%20the%20behaivor%20of%20solving%20times%20as%20the%20model%20gets%20larger%20and%20larger.%20In%20order%20to%20make%20this%20test%20same%20problem%20is%20solved%20with%2050%20images%20by%20using%20CVXPY%20and%20Fusion.%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(cvxpy_model%2C%20fusion_model%2C%20ones)%3A%0A%20%20%20%20n2%20%3D%2050%0A%20%20%20%20train_1_1%20%3D%20ones%5B%3An2%5D%0A%20%20%20%20cvxpy_model.reset()%0A%20%20%20%20fusion_model.reset()%0A%20%20%20%20return%20(train_1_1%2C)%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%20CVXPY%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(cvxpy_model%2C%20train_1_1)%3A%0A%20%20%20%20res_cvx%20%3D%20cvxpy_model.Wasserstein_BaryCenter(train_1_1)%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%20(res_cvx%2C)%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20plt%2C%20res_cvx)%3A%0A%20%20%20%20plt.style.use(%22default%22)%0A%20%20%20%20plt.figure(figsize%3D(3.3%2C3.3))%0A%20%20%20%20plt.imshow(np.reshape(res_cvx.squeeze()%2C%20(28%2C28))%2C%20cmap%3D'viridis')%0A%20%20%20%20plt.title('Barycenter')%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%23%20Fusion%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(fusion_model%2C%20train_1_1)%3A%0A%20%20%20%20res_f%20%3D%20fusion_model.Wasserstein_BaryCenter(train_1_1)%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('optimizerTime')))%0A%20%20%20%20print('The%20average%20Wasserstein%20distance%20between%20digits%20and%20the%20barycenter%3A%20%5Cn%20%7B0%7D'.format(fusion_model.objective))%0A%20%20%20%20return%20(res_f%2C)%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20plt%2C%20res_f)%3A%0A%20%20%20%20plt.figure(figsize%3D(3.3%2C3.3))%0A%20%20%20%20plt.imshow(np.reshape(res_f%2C(28%2C28))%2C%20cmap%3D'viridis')%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_model%2C%20plt)%3A%0A%20%20%20%20plt.figure(figsize%3D(8%2C3.5))%0A%0A%20%20%20%20total_t50%20%3D%20%5Bfusion_model.time%2C%20cvxpy_model.time%5D%0A%20%20%20%20solver_t50%20%3D%20%5Bfusion_model.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
57d9db3cb6a6d2a35ff5bb0dec382fa4edf2b8b5530840f6e08f244da46a6f78