Repository: WenjieDu/PyGrinder (language: python)

Issue:
Add MNAR and MAR patterns
### 1. Feature description

MNAR (missing not at random), missing values depend on both observed values and missing values;

MAR (missing at random), missing values only depend on observed parts.

### 2. Motivation

PyGrinder so far only supports MCAR pattern. MNAR and MAR both should be added as well.

### 3. Your contribution

Will create a PR to add them.
Add package source to PyPI
Only `.whl` file is available on PyPI. Please also make available the source (`*.tar.gz`) file.


Discussion / hints:

Thank you, @sugatoray 😃 It's on my list!
Sorry for the late update. It's already been solved by using [the auto-publishing workflow](https://github.com/WenjieDu/PyCorruptor/blob/main/.github/workflows/publish_to_PyPI.yml), just like PyPOTS and TSDB.

Fix patch:
diff --git a/pygrinder/__init__.py b/pygrinder/__init__.py
index 91956b2..78d0c21 100644
--- a/pygrinder/__init__.py
+++ b/pygrinder/__init__.py
@@ -25,6 +25,7 @@
 
 try:
     from pygrinder.mcar import mcar
+    from pygrinder.mnar import mnar_x, mnar_t
     from pygrinder.utils import (
         cal_missing_rate,
         masked_fill,
@@ -38,4 +39,6 @@
     "cal_missing_rate",
     "masked_fill",
     "mcar",
+    "mnar_x"
+    "mnar_t"
 ]
diff --git a/pygrinder/mar.py b/pygrinder/mar.py
index a5d6610..f37d57b 100644
--- a/pygrinder/mar.py
+++ b/pygrinder/mar.py
@@ -5,39 +5,152 @@
 # Created by Wenjie Du <wenjay.du@gmail.com>
 # License: GLP-v3
 
+from typing import Union, Tuple
 
-def mar(X, rate, nan=0):
-    """Create random missing values (MAR case).
+import numpy as np
+import torch
+from scipy import optimize
+
+
+def mar_logistic(
+    X: Union[torch.Tensor, np.ndarray],
+    obs_rate: float,
+    missing_rate: float,
+    nan: Union[float, int] = 0,
+) -> Union[
+    Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray],
+    Tuple[torch.Tensor, torch.Tensor, torch.Tensor, torch.Tensor],
+]:
+    """Create random missing values (MAR case) with a logistic model.
+    First, a subset of the variables without missing values is randomly selected.
+    Missing values will be introduced into the remaining variables according to a logistic model with random weights.
+    This implementation is inspired by the tutorial
+    https://rmisstastic.netlify.app/how-to/python/generate_html/how%20to%20generate%20missing%20values
 
     Parameters
     ----------
-    X : array,
-        Data vector. If X has any missing values, they should be numpy.nan.
+    X : shape of [n_steps, n_features]
+        A time series data vector without any missing data.
 
-    rate : float, in (0,1),
-        Artificially missing rate, rate of the observed values which will be artificially masked as missing.
+    obs_rate :
+        The proportion of variables without missing values that will be used for fitting the logistic masking model.
 
-        Note that,
-        `rate` = (number of artificially missing values) / np.sum(~np.isnan(self.data)),
-        not (number of artificially missing values) / np.product(self.data.shape),
-        considering that the given data may already contain missing values,
-        the latter way may be confusing because if the original missing rate >= `rate`,
-        the function will do nothing, i.e. it won't play the role it has to be.
+    missing_rate:
+        The proportion of missing values to generate for variables which will have missing values.
 
-    nan : int/float, optional, default=0
+    nan :
         Value used to fill NaN values.
 
     Returns
     -------
