{"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\nfrom sklearn.linear_model import Ridge\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-17T06:09:48.210786Z","iopub.execute_input":"2025-08-17T06:09:48.211207Z","iopub.status.idle":"2025-08-17T06:09:55.468681Z","shell.execute_reply.started":"2025-08-17T06:09:48.211183Z","shell.execute_reply":"2025-08-17T06:09:55.464535Z"}},"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-17T06:09:55.472608Z","iopub.execute_input":"2025-08-17T06:09:55.47311Z","iopub.status.idle":"2025-08-17T06:09:55.484568Z","shell.execute_reply.started":"2025-08-17T06:09:55.473082Z","shell.execute_reply":"2025-08-17T06:09:55.478497Z"}},"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-17T06:09:55.485924Z","iopub.execute_input":"2025-08-17T06:09:55.486149Z","iopub.status.idle":"2025-08-17T06:09:55.505277Z","shell.execute_reply.started":"2025-08-17T06:09:55.486127Z","shell.execute_reply":"2025-08-17T06:09:55.498906Z"}},"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-17T06:09:55.506025Z","iopub.execute_input":"2025-08-17T06:09:55.506241Z","iopub.status.idle":"2025-08-17T06:09:55.52094Z","shell.execute_reply.started":"2025-08-17T06:09:55.506219Z","shell.execute_reply":"2025-08-17T06:09:55.515974Z"}},"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-17T06:09:55.523582Z","iopub.execute_input":"2025-08-17T06:09:55.523813Z","iopub.status.idle":"2025-08-17T06:09:55.54053Z","shell.execute_reply.started":"2025-08-17T06:09:55.523792Z","shell.execute_reply":"2025-08-17T06:09:55.5368Z"}},"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-17T06:09:55.543313Z","iopub.execute_input":"2025-08-17T06:09:55.54361Z","iopub.status.idle":"2025-08-17T06:09:55.555586Z","shell.execute_reply.started":"2025-08-17T06:09:55.543588Z","shell.execute_reply":"2025-08-17T06:09:55.552053Z"}},"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-17T06:09:55.556642Z","iopub.execute_input":"2025-08-17T06:09:55.556845Z","iopub.status.idle":"2025-08-17T06:09:55.573929Z","shell.execute_reply.started":"2025-08-17T06:09:55.556825Z","shell.execute_reply":"2025-08-17T06:09:55.569248Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# df_train_final","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-17T06:09:55.57586Z","iopub.execute_input":"2025-08-17T06:09:55.576081Z","iopub.status.idle":"2025-08-17T06:09:55.595077Z","shell.execute_reply.started":"2025-08-17T06:09:55.576061Z","shell.execute_reply":"2025-08-17T06:09:55.58955Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# y_target","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-17T06:09:55.59737Z","iopub.execute_input":"2025-08-17T06:09:55.597579Z","iopub.status.idle":"2025-08-17T06:09:55.611057Z","shell.execute_reply.started":"2025-08-17T06:09:55.597559Z","shell.execute_reply":"2025-08-17T06:09:55.606316Z"}},"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-17T06:09:55.613566Z","iopub.execute_input":"2025-08-17T06:09:55.613788Z","iopub.status.idle":"2025-08-17T06:10:04.072376Z","shell.execute_reply.started":"2025-08-17T06:09:55.613767Z","shell.execute_reply":"2025-08-17T06:10:04.066541Z"}},"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 =30\n\nimport numpy as np\nfrom sklearn.linear_model import Lasso\n\ndef compute_a_b_from_etaj(y, etaj):\n    \"\"\"\n    Trả về a, b sao cho y(z) = a + b * z\n    Công thức: a = y - eta (eta^T y) / ||eta||^2\n             b = eta / ||eta||^2\n    y: vector kích thước (n,) hoặc (n,1)\n    etaj: vector (n,1) hoặc (n,)\n    \"\"\"\n    y_vec = y.reshape(-1)\n    eta = etaj.reshape(-1)\n    eta_norm2 = float(np.dot(eta, eta))\n    if eta_norm2 < 1e-12:\n        eta_norm2 = 1e-12\n    etaTy = float(np.dot(eta, y_vec))\n    a = y_vec - eta * (etaTy / eta_norm2)\n    b = eta / eta_norm2\n    return a, b\n\ndef _initial_lasso_active_set(X, yz, lambda_val):\n    \"\"\"\n    Fit one Lasso to get initial active set and signs.\n    lambda_val here should follow same convention as you use (i.e., the 'lambda' in your KKT).\n    sklearn.Lasso uses alpha = lambda_val / n for loss (if your earlier code uses BEST_ALPHA/n).\n    \"\"\"\n    n = X.shape[0]\n    clf = Lasso(alpha=lambda_val / n, fit_intercept=False, max_iter=100000, tol=1e-8)\n    clf.fit(X, yz)\n    bh = clf.coef_.copy()\n    # threshold tiny values\n    bh[np.abs(bh) < 1e-12] = 0.0\n    A = list(np.where(bh != 0)[0])\n    s = np.sign(bh[A]) if len(A) > 0 else np.array([], dtype=float)\n    return A, s, bh\n\ndef lasso_path_tracker(X, a, b, lambda_val, z0, z1, max_iter=10000, tol=1e-12, verbose=False):\n    \"\"\"\n    Path-following/homotopy tracker for Lasso with parametric response y(z) = a + b*z.\n    Returns:\n      zk_list: list of breakpoints z (ascending)\n      active_sets_list: list of active sets (list of indices) at each interval start zk[i]\n      betas_list: list of full beta arrays (length p) evaluated at zk[i]\n    Notes:\n      - lambda_val is the lambda used in KKT: condition |X_j^T r(z)| <= lambda_val\n      - The function finds the NEXT breakpoint exactly by solving linear eqns.\n    \"\"\"\n    n, p = X.shape\n    a = a.reshape(-1)\n    b = b.reshape(-1)\n    # Precompute Gram and X^T a, X^T b\n    G = X.T @ X        # (p,p)\n    Xt_a = X.T @ a     # (p,)\n    Xt_b = X.T @ b     # (p,)\n\n    # initial Lasso at z0 to get active set (fallback)\n    yz0 = a + b * z0\n    A, s_vec, bh0 = _initial_lasso_active_set(X, yz0, lambda_val)\n    s_vec = np.array(s_vec, dtype=float)\n\n    zk_list = [float(z0)]\n    active_sets_list = [list(A)]\n    betas_list = [bh0.copy()]\n\n    current_z = float(z0)\n    it = 0\n\n    while current_z < z1 - tol and it < max_iter:\n        it += 1\n        A = list(active_sets_list[-1])\n        k = len(A)\n\n        # compute u0,u1 for beta_A(z) = u0 + z * u1\n        if k == 0:\n            # beta empty vector; handle residual & correlations easily\n            u0 = np.zeros(0)\n            u1 = np.zeros(0)\n            XA_u0 = np.zeros(n)\n            XA_u1 = np.zeros(n)\n        else:\n            G_AA = G[np.ix_(A, A)]\n            # RHS: X_A^T a - lambda * s\n            rhs0 = Xt_a[A] - lambda_val * s_vec\n            rhs1 = Xt_b[A]\n            # Solve G_AA * u0 = rhs0 and G_AA * u1 = rhs1\n            # prefer solve; if singular use pinv fallback\n            try:\n                u0 = np.linalg.solve(G_AA, rhs0)\n                u1 = np.linalg.solve(G_AA, rhs1)\n            except np.linalg.LinAlgError:\n                # use pseudo-inverse\n                G_AA_pinv = np.linalg.pinv(G_AA)\n                u0 = G_AA_pinv @ rhs0\n                u1 = G_AA_pinv @ rhs1\n            # precompute XA @ u0 and XA @ u1 for residual\n            XA = X[:, A]\n            XA_u0 = XA @ u0\n            XA_u1 = XA @ u1\n\n        # residual r(z) = a + b*z - X_A beta_A(z) = (a - X_A u0) + z*(b - X_A u1)\n        r0 = a - XA_u0    # vector\n        r1 = b - XA_u1    # vector\n\n        # For each active coordinate k_i solve u0_i + z*u1_i = 0 => z = -u0_i / u1_i\n        candidate_z = []\n        candidate_events = []  # (\"drop\", index) or (\"add\", index, sign)\n        eps = 1e-12\n\n        if k > 0:\n            for idx_in_A, var_idx in enumerate(A):\n                ui0 = u0[idx_in_A]\n                ui1 = u1[idx_in_A]\n                if abs(ui1) > eps:\n                    z_hit = -ui0 / ui1\n                    # only future breakpoints\n                    if z_hit > current_z + 1e-13 and z_hit <= z1 + 1e-13:\n                        candidate_z.append(z_hit)\n                        candidate_events.append((\"drop\", var_idx))\n                # if ui1 == 0 then coefficient constant; it won't hit zero except if ui0==0 (already zero)\n                # ignore ui1 near zero\n\n        # For each inactive variable j, solve |c0_j + z*c1_j| = lambda_val\n        # where c0_j = X_j^T r0 = Xt_a[j] - G[j,A] @ u0\n        # and   c1_j = X_j^T r1 = Xt_b[j] - G[j,A] @ u1\n        inactive_idx = [j for j in range(p) if j not in A]\n        if len(inactive_idx) > 0:\n            # compute G_jA @ u0 and G_jA @ u1 efficiently by accessing rows of G\n            # vectorized:\n            if k == 0:\n                c0 = Xt_a[inactive_idx].copy()\n                c1 = Xt_b[inactive_idx].copy()\n            else:\n                G_jA = G[np.ix_(inactive_idx, A)]   # shape (p-k, k)\n                c0 = Xt_a[inactive_idx] - (G_jA @ u0)\n                c1 = Xt_b[inactive_idx] - (G_jA @ u1)\n            # solve c0 + z*c1 = ±lambda_val -> z = (±lambda_val - c0) / c1 when c1 != 0\n            for ii, j in enumerate(inactive_idx):\n                ci0 = float(c0[ii])\n                ci1 = float(c1[ii])\n                if abs(ci1) < eps:\n                    # then c0 must be checked: if abs(ci0) == lambda already at current_z -> immediate\n                    # but we'll ignore near-constant ones unless they already saturate (rare)\n                    continue\n                # two candidates\n                z_pos = ( lambda_val - ci0) / ci1\n                z_neg = (-lambda_val - ci0) / ci1\n                # we only care about the smallest z > current_z\n                if z_pos > current_z + 1e-13 and z_pos <= z1 + 1e-13:\n                    candidate_z.append(z_pos)\n                    # sign of entering variable at z_pos:\n                    sign_pos = np.sign(ci0 + z_pos * ci1)\n                    candidate_events.append((\"add\", j, int(sign_pos)))\n                if z_neg > current_z + 1e-13 and z_neg <= z1 + 1e-13:\n                    candidate_z.append(z_neg)\n                    sign_neg = np.sign(ci0 + z_neg * ci1)\n                    candidate_events.append((\"add\", j, int(sign_neg)))\n\n        if len(candidate_z) == 0:\n            # no more breakpoints until z1\n            # append z1 and beta at z1, then break\n            z_final = float(z1)\n            # compute beta at z_final\n            if k == 0:\n                beta_full = np.zeros(p)\n            else:\n                betaA_final = u0 + z_final * u1\n                beta_full = np.zeros(p)\n                beta_full[A] = betaA_final\n            zk_list.append(z_final)\n            active_sets_list.append(list(A))\n            betas_list.append(beta_full)\n            current_z = z_final\n            if verbose:\n                print(f\"[iter {it}] No more candidates -> finish at z={z_final:.6g}\")\n            break\n\n        # find smallest candidate z (with tie tolerance)\n        candidate_z = np.array(candidate_z, dtype=float)\n        min_idx = np.argmin(candidate_z)\n        z_next = float(candidate_z[min_idx])\n\n        # handle numerical tiny negative jump\n        if z_next <= current_z + 1e-13:\n            # avoid infinite loop; nudge forward tiny bit\n            z_next = current_z + 1e-10\n\n        # Evaluate event(s) occurring at z_next (there may be multiple => process all)\n        # collect all events with |z - z_next| <= tol_rel\n        tol_rel = 1e-9 + 1e-12 * max(1.0, abs(z_next))\n        events_at_z = []\n        for (zcand, ev) in zip(candidate_z, candidate_events):\n            if abs(zcand - z_next) <= tol_rel:\n                events_at_z.append(ev)\n\n        # compute beta at z_next\n        if k == 0:\n            beta_full = np.zeros(p)\n        else:\n            betaA_next = u0 + z_next * u1\n            beta_full = np.zeros(p)\n            beta_full[A] = betaA_next\n\n        # update active set based on events: first handle drops (coeff hit zero), then adds\n        A_new = set(A)\n        s_new = {}\n        for ev in events_at_z:\n            if ev[0] == \"drop\":\n                _, var_idx = ev\n                if var_idx in A_new:\n                    A_new.remove(var_idx)\n                    if var_idx in s_new: del s_new[var_idx]\n            elif ev[0] == \"add\":\n                _, var_idx, sign_var = ev\n                # only add if not present\n                if var_idx not in A_new:\n                    A_new.add(var_idx)\n                    s_new[var_idx] = sign_var\n        # build s_vec for new A in ordered list\n        A_new = sorted(list(A_new))\n        if len(A_new) == 0:\n            s_vec_new = np.array([], dtype=float)\n        else:\n            # if sign is not set by event (e.g., from previous A), try to recover sign from beta_full if present\n            s_vec_new = []\n            for j in A_new:\n                if j in s_new:\n                    s_vec_new.append(s_new[j])\n                else:\n                    bj = beta_full[j]\n                    # if bj nearly zero -> fallback to KKT sign via correlation\n                    if abs(bj) > 1e-14:\n                        s_vec_new.append(np.sign(bj))\n                    else:\n                        # compute cor = X_j^T r(z_next) to decide sign\n                        cor_j = (Xt_a[j] + z_next * Xt_b[j]) - (G[j, A] @ (u0 + z_next * u1)) if k > 0 else (Xt_a[j] + z_next * Xt_b[j])\n                        s_vec_new.append(int(np.sign(cor_j) if cor_j != 0 else 1))\n            s_vec_new = np.array(s_vec_new, dtype=float)\n\n        # append results\n        zk_list.append(z_next)\n        active_sets_list.append(list(A_new))\n        betas_list.append(beta_full.copy())\n\n        if verbose:\n            print(f\"[iter {it}] z {current_z:.6g} -> {z_next:.6g} | |A| {k} -> {len(A_new)} | events: {events_at_z}\")\n\n        # update for next iteration\n        current_z = z_next\n        # prepare s_vec for next loop; align ordering to A_new\n        s_vec = s_vec_new\n        # if A changed, next iteration will recompute u0/u1 from new A\n        # continue loop\n\n    return zk_list, active_sets_list, betas_list\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_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        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\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\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=0.00000000001, max_iter=2000,\n#                     lasso_tol=1e-3, lasso_max_iter=20000):\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-8)\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\nimport traceback\nimport traceback\nimport os\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=1e-11, max_iter=2000,\n                    lasso_tol=1e-3, lasso_max_iter=20000):\n    \"\"\"\n    Worker: dùng lasso_path_tracker để lấy toàn bộ breakpoints trên đoạn [z_start, z_end]\n    Trả về list các tuple: (zk, A, bhz)\n    \"\"\"\n\n    # ép kiểu ngay để tránh lỗi format trong logger\n    try:\n        z_start = float(z_start)\n        z_end = float(z_end)\n    except Exception:\n        # cố gắng lấy phần tử đầu nếu truyền ndarray\n        z_start = float(np.asarray(z_start).ravel()[0])\n        z_end = float(np.asarray(z_end).ravel()[0])\n\n    # giới hạn threads trong worker để tránh oversubscription\n    os.environ.setdefault(\"OMP_NUM_THREADS\", \"1\")\n    os.environ.setdefault(\"MKL_NUM_THREADS\", \"1\")\n    os.environ.setdefault(\"OPENBLAS_NUM_THREADS\", \"1\")\n\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:.6f}, {z_end:.6f}]\")\n\n    start_time = time.time()\n\n    # đảm bảo biến results luôn tồn tại\n    results = []\n\n    # biến handles để đóng trong finally\n    shmX = shmy = shmetaj = None\n\n    try:\n        # Attach shared arrays (có thể ném lỗi nếu name sai)\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        # ensure vectors are 1D numpy arrays\n        y_vec = np.asarray(y).reshape(-1)\n        etaj_vec = np.asarray(etaj).reshape(-1)\n\n        # compute a, b (y(z) = a + b*z)\n        a_vec, b_vec = compute_a_b_from_etaj(y_vec, etaj_vec)\n\n        # call path tracker (assume lasso_path_tracker is in global scope)\n        zk_list, active_sets_list, betas_list = lasso_path_tracker(\n            X, a_vec, b_vec, lambda_val=alpha,\n            z0=float(z_start), z1=float(z_end),\n            max_iter=max_iter, tol=1e-12, verbose=False\n        )\n\n        # build results in same structure as before: (zk, A, bhz)\n        thresh = 1e-9\n        for i, zk in enumerate(zk_list):\n            bhz = np.asarray(betas_list[i]).copy()\n            bhz[np.abs(bhz) <= thresh] = 0.0\n            A = list(np.where(np.abs(bhz) > thresh)[0])\n            results.append((float(zk), A, bhz))\n\n        end_time = time.time()\n        duration = end_time - start_time\n        logger.info(f\"KẾT THÚC xử lý đoạn [{z_start:.6f}, {z_end:.6f}] - breakpoints: {len(zk_list)} - Time: {duration:.2f}s\")\n\n        return results\n\n    except Exception as e:\n        # log traceback để dễ debug và re-raise để joblib in lỗi\n        logger.error(f\"LỖI trong worker đoạn [{z_start:.6f}, {z_end:.6f}]: {e}\")\n        logger.error(traceback.format_exc())\n        # đóng shm nếu đã attach\n        try:\n            if shmX is not None: shmX.close()\n            if shmy is not None: shmy.close()\n            if shmetaj is not None: shmetaj.close()\n        except Exception:\n            pass\n        raise\n\n    finally:\n        # đảm bảo đóng shared memory handles nếu vẫn mở\n        try:\n            if shmX is not None:\n                shmX.close()\n            if shmy is not None:\n                shmy.close()\n            if shmetaj is not None:\n                shmetaj.close()\n        except Exception:\n            # không raise trong finally\n            pass\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=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    overlap=0.0001\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 = 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    \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    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    if denominator < 1e-100: return 0\n    \n    pivot_val = float(numerator / denominator)\n    print(\"Done7\")\n    return 2 * min(pivot_val, 1 - pivot_val)\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            \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    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    \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    scaler_y = StandardScaler()\n    y_scaled = scaler_y.fit_transform(y_clean.values.reshape(-1, 1))\n    \n    n, p = X_scaled.shape\n    BEST_ALPHA = 9000\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-5).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    # 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    # --- 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        if (i==6):\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-17T06:58:43.196023Z","iopub.execute_input":"2025-08-17T06:58:43.196456Z"}},"outputs":[],"execution_count":null}]}