s5

import%20marimo%0A%0A__generated_with%20%3D%20%220.19.11%22%0Aapp%20%3D%20marimo.App(width%3D%22medium%22)%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20import%20marimo%20as%20mo%0A%0A%20%20%20%20return%20(mo%2C)%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20import%20numpy%20as%20np%0A%20%20%20%20rng%20%3D%20np.random.default_rng()%0A%20%20%20%20return%20np%2C%20rng%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20import%20scipy%20as%20si%0A%0A%20%20%20%20return%20(si%2C)%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20import%20matplotlib%0A%20%20%20%20import%20matplotlib.pyplot%20as%20plt%0A%20%20%20%20plt.ion()%3B%0A%20%20%20%20return%20(plt%2C)%0A%0A%0A%40app.cell%0Adef%20_(mo)%3A%0A%20%20%20%20from%20pathlib%20import%20Path%0A%20%20%20%20mo.pdf(src%3DPath(%22Homework5.pdf%22)%2C%20width%3D%22100%25%22%2C%20height%3D%2250vh%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%0A%20%20%20%20%23%23%20The%20Signal-to-Noise%20Ratio%20in%20Images%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(rng)%3A%0A%20%20%20%20xs%20%3D%20rng.poisson(lam%3D200%2C%20size%3D1000%20*%20500).reshape(1000%2C%20500)%0A%20%20%20%20xs.shape%0A%20%20%20%20return%20(xs%2C)%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20xs)%3A%0A%20%20%20%20np.mean(np.mean(xs%2C%20axis%3D0)%20%2F%20np.std(xs%2C%20axis%3D0))%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np)%3A%0A%20%20%20%20As%20%3D%20np.logspace(-6%2C%206%2C%20500)%0A%20%20%20%20As.shape%0A%20%20%20%20return%20(As%2C)%0A%0A%0A%40app.cell%0Adef%20_(As%2C%20np%2C%20xs)%3A%0A%20%20%20%20np.mean(np.mean(As%20*%20xs%2C%20axis%3D0)%20%2F%20np.std(As%20*%20xs%2C%20axis%3D0))%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(As%2C%20np%2C%20xs)%3A%0A%20%20%20%20ras%20%3D%20np.mean(xs%2C%20axis%3D0)%20%2F%20np.std(xs%2C%20axis%3D0)%20-%20np.mean(As%20*%20xs%2C%20axis%3D0)%20%2F%20np.std(As%20*%20xs%2C%20axis%3D0)%0A%20%20%20%20np.mean(ras)%0A%20%20%20%20return%20(ras%2C)%0A%0A%0A%40app.cell%0Adef%20_(plt%2C%20ras)%3A%0A%20%20%20%20fig%2C%20ax%20%3D%20plt.subplots()%0A%20%20%20%20ax.hist(ras%2C%20bins%3D%22auto%22)%0A%20%20%20%20fig%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%20Zero!%20To%20within%20typical%20machine%20error.%0A%20%20%20%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%0A%20%20%20%20%23%23%20Centroid%20Warmup%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np)%3A%0A%20%20%20%20i%20%3D%20np.arange(0%2C%201000)%0A%20%20%20%20I%20%3D%200.2%20*%20np.sqrt(i)%0A%20%20%20%20float(np.sum(i%20*%20I)%20%2F%20np.sum(I))%0A%20%20%20%20return%20(i%2C)%0A%0A%0A%40app.cell%0Adef%20_(i%2C%20np%2C%20rng)%3A%0A%20%20%20%20J%20%3D%2030%20%2F%20(i%20%2B%2020)%20%2B%200.05%20*%20rng.poisson(lam%3D10%2C%20size%3Di.size)%0A%20%20%20%20float(np.sum(i%20*%20J)%20%2F%20np.sum(J))%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%20Centroid%20Localization%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20si)%3A%0A%20%20%20%20def%20psfker(l%2C%20NA%3D0.9%2C%20side%3D1.0%2C%20N%3D210%2C%20center%3D(0.0%2C%200.0))%3A%0A%20%20%20%20%20%20%20%20xs%2C%20ys%20%3D%20np.meshgrid(np.linspace(-%20side%20%2F%202%2C%20side%20%2F%202%2C%20N)%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%20np.linspace(-%20side%20%2F%202%2C%20side%20%2F%202%2C%20N))%0A%20%20%20%20%20%20%20%20r2s%20%3D%20np.square(xs%20-%20center%5B0%5D)%20%2B%20np.square(ys%20-%20center%5B1%5D)%0A%20%20%20%20%20%20%20%20v%20%3D%202%20*%20np.pi%20*%20NA%20*%20np.sqrt(r2s)%20%2F%20l%0A%20%20%20%20%20%20%20%20psf%20%3D%204%20*%20np.square(si.special.j1(v)%20%2F%20v)%0A%20%20%20%20%20%20%20%20psf%5Bnp.isnan(psf)%5D%20%3D%201%0A%20%20%20%20%20%20%20%20scale%20%3D%20side%20%2F%20N%0A%20%20%20%20%20%20%20%20return%20psf%20%2F%20np.sum(psf)%2C%20scale%0A%0A%20%20%20%20return%20(psfker%2C)%0A%0A%0A%40app.cell%0Adef%20_(np)%3A%0A%20%20%20%20def%20pixelate(image%2C%20inscale%2C%20shape%3D(7%2C%207))%3A%0A%20%20%20%20%20%20%20%20outscale%20%3D%20image.shape%5B0%5D%20*%20inscale%20%2F%20shape%5B0%5D%0A%0A%20%20%20%20%20%20%20%20out%20%3D%20np.sum(image.reshape(shape%5B0%5D%2C%20int(image.shape%5B0%5D%20%2F%20shape%5B0%5D)%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%20shape%5B1%5D%2C%20int(image.shape%5B1%5D%20%2F%20shape%5B1%5D))%2C%20%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20axis%3D(1%2C%203))%0A%0A%20%20%20%20%20%20%20%20return%20out%2C%20outscale%0A%0A%20%20%20%20return%20(pixelate%2C)%0A%0A%0A%40app.cell%0Adef%20_(np)%3A%0A%20%20%20%20def%20fishify(image%2C%20N%3D1%2C%20rng%3Dnp.random.default_rng())%3A%0A%20%20%20%20%20%20%20%20return%20rng.poisson(lam%3DN%20*%20(image%20%2F%20np.sum(image)))%0A%0A%20%20%20%20return%20(fishify%2C)%0A%0A%0A%40app.cell%0Adef%20_(fishify%2C%20np%2C%20pixelate%2C%20psfker)%3A%0A%20%20%20%20def%20sim(l%3D0.51%2C%20NA%3D0.9%2C%20side%3D0.7%2C%20center%3D(0.0%2C%200.0)%2C%20Nphoton%3D500%2C%20Nfine%3D280%2C%20bgavg%3D10%2C%20shape%3D(7%2C%207))%3A%0A%20%20%20%20%20%20%20%20point%2C%20scale%20%3D%20psfker(l%2C%20NA%3DNA%2C%20side%3Dside%2C%20N%3DNfine%2C%20center%3Dcenter)%0A%20%20%20%20%20%20%20%20pxpt%2C%20scale%20%3D%20pixelate(point%2C%20scale%2C%20shape%3Dshape)%0A%20%20%20%20%20%20%20%20fg%20%3D%20fishify(pxpt%2C%20N%3DNphoton)%0A%20%20%20%20%20%20%20%20bg%20%3D%20fishify(np.ones(shape)%2C%20N%3Dbgavg%20*%20np.prod(shape))%0A%20%20%20%20%20%20%20%20return%20fg%20%2B%20bg%2C%20scale%0A%0A%20%20%20%20return%20(sim%2C)%0A%0A%0A%40app.cell%0Adef%20_(np)%3A%0A%20%20%20%20def%20centroid(im%2C%20scale%3DNone)%3A%0A%20%20%20%20%20%20%20%20ys%2C%20xs%20%3D%20np.indices(im.shape)%0A%20%20%20%20%20%20%20%20c%20%3D%20np.hstack(%5Bnp.sum(xs%20*%20im)%2C%20np.sum(ys%20*%20im)%5D)%20%2F%20np.sum(im)%0A%20%20%20%20%20%20%20%20c%20%3D%20(im%20*%20np.mgrid%5B0%3Aim.shape%5B0%5D%2C%200%3Aim.shape%5B1%5D%5D).sum(1).sum(1)%2Fim.sum()%0A%20%20%20%20%20%20%20%20return%20c%20if%20scale%20is%20None%20else%20scale%20*%20(c%20-%20np.array(im.shape)%20%2F%2F%202)%0A%0A%20%20%20%20return%20(centroid%2C)%0A%0A%0A%40app.cell%0Adef%20_(np)%3A%0A%20%20%20%20def%20rms(cs%2C%20center%3D(0.0%2C%200.0))%3A%0A%20%20%20%20%20%20%20%20return%20float(np.sqrt(np.mean(np.square(cs%5B%3A%2C0%5D%20-%20center%5B0%5D)%20%2B%20np.square(cs%5B%3A%2C1%5D%20-%20center%5B1%5D))))%0A%0A%20%20%20%20return%20(rms%2C)%0A%0A%0A%40app.cell%0Adef%20_(centroid%2C%20np%2C%20sim)%3A%0A%20%20%20%20sims1%20%3D%20%5Bsim(center%3D(0.0%2C%200.0))%20for%20m%20in%20range(100)%5D%0A%20%20%20%20cs1%20%3D%20np.array(%5Bcentroid(s%2C%20scale%3Dscale)%20for%20(s%2C%20scale)%20in%20sims1%5D)%0A%20%20%20%20return%20cs1%2C%20sims1%0A%0A%0A%40app.cell%0Adef%20_(centroid%2C%20cs1%2C%20plt%2C%20sims1)%3A%0A%20%20%20%20fig1%2C%20(ax1a%2C%20ax1b)%20%3D%20plt.subplots(1%2C%202%2C%20figsize%3D(8%2C%204))%0A%20%20%20%20ax1a.imshow(sims1%5B0%5D%5B0%5D%2C%20cmap%3D%22gray%22)%0A%20%20%20%20cx%2C%20cy%20%3D%20centroid(sims1%5B0%5D%5B0%5D)%0A%20%20%20%20ax1a.vlines(cx%2C%200%2C%206%2C%20color%3D%22red%22)%0A%20%20%20%20ax1a.hlines(cy%2C%200%2C%206%2C%20color%3D%22red%22)%0A%20%20%20%20ax1b.hist(cs1%5B%3A%2C0%5D%2C%20bins%3D%22auto%22)%0A%20%20%20%20fig1.tight_layout()%0A%20%20%20%20fig1%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(cs1%2C%20rms)%3A%0A%20%20%20%20rms(cs1)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(centroid%2C%20np%2C%20rms%2C%20sim)%3A%0A%20%20%20%20rmss%20%3D%20%5B%5D%0A%20%20%20%20for%20n%20in%20np.logspace(0%2C%205%2C%2015)%3A%0A%20%20%20%20%20%20%20%20sims%20%3D%20%5Bsim(Nphoton%3Dn)%20for%20m%20in%20range(100)%5D%0A%20%20%20%20%20%20%20%20cs%20%3D%20np.array(%5Bcentroid(s%2C%20scale%3Dsc)%20for%20(s%2C%20sc)%20in%20sims%5D)%0A%20%20%20%20%20%20%20%20rmss.append(rms(cs))%0A%20%20%20%20return%20(rmss%2C)%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20plt%2C%20rmss)%3A%0A%20%20%20%20fig2%2C%20ax2%20%3D%20plt.subplots()%0A%20%20%20%20ax2.scatter(np.logspace(0%2C%205%2C%2015)%2C%20rmss)%0A%20%20%20%20ax2.loglog(np.logspace(0%2C%205%2C%2010)%2C%201%20%2F%20np.sqrt(np.logspace(0%2C%205%2C%2010))%2C%20color%3D%22orange%22)%0A%20%20%20%20ax2.set_aspect('equal'%2C%20'box')%0A%20%20%20%20fig2.tight_layout()%0A%20%20%20%20fig2%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(centroid%2C%20np%2C%20sim)%3A%0A%20%20%20%20sims3a%20%3D%20%5Bsim(Nphoton%3D1000)%20for%20m%20in%20range(100)%5D%0A%20%20%20%20cs3a%20%3D%20np.array(%5Bcentroid(s%2C%20scale%3Dsc)%20for%20(s%2C%20sc)%20in%20sims3a%5D)%0A%20%20%20%20sims3b%20%3D%20%5Bsim(Nphoton%3D1000%2C%20center%3D(0.3%2C%200.0))%20for%20m%20in%20range(100)%5D%0A%20%20%20%20cs3b%20%3D%20np.array(%5Bcentroid(s%2C%20scale%3Dsc)%20for%20(s%2C%20sc)%20in%20sims3b%5D)%0A%20%20%20%20sims3c%20%3D%20%5Bsim(Nphoton%3D1000%2C%20center%3D(-0.3%2C%200.0))%20for%20m%20in%20range(100)%5D%0A%20%20%20%20cs3c%20%3D%20np.array(%5Bcentroid(s%2C%20scale%3Dsc)%20for%20(s%2C%20sc)%20in%20sims3c%5D)%0A%20%20%20%20return%20cs3a%2C%20cs3b%2C%20cs3c%0A%0A%0A%40app.cell%0Adef%20_(cs3a%2C%20cs3b%2C%20cs3c%2C%20plt)%3A%0A%20%20%20%20fig3%2C%20(ax3a%2C%20ax3b%2C%20ax3c)%20%3D%20plt.subplots(1%2C%203%2C%20figsize%3D(9%2C%203))%0A%20%20%20%20ax3a.hist(cs3a%5B%3A%2C%200%5D%20-%200.0%2C%20bins%3D%22auto%22)%0A%20%20%20%20ax3b.hist(cs3b%5B%3A%2C%200%5D%20-%200.3%2C%20bins%3D%22auto%22)%0A%20%20%20%20ax3c.hist(cs3c%5B%3A%2C%200%5D%20-%20(-0.3)%2C%20bins%3D%22auto%22)%0A%20%20%20%20fig3.tight_layout()%0A%20%20%20%20fig3%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(centroid%2C%20np%2C%20sim)%3A%0A%20%20%20%20Dxs%20%3D%20%5B%5D%0A%20%20%20%20for%20p%20in%20np.linspace(-7%2C%207%2C%2010)%3A%0A%20%20%20%20%20%20%20sms%20%3D%20%5Bsim(Nphoton%3D1000)%20for%20m%20in%20range(100)%5D%0A%20%20%20%20%20%20%20Dxs.append(np.mean(%5Bcentroid(s%2C%20scale%3Dsc)%5B0%5D%20-%20p%20for%20(s%2C%20sc)%20in%20sms%5D))%0A%20%20%20%20return%20(Dxs%2C)%0A%0A%0A%40app.cell%0Adef%20_(Dxs%2C%20np%2C%20plt)%3A%0A%20%20%20%20fig4%2C%20ax4%20%3D%20plt.subplots()%0A%20%20%20%20ax4.scatter(np.linspace(-0.5%2C%200.5%2C%2010)%2C%20Dxs)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_()%3A%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