+    X_intact : array,
+        Original data with missing values (nan) filled with given parameter `nan`, with observed values intact.
+        X_intact is for loss calculation in the masked imputation task.
+
+    X : array,
+        Original X with artificial missing values. X is for model input.
+        Both originally-missing and artificially-missing values are filled with given parameter `nan`.
+
+    missing_mask : array,
+        The mask indicates all missing values in X.
+        In it, 1 indicates observed values, and 0 indicates missing values.
+
+    indicating_mask : array,
+        The mask indicates the artificially-missing values in X, namely missing parts different from X_intact.
+        In it, 1 indicates artificially missing values,
+        and the other values (including originally observed/missing values) are indicated as 0.
 
     """
-    # TODO: Create missing values in MAR case
-    raise NotImplementedError("MAR case has not been implemented yet.")
+    if isinstance(X, list):
+        X = np.asarray(X)
+
+    if isinstance(X, np.ndarray) or isinstance(X, torch.Tensor):
+        X_intact, X, missing_mask, indicating_mask = _mar_logistic_torch(
+            X, obs_rate, missing_rate, nan
+        )
+    else:
+        raise TypeError(
+            "X must be type of list/numpy.ndarray/torch.Tensor, " f"but got {type(X)}"
+        )
+    return X_intact, X, missing_mask, indicating_mask
+
+
+def _mar_logistic_torch(X, rate, rate_obs, nan):
+    def pick_coefficients(X, idxs_obs=None, idxs_nas=None, self_mask=False):
+        n, d = X.shape
+        if self_mask:
+            coeffs = torch.randn(d)
+            Wx = X * coeffs
+            coeffs /= torch.std(Wx, 0)
+        else:
+            d_obs = len(idxs_obs)
+            d_na = len(idxs_nas)
+            coeffs = torch.randn(d_obs, d_na)
+            Wx = X[:, idxs_obs].mm(coeffs)
+            coeffs /= torch.std(Wx, 0, keepdim=True)
+        return coeffs
+
+    def fit_intercepts(X, coeffs, p, self_mask=False):
+        if self_mask:
+            d = len(coeffs)
+            intercepts = torch.zeros(d)
+            for j in range(d):
+
+                def f(x):
+                    return torch.sigmoid(X * coeffs[j] + x).mean().item() - p
+
+                intercepts[j] = optimize.bisect(f, -50, 50)
+        else:
+            d_obs, d_na = coeffs.shape
+            intercepts = torch.zeros(d_na)
+            for j in range(d_na):
+
+                def f(x):
+                    return torch.sigmoid(X.mv(coeffs[:, j]) + x).mean().item() - p
+
+                intercepts[j] = optimize.bisect(f, -50, 50)
+        return intercepts
+
+    assert len(X.shape) == 2, "X should be 2 dimensional"
+    n, d = X.shape
+
+    ori_type_is_np = isinstance(X, np.ndarray)
+    if ori_type_is_np:
+        X = torch.from_numpy(X).to(torch.float32)
+    else:
+        X = torch.clone(X).to(torch.float32)
+
+    assert (
+        torch.isnan(X).sum() == 0
+    ), "the input X of the mar_logistic() shouldn't containing originally missing data"
+
+    X_intact = torch.clone(X)
+    mask = torch.zeros(n, d).bool()
+
+    # number of variables that will have no missing values (at least one variable)
+    d_obs = max(int(rate_obs * d), 1)
+    d_na = d - d_obs  # number of variables that will have missing values
+
+    # Sample variables will all be observed, and the left will be with missing values
+    idxs_obs = np.random.choice(d, d_obs, replace=False)
+    idxs_nas = np.array([i for i in range(d) if i not in idxs_obs])
 
+    # Pick coefficients so that W^Tx has unit variance (avoids shrinking)
+    coeffs = pick_coefficients(X, idxs_obs, idxs_nas)
+    # Pick the intercepts to have a desired amount of missing values
+    intercepts = fit_intercepts(X[:, idxs_obs], coeffs, rate)
 
-def _mar_numpy(X, rate, nan=0):
-    pass
+    ps = torch.sigmoid(X[:, idxs_obs].mm(coeffs) + intercepts)
+    ber = torch.rand(n, d_na)
+    mask[:, idxs_nas] = ber < ps  # True means missing
 
+    X[mask] = torch.nan
+    X = torch.nan_to_num(X, nan)
+    missing_mask = (~mask).to(torch.float32)
+    indicating_mask = torch.clone(missing_mask)
 
-def _mar_torch(X, rate, nan=0):
-    pass
+    if ori_type_is_np:
+        X_intact, X, missing_mask, indicating_mask = [
+            i.numpy() for i in (X_intact, X, missing_mask, indicating_mask)
+        ]
+    return X_intact, X, missing_mask, indicating_mask
diff --git a/pygrinder/mnar.py b/pygrinder/mnar.py
index b005312..222483f 100644
--- a/pygrinder/mnar.py
+++ b/pygrinder/mnar.py
@@ -2,42 +2,268 @@
 Corrupt data by adding missing values to it with MNAR (missing not at random) pattern.
 """
 
