import%20marimo%0A%0A__generated_with%20%3D%20%220.15.5%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%20Optimization%20of%20cycles%20on%20surfaces%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%20this%20notebook%20we%20use%20MOSEK%20Fusion%20to%20solve%20some%20geometric%20optimization%20problems%20on%20surfaces.%20The%20code%20demonstrates%20setting%20up%20a%20linear%20optimization%20problem%2C%20extending%20the%20%60Model%60%20class%20and%20shows%20the%20practical%20difference%20between%20interior-point%20and%20basic%20solution.%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%23%20Requirements%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20%23%20'%25matplotlib%20notebook'%20command%20supported%20automatically%20in%20marimo%0A%20%20%20%20import%20mosek%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%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%20numpy%20as%20np%20%0A%20%20%20%20import%20sys%2C%20io%2C%20math%2C%20plyfile%0A%20%20%20%20import%20matplotlib.pyplot%20as%20plt%0A%0A%20%20%20%20%23Fetching%20examples%20predefined%20for%20this%20presentation%0A%20%20%20%20def%20getExample(name)%3A%0A%20%20%20%20%20%20%20%20return%20plyfile.PlyData.read(open(name%2C%20'rb'))%0A%20%20%20%20return%20Domain%2C%20Matrix%2C%20Model%2C%20ObjectiveSense%2C%20getExample%2C%20math%2C%20np%2C%20plt%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20%23%23%23%23%20References%0A%20%20%20%20%3Ca%20href%3D%22https%3A%2F%2Fpypi.python.org%2Fpypi%2Fplyfile%22%3E%60plyfile%60%3C%2Fa%3E%20is%20a%20package%20for%20reading%20the%20descriptions%20of%20surfaces%20in%20%3Ca%20href%3D%22https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FPLY_(file_format)%22%3EPLY%20format%3C%2Fa%3E.%20It%20can%20be%20installed%20with%20%60pip%20install%20plyfile%60.%20You%20will%20also%20need%20some%20data%20files%20with%20examples%20used%20in%20this%20presentation.%20The%20complete%20code%20and%20additional%20files%20are%20available%20from%20%3Ca%20href%3D%22https%3A%2F%2Fgithub.com%2FMOSEK%2FTutorials%2Ftree%2Fmaster%2Fsurfacecycles%22%3EGitHub%3C%2Fa%3E.%0A%0A%20%20%20%20One%20of%20the%20models%20used%20in%20this%20tutorial%20comes%20from%20%3Ca%20href%3D%22http%3A%2F%2Fpeople.sc.fsu.edu%2F~jburkardt%2Fdata%2Fply%2Fply.html%22%3EJohn%20Burkardt's%3C%2Fa%3E%20library%20of%20PLY%20files%20and%20the%20other%20ones%20were%20generated%20using%20%3Ca%20href%3D%22http%3A%2F%2Fwww.cgal.org%2F%22%3ECGAL%3C%2Fa%3E.%20A%20comprehensive%20survey%20of%20topological%20optimization%20problems%20can%20be%20found%20in%20the%20%3Ca%20href%3D%22http%3A%2F%2Fjeffe.cs.illinois.edu%2Fpubs%2Fpdf%2Foptcycles.pdf%22%3Esurvey%20paper%3C%2Fa%3E%20by%20Jeff%20Erickson.%20The%20model%20used%20in%20this%20tutorial%20comes%20from%20the%20paper%20%3Ca%20href%3D%22https%3A%2F%2Farxiv.org%2Fabs%2F1001.0338%22%3EOptimal%20Homologous%20Cycles%2C%20Total%20Unimodularity%2C%20and%20Linear%20Programming%3C%2Fa%3E%20by%20Dey%2C%20Hirani%20and%20Krishnamoorthy.%20The%20experts%20in%20the%20field%20are%20kindly%20asked%20to%20fill%20in%20some%20vague%20statements%20in%20the%20following%20text%20with%20the%20precise%20algebro-topological%20language%20on%20their%20own.%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%20%23%23%20A%20crash-course%20in%20homotopy%0A%20%20%20%20A%20triangulated%20surface%20is%20a%20mesh%20of%202D%20triangles%20in%203D%20that%20fit%20together%20to%20form%20a%20surface.%20The%20vertices%20and%20edges%20of%20a%20triangulation%20form%20a%20graph.%20The%20*cycles*%20in%20that%20graph%20can%20be%20used%20to%20detect%20interesting%20*topological%20features*.%20For%20an%20example%20of%20what%20that%20means%2C%20see%20the%20three%20loops%20on%20the%20surface%20below%3A%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%3Cimg%20src%3D%22https%3A%2F%2Fraw.githubusercontent.com%2FMOSEK%2FTutorials%2Fmaster%2Fsurfacecycles%2Fsurface1.png%22%20alt%3D%22surface1%22%20width%3D%22400%22%3E%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%20Now%20imagine%20that%20each%20loop%20is%20a%20rubber%20band%20that%20can%20be%20stretched%20in%20a%20continuous%20way%20and%20slide%20on%20the%20surface.%20Then%20the%20violet%20loop%20is%20not%20very%20interesting%20-%20it%20lies%20in%20a%20small%20patch%20of%20the%20surface%20and%20if%20we%20start%20shrinking%20it%2C%20it%20will%20collapse%20to%20a%20point%20(in%20other%20words%2C%20it%20is%20*contractible*).%20The%20yellow%20loop%20is%20different%20-%20it%20circles%20around%20a%20hole%20in%20the%20surface%2C%20and%20no%20matter%20how%20much%20we%20deform%20it%2C%20it%20will%20not%20disappear.%20We%20say%20this%20loop%20*detects%20a%20feature*%20-%20in%20this%20case%20the%20feature%20is%20one%20of%20the%20holes.%20The%20red%20loop%20has%20the%20same%20property.%20Moreover%2C%20the%20yellow%20and%20red%20loops%20are%20*homotopic*%20-%20each%20of%20them%20can%20be%20continuosly%20moulded%20into%20the%20second%20one%20-%20i.e.%20they%20detect%20the%20same%20feature%20(hole).%0A%0A%20%20%20%20In%20some%20applications%20it%20is%20useful%20to%20have%20*small*%20representations%20of%20topological%20features%2C%20in%20this%20case%20the%20shortest%20possible%20loops%20in%20a%20given%20homotopy%20class.%20There%20are%20at%20least%20two%20reasonable%20ways%20of%20measuring%20the%20length%20of%20a%20path%20on%20a%20surface%3A%20with%20the%20Euclidean%20metric%20on%20the%20edges%2C%20or%20just%20counting%20the%20number%20of%20edges.%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_(math)%3A%0A%20%20%20%20def%20euclideanM(x%2C%20y)%3A%0A%20%20%20%20%20%20%20%20return%20math.sqrt(sum(%5B(x%5B_i%5D%20-%20y%5B_i%5D)%20**%202%20for%20_i%20in%20range(3)%5D))%0A%0A%20%20%20%20def%20graphM(x%2C%20y)%3A%0A%20%20%20%20%20%20%20%20return%201%0A%20%20%20%20return%20(euclideanM%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20%23%23%20First%20problem%3A%20shortest%20equivalent%20representation%20of%20a%20cycle%0A%20%20%20%20Suppose%20the%20surface%20has%20%24v%24%20vertices%2C%20%24e%24%20edges%20and%20%24t%24%20triangles.%20It%20is%20convenient%20to%20represent%20it%20by%20*boundary%20matrices*%2C%20which%20are%20just%20variants%20of%20the%20%3Ca%20href%3D%22https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FIncidence_matrix%22%3Eincidence%20matrix%20of%20a%20graph%3C%2Fa%3E.%20Specifically%2C%20let%20%24V(i)%24%2C%20%24E(ij)%24%20and%20%24T(ijk)%24%20denote%20the%20vertex%2C%20edge%20and%20triangle%20with%20the%20given%20vertex%20set.%20Then%20we%20have%20a%20linear%20map%20%24D_1%3A%5Cmathbb%7BR%7D%5Ee%5Cto%5Cmathbb%7BR%7D%5Ev%24%20given%20on%20the%20basis%20vectors%20by%3A%0A%0A%20%20%20%20%24%24D_1(%5Cmathbf%7Be%7D_%7BE(ij)%7D)%3D%5Cmathbf%7Be%7D_%7BV(j)%7D-%5Cmathbf%7Be%7D_%7BV(i)%7D%24%24%0A%0A%20%20%20%20and%20we%20identify%20%24D_1%5Cin%20%5Cmathbb%7BR%7D%5E%7Bv%5Ctimes%20e%7D%24%20with%20the%20matrix%20of%20this%20map.%20Similarly%2C%20we%20can%20record%20the%20incidence%20between%20triangles%20and%20edges%20by%20%24D_2%3A%5Cmathbb%7BR%7D%5Et%5Cto%5Cmathbb%7BR%7D%5Ee%24%3A%0A%0A%20%20%20%20%24%24D_2(%5Cmathbf%7Be%7D_%7BT(ijk)%7D)%3D%5Cmathbf%7Be%7D_%7BE(jk)%7D-%5Cmathbf%7Be%7D_%7BE(ik)%7D%2B%5Cmathbf%7Be%7D_%7BE(ij)%7D%24%24%0A%0A%20%20%20%20and%20let%20%24D_2%5Cin%20%5Cmathbb%7BR%7D%5E%7Be%5Ctimes%20t%7D%24%20be%20its%20matrix.%20We%20refer%20to%20%3Ca%20href%3D%22https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FSimplicial_homology%22%3Eother%20sources%3C%2Fa%3E%20for%20more%20details.%20Note%20that%20matrices%20%24D_1%2CD_2%24%20are%20extremely%20sparse%20-%20they%20have%20just%20two%2C%20resp.%20three%2C%20nonzeros%20in%20each%20column.%0A%0A%20%20%20%20%23%23%23%20Shortest%20homologous%20cycle%20problem%0A%20%20%20%20A%20*chain*%20(more%20precisely%2C%20*1-chain*)%20is%20just%20a%20vector%20%24c%5Cin%20%5Cmathbb%7BR%7D%5Ee%24%2C%20that%20is%20an%20assignment%20of%20a%20real%20coefficient%20to%20every%20edge.%20If%20%24c%24%20is%20an%20actual%20path%20or%20cycle%20in%20the%20graph%2C%20then%20these%20coefficients%20are%20%24%5Cpm1%24%20for%20edges%20in%20the%20cycle%2C%20and%20%240%24%20otherwise%2C%20but%20more%20general%20chains%20are%20allowed.%20We%20will%20mostly%20concentrate%20on%20cycles.%0A%0A%20%20%20%20We%20can%20now%20formulate%20the%20following%20problem%3A%20given%20a%20cycle%20%24c%5Cin%5Cmathbb%7BR%7D%5Ee%24%2C%20find%20the%20shortest%20cycle%20that%20detects%20the%20same%20topological%20features%20as%20%24c%24%20(strictly%20speaking%2C%20the%20shortest%20cycle%20*homologous*%20to%20%24c%24).%20Suppose%20%24w_1%2C%5Cldots%2Cw_e%24%20are%20the%20weights%20(lengths)%20of%20the%20edges.%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cbegin%7Barray%7D%7Bll%7D%0A%20%20%20%20%5Ctextrm%7Bminimize%7D%20%20%20%26%20%5Csum_%7Bi%3D1%7D%5Ee%20w_i%7Cx_i%7C%20%5C%5C%0A%20%20%20%20%5Ctextrm%7Bsubject%20to%7D%20%26%20x-c%3DD_2y%20%5C%5C%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%26%20x%5Cin%5Cmathbb%7BR%7D%5Ee%2C%20y%5Cin%5Cmathbb%7BR%7D%5Et%0A%20%20%20%20%5Cend%7Barray%7D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20%23%23%23%20Implementation%0A%20%20%20%20This%20optimization%20problem%20can%20easily%20be%20formulated%20as%20a%20linear%20program.%20Below%20is%20the%20implementation%20of%20this%20program%20in%20a%20class%20%60Surface%60%2C%20which%20extends%20a%20Mosek%20Fusion%20%60Model%60.%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%20Matrix%2C%20Model%2C%20ObjectiveSense%2C%20plt)%3A%0A%20%20%20%20class%20Surface%3A%0A%0A%20%20%20%20%20%20%20%20def%20__init__(self%2C%20plydata%2C%20metric)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20self.v%20%3D%20plydata.elements%5B0%5D.count%0A%20%20%20%20%20%20%20%20%20%20%20%20self.t%20%3D%20plydata.elements%5B1%5D.count%0A%20%20%20%20%20%20%20%20%20%20%20%20self.coo%20%3D%20plydata%5B'vertex'%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20self.xc%2C%20self.yc%2C%20self.zc%20%3D%20zip(*self.coo)%0A%20%20%20%20%20%20%20%20%20%20%20%20self.tri%20%3D%20%5Btuple(sorted(f%5B0%5D))%20for%20f%20in%20plydata%5B'face'%5D%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20self.edg%20%3D%20list(set(%5Btuple(sorted(%5Bf%5B_i%5D%2C%20f%5Bj%5D%5D))%20for%20f%20in%20self.tri%20for%20_i%2C%20j%20in%20%5B(0%2C%201)%2C%20(0%2C%202)%2C%20(1%2C%202)%5D%5D))%0A%20%20%20%20%20%20%20%20%20%20%20%20self.e%20%3D%20len(self.edg)%0A%20%20%20%20%20%20%20%20%20%20%20%20self.wgh%20%3D%20%5Bmetric(self.coo%5Bf%5B0%5D%5D%2C%20self.coo%5Bf%5B1%5D%5D)%20for%20f%20in%20self.edg%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20d1%20%3D%20%5B(self.edg%5B_i%5D%5Bk%5D%2C%20_i%2C%20sgn)%20for%20_i%20in%20range(self.e)%20for%20k%2C%20sgn%20in%20%5B(1%2C%20%2B1)%2C%20(0%2C%20-1)%5D%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20s1%2C%20s2%2C%20val%20%3D%20zip(*d1)%0A%20%20%20%20%20%20%20%20%20%20%20%20self.D1%20%3D%20Matrix.sparse(self.v%2C%20self.e%2C%20list(s1)%2C%20list(s2)%2C%20list(val))%0A%20%20%20%20%20%20%20%20%20%20%20%20d2%20%3D%20%5B(self.edg.index((self.tri%5B_i%5D%5Bk%5D%2C%20self.tri%5B_i%5D%5Bl%5D))%2C%20_i%2C%20sgn)%20for%20_i%20in%20range(self.t)%20for%20k%2C%20l%2C%20sgn%20in%20%5B(1%2C%202%2C%20%2B1)%2C%20(0%2C%202%2C%20-1)%2C%20(0%2C%201%2C%20%2B1)%5D%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20s1%2C%20s2%2C%20val%20%3D%20zip(*d2)%0A%20%20%20%20%20%20%20%20%20%20%20%20self.D2%20%3D%20Matrix.sparse(self.e%2C%20self.t%2C%20list(s1)%2C%20list(s2)%2C%20list(val))%0A%0A%20%20%20%20%20%20%20%20def%20plot(self%2C%20chains%3D%5B%5D%2C%20colors%3D%5B%5D)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20fig%20%3D%20plt.figure()%0A%20%20%20%20%20%20%20%20%20%20%20%20ax%20%3D%20plt.axes(projection%3D'3d'%2C%20computed_zorder%3DFalse)%0A%20%20%20%20%20%20%20%20%20%20%20%20ax.plot_trisurf(self.xc%2C%20self.yc%2C%20self.zc%2C%20triangles%3Dself.tri%2C%20linewidth%3D0.5%2C%20alpha%3D0.5%2C%20edgecolor%3D'green')%0A%20%20%20%20%20%20%20%20%20%20%20%20_i%20%3D%200%0A%20%20%20%20%20%20%20%20%20%20%20%20for%20C%20in%20chains%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20for%20j%20in%20range(self.e)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20if%20abs(C%5Bj%5D)%20%3E%200.0001%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%20xc%2C%20yc%2C%20zc%20%3D%20zip(*%5Bself.coo%5Bself.edg%5Bj%5D%5B0%5D%5D%2C%20self.coo%5Bself.edg%5Bj%5D%5B1%5D%5D%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20ax.plot(xc%2C%20yc%2C%20zc%2C%20color%3Dcolors%5B_i%5D%2C%20alpha%3Dmin(1%2C%20abs(C%5Bj%5D))%2C%20zorder%3D1)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_i%20%2B%3D%201%0A%20%20%20%20%20%20%20%20%20%20%20%20plt.axis('off')%0A%20%20%20%20%20%20%20%20%20%20%20%20plt.show()%0A%0A%20%20%20%20%20%20%20%20def%20short(self%2C%20C)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20with%20Model('short')%20as%20M%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20x%20%3D%20M.variable(self.e%2C%20Domain.unbounded())%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20xabs%20%3D%20M.variable(self.e%2C%20Domain.unbounded())%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20y%20%3D%20M.variable(self.t%2C%20Domain.unbounded())%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20M.constraint(xabs%20%3E%3D%20x)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20M.constraint(xabs%20%3E%3D%20-x)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20M.constraint(x%20-%20C%20%3D%3D%20self.D2%20%40%20y)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20M.objective(ObjectiveSense.Minimize%2C%20xabs.T%20%40%20self.wgh)%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%20return%20x.level()%0A%20%20%20%20return%20(Surface%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20The%20initialization%20part%20processes%20the%20PLY%20description%20of%20the%20surface%20and%20initializes%20the%20sparse%20boundary%20matrices.%20Every%20time%20we%20optimize%20for%20the%20shortest%20cycle%20a%20new%20Fusion%20model%20is%20constructed%20and%20solved%20from%20the%20input%20data.%0A%0A%20%20%20%20%23%23%23%20Example%0A%20%20%20%20We%20predefined%20the%20yellow%20cycle%20on%20the%20mug%2C%20and%20Mosek%20finds%20the%20red%20cycle%20-%20its%20shortest%20topological%20equivalent.%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%3Cimg%20src%3D%22https%3A%2F%2Fraw.githubusercontent.com%2FMOSEK%2FTutorials%2Fmaster%2Fsurfacecycles%2Fmug.png%22%20alt%3D%22mug%22%20width%3D%22400%22%3E%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%22Indeed%2C%20if%20we%20again%20imagine%20the%20yellow%20loop%20is%20made%20out%20of%20rubber%2C%20then%20once%20it%20starts%20contracting%20it%20will%20reach%20its%20minimum%20length%20when%20it%20sits%20tightly%20around%20the%20mug's%20handle.%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20%23%23%20Detecting%20features%0A%20%20%20%20OK%2C%20but%20what%20if%20we%20don't%20have%20any%20predefined%20cycle%2C%20and%20we%20would%20simply%20like%20to%20visualize%20all%20topological%20features%20on%20a%20given%20input%20surface%20using%20shortest%20possible%20loops%3F%20In%20general%2C%20this%20is%20a%20hard%20problem%2C%20for%20which%20some%20algorithms%20are%20known%20on%20surfaces.%20However%2C%20we%20can%20use%20our%20model%20to%20try%20a%20more%20pedestrian%20approach%2C%20which%20is%20not%20perfect%2C%20but%20often%20good%20enough%20for%20visualization.%0A%0A%20%20%20%20Without%20going%20into%20too%20much%20algebra%2C%20the%20space%20of%20features%20on%20the%20surface%20can%20be%20identified%20with%20the%20nullspace%20of%20the%20matrix%20%24%24A%3D%5Cleft%5B%5Cbegin%7Barray%7D%7Bl%7D%20D_1%20%5C%5C%20D_2%5ET%5Cend%7Barray%7D%5Cright%5D%24%24%20(for%20the%20experts%2C%20we%20are%20now%20computing%20first%20homology%20with%20real%20coefficients).%20This%20can%20be%20easily%20found%20using%20singular%20value%20decomposition%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_(np)%3A%0A%20%20%20%20def%20nullspace(A)%3A%0A%20%20%20%20%20%20%20%20u%2C%20s%2C%20vh%20%3D%20np.linalg.svd(A)%0A%20%20%20%20%20%20%20%20nnz%20%3D%20(s%20%3E%3D%201e-13).sum()%0A%20%20%20%20%20%20%20%20ns%20%3D%20vh%5Bnnz%3A%5D.conj().T%0A%20%20%20%20%20%20%20%20return%20(ns%2C%20ns.shape%5B1%5D)%0A%0A%20%20%20%20def%20homology(S)%3A%0A%20%20%20%20%20%20%20%20return%20nullspace(np.reshape(np.concatenate((S.D1.getDataAsArray()%2C%20S.D2.transpose().getDataAsArray()))%2C%20%5BS.v%20%2B%20S.t%2C%20S.e%5D))%0A%20%20%20%20return%20(homology%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20For%20an%20orientable%20surface%20the%20dimension%20of%20the%20nullspace%20of%20%24A%24%20equals%20%242g%24%2C%20where%20%24g%24%20is%20the%20%3Ca%20href%3D%22https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FGenus_(mathematics)%22%3Egenus%3C%2Fa%3E%2C%20or%20the%20number%20of%20holes.%20For%20example%2C%20the%20torus%20has%20one%20hole%20and%20two%20essential%20features%3A%20the%20parallel%20and%20the%20meridian%20cycle.%0A%0A%20%20%20%20Since%20the%20SVD%20decomposition%20is%20oblivious%20to%20the%20underlying%20geometry%2C%20the%20chains%20it%20returns%20are%20typically%20dense%20vectors%20with%20real%20coordinates%2C%20not%20very%20useful%20for%20visualization.%20However%2C%20we%20can%20replace%20them%20with%20their%20shortest%20equivalents%20and%20see%20if%20they%20are%20any%20better.%20It%20turns%20out%20in%20practice%20they%20often%20are.%20To%20make%20this%20easier%20to%20see%20we%20also%20normalize%20so%20that%20the%20maximal%20coefficient%20of%20each%20chain%20is%201.%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_(Surface%2C%20euclideanM%2C%20getExample%2C%20homology)%3A%0A%20%20%20%20def%20normalize(C)%3A%0A%20%20%20%20%20%20%20%20return%20C%20%2F%20max(abs(C))%0A%20%20%20%20S%20%3D%20Surface(getExample('torus.ply')%2C%20euclideanM)%0A%20%20%20%20_generators%2C%20_betti1%20%3D%20homology(S)%0A%20%20%20%20_G%20%3D%20%5Bnormalize(_generators%5B%3A%2C%20_i%5D)%20for%20_i%20in%20range(_betti1)%5D%0A%20%20%20%20GS%20%3D%20%5Bnormalize(S.short(g))%20for%20g%20in%20_G%5D%0A%20%20%20%20for%20_i%20in%20range(_betti1)%3A%0A%20%20%20%20%20%20%20%20S.plot(%5B_G%5B_i%5D%2C%20GS%5B_i%5D%5D%2C%20%5B'yellow'%2C%20'red'%5D)%0A%20%20%20%20return%20(normalize%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20The%20two%20projections%20of%20the%20torus%20do%20not%20reveal%20where%20the%20holes%20are.%20In%20yellow%20we%20see%20the%20dense%20cycles%20found%20by%20SVD.%20%20Their%20shortest%20homologous%20red%20versions%20make%20this%20a%20bit%20more%20clear.%20For%20example%2C%20the%20coil-shaped%20part%20of%20the%20first%20cycle%20clearly%20detects%20the%20hole%20in%20the%20middle%20of%20the%20torus.%0A%0A%20%20%20%20Here%20is%20an%20example%20of%20this%20on%20a%20more%20complicated%20surface.%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_(Surface%2C%20euclideanM%2C%20getExample%2C%20homology%2C%20normalize)%3A%0A%20%20%20%20def%20_()%3A%0A%20%20%20%20%20%20%20%20S%20%3D%20Surface(getExample('surface.ply')%2C%20euclideanM)%0A%20%20%20%20%20%20%20%20_generators%2C%20_betti1%20%3D%20homology(S)%0A%20%20%20%20%20%20%20%20print('Number%20of%20features%3A%20%7B0%7D'.format(_betti1))%0A%20%20%20%20%20%20%20%20_G%20%3D%20%5Bnormalize(_generators%5B%3A%2C%20_i%5D)%20for%20_i%20in%20range(_betti1)%5D%0A%20%20%20%20%20%20%20%20for%20_i%20in%20%5B2%2C%205%5D%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20S.plot(%5B_G%5B_i%5D%2C%20normalize(S.short(_G%5B_i%5D))%5D%2C%20%5B'yellow'%2C%20'red'%5D)%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).%20%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
5d876c51a91480edc329b9f5389df5d63309cbe43d39496358d691b41ddfbc0b