+ A = 1 - G.sum(1)
+
+ # warnings.warn("harmonic recurrence", RuntimeWarning)
+ # har = torch.arange(t0, t1, device = G.device).float() + 1
+ # A = har / (har + 1)
+ # G = G / har
+
+ gated_V = torch.einsum("nhrt,nhtd->nrtd", G, V)
+ gated_K = torch.einsum("nhrt,nhtd->nrtd", G, K)
+
+ # We start from cached values, which matters in inference
+
+ init_rec_V = self.rec_V[:, :, t0 - L : t0]
+ init_rec_K = self.rec_K[:, :, t0 - L : t0]
+
+ #################################################################
+ # Associative scan
+
+ # Here there is a trick: Since the stack at position t is
+ # computed by updating that at position t-L, the parallel
+ # scan operates with a period of L. To do so we split the
+ # sequence indexing in two axes, the second of size L, and
+ # run the parallel scan using the first as the sequence index.
+
+ A = A.unflatten(2, (-1, L))
+ gated_V = gated_V.unflatten(2, (-1, L))
+ gated_K = gated_K.unflatten(2, (-1, L))
+
+ next_V = pscan_dim(A, gated_V, init_rec_V, dim=2)
+ next_K = pscan_dim(A, gated_K, init_rec_K, dim=2)
+
+ next_V = next_V.flatten(2, 3)
+ next_K = next_K.flatten(2, 3)
+
+ self.calibrator_rec_V.update(
+ next_V.permute(0, 1, 3, 2).reshape(-1, next_V.size(2))
+ )
+ self.calibrator_rec_K.update(
+ next_K.permute(0, 1, 3, 2).reshape(-1, next_K.size(2))
+ )
+
+ self.rec_V[:, :, t0:t1] = next_V
+ self.rec_K[:, :, t0:t1] = next_K