-# Created by Wenjie Du <wenjay.du@gmail.com>
+# Created by Jun Wang <jwangfx@connect.ust.hk> and Wenjie Du <wenjay.du@gmail.com>
 # License: GLP-v3
 
+from typing import Optional, Union, Tuple
 
-def mnar(X, rate, nan=0):
-    """Create not-random missing values (MNAR case).
+import numpy as np
+
+try:
+    import torch
+except ImportError:
+    pass
+
+
+def mnar_x(
+    X: Optional[Union[np.ndarray, torch.Tensor]],
+    offset: float = 0,
+    nan: Union[float, int] = 0,
+) -> Union[
+    Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray],
+    Tuple[torch.Tensor, torch.Tensor, torch.Tensor, torch.Tensor],
+]:
+    """Create not-random missing values related to values themselves (MNAR-x case ot self-masking MNAR case).
+    This case follows the setting in Ipsen et al. "not-MIWAE: Deep Generative Modelling with Missing Not at Random Data"
+    :cite`ipsen2021notmiwae`.
 
     Parameters
     ----------
+    X :
+        Data vector. If X has any missing values, they should be numpy.nan.
+
+    offset :
+        the weight of standard deviation. In MNAR-x case, for each time series,
+        the values larger than the mean of each time series plus offset*standard deviation will be missing
+
+    nan :
+        Value used to fill NaN values.
+
+    Returns
+    -------
+    X_intact : array,
+        Original data with missing values (nan) filled with given parameter `nan`, with observed values intact.
+        X_intact is for loss calculation in the masked imputation task.
+
     X : array,
