<jupyter_start><jupyter_text>Bayesian Temporal Tensor Factorization**Published**: October 29, 2020**Revised**: October 29, 2020**Download**: This Jupyter notebook is at our GitHub repository. If you want to evaluate the code, please download the notebook from the [**transdim**](https://github.com/xinychen/transdim/blob/master/imputer/BTTF.ipynb) repository.This notebook shows how to implement the Bayesian Temporal Tensor Factorization (BTTF), a fully Bayesian tensor factorization model, on some real-world data sets. To overcome the missing data problem in multidimensional time series, BTTF takes into account both low-rank tensor structure (i.e., CP decomposition) and time series autoregression. For an in-depth discussion of BTTF, please see [1].[1] <NAME>, <NAME> (2019). Bayesian temporal factorization for multidimensional time series prediction. arXiv:1910.06366. [PDF]<jupyter_code>import numpy as np
from numpy.linalg import inv as inv
from numpy.random import normal as normrnd
from scipy.linalg import khatri_rao as kr_prod
from scipy.stats import wishart
from scipy.stats import invwishart
from numpy.linalg import solve as solve
from numpy.linalg import cholesky as cholesky_lower
from scipy.linalg import cholesky as cholesky_upper
from scipy.linalg import solve_triangular as solve_ut
import matplotlib.pyplot as plt
%matplotlib inline
def mvnrnd_pre(mu, Lambda):
    src = normrnd(size = (mu.shape[0],))
    return solve_ut(cholesky_upper(Lambda, overwrite_a = True, check_finite = False), 
                    src, lower = False, check_finite = False, overwrite_b = True) + mu

def cov_mat(mat, mat_bar):
    mat = mat - mat_bar
    return mat.T @ mat

def ten2mat(tensor, mode):
    return np.reshape(np.moveaxis(tensor, mode, 0), (tensor.shape[mode], -1), order = 'F')
def sample_factor_u(tau_sparse_tensor, tau_ind, U, V, X, beta0 = 1):
    """Sampling M-by-R factor matrix U and its hyperparameters (mu_u, Lambda_u)."""
    
    dim1, rank = U.shape
    U_bar = np.mean(U, axis = 0)
    temp = dim1 / (dim1 + beta0)
    var_mu_hyper = temp * U_bar
    var_U_hyper = inv(np.eye(rank) + cov_mat(U, U_bar) + temp * beta0 * np.outer(U_bar, U_bar))
    var_Lambda_hyper = wishart.rvs(df = dim1 + rank, scale = var_U_hyper)
    var_mu_hyper = mvnrnd_pre(var_mu_hyper, (dim1 + beta0) * var_Lambda_hyper)

    var1 = kr_prod(X, V).T
    var2 = kr_prod(var1, var1)
    var3 = (var2 @ ten2mat(tau_ind, 0).T).reshape([rank, rank, dim1]) + var_Lambda_hyper[:, :, None]
    var4 = var1 @ ten2mat(tau_sparse_tensor, 0).T + (var_Lambda_hyper @ var_mu_hyper)[:, None]
    for i in range(dim1):
        U[i, :] = mvnrnd_pre(solve(var3[:, :, i], var4[:, i]), var3[:, :, i])
        
    return U
def sample_factor_v(tau_sparse_tensor, tau_ind, U, V, X, beta0 = 1):
    """Sampling N-by-R factor matrix V and its hyperparameters (mu_v, Lambda_v)."""
    
    dim2, rank = V.shape
    V_bar = np.mean(V, axis = 0)
    temp = dim2 / (dim2 + beta0)
    var_mu_hyper = temp * V_bar
    var_V_hyper = inv(np.eye(rank) + cov_mat(V, V_bar) + temp * beta0 * np.outer(V_bar, V_bar))
    var_Lambda_hyper = wishart.rvs(df = dim2 + rank, scale = var_V_hyper)
    var_mu_hyper = mvnrnd_pre(var_mu_hyper, (dim2 + beta0) * var_Lambda_hyper)

    var1 = kr_prod(X, U).T
    var2 = kr_prod(var1, var1)
    var3 = (var2 @ ten2mat(tau_ind, 1).T).reshape([rank, rank, dim2]) + var_Lambda_hyper[:, :, None]
    var4 = var1 @ ten2mat(tau_sparse_tensor, 1).T + (var_Lambda_hyper @ var_mu_hyper)[:, None]
    for j in range(dim2):
        V[j, :] = mvnrnd_pre(solve(var3[:, :, j], var4[:, j]), var3[:, :, j])
        
    return V
def mnrnd(M, U, V):
    """
    Generate matrix normal distributed random matrix.
    M is a m-by-n matrix, U is a m-by-m matrix, and V is a n-by-n matrix.
    """
    dim1, dim2 = M.shape
    X0 = np.random.randn(dim1, dim2)
    P = cholesky_lower(U)
    Q = cholesky_lower(V)
    
    return M + P @ X0 @ Q.T

def sample_var_coefficient(X, time_lags):
    dim, rank = X.shape
    d = time_lags.shape[0]
    tmax = np.max(time_lags)
    
    Z_mat = X[tmax : dim, :]
    Q_mat = np.zeros((dim - tmax, rank * d))
    for k in range(d):
        Q_mat[:, k * rank : (k + 1) * rank] = X[tmax - time_lags[k] : dim - time_lags[k], :]
    var_Psi0 = np.eye(rank * d) + Q_mat.T @ Q_mat
    var_Psi = inv(var_Psi0)
    var_M = var_Psi @ Q_mat.T @ Z_mat
    var_S = np.eye(rank) + Z_mat.T @ Z_mat - var_M.T @ var_Psi0 @ var_M
    Sigma = invwishart.rvs(df = rank + dim - tmax, scale = var_S)
    
    return mnrnd(var_M, var_Psi, Sigma), Sigma
def sample_factor_x(tau_sparse_tensor, tau_ind, time_lags, U, V, X, A, Lambda_x):
    """Sampling T-by-R factor matrix X."""

    dim3, rank = X.shape
    tmax = np.max(time_lags)
    tmin = np.min(time_lags)
    d = time_lags.shape[0]
    A0 = np.dstack([A] * d)
    for k in range(d):
        A0[k * rank : (k + 1) * rank, :, k] = 0
    mat0 = Lambda_x @ A.T
    mat1 = np.einsum('kij, jt -> kit', A.reshape([d, rank, rank]), Lambda_x)
    mat2 = np.einsum('kit, kjt -> ij', mat1, A.reshape([d, rank, rank]))
    
    var1 = kr_prod(V, U).T
    var2 = kr_prod(var1, var1)
    var3 = (var2 @ ten2mat(tau_ind, 2).T).reshape([rank, rank, dim3]) + Lambda_x[:, :, None]
    var4 = var1 @ ten2mat(tau_sparse_tensor, 2).T
    for t in range(dim3):
        Mt = np.zeros((rank, rank))
        Nt = np.zeros(rank)
        Qt = mat0 @ X[t - time_lags, :].reshape(rank * d)
        index = list(range(0, d))
        if t >= dim3 - tmax and t < dim3 - tmin:
            index = list(np.where(t + time_lags < dim3))[0]
        elif t < tmax:
            Qt = np.zeros(rank)
            index = list(np.where(t + time_lags >= tmax))[0]
        if t < dim3 - tmin:
            Mt = mat2.copy()
            temp = np.zeros((rank * d, len(index)))
            n = 0
            for k in index:
                temp[:, n] = X[t + time_lags[k] - time_lags, :].reshape(rank * d)
                n += 1
            temp0 = X[t + time_lags[index], :].T - np.einsum('ijk, ik -> jk', A0[:, :, index], temp)
            Nt = np.einsum('kij, jk -> i', mat1[index, :, :], temp0)
        
        var3[:, :, t] = var3[:, :, t] + Mt
        if t < tmax:
            var3[:, :, t] = var3[:, :, t] - Lambda_x + np.eye(rank)
        X[t, :] = mvnrnd_pre(solve(var3[:, :, t], var4[:, t] + Nt + Qt), var3[:, :, t])

    return X
def sample_precision_tau(sparse_tensor, tensor_hat, ind):
    var_alpha = 1e-6 + 0.5 * np.sum(ind, axis = 2)
    var_beta = 1e-6 + 0.5 * np.sum(((sparse_tensor - tensor_hat) ** 2) * ind, axis = 2)
    return np.random.gamma(var_alpha, 1 / var_beta)
def compute_mape(var, var_hat):
    return np.sum(np.abs(var - var_hat) / var) / var.shape[0]

def compute_rmse(var, var_hat):
    return  np.sqrt(np.sum((var - var_hat) ** 2) / var.shape[0])
def BTTF(dense_tensor, sparse_tensor, init, rank, time_lags, burn_iter, gibbs_iter, multi_steps = 1, vargin = 0):
    """Bayesian Temporal Tensor Factorization, BTTF."""
    
    dim1, dim2, dim3 = sparse_tensor.shape
    d = time_lags.shape[0]
    U = init["U"]
    V = init["V"]
    X = init["X"]
    if np.isnan(sparse_tensor).any() == False:
        ind = sparse_tensor != 0
        pos_obs = np.where(ind)
        pos_test = np.where((dense_tensor != 0) & (sparse_tensor == 0))
    elif np.isnan(sparse_tensor).any() == True:
        pos_test = np.where((dense_tensor != 0) & (np.isnan(sparse_tensor)))
        ind = ~np.isnan(sparse_tensor)
        pos_obs = np.where(ind)
        sparse_tensor[np.isnan(sparse_tensor)] = 0
    dense_test = dense_tensor[pos_test]
    del dense_tensor
    U_plus = np.zeros((dim1, rank))
    V_plus = np.zeros((dim2, rank))
    X_new_plus = np.zeros((dim3 + multi_steps, rank))
    A_plus = np.zeros((rank * d, rank))
    temp_hat = np.zeros(sparse_tensor.shape)
    show_iter = 200
    if vargin == 0: # scalar tau
        tau = 1
    elif vargin == 1: # matrix tau
        tau = np.ones((dim1, dim2))
    tensor_hat_plus = np.zeros(sparse_tensor.shape)
    for it in range(burn_iter + gibbs_iter):
        if vargin == 0: # scalar tau
            tau_ind = tau * ind
            tau_sparse_tensor = tau * sparse_tensor
            U = sample_factor_u(tau_sparse_tensor, tau_ind, U, V, X)
            V = sample_factor_v(tau_sparse_tensor, tau_ind, U, V, X)
            A, Sigma = sample_var_coefficient(X, time_lags)
            X = sample_factor_x(tau_sparse_tensor, tau_ind, time_lags, U, V, X, A, inv(Sigma))
            tensor_hat = np.einsum('is, js, ts -> ijt', U, V, X)
            tau = np.random.gamma(1e-6 + 0.5 * np.sum(ind), 
                                  1 / (1e-6 + 0.5 * np.sum(((sparse_tensor - tensor_hat) ** 2) * ind)))
        elif vargin == 1: # matrix tau
            tau_ind = tau[:, :, None] * ind
            tau_sparse_tensor = tau[:, :, None] * sparse_tensor
            U = sample_factor_u(tau_sparse_tensor, tau_ind, U, V, X)
            V = sample_factor_v(tau_sparse_tensor, tau_ind, U, V, X)
            A, Sigma = sample_var_coefficient(X, time_lags)
            X = sample_factor_x(tau_sparse_tensor, tau_ind, time_lags, U, V, X, A, inv(Sigma))
            tensor_hat = np.einsum('is, js, ts -> ijt', U, V, X)
            tau = sample_precision_tau(sparse_tensor, tensor_hat, ind)
        temp_hat += tensor_hat
        if (it + 1) % show_iter == 0 and it < burn_iter:
            temp_hat = temp_hat / show_iter
            print('Iter: {}'.format(it + 1))
            print('MAPE: {:.6}'.format(compute_mape(dense_test, temp_hat[pos_test])))
            print('RMSE: {:.6}'.format(compute_rmse(dense_test, temp_hat[pos_test])))
            temp_hat = np.zeros(sparse_tensor.shape)
            print()
        X_new = np.zeros((dim3 + multi_steps, rank))
        if it + 1 > burn_iter:
            U_plus += U
            V_plus += V
            A_plus += A
            X_new[: dim3, :] = X.copy()
            if multi_steps == 1:
                X_new[dim3, :] = A.T @ X_new[dim3 - time_lags, :].reshape(rank * d)
            elif multi_steps > 1:
                for t0 in range(multi_steps):
                    X_new[dim3 + t0, :] = A.T @ X_new[dim3 + t0 - time_lags, :].reshape(rank * d)
            X_new_plus += X_new
            tensor_hat_plus += tensor_hat
    tensor_hat = tensor_hat_plus / gibbs_iter
    U = U_plus / gibbs_iter
    V = V_plus / gibbs_iter
    X_new = X_new_plus / gibbs_iter
    A = A_plus / gibbs_iter
    print('Imputation MAPE: {:.6}'.format(compute_mape(dense_test, tensor_hat[pos_test])))
    print('Imputation RMSE: {:.6}'.format(compute_rmse(dense_test, tensor_hat[pos_test])))
    print()
    
    return tensor_hat, U, V, X_new, A<jupyter_output><empty_output><jupyter_text>Evaluation on New York Taxi Data**Scenario setting**:- Tensor size: $30\times 30\times 1464$ (origin, destination, time)- Random missing (RM)- 40% missing rate<jupyter_code>import scipy.io

dense_tensor = scipy.io.loadmat('../datasets/NYC-data-set/tensor.mat')['tensor'].astype(np.float32)
rm_tensor = scipy.io.loadmat('../datasets/NYC-data-set/rm_tensor.mat')['rm_tensor']
missing_rate = 0.4

## Random missing (RM)
binary_tensor = np.round(rm_tensor + 0.5 - missing_rate)
sparse_tensor = dense_tensor.copy()
sparse_tensor[binary_tensor == 0] = np.nan<jupyter_output><empty_output><jupyter_text>**Model setting**:- Low rank: 30- Time lags: {1, 2, 24}- The number of burn-in iterations: 1000- The number of Gibbs iterations: 200<jupyter_code>import time
start = time.time()
dim1, dim2, dim3 = sparse_tensor.shape
rank = 30
time_lags = np.array([1, 2, 24])
init = {"U": 0.1 * np.random.randn(dim1, rank), "V": 0.1 * np.random.randn(dim2, rank),
        "X": 0.1 * np.random.randn(dim3, rank)}
burn_iter = 1000
gibbs_iter = 200
tensor_hat, U, V, X_new, A = BTTF(dense_tensor, sparse_tensor, init, rank, time_lags, burn_iter, gibbs_iter)
end = time.time()
print('Running time: %d seconds'%(end - start))<jupyter_output>Iter: 200
MAPE: 0.479829
RMSE: 4.77438

Iter: 400
MAPE: 0.478661
RMSE: 4.83144

Iter: 600
MAPE: 0.477919
RMSE: 4.83232

Iter: 800
MAPE: 0.478
RMSE: 4.8425

Iter: 1000
MAPE: 0.478093
RMSE: 4.84608

Imputation MAPE: 0.478035
Imputation RMSE: 4.85489

Running time: 1103 seconds<jupyter_text>**Scenario setting**:- Tensor size: $30\times 30\times 1464$ (origin, destination, time)- Random missing (RM)- 60% missing rate<jupyter_code>import scipy.io

dense_tensor = scipy.io.loadmat('../datasets/NYC-data-set/tensor.mat')['tensor'].astype(np.float32)
rm_tensor = scipy.io.loadmat('../datasets/NYC-data-set/rm_tensor.mat')['rm_tensor']
missing_rate = 0.6

## Random missing (RM)
binary_tensor = np.round(rm_tensor + 0.5 - missing_rate)
sparse_tensor = dense_tensor.copy()
sparse_tensor[binary_tensor == 0] = np.nan<jupyter_output><empty_output><jupyter_text>**Model setting**:- Low rank: 30- Time lags: {1, 2, 24}- The number of burn-in iterations: 1000- The number of Gibbs iterations: 200<jupyter_code>import time
start = time.time()
dim1, dim2, dim3 = sparse_tensor.shape
rank = 30
time_lags = np.array([1, 2, 24])
init = {"U": 0.1 * np.random.randn(dim1, rank), "V": 0.1 * np.random.randn(dim2, rank),
        "X": 0.1 * np.random.randn(dim3, rank)}
burn_iter = 1000
gibbs_iter = 200
tensor_hat, U, V, X_new, A = BTTF(dense_tensor, sparse_tensor, init, rank, time_lags, burn_iter, gibbs_iter)
end = time.time()
print('Running time: %d seconds'%(end - start))<jupyter_output>Iter: 200
MAPE: 0.485985
RMSE: 4.94286

Iter: 400
MAPE: 0.48869
RMSE: 5.04894

Iter: 600
MAPE: 0.488934
RMSE: 5.05074

Iter: 800
MAPE: 0.488713
RMSE: 5.06594

Iter: 1000
MAPE: 0.489163
RMSE: 5.07669

Imputation MAPE: 0.48909
Imputation RMSE: 5.08303

Running time: 1102 seconds<jupyter_text>**Scenario setting**:- Tensor size: $30\times 30\times 1464$ (origin, destination, time)- Non-random missing (NM)- 40% missing rate<jupyter_code>import scipy.io

dense_tensor = scipy.io.loadmat('../datasets/NYC-data-set/tensor.mat')['tensor'].astype(np.float32)
nm_tensor = scipy.io.loadmat('../datasets/NYC-data-set/nm_tensor.mat')['nm_tensor']
missing_rate = 0.4

## Non-random missing (NM)
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 = dense_tensor.copy()
sparse_tensor[binary_tensor == 0] = np.nan<jupyter_output><empty_output><jupyter_text>**Model setting**:- Low rank: 30- Time lags: {1, 2, 24}- The number of burn-in iterations: 1000- The number of Gibbs iterations: 200<jupyter_code>import time
start = time.time()
dim1, dim2, dim3 = sparse_tensor.shape
rank = 30
time_lags = np.array([1, 2, 24])
init = {"U": 0.1 * np.random.randn(dim1, rank), "V": 0.1 * np.random.randn(dim2, rank),
        "X": 0.1 * np.random.randn(dim3, rank)}
burn_iter = 1000
gibbs_iter = 200
tensor_hat, U, V, X_new, A = BTTF(dense_tensor, sparse_tensor, init, rank, time_lags, burn_iter, gibbs_iter)
end = time.time()
print('Running time: %d seconds'%(end - start))
import numpy as np
import scipy.io
np.random.seed(1000)

dense_tensor = scipy.io.loadmat('../datasets/NYC-data-set/tensor.mat')['tensor'].astype(np.float32)
binary_tensor = np.zeros(dense_tensor.shape)
dim1, dim2, dim3 = dense_tensor.shape
missing_rate = 0.2

## Fiber missing (FM)
binary = np.round(np.random.rand(dim1, dim2) + 0.5 - missing_rate)
for i in range(dim1):
    dense_tensor[i, i, :] = 0
    binary_tensor[i, i, :] = 0
binary_tensor = binary[:, :, None] * np.ones((dim1, dim2, dim3))
sparse_tensor = dense_tensor.copy()
sparse_tensor[binary_tensor == 0] = np.nan
import time
start = time.time()
dim1, dim2, dim3 = sparse_tensor.shape
rank = 10
time_lags = np.array([1, 2, 24])
init = {"U": 0.1 * np.random.randn(dim1, rank), "V": 0.1 * np.random.randn(dim2, rank),
        "X": 0.1 * np.random.randn(dim3, rank)}
burn_iter = 1000
gibbs_iter = 200
tensor_hat, U, V, X_new, A = BTTF(dense_tensor, sparse_tensor, init, rank, time_lags, burn_iter, gibbs_iter)
end = time.time()
print('Running time: %d seconds'%(end - start))<jupyter_output>Iter: 200
MAPE: 2.62378
RMSE: 36.1879

Iter: 400
MAPE: 4.38537
RMSE: 75.5418

Iter: 600
MAPE: 11.6644
RMSE: 190.616

Iter: 800
MAPE: 20.1019
RMSE: 314.556

Iter: 1000
MAPE: 24.7993
RMSE: 394.532

Imputation MAPE: 27.1575
Imputation RMSE: 435.979

Running time: 452 seconds<jupyter_text>Evaluation on Pacific Surface Temperature Data**Scenario setting**:- Tensor size: $30\times 84\times 396$ (location x, location y, month)- Random missing (RM)- 40% missing rate<jupyter_code>import numpy as np
np.random.seed(1000)

dense_tensor = np.load('../datasets/Temperature-data-set/tensor.npy').astype(np.float32)
pos = np.where(dense_tensor[:, 0, :] > 50)
dense_tensor[pos[0], :, pos[1]] = 0
random_tensor = np.random.rand(dense_tensor.shape[0], dense_tensor.shape[1], dense_tensor.shape[2])
missing_rate = 0.4

## Random missing (RM)
binary_tensor = np.round(random_tensor + 0.5 - missing_rate)
sparse_tensor = dense_tensor.copy()
sparse_tensor[binary_tensor == 0] = np.nan
sparse_tensor[sparse_tensor == 0] = np.nan<jupyter_output><empty_output><jupyter_text>**Model setting**:- Low rank: 20- Time lags: {1, 2, 12}- The number of burn-in iterations: 1000- The number of Gibbs iterations: 200<jupyter_code>import time
start = time.time()
dim1, dim2, dim3 = sparse_tensor.shape
rank = 30
time_lags = np.array([1, 2, 12])
init = {"U": 0.1 * np.random.randn(dim1, rank), "V": 0.1 * np.random.randn(dim2, rank),
        "X": 0.1 * np.random.randn(dim3, rank)}
burn_iter = 1000
gibbs_iter = 200
tensor_hat, U, V, X_new, A = BTTF(dense_tensor, sparse_tensor, init, rank, time_lags, burn_iter, gibbs_iter)
end = time.time()
print('Running time: %d seconds'%(end - start))<jupyter_output>Iter: 200
MAPE: 0.015674
RMSE: 0.524397

Iter: 400
MAPE: 0.0152527
RMSE: 0.509142

Iter: 600
MAPE: 0.0151012
RMSE: 0.503865

Iter: 800
MAPE: 0.0150042
RMSE: 0.500957

Iter: 1000
MAPE: 0.0148908
RMSE: 0.496983

Imputation MAPE: 0.014815
Imputation RMSE: 0.494163

Running time: 483 seconds<jupyter_text>**Scenario setting**:- Tensor size: $30\times 84\times 396$ (location x, location y, month)- Random missing (RM)- 60% missing rate<jupyter_code>import numpy as np
np.random.seed(1000)

dense_tensor = np.load('../datasets/Temperature-data-set/tensor.npy').astype(np.float32)
pos = np.where(dense_tensor[:, 0, :] > 50)
dense_tensor[pos[0], :, pos[1]] = 0
random_tensor = np.random.rand(dense_tensor.shape[0], dense_tensor.shape[1], dense_tensor.shape[2])
missing_rate = 0.6

## Random missing (RM)
binary_tensor = np.round(random_tensor + 0.5 - missing_rate)
sparse_tensor = dense_tensor.copy()
sparse_tensor[binary_tensor == 0] = np.nan
sparse_tensor[sparse_tensor == 0] = np.nan<jupyter_output><empty_output><jupyter_text>**Model setting**:- Low rank: 20- Time lags: {1, 2, 12}- The number of burn-in iterations: 1000- The number of Gibbs iterations: 200<jupyter_code>import time
start = time.time()
dim1, dim2, dim3 = sparse_tensor.shape
rank = 30
time_lags = np.array([1, 2, 12])
init = {"U": 0.1 * np.random.randn(dim1, rank), "V": 0.1 * np.random.randn(dim2, rank),
        "X": 0.1 * np.random.randn(dim3, rank)}
burn_iter = 1000
gibbs_iter = 200
tensor_hat, U, V, X_new, A = BTTF(dense_tensor, sparse_tensor, init, rank, time_lags, burn_iter, gibbs_iter)
end = time.time()
print('Running time: %d seconds'%(end - start))<jupyter_output>Iter: 200
MAPE: 0.0161395
RMSE: 0.542084

Iter: 400
MAPE: 0.0155541
RMSE: 0.521575

Iter: 600
MAPE: 0.0152756
RMSE: 0.51217

Iter: 800
MAPE: 0.0150673
RMSE: 0.505096

Iter: 1000
MAPE: 0.0149478
RMSE: 0.501089

Imputation MAPE: 0.014884
Imputation RMSE: 0.499067

Running time: 478 seconds<jupyter_text>**Scenario setting**:- Tensor size: $30\times 84\times 396$ (location x, location y, month)- Non-random missing (NM)- 40% missing rate<jupyter_code>import numpy as np
np.random.seed(1000)

dense_tensor = np.load('../datasets/Temperature-data-set/tensor.npy').astype(np.float32)
pos = np.where(dense_tensor[:, 0, :] > 50)
dense_tensor[pos[0], :, pos[1]] = 0
random_tensor = np.random.rand(dense_tensor.shape[0], dense_tensor.shape[1], int(dense_tensor.shape[2] / 3))
missing_rate = 0.4

## Non-random missing (NM)
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(int(dense_tensor.shape[2] / 3)):
            binary_tensor[i1, i2, i3 * 3 : (i3 + 1) * 3] = np.round(random_tensor[i1, i2, i3] + 0.5 - missing_rate)
sparse_tensor = dense_tensor.copy()
sparse_tensor[binary_tensor == 0] = np.nan
sparse_tensor[sparse_tensor == 0] = np.nan
import time
start = time.time()
dim1, dim2, dim3 = sparse_tensor.shape
rank = 30
time_lags = np.array([1, 2, 12])
init = {"U": 0.1 * np.random.randn(dim1, rank), "V": 0.1 * np.random.randn(dim2, rank),
        "X": 0.1 * np.random.randn(dim3, rank)}
burn_iter = 1000
gibbs_iter = 200
tensor_hat, U, V, X_new, A = BTTF(dense_tensor, sparse_tensor, init, rank, time_lags, burn_iter, gibbs_iter)
end = time.time()
print('Running time: %d seconds'%(end - start))<jupyter_output>Iter: 200
MAPE: 0.0158287
RMSE: 0.531278

Iter: 400
MAPE: 0.0153527
RMSE: 0.514709

Iter: 600
MAPE: 0.0150346
RMSE: 0.501843

Iter: 800
MAPE: 0.0148943
RMSE: 0.496861

Iter: 1000
MAPE: 0.014797
RMSE: 0.493559

Imputation MAPE: 0.0147441
Imputation RMSE: 0.491632

Running time: 729 seconds
{'fraction_non_alphanumeric': 0.14843859403444895, 'fraction_numerical': 0.058581897960136305, 'mean_word_length': 3.5179249261914802, 'pattern_counts': {'":': 21, '<': 51, '<?xml version=': 0, '>': 58, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '99824', 'license': 'MIT', 'n_tokens_mistral': 9324, 'n_tokens_neox': 8380, 'n_words': 2037, 'path': 'imputer/BTTF.ipynb', 'repo_name': 'shawnwang-tech/transdim', 'size': '21430'}