ellipsoid peeling

%23%20marimo%20notebook%3A%20ellipsoid%20peeling%20with%20convex%20hull%0A%23%20Run%3A%20%20marimo%20run%20app.py%0A%0Aimport%20marimo%0A%0A__generated_with%20%3D%20%220.20.4%22%0Aapp%20%3D%20marimo.App(width%3D%22full%22)%0A%0A%0A%40app.cell%0Aasync%20def%20_()%3A%0A%20%20%20%20import%20sys%0A%0A%20%20%20%20%23%20Only%20run%20this%20in%20Pyodide%2FWASM%0A%20%20%20%20if%20sys.platform%20%3D%3D%20%22emscripten%22%3A%0A%0A%20%20%20%20%20%20%20%20import%20micropip%0A%0A%20%20%20%20%20%20%20%20%23%20Minimal%20deps%20for%20this%20notebook%20in-browser%0A%20%20%20%20%20%20%20%20%23%20NOTE%3A%20cvxpy%20in%20Pyodide%20is%20published%20as%20%22cvxpy-base%22%0A%20%20%20%20%20%20%20%20await%20micropip.install(%5B%22cvxpy-base%22%2C%20%22clarabel%22%2C%20%22numpy%22%2C%20%22matplotlib%22%2C%20%22anywidget%22%2C%20%22wigglystuff%22%5D)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20import%20marimo%20as%20mo%0A%20%20%20%20import%20numpy%20as%20np%0A%20%20%20%20import%20cvxpy%20as%20cp%0A%20%20%20%20import%20pandas%20as%20pd%0A%20%20%20%20import%20matplotlib.pyplot%20as%20plt%0A%20%20%20%20from%20wigglystuff%20import%20ChartPuck%0A%0A%20%20%20%20return%20ChartPuck%2C%20cp%2C%20mo%2C%20np%0A%0A%0A%40app.cell%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%20Ellipsoid%20peeling%20using%20dual%20values%0A%20%20%20%20%23%23%20Minimum%20volume%20ellipsoid%20around%20a%20set%20of%20points%0A%20%20%20%20*%20Ellipsoid%20parametarized%20as%20%24%5Cvarepsilon%3D%5C%7B%5Cnu%20%7C%20%7C%7CA%5Cnu%20%2Bb%7C%7C_2%20%5Cleq%201%5C%7D%24.%0A%20%20%20%20*%20Finite%20set%20of%20points%20%24%5C%7Bx_1%2C%20x_2%2C%20...%2C%20x_m%5C%7D%24.%0A%20%20%20%20*%20Volume%20is%20proportional%20to%20%24%5Cdet%20A%5E%7B-1%7D%24%2C%20so%20can%20find%20%24%5Cvarepsilon%24%20by%20solving%20the%20convex%20problem%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Ctext%7Bminimize%20(over%20%7D%20A%2Cb)%20%5C%3B%5C%3B%5C%3B%5C%3B%20%5Clog%5Cdet%20A%5E%7B-1%7D%20%5C%5C%0A%20%20%20%20%5Ctext%7Bsubject%20to%20%7D%20%7C%7CAx_i%2Bb%7C%7C_2%20%5Cleq%201%2C%20%5C%3B%5C%3B%5C%3B%5C%3B%20i%3D1%2C...%2C%20m%0A%20%20%20%20%24%24%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(ChartPuck%2C%20cp%2C%20mo%2C%20np)%3A%0A%20%20%20%20n%20%3D%207%0A%0A%0A%20%20%20%20def%20draw_func(ax%2C%20widget)%3A%0A%20%20%20%20%20%20%20%20positions%20%3D%20np.array(%5B%5Bwidget.x%5Bi%5D%2C%20widget.y%5Bi%5D%5D%20for%20i%20in%20range(n)%5D)%0A%20%20%20%20%20%20%20%20A%20%3D%20cp.Variable((2%2C%202)%2C%20PSD%3DTrue)%0A%20%20%20%20%20%20%20%20b%20%3D%20cp.Variable(2)%0A%0A%20%20%20%20%20%20%20%20cons%20%3D%20%5B%5D%0A%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%20cons%20%2B%3D%20%5Bcp.norm2(A%20%40%20positions%5Bj%5D%20%2B%20b)%20%3C%3D%201%5D%0A%0A%20%20%20%20%20%20%20%20prob%20%3D%20cp.Problem(cp.Maximize(cp.log_det(A))%2C%20cons)%0A%20%20%20%20%20%20%20%20prob.solve(solver%3Dcp.CLARABEL)%0A%0A%0A%20%20%20%20%20%20%20%20Ahat%20%3D%20A.value%0A%20%20%20%20%20%20%20%20bhat%20%3D%20b.value%0A%0A%20%20%20%20%20%20%20%20if%20Ahat%20is%20None%20or%20bhat%20is%20None%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20raise%20ValueError(%22No%20values%20set%20for%20A%20or%20b%22)%0A%0A%20%20%20%20%20%20%20%20%23%20parameterize%20boundary%0A%20%20%20%20%20%20%20%20theta%20%3D%20np.linspace(0%2C%202*np.pi%2C%20400)%0A%20%20%20%20%20%20%20%20U%20%3D%20np.vstack(%5Bnp.cos(theta)%2C%20np.sin(theta)%5D)%20%20%20%20%20%20%20%20%20%20%23%20Pick%20points%20on%20unit%20circle%0A%20%20%20%20%20%20%20%20X%20%3D%20np.linalg.solve(Ahat%2C%20U%20-%20bhat.reshape(2%2C1))%20%20%20%20%20%20%20%23%20Find%20pre-image%20of%20points%0A%0A%0A%20%20%20%20%20%20%20%20ax.scatter(positions%5B%3A%2C%200%5D%2C%20positions%5B%3A%2C%201%5D%2C%20s%3D20%2C%20label%3D%22points%22)%0A%20%20%20%20%20%20%20%20ax.plot(X%5B0%2C%20%3A%5D%2C%20X%5B1%2C%20%3A%5D%2C%20linewidth%3D2%2C%20label%3D%22ellipsoid%20boundary%22)%0A%0A%20%20%20%20%20%20%20%20%23%20Duals%20calculation%0A%20%20%20%20%20%20%20%20lambdas%20%3D%20np.array(%5Bc.dual_value%20for%20c%20in%20cons%5D%2C%20dtype%3Dfloat)%20%20%23%20shape%20(n%2C)%0A%20%20%20%20%20%20%20%20lambdas%20%3D%20np.maximum(lambdas%2C%200.0)%0A%0A%20%20%20%20%20%20%20%20Y%20%3D%20(Ahat%20%40%20positions.T%20%2B%20bhat.reshape(2%2C%201))%0A%20%20%20%20%20%20%20%20Ynorm%20%3D%20np.linalg.norm(Y%2C%20axis%3D0)%20%2B%201e-12%0A%20%20%20%20%20%20%20%20N%20%3D%20(Ahat.T%20%40%20(Y%20%2F%20Ynorm))%0A%0A%20%20%20%20%20%20%20%20if%20np.max(lambdas)%20%3E%200%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20lam_scaled%20%3D%20lambdas%20%2F%20np.max(lambdas)%0A%20%20%20%20%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20lam_scaled%20%3D%20lambdas%0A%0A%20%20%20%20%20%20%20%20base%20%3D%202.0%0A%20%20%20%20%20%20%20%20Ux%20%3D%20-N%5B0%2C%20%3A%5D%20*%20(base%20*%20lam_scaled)%0A%20%20%20%20%20%20%20%20Uy%20%3D%20-N%5B1%2C%20%3A%5D%20*%20(base%20*%20lam_scaled)%0A%0A%20%20%20%20%20%20%20%20ax.quiver(%0A%20%20%20%20%20%20%20%20%20%20%20%20positions%5B%3A%2C%200%5D%2C%20positions%5B%3A%2C%201%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20Ux%2C%20Uy%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20angles%3D%22xy%22%2C%20scale_units%3D%22xy%22%2C%20scale%3D1%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20width%3D0.006%2C%20alpha%3D0.9%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20label%3D%22dual%20normal%20force%20(scaled)%22%0A%20%20%20%20%20%20%20%20)%0A%0A%20%20%20%20%20%20%20%20ax.set_xlim(-4%2C%204)%0A%20%20%20%20%20%20%20%20ax.set_ylim(-4%2C%204)%0A%20%20%20%20%20%20%20%20ax.set_xlabel(%22X%22)%0A%20%20%20%20%20%20%20%20ax.set_ylabel(%22Y%22)%0A%20%20%20%20%20%20%20%20ax.set_title(%22Drag%20points%20to%20move%20them%22)%0A%20%20%20%20%20%20%20%20ax.grid(True%2C%20alpha%3D0.3)%0A%0A%20%20%20%20multi_puck%20%3D%20ChartPuck.from_callback(%0A%20%20%20%20%20%20%20%20draw_fn%3Ddraw_func%2C%0A%20%20%20%20%20%20%20%20figsize%3D(8%2C%208)%2C%0A%20%20%20%20%20%20%20%20x_bounds%3D(-4%2C%204)%2C%0A%20%20%20%20%20%20%20%20y_bounds%3D(-4%2C%204)%2C%0A%20%20%20%20%20%20%20%20x%3D%5B-1.5%5D%20%2B%20%5B0%5D%20*%20(n-2)%20%2B%20%5B1.5%5D%2C%0A%20%20%20%20%20%20%20%20y%3D%5B0%5D%20%2B%20%5B1.5%5D%20*%20(n-1)%2C%0A%20%20%20%20%20%20%20%20puck_color%3D%22red%22%2C%0A%20%20%20%20%20%20%20%20throttle%3D100%0A%20%20%20%20)%0A%0A%20%20%20%20multi_widget%20%3D%20mo.ui.anywidget(multi_puck)%0A%20%20%20%20return%20(multi_widget%2C)%0A%0A%0A%40app.cell%0Adef%20_(multi_widget)%3A%0A%20%20%20%20multi_widget%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Observations%3A%0A%20%20%20%20*%20Dual%20corresponding%20to%20each%20point%20demonstrates%20how%20much%20the%20volume%20of%20the%20elliposid%20would%20shrink%20if%20we%20were%20to%20remove%20that%20point.%0A%20%20%20%20*%20If%20the%20ellipsoid%20were%20an%20ellastic%20membrane%2C%20we%20could%20think%20of%20this%20as%20a%20%22normal%20force%22%20acting%20on%20each%20point.%0A%20%20%20%20*%20The%20points%20not%20touching%20the%20boundary%20have%20no%20%22force%22%20exerted%20on%20them.%20Their%20corresponding%20dual%20is%20zero%20(complementary%20slackness)%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0Aif%20__name__%20%3D%3D%20%22__main__%22%3A%0A%20%20%20%20app.run()%0A