{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.12.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":112509,"databundleVersionId":14254895,"sourceType":"competition"}],"dockerImageVersionId":31234,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"\"\"\"\nCSIRO Image2Biomass - Tools for EDA, data quality checks, and model validation\n\"\"\"\n\nimport pandas as pd\nimport numpy as np\nfrom pathlib import Path\nimport warnings\nwarnings.filterwarnings('ignore')\n\n# Core ML & Stats\nfrom sklearn.ensemble import IsolationForest, RandomForestRegressor, RandomForestClassifier\nfrom sklearn.neighbors import LocalOutlierFactor\nfrom sklearn.cluster import KMeans, DBSCAN\nfrom sklearn.preprocessing import StandardScaler, LabelEncoder\nfrom sklearn.decomposition import PCA\nfrom sklearn.metrics import silhouette_score, davies_bouldin_score\nfrom scipy import stats\n\n# Visualization\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport plotly.express as px\nimport plotly.graph_objects as go\nfrom plotly.subplots import make_subplots\n\nPLOTLY_AVAILABLE = True","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-12-28T15:31:17.415540Z","iopub.execute_input":"2025-12-28T15:31:17.415811Z","iopub.status.idle":"2025-12-28T15:31:17.422098Z","shell.execute_reply.started":"2025-12-28T15:31:17.415793Z","shell.execute_reply":"2025-12-28T15:31:17.421050Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def load_and_pivot_biomass_data(csv_path):\n    \"\"\"\n    Load CSIRO train.csv and convert from long to wide format.\n    \n    Args:\n        csv_path: Path to train.csv\n        \n    Returns:\n        DataFrame in wide format (1 row per image)\n    \"\"\"\n    print(\"\\nLoading CSIRO Biomass Dataset...\")\n    print(\"=\" * 60)\n    \n    # Load data\n    df_long = pd.read_csv(csv_path)\n    print(f\"Loaded: {len(df_long)} rows × {len(df_long.columns)} columns (long format)\")\n    \n    # Extract image ID from sample_id\n    df_long['image_id'] = df_long['sample_id'].str.split('__').str[0]\n    \n    # Pivot targets from long to wide\n    df_wide = df_long.pivot_table(\n        index=['image_id', 'image_path', 'Sampling_Date', 'State', 'Species', \n               'Pre_GSHH_NDVI', 'Height_Ave_cm'],\n        columns='target_name',\n        values='target',\n        aggfunc='first'\n    ).reset_index()\n    \n    print(f\"Pivoted: {len(df_wide)} images × {len(df_wide.columns)} columns (wide format)\")\n    \n    # Feature Engineering: Date features\n    print(\"\\nEngineering temporal features...\")\n    df_wide['Sampling_Date'] = pd.to_datetime(df_wide['Sampling_Date'], format='%m/%d/%Y', errors='coerce')\n    df_wide['Month'] = df_wide['Sampling_Date'].dt.month\n    df_wide['DayOfYear'] = df_wide['Sampling_Date'].dt.dayofyear\n    df_wide['Season'] = pd.cut(df_wide['Month'], \n                                bins=[0, 3, 6, 9, 12], \n                                labels=['Summer', 'Autumn', 'Winter', 'Spring'])\n    \n    # Cyclical encoding for seasonality\n    df_wide['Month_sin'] = np.sin(2 * np.pi * df_wide['Month'] / 12)\n    df_wide['Month_cos'] = np.cos(2 * np.pi * df_wide['Month'] / 12)\n    \n    # Feature Engineering: Species complexity\n    df_wide['Species_Count'] = df_wide['Species'].str.count('_') + 1\n    df_wide['Has_Clover'] = df_wide['Species'].str.contains('Clover|clover', na=False).astype(int)\n    df_wide['Has_Ryegrass'] = df_wide['Species'].str.contains('Ryegrass|ryegrass', na=False).astype(int)\n    \n    # Encode State as numeric\n    le_state = LabelEncoder()\n    df_wide['State_Encoded'] = le_state.fit_transform(df_wide['State'])\n    \n    # Derived biomass ratios (as mentioned in research)\n    df_wide['Dead_Ratio'] = df_wide['Dry_Dead_g'] / (df_wide['Dry_Total_g'] + 1e-6)\n    df_wide['Green_Ratio'] = df_wide['Dry_Green_g'] / (df_wide['Dry_Total_g'] + 1e-6)\n    df_wide['Clover_Ratio'] = df_wide['Dry_Clover_g'] / (df_wide['Dry_Total_g'] + 1e-6)\n    \n    # NDVI-Height interaction (mentioned in research as strong predictor)\n    df_wide['NDVI_Height_Interaction'] = df_wide['Pre_GSHH_NDVI'] * df_wide['Height_Ave_cm']\n    \n    # Biomass density\n    df_wide['Biomass_per_Height'] = df_wide['Dry_Total_g'] / (df_wide['Height_Ave_cm'] + 1e-6)\n    \n    print(f\"Final dataset: {len(df_wide)} samples x {len(df_wide.columns)} features\")\n    \n    # Data quality report\n    print(\"\\nData Quality Summary:\")\n    print(f\"  - Missing values: {df_wide.isnull().sum().sum()}\")\n    print(f\"  - States: {df_wide['State'].unique().tolist()}\")\n    print(f\"  - Date range: {df_wide['Sampling_Date'].min()} to {df_wide['Sampling_Date'].max()}\")\n    \n    # Target statistics\n    target_cols = ['Dry_Green_g', 'Dry_Dead_g', 'Dry_Clover_g', 'GDM_g', 'Dry_Total_g']\n    print(\"\\nTarget Statistics:\")\n    for col in target_cols:\n        zero_pct = (df_wide[col] == 0).sum() / len(df_wide) * 100\n        print(f\"  - {col:15s}: mean={df_wide[col].mean():6.2f}g, \"\n              f\"zeros={zero_pct:5.1f}%, CV={df_wide[col].std() / (df_wide[col].mean() + 1e-6):.2f}\")\n    \n    return df_wide","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-28T15:31:17.429617Z","iopub.execute_input":"2025-12-28T15:31:17.429923Z","iopub.status.idle":"2025-12-28T15:31:17.451035Z","shell.execute_reply.started":"2025-12-28T15:31:17.429899Z","shell.execute_reply":"2025-12-28T15:31:17.449856Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"class AnomalyDetector:\n    \"\"\"\n    Multi-algorithm anomaly detection with ensemble scoring.\n    Combines Isolation Forest, LOF, and statistical methods.\n    \"\"\"\n\n    def __init__(self, contamination=0.05, n_neighbors=20):\n        self.contamination = contamination\n        self.n_neighbors = n_neighbors\n        self.models = {}\n        self.results = {}\n\n    def detect_anomalies(self, df, numeric_cols):\n        \"\"\"\n        Run multiple anomaly detection algorithms and ensemble results.\n        \"\"\"\n        print(\"\\nANOMALY DETECTION\")\n        print(\"=\" * 60)\n\n        X = df[numeric_cols].copy()\n        X = X.fillna(X.median())\n\n        # Standardize for distance-based methods\n        scaler = StandardScaler()\n        X_scaled = scaler.fit_transform(X)\n\n        # Handle Zero Variance columns\n        # StandardScaler creates NaNs if a column has 0 variance (division by zero) => replace these NaNs with 0\n        X_scaled = np.nan_to_num(X_scaled)\n\n        # 1. Isolation Forest (tree-based)\n        print(\"Running Isolation Forest...\")\n        iso_forest = IsolationForest(\n            contamination=self.contamination, \n            random_state=42, \n            n_estimators=200\n        )\n        iso_scores = iso_forest.fit_predict(X_scaled)\n        iso_anomaly_scores = iso_forest.score_samples(X_scaled)\n        \n        # 2. Local Outlier Factor (density-based)\n        print(\"Running Local Outlier Factor...\")\n        lof = LocalOutlierFactor(\n            n_neighbors=self.n_neighbors, \n            contamination=self.contamination\n        )\n        lof_scores = lof.fit_predict(X_scaled)\n        lof_anomaly_scores = lof.negative_outlier_factor_\n\n        # 3. Statistical Z-score method (multivariate)\n        print(\"Running Statistical Z-score...\")\n        z_scores = np.abs(stats.zscore(X_scaled, axis=0))\n        z_anomalies = (z_scores > 3).any(axis=1).astype(int)\n        z_anomalies = np.where(z_anomalies == 1, -1, 1)\n\n        # -----------\n        # Ensemble: Majority voting\n        ensemble_votes = (iso_scores + lof_scores + z_anomalies) / 3\n        is_anomaly = ensemble_votes < 0\n\n        # Composite anomaly score (normalized)\n        composite_score = (\n            -iso_anomaly_scores / np.abs(iso_anomaly_scores).max() +\n            -lof_anomaly_scores / np.abs(lof_anomaly_scores).max() +\n            z_scores.max(axis=1) / 3\n        ) / 3\n\n        self.results = {\n            'is_anomaly': is_anomaly,\n            'anomaly_score': composite_score,\n            'isolation_forest': iso_scores,\n            'lof': lof_scores,\n            'z_score': z_anomalies\n        }\n        \n        # Report findings\n        n_anomalies = is_anomaly.sum()\n        print(f\"\\nDetected {n_anomalies} anomalies ({n_anomalies/len(df)*100:.2f}%)\")\n        \n        if n_anomalies > 0:\n            anomaly_indices = np.where(is_anomaly)[0]\n            print(f\"\\nTop 5 anomaly indices: {anomaly_indices[:5].tolist()}\")\n            \n            # Show which states have anomalies\n            if 'State' in df.columns:\n                anomaly_states = df.iloc[anomaly_indices]['State'].value_counts()\n                print(f\"\\nAnomalies by State:\")\n                for state, count in anomaly_states.items():\n                    print(f\"  {state}: {count}\")\n        \n        return self.results","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-28T15:31:17.452521Z","iopub.execute_input":"2025-12-28T15:31:17.452820Z","iopub.status.idle":"2025-12-28T15:31:17.484881Z","shell.execute_reply.started":"2025-12-28T15:31:17.452801Z","shell.execute_reply":"2025-12-28T15:31:17.483574Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"class MultivariateAnalyzer:\n    \"\"\"\n    Discovers interactions, correlations, and generates engineered features.\n    \"\"\"\n    def __init__(self, correlation_threshold=0.7, interaction_threshold=0.3):\n        self.correlation_threshold = correlation_threshold\n        self.interaction_threshold = interaction_threshold\n        self.interaction_features = []\n\n    def analyze(self, df, numeric_cols, target_col=None):\n        \"\"\"\n        Comprehensive multivariate analysis with automatic feature engineering.\n        \"\"\"\n        print(\"\\nMULTIVARIATE INTERACTION ANALYSIS\")\n        print(\"=\" * 60)\n\n        X = df[numeric_cols].copy()\n\n        # 1. Correlation Analysis\n        print(\"Computing correlation matrix...\")\n        corr_matrix = X.corr()\n\n        # Find strong correlations\n        strong_corr = []\n        for i in range(len(corr_matrix.columns)):\n            for j in range(i+1, len(corr_matrix.columns)):\n                corr_val = corr_matrix.iloc[i, j]\n                if abs(corr_val) > self.correlation_threshold:\n                    strong_corr.append((corr_matrix.columns[i], corr_matrix.columns[j], corr_val))\n        \n        if strong_corr:\n            print(f\"\\nFound {len(strong_corr)} strong correlations (|r| > {self.correlation_threshold}):\")\n            for col1, col2, corr in sorted(strong_corr, key=lambda x: abs(x[2]), reverse=True)[:15]:\n                print(f\"  {col1:25s} <-> {col2:25s}: r={corr:6.3f}\")\n\n        # 2. Automatic Interaction Feature Generation\n        print(\"\\nGenerating interaction features...\")\n        new_features = {}\n\n        # Generate multiplicative and ratio features for correlated pairs\n        for col1, col2, corr in strong_corr[:20]:\n            # Multiplicative interaction\n            interaction_name = f\"{col1}_x_{col2}\"\n            new_features[interaction_name] = df[col1] * df[col2]\n\n            # Ratio interaction (if no zeros)\n            if (df[col2] != 0).all():\n                ratio_name = f\"{col1}_div_{col2}\"\n                new_features[ratio_name] = df[col1] / df[col2]\n\n        # 3. Polynomial features for highly skewed columns\n        for col in numeric_cols[:10]:\n            if df[col].std() > 0:\n                skew = df[col].skew()\n                if abs(skew) > 1:\n                    new_features[f\"{col}_squared\"] = df[col] ** 2\n                    new_features[f\"{col}_log\"] = np.log1p(df[col].clip(lower=0))\n\n        self.interaction_features = list(new_features.keys())\n        print(f\"Generated {len(new_features)} interaction features\")\n\n        # 4. If target provided, rank interactions by correlation with target\n        if target_col and target_col in df.columns:\n            print(f\"\\nRanking features by correlation with '{target_col}':\")\n            feature_df = pd.DataFrame(new_features)\n            feature_corr = feature_df.corrwith(df[target_col]).abs().sort_values(ascending=False)\n            \n            print(\"Top 10 interaction features:\")\n            for feat, corr in feature_corr.head(10).items():\n                print(f\"  {feat:40s}: r={corr:.3f}\")\n        \n        return new_features, corr_matrix","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-28T15:31:17.486135Z","iopub.execute_input":"2025-12-28T15:31:17.486395Z","iopub.status.idle":"2025-12-28T15:31:17.506095Z","shell.execute_reply.started":"2025-12-28T15:31:17.486376Z","shell.execute_reply":"2025-12-28T15:31:17.504982Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"class ClusteringEngine:\n    \"\"\"\n    Unsupervised pattern discovery with automatic cluster optimization.\n    \"\"\"\n    \n    def __init__(self, max_clusters=8):\n        self.max_clusters = max_clusters\n        self.best_model = None\n        self.cluster_labels = None\n        \n    def discover_patterns(self, df, numeric_cols):\n        \"\"\"\n        Automatic clustering with K-Means and DBSCAN, selecting optimal k.\n        \"\"\"\n        print(\"\\nUNSUPERVISED PATTERN DISCOVERY\")\n        print(\"=\" * 60)\n        \n        X = df[numeric_cols].copy()\n        X = X.fillna(X.median())\n\n        \n        # Standardize\n        scaler = StandardScaler()\n        X_scaled = scaler.fit_transform(X)\n\n        X_scaled = np.nan_to_num(X_scaled)\n        \n        # 1. Determine optimal k using elbow method + silhouette\n        print(\"Finding optimal number of clusters...\")\n        inertias = []\n        silhouettes = []\n        \n        max_k = min(self.max_clusters + 1, len(X) // 10)\n        \n        for k in range(2, max_k):\n            kmeans = KMeans(n_clusters=k, random_state=42, n_init=10)\n            labels = kmeans.fit_predict(X_scaled)\n            inertias.append(kmeans.inertia_)\n            silhouettes.append(silhouette_score(X_scaled, labels))\n\n        # Find elbow using second derivative\n        if len(silhouettes) > 2:\n            inertia_diff = np.diff(inertias, 2)\n            optimal_k = np.argmax(inertia_diff) + 2\n        else:\n            optimal_k = 2\n                \n        # Validate with silhouette\n        sil_optimal_k = np.argmax(silhouettes) + 2\n        \n        # Use average of two methods\n        final_k = int(np.mean([optimal_k, sil_optimal_k]))\n        \n        print(f\"Optimal clusters: {final_k} (elbow: {optimal_k}, silhouette: {sil_optimal_k})\")\n        \n        # 2. Fit final K-Means model\n        kmeans = KMeans(\n            n_clusters=final_k, \n            random_state=42, \n            n_init=20\n        )\n        kmeans_labels = kmeans.fit_predict(X_scaled)\n        \n        # 3. Try DBSCAN for density-based clustering\n        print(\"\\nAttempting DBSCAN for density-based patterns...\")\n        dbscan = DBSCAN(eps=0.5, min_samples=5)\n        dbscan_labels = dbscan.fit_predict(X_scaled)\n        n_dbscan_clusters = len(set(dbscan_labels)) - (1 if -1 in dbscan_labels else 0)\n        \n        print(f\"DBSCAN found {n_dbscan_clusters} clusters + {(dbscan_labels == -1).sum()} noise points\")\n        \n        # Use K-Means as primary, DBSCAN for validation\n        self.cluster_labels = kmeans_labels\n        self.best_model = kmeans\n        \n        # 4. Profile each cluster\n        print(f\"\\nCLUSTER PROFILES (K={final_k}):\")\n        cluster_df = X.copy()\n        cluster_df['Cluster'] = kmeans_labels\n        \n        for cluster_id in range(final_k):\n            cluster_data = cluster_df[cluster_df['Cluster'] == cluster_id]\n            size = len(cluster_data)\n            pct = size / len(df) * 100\n            \n            print(f\"\\nCluster {cluster_id} ({size} samples, {pct:.1f}%):\")\n            \n            # Find distinguishing features (highest z-score deviation from mean)\n            cluster_mean = cluster_data[numeric_cols].mean()\n            global_mean = X[numeric_cols].mean()\n            global_std = X[numeric_cols].std()\n            \n            z_scores = (cluster_mean - global_mean) / (global_std + 1e-10)\n            top_features = z_scores.abs().sort_values(ascending=False).head(3)\n            \n            for feat, z in top_features.items():\n                direction = \"higher\" if z_scores[feat] > 0 else \"lower\"\n                print(f\"- {feat:25s}: {direction} than average (z={z:.2f})\")\n        \n        return {\n            'labels': kmeans_labels,\n            'n_clusters': final_k,\n            'silhouette': silhouettes[final_k - 2] if final_k - 2 < len(silhouettes) else silhouettes[-1],\n            'dbscan_labels': dbscan_labels,\n            'cluster_centers': kmeans.cluster_centers_\n        }","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-28T15:31:17.508171Z","iopub.execute_input":"2025-12-28T15:31:17.508407Z","iopub.status.idle":"2025-12-28T15:31:17.539148Z","shell.execute_reply.started":"2025-12-28T15:31:17.508386Z","shell.execute_reply":"2025-12-28T15:31:17.538142Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"class PredictiveSignalExtractor:\n    \"\"\"\n    Feature importance analysis using ensemble methods.\n    \"\"\"\n    \n    def __init__(self, task_type='auto'):\n        self.task_type = task_type\n        self.model = None\n        self.feature_importance = None\n        \n    def extract_signals(self, df, numeric_cols, target_col):\n        \"\"\"\n        Use Random Forest to determine feature importance.\n        \"\"\"\n        print(\"\\nPREDICTIVE SIGNAL EXTRACTION\")\n        print(\"=\" * 60)\n        \n        if target_col not in df.columns:\n            print(f\"Target column '{target_col}' not found\")\n            return None\n        \n        X = df[numeric_cols].copy()\n        X = X.fillna(X.median())\n        y = df[target_col].copy()\n        \n        # Remove target from features if present\n        if target_col in X.columns:\n            X = X.drop(columns=[target_col])\n        \n        # Auto-detect task type\n        if self.task_type == 'auto':\n            n_unique = y.nunique()\n            task = 'classification' if n_unique <= 10 else 'regression'\n        else:\n            task = self.task_type\n        \n        print(f\"Task type: {task.upper()}\")\n        print(f\"Features: {len(X.columns)}, Samples: {len(X)}\")\n        \n        # Train model\n        if task == 'classification':\n            model = RandomForestClassifier(n_estimators=200, max_depth=10, random_state=42)\n        else:\n            model = RandomForestRegressor(n_estimators=200, max_depth=10, random_state=42)\n        \n        model.fit(X, y)\n        self.model = model\n        \n        # Extract feature importance\n        importance_df = pd.DataFrame({\n            'feature': X.columns,\n            'importance': model.feature_importances_\n        }).sort_values('importance', ascending=False)\n        \n        self.feature_importance = importance_df\n        \n        # Calculate cumulative importance\n        importance_df['cumulative'] = importance_df['importance'].cumsum()\n        \n        # Find features that contribute to 80% of importance\n        n_features_80 = (importance_df['cumulative'] <= 0.8).sum()\n        \n        print(f\"\\nFEATURE IMPORTANCE RANKINGS:\")\n        print(f\"Top {n_features_80} features explain 80% of predictive power\\n\")\n        \n        for idx, row in importance_df.head(15).iterrows():\n            bar = '█' * int(row['importance'] * 50)\n            print(f\"{row['feature']:30s} {bar} {row['importance']:.4f}\")\n        \n        # Identify low-importance features for potential removal\n        low_importance = importance_df[importance_df['importance'] < 0.01]\n        if len(low_importance) > 0:\n            print(f\"\\n{len(low_importance)} features have <1% importance (consider removal)\")\n        \n        return importance_df","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-28T15:31:17.540719Z","iopub.execute_input":"2025-12-28T15:31:17.540995Z","iopub.status.idle":"2025-12-28T15:31:17.573188Z","shell.execute_reply.started":"2025-12-28T15:31:17.540977Z","shell.execute_reply":"2025-12-28T15:31:17.572088Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"class Visualizer:\n    \"\"\"\n    Interactive and static visualizations.\n    \"\"\"\n    \n    def __init__(self, output_dir='insights_output'):\n        self.output_dir = Path(output_dir)\n        self.output_dir.mkdir(exist_ok=True)\n\n    def plot_correlation_heatmap(self, corr_matrix, title=\"Correlation Heatmap\"):\n        plt.figure(figsize=(16, 14))\n        \n        sns.heatmap(corr_matrix, annot=False, fmt='.2f', cmap='RdBu_r', \n                    center=0, square=True, linewidths=0.5, cbar_kws={\"shrink\": 0.8})\n        plt.title(title, fontsize=16, fontweight='bold')\n        plt.tight_layout()\n        \n        filename = self.output_dir / 'correlation_heatmap.png'\n        plt.savefig(filename, dpi=300, bbox_inches='tight')\n        print(f\"Saved: {filename}\")\n        plt.close()\n        \n    def plot_pairplot(self, df, numeric_cols, hue_col=None, max_vars=8):\n        \"\"\"High-dimensional pair plot with KDE.\"\"\"\n        plot_cols = numeric_cols[:max_vars]\n        \n        if hue_col and hue_col in df.columns:\n            plot_df = df[plot_cols + [hue_col]].copy()\n        else:\n            plot_df = df[plot_cols].copy()\n            hue_col = None\n        \n        g = sns.pairplot(plot_df, hue=hue_col, diag_kind='kde', plot_kws={'alpha': 0.6})\n        g.fig.suptitle('Multivariate Pair Plot', y=1.02, fontsize=16, fontweight='bold')\n        \n        filename = self.output_dir / 'pairplot.png'\n        plt.savefig(filename, dpi=300, bbox_inches='tight')\n        print(f\"Saved: {filename}\")\n        plt.close()\n        \n\n    def plot_state_analysis(self, df, output_dir):\n        \"\"\"Biomass-specific: Analyze by State (WA anomaly)\"\"\"\n        fig, axes = plt.subplots(2, 2, figsize=(14, 10))\n        \n        # Dead biomass by state\n        df.boxplot(column='Dry_Dead_g', by='State', ax=axes[0, 0])\n        axes[0, 0].set_title('Dead Biomass by State')\n        axes[0, 0].set_ylabel('Dry Dead (g)')\n        \n        # Clover by state\n        df.boxplot(column='Dry_Clover_g', by='State', ax=axes[0, 1])\n        axes[0, 1].set_title('Clover Biomass by State')\n        axes[0, 1].set_ylabel('Dry Clover (g)')\n        \n        # Total biomass by state\n        df.boxplot(column='Dry_Total_g', by='State', ax=axes[1, 0])\n        axes[1, 0].set_title('Total Biomass by State')\n        axes[1, 0].set_ylabel('Dry Total (g)')\n        \n        # NDVI by state\n        df.boxplot(column='Pre_GSHH_NDVI', by='State', ax=axes[1, 1])\n        axes[1, 1].set_title('NDVI by State')\n        axes[1, 1].set_ylabel('NDVI')\n        \n        plt.suptitle('State-Level Analysis (WA Anomaly Investigation)', fontsize=14, fontweight='bold')\n        plt.tight_layout()\n        \n        filename = output_dir / 'state_analysis.png'\n        plt.savefig(filename, dpi=300, bbox_inches='tight')\n        print(f\"Saved: {filename}\")\n        plt.close()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-28T15:31:17.574338Z","iopub.execute_input":"2025-12-28T15:31:17.574639Z","iopub.status.idle":"2025-12-28T15:31:17.606173Z","shell.execute_reply.started":"2025-12-28T15:31:17.574609Z","shell.execute_reply":"2025-12-28T15:31:17.605120Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"class Biomass:\n    \"\"\"\n    Main orchestrator\n    \"\"\"\n    \n    def __init__(self, output_dir='insights_output'):\n        self.anomaly_detector = AnomalyDetector(contamination=0.05)\n        self.multivariate_analyzer = MultivariateAnalyzer(correlation_threshold=0.6)\n        self.clustering_engine = ClusteringEngine(max_clusters=8)\n        self.signal_extractor = PredictiveSignalExtractor()\n        self.visualizer = Visualizer(output_dir)\n        self.insights = {}\n        \n    def analyze(self, df):\n        \"\"\"\n        Run complete analysis pipeline.\n        \n        Args:\n            df: Input DataFrame\n            target_col: Target variable for supervised analysis\n            id_cols: Columns to exclude from analysis (IDs, text, etc.)\n        \"\"\"\n        print(\"\\n\" + \"=\"*60)\n        print(\"ANALYZING\")\n        print(\"=\"*60)\n        \n        # Identify numeric columns\n        target_cols = ['Dry_Green_g', 'Dry_Dead_g', 'Dry_Clover_g', 'GDM_g', 'Dry_Total_g']\n        id_cols = ['image_id', 'image_path']\n        \n        # Numeric features for analysis\n        numeric_cols = df.select_dtypes(include=[np.number]).columns.tolist()\n        feature_cols = [col for col in numeric_cols if col not in target_cols]\n        \n        print(f\"\\nAnalysis Configuration:\")\n        print(f\"  - Total samples: {len(df)}\")\n        print(f\"  - Feature columns: {len(feature_cols)}\")\n        print(f\"  - Target columns: {len(target_cols)}\")\n        \n        # Module 1: Anomaly Detection\n        anomaly_results = self.anomaly_detector.detect_anomalies(df, feature_cols)\n        df['is_anomaly'] = anomaly_results['is_anomaly']\n        df['anomaly_score'] = anomaly_results['anomaly_score']\n        self.insights['anomalies'] = anomaly_results\n        \n        # Module 2: Multivariate Analysis (using Dry_Total_g as primary target)\n        interaction_features, corr_matrix = self.multivariate_analyzer.analyze(\n            df, feature_cols, target_col='Dry_Total_g'\n        )\n        self.insights['interactions'] = interaction_features\n        self.insights['correlations'] = corr_matrix\n        \n        # Module 3: Clustering\n        cluster_results = self.clustering_engine.discover_patterns(df, feature_cols)\n        df['Cluster'] = cluster_results['labels']\n        self.insights['clusters'] = cluster_results\n        \n        # Module 4: Predictive Signals\n        print(\"\\n\" + \"=\"*60)\n        print(\"PREDICTIVE SIGNALS FOR EACH TARGET\")\n        print(\"=\"*60)\n        \n        importance_results = {}\n        for target in target_cols:\n            print(f\"\\n--- Analyzing: {target} ---\")\n            importance_df = self.signal_extractor.extract_signals(df, feature_cols, target)\n            importance_results[target] = importance_df\n        \n        self.insights['feature_importance'] = importance_results\n        \n        # Module 5: Visualizations\n        print(\"\\nGENERATING VISUALIZATIONS\")\n        print(\"=\" * 60)\n\n        self.visualizer.plot_correlation_heatmap(corr_matrix, \"Biomass Features Correlation Matrix\")\n        self.visualizer.plot_pairplot(df, feature_cols[:8], \n                                      hue_col='Cluster' if 'Cluster' not in feature_cols else None)\n        self.visualizer.plot_state_analysis(df, self.visualizer.output_dir)\n\n        # Save enhanced dataset\n        output_csv = self.visualizer.output_dir / 'train_with_insights.csv'\n        df.to_csv(output_csv, index=False)\n        print(f\"Saved enhanced dataset: {output_csv}\")\n        \n        print(\"\\n\" + \"=\"*60)\n        print(\"ANALYSIS COMPLETE\")\n        print(\"=\"*60)\n        print(f\"\\nAll outputs saved to: {self.visualizer.output_dir}\")\n        \n        return df, self.insights\n    \n    def generate_biomass_report(self, df, output_file='biomass_report.txt'):\n        \"\"\"Generate a comprehensive text report of findings.\"\"\"\n        report_path = self.visualizer.output_dir / output_file\n        \n        with open(report_path, 'w') as f:\n            f.write(\"=\"*70 + \"\\n\")\n            f.write(\"CSIRO BIOMASS COMPETITION - REPORT\\n\")\n            f.write(\"=\"*70 + \"\\n\\n\")\n            \n            # Executive Summary\n            f.write(\"EXECUTIVE SUMMARY\\n\")\n            f.write(\"-\"*70 + \"\\n\")\n            f.write(f\"Dataset: {len(df)} pasture images analyzed\\n\")\n            f.write(f\"States: {df['State'].unique().tolist()}\\n\")\n            f.write(f\"Date Range: {df['Sampling_Date'].min()} to {df['Sampling_Date'].max()}\\n\\n\")\n            \n            # Key Finding 1: WA Dead Matter Anomaly\n            f.write(\"KEY FINDING 1: Western Australia Dead Matter Anomaly\\n\")\n            f.write(\"-\"*70 + \"\\n\")\n            wa_dead = df[df['State'] == 'WA']['Dry_Dead_g']\n            other_dead = df[df['State'] != 'WA']['Dry_Dead_g']\n            f.write(f\"WA Dead Biomass:    mean={wa_dead.mean():.2f}g, zeros={(wa_dead==0).sum()/len(wa_dead)*100:.1f}%\\n\")\n            f.write(f\"Other Dead Biomass: mean={other_dead.mean():.2f}g, zeros={(other_dead==0).sum()/len(other_dead)*100:.1f}%\\n\")\n            f.write(\"RECOMMENDATION: Apply state-conditional post-processing for Dead predictions\\n\\n\")\n            \n            # Anomalies\n            f.write(\"ANOMALY DETECTION RESULTS\\n\")\n            f.write(\"-\"*70 + \"\\n\")\n            if 'anomalies' in self.insights:\n                n_anomalies = self.insights['anomalies']['is_anomaly'].sum()\n                f.write(f\"Detected: {n_anomalies} anomalous samples ({n_anomalies/len(df)*100:.2f}%)\\n\\n\")\n            \n            # Clustering\n            f.write(\"CLUSTER ANALYSIS\\n\")\n            f.write(\"-\"*70 + \"\\n\")\n            if 'clusters' in self.insights:\n                f.write(f\"Optimal clusters: {self.insights['clusters']['n_clusters']}\\n\")\n                f.write(f\"Silhouette score: {self.insights['clusters']['silhouette']:.3f}\\n\\n\")\n            \n            # Feature Importance Summary\n            f.write(\"TOP PREDICTIVE FEATURES (for Dry_Total_g)\\n\")\n            f.write(\"-\"*70 + \"\\n\")\n            if 'feature_importance' in self.insights:\n                importance_df = self.insights['feature_importance']['Dry_Total_g']\n                for _, row in importance_df.head(15).iterrows():\n                    f.write(f\"  {row['feature']:40s} {row['importance']:.4f}\\n\")\n        \n        print(f\"Report saved: {report_path}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-28T15:31:17.607289Z","iopub.execute_input":"2025-12-28T15:31:17.607482Z","iopub.status.idle":"2025-12-28T15:31:17.632082Z","shell.execute_reply.started":"2025-12-28T15:31:17.607464Z","shell.execute_reply":"2025-12-28T15:31:17.630955Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def main():\n    \"\"\"Main execution function.\"\"\"\n    \n    # Configuration\n    TRAIN_CSV_PATH = '/kaggle/input/csiro-biomass/train.csv'\n    OUTPUT_DIR = 'biomass_insights'\n    \n    print(\"CSIRO Biomass Competition - Insight Analysis\")\n    print(\"=\"*60)\n    \n    # Step 1: Load and preprocess data\n    df_wide = load_and_pivot_biomass_data(TRAIN_CSV_PATH)\n    \n    # Step 2: Run comprehensive analysis\n    engine = Biomass(output_dir=OUTPUT_DIR)\n    df_enhanced, insights = engine.analyze(df_wide)\n    \n    # Step 3: Generate report\n    engine.generate_biomass_report(df_enhanced)\n    \n    print(\"\\nAnalysis Complete!\")\n    print(f\"All outputs saved to: {OUTPUT_DIR}/\")\n\nif __name__ == \"__main__\":\n    main()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-28T15:31:17.633442Z","iopub.execute_input":"2025-12-28T15:31:17.633740Z","iopub.status.idle":"2025-12-28T15:31:31.808947Z","shell.execute_reply.started":"2025-12-28T15:31:17.633722Z","shell.execute_reply":"2025-12-28T15:31:31.807864Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}