# Public copy: internal compute identifiers and local environment paths were removed or generalized. Scientific values are retained. Paths may need adapting and some original run inputs are not included. """Animated explainer: why a geophysical inversion is non-unique. Data provenance --------------- All numbers, the convergence curves, the recovered-model sections, and the observed-TMI profile are taken from the real SimPEG run harvested from the compute environment (anim_assets/curves.json + the baked section PNGs). The buried-body cartoon in part 1 is a SCHEMATIC: the polygons are drawn by hand to introduce the geometry, and are not the true model array. Everything that carries a number on screen is measured. """ import json import numpy as np import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt from manim import * config.background_color = "#161616" BLUE_3B1B = "#58C4DD" YELLOW_3B1B = "#FFFF00" RED_3B1B = "#E86A6A" GREEN_3B1B = "#39FF88" GREY_3B1B = "#888888" D = json.load(open("anim_assets/curves.json")) ND = D["nD"] FIN_L2 = D["final"]["L2"] FIN_SP = D["final"]["sparse"] def tex_png(expr, path, color="white", fs=64, dpi=300): """Math expression -> transparent PNG via matplotlib mathtext (no LaTeX).""" fig = plt.figure(figsize=(0.01, 0.01)) fig.patch.set_alpha(0) fig.text(0, 0, f"${expr}$", color=color, fontsize=fs) fig.savefig(path, dpi=dpi, transparent=True, bbox_inches="tight", pad_inches=0.05) plt.close(fig) return path class NonUniqueness(Scene): def construct(self): self.part1_setup() self.part2_forward() self.part3_misfit() self.part4_two_models() self.part5_closing() # ------------------------------------------------------------------ part 1 def part1_setup(self): title = Text("Geophysical inversion", font_size=52, color=BLUE_3B1B) sub = Text("recovering the subsurface from surface data", font_size=25, color=GREY_3B1B).next_to(title, DOWN, buff=0.3) self.play(Write(title, run_time=1.4)) self.play(FadeIn(sub, shift=UP * 0.25), run_time=0.9) self.wait(1.0) self.play(FadeOut(title), FadeOut(sub), run_time=0.6) # ground line + two buried bodies (schematic, matching the real model) ground = Line(LEFT * 5.9, RIGHT * 5.9, color=GREY_3B1B).shift(UP * 1.15) gl = Text("surface", font_size=17, color=GREY_3B1B)\ .next_to(ground.get_left(), UP, buff=0.09).shift(RIGHT * 0.35) body1 = Polygon( [-3.7, -0.25, 0], [-1.6, 0.55, 0], [-1.1, -0.05, 0], [-3.2, -0.85, 0], color=RED_3B1B, fill_opacity=0.85, stroke_width=0, ) body2 = Polygon( [1.2, 0.25, 0], [2.8, -0.65, 0], [3.3, -0.10, 0], [1.7, 0.75, 0], color="#F0A0A0", fill_opacity=0.75, stroke_width=0, ) b1l = Text("0.06 SI", font_size=17, color=RED_3B1B).next_to(body1, DOWN, buff=0.12) b2l = Text("0.035 SI", font_size=17, color="#F0A0A0").next_to(body2, DOWN, buff=0.12) self.play(Create(ground), FadeIn(gl), run_time=0.9) self.play(FadeIn(body1, scale=0.85), FadeIn(body2, scale=0.85), run_time=1.1) self.play(FadeIn(b1l), FadeIn(b2l), run_time=0.6) cap = Text("Two magnetized bodies we cannot see", font_size=24, color=YELLOW_3B1B).to_edge(DOWN, buff=0.5) self.play(FadeIn(cap, shift=UP * 0.2), run_time=0.8) self.wait(1.2) self.play(FadeOut(cap), run_time=0.5) self.ground, self.body1, self.body2 = ground, body1, body2 self.gl = gl self.blabels = VGroup(b1l, b2l) # ------------------------------------------------------------------ part 2 def part2_forward(self): # flight line of sensors just above the surface sensors = VGroup(*[ Triangle(color=BLUE_3B1B, fill_opacity=1).scale(0.052).rotate(PI) .move_to([x, 1.72, 0]) for x in np.linspace(-5.2, 5.2, 21) ]) slab = Text("441 airborne magnetometer stations", font_size=21, color=BLUE_3B1B).to_edge(DOWN, buff=0.5) self.play(LaggedStart(*[FadeIn(s, shift=DOWN * 0.2) for s in sensors], lag_ratio=0.04), run_time=1.6) self.play(FadeIn(slab), run_time=0.6) self.wait(0.5) # the measured anomaly curve, drawn as the field is "sensed" # REAL measured TMI: the observed survey line through the peak anomaly # (dobs.npy, 21 stations). Plotted as the actual data points joined by # line segments -- no fitted or invented shape. px = np.array(D["tmi_profile"]["x_m"], dtype=float) pv = np.array(D["tmi_profile"]["tmi_nT"], dtype=float) ax = Axes(x_range=[px.min(), px.max(), 250], y_range=[-60, 420, 200], x_length=10.4, y_length=1.45, axis_config={"color": GREY_3B1B, "stroke_width": 2, "include_ticks": False}, tips=False).shift(UP * 2.95) curve = VMobject(color=YELLOW_3B1B, stroke_width=4) curve.set_points_as_corners([ax.c2p(a, b) for a, b in zip(px, pv)]) dlab = Text("observed TMI", font_size=19, color=YELLOW_3B1B)\ .to_corner(UL, buff=0.28) rng = Text(f"{pv.min():.0f} to {pv.max():.0f} nT along this line", font_size=15, color=GREY_3B1B)\ .next_to(dlab, DOWN, buff=0.08).align_to(dlab, LEFT) self.play(Create(curve, run_time=2.2, rate_func=linear)) self.play(FadeIn(dlab), FadeIn(rng), run_time=0.6) self.wait(0.8) fwd = Text("Forward problem: model \u2192 data (one answer)", font_size=25, color=BLUE_3B1B).to_edge(DOWN, buff=0.5) self.play(FadeOut(slab), FadeIn(fwd, shift=UP * 0.2), run_time=0.8) arr = Arrow(self.body1.get_top() + UP * 0.05, ax.c2p(px[5], pv[5]), color=BLUE_3B1B, stroke_width=3, buff=0.18, max_tip_length_to_length_ratio=0.05) self.play(GrowArrow(arr), run_time=0.9) self.wait(1.3) inv = Text("Inverse problem: data \u2192 model (many answers)", font_size=25, color=RED_3B1B).to_edge(DOWN, buff=0.5) arr2 = Arrow(ax.c2p(px[5], pv[5]), self.body1.get_top() + UP * 0.05, color=RED_3B1B, stroke_width=3, buff=0.18, max_tip_length_to_length_ratio=0.05) self.play(ReplacementTransform(fwd, inv), ReplacementTransform(arr, arr2), run_time=1.2) self.wait(1.4) self.play(FadeOut(VGroup(sensors, arr2, inv, dlab, rng)), FadeOut(self.blabels), run_time=0.8) self.dcurve, self.dax = curve, ax # ------------------------------------------------------------------ part 3 def part3_misfit(self): schematic = VGroup(self.ground, self.gl, self.body1, self.body2) self.play(FadeOut(schematic), FadeOut(self.dcurve), FadeOut(self.dax), run_time=0.8) head = Text("How well does a model fit the data?", font_size=33, color=BLUE_3B1B).to_edge(UP, buff=0.55) self.play(Write(head), run_time=1.1) tex_png(r"\phi_d = \left\| W_d \,(d_{pred} - d_{obs}) \right\|^2", "anim_assets/eq_phid.png", color="#FFFFFF", fs=54) eq = ImageMobject("anim_assets/eq_phid.png") eq.set(width=6.4).shift(UP * 1.15) self.play(FadeIn(eq, shift=UP * 0.25), run_time=1.0) self.wait(0.9) note = Text("Fit the data to its noise level, and \u03c6d lands near the data count", font_size=21, color=GREY_3B1B).shift(DOWN * 0.25) self.play(FadeIn(note), run_time=0.8) tex_png(r"\phi_d^{target} = n_D = %d" % ND, "anim_assets/eq_target.png", color=YELLOW_3B1B, fs=54) tgt = ImageMobject("anim_assets/eq_target.png") tgt.set(width=3.6).next_to(note, DOWN, buff=0.55) self.play(FadeIn(tgt, scale=1.06), run_time=1.0) self.wait(1.6) self.play(FadeOut(eq), FadeOut(note), FadeOut(tgt), FadeOut(head), run_time=0.8) # real convergence curves head2 = Text("Both inversions descend to that target", font_size=31, color=BLUE_3B1B).to_edge(UP, buff=0.4) self.play(Write(head2), run_time=1.0) itL, pdL = D["L2"]["it"], D["L2"]["phi_d"] itS, pdS = D["sparse"]["it"], D["sparse"]["phi_d"] # NOTE: do NOT pass include_numbers=True — Axes renders numbers through # DecimalNumber, which is LaTeX-backed and fails on this platform. # Tick labels are placed manually as Text below. ax = Axes(x_range=[0, 36, 6], y_range=[2, 6, 1], x_length=8.6, y_length=3.7, axis_config={"color": GREY_3B1B, "font_size": 18}, tips=False).shift(DOWN * 0.55 + RIGHT * 0.35) xticks = VGroup(*[ Text(str(v), font_size=15, color=GREY_3B1B) .next_to(ax.c2p(v, 2), DOWN, buff=0.14) for v in (0, 6, 12, 18, 24, 30, 36) ]) yticks = VGroup(*[ Text(f"10{'\u2070\u00b9\u00b2\u00b3\u2074\u2075\u2076'[v]}", font_size=15, color=GREY_3B1B) .next_to(ax.c2p(0, v), LEFT, buff=0.14) for v in (2, 3, 4, 5, 6) ]) xl = Text("iteration", font_size=19, color=GREY_3B1B).next_to(xticks, DOWN, buff=0.10) yl = Text("\u03c6d", font_size=19, color=GREY_3B1B)\ .rotate(PI / 2).next_to(yticks, LEFT, buff=0.12) self.play(Create(ax), FadeIn(xticks), FadeIn(yticks), FadeIn(xl), FadeIn(yl), run_time=1.2) tline = DashedLine(ax.c2p(0, np.log10(ND)), ax.c2p(36, np.log10(ND)), color=YELLOW_3B1B, stroke_width=3, dash_length=0.12) tl_lab = Text(f"target \u03c6d = {ND}", font_size=18, color=YELLOW_3B1B)\ .next_to(ax.c2p(24, np.log10(ND)), UP, buff=0.12) self.play(Create(tline), FadeIn(tl_lab), run_time=1.0) def poly(its, pds, color): pts = [ax.c2p(i, np.log10(v)) for i, v in zip(its, pds)] return VMobject(color=color, stroke_width=4).set_points_as_corners(pts) cL = poly(itL, pdL, BLUE_3B1B) cS = poly(itS, pdS, RED_3B1B) lL = Text("smooth (L2)", font_size=19, color=BLUE_3B1B)\ .move_to(ax.c2p(15, 5.35)) lS = Text("compact (sparse L0)", font_size=19, color=RED_3B1B)\ .move_to(ax.c2p(21, 4.65)) self.play(Create(cL, run_time=2.0, rate_func=linear), FadeIn(lL)) self.play(Create(cS, run_time=2.4, rate_func=linear), FadeIn(lS)) self.wait(0.6) chis = Text(f"\u03c7 = {FIN_L2['achieved_chi']:.2f} and " f"\u03c7 = {FIN_SP['achieved_chi']:.2f} \u2014 both fit the data", font_size=23, color=GREEN_3B1B).to_edge(DOWN, buff=0.16) self.play(FadeIn(chis, shift=UP * 0.2), run_time=0.9) self.wait(1.8) self.play(FadeOut(VGroup(ax, xticks, yticks, xl, yl, tline, tl_lab, cL, cS, lL, lS, chis, head2)), run_time=0.9) # ------------------------------------------------------------------ part 4 def part4_two_models(self): head = Text("Same data. Same fit. Different Earth.", font_size=36, color=YELLOW_3B1B).to_edge(UP, buff=0.45) self.play(Write(head), run_time=1.3) self.wait(0.4) secL = ImageMobject("anim_assets/sec_L2.png").scale(0.92) secS = ImageMobject("anim_assets/sec_sparse.png").scale(0.92) secL.shift(LEFT * 3.35 + DOWN * 0.35) secS.shift(RIGHT * 3.35 + DOWN * 0.35) fL = SurroundingRectangle(secL, color=BLUE_3B1B, buff=0.0, stroke_width=2) fS = SurroundingRectangle(secS, color=RED_3B1B, buff=0.0, stroke_width=2) tL = Text("smooth (L2)", font_size=22, color=BLUE_3B1B).next_to(fL, UP, buff=0.14) tS = Text("compact (sparse L0)", font_size=22, color=RED_3B1B).next_to(fS, UP, buff=0.14) self.play(FadeIn(secL, shift=RIGHT * 0.3), Create(fL), FadeIn(tL), run_time=1.1) self.play(FadeIn(secS, shift=LEFT * 0.3), Create(fS), FadeIn(tS), run_time=1.1) gl = Text("dashed green = true body outline", font_size=18, color=GREEN_3B1B)\ .to_edge(DOWN, buff=0.62) self.play(FadeIn(gl), run_time=0.7) self.wait(1.6) # contrast the two statistics that differ s1 = Text(f"peak {FIN_L2['chi_max_SI']:.3f} SI", font_size=20, color=BLUE_3B1B)\ .next_to(fL, DOWN, buff=0.14) s2 = Text(f"peak {FIN_SP['chi_max_SI']:.3f} SI", font_size=20, color=RED_3B1B)\ .next_to(fS, DOWN, buff=0.14) self.play(FadeIn(s1), FadeIn(s2), run_time=0.8) self.wait(1.0) self.play(Indicate(s1, color=YELLOW_3B1B, scale_factor=1.15), run_time=0.9) self.play(Indicate(s2, color=YELLOW_3B1B, scale_factor=1.15), run_time=0.9) self.wait(1.2) punch = Text("The regularization chose the shape \u2014 not the data", font_size=27, color=YELLOW_3B1B).to_edge(DOWN, buff=0.2) self.play(FadeOut(gl), FadeIn(punch, shift=UP * 0.25), run_time=1.0) self.wait(2.0) self.play(FadeOut(VGroup(fL, fS, tL, tS, s1, s2, punch, head)), FadeOut(secL), FadeOut(secS), run_time=1.0) # ------------------------------------------------------------------ part 5 def part5_closing(self): l1 = Text("An inversion returns one of many models", font_size=32, color="#FFFFFF") l2 = Text("that fit the data equally well.", font_size=32, color="#FFFFFF") grp = VGroup(l1, l2).arrange(DOWN, buff=0.28).shift(UP * 0.75) self.play(Write(l1, run_time=1.2)) self.play(Write(l2, run_time=1.1)) self.wait(1.0) l3 = Text("Always state the regularization with the result.", font_size=27, color=YELLOW_3B1B).next_to(grp, DOWN, buff=0.6) self.play(FadeIn(l3, shift=UP * 0.25), run_time=1.0) self.wait(1.6) cred = Text( f"SimPEG {'0.24.0'} \u00b7 {D['mesh']['n_active']:,} active octree cells" f" \u00b7 {ND} TMI data", font_size=19, color=GREY_3B1B).to_edge(DOWN, buff=0.42) self.play(FadeIn(cred), run_time=0.9) self.wait(2.2) self.play(FadeOut(VGroup(grp, l3, cred)), run_time=1.0)