+        Original X with artificial missing values. X is for model input.
+        Both originally-missing and artificially-missing values are filled with given parameter `nan`.
+
+    missing_mask : array,
+        The mask indicates all missing values in X.
+        In it, 1 indicates observed values, and 0 indicates missing values.
+
+    indicating_mask : array,
+        The mask indicates the artificially-missing values in X, namely missing parts different from X_intact.
+        In it, 1 indicates artificially missing values,
+        and the other values (including originally observed/missing values) are indicated as 0.
+
+    """
+    if isinstance(X, list):
+        X = np.asarray(X)
+
+    if isinstance(X, np.ndarray):
+        X_intact, X, missing_mask, indicating_mask = _mnar_x_numpy(X, offset, nan)
+    elif isinstance(X, torch.Tensor):
+        X_intact, X, missing_mask, indicating_mask = _mnar_x_torch(X, offset, nan)
+    else:
+        raise TypeError(
+            "X must be type of list/numpy.ndarray/torch.Tensor, " f"but got {type(X)}"
+        )
+    return X_intact, X, missing_mask, indicating_mask
+
+
+def _mnar_x_numpy(
+    X: np.ndarray,
+    offset: float = 0,
+    nan: Union[float, int] = 0,
+) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
+    # clone X to ensure values of X out of this function not being affected
+    X = np.copy(X)
+
+    X_intact = np.copy(X)  # keep a copy of originally observed values in X_intact
+
+    n_s, n_l, n_c = X.shape
+
+    ori_mask = ~np.isnan(X)
+    X = np.nan_to_num(X, nan=nan)
+    X_intact = np.nan_to_num(X_intact, nan=nan)
+
+    mask_sum = ori_mask.sum(1)
+    mask_sum[mask_sum == 0] = 1
+    X_mean = np.repeat(
+        ((X * ori_mask).sum(1) / mask_sum).reshape(n_s, 1, n_c), n_l, axis=1
+    )
+    X_std = np.repeat(
+        np.sqrt(np.square((X - X_mean) * ori_mask).sum(1) / mask_sum).reshape(
+            n_s, 1, n_c
+        ),
+        n_l,
+        axis=1,
+    )
+
+    mnar_missing_mask = np.zeros_like(X)
+    mnar_missing_mask[X <= (X_mean + offset * X_std)] = 1
+    missing_mask = ori_mask * mnar_missing_mask
+    indicating_mask = ori_mask - missing_mask
+    X[missing_mask == 0] = nan
+    return X_intact, X, missing_mask, indicating_mask
+
+
+def _mnar_x_torch(
+    X: torch.Tensor,
+    offset: float = 0,
+    nan: Union[float, int] = 0,
+) -> Tuple[torch.Tensor, torch.Tensor, torch.Tensor, torch.Tensor]:
+    # clone X to ensure values of X out of this function not being affected
+    X = torch.clone(X)
+
+    X_intact = torch.clone(X)  # keep a copy of originally observed values in X_intact
+    n_s, n_l, n_c = X.shape
+
+    ori_mask = (~torch.isnan(X)).type(torch.float32)
+    X = torch.nan_to_num(X, nan=nan)
+    X_intact = torch.nan_to_num(X_intact, nan=nan)
+
+    mask_sum = ori_mask.sum(1)
+    mask_sum[mask_sum == 0] = 1
+    X_mean = ((X * ori_mask).sum(1) / mask_sum).reshape(n_s, 1, n_c).repeat(1, n_l, 1)
+    X_std = (
+        (((X - X_mean) * ori_mask).pow(2).sum(1) / mask_sum)
+        .sqrt()
+        .reshape(n_s, 1, n_c)
+        .repeat(1, n_l, 1)
+    )
+
+    mnar_missing_mask = torch.zeros_like(X)
+    mnar_missing_mask[X <= (X_mean + offset * X_std)] = 1
+    missing_mask = ori_mask * mnar_missing_mask
+    indicating_mask = ori_mask - missing_mask
+    X[missing_mask == 0] = nan
+    return X_intact, X, missing_mask, indicating_mask
+
+
+def mnar_t(
+    X: Optional[Union[np.ndarray, torch.Tensor]],
+    cycle: float = 20,
+    pos: float = 10,
+    scale: float = 3,
+    nan: Union[float, int] = 0,
+) -> Union[
+    Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray],
+    Tuple[torch.Tensor, torch.Tensor, torch.Tensor, torch.Tensor],
+]:
+    """Create not-random missing values related to temporal dynamics (MNAR-t case).
+    In particular, the missingness is generated by an intensity function f(t) = exp(3*torch.sin(cycle*t + pos)).
+    This case mainly follows the setting in https://hawkeslib.readthedocs.io/en/latest/tutorial.html.
+
+    Parameters
+    ----------
+    X :
         Data vector. If X has any missing values, they should be numpy.nan.
 
-    rate : float, in (0,1),
-        Artificially missing rate, rate of the observed values which will be artificially masked as missing.
+    cycle :
+        The cycle of the used intensity function
+
+    pos :
+        The displacement of the used intensity function
 
