{"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"}],"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\nfrom sklearn.linear_model import Ridge\nfrom sklearn.ensemble import StackingRegressor\nwarnings.filterwarnings(\"ignore\")\n\nfrom itertools import product\n\nhours = range(0, 6, 1)\nminutes = range(0, 60, 5)\n\ntarget_col = \"bg+1-00\"\ngroup_col = \"p_num\"\ndate_col = \"time\"\n\nbg_cols = [f\"bg-{i}-{j:02d}\" for i, j in product(hours, minutes)]\ninsu_cols = [f\"insulin-{i}-{j:02d}\" for i, j in product(hours, minutes)]\ncarb_cols = [f\"carbs-{i}-{j:02d}\" for i, j in product(hours, minutes)]\nhr_cols = [f\"hr-{i}-{j:02d}\" for i, j in product(hours, minutes)]\nstep_cols = [f\"steps-{i}-{j:02d}\" for i, j in product(hours, minutes)]\ncals_cols = [f\"cals-{i}-{j:02d}\" for i, j in product(hours, minutes)]\n\nfeature_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-11T14:23:33.766453Z","iopub.execute_input":"2025-08-11T14:23:33.766808Z","iopub.status.idle":"2025-08-11T14:23:40.412166Z","shell.execute_reply.started":"2025-08-11T14:23:33.766785Z","shell.execute_reply":"2025-08-11T14:23:40.40799Z"}},"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\ndf_train.columns = df_train.columns.str.replace(':', '-')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-11T14:23:40.415418Z","iopub.execute_input":"2025-08-11T14:23:40.415837Z","iopub.status.idle":"2025-08-11T14:23:51.212603Z","shell.execute_reply.started":"2025-08-11T14:23:40.415811Z","shell.execute_reply":"2025-08-11T14:23:51.207751Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import random\nimport os\nimport numpy as np\nimport torch\n\n\ndef 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\ndf_train.columns = df_train.columns.str.replace(\":\", \"-\")\n\n\nfor 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\nfeature_columns = {\n    'bg': bg_cols,\n    'hr': hr_cols,\n    'steps': step_cols\n}\n\ndef 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\nfor feature, columns in feature_columns.items():\n    df_train = impute_time_series_features(df_train, columns, method='median')\n\nall_imputed_columns = bg_cols + hr_cols + step_cols\nprint(\"Remaining null values in training set:\", df_train[all_imputed_columns].isnull().sum().sum())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-11T14:23:51.215278Z","iopub.execute_input":"2025-08-11T14:23:51.215538Z","iopub.status.idle":"2025-08-11T14:24:50.224855Z","shell.execute_reply.started":"2025-08-11T14:23:51.215515Z","shell.execute_reply":"2025-08-11T14:24:50.219028Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\nfeature_columns_insulin_cals = {\n    'insulin': insu_cols,\n    'cals': cals_cols\n}\n\ndef 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\nfor feature, columns in feature_columns_insulin_cals.items():\n    df_train = impute_non_temporal_features(df_train, columns, method='median')\n\nall_imputed_columns_insulin_cals = insu_cols + cals_cols\nprint(\"Remaining null values in training set:\", df_train[all_imputed_columns_insulin_cals].isnull().sum().sum())\n\n\n\nfeature_columns_insulin_cals = {\n    'insulin': insu_cols,\n    'cals': cals_cols\n}\n\ndef 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\nfor feature, columns in feature_columns_insulin_cals.items():\n    df_train = impute_non_temporal_features(df_train, columns, method='median')\n\nall_imputed_columns_insulin_cals = insu_cols + cals_cols\nprint(\"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-11T14:24:50.225669Z","iopub.execute_input":"2025-08-11T14:24:50.226059Z","iopub.status.idle":"2025-08-11T14:25:08.781727Z","shell.execute_reply.started":"2025-08-11T14:24:50.226033Z","shell.execute_reply":"2025-08-11T14:25:08.777034Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"activity_cols = [col for col in df_train.columns if col.startswith('activity')]\ndominant_activity_train = df_train[activity_cols].mode(axis=1)[0]\ndf_train['dominant_activity'] = dominant_activity_train\nprint(df_train['dominant_activity'].value_counts())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-11T14:25:08.783754Z","iopub.execute_input":"2025-08-11T14:25:08.784088Z","iopub.status.idle":"2025-08-11T14:25:49.515291Z","shell.execute_reply.started":"2025-08-11T14:25:08.784065Z","shell.execute_reply":"2025-08-11T14:25:49.510874Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"activity_replace = {\n    'Walk': 'Walking',\n    'Run': 'Running',\n    'Swim': 'Swimming',\n}\ndf_train['dominant_activity'] = df_train['dominant_activity'].replace(activity_replace)\nprint(df_train['dominant_activity'])","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-11T14:25:49.518495Z","iopub.execute_input":"2025-08-11T14:25:49.518742Z","iopub.status.idle":"2025-08-11T14:25:49.547822Z","shell.execute_reply.started":"2025-08-11T14:25:49.51872Z","shell.execute_reply":"2025-08-11T14:25:49.543424Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"df_train['dominant_activity'].fillna('No Activity', inplace=True)\n\nfrom sklearn.preprocessing import LabelEncoder\n\nlabel_encoder = LabelEncoder()\n\ncombined_activities = pd.concat([\n    df_train['dominant_activity']\n])\n\nlabel_encoder.fit(combined_activities)\n\ndf_train['dominant_activity_encoded'] = label_encoder.transform(df_train['dominant_activity'])\n\ndf_train = df_train.drop(['dominant_activity'], axis = 1)\n\nfeature_cols = feature_cols + ['dominant_activity_encoded']\n\n\nimputer = SimpleImputer()\n\ndf_train[feature_cols] = imputer.fit_transform(df_train[feature_cols])\n\ndf_train['time_hour'] = pd.to_datetime(df_train['time']).dt.hour\n\nfeature_cols.extend([\"time_hour\"])\n\ndf_train_final = df_train[feature_cols]\n\ny_target = df_train[[target_col]]\n\n\nnp.isnan(df_train_final).sum()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-11T14:25:49.550334Z","iopub.execute_input":"2025-08-11T14:25:49.550621Z","iopub.status.idle":"2025-08-11T14:25:53.170964Z","shell.execute_reply.started":"2025-08-11T14:25:49.550581Z","shell.execute_reply":"2025-08-11T14:25:53.166398Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"df_train_final","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-11T14:25:53.173551Z","iopub.execute_input":"2025-08-11T14:25:53.173804Z","iopub.status.idle":"2025-08-11T14:25:53.372472Z","shell.execute_reply.started":"2025-08-11T14:25:53.173781Z","shell.execute_reply":"2025-08-11T14:25:53.366872Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"y_target","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-11T14:25:53.37428Z","iopub.execute_input":"2025-08-11T14:25:53.374507Z","iopub.status.idle":"2025-08-11T14:25:53.392116Z","shell.execute_reply.started":"2025-08-11T14:25:53.374486Z","shell.execute_reply":"2025-08-11T14:25:53.386127Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"#1.SI\n","metadata":{}},{"cell_type":"code","source":"!pip install statsmodels","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-11T14:25:53.393811Z","iopub.execute_input":"2025-08-11T14:25:53.394038Z","iopub.status.idle":"2025-08-11T14:26:01.901145Z","shell.execute_reply.started":"2025-08-11T14:25:53.394017Z","shell.execute_reply":"2025-08-11T14:26:01.894291Z"}},"outputs":[],"execution_count":null},{"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\nfrom sklearn.linear_model 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 = 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)[0, 0]\n\n    return etaj, etajTy\n\n\n\ndef compute_yz(y, etaj, zk, n, sigma_squared):\n\n    sq_norm = sigma_squared*(etaj.T @ etaj)\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)\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        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\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\ndef parametric_lasso(X, yz, lamda, b, n, p,eta_j):\n    yz_flatten = yz.flatten()\n    clf = linear_model.Lasso(alpha=lamda/n, fit_intercept=False, tol=1e-10)\n    clf.fit(X, yz_flatten)\n    bhz = clf.coef_\n\n    A, XA = construct_A_XA_Ac_XAc_bhA(X, bhz)\n    s=np.sign(bhz[A])\n    A1,b1= compute_polyhedron(X,A,s, lamda)\n    V_l,V_r=compute_truncation_bounds(A1, b1, eta_j, yz_flatten)\n    print(\"Done5\")\n    print(\"V_l: \",V_l)\n    print(\"V_r: \",V_r)\n\n    return V_l,V_r, A, bhz\n\n# SỬA ĐỔI 2: Hàm run_... giờ nhận sigma_squared\n# def run_parametric_lasso(X, y, alpha, etaj, threshold, sigma_squared):\n#     n, p = X.shape\n#     zk = -threshold\n#     list_zk = [zk]\n#     list_active_set = []\n#     list_bhz = []\n    \n#     while zk < threshold:\n#         yz, b = compute_yz(y, etaj, zk, n,sigma_squared)\n#         V_l,V_r, Akz, bhkz = parametric_lasso(X, yz, alpha, b, n, p,etaj)\n#         zk = V_r + 0.0001\n\n\n#         if zk <= threshold:\n#             list_zk.append(zk-0.0001)\n#         else:\n#             list_zk.append(threshold)\n#         list_active_set.append(Akz)\n#         list_bhz.append(bhkz)\n\n#     print(\"DOne6\")\n#     return list_zk, list_bhz, list_active_set\n\n\n\n\n\n\n\n\n\n\n# ============================\n# Hàm tạo/attach shared memory\n# ============================\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# def 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, small_eps=1e-6, max_iter=2000,\n#                     lasso_tol=1e-10, lasso_max_iter=20000):\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, tol=lasso_tol, warm_start=True, max_iter=lasso_max_iter)\n\n#     zk = float(z_start)\n#     it = 0\n#     while zk < z_end and it < max_iter:\n#         it += 1\n#         yz, b_vec = compute_yz(y, etaj, zk, n, sigma_sq)\n#         yz_flat = yz.flatten()\n\n#         clf.fit(X, yz_flat)\n#         bhz = clf.coef_.copy()\n\n#         A, XA= 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#             try:\n#                 V_l, V_r = compute_truncation_bounds(A1, b1, etaj, yz_flat)\n#             except Exception:\n#                 V_l, V_r = zk, zk + small_eps\n\n    #     results.append((zk, A, bhz))\n\n    #     if not np.isfinite(V_r) or V_r <= zk + 1e-12:\n    #         zk = zk + small_eps\n    #     else:\n    #         zk = min(V_r + small_eps, z_end)\n\n    # shmX.close()\n    # shmy.close()\n    # shmetaj.close()\n\n    # return results\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, small_eps=0.001, max_iter=2000,\n                    lasso_tol=1e-4, lasso_max_iter=20000):\n    \n    # --- SỬA ĐỔI 1: Thiết lập logger và log khi bắt đầu ---\n    worker_id = os.getpid()\n    logger = setup_worker_logging(worker_id)\n    logger.info(f\"BẮT ĐẦU xử lý đoạn [{z_start:.4f}, {z_end:.4f}]\")\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, tol=lasso_tol, warm_start=True, max_iter=lasso_max_iter)\n\n    zk = float(z_start)\n    it = 0\n    # --- SỬA ĐỔI 2: Tạo biến để kiểm soát tần suất log, tránh làm ngập màn hình ---\n    log_interval = 20  # Chỉ in log sau mỗi 250 vòng lặp\n\n    while zk < z_end and it < max_iter:\n        it += 1\n        \n        # --- SỬA ĐỔI 3: Log \"nhịp tim\" của worker ---\n        # Điều này cho bạn biết worker có bị \"treo\" hay không.\n        if it == 1 or it % log_interval == 0:\n            logger.info(f\"Vòng lặp #{it}... zk hiện tại = {zk:.6f}\")\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\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            try:\n                V_l, V_r = compute_truncation_bounds(A1, b1, etaj, yz_flat)\n            except Exception as e:\n                logger.error(f\"Lỗi tại vòng lặp {it}, zk={zk:.4f}: {e}. Dùng bước nhảy nhỏ.\")\n                V_l, V_r = zk, zk + small_eps\n\n        results.append((zk, A, bhz))\n        \n        # --- SỬA ĐỔI 4 (Tùy chọn): Log kết quả của mỗi bước nhảy ---\n        # Bỏ comment dòng dưới nếu bạn muốn xem chi tiết tuyệt đối, nhưng nó sẽ tạo ra RẤT NHIỀU log.\n        # logger.debug(f\"Loop {it}: zk={zk:.4f} -> new V_r={V_r:.4f}. Active set size: {len(A)}\")\n        # Lưu ý: Phải đổi logger.setLevel(logging.DEBUG) ở trên để thấy log này.\n\n        if not np.isfinite(V_r) or V_r <= zk + 1e-12:\n            zk = zk + small_eps\n        else:\n            zk = min(V_r + small_eps, z_end)\n\n    shmX.close()\n    shmy.close()\n    shmetaj.close()\n    \n    # --- SỬA ĐỔI 5: Log khi kết thúc và báo cáo tổng kết ---\n    end_time = time.time()\n    duration = end_time - start_time\n    \n    if it >= max_iter:\n        logger.warning(f\"ĐÃ ĐẠT GIỚI HẠN {max_iter} vòng lặp. Worker có thể chưa xử lý xong.\")\n\n    logger.info(f\"KẾT THÚC xử lý đoạn [{z_start:.4f}, {z_end:.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=5, overlap=1e-8, **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    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                **worker_kwargs\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 = 1e-8\n    for zk, A, bhz in all_entries:\n        if merged_zk and abs(zk - merged_zk[-1]) < tol:\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\n# SỬA ĐỔI 3: Hàm p_value nhận z_obs và đã sửa lỗi logic\ndef p_value(A_obs, list_active_set, list_zk, etaj, z_obs, cov_matrix):\n    tn_sigma = np.sqrt((etaj.T @ cov_matrix @ etaj))[0, 0]\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    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            \n    if denominator < 1e-100: return None\n    \n    pivot_val = float(numerator / denominator)\n    print(\"Done7\")\n    return 2 * min(pivot_val, 1 - pivot_val)\n\n# =============================================================================\n# HÀM CHẠY CHÍNH: QUY TRÌNH PHÂN TÍCH DỮ LIỆU THỰC\n# =============================================================================\n\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    X_clean =df_train_final\n    y_clean =y_target\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, p = X_scaled.shape\n    BEST_ALPHA = 6000\n    print(f\"Sử dụng alpha cố định: {BEST_ALPHA:.4f}\")\n    lasso_model = Lasso(alpha=BEST_ALPHA/n).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    # SỬA ĐỔI 4: Ước lượng phương sai nhiễu (σ²)\n\n    \n    y_pred = lasso_model.predict(X_scaled)\n    df_residuals = n - len(A_obs)\n    if df_residuals <= 0:\n        raise ValueError(\"Không thể ước lượng sigma^2.\")\n    # sigma_squared_hat = np.sum((y_clean.values - y_pred)**2) / df_residuals\n    sigma_squared_hat = 1\n    print(f\"Ước lượng phương sai nhiễu (σ^2): {sigma_squared_hat:.4f}\")\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.values.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}'...\")\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        # Tối ưu hóa việc tính a_vec, tránh tạo ma trận n x n\n        # a_vec = y_col_vector - (etaj @ etaj.T @ y_col_vector) / sq_norm_eta # Cách không hiệu quả\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        sigma_z = np.sqrt(sq_norm_eta * sigma_squared_hat)\n        # threshold = 20 * sigma_z\n        threshold = 20 * 1\n        \n        # list_zk, _, list_active_set = run_parametric_lasso(\n        #     X_scaled, y_col_vector, BEST_ALPHA, etaj, threshold, sigma_squared_hat\n        # )\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=20.0, sigma_sq=1.0,\n            n_segments=n_segments,  # CHỈ THÊM DÒNG NÀY\n            small_eps=0.001, max_iter=2000\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-11T14:26:01.903079Z","iopub.execute_input":"2025-08-11T14:26:01.903382Z"}},"outputs":[],"execution_count":null}]}