<jupyter_start><jupyter_text>Quick RunThis notebook is publicly available for any usage at our data imputation project. Please click [**transdim**](https://github.com/xinychen/transdim). We start by importing the necessary dependencies. We will make use of `numpy` and `scipy`.<jupyter_code>import numpy as np
from numpy.linalg import inv as inv<jupyter_output><empty_output><jupyter_text>Part 1: Matrix Computation Concepts 1) Kronecker product- **Definition**:Given two matrices $A\in\mathbb{R}^{m_1\times n_1}$ and $B\in\mathbb{R}^{m_2\times n_2}$, then, the **Kronecker product** between these two matrices is defined as$$A\otimes B=\left[ \begin{array}{cccc} a_{11}B & a_{12}B & \cdots & a_{1m_2}B \\ a_{21}B & a_{22}B & \cdots & a_{2m_2}B \\ \vdots & \vdots & \ddots & \vdots \\ a_{m_11}B & a_{m_12}B & \cdots & a_{m_1m_2}B \\ \end{array} \right]$$where the symbol $\otimes$ denotes Kronecker product, and the size of resulted $A\otimes B$ is $(m_1m_2)\times (n_1n_2)$ (i.e., $m_1\times m_2$ columns and $n_1\times n_2$ rows).- **Example**:If $A=\left[ \begin{array}{cc} 1 & 2 \\ 3 & 4 \\ \end{array} \right]$ and $B=\left[ \begin{array}{ccc} 5 & 6 & 7\\ 8 & 9 & 10 \\ \end{array} \right]$, then, we have$$A\otimes B=\left[ \begin{array}{cc} 1\times \left[ \begin{array}{ccc} 5 & 6 & 7\\ 8 & 9 & 10\\ \end{array} \right] & 2\times \left[ \begin{array}{ccc} 5 & 6 & 7\\ 8 & 9 & 10\\ \end{array} \right] \\ 3\times \left[ \begin{array}{ccc} 5 & 6 & 7\\ 8 & 9 & 10\\ \end{array} \right] & 4\times \left[ \begin{array}{ccc} 5 & 6 & 7\\ 8 & 9 & 10\\ \end{array} \right] \\ \end{array} \right]$$$$=\left[ \begin{array}{cccccc} 5 & 6 & 7 & 10 & 12 & 14 \\ 8 & 9 & 10 & 16 & 18 & 20 \\ 15 & 18 & 21 & 20 & 24 & 28 \\ 24 & 27 & 30 & 32 & 36 & 40 \\ \end{array} \right]\in\mathbb{R}^{4\times 6}.$$ 2) Khatri-Rao product (`kr_prod`)- **Definition**:Given two matrices $A=\left( \boldsymbol{a}_1,\boldsymbol{a}_2,...,\boldsymbol{a}_r \right)\in\mathbb{R}^{m\times r}$ and $B=\left( \boldsymbol{b}_1,\boldsymbol{b}_2,...,\boldsymbol{b}_r \right)\in\mathbb{R}^{n\times r}$ with same number of columns, then, the **Khatri-Rao product** (or **column-wise Kronecker product**) between $A$ and $B$ is given as follows,$$A\odot B=\left( \boldsymbol{a}_1\otimes \boldsymbol{b}_1,\boldsymbol{a}_2\otimes \boldsymbol{b}_2,...,\boldsymbol{a}_r\otimes \boldsymbol{b}_r \right)\in\mathbb{R}^{(mn)\times r},$$where the symbol $\odot$ denotes Khatri-Rao product, and $\otimes$ denotes Kronecker product.- **Example**:If $A=\left[ \begin{array}{cc} 1 & 2 \\ 3 & 4 \\ \end{array} \right]=\left( \boldsymbol{a}_1,\boldsymbol{a}_2 \right) $ and $B=\left[ \begin{array}{cc} 5 & 6 \\ 7 & 8 \\ 9 & 10 \\ \end{array} \right]=\left( \boldsymbol{b}_1,\boldsymbol{b}_2 \right) $, then, we have$$A\odot B=\left( \boldsymbol{a}_1\otimes \boldsymbol{b}_1,\boldsymbol{a}_2\otimes \boldsymbol{b}_2 \right) $$$$=\left[ \begin{array}{cc} \left[ \begin{array}{c} 1 \\ 3 \\ \end{array} \right]\otimes \left[ \begin{array}{c} 5 \\ 7 \\ 9 \\ \end{array} \right] & \left[ \begin{array}{c} 2 \\ 4 \\ \end{array} \right]\otimes \left[ \begin{array}{c} 6 \\ 8 \\ 10 \\ \end{array} \right] \\ \end{array} \right]$$$$=\left[ \begin{array}{cc} 5 & 12 \\ 7 & 16 \\ 9 & 20 \\ 15 & 24 \\ 21 & 32 \\ 27 & 40 \\ \end{array} \right]\in\mathbb{R}^{6\times 2}.$$<jupyter_code>def kr_prod(a, b):
    return np.einsum('ir, jr -> ijr', a, b).reshape(a.shape[0] * b.shape[0], -1)