-        Note that,
-        `rate` = (number of artificially missing values) / np.sum(~np.isnan(self.data)),
-        not (number of artificially missing values) / np.product(self.data.shape),
-        considering that the given data may already contain missing values,
-        the latter way may be confusing because if the original missing rate >= `rate`,
-        the function will do nothing, i.e. it won't play the role it has to be.
+    scale :
+        The scale number to control the missing rate
 
-    nan : int/float, optional, default=0
+    nan :
         Value used to fill NaN values.
 
     Returns
     -------
+    X_intact : array,
+        Original data with missing values (nan) filled with given parameter `nan`, with observed values intact.
+        X_intact is for loss calculation in the masked imputation task.
+
+    X : array,
+        Original X with artificial missing values. X is for model input.
+        Both originally-missing and artificially-missing values are filled with given parameter `nan`.
+
+    missing_mask : array,
+        The mask indicates all missing values in X.
+        In it, 1 indicates observed values, and 0 indicates missing values.
+
+    indicating_mask : array,
+        The mask indicates the artificially-missing values in X, namely missing parts different from X_intact.
+        In it, 1 indicates artificially missing values,
+        and the other values (including originally observed/missing values) are indicated as 0.
 
     """
-    # TODO: Create missing values in MNAR case
-    raise NotImplementedError("MNAR case has not been implemented yet.")
+    if isinstance(X, list):
+        X = np.asarray(X)
 
+    if isinstance(X, np.ndarray):
+        X_intact, X, missing_mask, indicating_mask = _mnar_t_numpy(
+            X, cycle, pos, scale, nan
+        )
+    elif isinstance(X, torch.Tensor):
+        X_intact, X, missing_mask, indicating_mask = _mnar_t_torch(
+            X, cycle, pos, scale, nan
+        )
+    else:
+        raise TypeError(
+            "X must be type of list/numpy.ndarray/torch.Tensor, " f"but got {type(X)}"
+        )
+    return X_intact, X, missing_mask, indicating_mask
 
-def _mnar_numpy(X, rate, nan=0):
-    pass
 
+def _mnar_t_numpy(
+    X: np.ndarray,
+    cycle: float = 20,
+    pos: float = 10,
+    scale: float = 3,
+    nan: Union[float, int] = 0,
+) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
+    # clone X to ensure values of X out of this function not being affected
+    X = np.copy(X)
 
-def _mnar_torch(X, rate, nan=0):
-    pass
+    X_intact = np.copy(X)  # keep a copy of originally observed values in X_intact
+
+    n_s, n_l, n_c = X.shape
+
+    ori_mask = (~np.isnan(X)).astype(np.float32)
+    X = np.nan_to_num(X, nan=nan)
+    X_intact = np.nan_to_num(X_intact, nan=nan)
+
+    ts = np.linspace(0, 1, n_l).reshape(1, n_l, 1)
+    ts = np.repeat(ts, n_s, axis=0)
+    ts = np.repeat(ts, n_c, axis=2)
+    intensity = np.exp(3 * np.sin(cycle * ts + pos))
+    mnar_missing_mask = (np.random.rand(n_s, n_l, n_c) * scale) < intensity
+
+    missing_mask = ori_mask * mnar_missing_mask
+    indicating_mask = ori_mask - missing_mask
+    X[missing_mask == 0] = nan
+    return X_intact, X, missing_mask, indicating_mask
+
+
+def _mnar_t_torch(
+    X: torch.Tensor,
+    cycle: float = 20,
+    pos: float = 10,
+    scale: float = 3,
+    nan: Union[float, int] = 0,
+) -> Tuple[torch.Tensor, torch.Tensor, torch.Tensor, torch.Tensor]:
+    # clone X to ensure values of X out of this function not being affected
+    X = torch.clone(X)
+
+    X_intact = torch.clone(X)  # keep a copy of originally observed values in X_intact
+    n_s, n_l, n_c = X.shape
+
+    ori_mask = (~torch.isnan(X)).type(torch.float32)
+    X = torch.nan_to_num(X, nan=nan)
+    X_intact = torch.nan_to_num(X_intact, nan=nan)
+
+    ts = torch.linspace(0, 1, n_l).reshape(1, n_l, 1).repeat(n_s, 1, n_c)
+    intensity = torch.exp(3 * torch.sin(cycle * ts + pos))
+    mnar_missing_mask = (torch.randn(X.size()).uniform_(0, 1) * scale) < intensity
+
+    missing_mask = ori_mask * mnar_missing_mask
+    indicating_mask = ori_mask - missing_mask
+    X[missing_mask == 0] = nan
+    return X_intact, X, missing_mask, indicating_mask


Test patch:
diff --git a/tests/test_pygrinder.py b/tests/test_pygrinder.py
index ee9aaa9..bd69795 100644
--- a/tests/test_pygrinder.py
+++ b/tests/test_pygrinder.py
@@ -9,7 +9,12 @@
 
 import numpy as np
 
-from pygrinder import mcar
+from pygrinder import (
+    mcar,
+    mar_logistic,
+    mnar_x,
+    mnar_t,
+)
 
 DEFAULT_MISSING_RATE = 0.1
 NAN = 1
@@ -17,15 +22,49 @@
 
 class TestPyGrinder(unittest.TestCase):
     def test_0_mcar(self):
-        d = np.random.randn(128, 10, 36)
-        d_intact, d_with_missing, missing_mask, indicating_mask = mcar(
-            d, p=DEFAULT_MISSING_RATE, nan=NAN
+        X = np.random.randn(128, 10, 36)
+        X_intact, X_with_missing, missing_mask, indicating_mask = mcar(
+            X, p=DEFAULT_MISSING_RATE, nan=NAN
         )
-        shape_product = np.product(d.shape)
+        shape_product = np.product(X_intact.shape)
         actual_missing_rate = (shape_product - missing_mask.sum()) / shape_product
         assert (
             round(actual_missing_rate, 1) == DEFAULT_MISSING_RATE
         ), f"Actual missing rate is {actual_missing_rate}, not given {DEFAULT_MISSING_RATE}"
-        assert np.sum(d_with_missing[(1 - missing_mask).astype(bool)]) == NAN * np.sum(
+        assert np.sum(X_with_missing[(1 - missing_mask).astype(bool)]) == NAN * np.sum(
             1 - missing_mask
         )
+
+    def test_1_mar(self):
+        X = np.random.randn(128, 36)
+        X_intact, X_with_missing, missing_mask, indicating_mask = mar_logistic(
+            X, obs_rate=0.1, missing_rate=0.2
+        )
+        shape_product = np.product(X_intact.shape)
+        actual_missing_rate = (shape_product - missing_mask.sum()) / shape_product
+        assert (
+            round(actual_missing_rate, 1) > 0
+        ), f"Actual missing rate is {actual_missing_rate}"
+
+    def test_2_mnar(self):
+        X = np.random.randn(128, 10, 36)
+
+        # mnar_x
+        X_intact, X_with_missing, missing_mask, indicating_mask = mnar_x(
+            X, offset=0, nan=NAN
+        )
+        shape_product = np.product(X_intact.shape)
+        actual_missing_rate = (shape_product - missing_mask.sum()) / shape_product
+        assert (
+            round(actual_missing_rate, 1) > 0
+        ), f"Actual missing rate is {actual_missing_rate}"
+
+        # mnar_t
+        X_intact, X_with_missing, missing_mask, indicating_mask = mnar_t(
+            X, cycle=20, pos=10, scale=3, nan=NAN
+        )
+        shape_product = np.product(X_intact.shape)
+        actual_missing_rate = (shape_product - missing_mask.sum()) / shape_product
+        assert (
+            round(actual_missing_rate, 1) > 0
+        ), f"Actual missing rate is {actual_missing_rate}"