da l63 tutorial

%23%20%2F%2F%2F%20script%0A%23%20requires-python%20%3D%20%22%3E%3D3.12%22%0A%23%20dependencies%20%3D%20%5B%0A%23%20%20%20%20%20%22marimo%3E%3D0.23.13%22%2C%0A%23%20%5D%0A%23%20%2F%2F%2F%0A%0Aimport%20marimo%0A%0A__generated_with%20%3D%20%220.23.13%22%0Aapp%20%3D%20marimo.App(width%3D%22medium%22)%0A%0A%0A%40app.cell(hide_code%3DTrue)%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(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%20Data%20Assimilation%3A%20From%20Chaos%20to%20Analysis%0A%20%20%20%20%23%23%203DVAR%20%C2%B7%204DVAR%20%C2%B7%20Ensemble%20Kalman%20Filter%20on%20the%20Lorenz%2063%20System%0A%0A%20%20%20%20**FERS%20Summer%20School%20%E2%80%94%20AI%20and%20Machine%20Learning%20for%20Earth%20System%20Modeling%20and%20Prediction**%0A%0A%20%20%20%20---%0A%0A%20%20%20%20%3E%20*%22All%20models%20are%20wrong%2C%20but%20some%20are%20useful%20%E2%80%94%20and%20with%20data%20assimilation%2C%20we%20can%20make%20them%20more%20useful.%22*%0A%0A%20%20%20%20Data%20assimilation%20(DA)%20is%20the%20discipline%20of%20**combining%20a%20dynamical%20model%20with%20noisy%20observations**%0A%20%20%20%20to%20obtain%20the%20best%20possible%20estimate%20of%20the%20true%20state%20of%20a%20system.%0A%20%20%20%20It%20sits%20at%20the%20heart%20of%20numerical%20weather%20prediction%20(NWP)%2C%20ocean%20reanalysis%2C%20and%20increasingly%2C%0A%20%20%20%20hybrid%20AI%2Fphysics%20forecasting%20systems.%0A%0A%20%20%20%20In%20this%20tutorial%20we%20explore%20three%20foundational%20DA%20algorithms%20using%20the%20**Lorenz%2063%20system**%20as%0A%20%20%20%20a%20low-dimensional%20but%20chaotic%20testbed%3A%0A%0A%20%20%20%20%7C%20Algorithm%20%7C%20Temporal%20scope%20%7C%20Background%20covariance%20%7C%0A%20%20%20%20%7C-----------%7C---------------%7C----------------------%7C%0A%20%20%20%20%7C%20**3DVAR**%20%7C%20Single%20analysis%20time%20%7C%20Static%20%24%5Cmathbf%7BB%7D%24%20%7C%0A%20%20%20%20%7C%20**4DVAR**%20%7C%20Time%20window%20%24%5Bt_0%2C%20t_N%5D%24%20%7C%20Static%20%24%5Cmathbf%7BB%7D%24%2C%20implicit%20flow-dependence%20%7C%0A%20%20%20%20%7C%20**EnKF**%20%20%7C%20Sequential%20%2F%20cycling%20%7C%20Flow-dependent%20from%20ensemble%20%7C%0A%0A%20%20%20%20%23%23%23%20Prerequisites%20assumed%0A%20%20%20%20-%20Familiarity%20with%20ODEs%20and%20chaos%20(you%20know%20what%20a%20strange%20attractor%20is)%0A%20%20%20%20-%20Basic%20linear%20algebra%20(matrix%20inversion%2C%20eigenvalues)%0A%20%20%20%20-%20Some%20exposure%20to%20Bayesian%20reasoning%20or%20least-squares%20estimation%0A%20%20%20%20-%20Python%20%2F%20NumPy%20comfort%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---%0A%20%20%20%20%23%23%201%20%20The%20Lorenz%2063%20System%0A%0A%20%20%20%20Edward%20Lorenz%20(1963)%20derived%20a%20minimal%203-variable%20ODE%20system%20from%20the%20Navier%E2%80%93Stokes%0A%20%20%20%20equations%20for%20Rayleigh%E2%80%93B%C3%A9nard%20convection%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cfrac%7Bdx%7D%7Bdt%7D%20%3D%20%5Csigma(y%20-%20x)%0A%20%20%20%20%5Cqquad%0A%20%20%20%20%5Cfrac%7Bdy%7D%7Bdt%7D%20%3D%20x(%5Crho%20-%20z)%20-%20y%0A%20%20%20%20%5Cqquad%0A%20%20%20%20%5Cfrac%7Bdz%7D%7Bdt%7D%20%3D%20xy%20-%20%5Cbeta%20z%0A%20%20%20%20%24%24%0A%0A%20%20%20%20With%20the%20classic%20parameters%20%24%5Csigma%3D10%2C%5C%3B%20%5Crho%3D28%2C%5C%3B%20%5Cbeta%3D8%2F3%24%20the%20system%20is%20**chaotic**%3A%0A%20%20%20%20trajectories%20starting%20from%20nearly%20identical%20initial%20conditions%20diverge%20exponentially%2C%0A%20%20%20%20with%20a%20Lyapunov%20time%20%24%5Ctau_%5CLambda%20%5Capprox%201.1%24%20time%20units.%0A%0A%20%20%20%20This%20makes%20L63%20an%20ideal%20DA%20testbed%20%E2%80%94%20the%20chaotic%20divergence%20*is*%20the%20problem%20DA%20is%20trying%20to%20solve.%0A%20%20%20%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%20numpy%20as%20np%0A%20%20%20%20from%20scipy.integrate%20import%20solve_ivp%0A%20%20%20%20from%20scipy.optimize%20import%20minimize%0A%20%20%20%20import%20matplotlib.pyplot%20as%20plt%0A%20%20%20%20import%20matplotlib.gridspec%20as%20gridspec%0A%20%20%20%20from%20matplotlib.colors%20import%20Normalize%0A%20%20%20%20import%20warnings%0A%20%20%20%20warnings.filterwarnings(%22ignore%22)%0A%0A%20%20%20%20%23%20reproducibility%0A%20%20%20%20rng%20%3D%20np.random.default_rng(42)%0A%20%20%20%20return%20minimize%2C%20np%2C%20plt%2C%20rng%2C%20solve_ivp%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20solve_ivp)%3A%0A%20%20%20%20def%20lorenz63(t%2C%20state%2C%20sigma%3D10.0%2C%20rho%3D28.0%2C%20beta%3D8.0%20%2F%203.0)%3A%0A%20%20%20%20%20%20%20%20x%2C%20y%2C%20z%20%3D%20state%0A%20%20%20%20%20%20%20%20return%20%5B%0A%20%20%20%20%20%20%20%20%20%20%20%20sigma%20*%20(y%20-%20x)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20x%20*%20(rho%20-%20z)%20-%20y%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20x%20*%20y%20-%20beta%20*%20z%2C%0A%20%20%20%20%20%20%20%20%5D%0A%0A%20%20%20%20def%20integrate_l63(%0A%20%20%20%20%20%20%20%20x0%2C%20t_span%2C%20t_eval%3DNone%2C%20sigma%3D10.0%2C%20rho%3D28.0%2C%20beta%3D8.0%20%2F%203.0%0A%20%20%20%20)%3A%0A%20%20%20%20%20%20%20%20%23%20Clip%20t_eval%20to%20t_span%20to%20guard%20against%20floating-point%20overshoot%0A%20%20%20%20%20%20%20%20if%20t_eval%20is%20not%20None%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20mask%20%20%20%3D%20(t_eval%20%3E%3D%20t_span%5B0%5D%20-%201e-10)%20%26%20(t_eval%20%3C%3D%20t_span%5B1%5D%20%2B%201e-10)%0A%20%20%20%20%20%20%20%20%20%20%20%20t_eval%20%3D%20np.clip(t_eval%5Bmask%5D%2C%20t_span%5B0%5D%2C%20t_span%5B1%5D)%0A%20%20%20%20%20%20%20%20sol%20%3D%20solve_ivp(%0A%20%20%20%20%20%20%20%20%20%20%20%20lorenz63%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20t_span%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20x0%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20t_eval%3Dt_eval%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20args%3D(sigma%2C%20rho%2C%20beta)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20method%3D%22RK45%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20rtol%3D1e-9%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20atol%3D1e-12%2C%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20return%20sol.y.T%20%20%23%20shape%20(n_steps%2C%203)%0A%0A%20%20%20%20return%20(integrate_l63%2C)%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%23%201.1%20%20Visualising%20the%20attractor%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(integrate_l63%2C%20np%2C%20plt)%3A%0A%20%20%20%20%23%20Spin-up%20trajectory%20to%20land%20on%20the%20attractor%2C%20then%20visualise%0A%20%20%20%20_t_spinup%20%3D%20np.linspace(0%2C%2010%2C%205000)%0A%20%20%20%20_x0_spinup%20%3D%20np.array(%5B1.0%2C%200.0%2C%200.0%5D)%0A%20%20%20%20_spinup%20%3D%20integrate_l63(_x0_spinup%2C%20(0%2C%2010)%2C%20_t_spinup)%0A%0A%20%20%20%20_t_attr%20%3D%20np.linspace(0%2C%2030%2C%2015000)%0A%20%20%20%20_traj%20%3D%20integrate_l63(_spinup%5B-1%5D%2C%20(0%2C%2030)%2C%20_t_attr)%0A%0A%20%20%20%20_fig%20%3D%20plt.figure(figsize%3D(13%2C%204))%0A%20%20%20%20_ax1%20%3D%20_fig.add_subplot(131%2C%20projection%3D%223d%22)%0A%20%20%20%20_ax1.plot(*_traj.T%2C%20lw%3D0.4%2C%20alpha%3D0.7%2C%20color%3D%22steelblue%22)%0A%20%20%20%20_ax1.set_title(%22Lorenz%20Attractor%20(3D)%22)%0A%20%20%20%20_ax1.set_xlabel(%22x%22)%3B%20_ax1.set_ylabel(%22y%22)%3B%20_ax1.set_zlabel(%22z%22)%0A%20%20%20%20_ax1.tick_params(labelsize%3D7)%0A%0A%20%20%20%20_ax2%20%3D%20_fig.add_subplot(132)%0A%20%20%20%20for%20_i%2C%20(_lbl%2C%20_col)%20in%20enumerate(%0A%20%20%20%20%20%20%20%20zip(%5B%22x%22%2C%20%22y%22%2C%20%22z%22%5D%2C%20%5B%22tab%3Ablue%22%2C%20%22tab%3Aorange%22%2C%20%22tab%3Agreen%22%5D)%0A%20%20%20%20)%3A%0A%20%20%20%20%20%20%20%20_ax2.plot(_t_attr%2C%20_traj%5B%3A%2C%20_i%5D%2C%20lw%3D0.6%2C%20label%3D_lbl%2C%20color%3D_col)%0A%20%20%20%20_ax2.set_xlabel(%22time%22)%3B%20_ax2.set_title(%22Time%20series%22)%3B%20_ax2.legend(fontsize%3D9)%0A%0A%20%20%20%20_ax3%20%3D%20_fig.add_subplot(133)%0A%20%20%20%20_ax3.scatter(%0A%20%20%20%20%20%20%20%20_traj%5B%3A%3A10%2C%200%5D%2C%20_traj%5B%3A%3A10%2C%202%5D%2C%0A%20%20%20%20%20%20%20%20s%3D0.5%2C%20alpha%3D0.4%2C%20c%3D_t_attr%5B%3A%3A10%5D%2C%20cmap%3D%22viridis%22%2C%0A%20%20%20%20)%0A%20%20%20%20_ax3.set_xlabel(%22x%22)%3B%20_ax3.set_ylabel(%22z%22)%3B%20_ax3.set_title(%22x%E2%80%93z%20projection%22)%0A%0A%20%20%20%20_fig.tight_layout()%0A%20%20%20%20_fig%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---%0A%20%20%20%20%23%23%202%20%20The%20Twin-Experiment%20Framework%0A%0A%20%20%20%20We%20work%20in%20a%20**perfect-model%2C%20twin-experiment**%20setup%3A%0A%0A%20%20%20%201.%20**Nature%20run**%20(truth)%3A%20integrate%20L63%20from%20a%20known%20IC%20for%20time%20%24T%24.%0A%20%20%20%202.%20**Observations**%3A%20sample%20the%20truth%20at%20every%20%24%5CDelta%20t_%7B%5Crm%20obs%7D%24%20and%20add%20Gaussian%20noise%20%24%5Cmathcal%7BN%7D(0%2C%20%5Cmathbf%7BR%7D)%24.%0A%20%20%20%203.%20**Background**%20%24%5Cmathbf%7Bx%7D%5Eb%24%3A%20a%20perturbed%20IC%20that%20represents%20our%20prior%20%2F%20first-guess.%0A%20%20%20%204.%20**Analysis**%20%24%5Cmathbf%7Bx%7D%5Ea%24%3A%20the%20DA%20estimate%20we%20want%20to%20evaluate.%0A%0A%20%20%20%20**Error%20metric**%3A%20Root-Mean-Square%20Error%20(RMSE)%20between%20analysis%2Fbackground%20and%20truth.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(integrate_l63%2C%20np%2C%20rng)%3A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20Experiment%20parameters%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20DT_MODEL%20%20%3D%200.01%20%20%20%20%20%20%20%20%20%20%23%20model%20time%20step%20(non-dimensional)%0A%20%20%20%20T_SPINUP%20%20%3D%205.0%20%20%20%20%20%20%20%20%20%20%20%23%20spin-up%20to%20attractor%0A%20%20%20%20T_ASSIM%20%20%20%3D%205.0%20%20%20%20%20%20%20%20%20%20%20%23%20assimilation%20window%20length%0A%20%20%20%20DT_OBS%20%20%20%20%3D%200.2%20%20%20%20%20%20%20%20%20%20%20%23%20observation%20interval%0A%20%20%20%20OBS_SIGMA%20%3D%202.0%20%20%20%20%20%20%20%20%20%20%20%23%20observation%20error%20std%20(each%20component)%0A%20%20%20%20BG_SIGMA%20%20%3D%203.0%20%20%20%20%20%20%20%20%20%20%20%23%20background%20error%20std%20(each%20component)%0A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20Spin-up%3A%20land%20on%20the%20attractor%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20t_spinup%20%3D%20np.linspace(0%2C%20T_SPINUP%2C%20int(round(T_SPINUP%20%2F%20DT_MODEL))%20%2B%201)%0A%20%20%20%20x0_truth%20%3D%20np.array(%5B1.0%2C%201.0%2C%2020.0%5D)%0A%20%20%20%20spinup%20%20%20%3D%20integrate_l63(x0_truth%2C%20(0%2C%20T_SPINUP)%2C%20t_spinup)%0A%20%20%20%20x_start%20%20%3D%20spinup%5B-1%5D%20%20%23%20point%20on%20attractor%0A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20Nature%20run%20(truth)%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20n_steps%20%20%3D%20int(round(T_ASSIM%20%2F%20DT_MODEL))%20%2B%201%0A%20%20%20%20t_assim%20%20%3D%20np.linspace(0%2C%20T_ASSIM%2C%20n_steps)%0A%20%20%20%20truth%20%20%20%20%3D%20integrate_l63(x_start%2C%20(0%2C%20T_ASSIM)%2C%20t_assim)%20%20%23%20(n_steps%2C%203)%0A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20Synthetic%20observations%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20%23%20Observe%20all%203%20components%20(H%20%3D%20I_3).%20%20Can%20be%20made%20partial%20easily.%0A%20%20%20%20obs_times_idx%20%3D%20np.arange(0%2C%20n_steps%2C%20int(DT_OBS%20%2F%20DT_MODEL))%0A%20%20%20%20obs_times%20%20%20%20%20%3D%20t_assim%5Bobs_times_idx%5D%0A%20%20%20%20n_obs_times%20%20%20%3D%20len(obs_times_idx)%0A%20%20%20%20R_diag%20%20%20%20%20%20%20%20%3D%20OBS_SIGMA**2%20*%20np.ones(3)%0A%20%20%20%20R%20%20%20%20%20%20%20%20%20%20%20%20%20%3D%20np.diag(R_diag)%0A%20%20%20%20R_inv%20%20%20%20%20%20%20%20%20%3D%20np.diag(1.0%20%2F%20R_diag)%0A%20%20%20%20observations%20%20%3D%20truth%5Bobs_times_idx%5D%20%2B%20rng.normal(0%2C%20OBS_SIGMA%2C%20(n_obs_times%2C%203))%0A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20Background%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20%23%20Perturb%20the%20true%20IC%20to%20create%20an%20%22imperfect%22%20background%0A%20%20%20%20x_bg%20%20%20%3D%20x_start%20%2B%20rng.normal(0%2C%20BG_SIGMA%2C%203)%0A%20%20%20%20B_diag%20%3D%20BG_SIGMA**2%20*%20np.ones(3)%0A%20%20%20%20B%20%20%20%20%20%20%3D%20np.diag(B_diag)%0A%20%20%20%20B_inv%20%20%3D%20np.diag(1.0%20%2F%20B_diag)%0A%0A%20%20%20%20print(f%22Truth%20start%20%20%20%20%3A%20%7Bx_start%7D%22)%0A%20%20%20%20print(f%22Background%20start%3A%20%7Bx_bg%7D%22)%0A%20%20%20%20print(f%22Initial%20bg%20error%3A%20%7Bnp.linalg.norm(x_bg%20-%20x_start)%3A.3f%7D%22)%0A%20%20%20%20return%20(%0A%20%20%20%20%20%20%20%20B_diag%2C%0A%20%20%20%20%20%20%20%20B_inv%2C%0A%20%20%20%20%20%20%20%20DT_MODEL%2C%0A%20%20%20%20%20%20%20%20R%2C%0A%20%20%20%20%20%20%20%20R_inv%2C%0A%20%20%20%20%20%20%20%20obs_times%2C%0A%20%20%20%20%20%20%20%20obs_times_idx%2C%0A%20%20%20%20%20%20%20%20observations%2C%0A%20%20%20%20%20%20%20%20t_assim%2C%0A%20%20%20%20%20%20%20%20truth%2C%0A%20%20%20%20%20%20%20%20x_bg%2C%0A%20%20%20%20%20%20%20%20x_start%2C%0A%20%20%20%20)%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---%0A%20%20%20%20%23%23%203%20%203DVAR%20%E2%80%94%20Three-Dimensional%20Variational%20Assimilation%0A%0A%20%20%20%20%23%23%23%20Theory%0A%0A%20%20%20%203DVAR%20finds%20the%20analysis%20%24%5Cmathbf%7Bx%7D%5Ea%24%20by%20minimising%20the%20**cost%20function**%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cmathcal%7BJ%7D(%5Cmathbf%7Bx%7D)%20%3D%0A%20%20%20%20%5Cunderbrace%7B%5Cfrac%7B1%7D%7B2%7D(%5Cmathbf%7Bx%7D%20-%20%5Cmathbf%7Bx%7D%5Eb)%5ET%20%5Cmathbf%7BB%7D%5E%7B-1%7D(%5Cmathbf%7Bx%7D%20-%20%5Cmathbf%7Bx%7D%5Eb)%7D_%7B%5Cmathcal%7BJ%7D_b%20%5Ctext%7B%20%E2%80%94%20background%20term%7D%7D%0A%20%20%20%20%2B%0A%20%20%20%20%5Cunderbrace%7B%5Cfrac%7B1%7D%7B2%7D(%5Cmathbf%7By%7D%20-%20H%5Cmathbf%7Bx%7D)%5ET%20%5Cmathbf%7BR%7D%5E%7B-1%7D(%5Cmathbf%7By%7D%20-%20H%5Cmathbf%7Bx%7D)%7D_%7B%5Cmathcal%7BJ%7D_o%20%5Ctext%7B%20%E2%80%94%20observation%20term%7D%7D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20-%20%24%5Cmathbf%7Bx%7D%5Eb%20%5Cin%20%5Cmathbb%7BR%7D%5En%24%20%E2%80%94%20background%20(prior)%20state%0A%20%20%20%20-%20%24%5Cmathbf%7BB%7D%20%5Cin%20%5Cmathbb%7BR%7D%5E%7Bn%20%5Ctimes%20n%7D%24%20%E2%80%94%20background%20error%20covariance%0A%20%20%20%20-%20%24%5Cmathbf%7By%7D%20%5Cin%20%5Cmathbb%7BR%7D%5Em%24%20%E2%80%94%20observation%20vector%0A%20%20%20%20-%20%24H%3A%20%5Cmathbb%7BR%7D%5En%20%5Cto%20%5Cmathbb%7BR%7D%5Em%24%20%E2%80%94%20observation%20operator%0A%20%20%20%20-%20%24%5Cmathbf%7BR%7D%20%5Cin%20%5Cmathbb%7BR%7D%5E%7Bm%20%5Ctimes%20m%7D%24%20%E2%80%94%20observation%20error%20covariance%0A%0A%20%20%20%20**Analytic%20solution**%20(linear%20%24H%24)%3A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cmathbf%7Bx%7D%5Ea%20%3D%20%5Cmathbf%7Bx%7D%5Eb%20%2B%20%5Cmathbf%7BK%7D(%5Cmathbf%7By%7D%20-%20H%5Cmathbf%7Bx%7D%5Eb)%2C%20%5Cquad%0A%20%20%20%20%5Cmathbf%7BK%7D%20%3D%20%5Cmathbf%7BB%7D%20H%5ET%20(H%20%5Cmathbf%7BB%7D%20H%5ET%20%2B%20%5Cmathbf%7BR%7D)%5E%7B-1%7D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20The%20gradient%20needed%20by%20the%20minimiser%20is%3A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cnabla_%7B%5Cmathbf%7Bx%7D%7D%20%5Cmathcal%7BJ%7D%20%3D%20%5Cmathbf%7BB%7D%5E%7B-1%7D(%5Cmathbf%7Bx%7D%20-%20%5Cmathbf%7Bx%7D%5Eb)%20-%20H%5ET%20%5Cmathbf%7BR%7D%5E%7B-1%7D(%5Cmathbf%7By%7D%20-%20H%5Cmathbf%7Bx%7D)%0A%20%20%20%20%24%24%0A%0A%20%20%20%20%23%23%23%20Key%20limitations%0A%0A%20%20%20%20*%20**Atemporal**%3A%20uses%20observations%20at%20a%20single%20time%20%E2%80%94%20no%20memory%20of%20the%20trajectory.%0A%20%20%20%20*%20**Static%20%24%5Cmathbf%7BB%7D%24**%3A%20does%20not%20adapt%20to%20the%20flow%20of%20the%20day.%0A%20%20%20%20*%20Simple%20to%20implement%20and%20fast%20%E2%80%94%20still%20used%20operationally%20in%20many%20regional%20NWP%20systems.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20B_inv%2C%0A%20%20%20%20DT_MODEL%2C%0A%20%20%20%20R_inv%2C%0A%20%20%20%20integrate_l63%2C%0A%20%20%20%20minimize%2C%0A%20%20%20%20np%2C%0A%20%20%20%20obs_times_idx%2C%0A%20%20%20%20observations%2C%0A%20%20%20%20t_assim%2C%0A%20%20%20%20truth%2C%0A%20%20%20%20x_bg%2C%0A)%3A%0A%20%20%20%20def%20var3d_cost_and_grad(x%2C%20xb%2C%20B_inv%2C%20y_obs%2C%20H%2C%20R_inv)%3A%0A%20%20%20%20%20%20%20%20%22%22%223DVAR%20cost%20function%20J%20and%20its%20gradient%20%E2%88%87J.%22%22%22%0A%20%20%20%20%20%20%20%20innov%20%3D%20y_obs%20-%20H%20%40%20x%0A%20%20%20%20%20%20%20%20bg_departure%20%3D%20x%20-%20xb%0A%20%20%20%20%20%20%20%20J%20%20%3D%200.5%20*%20bg_departure%20%40%20B_inv%20%40%20bg_departure%20%2B%200.5%20*%20innov%20%40%20R_inv%20%40%20innov%0A%20%20%20%20%20%20%20%20dJ%20%3D%20B_inv%20%40%20bg_departure%20-%20H.T%20%40%20R_inv%20%40%20innov%0A%20%20%20%20%20%20%20%20return%20J%2C%20dJ%0A%0A%20%20%20%20def%20run_3dvar(xb%2C%20B_inv%2C%20obs_seq%2C%20obs_idx%2C%20R_inv%2C%20t_grid%2C%20dt)%3A%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20Cycling%203DVAR%3A%20at%20each%20observation%20time%2C%20assimilate%20the%20observation%2C%0A%20%20%20%20%20%20%20%20then%20propagate%20the%20analysis%20forward%20to%20the%20next%20observation%20time.%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20H%20%3D%20np.eye(3)%20%20%20%23%20observe%20all%20components%0A%20%20%20%20%20%20%20%20x_curr%20%3D%20xb.copy()%0A%20%20%20%20%20%20%20%20analyses%20%20%3D%20%5B%5D%0A%20%20%20%20%20%20%20%20fcst_traj%20%3D%20%5B%5D%0A%0A%20%20%20%20%20%20%20%20for%20k%2C%20idx%20in%20enumerate(obs_idx)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20y_k%20%3D%20obs_seq%5Bk%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20t_start%20%3D%20t_grid%5Bidx%5D%20if%20k%20%3D%3D%200%20else%20t_grid%5Bobs_idx%5Bk%20-%201%5D%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20t_end%20%20%20%3D%20t_grid%5Bidx%5D%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20Forecast%20to%20obs%20time%20(except%20at%20k%3D0%20we%20are%20already%20there)%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20k%20%3E%200%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_dt_seg%20%3D%20t_end%20-%20t_start%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20t_seg%20%20%20%3D%20np.linspace(0%2C%20_dt_seg%2C%20max(2%2C%20int(round(_dt_seg%20%2F%20dt))%20%2B%201))%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20segment%20%3D%20integrate_l63(x_curr%2C%20(0%2C%20_dt_seg)%2C%20t_seg)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20fcst_traj.append((t_start%2C%20t_end%2C%20t_seg%2C%20segment))%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20x_curr%20%20%3D%20segment%5B-1%5D%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20Minimise%20J%0A%20%20%20%20%20%20%20%20%20%20%20%20res%20%3D%20minimize(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20var3d_cost_and_grad%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20x_curr%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20jac%3DTrue%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20args%3D(x_curr%2C%20B_inv%2C%20y_k%2C%20H%2C%20R_inv)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20method%3D%22L-BFGS-B%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20options%3D%7B%22maxiter%22%3A%20200%2C%20%22ftol%22%3A%201e-12%2C%20%22gtol%22%3A%201e-8%7D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20%20%20%20%20x_curr%20%3D%20res.x%0A%20%20%20%20%20%20%20%20%20%20%20%20analyses.append((t_grid%5Bidx%5D%2C%20x_curr.copy()))%0A%0A%20%20%20%20%20%20%20%20return%20analyses%2C%20fcst_traj%0A%0A%20%20%20%20analyses_3dvar%2C%20_%20%3D%20run_3dvar(%0A%20%20%20%20%20%20%20%20x_bg%2C%20B_inv%2C%20observations%2C%20obs_times_idx%2C%20R_inv%2C%20t_assim%2C%20DT_MODEL%0A%20%20%20%20)%0A%0A%20%20%20%20%23%20Reconstruct%20full%20trajectory%20from%20analyses%0A%20%20%20%20%23%20Propagate%20each%20analysis%20to%20the%20next%20obs%20time%0A%20%20%20%20traj_3dvar%20%3D%20np.full_like(truth%2C%20np.nan)%0A%20%20%20%20traj_3dvar%5B0%5D%20%3D%20x_bg%0A%20%20%20%20for%20_k%20in%20range(len(analyses_3dvar)%20-%201)%3A%0A%20%20%20%20%20%20%20%20_t0%2C%20_xa%20%3D%20analyses_3dvar%5B_k%5D%0A%20%20%20%20%20%20%20%20_t1%2C%20_%20%20%20%3D%20analyses_3dvar%5B_k%20%2B%201%5D%0A%20%20%20%20%20%20%20%20_ddt%20%20%20%20%20%3D%20_t1%20-%20_t0%0A%20%20%20%20%20%20%20%20_seg_t%20%20%20%3D%20np.linspace(0%2C%20_ddt%2C%20max(2%2C%20int(round(_ddt%20%2F%20DT_MODEL))%20%2B%201))%0A%20%20%20%20%20%20%20%20_seg%20%20%20%20%20%3D%20integrate_l63(_xa%2C%20(0%2C%20_ddt)%2C%20_seg_t)%0A%20%20%20%20%20%20%20%20_i0%20%20%20%20%20%20%3D%20obs_times_idx%5B_k%5D%0A%20%20%20%20%20%20%20%20_i1%20%20%20%20%20%20%3D%20obs_times_idx%5B_k%20%2B%201%5D%0A%20%20%20%20%20%20%20%20_n%20%20%20%20%20%20%20%3D%20min(len(_seg)%2C%20_i1%20-%20_i0%20%2B%201)%0A%20%20%20%20%20%20%20%20traj_3dvar%5B_i0%3A_i0%20%2B%20_n%5D%20%3D%20_seg%5B%3A_n%5D%0A%0A%20%20%20%20%23%20RMSE%20at%20analysis%20times%0A%20%20%20%20rmse_3dvar_bg%20%3D%20np.sqrt(np.mean(%0A%20%20%20%20%20%20%20%20%5B(truth%5Bobs_times_idx%5Bk%5D%5D%20-%20x_bg%20if%20k%20%3D%3D%200%20else%0A%20%20%20%20%20%20%20%20%20%20truth%5Bobs_times_idx%5Bk%5D%5D%20-%20analyses_3dvar%5Bk%20-%201%5D%5B1%5D)**2%0A%20%20%20%20%20%20%20%20%20for%20k%20in%20range(len(analyses_3dvar))%5D%0A%20%20%20%20))%0A%20%20%20%20rmse_3dvar_an%20%3D%20np.sqrt(np.mean(%0A%20%20%20%20%20%20%20%20%5B(truth%5Bobs_times_idx%5Bk%5D%5D%20-%20xa)**2%0A%20%20%20%20%20%20%20%20%20for%20k%2C%20(_%2C%20xa)%20in%20enumerate(analyses_3dvar)%5D%0A%20%20%20%20))%0A%20%20%20%20print(f%223DVAR%20%7C%20Background%20RMSE%3A%20%7Brmse_3dvar_bg%3A.3f%7D%20%20Analysis%20RMSE%3A%20%7Brmse_3dvar_an%3A.3f%7D%22)%0A%20%20%20%20return%20analyses_3dvar%2C%20rmse_3dvar_an%2C%20run_3dvar%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%23%203DVAR%20%E2%80%94%20Results%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(analyses_3dvar%2C%20obs_times%2C%20observations%2C%20plt%2C%20t_assim%2C%20truth%2C%20x_bg)%3A%0A%20%20%20%20_fig%2C%20_axes%20%3D%20plt.subplots(3%2C%201%2C%20figsize%3D(12%2C%207)%2C%20sharex%3DTrue)%0A%20%20%20%20_labels%20%3D%20%5B%22x%22%2C%20%22y%22%2C%20%22z%22%5D%0A%20%20%20%20_colors%20%3D%20%5B%22tab%3Ablue%22%2C%20%22tab%3Aorange%22%2C%20%22tab%3Agreen%22%5D%0A%0A%20%20%20%20for%20_i%20in%20range(3)%3A%0A%20%20%20%20%20%20%20%20_ax%20%3D%20_axes%5B_i%5D%0A%20%20%20%20%20%20%20%20_ax.plot(t_assim%2C%20truth%5B%3A%2C%20_i%5D%2C%20%22k-%22%2C%20lw%3D1.5%2C%20label%3D%22Truth%22)%0A%20%20%20%20%20%20%20%20_ax.plot(%0A%20%20%20%20%20%20%20%20%20%20%20%20obs_times%2C%20observations%5B%3A%2C%20_i%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%22.%22%2C%20color%3D%22gray%22%2C%20ms%3D4%2C%20alpha%3D0.6%2C%20label%3D%22Obs%22%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20%23%20background%20trajectory%20(free%20forecast%20from%20perturbed%20IC)%0A%20%20%20%20%20%20%20%20_ax.axhline(x_bg%5B_i%5D%2C%20ls%3D%22--%22%2C%20color%3D%22salmon%22%2C%20lw%3D1%2C%20label%3D%22Background%20IC%22)%0A%20%20%20%20%20%20%20%20%23%20analysis%20points%0A%20%20%20%20%20%20%20%20_an_times%20%3D%20%5Bt%20for%20t%2C%20_%20in%20analyses_3dvar%5D%0A%20%20%20%20%20%20%20%20_an_vals%20%20%3D%20%5Bxa%5B_i%5D%20for%20_%2C%20xa%20in%20analyses_3dvar%5D%0A%20%20%20%20%20%20%20%20_ax.plot(_an_times%2C%20_an_vals%2C%20%22o%22%2C%20color%3D_colors%5B_i%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20ms%3D5%2C%20zorder%3D5%2C%20label%3D%223DVAR%20analysis%22)%0A%20%20%20%20%20%20%20%20_ax.set_ylabel(_labels%5B_i%5D%2C%20fontsize%3D12)%0A%20%20%20%20%20%20%20%20_ax.legend(loc%3D%22upper%20right%22%2C%20fontsize%3D7%2C%20ncol%3D4)%0A%20%20%20%20_axes%5B-1%5D.set_xlabel(%22time%22)%0A%20%20%20%20_fig.suptitle(%223DVAR%3A%20Cycling%20Analysis%20vs%20Truth%22%2C%20fontsize%3D13)%0A%20%20%20%20_fig.tight_layout()%0A%20%20%20%20_fig%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---%0A%20%20%20%20%23%23%204%20%204DVAR%20%E2%80%94%20Four-Dimensional%20Variational%20Assimilation%0A%0A%20%20%20%20%23%23%23%20Theory%0A%0A%20%20%20%204DVAR%20extends%203DVAR%20to%20**a%20time%20window**%20%24%5Bt_0%2C%20t_N%5D%24.%0A%20%20%20%20We%20seek%20the%20initial%20condition%20%24%5Cmathbf%7Bx%7D_0%24%20that%2C%20when%20propagated%20by%20the%20model%0A%20%20%20%20%24%5Cmathcal%7BM%7D%24%2C%20best%20fits%20*all*%20observations%20in%20the%20window%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cmathcal%7BJ%7D(%5Cmathbf%7Bx%7D_0)%20%3D%0A%20%20%20%20%5Cfrac%7B1%7D%7B2%7D(%5Cmathbf%7Bx%7D_0%20-%20%5Cmathbf%7Bx%7D_0%5Eb)%5ET%20%5Cmathbf%7BB%7D%5E%7B-1%7D(%5Cmathbf%7Bx%7D_0%20-%20%5Cmathbf%7Bx%7D_0%5Eb)%0A%20%20%20%20%2B%0A%20%20%20%20%5Csum_%7Bk%3D0%7D%5E%7BN%7D%0A%20%20%20%20%5Cfrac%7B1%7D%7B2%7D(%5Cmathbf%7By%7D_k%20-%20H%5Cmathbf%7Bx%7D_k)%5ET%20%5Cmathbf%7BR%7D%5E%7B-1%7D(%5Cmathbf%7By%7D_k%20-%20H%5Cmathbf%7Bx%7D_k)%0A%20%20%20%20%24%24%0A%0A%20%20%20%20where%20%24%5Cmathbf%7Bx%7D_k%20%3D%20%5Cmathcal%7BM%7D_%7B0%20%5Cto%20t_k%7D(%5Cmathbf%7Bx%7D_0)%24.%0A%0A%20%20%20%20**Gradient**%20requires%20the%20**adjoint%20model**%20%24%5Cmathcal%7BM%7D%5E%7B%5Cast%7D%24%20(backward%20integration)%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cnabla_%7B%5Cmathbf%7Bx%7D_0%7D%5Cmathcal%7BJ%7D%20%3D%20%5Cmathbf%7BB%7D%5E%7B-1%7D(%5Cmathbf%7Bx%7D_0%20-%20%5Cmathbf%7Bx%7D_0%5Eb)%20%2B%0A%20%20%20%20%5Cmathcal%7BM%7D%5E%7B%5Cast%7D_%7Bt_N%20%5Cto%20t_0%7D%5Cleft%5B%0A%20%20%20%20%5Csum_k%20H%5ET%20%5Cmathbf%7BR%7D%5E%7B-1%7D(%5Cmathbf%7By%7D_k%20-%20H%5Cmathbf%7Bx%7D_k)%5Cdelta(t%20-%20t_k)%0A%20%20%20%20%5Cright%5D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20Here%20we%20approximate%20the%20gradient%20via%20**finite%20differences**%20(also%20called%0A%20%20%20%20implicit-differentiation%20%2F%20ensemble%20gradient%20in%20some%20contexts).%20For%20production%0A%20%20%20%20systems%2C%20an%20actual%20adjoint%20code%20is%20used.%0A%0A%20%20%20%20%23%23%23%20Advantage%20over%203DVAR%0A%20%20%20%20*%20Uses%20observations%20distributed%20over%20a%20**time%20window**%20%E2%86%92%20the%20model%20trajectory%20is%0A%20%20%20%20%20%20simultaneously%20fit%20to%20multiple%20obs%20times.%0A%20%20%20%20*%20Implicitly%20captures%20**flow-dependent**%20error%20growth%20during%20the%20window%20via%20the%0A%20%20%20%20%20%20model's%20nonlinear%20dynamics.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20B_inv%2C%0A%20%20%20%20DT_MODEL%2C%0A%20%20%20%20R_inv%2C%0A%20%20%20%20integrate_l63%2C%0A%20%20%20%20minimize%2C%0A%20%20%20%20np%2C%0A%20%20%20%20obs_times_idx%2C%0A%20%20%20%20observations%2C%0A%20%20%20%20t_assim%2C%0A%20%20%20%20truth%2C%0A%20%20%20%20x_bg%2C%0A)%3A%0A%20%20%20%20def%20_linspace_seg(t_end%2C%20dt)%3A%0A%20%20%20%20%20%20%20%20n%20%3D%20max(2%2C%20int(round(t_end%20%2F%20dt))%20%2B%201)%0A%20%20%20%20%20%20%20%20return%20np.linspace(0.0%2C%20t_end%2C%20n)%0A%0A%20%20%20%20def%20var4d_cost_and_grad(x0%2C%20xb%2C%20B_inv%2C%20obs_window%2C%20obs_idx_window%2C%20t_grid%2C%20dt%2C%20R_inv%2C%20eps%3D1e-5)%3A%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%204DVAR%20cost%20J%20and%20gradient%20%E2%88%87J%20w.r.t.%20x0.%0A%20%20%20%20%20%20%20%20Gradient%20computed%20via%20finite%20differences%20(approximates%20adjoint).%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20H%20%3D%20np.eye(3)%0A%20%20%20%20%20%20%20%20n%20%3D%20len(x0)%0A%0A%20%20%20%20%20%20%20%20def%20_J(x0_)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20t_end%20%3D%20t_grid%5Bobs_idx_window%5B-1%5D%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20Jb%20%3D%200.5%20*%20(x0_%20-%20xb)%20%40%20B_inv%20%40%20(x0_%20-%20xb)%0A%20%20%20%20%20%20%20%20%20%20%20%20Jo%20%3D%200.0%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20t_end%20%3C%20dt%20*%200.5%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%23%20Degenerate%20single-obs%20window%20at%20t%3D0%3A%20no%20model%20propagation%20needed%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20for%20k%2C%20idx%20in%20enumerate(obs_idx_window)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20innov%20%3D%20obs_window%5Bk%5D%20-%20H%20%40%20x0_%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20Jo%20%20%20%2B%3D%200.5%20*%20innov%20%40%20R_inv%20%40%20innov%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%20t_seg%20%3D%20_linspace_seg(t_end%2C%20dt)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20traj%20%20%3D%20integrate_l63(x0_%2C%20(0.0%2C%20t_end)%2C%20t_seg)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20for%20k%2C%20idx%20in%20enumerate(obs_idx_window)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%23%20map%20obs%20time%20to%20nearest%20traj%20index%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20ti%20%3D%20t_grid%5Bidx%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20li%20%3D%20min(int(round(ti%20%2F%20dt))%2C%20len(traj)%20-%201)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20innov%20%3D%20obs_window%5Bk%5D%20-%20H%20%40%20traj%5Bli%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20Jo%20%20%20%2B%3D%200.5%20*%20innov%20%40%20R_inv%20%40%20innov%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20Jb%20%2B%20Jo%0A%0A%20%20%20%20%20%20%20%20J0%20%3D%20_J(x0)%0A%20%20%20%20%20%20%20%20grad%20%3D%20np.zeros(n)%0A%20%20%20%20%20%20%20%20for%20_j%20in%20range(n)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20x_p%20%20%20%20%20%20%3D%20x0.copy()%3B%20x_p%5B_j%5D%20%2B%3D%20eps%0A%20%20%20%20%20%20%20%20%20%20%20%20grad%5B_j%5D%20%3D%20(_J(x_p)%20-%20J0)%20%2F%20eps%0A%20%20%20%20%20%20%20%20return%20J0%2C%20grad%0A%0A%20%20%20%20def%20run_4dvar(xb%2C%20B_inv%2C%20obs_seq%2C%20obs_idx%2C%20R_inv%2C%20t_grid%2C%20dt%2C%20window_size%3D5)%3A%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20Cycling%204DVAR%3A%20assimilate%20'window_size'%20observations%20at%20a%20time%2C%0A%20%20%20%20%20%20%20%20then%20advance%20the%20analysis%20to%20the%20start%20of%20the%20next%20window.%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20x_curr%20%20%20%3D%20xb.copy()%0A%20%20%20%20%20%20%20%20analyses%20%3D%20%5B%5D%0A%20%20%20%20%20%20%20%20n_obs%20%20%20%20%3D%20len(obs_idx)%0A%0A%20%20%20%20%20%20%20%20for%20k%20in%20range(0%2C%20n_obs%2C%20window_size)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20win_end%20%20%20%20%20%3D%20min(k%20%2B%20window_size%2C%20n_obs)%0A%20%20%20%20%20%20%20%20%20%20%20%20win_obs_idx%20%3D%20obs_idx%5Bk%3Awin_end%5D%20-%20obs_idx%5Bk%5D%20%20%20%23%20relative%20indices%0A%20%20%20%20%20%20%20%20%20%20%20%20win_obs%20%20%20%20%20%3D%20obs_seq%5Bk%3Awin_end%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20t_start%20%20%20%20%20%3D%20t_grid%5Bobs_idx%5Bk%5D%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20t_win_end%20%20%20%3D%20t_grid%5Bobs_idx%5Bwin_end%20-%201%5D%5D%20-%20t_start%0A%20%20%20%20%20%20%20%20%20%20%20%20t_win%20%20%20%20%20%20%20%3D%20_linspace_seg(t_win_end%2C%20dt)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20res%20%3D%20minimize(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20var4d_cost_and_grad%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20x_curr%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20jac%3DTrue%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20args%3D(x_curr%2C%20B_inv%2C%20win_obs%2C%20win_obs_idx%2C%20t_win%2C%20dt%2C%20R_inv)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20method%3D%22L-BFGS-B%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20options%3D%7B%22maxiter%22%3A%20100%2C%20%22ftol%22%3A%201e-10%2C%20%22gtol%22%3A%201e-6%7D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20%20%20%20%20x0_opt%20%3D%20res.x%0A%20%20%20%20%20%20%20%20%20%20%20%20analyses.append((t_start%2C%20x0_opt.copy()))%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20Propagate%20to%20end%20of%20window%20%E2%86%92%20becomes%20background%20for%20next%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20win_end%20%3C%20n_obs%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20t_advance%20%3D%20t_grid%5Bobs_idx%5Bwin_end%5D%5D%20-%20t_start%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20seg%20%20%20%20%20%20%20%3D%20integrate_l63(x0_opt%2C%20(0.0%2C%20t_advance)%2C%20_linspace_seg(t_advance%2C%20dt))%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20x_curr%20%20%20%20%3D%20seg%5B-1%5D%0A%0A%20%20%20%20%20%20%20%20return%20analyses%0A%0A%20%20%20%20analyses_4dvar%20%3D%20run_4dvar(%0A%20%20%20%20%20%20%20%20x_bg%2C%20B_inv%2C%20observations%2C%20obs_times_idx%2C%20R_inv%2C%20t_assim%2C%20DT_MODEL%2C%20window_size%3D5%0A%20%20%20%20)%0A%0A%20%20%20%20%23%20Reconstruct%20trajectory%20from%204DVAR%20analyses%0A%20%20%20%20traj_4dvar%20%3D%20np.full_like(truth%2C%20np.nan)%0A%20%20%20%20traj_4dvar%5B0%5D%20%3D%20x_bg%0A%20%20%20%20for%20_k%20in%20range(len(analyses_4dvar)%20-%201)%3A%0A%20%20%20%20%20%20%20%20_t0%2C%20_xa%20%3D%20analyses_4dvar%5B_k%5D%0A%20%20%20%20%20%20%20%20_t1%2C%20_%20%20%20%3D%20analyses_4dvar%5B_k%20%2B%201%5D%0A%20%20%20%20%20%20%20%20_dt%20%20%20%20%20%20%3D%20_t1%20-%20_t0%0A%20%20%20%20%20%20%20%20_seg%20%20%20%20%20%3D%20integrate_l63(_xa%2C%20(0.0%2C%20_dt)%2C%20_linspace_seg(_dt%2C%20DT_MODEL))%0A%20%20%20%20%20%20%20%20_i0%20%20%20%20%20%20%3D%20int(round(_t0%20%2F%20DT_MODEL))%0A%20%20%20%20%20%20%20%20_n%20%20%20%20%20%20%20%3D%20min(len(_seg)%2C%20int(round(_dt%20%2F%20DT_MODEL))%20%2B%201)%0A%20%20%20%20%20%20%20%20traj_4dvar%5B_i0%3A_i0%20%2B%20_n%5D%20%3D%20_seg%5B%3A_n%5D%0A%0A%20%20%20%20rmse_4dvar%20%3D%20np.sqrt(np.mean(%0A%20%20%20%20%20%20%20%20%5B(truth%5Bobs_times_idx%5Bk%20*%205%5D%5D%20-%20xa)**2%0A%20%20%20%20%20%20%20%20%20for%20k%2C%20(_%2C%20xa)%20in%20enumerate(analyses_4dvar)%0A%20%20%20%20%20%20%20%20%20if%20k%20*%205%20%3C%20len(obs_times_idx)%5D%0A%20%20%20%20))%0A%20%20%20%20print(f%224DVAR%20%7C%20Analysis%20RMSE%20at%20window%20starts%3A%20%7Brmse_4dvar%3A.3f%7D%22)%0A%20%20%20%20return%20analyses_4dvar%2C%20rmse_4dvar%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%23%204DVAR%20%E2%80%94%20Results%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(analyses_4dvar%2C%20obs_times%2C%20observations%2C%20plt%2C%20t_assim%2C%20truth%2C%20x_bg)%3A%0A%20%20%20%20_fig%2C%20_axes%20%3D%20plt.subplots(3%2C%201%2C%20figsize%3D(12%2C%207)%2C%20sharex%3DTrue)%0A%20%20%20%20_labels%20%3D%20%5B%22x%22%2C%20%22y%22%2C%20%22z%22%5D%0A%20%20%20%20_colors%20%3D%20%5B%22tab%3Ablue%22%2C%20%22tab%3Aorange%22%2C%20%22tab%3Agreen%22%5D%0A%20%20%20%20_win%20%3D%205%20%20%20%23%20window_size%20used%20above%0A%0A%20%20%20%20for%20_i%20in%20range(3)%3A%0A%20%20%20%20%20%20%20%20_ax%20%3D%20_axes%5B_i%5D%0A%20%20%20%20%20%20%20%20_ax.plot(t_assim%2C%20truth%5B%3A%2C%20_i%5D%2C%20%22k-%22%2C%20lw%3D1.5%2C%20label%3D%22Truth%22)%0A%20%20%20%20%20%20%20%20_ax.plot(obs_times%2C%20observations%5B%3A%2C%20_i%5D%2C%20%22.%22%2C%20color%3D%22gray%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20ms%3D4%2C%20alpha%3D0.5%2C%20label%3D%22Obs%22)%0A%20%20%20%20%20%20%20%20_ax.axhline(x_bg%5B_i%5D%2C%20ls%3D%22--%22%2C%20color%3D%22salmon%22%2C%20lw%3D1%2C%20label%3D%22Background%20IC%22)%0A%20%20%20%20%20%20%20%20_an_times%20%3D%20%5Bt%20for%20t%2C%20_%20in%20analyses_4dvar%5D%0A%20%20%20%20%20%20%20%20_an_vals%20%20%3D%20%5Bxa%5B_i%5D%20for%20_%2C%20xa%20in%20analyses_4dvar%5D%0A%20%20%20%20%20%20%20%20_ax.plot(_an_times%2C%20_an_vals%2C%20%22s%22%2C%20color%3D_colors%5B_i%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20ms%3D6%2C%20zorder%3D5%2C%20label%3D%224DVAR%20analysis%20(window%20start)%22)%0A%20%20%20%20%20%20%20%20_ax.set_ylabel(_labels%5B_i%5D%2C%20fontsize%3D12)%0A%20%20%20%20%20%20%20%20_ax.legend(loc%3D%22upper%20right%22%2C%20fontsize%3D7%2C%20ncol%3D4)%0A%20%20%20%20_axes%5B-1%5D.set_xlabel(%22time%22)%0A%20%20%20%20_fig.suptitle(%224DVAR%3A%20Cycling%20Analysis%20vs%20Truth%22%2C%20fontsize%3D13)%0A%20%20%20%20_fig.tight_layout()%0A%20%20%20%20_fig%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---%0A%20%20%20%20%23%23%205%20%20Ensemble%20Kalman%20Filter%20(EnKF)%0A%0A%20%20%20%20%23%23%23%20Theory%0A%0A%20%20%20%20The%20EnKF%20(Evensen%201994)%20replaces%20the%20static%20**B**%20with%20a%20**flow-dependent**%0A%20%20%20%20background%20error%20covariance%20estimated%20from%20an%20ensemble%20of%20%24N%24%20model%20trajectories%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cmathbf%7BB%7D%5Ef%20%5Capprox%20%5Cmathbf%7BP%7D%5Ef%20%3D%20%5Cfrac%7B1%7D%7BN-1%7D%5Csum_%7Bi%3D1%7D%5E%7BN%7D%0A%20%20%20%20(%5Cmathbf%7Bx%7D%5Ef_i%20-%20%5Cbar%7B%5Cmathbf%7Bx%7D%7D%5Ef)(%5Cmathbf%7Bx%7D%5Ef_i%20-%20%5Cbar%7B%5Cmathbf%7Bx%7D%7D%5Ef)%5ET%0A%20%20%20%20%24%24%0A%0A%20%20%20%20**Forecast%20step**%3A%20propagate%20each%20ensemble%20member%20through%20the%20model%3A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cmathbf%7Bx%7D%5Ef_i(t%2B1)%20%3D%20%5Cmathcal%7BM%7D(%5Cmathbf%7Bx%7D%5Ea_i(t))%0A%20%20%20%20%24%24%0A%0A%20%20%20%20**Analysis%20step**%20(perturbed-observation%20EnKF)%3A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cmathbf%7Bx%7D%5Ea_i%20%3D%20%5Cmathbf%7Bx%7D%5Ef_i%20%2B%20%5Cmathbf%7BK%7D(%5Cmathbf%7By%7D_i%20-%20H%5Cmathbf%7Bx%7D%5Ef_i)%2C%20%5Cquad%0A%20%20%20%20%5Cmathbf%7By%7D_i%20%3D%20%5Cmathbf%7By%7D%20%2B%20%5Cboldsymbol%7B%5Cepsilon%7D_i%2C%20%5Cquad%0A%20%20%20%20%5Cboldsymbol%7B%5Cepsilon%7D_i%20%5Csim%20%5Cmathcal%7BN%7D(0%2C%20%5Cmathbf%7BR%7D)%0A%20%20%20%20%24%24%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cmathbf%7BK%7D%20%3D%20%5Cmathbf%7BP%7D%5Ef%20H%5ET%20(H%20%5Cmathbf%7BP%7D%5Ef%20H%5ET%20%2B%20%5Cmathbf%7BR%7D)%5E%7B-1%7D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20%23%23%23%20Key%20strengths%0A%20%20%20%20*%20Fully%20nonlinear%20forecast%20step%20%E2%80%94%20no%20adjoint%20needed.%0A%20%20%20%20*%20%24%5Cmathbf%7BP%7D%5Ef%24%20adapts%20to%20the%20local%20geometry%20of%20the%20attractor%20(e.g.%2C%20collapses%20onto%0A%20%20%20%20%20%20the%20unstable%20manifold).%0A%20%20%20%20*%20Natural%20framework%20for%20ensemble%20forecasting.%0A%0A%20%20%20%20%23%23%23%20Practical%20notes%0A%20%20%20%20*%20**Ensemble%20size**%20%24N%24%20is%20critical%3A%20too%20small%20%E2%86%92%20**rank%20deficiency**%20and%0A%20%20%20%20%20%20**spurious%20long-range%20correlations**%20%E2%86%92%20filter%20divergence.%0A%20%20%20%20*%20Remedies%3A%20**inflation**%20(multiply%20%24%5Cmathbf%7BP%7D%5Ef%24%20by%20%241%20%2B%20%5Cdelta%24)%20and%0A%20%20%20%20%20%20**localisation**%20(Gaspari%E2%80%93Cohn%20taper).%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%20N_ens%20%3D%20mo.ui.slider(start%3D10%2C%20stop%3D100%2C%20step%3D10%2C%20value%3D20%2C%20label%3D%22Ensemble%20size%20N%22)%0A%20%20%20%20inflation%20%3D%20mo.ui.slider(start%3D1.0%2C%20stop%3D1.5%2C%20step%3D0.05%2C%20value%3D1.05%2C%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%20label%3D%22Inflation%20factor%20%CE%B1%22)%0A%20%20%20%20mo.vstack(%5B%0A%20%20%20%20%20%20%20%20mo.md(%22**Interactive%20controls**%20%E2%80%94%20change%20these%20and%20the%20EnKF%20results%20below%20will%20update%3A%22)%2C%0A%20%20%20%20%20%20%20%20mo.hstack(%5BN_ens%2C%20inflation%5D)%2C%0A%20%20%20%20%5D)%0A%20%20%20%20return%20N_ens%2C%20inflation%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20B_diag%2C%0A%20%20%20%20DT_MODEL%2C%0A%20%20%20%20N_ens%2C%0A%20%20%20%20R%2C%0A%20%20%20%20inflation%2C%0A%20%20%20%20integrate_l63%2C%0A%20%20%20%20np%2C%0A%20%20%20%20obs_times_idx%2C%0A%20%20%20%20observations%2C%0A%20%20%20%20rng%2C%0A%20%20%20%20t_assim%2C%0A%20%20%20%20truth%2C%0A%20%20%20%20x_start%2C%0A)%3A%0A%20%20%20%20def%20run_enkf(x_start%2C%20N%2C%20B_diag%2C%20obs_seq%2C%20obs_idx%2C%20R%2C%20t_grid%2C%20dt%2C%20alpha%3D1.05)%3A%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20Perturbed-observation%20EnKF%20with%20multiplicative%20inflation%20%CE%B1.%0A%20%20%20%20%20%20%20%20Returns%20ensemble%20mean%20trajectory%20and%20spread%20at%20each%20obs%20time.%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20n_state%20%3D%203%0A%20%20%20%20%20%20%20%20H%20%20%20%20%20%20%20%3D%20np.eye(n_state)%0A%0A%20%20%20%20%20%20%20%20%23%20Initialise%20ensemble%20around%20true%20IC%20with%20background%20uncertainty%0A%20%20%20%20%20%20%20%20ensemble%20%3D%20x_start%20%2B%20rng.normal(0%2C%20np.sqrt(B_diag)%2C%20(N%2C%20n_state))%0A%20%20%20%20%20%20%20%20means%2C%20spreads%2C%20analyses%20%3D%20%5B%5D%2C%20%5B%5D%2C%20%5B%5D%0A%20%20%20%20%20%20%20%20prev_idx%20%3D%200%0A%0A%20%20%20%20%20%20%20%20for%20k%2C%20idx%20in%20enumerate(obs_idx)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20%E2%94%80%E2%94%80%20Forecast%3A%20advance%20each%20member%20from%20prev_idx%20to%20idx%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20%20%20%20%20%20%20%20%20dt_seg%20%20%20%3D%20t_grid%5Bidx%5D%20-%20t_grid%5Bprev_idx%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20dt_seg%20%3E%200%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20t_seg%20%3D%20np.linspace(0%2C%20dt_seg%2C%20max(2%2C%20int(round(dt_seg%20%2F%20dt))%20%2B%201))%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20new_ens%20%3D%20np.array(%5B%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20integrate_l63(ensemble%5Bi%5D%2C%20(0%2C%20dt_seg)%2C%20t_seg)%5B-1%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20for%20i%20in%20range(N)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%5D)%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%20new_ens%20%3D%20ensemble.copy()%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20%E2%94%80%E2%94%80%20Inflate%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20%20%20%20%20%20%20%20%20xf_mean%20%3D%20new_ens.mean(axis%3D0)%0A%20%20%20%20%20%20%20%20%20%20%20%20new_ens%20%3D%20xf_mean%20%2B%20alpha%20*%20(new_ens%20-%20xf_mean)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20%E2%94%80%E2%94%80%20Compute%20Pf%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20%20%20%20%20%20%20%20%20Xf%20%20%3D%20new_ens%20-%20xf_mean%20%20%20%20%20%20%20%20%20%20%20%23%20(N%2C%203)%20anomalies%0A%20%20%20%20%20%20%20%20%20%20%20%20Pf%20%20%3D%20(Xf.T%20%40%20Xf)%20%2F%20(N%20-%201)%20%20%20%20%20%20%23%20(3%2C%203)%20sample%20covariance%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20%E2%94%80%E2%94%80%20Kalman%20gain%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20%20%20%20%20%20%20%20%20S%20%20%20%3D%20H%20%40%20Pf%20%40%20H.T%20%2B%20R%20%20%20%20%20%20%20%20%20%20%20%20%23%20(3%2C%203)%0A%20%20%20%20%20%20%20%20%20%20%20%20K%20%20%20%3D%20Pf%20%40%20H.T%20%40%20np.linalg.inv(S)%20%23%20(3%2C%203)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20%E2%94%80%E2%94%80%20Analysis%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20%20%20%20%20%20%20%20%20y_pert%20%20%20%3D%20obs_seq%5Bk%5D%20%2B%20rng.multivariate_normal(np.zeros(3)%2C%20R%2C%20N)%0A%20%20%20%20%20%20%20%20%20%20%20%20ensemble%20%3D%20new_ens%20%2B%20(y_pert%20-%20(H%20%40%20new_ens.T).T)%20%40%20K.T%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20xa_mean%20%3D%20ensemble.mean(axis%3D0)%0A%20%20%20%20%20%20%20%20%20%20%20%20spread%20%20%3D%20np.sqrt(np.mean(np.var(ensemble%2C%20axis%3D0)))%0A%20%20%20%20%20%20%20%20%20%20%20%20means.append(xa_mean)%0A%20%20%20%20%20%20%20%20%20%20%20%20spreads.append(spread)%0A%20%20%20%20%20%20%20%20%20%20%20%20analyses.append((t_grid%5Bidx%5D%2C%20xa_mean.copy()%2C%20ensemble.copy()))%0A%20%20%20%20%20%20%20%20%20%20%20%20prev_idx%20%3D%20idx%0A%0A%20%20%20%20%20%20%20%20return%20analyses%2C%20means%2C%20spreads%0A%0A%20%20%20%20enkf_results%2C%20enkf_means%2C%20enkf_spreads%20%3D%20run_enkf(%0A%20%20%20%20%20%20%20%20x_start%2C%20N_ens.value%2C%20B_diag%2C%20observations%2C%20obs_times_idx%2C%0A%20%20%20%20%20%20%20%20R%2C%20t_assim%2C%20DT_MODEL%2C%20alpha%3Dinflation.value%0A%20%20%20%20)%0A%0A%20%20%20%20%23%20Evaluate%20RMSE%0A%20%20%20%20rmse_enkf%20%3D%20np.sqrt(np.mean(%0A%20%20%20%20%20%20%20%20%5B(truth%5Bobs_times_idx%5Bk%5D%5D%20-%20m)**2%0A%20%20%20%20%20%20%20%20%20for%20k%2C%20m%20in%20enumerate(enkf_means)%5D%0A%20%20%20%20))%0A%20%20%20%20print(f%22EnKF%20(N%3D%7BN_ens.value%7D%2C%20%CE%B1%3D%7Binflation.value%7D)%20%7C%20Analysis%20RMSE%3A%20%7Brmse_enkf%3A.3f%7D%22)%0A%20%20%20%20return%20enkf_results%2C%20enkf_spreads%2C%20rmse_enkf%2C%20run_enkf%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%23%20EnKF%20%E2%80%94%20Results%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(enkf_results%2C%20np%2C%20obs_times%2C%20observations%2C%20plt%2C%20t_assim%2C%20truth)%3A%0A%20%20%20%20_fig%2C%20_axes%20%3D%20plt.subplots(3%2C%201%2C%20figsize%3D(12%2C%207)%2C%20sharex%3DTrue)%0A%20%20%20%20_labels%20%3D%20%5B%22x%22%2C%20%22y%22%2C%20%22z%22%5D%0A%20%20%20%20_colors%20%3D%20%5B%22tab%3Ablue%22%2C%20%22tab%3Aorange%22%2C%20%22tab%3Agreen%22%5D%0A%0A%20%20%20%20for%20_i%20in%20range(3)%3A%0A%20%20%20%20%20%20%20%20_ax%20%3D%20_axes%5B_i%5D%0A%20%20%20%20%20%20%20%20_ax.plot(t_assim%2C%20truth%5B%3A%2C%20_i%5D%2C%20%22k-%22%2C%20lw%3D1.5%2C%20label%3D%22Truth%22)%0A%20%20%20%20%20%20%20%20_ax.plot(obs_times%2C%20observations%5B%3A%2C%20_i%5D%2C%20%22.%22%2C%20color%3D%22gray%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20ms%3D4%2C%20alpha%3D0.5%2C%20label%3D%22Obs%22)%0A%0A%20%20%20%20%20%20%20%20_an_times%20%3D%20%5Bt%20for%20t%2C%20_%2C%20_%20in%20enkf_results%5D%0A%20%20%20%20%20%20%20%20_an_means%20%3D%20%5Bxa%5B_i%5D%20for%20_%2C%20xa%2C%20_%20in%20enkf_results%5D%0A%20%20%20%20%20%20%20%20_an_spreads%20%3D%20%5B%0A%20%20%20%20%20%20%20%20%20%20%20%20np.std(%5Bens%5Bj%2C%20_i%5D%20for%20j%20in%20range(ens.shape%5B0%5D)%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20for%20_%2C%20_%2C%20ens%20in%20enkf_results%0A%20%20%20%20%20%20%20%20%5D%0A%20%20%20%20%20%20%20%20_ax.plot(_an_times%2C%20_an_means%2C%20%22%5E%22%2C%20color%3D_colors%5B_i%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20ms%3D5%2C%20zorder%3D5%2C%20label%3D%22EnKF%20mean%22)%0A%20%20%20%20%20%20%20%20_ax.fill_between(%0A%20%20%20%20%20%20%20%20%20%20%20%20_an_times%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20np.array(_an_means)%20-%20np.array(_an_spreads)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20np.array(_an_means)%20%2B%20np.array(_an_spreads)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20alpha%3D0.2%2C%20color%3D_colors%5B_i%5D%2C%20label%3D%22%C2%B11%20std%22%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20_ax.set_ylabel(_labels%5B_i%5D%2C%20fontsize%3D12)%0A%20%20%20%20%20%20%20%20_ax.legend(loc%3D%22upper%20right%22%2C%20fontsize%3D7%2C%20ncol%3D5)%0A%20%20%20%20_axes%5B-1%5D.set_xlabel(%22time%22)%0A%20%20%20%20_fig.suptitle(%22EnKF%3A%20Ensemble%20Mean%20%C2%B1%20Spread%20vs%20Truth%22%2C%20fontsize%3D13)%0A%20%20%20%20_fig.tight_layout()%0A%20%20%20%20_fig%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---%0A%20%20%20%20%23%23%206%20%20Side-by-Side%20Comparison%0A%0A%20%20%20%20Let's%20put%20all%20three%20methods%20on%20the%20same%20plot%20and%20compare%20their%20RMSE.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(%0A%20%20%20%20analyses_3dvar%2C%0A%20%20%20%20analyses_4dvar%2C%0A%20%20%20%20enkf_results%2C%0A%20%20%20%20obs_times%2C%0A%20%20%20%20observations%2C%0A%20%20%20%20plt%2C%0A%20%20%20%20rmse_3dvar_an%2C%0A%20%20%20%20rmse_4dvar%2C%0A%20%20%20%20rmse_enkf%2C%0A%20%20%20%20t_assim%2C%0A%20%20%20%20truth%2C%0A)%3A%0A%20%20%20%20_fig%20%3D%20plt.figure(figsize%3D(14%2C%2010))%0A%20%20%20%20_gs%20%20%3D%20_fig.add_gridspec(4%2C%203%2C%20hspace%3D0.45%2C%20wspace%3D0.35)%0A%0A%20%20%20%20_labels%20%3D%20%5B%22x%22%2C%20%22y%22%2C%20%22z%22%5D%0A%20%20%20%20_methods%20%3D%20%5B%223DVAR%22%2C%20%224DVAR%22%2C%20%22EnKF%22%5D%0A%20%20%20%20_rmses%20%20%20%3D%20%5Brmse_3dvar_an%2C%20rmse_4dvar%2C%20rmse_enkf%5D%0A%20%20%20%20_colors%20%20%3D%20%5B%22tab%3Ared%22%2C%20%22tab%3Apurple%22%2C%20%22tab%3Acyan%22%5D%0A%0A%20%20%20%20def%20_get_vals(method%2C%20comp)%3A%0A%20%20%20%20%20%20%20%20if%20method%20%3D%3D%20%223DVAR%22%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%5Bt%20for%20t%2C%20_%20in%20analyses_3dvar%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%5Bxa%5Bcomp%5D%20for%20_%2C%20xa%20in%20analyses_3dvar%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20elif%20method%20%3D%3D%20%224DVAR%22%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%5Bt%20for%20t%2C%20_%20in%20analyses_4dvar%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%5Bxa%5Bcomp%5D%20for%20_%2C%20xa%20in%20analyses_4dvar%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%5Bt%20for%20t%2C%20_%2C%20_%20in%20enkf_results%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%5Bxa%5Bcomp%5D%20for%20_%2C%20xa%2C%20_%20in%20enkf_results%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20)%0A%0A%20%20%20%20for%20_col%2C%20(_method%2C%20_col_color%2C%20_rmse)%20in%20enumerate(zip(_methods%2C%20_colors%2C%20_rmses))%3A%0A%20%20%20%20%20%20%20%20for%20_row%2C%20_lbl%20in%20enumerate(_labels)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_ax%20%3D%20_fig.add_subplot(_gs%5B_row%2C%20_col%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20_ax.plot(t_assim%2C%20truth%5B%3A%2C%20_row%5D%2C%20%22k-%22%2C%20lw%3D1.2%2C%20alpha%3D0.7)%0A%20%20%20%20%20%20%20%20%20%20%20%20_ax.plot(obs_times%2C%20observations%5B%3A%2C%20_row%5D%2C%20%22.%22%2C%20color%3D%22lightgray%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20ms%3D3%2C%20alpha%3D0.5)%0A%20%20%20%20%20%20%20%20%20%20%20%20_t%2C%20_v%20%3D%20_get_vals(_method%2C%20_row)%0A%20%20%20%20%20%20%20%20%20%20%20%20_ax.plot(_t%2C%20_v%2C%20%22o%22%2C%20color%3D_col_color%2C%20ms%3D4%2C%20label%3D_method)%0A%20%20%20%20%20%20%20%20%20%20%20%20_ax.set_ylabel(_lbl%2C%20fontsize%3D10)%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20_row%20%3D%3D%200%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_ax.set_title(f%22%7B_method%7D%5CnRMSE%3D%7B_rmse%3A.3f%7D%22%2C%20fontsize%3D10%2C%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%20color%3D_col_color%2C%20fontweight%3D%22bold%22)%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20_row%20%3D%3D%202%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_ax.set_xlabel(%22time%22%2C%20fontsize%3D9)%0A%0A%20%20%20%20%23%20RMSE%20bar%20chart%0A%20%20%20%20_ax_bar%20%3D%20_fig.add_subplot(_gs%5B3%2C%20%3A%5D)%0A%20%20%20%20_bars%20%20%20%3D%20_ax_bar.bar(_methods%2C%20_rmses%2C%20color%3D_colors%2C%20edgecolor%3D%22k%22%2C%20width%3D0.4)%0A%20%20%20%20for%20_bar%2C%20_r%20in%20zip(_bars%2C%20_rmses)%3A%0A%20%20%20%20%20%20%20%20_ax_bar.text(_bar.get_x()%20%2B%20_bar.get_width()%20%2F%202%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_bar.get_height()%20%2B%200.005%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20f%22%7B_r%3A.3f%7D%22%2C%20ha%3D%22center%22%2C%20va%3D%22bottom%22%2C%20fontsize%3D11)%0A%20%20%20%20_ax_bar.set_ylabel(%22Analysis%20RMSE%22%2C%20fontsize%3D11)%0A%20%20%20%20_ax_bar.set_title(%22RMSE%20Comparison%22%2C%20fontsize%3D11)%0A%20%20%20%20_ax_bar.set_ylim(0%2C%20max(_rmses)%20*%201.3)%0A%0A%20%20%20%20_fig.suptitle(%22Data%20Assimilation%20Methods%20%E2%80%94%20Lorenz%2063%20Testbed%22%2C%20fontsize%3D14%2C%20y%3D1.01)%0A%20%20%20%20_fig%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---%0A%20%20%20%20%23%23%207%20%20Ensemble%20Spread%20vs%20Error%20%E2%80%94%20Reliability%20Diagram%0A%0A%20%20%20%20A%20well-calibrated%20ensemble%20should%20have%20**spread%20%E2%89%88%20RMSE**%20at%20each%20analysis%20time.%0A%20%20%20%20This%20is%20the%20cornerstone%20of%20probabilistic%20DA%20and%20ensemble%20NWP.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(enkf_results%2C%20enkf_spreads%2C%20np%2C%20obs_times%2C%20obs_times_idx%2C%20plt%2C%20truth)%3A%0A%20%20%20%20_errors%20%3D%20%5Bnp.linalg.norm(truth%5Bobs_times_idx%5Bk%5D%5D%20-%20xa)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20for%20k%2C%20(_%2C%20xa%2C%20_)%20in%20enumerate(enkf_results)%5D%0A%0A%20%20%20%20_fig%2C%20_axes%20%3D%20plt.subplots(1%2C%202%2C%20figsize%3D(11%2C%204))%0A%0A%20%20%20%20_axes%5B0%5D.plot(obs_times%2C%20_errors%2C%20%22k-o%22%2C%20ms%3D4%2C%20label%3D%22Analysis%20error%20%E2%80%96e%E2%82%90%E2%80%96%22)%0A%20%20%20%20_axes%5B0%5D.plot(obs_times%2C%20enkf_spreads%2C%20%22b--s%22%2C%20ms%3D4%2C%20label%3D%22Ensemble%20spread%20%CF%83%22)%0A%20%20%20%20_axes%5B0%5D.set_xlabel(%22time%22)%3B%20_axes%5B0%5D.set_ylabel(%22magnitude%22)%0A%20%20%20%20_axes%5B0%5D.set_title(%22Spread%20vs%20Error%20over%20time%22)%0A%20%20%20%20_axes%5B0%5D.legend()%0A%0A%20%20%20%20_axes%5B1%5D.scatter(_errors%2C%20enkf_spreads%2C%20c%3Dobs_times%2C%20cmap%3D%22viridis%22%2C%20s%3D30%2C%20alpha%3D0.8)%0A%20%20%20%20_lim%20%3D%20max(max(_errors)%2C%20max(enkf_spreads))%20*%201.05%0A%20%20%20%20_axes%5B1%5D.plot(%5B0%2C%20_lim%5D%2C%20%5B0%2C%20_lim%5D%2C%20%22r--%22%2C%20lw%3D1.5%2C%20label%3D%22Perfect%20calibration%22)%0A%20%20%20%20_axes%5B1%5D.set_xlabel(%22Analysis%20error%22)%3B%20_axes%5B1%5D.set_ylabel(%22Ensemble%20spread%22)%0A%20%20%20%20_axes%5B1%5D.set_title(%22Reliability%3A%20spread%20vs%20error%22)%0A%20%20%20%20_axes%5B1%5D.legend()%0A%0A%20%20%20%20_fig.tight_layout()%0A%20%20%20%20_fig%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---%0A%20%20%20%20%23%23%208%20%20Observation%20Impact%20%E2%80%94%20Sensitivity%20to%20Observation%20Frequency%0A%0A%20%20%20%20How%20does%20the%20analysis%20quality%20degrade%20as%20we%20observe%20less%20frequently%3F%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20B_diag%2C%0A%20%20%20%20B_inv%2C%0A%20%20%20%20DT_MODEL%2C%0A%20%20%20%20R%2C%0A%20%20%20%20R_inv%2C%0A%20%20%20%20np%2C%0A%20%20%20%20plt%2C%0A%20%20%20%20rng%2C%0A%20%20%20%20run_3dvar%2C%0A%20%20%20%20run_enkf%2C%0A%20%20%20%20t_assim%2C%0A%20%20%20%20truth%2C%0A%20%20%20%20x_bg%2C%0A%20%20%20%20x_start%2C%0A)%3A%0A%20%20%20%20_dt_obs_list%20%20%3D%20%5B0.1%2C%200.2%2C%200.3%2C%200.5%2C%200.8%5D%0A%20%20%20%20_rmse_3dvar_v%20%3D%20%5B%5D%0A%20%20%20%20_rmse_enkf_v%20%20%3D%20%5B%5D%0A%0A%20%20%20%20for%20_dt_obs%20in%20_dt_obs_list%3A%0A%20%20%20%20%20%20%20%20_step%20%20%20%3D%20max(1%2C%20int(_dt_obs%20%2F%20DT_MODEL))%0A%20%20%20%20%20%20%20%20_oi%20%20%20%20%20%3D%20np.arange(0%2C%20len(t_assim)%2C%20_step)%0A%20%20%20%20%20%20%20%20_ot%20%20%20%20%20%3D%20t_assim%5B_oi%5D%0A%20%20%20%20%20%20%20%20_obs%20%20%20%20%3D%20truth%5B_oi%5D%20%2B%20rng.normal(0%2C%202.0%2C%20(len(_oi)%2C%203))%0A%0A%20%20%20%20%20%20%20%20_an3%20%20%3D%20run_3dvar(x_bg%2C%20B_inv%2C%20_obs%2C%20_oi%2C%20R_inv%2C%20t_assim%2C%20DT_MODEL)%0A%20%20%20%20%20%20%20%20_r3%20%20%20%3D%20np.sqrt(np.mean(%5B(truth%5B_oi%5Bk%5D%5D%20-%20xa)**2%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%20for%20k%2C%20(_%2C%20xa)%20in%20enumerate(_an3%5B0%5D)%5D))%0A%20%20%20%20%20%20%20%20_rmse_3dvar_v.append(_r3)%0A%0A%20%20%20%20%20%20%20%20_ank%2C%20_%2C%20_%20%3D%20run_enkf(x_start%2C%2020%2C%20B_diag%2C%20_obs%2C%20_oi%2C%20R%2C%20t_assim%2C%20DT_MODEL)%0A%20%20%20%20%20%20%20%20_re%20%20%3D%20np.sqrt(np.mean(%5B(truth%5B_oi%5Bk%5D%5D%20-%20xa)**2%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%20for%20k%2C%20(_%2C%20xa%2C%20_)%20in%20enumerate(_ank)%5D))%0A%20%20%20%20%20%20%20%20_rmse_enkf_v.append(_re)%0A%0A%20%20%20%20_fig%2C%20_ax%20%3D%20plt.subplots(figsize%3D(7%2C%204))%0A%20%20%20%20_ax.plot(_dt_obs_list%2C%20_rmse_3dvar_v%2C%20%22o-%22%2C%20color%3D%22tab%3Ared%22%2C%20%20label%3D%223DVAR%22)%0A%20%20%20%20_ax.plot(_dt_obs_list%2C%20_rmse_enkf_v%2C%20%20%22%5E-%22%2C%20color%3D%22tab%3Acyan%22%2C%20label%3D%22EnKF%20(N%3D20)%22)%0A%20%20%20%20_ax.set_xlabel(%22Observation%20interval%20%CE%94t_obs%22)%0A%20%20%20%20_ax.set_ylabel(%22Analysis%20RMSE%22)%0A%20%20%20%20_ax.set_title(%22RMSE%20vs%20Observation%20Frequency%22)%0A%20%20%20%20_ax.legend()%0A%20%20%20%20_ax.grid(alpha%3D0.3)%0A%20%20%20%20_fig%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---%0A%20%20%20%20%23%23%209%20%20Key%20Takeaways%0A%0A%20%20%20%20%7C%20Feature%20%7C%203DVAR%20%7C%204DVAR%20%7C%20EnKF%20%7C%0A%20%20%20%20%7C---------%7C-------%7C-------%7C------%7C%0A%20%20%20%20%7C%20**Background%20covariance**%20%7C%20Static%20%24%5Cmathbf%7BB%7D%24%20%7C%20Static%20%24%5Cmathbf%7BB%7D%24%20%7C%20Flow-dependent%20%24%5Cmathbf%7BP%7D%5Ef%24%20%7C%0A%20%20%20%20%7C%20**Temporal%20scope**%20%7C%20Single%20time%20%7C%20Time%20window%20%7C%20Sequential%20%7C%0A%20%20%20%20%7C%20**Model%20adjoint**%20%7C%20Not%20required%20%7C%20Required%20(or%20FD)%20%7C%20Not%20required%20%7C%0A%20%20%20%20%7C%20**Ensemble**%20%7C%20No%20%7C%20No%20%7C%20Yes%20%7C%0A%20%20%20%20%7C%20**Parallelism**%20%7C%20Low%20%7C%20Low%20%7C%20High%20(member%20level)%20%7C%0A%20%20%20%20%7C%20**Typical%20use**%20%7C%20Regional%20NWP%2C%20fast%20cycling%20%7C%20Global%20NWP%20(ECMWF)%20%7C%20Ensemble%20NWP%2C%20ocean%2C%20land%20%7C%0A%0A%20%20%20%20%23%23%23%20Connections%20to%20ML%20you%20already%20know%0A%0A%20%20%20%20*%20**3DVAR**%20is%20Tikhonov-regularised%20least%20squares%20%E2%80%94%20the%20same%20as%20Ridge%20Regression%0A%20%20%20%20%20%20with%20%24%5Cmathbf%7BB%7D%5E%7B-1%7D%24%20as%20the%20regularisation%20matrix.%0A%20%20%20%20*%20**4DVAR**%20is%20**backpropagation%20through%20time**%20(BPTT)%20for%20a%20physics%20model%20%E2%80%94%0A%20%20%20%20%20%20the%20adjoint%20is%20exactly%20the%20reverse-mode%20AD%20graph.%0A%20%20%20%20*%20**EnKF**%20is%20a%20**Monte-Carlo%20Kalman%20Filter**%20%E2%80%94%20analogous%20to%20particle%20methods%0A%20%20%20%20%20%20in%20Bayesian%20deep%20learning.%0A%20%20%20%20*%20Modern%20hybrid%20systems%20(e.g.%20ECMWF's%20IFS-AIFS)%20combine%204DVAR%20with%20neural%0A%20%20%20%20%20%20emulators%20for%20the%20forecast%20model%2C%20making%20the%20adjoint%20tractable%20via%20AD%20frameworks.%0A%0A%20%20%20%20%23%23%23%20Further%20reading%0A%0A%20%20%20%20*%20Kalnay%2C%20E.%20(2003).%20*Atmospheric%20Modeling%2C%20Data%20Assimilation%20and%20Predictability*.%20Cambridge.%0A%20%20%20%20*%20Evensen%2C%20G.%20(2009).%20*Data%20Assimilation%3A%20The%20Ensemble%20Kalman%20Filter*.%20Springer.%0A%20%20%20%20*%20Bocquet%2C%20M.%20et%20al.%20(2023).%20*A%20guide%20to%20ensemble%20Kalman%20methods%20with%20implementation%0A%20%20%20%20%20%20in%20Python*.%20arXiv%3A2305.00087.%0A%20%20%20%20*%20Bauer%2C%20P.%20et%20al.%20(2015).%20The%20quiet%20revolution%20of%20NWP.%20*Nature*%2C%20525%2C%2047%E2%80%9355.%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