{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"tpu1vmV38","dataSources":[{"sourceId":82611,"databundleVersionId":9553358,"sourceType":"competition"},{"sourceId":12784398,"sourceType":"datasetVersion","datasetId":8082575},{"sourceId":12784401,"sourceType":"datasetVersion","datasetId":8082577}],"dockerImageVersionId":31091,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"from tqdm import tqdm\nfrom itertools import product\n\nimport numpy as np\nimport pandas as pd\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom sklearn.model_selection import GroupKFold\nfrom sklearn.impute import SimpleImputer, KNNImputer\nimport warnings\n\n\nfrom sklearn.ensemble import StackingRegressor\nwarnings.filterwarnings(\"ignore\")\n\nfrom itertools import product\nfrom scipy.stats import norm\nfrom scipy.special import logsumexp\n\n\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# df_train = pd.read_csv(\n#     '/kaggle/input/brist1d/train.csv',\n#     index_col='id',\n#     parse_dates=['time'],\n# )\n\n# df_train.columns = df_train.columns.str.replace(':', '-')","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# import random\n# import os\n# import numpy as np\n# import torch\n\n\n# def seed_everything(seed):\n#     random.seed(seed)\n#     os.environ[\"PYTHONHASHSEED\"] = str(seed)\n#     np.random.seed(seed)\n#     torch.manual_seed(seed)\n#     torch.cuda.manual_seed(seed)\n#     torch.backends.cudnn.deterministic = True\n#     torch.backends.cudnn.benchmark = True\n\n\n# df_train.columns = df_train.columns.str.replace(\":\", \"-\")\n\n\n# for colset in [bg_cols, insu_cols, carb_cols, hr_cols, step_cols, cals_cols]:\n#     df_train[colset] = (\n#         df_train[colset]\n#         .interpolate(axis=1)\n#         .fillna(method=\"bfill\", axis=1)\n#         .fillna(method=\"ffill\", axis=1)\n#     )\n\n# feature_columns = {\n#     'bg': bg_cols,\n#     'hr': hr_cols,\n#     'steps': step_cols\n# }\n\n# def impute_time_series_features(df, columns, method='median', window=5):\n#     for col in columns:\n#         df[col] = df[col].fillna(df[col].rolling(window=window, min_periods=1).median())\n\n#         df[col] = df[col].fillna(method='ffill')\n#         df[col] = df[col].fillna(method='bfill')\n\n#         df[col] = df.groupby('p_num')[col].transform(lambda x: x.fillna(x.median()))\n#     return df\n\n# for feature, columns in feature_columns.items():\n#     df_train = impute_time_series_features(df_train, columns, method='median')\n\n# all_imputed_columns = bg_cols + hr_cols + step_cols\n# print(\"Remaining null values in training set:\", df_train[all_imputed_columns].isnull().sum().sum())","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\n# feature_columns_insulin_cals = {\n#     'insulin': insu_cols,\n#     'cals': cals_cols\n# }\n\n# def impute_non_temporal_features(df, columns, method='median', window=5):\n#     for col in columns:\n#         df[col] = df.groupby('p_num')[col].transform(lambda x: x.fillna(x.median()))\n\n#         df[col] = df[col].fillna(df[col].rolling(window=window, min_periods=1).median())\n#     return df\n\n# for feature, columns in feature_columns_insulin_cals.items():\n#     df_train = impute_non_temporal_features(df_train, columns, method='median')\n\n# all_imputed_columns_insulin_cals = insu_cols + cals_cols\n# print(\"Remaining null values in training set:\", df_train[all_imputed_columns_insulin_cals].isnull().sum().sum())\n\n\n\n# feature_columns_insulin_cals = {\n#     'insulin': insu_cols,\n#     'cals': cals_cols\n# }\n\n# def impute_non_temporal_features(df, columns, method='median', window=5):\n#     for col in columns:\n#         df[col] = df.groupby('p_num')[col].transform(lambda x: x.fillna(x.median()))\n\n#         df[col] = df[col].fillna(df[col].rolling(window=window, min_periods=1).median())\n#     return df\n\n# for feature, columns in feature_columns_insulin_cals.items():\n#     df_train = impute_non_temporal_features(df_train, columns, method='median')\n\n# all_imputed_columns_insulin_cals = insu_cols + cals_cols\n# print(\"Remaining null values in training set:\", df_train[all_imputed_columns_insulin_cals].isnull().sum().sum())\n","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# activity_cols = [col for col in df_train.columns if col.startswith('activity')]\n# dominant_activity_train = df_train[activity_cols].mode(axis=1)[0]\n# df_train['dominant_activity'] = dominant_activity_train\n# print(df_train['dominant_activity'].value_counts())","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# activity_replace = {\n#     'Walk': 'Walking',\n#     'Run': 'Running',\n#     'Swim': 'Swimming',\n# }\n# df_train['dominant_activity'] = df_train['dominant_activity'].replace(activity_replace)\n# print(df_train['dominant_activity'])","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# df_train['dominant_activity'].fillna('No Activity', inplace=True)\n\n# from sklearn.preprocessing import LabelEncoder\n\n# label_encoder = LabelEncoder()\n\n# combined_activities = pd.concat([\n#     df_train['dominant_activity']\n# ])\n\n# label_encoder.fit(combined_activities)\n\n# df_train['dominant_activity_encoded'] = label_encoder.transform(df_train['dominant_activity'])\n\n# df_train = df_train.drop(['dominant_activity'], axis = 1)\n\n# feature_cols = feature_cols + ['dominant_activity_encoded']\n\n\n# imputer = SimpleImputer()\n\n# df_train[feature_cols] = imputer.fit_transform(df_train[feature_cols])\n\n# df_train['time_hour'] = pd.to_datetime(df_train['time']).dt.hour\n\n# feature_cols.extend([\"time_hour\"])\n\n# df_train_final = df_train[feature_cols]\n\n# y_target = df_train[[target_col]]\n\n\n# np.isnan(df_train_final).sum()\n","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# y_clean.shape","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# X_clean.shape","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# X_clean =pd.read_csv(\"/kaggle/input/x-clean/X_clean.csv\")\n# y_clean=pd.read_csv(\"/kaggle/input/y-clean/y_clean.csv\")\n\n","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# X=pd.concat([X_clean,y_clean],axis=1)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# X.shape","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# import pandas as pd\n# import matplotlib.pyplot as plt\n# import seaborn as sns\n\n\n# # Nếu có quá nhiều đặc trưng, hãy vẽ cho từng nhóm nhỏ\n# for i in range(0, 434, 20):\n#     plt.figure(figsize=(15, 6))\n#     sns.boxplot(data=X.iloc[:, i:i+20])\n#     plt.title(f\"Box Plot cho đặc trưng {i} đến {i+19}\")\n#     plt.show()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# import pandas as pd\n# from collections import Counter\n\n# outlier_indices = []\n\n# for column in X.columns[-1]:\n#     Q1 = X[column].quantile(0.25)\n#     Q3 = X[column].quantile(0.75)\n#     IQR_val = Q3 - Q1\n#     lower_bound = Q1 - 1.5 * IQR_val\n#     upper_bound = Q3 + 1.5 * IQR_val\n#     column_outliers = X[(X[column] < lower_bound) | (X[column] > upper_bound)].index\n#     outlier_indices.extend(column_outliers)\n\n# outlier_counts = Counter(outlier_indices)\n\n# feature_threshold = 1\n\n# multiple_outliers = [i for i, count in outlier_counts.items() if count > feature_threshold]\n\n# print(f\"Số hàng ban đầu: {len(X)}\")\n# print(f\"Số hàng bị coi là outlier nghiêm trọng: {len(multiple_outliers)}\")\n\n# X_cleaned = X.drop(multiple_outliers)\n# if 'y' in locals():\n#     y_cleaned = y.drop(multiple_outliers)\n\n# print(f\"Số hàng sau khi xóa: {len(X_cleaned)}\")\n","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# import pandas as pd\n# import matplotlib.pyplot as plt\n# import seaborn as sns\n\n\n# # Nếu có quá nhiều đặc trưng, hãy vẽ cho từng nhóm nhỏ\n# for i in range(0, 434, 20):\n#     plt.figure(figsize=(15, 6))\n#     sns.boxplot(data=X.iloc[:, i:i+20])\n#     plt.title(f\"Box Plot cho đặc trưng {i} đến {i+19}\")\n#     plt.show()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"!pip install skglm\n!pip install statsmodels\n# !pip install celer\n!pip install -U scikit-learn-intelex skglm\n","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"#1.SI\n","metadata":{}},{"cell_type":"code","source":"# def _log_cdf_diff(V_r, V_l):\n#     \"\"\"\n#     Tính log(P(V_l < Z < V_r)) một cách ổn định về mặt số học.\n#     Z ~ N(0, 1).\n#     Sử dụng công thức: log(a - b) = log(a) + log(1 - exp(log(b) - log(a)))\n#     \"\"\"\n#     if V_l >= V_r:\n#         return -np.inf  # log(0)\n\n#     log_cdf_r = norm.logcdf(V_r)\n#     log_cdf_l = norm.logcdf(V_l)\n\n#     # np.log1p(x) tính log(1 + x) một cách chính xác hơn cho x nhỏ.\n#     # Ta có log(1 - exp(y)) = log1p(-exp(y))\n#     # y = log_cdf_l - log_cdf_r luôn âm vì V_l < V_r\n#     return log_cdf_r + np.log1p(-np.exp(log_cdf_l - log_cdf_r))\n\n# def _calculate_pivot_under_null(V_intervals, V_obs):\n#     \"\"\"\n#     Tính giá trị của đại lượng chốt (CDF của phân phối cụt) dưới H0.\n#     Giá trị trả về tuân theo phân phối Uniform(0, 1).\n#     \"\"\"\n#     # 1. Tính log(denominator) - Tổng khối xác suất của vùng cắt cụt\n#     log_probs_denom = [_log_cdf_diff(V_r, V_l) for V_l, V_r in V_intervals]\n#     log_denominator = logsumexp(log_probs_denom)\n#     print(\"log_denominator:\",log_denominator)\n    \n\n\n#     # 2. Tính log(numerator) - Khối xác suất từ -inf đến z_obs trong vùng cắt cụt\n#     log_probs_num = []\n#     for V_l, V_r in V_intervals:\n#         # Khoảng tích phân là [V_l, min(V_r, V_obs)]\n#         integration_end = min(V_r, V_obs)\n#         if V_l < integration_end:\n#             log_probs_num.append(_log_cdf_diff(integration_end, V_l))\n    \n    \n#     log_numerator = logsumexp(log_probs_num)\n#     print(\"log_numerator:\",log_numerator)\n\n#     # 3. Tính log của pivot\n#     log_pivot = log_numerator - log_denominator\n    \n#     return np.exp(log_pivot)\n\n# def p_value(A_obs, list_active_set, list_zk, etaj, z_obs, cov_matrix):\n#     \"\"\"\n#     Tính p-value cho suy luận chọn lọc bằng cách trước tiên tính đại lượng chốt (pivot).\n#     \"\"\"\n#     print(\"    Bắt đầu tính p-value thông qua pivot trên log-space...\")\n#     start_time_pval = time.time()\n    \n#     # --- Bước A: Xây dựng vùng cắt cụt z_intervals ---\n#     z_intervals = []\n#     set_A_obs = set(A_obs)\n#     for i in range(len(list_zk) - 1):\n#         if set(list_active_set[i]) == set_A_obs:\n#             z_intervals.append((list_zk[i], list_zk[i+1]))\n\n#     if not z_intervals:\n#         print(\"    CẢNH BÁO: Không tìm thấy khoảng khớp. Trả về p-value = 1.0\")\n#         return 1.0\n\n#     # --- Bước B: Tính tham số phân phối ---\n#     sq_norm_eta = (etaj.T @ etaj).item()\n#     sigma_squared_hat = cov_matrix[0, 0]\n#     tn_sigma = np.sqrt(sq_norm_eta * sigma_squared_hat)\n    \n#     if tn_sigma < 1e-12: tn_sigma = 1e-12\n#     print(f\"    tn_sigma tính được: {tn_sigma:.6f}\")\n\n#     # --- Bước C: Chuẩn hóa sang không gian V ---\n#     V_intervals = [(z_l / tn_sigma, z_r / tn_sigma) for z_l, z_r in z_intervals]\n#     V_obs = z_obs / tn_sigma\n\n#     # --- Bước D: Tính giá trị pivot dưới H0 ---\n#     pivot_value = _calculate_pivot_under_null(V_intervals, V_obs)\n#     print(f\"    Giá trị Pivot (CDF) tính được: {pivot_value:.6f}\")\n    \n#     # --- Bước E: Chuyển pivot thành p-value hai phía ---\n#     final_p_value = 2 * min(pivot_value, 1 - pivot_value)\n    \n#     print(f\"    Thời gian tính p-value: {(time.time() - start_time_pval):.4f}s\")\n\n#     return final_p_value","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"!pip install skglm\nfrom tqdm import tqdm\nfrom itertools import product\n\nimport numpy as np\nimport pandas as pd\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom sklearn.model_selection import GroupKFold\nfrom sklearn.impute import SimpleImputer, KNNImputer\nimport warnings\n\n\nfrom sklearn.ensemble import StackingRegressor\nwarnings.filterwarnings(\"ignore\")\n\nfrom itertools import product\nfrom scipy.stats import norm\nfrom scipy.special import logsumexp\n\n\n\nimport numpy as np\nimport logging\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport statsmodels.api as sm\nfrom sklearn.preprocessing import StandardScaler\n# from celer import Lasso\nfrom skglm import Lasso\nfrom sklearn import linear_model\nfrom statsmodels.stats.multitest import multipletests\nfrom mpmath import mp\nimport time\n\nfrom joblib import Parallel, delayed\nfrom multiprocessing import shared_memory\n\nimport os\n\nmp.dps =400\n\ndef construct_A_XA_Ac_XAc_bhA(X, bh):\n    p = X.shape[1]\n    A = sorted(list(np.where(bh != 0)[0]))\n    XA = X[:, A] if A else np.empty((X.shape[0], 0))\n\n    return A, XA\n\ndef construct_test_statistic(j, XA, y, A):\n    ej = np.zeros((len(A), 1))\n    ej[A.index(j)] = 1\n    inv = np.linalg.pinv(XA.T @ XA)\n    etaj = XA @ inv @ ej\n    etajTy = (etaj.T @ y).item()\n\n    return etaj, etajTy\n\n\n\ndef compute_yz(y, etaj, zk, n, sigma_squared):\n\n    sq_norm_val =(etaj.T @ etaj).item()\n    sq_norm = sigma_squared * sq_norm_val\n    \n    \n    if sq_norm < 1e-12: sq_norm = 1e-12\n\n    # --- PHẦN TỐI ƯU HÓA ---\n    # Thay vì tạo ma trận proj_matrix (n x n), chúng ta tính toán trực tiếp.\n    # a = (I - (ηηᵀ)/||η||²)y = y - η(ηᵀy)/||η||²\n    \n    # 1. Tính tích trong (inner product) trước -> ra một số vô hướng\n    # etajTy_scalar = (etaj.T @ y)[0, 0]\n    etajTy_scalar = (etaj.T @ y).item()\n\n    \n    # 2. Tính vector \"a\" bằng các phép toán trên vector\n    a = y - etaj * (etajTy_scalar / sq_norm)*sigma_squared\n    # -----------------------\n    \n    b = (etaj / sq_norm) * sigma_squared\n    yz = a + b * zk\n   \n    return yz, b\n\n\n\n\ndef compute_polyhedron(X, M, s, lambda_val):\n    n, p = X.shape\n    if not M: return np.empty((0, n)), np.empty(0)\n    \n    X_M = X[:, M]\n    X_notM_indices = np.setdiff1d(np.arange(p), M)\n    X_notM = X[:, X_notM_indices]\n    \n    try:\n        XATXA_inv = np.linalg.pinv(X_M.T @ X_M)\n    except np.linalg.LinAlgError:\n        return np.empty((0, n)), np.empty(0)\n\n    A0, b0 = np.empty((0, n)), np.empty(0)\n    if X_notM.shape[1] > 0:\n        # TỐI ƯU HÓA TỐT HƠN NỮA\n        # Tính X_notM.T @ P_M mà không cần tạo P_M (n x n)\n        # X_notM.T @ P_M = X_notM.T @ (X_M @ inv(X_M.T @ X_M) @ X_M.T)\n        #                = (X_notM.T @ X_M) @ inv(X_M.T @ X_M) @ X_M.T\n        term_k_x_k = X_notM.T @ X_M @ XATXA_inv\n        X_notM_T_P_M = term_k_x_k @ X_M.T\n        \n        # A0_T = (1/lambda_val) * (X_notM.T - X_notM.T @ P_M)\n        A0_T = (1/lambda_val) * (X_notM.T - X_notM_T_P_M)\n        \n        A0 = np.vstack([A0_T, -A0_T])\n        \n        XT_M_pinv = np.linalg.pinv(X_M.T)\n        s_term = X_notM.T @ XT_M_pinv @ s\n        s_term = s_term.flatten()\n        b0 = np.hstack([np.ones(X_notM.shape[1]) - s_term, \n                        np.ones(X_notM.shape[1]) + s_term])\n    \n    A1 = -np.diag(s) @ XATXA_inv @ X_M.T\n    b1 = -lambda_val * np.diag(s) @ XATXA_inv @ s\n    \n    return np.vstack([A0, A1]), np.hstack([b0, b1])\n\n\n\n\n\ndef compute_truncation_bounds(A, b, eta, y):\n    eta = eta.ravel()\n    y = y.ravel()\n    eta_norm = np.dot(eta, eta)\n    if eta_norm < 1e-10:\n        raise ValueError(\"eta has near-zero norm\")\n    c = eta / eta_norm\n    z = y - np.dot(eta, y) * c\n    Ac = A @ c\n    Az = A @ z\n    V_minus = None\n    V_plus = None\n    V_zero = np.inf\n    for j in range(A.shape[0]):\n        if abs(Ac[j]) == 0:\n            V_zero = min(V_zero, b[j] - Az[j])\n        elif Ac[j] < 0:\n            if V_minus is None or (b[j] - Az[j]) / Ac[j] > V_minus:\n                V_minus = (b[j] - Az[j]) / Ac[j]\n        elif Ac[j] > 0:\n            if V_plus is None or (b[j] - Az[j]) / Ac[j] < V_plus:\n                V_plus = (b[j] - Az[j]) / Ac[j]\n    if V_minus is None:\n        V_minus = -np.inf\n    if V_plus is None:\n        V_plus = np.inf\n    if V_zero < 0:\n        raise ValueError(\"Invalid truncation bounds: V_zero < 0\")\n\n \n\n    return V_minus, V_plus\n\ndef create_shared_array(arr: np.ndarray):\n    shm = shared_memory.SharedMemory(create=True, size=arr.nbytes)\n    shm_arr = np.ndarray(arr.shape, dtype=arr.dtype, buffer=shm.buf)\n    shm_arr[:] = arr[:]\n    return shm, arr.shape, arr.dtype.str\n\ndef attach_shared_array(shm_name, shape, dtype_str):\n    shm = shared_memory.SharedMemory(name=shm_name)\n    arr = np.ndarray(shape, dtype=np.dtype(dtype_str), buffer=shm.buf)\n    return shm, arr\n\n# ============================\n# Worker xử lý 1 đoạn zk\n# ============================\n\n\ndef setup_worker_logging(worker_id):\n    \"\"\"Cấu hình logger cho mỗi tiến trình con.\"\"\"\n    logger = logging.getLogger(f\"worker_{worker_id}\")\n    logger.setLevel(logging.INFO)  # Chỉ ghi log từ mức INFO trở lên\n\n    # Tạo handler để ghi ra console\n    handler = logging.StreamHandler()\n    \n    # Tạo formatter để định dạng log message\n    # Bao gồm: Thời gian - ID Tiến trình - Tên Logger - Mức độ - Nội dung\n    formatter = logging.Formatter(\n        '%(asctime)s - PID:%(process)d - %(name)s - %(levelname)s - %(message)s'\n    )\n    handler.setFormatter(formatter)\n\n    # Thêm handler vào logger (nếu chưa có)\n    if not logger.handlers:\n        logger.addHandler(handler)\n        \n    return logger\n\n\n\n\n\ndef process_segment(z_start, z_end,\n                    X_shm_name, X_shape, X_dtype,\n                    y_shm_name, y_shape, y_dtype,\n                    etaj_shm_name, etaj_shape, etaj_dtype,\n                    alpha, sigma_sq):\n\n\n    small_eps=0.000001\n    \n\n    worker_id = os.getpid()\n    logger = setup_worker_logging(worker_id)\n    logger.info(f\"BẮT ĐẦU xử lý đoạn [{float(z_start):.4f}, {float(z_end):.4f}]\")\n\n    \n    start_time = time.time()\n\n    # Attach shared arrays\n    shmX, X = attach_shared_array(X_shm_name, X_shape, X_dtype)\n    shmy, y = attach_shared_array(y_shm_name, y_shape, y_dtype)\n    shmetaj, etaj = attach_shared_array(etaj_shm_name, etaj_shape, etaj_dtype)\n\n    n = X.shape[0]\n    results = []\n    clf = Lasso(alpha=alpha / n, fit_intercept=False, warm_start=True,max_iter=80000,tol=1e-9)\n\n    zk = float(z_start)\n    it = 0\n   \n    log_interval = 250\n\n    while zk < z_end:\n        it += 1\n        if it == 1 or it % log_interval == 0:\n            logger.info(f\"Vòng lặp #{it}... zk hiện tại = {float(zk):.6f}\")\n\n\n        yz, b_vec = compute_yz(y, etaj, zk, n, sigma_sq)\n        yz_flat = yz.flatten()\n\n    \n        \n       \n        clf.fit(X, yz_flat)\n       \n        bhz = clf.coef_.copy()\n        # Bước \"làm sạch\"\n        coeff_threshold = 1e-9 \n        bhz[np.abs(bhz) < coeff_threshold] = 0\n\n        A, _ = construct_A_XA_Ac_XAc_bhA(X, bhz)\n\n        if len(A) == 0:\n            V_l, V_r = zk, z_end\n        else:\n            s = np.sign(bhz[A])\n            A1, b1 = compute_polyhedron(X, A, s, alpha)\n        \n            V_l, V_r = compute_truncation_bounds(A1, b1, etaj, yz_flat)\n            \n        print(f\"Vòng lặp #{it}: V_r={V_r}, V_l={V_l}, V_r-V_l={V_r-V_l}\") # In ra để debug\n        results.append((zk, A, bhz))\n        \n        # KIỂM TRA TÍNH HỢP LỆ CỦA KHOẢNG TRƯỚC KHI CẬP NHẬT\n        if np.isfinite(V_r) and np.isfinite(V_l) and V_r > V_l and V_r > zk + 1e-12:\n           \n            zk = min(V_r + small_eps, z_end)\n            \n        else:\n\n            logger.warning(f\"Vòng lặp #{it}: Đa diện không hợp lệ hoặc không tiến lên (V_r={V_r}, V_l={V_l}). Dùng bước nhảy nhỏ.\")\n            zk = zk + small_eps\n            \n   \n    shmX.close()\n    shmy.close()\n    shmetaj.close()\n    end_time = time.time()\n    duration = end_time - start_time\n    \n   \n\n    logger.info(f\"KẾT THÚC xử lý đoạn [{z_start.item():.4f}, {z_end.item():.4f}]. \"\n            f\"Tổng cộng: {it} vòng lặp. Thời gian: {duration:.2f} giây.\")\n\n    return results\n\n\n\n\n# ============================\n# Hàm song song hóa chính\n# ============================\n\ndef run_parametric_lasso(X, y, alpha, etaj, threshold, sigma_sq,\n                                 n_segments, overlap=0, **worker_kwargs):\n    shmX, X_shape, X_dtype = create_shared_array(X)\n    shmy, y_shape, y_dtype = create_shared_array(y.reshape(-1))\n    shmetaj, etaj_shape, etaj_dtype = create_shared_array(etaj.reshape(-1))\n\n    edges = np.linspace(-threshold, threshold, n_segments + 1)\n    segments = []\n    overlap=0\n    for i in range(n_segments):\n        z0 = edges[i] - (overlap if i > 0 else 0.0)\n        z1 = edges[i+1] + (overlap if i < n_segments - 1 else 0.0)\n        z0 = max(z0, -threshold)\n        z1 = min(z1, threshold)\n        segments.append((z0, z1))\n\n    parallel = Parallel(n_jobs=n_segments)\n    jobs = (delayed(process_segment)(\n                z0, z1,\n                shmX.name, X_shape, X_dtype,\n                shmy.name, y_shape, y_dtype,\n                shmetaj.name, etaj_shape, etaj_dtype,\n                alpha, sigma_sq\n            ) for (z0, z1) in segments)\n    start=time.time()\n    segment_results = parallel(jobs)\n    print(\"Thời gian cho ra segment_results: \",time.time()-start)\n\n    shmX.close(); shmX.unlink()\n    shmy.close(); shmy.unlink()\n    shmetaj.close(); shmetaj.unlink()\n\n    all_entries = []\n    for seg in segment_results:\n        all_entries.extend(seg)\n    all_entries.sort(key=lambda t: t[0])\n\n    merged_zk = []\n    merged_active_sets = []\n    merged_bhz = []\n    tol = 0.00001\n    for zk, A, bhz in all_entries:\n        if merged_zk and abs(zk - merged_zk[-1]) < tol:\n            continue\n        if merged_active_sets and (merged_active_sets[-1]==A):\n            continue\n        merged_zk.append(float(zk))\n        merged_active_sets.append(A)\n        merged_bhz.append(bhz)\n\n    if len(merged_zk) == 0 or merged_zk[-1] < threshold - 1e-12:\n        merged_zk.append(float(threshold))\n\n    return merged_zk, merged_active_sets, merged_bhz\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n# def p_value(A_obs, list_active_set, list_zk, etaj, z_obs, cov_matrix):\n#     start=time.time()\n#     tn_sigma = (np.sqrt((etaj.T @ cov_matrix @ etaj))).item()\n#     print(\"tn_sigma:\",tn_sigma)\n#     if tn_sigma < 1e-9: return None\n    \n#     z_intervals = []\n#     for i, active_set in enumerate(list_active_set):\n#         if set(A_obs) == set(active_set):\n#             z_intervals.append([list_zk[i], list_zk[i+1]])\n\n#     print(\"z_intervals: \",z_intervals)\n\n#     if not z_intervals: return None\n\n#     tn_mu = 0\n#     numerator = mp.mpf(0)\n#     denominator = mp.mpf(0)\n    \n#     for al, ar in z_intervals:\n#         prob_interval = mp.ncdf((ar - tn_mu)/tn_sigma) - mp.ncdf((al - tn_mu)/tn_sigma)\n#         denominator += prob_interval\n#         if z_obs >= ar:\n#             numerator += prob_interval\n#         elif z_obs > al:\n#             numerator += mp.ncdf((z_obs - tn_mu)/tn_sigma) - mp.ncdf((al - tn_mu)/tn_sigma)\n#     print(\"numerator: \",numerator)       \n#     print(\"denominator:\",denominator)\n   \n#     # if denominator<1e-100:\n    \n#     #     return 0\n    \n#     pivot_val = float(numerator / denominator)\n#     print(\"Thơi gian tính p value: \",time.time()-start)\n#     return 2 * min(pivot_val, 1 - pivot_val)\n\nimport mpmath as mp\nimport time\n\nmp.mp.dps = 2000  # độ chính xác (có thể chỉnh)\n\n# ====== Các hàm hỗ trợ ======\n\ndef log_survival(t):\n    return mp.log(0.5 * mp.erfc(t / mp.sqrt(2)))\n\ndef log_interval_mass(la, lb):\n    delta = lb - la\n    if delta <= -50:\n        return la\n    return la + mp.log1p(- mp.e**(delta))\n\ndef log_sum_exp(log_vals):\n    if not log_vals:\n        return mp.ninf\n    L = max(log_vals)\n    if L == mp.ninf:\n        return mp.ninf\n    s = mp.mpf('0')\n    for lv in log_vals:\n        s += mp.e**(lv - L)\n    return L + mp.log(s)\n\n# ====== Hàm p_value ổn định ======\n\ndef p_value(A_obs, list_active_set, list_zk, etaj, z_obs, cov_matrix):\n    start = time.time()\n    \n    # 1) Tính sigma\n    tn_sigma = (np.sqrt((etaj.T @ cov_matrix @ etaj))).item()\n    print(\"tn_sigma:\", tn_sigma)\n    if tn_sigma < 1e-9:\n        return None\n\n    # 2) Lấy các khoảng cắt cụt\n    z_intervals = []\n    for i, active_set in enumerate(list_active_set):\n        if set(A_obs) == set(active_set):\n            z_intervals.append((list_zk[i], list_zk[i+1]))\n\n    print(\"z_intervals:\", z_intervals)\n    if not z_intervals:\n        return None\n\n    tn_mu = 0\n\n    # 3) Chuẩn hóa z và các khoảng\n    z_std = (z_obs - tn_mu) / tn_sigma\n    list_intervals_std = [((al - tn_mu)/tn_sigma, (ar - tn_mu)/tn_sigma) for al, ar in z_intervals]\n\n    # 4) Tính logN và logD\n    log_masses = []\n    log_parts = []\n    lz = log_survival(z_std)\n\n    for (a_std, b_std) in list_intervals_std:\n        la = log_survival(a_std)\n        lb = log_survival(b_std)\n\n        # toàn bộ khối lượng khoảng\n        log_mk = log_interval_mass(la, lb)\n        log_masses.append(log_mk)\n\n        # phần cho numerator\n        if z_std <= a_std:\n            pass  # không đóng góp\n        elif z_std >= b_std:\n            log_parts.append(log_mk)\n        else:\n            log_pk = log_interval_mass(la, lz)\n            log_parts.append(log_pk)\n\n    logD = log_sum_exp(log_masses)\n    logN = log_sum_exp(log_parts)\n\n    if logD == mp.ninf:\n        return None  # denominator = 0\n\n    pivot = mp.e**(logN - logD)\n    p_value = 2 * mp.mpf(min(pivot, 1 - pivot))\n\n    print(\"pivot:\", pivot)\n    print(\"p-value:\", p_value)\n    print(\"Thời gian tính p value:\", time.time() - start)\n\n    return float(p_value)\n\n\n\n\n# =============================================================================\n# HÀM CHẠY CHÍNH: QUY TRÌNH PHÂN TÍCH DỮ LIỆU THỰC\n# =============================================================================\nimport numpy as np\nimport random\nfrom sklearn.preprocessing import RobustScaler\nif __name__ == '__main__':\n    # --- BƯỚC 0: CHUẨN BỊ DỮ LIỆU ĐẦU VÀO ---\n    \n    print(\"Sử dụng dữ liệu X_clean và y_clean...\")\n\n            \n    X_clean =pd.read_csv(\"/kaggle/input/x-clean/X_clean.csv\")\n    y_clean=pd.read_csv(\"/kaggle/input/y-clean/y_clean.csv\")\n \n\n    \n    \n    print(f\"Dữ liệu đầu vào: X_clean.shape={X_clean.shape}, y_clean.shape={y_clean.shape}\")\n\n    # --- GIAI ĐOẠN 1: LỰA CHỌN MÔ HÌNH BẰNG LASSO (KHÔNG DÙNG CV) ---\n    print(\"\\n--- Giai đoạn 1: Lựa chọn Mô hình bằng Lasso ---\")\n\n    # 1.1. Chuẩn hóa X\n    print(\"Đang chuẩn hóa X_clean...\")\n    scaler = StandardScaler()\n\n\n    X_columns = X_clean.columns\n    X_scaled = scaler.fit_transform(X_clean)\n   \n    y_scaler = StandardScaler(with_std=False) \n    y_clean = y_scaler.fit_transform(y_clean)\n    n, p = X_scaled.shape\n    BEST_ALPHA = 10000\n    print(f\"Sử dụng alpha cố định: {BEST_ALPHA:.4f}\")\n    lasso_model = Lasso(alpha=BEST_ALPHA/n, fit_intercept=False,max_iter=80000,tol=1e-9).fit(X_scaled, y_clean)\n    A_obs, XA_obs = construct_A_XA_Ac_XAc_bhA(X_scaled, lasso_model.coef_)\n\n    if not A_obs:\n        print(\"Lasso không chọn đặc trưng nào. Kết thúc.\")\n        exit()\n\n    print(f\"\\nLasso đã chọn {len(A_obs)} đặc trưng.\")\n    print(\"Các đặc trưng được chọn:\", X_columns[A_obs].tolist())\n\n\n    print(\"Ước lượng sigma^2 từ phần dư của Lasso với bậc tự do điều chỉnh (phương pháp Reid et al.)...\")\n    \n    # 1. Lấy y_pred trực tiếp từ mô hình Lasso đã fit\n    y_pred_lasso = lasso_model.predict(X_scaled)\n    \n    # 2. Tính phần dư (residuals)\n    # Đảm bảo cả hai đều là mảng 1D để trừ\n    if isinstance(y_clean, pd.Series) or isinstance(y_clean, pd.DataFrame):\n        y_true_flat = y_clean.values.flatten()\n    else:\n        y_true_flat = y_clean.flatten()\n        \n    residuals_lasso = y_true_flat - y_pred_lasso.flatten()\n\n\n\n    # 3. Tính tổng bình phương phần dư (RSS)\n    rss_lasso = np.sum(residuals_lasso**2)\n    \n    # 4. Xác định bậc tự do hiệu dụng (effective degrees of freedom)\n    # Theo lý thuyết, bậc tự do hiệu dụng của Lasso chính là số lượng biến khác 0\n    df_residuals = n - len(A_obs) \n    if df_residuals <= 0:\n        raise ValueError(f\"Không thể ước lượng sigma^2 vì bậc tự do (n - |A_obs|) = ({n} - {len(A_obs)}) <= 0.\")\n    \n    # 5. Tính ước lượng phương sai nhiễu\n    sigma_squared_hat = rss_lasso / df_residuals\n    print(f\"Ước lượng phương sai nhiễu (σ^2) : {sigma_squared_hat:.4f}\")\n\n    \n\n\n        \n\n    # --- GIAI ĐOẠN 2: SUY LUẬN CHỌN LỌC ---\n    print(f\"\\n--- Giai đoạn 2: Bắt đầu Suy luận Chọn lọc cho {len(A_obs)} đặc trưng ---\")\n    p_values_dict = {}\n    \n    y_col_vector = y_clean.reshape(-1, 1)\n\n    for i, j_selected in enumerate(A_obs):\n        if (i==0):\n        \n            feature_name = X_columns[j_selected]\n            start_time = time.time()\n            print(f\"  [{i+1}/{len(A_obs)}] Đang phân tích '{feature_name}'...\")\n            etaj, etajTy_orig = construct_test_statistic(j_selected, XA_obs, y_col_vector, A_obs)\n            sq_norm_eta = (etaj.T @ etaj)\n            if sq_norm_eta < 1e-12: sq_norm_eta = 1e-12\n            \n      \n            etajTy_scalar = (etaj.T @ y_col_vector)\n            a_vec = y_col_vector - etaj * (etajTy_scalar / sq_norm_eta) # Cách hiệu quả\n            z_obs = etajTy_orig \n            print(\"z_obs: \",z_obs)\n            sigma_z = np.sqrt(sq_norm_eta * sigma_squared_hat)\n            \n            threshold = 20*sigma_z\n            if (abs(z_obs)>threshold):\n                print(\"Vì z_obs>20*sigma_z nên threshold=z_obs*1.5 \")\n                threshold=abs(z_obs)*1.5\n            print(\"threshold: \",threshold)\n      \n            \n            n_segments = os.cpu_count()  # số core logic của CPU\n            \n            list_zk, list_active_set, list_bhz = run_parametric_lasso(\n                X_scaled, y_col_vector, BEST_ALPHA, etaj,\n                threshold=threshold, sigma_sq=sigma_squared_hat,\n                n_segments=n_segments,\n                small_eps=0.001, max_iter=10000\n            )\n        \n            cov_matrix = np.identity(n) * sigma_squared_hat\n        \n            p_val = p_value(A_obs, list_active_set, list_zk, etaj, z_obs, cov_matrix)\n            \n            \n            p_values_dict[feature_name] = p_val\n            elapsed_time = time.time() - start_time\n            print(f\"    -> p-value = {p_val if p_val is not None else 'N/A'} (Thời gian: {elapsed_time:.2f}s)\")\n            break\n        \n    # --- GIAI ĐOẠN 3: PHÂN TÍCH KẾT QUẢ ---\n\n    print(\"\\n--- Giai đoạn 3: Phân tích Kết quả Cuối cùng ---\")\n    results_df = pd.DataFrame.from_dict(p_values_dict, orient='index', columns=['raw_p_value']).dropna()\n    results_df.index.name = 'feature'\n    results_df.sort_values('raw_p_value', inplace=True)\n\n    if not results_df.empty:\n        reject, p_adjusted, _, _ = multipletests(results_df['raw_p_value'], alpha=0.05, method='fdr_bh')\n        results_df['adjusted_p_value'] = p_adjusted\n        results_df['is_significant'] = reject\n        \n        lasso_coef_df = pd.DataFrame({\n            'feature': X_columns[A_obs],\n            'lasso_coefficient': lasso_model.coef_[A_obs]\n        })\n        results_df = pd.merge(results_df, lasso_coef_df, left_index=True, right_on='feature').set_index('feature')\n        results_df = results_df.sort_values('adjusted_p_value')\n\n        print(\"\\nBảng kết quả Suy luận Chọn lọc:\")\n        print(results_df)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# def generate_data(n, p, true_beta, noise_std=1):\n#     \"\"\"Hàm tạo dữ liệu mô phỏng.\"\"\"\n#     print(f\"Đang tạo dữ liệu với n={n}, p={p}...\")\n#     X = np.random.normal(loc=0, scale=1, size=(n, p))\n#     true_beta_reshaped = np.reshape(true_beta, (p, 1))\n#     noise = np.random.normal(loc=0, scale=noise_std, size=(n, 1))\n#     y = X @ true_beta_reshaped + noise\n#     print(\"Tạo dữ liệu xong.\")\n#     return X, y\n\n# if __name__ == '__main__':\n#     # --- BƯỚC 0: CHUẨN BỊ DỮ LIỆU MÔ PHỎNG ---\n    \n#     # 0.1. Thiết lập các tham số cho mô phỏng\n#     n = 170000\n#     p = 434\n#     noise_std = 1\n#     num_true_features = 25 # Số lượng đặc trưng thực sự có ảnh hưởng\n    \n#     # 0.2. Tạo vector hệ số beta thực\n#     true_beta = np.zeros(p)\n#     # Chọn ngẫu nhiên các chỉ số cho các đặc trưng có ảnh hưởng\n#     true_feature_indices = np.random.choice(p, num_true_features, replace=False)\n#     # Gán giá trị cho các hệ số beta thực (ví dụ: 5)\n#     true_beta[true_feature_indices] = 1\n    \n#     print(f\"Các đặc trưng có ảnh hưởng thực sự (chỉ số): {sorted(list(true_feature_indices))}\")\n    \n#     # 0.3. Tạo dữ liệu X và y\n#     X_sim, y_sim = generate_data(n, p, true_beta, noise_std)\n    \n#     # 0.4 Tạo tên cột giả để tương thích với mã cũ\n#     X_columns = [f'feature_{i}' for i in range(p)]\n\n#     print(f\"Dữ liệu mô phỏng: X.shape={X_sim.shape}, y.shape={y_sim.shape}\")\n\n#     # --- GIAI ĐOẠN 1: LỰA CHỌN MÔ HÌNH BẰNG LASSO ---\n#     print(\"\\n--- Giai đoạn 1: Lựa chọn Mô hình bằng Lasso ---\")\n\n#     # 1.1. Chuẩn hóa X và y\n#     print(\"Đang chuẩn hóa X và y...\")\n#     scaler_X = StandardScaler()\n#     X_scaled = X_sim\n   \n#     scaler_y = StandardScaler(with_std=False) # Chỉ trừ trung bình, không chia cho std dev\n#     y_scaled = y_sim\n    \n#     n, p = X_scaled.shape\n    \n\n#     BEST_ALPHA = 6000\n  \n#     lasso_model = Lasso(alpha=BEST_ALPHA/n, fit_intercept=False, max_iter=80000, tol=1e-9).fit(X_scaled, y_scaled.ravel())\n#     A_obs, XA_obs = construct_A_XA_Ac_XAc_bhA(X_scaled, lasso_model.coef_)\n\n#     if not A_obs:\n#         print(\"Lasso không chọn đặc trưng nào. Kết thúc.\")\n#         exit()\n\n#     print(f\"\\nLasso đã chọn {len(A_obs)} đặc trưng.\")\n#     print(\"Các đặc trưng được chọn (chỉ số):\", A_obs)\n\n#     # 1.4. Ước lượng sigma^2\n#     print(\"Ước lượng sigma^2 từ phần dư của Lasso...\")\n#     y_pred_lasso = lasso_model.predict(X_scaled)\n#     residuals_lasso = y_scaled.flatten() - y_pred_lasso.flatten()\n#     rss_lasso = np.sum(residuals_lasso**2)\n#     df_residuals = n - len(A_obs) \n#     if df_residuals <= 0:\n#         raise ValueError(f\"Không thể ước lượng sigma^2 vì bậc tự do (n - |A_obs|) = ({n} - {len(A_obs)}) <= 0.\")\n#     sigma_squared_hat = rss_lasso / df_residuals\n#     print(f\"Ước lượng phương sai nhiễu (σ^2) : {sigma_squared_hat:.4f} (Giá trị thực tế là {noise_std**2})\")\n\n#     # --- GIAI ĐOẠN 2: SUY LUẬN CHỌN LỌC ---\n#     print(f\"\\n--- Giai đoạn 2: Bắt đầu Suy luận Chọn lọc cho {len(A_obs)} đặc trưng ---\")\n#     p_values_dict = {}\n    \n#     y_col_vector = y_scaled.reshape(-1, 1)\n\n#     for i, j_selected in enumerate(A_obs):\n#         feature_name = X_columns[j_selected]\n#         start_time = time.time()\n#         print(f\"  [{i+1}/{len(A_obs)}] Đang phân tích '{feature_name}' (chỉ số {j_selected})...\")\n        \n#         etaj, etajTy_orig = construct_test_statistic(j_selected, XA_obs, y_col_vector, A_obs)\n#         sq_norm_eta = (etaj.T @ etaj).item()\n#         if sq_norm_eta < 1e-12: sq_norm_eta = 1e-12\n        \n#         z_obs = etajTy_orig\n#         print(f\"    z_obs = {z_obs:.4f}\")\n        \n#         sigma_z = np.sqrt(sq_norm_eta * sigma_squared_hat)\n#         threshold = 20 * sigma_z\n#         if abs(z_obs) > threshold:\n#             print(f\"    z_obs ({abs(z_obs):.2f}) vượt ngưỡng 20*sigma_z ({threshold:.2f}). Đặt threshold = |z_obs|*1.5\")\n#             threshold = abs(z_obs) * 1.5\n#         print(f\"    Sử dụng threshold = {threshold:.4f}\")\n        \n#         # Sử dụng số core logic của CPU làm số đoạn xử lý\n#         n_segments = os.cpu_count() or 4 # Mặc định là 4 nếu không phát hiện được\n#         print(f\"    Bắt đầu `run_parametric_lasso` với {n_segments} tiến trình...\")\n        \n        \n#         list_zk, list_active_set, list_bhz = run_parametric_lasso(\n#             X_scaled, y_col_vector, BEST_ALPHA, etaj,\n#             threshold=threshold, sigma_sq=sigma_squared_hat,\n#             n_segments=n_segments\n#         )\n    \n#         print(f\"    `run_parametric_lasso` hoàn tất. Bắt đầu tính p-value...\")\n#         cov_matrix = np.identity(n) * sigma_squared_hat\n#         p_val = p_value(A_obs, list_active_set, list_zk, etaj, z_obs, cov_matrix)\n        \n#         p_values_dict[feature_name] = p_val\n#         elapsed_time = time.time() - start_time\n#         print(f\"    -> p-value = {p_val if p_val is not None else 'N/A'} (Thời gian: {elapsed_time:.2f}s)\")\n        \n#         # LƯU Ý: Lệnh break này sẽ dừng vòng lặp sau khi phân tích xong đặc trưng đầu tiên.\n#         # Xóa nó đi nếu bạn muốn phân tích TẤT CẢ các đặc trưng Lasso đã chọn.\n#         break\n        \n#     # --- GIAI ĐOẠN 3: PHÂN TÍCH KẾT QUẢ ---\n#     print(\"\\n--- Giai đoạn 3: Phân tích Kết quả Cuối cùng ---\")\n    \n#     # Lấy các hệ số Lasso cho các đặc trưng được chọn\n#     lasso_coef_df = pd.DataFrame({\n#         'feature': [X_columns[i] for i in A_obs],\n#         'lasso_coefficient': lasso_model.coef_[A_obs]\n#     })\n    \n#     # Tạo DataFrame kết quả từ p-values\n#     results_df = pd.DataFrame.from_dict(p_values_dict, orient='index', columns=['raw_p_value']).dropna()\n#     results_df.index.name = 'feature'\n    \n#     if not results_df.empty:\n#         # Hiệu chỉnh p-value cho đa kiểm định\n#         reject, p_adjusted, _, _ = multipletests(results_df['raw_p_value'], alpha=0.05, method='fdr_bh')\n#         results_df['adjusted_p_value'] = p_adjusted\n#         results_df['is_significant'] = reject\n        \n#         # Kết hợp với DataFrame hệ số Lasso\n#         results_df = pd.merge(results_df.reset_index(), lasso_coef_df, on='feature').set_index('feature')\n        \n#         # Thêm cột cho biết đây có phải là đặc trưng \"thật\" không\n#         true_feature_names = {X_columns[i] for i in true_feature_indices}\n#         results_df['is_true_feature'] = results_df.index.isin(true_feature_names)\n        \n#         # Sắp xếp lại\n#         results_df = results_df.sort_values('adjusted_p_value')\n\n#         print(\"\\nBảng kết quả Suy luận Chọn lọc:\")\n#         print(results_df[['lasso_coefficient', 'raw_p_value', 'adjusted_p_value', 'is_significant', 'is_true_feature']])\n#     else:\n#         print(\"Không có p-value nào được tính toán.\")","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}