cost functional (w-h(u))^T * R^{-1} * (w-h(u))
def loss(rhs,ObsOp,t,ind_m,u0,w,R,opt,*args):
n = len(u0)
#determine the assimilation window
t = t[:ind_m[-1]+1] #cut the time till the last observation point
nt = len(t)-1
dt = t[1] - t[0]
u = np.zeros([n,nt+1]) #trajectory
u[:,0] = u0
```

```
Ri = np.linalg.inv(R)
floss = 0
km = 0
nt_m = len(ind_m)
if opt == 0: #Euler
#forward model
for k in range(nt):
u[:,k+1] = euler(rhs,u[:,k],dt,*args)
if (km<nt_m) and (k+1==ind_m[km]):
tmp = w[:,km] - ObsOp(u[:,k+1])
tmp = tmp.reshape(-1,1)
floss = floss + np.linalg.multi_dot(( tmp.T, Ri , tmp ))
km = km + 1
elif opt == 1: #RK4