code: def _twofilter_smoothing_ON(self, t, ti, info, phi, lwinfo, return_ess,
                                 modif_forward, modif_info):
          
          if modif_info is not None:
              lwinfo += modif_info
          Winfo = rs.exp_and_normalise(lwinfo)
          I = rs.multinomial(Winfo)
          if modif_forward is not None:
              lw = self.wgt[t].lw + modif_forward
              W = rs.exp_and_normalise(lw)
          else:
              W = self.wgt[t].W
          J = rs.multinomial(W)
          log_omega = self.model.logpt(t + 1, self.X[t][J], info.hist.X[ti][I])
          if modif_forward is not None:
              log_omega -= modif_forward[J]
          if modif_info is not None:
              log_omega -= modif_info[I]
          Om = rs.exp_and_normalise(log_omega)
          est = np.average(phi(self.X[t][J], info.hist.X[ti][I]), axis=0,
                           weights=Om)
          if return_ess:
              return (est, 1. / np.sum(Om**2))
          else:
              return est
docs: O(N) version of two-filter smoothing.
  
          This method should not be called directly, see twofilter_smoothing.