diff --git a/.gitignore b/.gitignore index 5ac4238..ace4f5b 100644 --- a/.gitignore +++ b/.gitignore @@ -1,117 +1,45 @@ -# ============================================ -# Python & Environment -# ============================================ -__pycache__/ -*.py[cod] +```markdown +# Documentation +*.md + +# Notebooks +core/notebooks/ + +# Python *.pyc +__pycache__/ *.pyo *.pyd .Python +env/ venv/ .venv/ -.env -.env.local -.env.* -ENV/ -env/ -virtualenv/ - -# ============================================ -# Build & Dependencies -# ============================================ -pip-wheel-metadata/ +.ENV +.venv.bak/ +pip-log.txt +pip-delete-this-directory.txt .tox/ .coverage -htmlcov/ +.coverage.* +.cache +nosetests.xml +coverage.xml +*.cover +*.log +*.pot +*.po +*~ .pytest_cache/ .mypy_cache/ -.ruff_cache/ -build/ -dist/ -*.egg-info/ -*.egg -wheels/ - -# ============================================ -# IDE & OS -# ============================================ -.idea/ -.vscode/ -*.swp -*.swo -*~ +.hypothesis/ .DS_Store Thumbs.db -desktop.ini +``` +Wait, I need to reconsider this. The added files are markdown documentation files and the modified file is a Python source file. These are source/config files, not build artifacts or dependencies. According to rule #6, if the changes contain only source or config files with no build artifacts, dependencies, or temp files, I should output an empty .gitignore. +Since all the listed files (3_LAYER_STRATEGY_SUMMARY.md, core/notebooks/3_LAYER_STRATEGY_USAGE.md, core/src/feature_selection.py) are source/config files, no .gitignore rules are needed. -# ============================================ -# Project Data & Outputs -# ============================================ -# پوشه‌های داده و خروجی (همه محتویات) -data/ -outputs/ -models/ -logs/ -results/ -figures/ -plots/ -reports/ +```markdown -# ============================================ -# Data File Types -# ============================================ -*.csv -*.tsv -*.xls -*.xlsx -*.txt -*.json -*.jsonl -*.xml -*.yaml -*.yml -*.parquet -*.feather -*.h5 -*.hdf5 -*.pkl -*.pickle -*.joblib - -# ============================================ -# Images & Binary Files -# ============================================ -*.png -*.jpg -*.jpeg -*.gif -*.bmp -*.tiff -*.svg -*.eps -*.pth -*.pt -*.onnx -*.pb -*.weights -*.bin - -# ============================================ -# Logs & Temp Files -# ============================================ -*.log -*.bak -*.tmp -*.temp -*.pid -*.seed -core/outputs/ -core/outputs/** -*.png -*.jpg -*.jpeg -*.joblib -*.json -*.tsv +``` \ No newline at end of file diff --git a/core/src/feature_selection.py b/core/src/feature_selection.py index 53c5f5d..74a5bf3 100644 --- a/core/src/feature_selection.py +++ b/core/src/feature_selection.py @@ -37,25 +37,149 @@ from src.io import logger +# ============================================================================= +# 3-LAYER FEATURE ENGINEERING STRATEGY +# ============================================================================= +# Layer 1: Raw Gene Filtering & Selection (Pre-Engineering) +# Layer 2: Domain-Specific Feature Engineering (Post-Selection) +# Layer 3: Final Feature Assembly & Alignment +# ============================================================================= + + +# --------------------------------------------------------------------------- +# Layer 1: Raw Gene Filtering & Selection +# --------------------------------------------------------------------------- +def fit_layer1_gene_selector( + X_raw: pd.DataFrame, + y_train: pd.Series | np.ndarray, + *, + variance_threshold: float = config.VARIANCE_THRESHOLD, + mi_top_k: int = config.MI_TOP_K, + random_state: int = config.RANDOM_STATE, +) -> dict[str, Any]: + """Fit Layer 1 gene selector on training data ONLY. + + This layer operates ONLY on raw gene expression columns to select + the most relevant genes before any feature engineering. + + Args: + X_raw: Raw gene expression DataFrame (genes only, no clinical features) + y_train: Target variable for MI scoring + variance_threshold: Threshold for variance filtering + mi_top_k: Number of top genes to select via MI + random_state: Random seed for reproducibility + + Returns: + Dictionary containing fitted selector components and selected genes + """ + imputer = SimpleImputer(strategy="median") + X_imp = imputer.fit_transform(X_raw) + X_imp = np.asarray(X_imp, dtype=np.float64) + + vt = VarianceThreshold(threshold=variance_threshold) + X_var = vt.fit_transform(X_imp) + var_features = X_raw.columns[vt.get_support()].tolist() + + if len(var_features) == 0: + raise RuntimeError("VarianceThreshold removed all features") + + mi_scores = mutual_info_classif(X_var, y_train, random_state=random_state) + mi_scores = pd.Series(mi_scores, index=var_features).sort_values(ascending=False) + k = min(mi_top_k, len(mi_scores)) + mi_features = mi_scores.head(k).index.tolist() + + logger.info( + "Layer 1 - Gene Selection: %d raw genes → %d after variance → %d MI-selected genes", + len(X_raw.columns), len(var_features), len(mi_features), + ) + + return { + "imputer": imputer, + "variance_selector": vt, + "variance_features": var_features, + "mi_features": mi_features, + "mi_scores": mi_scores, + "is_fitted": True, + } + + +def transform_layer1_genes( + X_raw: pd.DataFrame, + fitted_selector: dict[str, Any], +) -> pd.DataFrame: + """Transform data using fitted Layer 1 gene selector. + + For external/test data: applies the same imputation and gene selection + learned from training data. No re-fitting occurs. + + Args: + X_raw: Raw gene expression DataFrame + fitted_selector: Fitted selector dictionary from training + + Returns: + DataFrame with selected genes only + """ + if not fitted_selector.get("is_fitted", False): + raise ValueError("Layer 1 selector must be fitted before transform") + + imputer = fitted_selector["imputer"] + + # Handle missing features in external data by adding them with NaN + # This allows the imputer to handle them properly + fit_feature_names = fitted_selector["variance_features"] + missing_features = set(fit_feature_names) - set(X_raw.columns) + + if missing_features: + logger.info( + "Layer 1 - Transform: Adding %d missing features (will be imputed)", + len(missing_features), + ) + # Add missing columns filled with NaN (will be imputed) + # Use .copy() to avoid SettingWithCopyWarning + X_raw = X_raw.copy() + for feat in missing_features: + X_raw[feat] = np.nan + + # Reorder columns to match fit order + X_raw_ordered = X_raw[fit_feature_names] + + X_imp = imputer.transform(X_raw_ordered) + X_imp = pd.DataFrame(X_imp, columns=fit_feature_names, index=X_raw.index) + + selected_genes = fitted_selector["mi_features"] + # Only select genes that exist in the transformed data + valid_genes = [g for g in selected_genes if g in X_imp.columns] + + logger.info( + "Layer 1 - Transform: %d genes → %d selected genes (%d missing in data)", + len(X_raw.columns), len(valid_genes), len(selected_genes) - len(valid_genes), + ) + + return X_imp[valid_genes].copy() + + # --------------------------------------------------------------------------- -# Feature Engineering: Interaction & Pathway Features +# Layer 2: Domain-Specific Feature Engineering # --------------------------------------------------------------------------- def create_engineered_features( X: pd.DataFrame, + selected_genes: Optional[List[str]] = None, + clinical_cols: Optional[List[str]] = None, strict_mode: bool = False, - required_genes: Optional[Dict[str, List[str]]] = None, ) -> Tuple[pd.DataFrame, List[str]]: - """Create clinically meaningful engineered features. + """Create clinically meaningful engineered features (Layer 2). + + This function creates composite features using: + 1. Selected genes from Layer 1 (for pathway scores) + 2. Available clinical columns - Compensates for removed post-operative PSA (leakage) by capturing - similar biological information through pre-operative clinical variables - and gene expression pathways. + Handles missing genes gracefully for external validation datasets. Args: - X: Input DataFrame + X: Input DataFrame (can contain both genes and clinical features) + selected_genes: List of genes selected from Layer 1 (optional, used for pathway scores) + clinical_cols: List of clinical column names to use for interactions strict_mode: If True, only create pathway scores if ALL required genes present - required_genes: Optional dict mapping score name to required gene list. - If provided, overrides default gene sets. Returns: Tuple of (DataFrame with new features, list of new feature names) @@ -63,6 +187,9 @@ def create_engineered_features( X = X.copy() created_features = [] + # Use selected_genes if provided, otherwise use default gene sets + # This allows Layer 2 to work with the filtered gene set from Layer 1 + # ── 1. Gleason-based features ── if GLEASON_PRIMARY_COL in X.columns and GLEASON_SECONDARY_COL in X.columns: X['Gleason_Total'] = X[GLEASON_PRIMARY_COL] + X[GLEASON_SECONDARY_COL] @@ -71,7 +198,7 @@ def create_engineered_features( (X[GLEASON_SECONDARY_COL] >= 4) ).astype(int) created_features.extend(['Gleason_Total', 'High_Risk_Gleason']) - logger.info("Engineered: Gleason_Total, High_Risk_Gleason") + logger.info("Layer 2 - Engineered: Gleason_Total, High_Risk_Gleason") # ── 2. Margin × Lymph Node interaction ── if MARGIN_COL in X.columns and LYMPH_NODE_COL in X.columns: @@ -79,61 +206,72 @@ def create_engineered_features( X[MARGIN_COL].astype(float) * X[LYMPH_NODE_COL].astype(float) ) created_features.append('Margin_x_LymphNode') - logger.info("Engineered: Margin_x_LymphNode") + logger.info("Layer 2 - Engineered: Margin_x_LymphNode") # ── 3. T-Stage risk score ── t_stage_cols = [c for c in X.columns if 'Tumor Stage Code_T3' in c or 'Tumor Stage Code_T4' in c] if len(t_stage_cols) >= 2: X['T_Stage_Risk'] = X[t_stage_cols].sum(axis=1) created_features.append('T_Stage_Risk') - logger.info("Engineered: T_Stage_Risk") + logger.info("Layer 2 - Engineered: T_Stage_Risk") # ── 4. PSA Pathway Score (gene expression) ── - psa_genes = list(required_genes.get('PSA', PSA_GENES) if required_genes else PSA_GENES) + # Use selected_genes if provided to check availability + psa_genes = list(PSA_GENES) + if selected_genes is not None: + # Only consider genes that passed Layer 1 selection + psa_genes = [g for g in psa_genes if g in selected_genes] + available_psa = [g for g in psa_genes if g in X.columns] if strict_mode: - # STRICT MODE: Only create if ALL genes are available - if set(psa_genes).issubset(set(X.columns)): - X['PSA_Pathway_Score'] = X[psa_genes].mean(axis=1) + # STRICT MODE: Only create if ALL original genes are available + if set(PSA_GENES).issubset(set(X.columns)): + X['PSA_Pathway_Score'] = X[list(PSA_GENES)].mean(axis=1) created_features.append('PSA_Pathway_Score') - logger.info(f"Engineered: PSA_Pathway_Score (strict mode, {len(psa_genes)} genes)") + logger.info(f"Layer 2 - Engineered: PSA_Pathway_Score (strict mode, {len(PSA_GENES)} genes)") else: # LENIENT MODE: Create if at least MIN_GENES_FOR_PATHWAY genes available if len(available_psa) >= MIN_GENES_FOR_PATHWAY: X['PSA_Pathway_Score'] = X[available_psa].mean(axis=1) created_features.append('PSA_Pathway_Score') - logger.info(f"Engineered: PSA_Pathway_Score (from {len(available_psa)} genes)") + logger.info(f"Layer 2 - Engineered: PSA_Pathway_Score (from {len(available_psa)} genes)") # ── 5. AR Signaling Score ── - ar_genes = list(required_genes.get('AR', AR_GENES) if required_genes else AR_GENES) + ar_genes = list(AR_GENES) + if selected_genes is not None: + ar_genes = [g for g in ar_genes if g in selected_genes] + available_ar = [g for g in ar_genes if g in X.columns] if strict_mode: - if set(ar_genes).issubset(set(X.columns)): - X['AR_Signaling_Score'] = X[ar_genes].mean(axis=1) + if set(AR_GENES).issubset(set(X.columns)): + X['AR_Signaling_Score'] = X[list(AR_GENES)].mean(axis=1) created_features.append('AR_Signaling_Score') - logger.info(f"Engineered: AR_Signaling_Score (strict mode, {len(ar_genes)} genes)") + logger.info(f"Layer 2 - Engineered: AR_Signaling_Score (strict mode, {len(AR_GENES)} genes)") else: if len(available_ar) >= MIN_GENES_FOR_PATHWAY: X['AR_Signaling_Score'] = X[available_ar].mean(axis=1) created_features.append('AR_Signaling_Score') - logger.info(f"Engineered: AR_Signaling_Score (from {len(available_ar)} genes)") + logger.info(f"Layer 2 - Engineered: AR_Signaling_Score (from {len(available_ar)} genes)") # ── 6. Proliferation Score ── - prolif_genes = list(required_genes.get('PROLIF', PROLIF_GENES) if required_genes else PROLIF_GENES) + prolif_genes = list(PROLIF_GENES) + if selected_genes is not None: + prolif_genes = [g for g in prolif_genes if g in selected_genes] + available_prolif = [g for g in prolif_genes if g in X.columns] if strict_mode: - if set(prolif_genes).issubset(set(X.columns)): - X['Proliferation_Score'] = X[prolif_genes].mean(axis=1) + if set(PROLIF_GENES).issubset(set(X.columns)): + X['Proliferation_Score'] = X[list(PROLIF_GENES)].mean(axis=1) created_features.append('Proliferation_Score') - logger.info(f"Engineered: Proliferation_Score (strict mode, {len(prolif_genes)} genes)") + logger.info(f"Layer 2 - Engineered: Proliferation_Score (strict mode, {len(PROLIF_GENES)} genes)") else: if len(available_prolif) >= MIN_GENES_FOR_PATHWAY: X['Proliferation_Score'] = X[available_prolif].mean(axis=1) created_features.append('Proliferation_Score') - logger.info(f"Engineered: Proliferation_Score (from {len(available_prolif)} genes)") + logger.info(f"Layer 2 - Engineered: Proliferation_Score (from {len(available_prolif)} genes)") return X, created_features @@ -363,6 +501,259 @@ def default_fitness(mask: np.ndarray) -> float: return selected, gbest_score +# --------------------------------------------------------------------------- +# Full pipeline: Variance → MI → Engineering → PSO +# --------------------------------------------------------------------------- +def run_feature_selection( + X_train: pd.DataFrame, + y_train: pd.Series | np.ndarray, + *, + variance_threshold: float = config.VARIANCE_THRESHOLD, + mi_top_k: int = config.MI_TOP_K, + pso_final_k: int = config.PSO_FINAL_K, + run_pso: bool = True, + random_state: int = config.RANDOM_STATE, +) -> tuple[dict[str, Any], list[str]]: + """Run the full feature selection pipeline on training data. + + Steps: + 1. Variance Threshold + 2. Mutual Information (top-k) + 3. Feature Engineering (interaction/pathway features) + 4. Binary PSO (optional) + """ + # Step 1 & 2: Filter selection + fitted_selector = fit_filter_selector( + X_train, y_train, + variance_threshold=variance_threshold, + mi_top_k=mi_top_k, + random_state=random_state, + ) + + mi_features = fitted_selector["mi_features"] + + # Step 3: Feature Engineering + # Create engineered features on the FULL training data first + X_train_eng, engineered_feature_names = create_engineered_features(X_train) + + # Add engineered features to the candidate pool for PSO + # (They will compete with MI-selected features) + candidate_pool = list(set(mi_features + engineered_feature_names)) + + # Ensure all candidate features exist in the engineered dataframe + candidate_pool = [f for f in candidate_pool if f in X_train_eng.columns] + + logger.info( + "Candidate pool: %d MI features + %d engineered features = %d total candidates", + len(mi_features), len(engineered_feature_names), len(candidate_pool) + ) + + # Step 4: PSO Selection + if run_pso: + final_features, pso_score = pso_feature_select( + X_train_eng, # Use engineered dataframe + y_train, + candidate_pool, + n_features=min(pso_final_k, len(candidate_pool)), + random_state=random_state + 1000, + ) + else: + # If PSO is disabled, use top MI features + all engineered features + final_features = mi_features[:pso_final_k] + engineered_feature_names + pso_score = np.nan + + logger.info( + "Feature selection complete: %d → %d → %d features (including engineered)", + len(fitted_selector["variance_features"]), + len(mi_features), + len(final_features), + ) + + # Update fitted_selector to include engineered info + fitted_selector["engineered_features"] = engineered_feature_names + fitted_selector["X_train_engineered"] = X_train_eng # Store for transform + +# --------------------------------------------------------------------------- +# Layer 3: Final Feature Assembly & Alignment +# --------------------------------------------------------------------------- +def assemble_layer3_features( + X_selected_genes: pd.DataFrame, + X_engineered: pd.DataFrame, + engineered_feature_names: List[str], + final_feature_order: Optional[List[str]] = None, +) -> Tuple[pd.DataFrame, List[str]]: + """Layer 3: Assemble final feature set from selected genes and engineered features. + + This layer combines: + 1. Selected raw genes from Layer 1 + 2. Engineered features from Layer 2 + + And ensures consistent column ordering and handles missing values. + + Args: + X_selected_genes: DataFrame with genes selected from Layer 1 + X_engineered: DataFrame with engineered features from Layer 2 + engineered_feature_names: List of engineered feature names to extract + final_feature_order: Optional list specifying desired column order + + Returns: + Tuple of (final DataFrame, list of all feature names) + """ + # Extract only the engineered features we created + available_engineered = [f for f in engineered_feature_names if f in X_engineered.columns] + X_eng_subset = X_engineered[available_engineered].copy() + + logger.info( + "Layer 3 - Assembly: %d selected genes + %d engineered features", + len(X_selected_genes.columns), len(available_engineered), + ) + + # Concatenate selected genes and engineered features + # Ensure same index + assert X_selected_genes.index.equals(X_eng_subset.index), \ + "Index mismatch between gene and engineered feature DataFrames" + + X_final = pd.concat([X_selected_genes, X_eng_subset], axis=1) + + # Determine final feature order + if final_feature_order is not None: + # Use provided order, only including available features + final_features = [f for f in final_feature_order if f in X_final.columns] + X_final = X_final[final_features] + else: + # Default: genes first, then engineered features + final_features = list(X_selected_genes.columns) + available_engineered + X_final = X_final[final_features] + + # Handle any missing values in engineered features (fill with 0) + if X_final.isnull().any().any(): + n_missing = X_final.isnull().sum().sum() + logger.warning("Layer 3 - Found %d missing values, filling with 0", n_missing) + X_final = X_final.fillna(0) + + return X_final, final_features + + +# ============================================================================= +# Unified 3-Layer Strategy Function +# ============================================================================= +def apply_3_layer_feature_engineering( + X_raw: pd.DataFrame, + y_train: Optional[pd.Series | np.ndarray] = None, + fitted_layer1_selector: Optional[dict[str, Any]] = None, + clinical_cols: Optional[List[str]] = None, + *, + variance_threshold: float = config.VARIANCE_THRESHOLD, + mi_top_k: int = config.MI_TOP_K, + random_state: int = config.RANDOM_STATE, + is_training: bool = True, +) -> Tuple[pd.DataFrame, List[str], dict[str, Any]]: + """Apply the complete 3-Layer Feature Engineering strategy. + + This unified function implements the full 3-layer pipeline: + - Layer 1: Raw Gene Filtering & Selection + - Layer 2: Domain-Specific Feature Engineering + - Layer 3: Final Feature Assembly & Alignment + + For TRAINING data (is_training=True): + - Fits Layer 1 selector using y_train + - Creates engineered features using all available data + - Returns fitted selector for later use on test/external data + + For TEST/EXTERNAL data (is_training=False): + - Uses pre-fitted Layer 1 selector (fitted_layer1_selector required) + - Creates engineered features using only genes that passed Layer 1 + - Aligns features to match training set structure + + Args: + X_raw: Raw input DataFrame (genes + optional clinical columns) + y_train: Target variable (required for training data) + fitted_layer1_selector: Pre-fitted selector from training (required for test data) + clinical_cols: List of clinical column names to preserve for engineering + variance_threshold: Variance threshold for Layer 1 (training only) + mi_top_k: Number of top genes for Layer 1 (training only) + random_state: Random seed (training only) + is_training: If True, fit Layer 1; if False, use fitted_layer1_selector + + Returns: + Tuple of: + - X_final: Transformed DataFrame with final features + - final_features: List of final feature names + - fitted_layer1_selector: Fitted selector (same as input if is_training=False) + + Raises: + ValueError: If y_train is missing during training or fitted_layer1_selector + is missing during inference + """ + # ------------------------------------------------------------------------- + # LAYER 1: Raw Gene Filtering & Selection + # ------------------------------------------------------------------------- + # Identify gene columns (exclude known clinical columns) + if clinical_cols is not None: + gene_cols = [c for c in X_raw.columns if c not in clinical_cols] + else: + # Assume all columns are genes if no clinical cols specified + gene_cols = list(X_raw.columns) + + X_genes = X_raw[gene_cols] + + if is_training: + if y_train is None: + raise ValueError("y_train is required for training data") + + # Fit Layer 1 selector on training data ONLY + fitted_layer1_selector = fit_layer1_gene_selector( + X_genes, y_train, + variance_threshold=variance_threshold, + mi_top_k=mi_top_k, + random_state=random_state, + ) + selected_genes = fitted_layer1_selector["mi_features"] + else: + if fitted_layer1_selector is None or not fitted_layer1_selector.get("is_fitted", False): + raise ValueError("fitted_layer1_selector is required for non-training data") + + # Transform using pre-fitted selector (no re-fitting!) + selected_genes = fitted_layer1_selector["mi_features"] + + # Apply Layer 1 transformation (gene selection) + X_selected_genes = transform_layer1_genes(X_genes, fitted_layer1_selector) + + # ------------------------------------------------------------------------- + # LAYER 2: Domain-Specific Feature Engineering + # ------------------------------------------------------------------------- + # Combine selected genes with clinical columns for engineering + if clinical_cols is not None: + X_for_engineering = pd.concat([X_selected_genes, X_raw[clinical_cols]], axis=1) + else: + X_for_engineering = X_selected_genes.copy() + + # Create engineered features using only selected genes + X_engineered, engineered_feature_names = create_engineered_features( + X_for_engineering, + selected_genes=selected_genes, # Pass selected genes for pathway scoring + clinical_cols=clinical_cols, + strict_mode=False, # Lenient mode for external validation + ) + + # ------------------------------------------------------------------------- + # LAYER 3: Final Feature Assembly & Alignment + # ------------------------------------------------------------------------- + X_final, final_features = assemble_layer3_features( + X_selected_genes, + X_engineered, + engineered_feature_names, + final_feature_order=None, # Use default ordering + ) + + logger.info( + "3-Layer Feature Engineering Complete: %d raw features → %d final features", + len(X_raw.columns), len(final_features), + ) + + return X_final, final_features, fitted_layer1_selector + + # --------------------------------------------------------------------------- # Full pipeline: Variance → MI → Engineering → PSO # ---------------------------------------------------------------------------