A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8], [9, 10]])
print(kr_prod(A, B))<jupyter_output>[[ 5 12]
 [ 7 16]
 [ 9 20]
 [15 24]
 [21 32]
 [27 40]]<jupyter_text>3) CP decomposition (`cp_combine`)- **Definition**:The CP decomposition factorizes a tensor into a sum of outer products of vectors. For example, for a third-order tensor $\mathcal{Y}\in\mathbb{R}^{m\times n\times f}$, the CP decomposition can be written as$$\hat{\mathcal{Y}}=\sum_{s=1}^{r}\boldsymbol{u}_{s}\circ\boldsymbol{v}_{s}\circ\boldsymbol{x}_{s},$$or element-wise,$$\hat{y}_{ijt}=\sum_{s=1}^{r}u_{is}v_{js}x_{ts},\forall (i,j,t),$$where vectors $\boldsymbol{u}_{s}\in\mathbb{R}^{m},\boldsymbol{v}_{s}\in\mathbb{R}^{n},\boldsymbol{x}_{s}\in\mathbb{R}^{f}$ are columns of factor matrices $U\in\mathbb{R}^{m\times r},V\in\mathbb{R}^{n\times r},X\in\mathbb{R}^{f\times r}$, respectively. The symbol $\circ$ denotes vector outer product.- **Example**:Given matrices $U=\left[ \begin{array}{cc} 1 & 2 \\ 3 & 4 \\ \end{array} \right]\in\mathbb{R}^{2\times 2}$, $V=\left[ \begin{array}{cc} 1 & 2 \\ 3 & 4 \\ 5 & 6 \\ \end{array} \right]\in\mathbb{R}^{3\times 2}$ and $X=\left[ \begin{array}{cc} 1 & 5 \\ 2 & 6 \\ 3 & 7 \\ 4 & 8 \\ \end{array} \right]\in\mathbb{R}^{4\times 2}$, then if $\hat{\mathcal{Y}}=\sum_{s=1}^{r}\boldsymbol{u}_{s}\circ\boldsymbol{v}_{s}\circ\boldsymbol{x}_{s}$, then, we have$$\hat{Y}_1=\hat{\mathcal{Y}}(:,:,1)=\left[ \begin{array}{ccc} 31 & 42 & 65 \\ 63 & 86 & 135 \\ \end{array} \right],$$$$\hat{Y}_2=\hat{\mathcal{Y}}(:,:,2)=\left[ \begin{array}{ccc} 38 & 52 & 82 \\ 78 & 108 & 174 \\ \end{array} \right],$$$$\hat{Y}_3=\hat{\mathcal{Y}}(:,:,3)=\left[ \begin{array}{ccc} 45 & 62 & 99 \\ 93 & 130 & 213 \\ \end{array} \right],$$$$\hat{Y}_4=\hat{\mathcal{Y}}(:,:,4)=\left[ \begin{array}{ccc} 52 & 72 & 116 \\ 108 & 152 & 252 \\ \end{array} \right].$$<jupyter_code>def cp_combine(U, V, X):
    return np.einsum('is, js, ts -> ijt', U, V, X)
