{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.18","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\n\n# hours = range(0, 6, 1)\n# minutes = range(0, 60, 5)\n\n# target_col = \"bg+1-00\"\n# group_col = \"p_num\"\n# date_col = \"time\"\n\n# bg_cols = [f\"bg-{i}-{j:02d}\" for i, j in product(hours, minutes)]\n# insu_cols = [f\"insulin-{i}-{j:02d}\" for i, j in product(hours, minutes)]\n# carb_cols = [f\"carbs-{i}-{j:02d}\" for i, j in product(hours, minutes)]\n# hr_cols = [f\"hr-{i}-{j:02d}\" for i, j in product(hours, minutes)]\n# step_cols = [f\"steps-{i}-{j:02d}\" for i, j in product(hours, minutes)]\n# cals_cols = [f\"cals-{i}-{j:02d}\" for i, j in product(hours, minutes)]\n\n# feature_cols = bg_cols + insu_cols + carb_cols + hr_cols + step_cols + cals_cols\n\n\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-08-23T14:25:47.117203Z","iopub.execute_input":"2025-08-23T14:25:47.117755Z","iopub.status.idle":"2025-08-23T14:25:53.483039Z","shell.execute_reply.started":"2025-08-23T14:25:47.117725Z","shell.execute_reply":"2025-08-23T14:25:53.479314Z"}},"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,"execution":{"iopub.status.busy":"2025-08-23T14:25:53.486726Z","iopub.execute_input":"2025-08-23T14:25:53.487093Z","iopub.status.idle":"2025-08-23T14:25:53.495514Z","shell.execute_reply.started":"2025-08-23T14:25:53.487068Z","shell.execute_reply":"2025-08-23T14:25:53.490945Z"}},"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,"execution":{"iopub.status.busy":"2025-08-23T14:25:53.497339Z","iopub.execute_input":"2025-08-23T14:25:53.497558Z","iopub.status.idle":"2025-08-23T14:25:54.741274Z","shell.execute_reply.started":"2025-08-23T14:25:53.497538Z","shell.execute_reply":"2025-08-23T14:25:54.735201Z"}},"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,"execution":{"iopub.status.busy":"2025-08-23T14:25:54.742185Z","iopub.execute_input":"2025-08-23T14:25:54.742544Z","iopub.status.idle":"2025-08-23T14:25:54.759478Z","shell.execute_reply.started":"2025-08-23T14:25:54.742515Z","shell.execute_reply":"2025-08-23T14:25:54.754239Z"}},"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,"execution":{"iopub.status.busy":"2025-08-23T14:25:54.760509Z","iopub.execute_input":"2025-08-23T14:25:54.760727Z","iopub.status.idle":"2025-08-23T14:25:54.776462Z","shell.execute_reply.started":"2025-08-23T14:25:54.760706Z","shell.execute_reply":"2025-08-23T14:25:54.772291Z"}},"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,"execution":{"iopub.status.busy":"2025-08-23T14:25:54.778527Z","iopub.execute_input":"2025-08-23T14:25:54.778938Z","iopub.status.idle":"2025-08-23T14:25:54.786834Z","shell.execute_reply.started":"2025-08-23T14:25:54.778917Z","shell.execute_reply":"2025-08-23T14:25:54.782469Z"}},"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,"execution":{"iopub.status.busy":"2025-08-23T14:25:54.788336Z","iopub.execute_input":"2025-08-23T14:25:54.788537Z","iopub.status.idle":"2025-08-23T14:25:54.799831Z","shell.execute_reply.started":"2025-08-23T14:25:54.788518Z","shell.execute_reply":"2025-08-23T14:25:54.79519Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# df_train_final","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-23T14:25:54.803154Z","iopub.execute_input":"2025-08-23T14:25:54.803491Z","iopub.status.idle":"2025-08-23T14:25:54.815739Z","shell.execute_reply.started":"2025-08-23T14:25:54.803467Z","shell.execute_reply":"2025-08-23T14:25:54.810059Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"!pip install skglm\n!pip install statsmodels\n# !pip install celer\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-23T14:25:54.817908Z","iopub.execute_input":"2025-08-23T14:25:54.818344Z","iopub.status.idle":"2025-08-23T14:26:13.05831Z","shell.execute_reply.started":"2025-08-23T14:25:54.818321Z","shell.execute_reply":"2025-08-23T14:26:13.053874Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"#1.SI\n","metadata":{}},{"cell_type":"code","source":"import 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\nfrom scipy.stats import norm\n\nnp.random.seed(2025)\n\ndef generate_data(n, p, true_beta):\n\n    X = np.random.normal(loc=0, scale=1, size=(n, p))\n    true_beta = np.reshape(true_beta, (p, 1))\n    noise = np.random.normal(loc=0, scale=1, size=(n, 1))\n    y = X @ true_beta + noise\n    return X, y\n\nmp.dps =50\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        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    segment_results = parallel(jobs)\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 (set(merged_active_sets[-1])==set(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# def p_value(A_obs, list_active_set, list_zk, etaj, z_obs, cov_matrix):\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(\"Done7\")\n    # return 2 * min(pivot_val, 1 - pivot_val)\n\nimport numpy as np\n\ndef p_value(A_obs, list_active_set, list_zk, etaj, z_obs, cov_matrix,\n            extreme_threshold=8.0, log_underflow_thresh=-745.0):\n    \"\"\"\n    Robust p-value calculation:\n    - compute intervals where active_set == A_obs\n    - work in log-space with logsumexp\n    - for extreme tails (a_s > extreme_threshold) use approximation:\n        log P(a<b) ≈ log phi(a_s) + log(b_s - a_s)\n    - fallback to mpmath only if needed\n    \"\"\"\n    from math import log, pi\n    tn_sigma = (np.sqrt((etaj.T @ cov_matrix @ etaj))).item()\n    print(\"tn_sigma:\", tn_sigma)\n    if tn_sigma < 1e-12:\n        return None\n\n    # collect intervals (float)\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((float(list_zk[i]), float(list_zk[i+1])))\n    print(\"z_intervals:\", z_intervals)\n    if not z_intervals:\n        return None\n\n    # try scipy logcdf first\n    try:\n        from scipy.stats import norm\n        from scipy.special import logsumexp\n        use_scipy = True\n    except Exception:\n        use_scipy = False\n\n    log_probs = []\n    log_partial = []\n\n    def log_phi(x):\n        # log density of standard normal at x\n        return -0.5 * (x**2) - 0.5 * np.log(2.0 * np.pi)\n\n    if use_scipy:\n        def log_cdf_diff_scipy(a, b):\n            # robust log(Φ(b) - Φ(a)) using logcdf\n            # assume a <= b\n            if np.isneginf(a) and np.isposinf(b):\n                return 0.0\n            if np.isneginf(a):\n                return norm.logcdf(b)\n            if np.isposinf(b):\n                la = norm.logcdf(a)\n                return np.log1p(-np.exp(la)) if la < 0 else np.log(1.0 - np.exp(la))\n            La = norm.logcdf(a)\n            Lb = norm.logcdf(b)\n            if Lb < La:\n                La, Lb = Lb, La\n            # Lb + log(1 - exp(La-Lb))\n            return Lb + np.log1p(-np.exp(La - Lb))\n\n        for (al, ar) in z_intervals:\n            a_s = al / tn_sigma\n            b_s = ar / tn_sigma\n            width_s = b_s - a_s\n\n            # if in extreme tail, use density*width approx (in log-space) to avoid underflow\n            if a_s > extreme_threshold:\n                # log prob ≈ log phi(a_s) + log(width_s)\n                lp = log_phi(a_s) + np.log(width_s)\n            else:\n                # normal route with logcdf diff\n                lp = log_cdf_diff_scipy(a_s, b_s)\n                # if this gave -inf or nan, fallback to density approx if width small\n                if not np.isfinite(lp):\n                    lp = log_phi(a_s) + np.log(max(width_s, 1e-300))\n\n            log_probs.append(lp)\n\n            # numerator piece\n            if z_obs >= ar:\n                log_partial.append(lp)\n            elif z_obs > al:\n                z_s = float(z_obs) / tn_sigma\n                if a_s > extreme_threshold:\n                    # partial width: (z_s - a_s)\n                    part_width = z_s - a_s\n                    lp_part = log_phi(a_s) + np.log(max(part_width, 1e-300))\n                else:\n                    lp_part = log_cdf_diff_scipy(a_s, z_s)\n                    if not np.isfinite(lp_part):\n                        lp_part = log_phi(a_s) + np.log(max(z_s - a_s, 1e-300))\n                log_partial.append(lp_part)\n            else:\n                # no contrib\n                pass\n\n        # aggregate\n        finite_logs = [v for v in log_probs if np.isfinite(v)]\n        if len(finite_logs) == 0:\n            use_scipy = False\n        else:\n            log_den = logsumexp(log_probs)\n            log_num = logsumexp([v for v in log_partial if np.isfinite(v)]) if any(np.isfinite(x) for x in log_partial) else -np.inf\n            print(\"log_den (scipy-approx):\", log_den, \"log_num:\", log_num)\n            # if denominator extremely small -> fallback\n            # Nếu cả hai log đều hữu hạn thì ta có thể tính pivot ngay bằng exp(log_num - log_den)\n            if np.isfinite(log_den) and np.isfinite(log_num):\n                pivot_val = float(np.exp(log_num - log_den))\n                pivot_val = max(0.0, min(1.0, pivot_val))\n                pval = 2.0 * min(pivot_val, 1.0 - pivot_val)\n                return pval\n            # Nếu một trong hai không hữu hạn, mới fallback sang mpmath\n            use_scipy = False\n\n    # fallback to high-precision mpmath if needed\n    try:\n        import mpmath as mp\n        if mp.mp.dps < 80:\n            mp.mp.dps = 80\n    except Exception:\n        mp = None\n\n    if mp is None:\n        # no high-precision available; return extremely small p-value (degenerate)\n        # better to return 0.0 than crash; caller can inspect logs for log_den\n        return 0.0\n\n    # mpmath computations\n    def mp_norm_cdf(x):\n        xm = mp.mpf(x)\n        return mp.mpf('0.5') * (1 + mp.erf(xm / mp.sqrt(2)))\n\n    denom_mp = mp.mpf('0')\n    numer_mp = mp.mpf('0')\n    for (al, ar) in z_intervals:\n        a_s = mp.mpf(al) / mp.mpf(tn_sigma)\n        b_s = mp.mpf(ar) / mp.mpf(tn_sigma)\n        prob = mp_norm_cdf(b_s) - mp_norm_cdf(a_s)\n        denom_mp += prob\n        if mp.mpf(z_obs) >= mp.mpf(ar):\n            numer_mp += prob\n        elif mp.mpf(z_obs) > mp.mpf(al):\n            z_s = mp.mpf(z_obs) / mp.mpf(tn_sigma)\n            numer_mp += mp_norm_cdf(z_s) - mp_norm_cdf(a_s)\n\n    if denom_mp == 0:\n        return 0.0\n\n    pivot_val_mp = numer_mp / denom_mp\n    try:\n        pivot_val = float(pivot_val_mp)\n    except Exception:\n        pivot_val = 0.0 if pivot_val_mp == 0 else 1.0\n    pivot_val = max(0.0, min(1.0, pivot_val))\n    pval = 2.0 * min(pivot_val, 1.0 - pivot_val)\n    print(\"Done (mpmath fallback).\")\n    return pval\n\ndef ensure_series(y, name=\"y\"):\n    import numpy as np, pandas as pd\n    # DataFrame 1 cột -> chuyển thành Series\n    if isinstance(y, pd.DataFrame):\n        if y.shape[1] == 1:\n            return y.iloc[:, 0].copy()\n        else:\n            raise ValueError(\"y là DataFrame nhiều cột — cần chỉ rõ cột target.\")\n    # ndarray\n    if isinstance(y, np.ndarray):\n        if y.ndim == 1:\n            return pd.Series(y, name=name)\n        if y.ndim == 2:\n            # trường hợp (n,1) hoặc (1,n)\n            if y.shape[1] == 1 or y.shape[0] == 1:\n                return pd.Series(y.ravel(), name=name)\n            else:\n                raise ValueError(f\"ndarray không phải vector 1D hay cột 1: shape={y.shape}\")\n    # nếu là list, tuple, pd.Series, v.v.\n    return pd.Series(y, name=name)\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\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 = 177000\n    p = 433\n    true_beta = [0.25 if i < 25 else 0 for i in range(p)]\n\n    X_clean, y_clean = generate_data(n, p, true_beta)  \n    # if isinstance(X_clean, np.ndarray):\n    #     col_names = [f\"x{i}\" for i in range(X_clean.shape[1])]  # hoặc dùng tên khác\n    #     X_clean = pd.DataFrame(X_clean, columns=col_names)\n    # else:\n    #     # nếu đã là DataFrame, copy để tránh thay đổi ngoài ý muốn\n    #     X_clean = X_clean.copy()\n\n    # y_clean = ensure_series(y_clean, name=\"y\")\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    import pandas as pd\n\n    # Nếu X_clean là numpy array 2D\n    X_clean = pd.DataFrame(X_clean, columns=[f\"x{i}\" for i in range(X_clean.shape[1])])\n    \n    # Nếu y_clean là numpy array 1D hoặc 2D\n    y_clean = pd.DataFrame(y_clean, columns=[\"y\"])\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    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==5):\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,"execution":{"iopub.status.busy":"2025-08-23T14:26:13.061151Z","iopub.execute_input":"2025-08-23T14:26:13.061397Z","iopub.status.idle":"2025-08-23T14:36:03.262064Z","shell.execute_reply.started":"2025-08-23T14:26:13.061372Z","shell.execute_reply":"2025-08-23T14:36:03.254728Z"}},"outputs":[],"execution_count":null}]}