U = np.array([[1, 2], [3, 4]])
V = np.array([[1, 3], [2, 4], [5, 6]])
X = np.array([[1, 5], [2, 6], [3, 7], [4, 8]])
print(cp_combine(U, V, X))
print()
print('tensor size:')
print(cp_combine(U, V, X).shape)<jupyter_output>[[[ 31  38  45  52]
  [ 42  52  62  72]
  [ 65  82  99 116]]

 [[ 63  78  93 108]
  [ 86 108 130 152]
  [135 174 213 252]]]

tensor size:
(2, 3, 4)<jupyter_text>4) Tensor Unfolding (`ten2mat`) and Matrix Folding (`mat2ten`)Using numpy reshape to perform 3rd rank tensor unfold operation. [[**link**](https://stackoverflow.com/questions/49970141/using-numpy-reshape-to-perform-3rd-rank-tensor-unfold-operation)]<jupyter_code>import numpy as np
def ten2mat(tensor, mode):
    return np.reshape(np.moveaxis(tensor, mode, 0), (tensor.shape[mode], -1), order = 'F')
X = np.array([[[1, 2, 3, 4], [3, 4, 5, 6]], [[5, 6, 7, 8], [7, 8, 9, 10]], [[9, 10, 11, 12], [11, 12, 13, 14]]])
print('tensor size:')
print(X.shape)
print('original tensor:')
print(X)
print()
print('(1) mode-1 tensor unfolding:')
print(ten2mat(X, 0))
print()
print('(2) mode-2 tensor unfolding:')
print(ten2mat(X, 1))
print()
print('(3) mode-3 tensor unfolding:')
print(ten2mat(X, 2))
def mat2ten(mat, tensor_size, mode):
    index = list()
    index.append(mode)
    for i in range(tensor_size.shape[0]):
        if i != mode:
            index.append(i)
    return np.moveaxis(np.reshape(mat, list(tensor_size[index]), order = 'F'), 0, mode)<jupyter_output><empty_output><jupyter_text>Part 2: Temporal Regularized Tensor Factorization (TRTF)<jupyter_code>def TRTF(dense_tensor, sparse_tensor, U, V, X, theta, time_lags, multi_steps,
         lambda_u, lambda_v, lambda_ar, eta, lambda_theta, maxiter):
    dim1, dim2, dim3 = dense_tensor.shape
    binary_tensor = np.zeros((dim1, dim2, dim3))
    position = np.where(sparse_tensor > 0)
    binary_tensor[position] = 1
    pos = np.where((dense_tensor > 0) & (sparse_tensor == 0))
    d = len(time_lags)
    rank = U.shape[1]
    
    for iters in range(maxiter):
        var1 = kr_prod(X, V).T
        var2 = kr_prod(var1, var1)
        var3 = np.matmul(var2, ten2mat(binary_tensor, 0).T).reshape([rank, rank, dim1]) + np.dstack([lambda_u * np.eye(rank)] * dim1)
        var4 = np.matmul(var1, ten2mat(sparse_tensor, 0).T)
        for i in range(dim1):
            var_Lambda1 = var3[ :, :, i]
            inv_var_Lambda1 = np.linalg.inv((var_Lambda1 + var_Lambda1.T)/2)
            U[i, :] = np.matmul(inv_var_Lambda1, var4[:, i])

        var1 = kr_prod(X, U).T
        var2 = kr_prod(var1, var1)
        var3 = np.matmul(var2, ten2mat(binary_tensor, 1).T).reshape([rank, rank, dim2]) + np.dstack([lambda_v * np.eye(rank)] * dim2)
        var4 = np.matmul(var1, ten2mat(sparse_tensor, 1).T)
        for j in range(dim2):
            var_Lambda1 = var3[ :, :, j]
            inv_var_Lambda1 = np.linalg.inv((var_Lambda1 + var_Lambda1.T)/2)
            V[j, :] = np.matmul(inv_var_Lambda1, var4[:, j])

        var1 = kr_prod(V, U).T
        var2 = kr_prod(var1, var1)
        var3 = np.matmul(var2, ten2mat(binary_tensor, 2).T).reshape([rank, rank, dim3])
        var4 = np.matmul(var1, ten2mat(sparse_tensor, 2).T)
        for t in range(dim3):
            Mt = np.zeros((rank, rank))
            Nt = np.zeros(rank)
            if t < max(time_lags):
                Pt = np.zeros((rank, rank))
                Qt = np.zeros(rank)
            else:
                Pt = np.eye(rank)
                Qt = np.einsum('ij, ij -> j', theta, X[t - time_lags, :])
            if t < dim3 - np.min(time_lags):
                if t >= np.max(time_lags) and t < dim3 - np.max(time_lags):
                    index = list(range(0, d))
                else:
                    index = list(np.where((t + time_lags >= np.max(time_lags)) & (t + time_lags < dim3)))[0]
                for k in index:
                    theta0 = theta.copy()
                    theta0[k, :] = 0
                    Mt = Mt + np.diag(theta[k, :] ** 2);
                    Nt = Nt + np.multiply(theta[k, :], (X[t + time_lags[k], :] 
                                                        - np.einsum('ij, ij -> j', theta0,
                                                                    X[t + time_lags[k] - time_lags, :])))
                X[t, :] = np.matmul(np.linalg.inv(var3[:, :, t]
                                                  + lambda_ar * Pt + lambda_ar * Mt + lambda_ar * eta * np.eye(rank)),
                                    (var4[:, t] + lambda_ar * Qt + lambda_ar * Nt))
            elif t >= dim3 - np.min(time_lags):
                X[t, :] = np.matmul(np.linalg.inv(var3[:, :, t]
                                                  + lambda_ar * Pt +
                                                  lambda_ar * eta * np.eye(rank)), (var4[:, t] + Qt))

        for k in range(d):
            var1 = X[np.max(time_lags) - time_lags[k] : dim3 - time_lags[k], :]
            var2 = np.linalg.inv(np.diag(np.einsum('ij, ij -> j', var1, var1))
                                 + (lambda_theta / lambda_ar) * np.eye(rank))
            var3 = np.zeros(rank)
            for t in range(np.max(time_lags) - time_lags[k], dim3 - time_lags[k]):
                var3 += np.multiply(X[t, :], (X[t + time_lags[k], :] 
                                              - np.einsum('ij, ij -> j', theta, X[t + time_lags[k] - time_lags, :])
                                              + np.multiply(theta[k, :], X[t, :])))
            theta[k, :] = np.matmul(var2, var3)
            
        tensor_hat = cp_combine(U, V, X)
        mape = np.sum(np.abs(dense_tensor[pos] - 
                             tensor_hat[pos])/dense_tensor[pos])/dense_tensor[pos].shape[0]
        rmse = np.sqrt(np.sum((dense_tensor[pos] -
                               tensor_hat[pos])**2)/dense_tensor[pos].shape[0])
        
        if (iters + 1) % 100 == 0:
            print('Iter: {}'.format(iters + 1))
            print('MAPE: {:.6}'.format(mape))
            print('RMSE: {:.6}'.format(rmse))
            print()
            
    X_new = np.zeros((dim3 + multi_steps, rank))
    X_new[0 : dim3, :] = X.copy()
    for t0 in range(multi_steps):
        X_new[dim3 + t0, :] = np.einsum('ij, ij -> j', theta, X_new[dim3 + t0 - time_lags, :])        
    return cp_combine(U, V, X_new[dim3 : dim3 + multi_steps, :]), U, V, X_new, theta<jupyter_output><empty_output><jupyter_text>Multi-step predictionIn the multi-step prediction task, to enable training data for each rolling step informative, we do not apply an online implementation anymore.Involving rolling prediction tasks, there are two crucial inputs:- **`pred_time_steps`**: the number of steps we should forecast, e.g., if we want to forecast time series within 5 days (144 time slots/steps per day) in advance, then the `pred_time_steps` is $5\times 144=720$;- **`multi_steps`**: the number of steps we should forecast at the current step, e.g., if we want to forecast time series within 2 hours (7 time slots/steps per hour) in advance, then the `multi_stpes` is $2\times 6=12$.<jupyter_code>def multi_prediction(dense_tensor, sparse_tensor, pred_time_steps, rank, time_lags, multi_steps,
                     lambda_u, lambda_v, lambda_ar, eta, lambda_theta, maxiter):
    T = dense_tensor.shape[2]
    start_time = T - pred_time_steps
    dim1 = dense_tensor.shape[0]
    dim2 = dense_tensor.shape[1]
    d = time_lags.shape[0]
    tensor_hat = np.zeros((dim1, dim2, pred_time_steps))
    
    for t in range(int(pred_time_steps/multi_steps)):
        if t == 0:
            ten, U, V, X, theta = TRTF(dense_tensor[:, :, 0 : start_time], sparse_tensor[:, :, 0 : start_time],
                                       0.1 * np.random.rand(dim1, rank), 0.1 * np.random.rand(dim2, rank),
                                       0.1 * np.random.rand(start_time, rank), 0.1 * np.random.rand(d, rank),
                                       time_lags, multi_steps,
                                       lambda_u, lambda_v, lambda_ar, eta, lambda_theta, maxiter[0])
        else:
            ten, U, V, X, theta = TRTF(dense_tensor[:, :, 0 : start_time + t * multi_steps], 
                                       sparse_tensor[:, :, 0 : start_time + t * multi_steps],
                                       U, V, X, theta, time_lags, multi_steps,
                                       lambda_u, lambda_v, lambda_ar, eta, lambda_theta, maxiter[1])
        tensor_hat[:, :, t * multi_steps : (t + 1) * multi_steps] = ten[:, :, ten.shape[2] - multi_steps : ten.shape[2]]

    small_dense_tensor = dense_tensor[:, :, start_time : dense_tensor.shape[2]]
    pos = np.where(small_dense_tensor != 0)
    final_mape = np.sum(np.abs(small_dense_tensor[pos] - 
                               tensor_hat[pos])/small_dense_tensor[pos])/small_dense_tensor[pos].shape[0]
    final_rmse = np.sqrt(np.sum((small_dense_tensor[pos] - 
                                 tensor_hat[pos]) ** 2)/small_dense_tensor[pos].shape[0])
    print('Final MAPE: {:.6}'.format(final_mape))
    print('Final RMSE: {:.6}'.format(final_rmse))
    print()
    return tensor_hat<jupyter_output><empty_output><jupyter_text>Part 3: Experiments on New York Data Set<jupyter_code>import scipy.io

tensor = scipy.io.loadmat('../datasets/NYC-data-set/tensor.mat')
dense_tensor = tensor['tensor']
rm_tensor = scipy.io.loadmat('../datasets/NYC-data-set/rm_tensor.mat')
rm_tensor = rm_tensor['rm_tensor']
nm_tensor = scipy.io.loadmat('../datasets/NYC-data-set/nm_tensor.mat')
nm_tensor = nm_tensor['nm_tensor']

missing_rate = 0.0

# =============================================================================
### Random missing (RM) scenario
### Set the RM scenario by:
binary_tensor = np.round(rm_tensor + 0.5 - missing_rate)
# =============================================================================

sparse_tensor = np.multiply(dense_tensor, binary_tensor)
import time
start = time.time()
pred_time_steps = 24 * 7
multi_steps = 24
rank = 10
time_lags = np.array([1, 2, 3, 24, 24+1, 24+2, 7*24, 7*24+1, 7*24+2])
maxiter = np.array([200, 20])
theta = 0.1 * np.random.rand(time_lags.shape[0], rank)
lambda_u = 500
lambda_v = 500
lambda_ar = 500
eta = 2e-2
lambda_theta = 100
tensor_hat = multi_prediction(dense_tensor, sparse_tensor, pred_time_steps, rank, time_lags, multi_steps,
                              lambda_u, lambda_v, lambda_ar, eta, lambda_theta, maxiter)
end = time.time()
print('Running time: %d seconds'%(end - start))
import scipy.io

tensor = scipy.io.loadmat('../datasets/NYC-data-set/tensor.mat')
dense_tensor = tensor['tensor']
rm_tensor = scipy.io.loadmat('../datasets/NYC-data-set/rm_tensor.mat')
rm_tensor = rm_tensor['rm_tensor']
nm_tensor = scipy.io.loadmat('../datasets/NYC-data-set/nm_tensor.mat')
nm_tensor = nm_tensor['nm_tensor']

missing_rate = 0.1

# =============================================================================
### Random missing (RM) scenario
### Set the RM scenario by:
binary_tensor = np.round(rm_tensor + 0.5 - missing_rate)
# =============================================================================

sparse_tensor = np.multiply(dense_tensor, binary_tensor)
import time
start = time.time()
pred_time_steps = 24 * 7
multi_steps = 24
rank = 10
time_lags = np.array([1, 2, 3, 24, 24+1, 24+2, 7*24, 7*24+1, 7*24+2])
maxiter = np.array([200, 20])
theta = 0.1 * np.random.rand(time_lags.shape[0], rank)
lambda_u = 500
lambda_v = 500
lambda_ar = 500
eta = 2e-2
lambda_theta = 100
tensor_hat = multi_prediction(dense_tensor, sparse_tensor, pred_time_steps, rank, time_lags, multi_steps,
                              lambda_u, lambda_v, lambda_ar, eta, lambda_theta, maxiter)
end = time.time()
print('Running time: %d seconds'%(end - start))
import scipy.io

tensor = scipy.io.loadmat('../datasets/NYC-data-set/tensor.mat')
dense_tensor = tensor['tensor']
rm_tensor = scipy.io.loadmat('../datasets/NYC-data-set/rm_tensor.mat')
rm_tensor = rm_tensor['rm_tensor']
nm_tensor = scipy.io.loadmat('../datasets/NYC-data-set/nm_tensor.mat')
nm_tensor = nm_tensor['nm_tensor']

missing_rate = 0.3

# =============================================================================
### Random missing (RM) scenario
### Set the RM scenario by:
binary_tensor = np.round(rm_tensor + 0.5 - missing_rate)
# =============================================================================

sparse_tensor = np.multiply(dense_tensor, binary_tensor)
import time
start = time.time()
pred_time_steps = 24 * 7
multi_steps = 24
rank = 10
time_lags = np.array([1, 2, 3, 24, 24+1, 24+2, 7*24, 7*24+1, 7*24+2])
maxiter = np.array([200, 20])
theta = 0.1 * np.random.rand(time_lags.shape[0], rank)
lambda_u = 500
lambda_v = 500
lambda_ar = 500
eta = 2e-2
lambda_theta = 100
tensor_hat = multi_prediction(dense_tensor, sparse_tensor, pred_time_steps, rank, time_lags, multi_steps,
                              lambda_u, lambda_v, lambda_ar, eta, lambda_theta, maxiter)
end = time.time()
print('Running time: %d seconds'%(end - start))
import scipy.io

tensor = scipy.io.loadmat('../datasets/NYC-data-set/tensor.mat')
dense_tensor = tensor['tensor']
rm_tensor = scipy.io.loadmat('../datasets/NYC-data-set/rm_tensor.mat')
rm_tensor = rm_tensor['rm_tensor']
nm_tensor = scipy.io.loadmat('../datasets/NYC-data-set/nm_tensor.mat')
nm_tensor = nm_tensor['nm_tensor']

missing_rate = 0.1

# =============================================================================
### Non-random missing (NM) scenario
### Set the NM scenario by:
binary_tensor = np.zeros(dense_tensor.shape)
for i1 in range(dense_tensor.shape[0]):
    for i2 in range(dense_tensor.shape[1]):
        for i3 in range(61):
            binary_tensor[i1, i2, i3 * 24 : (i3 + 1) * 24] = np.round(nm_tensor[i1, i2, i3] 
                                                                      + 0.5 - missing_rate)
# =============================================================================

sparse_tensor = np.multiply(dense_tensor, binary_tensor)
import time
start = time.time()
pred_time_steps = 24 * 7
multi_steps = 24
rank = 10
time_lags = np.array([1, 2, 3, 24, 24+1, 24+2, 7*24, 7*24+1, 7*24+2])
maxiter = np.array([200, 20])
theta = 0.1 * np.random.rand(time_lags.shape[0], rank)
lambda_u = 500
lambda_v = 500
lambda_ar = 500
eta = 2e-2
lambda_theta = 100
tensor_hat = multi_prediction(dense_tensor, sparse_tensor, pred_time_steps, rank, time_lags, multi_steps,
                              lambda_u, lambda_v, lambda_ar, eta, lambda_theta, maxiter)
end = time.time()
print('Running time: %d seconds'%(end - start))
import scipy.io

tensor = scipy.io.loadmat('../datasets/NYC-data-set/tensor.mat')
dense_tensor = tensor['tensor']
rm_tensor = scipy.io.loadmat('../datasets/NYC-data-set/rm_tensor.mat')
rm_tensor = rm_tensor['rm_tensor']
nm_tensor = scipy.io.loadmat('../datasets/NYC-data-set/nm_tensor.mat')
nm_tensor = nm_tensor['nm_tensor']

missing_rate = 0.3

# =============================================================================
### Non-random missing (NM) scenario
### Set the NM scenario by:
binary_tensor = np.zeros(dense_tensor.shape)
for i1 in range(dense_tensor.shape[0]):
    for i2 in range(dense_tensor.shape[1]):
        for i3 in range(61):
            binary_tensor[i1, i2, i3 * 24 : (i3 + 1) * 24] = np.round(nm_tensor[i1, i2, i3] 
                                                                      + 0.5 - missing_rate)
# =============================================================================

sparse_tensor = np.multiply(dense_tensor, binary_tensor)
import time
start = time.time()
pred_time_steps = 24 * 7
multi_steps = 24
rank = 10
time_lags = np.array([1, 2, 3, 24, 24+1, 24+2, 7*24, 7*24+1, 7*24+2])
maxiter = np.array([200, 20])
theta = 0.1 * np.random.rand(time_lags.shape[0], rank)
lambda_u = 500
lambda_v = 500
lambda_ar = 500
eta = 2e-2
lambda_theta = 100
tensor_hat = multi_prediction(dense_tensor, sparse_tensor, pred_time_steps, rank, time_lags, multi_steps,
                              lambda_u, lambda_v, lambda_ar, eta, lambda_theta, maxiter)
end = time.time()
print('Running time: %d seconds'%(end - start))<jupyter_output>Iter: 100
MAPE: 0.861796
RMSE: 7.08837

Iter: 200
MAPE: 0.855741
RMSE: 7.10742

Final MAPE: 0.860435
Final RMSE: 7.22303

Running time: 171 seconds
{'fraction_non_alphanumeric': 0.20530752130418672, 'fraction_numerical': 0.045526120785476104, 'mean_word_length': 3.0098421541318476, 'pattern_counts': {'":': 0, '<': 30, '<?xml version=': 0, '>': 38, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 37, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '147530', 'license': 'MIT', 'n_tokens_mistral': 8527, 'n_tokens_neox': 7686, 'n_words': 1759, 'path': 'experiments/Prediction-Multi-TRTF.ipynb', 'repo_name': 'shawnwang-tech/transdim', 'size': '21592'}