{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":114397,"databundleVersionId":13696770,"sourceType":"competition"}],"dockerImageVersionId":31089,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport os\nfrom PIL import Image\nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\nimport torchvision.transforms as transforms\nfrom torchvision import models\nfrom torch.utils.data import Dataset, DataLoader\nfrom sklearn.cluster import KMeans, AgglomerativeClustering, SpectralClustering\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.decomposition import PCA\nfrom sklearn.metrics import silhouette_score, calinski_harabasz_score, davies_bouldin_score\nfrom sklearn.mixture import GaussianMixture\nfrom scipy.cluster.hierarchy import linkage, fcluster\nfrom scipy.spatial.distance import pdist, squareform\nfrom tqdm import tqdm\nimport warnings\nimport pickle\nimport gc\nfrom typing import Dict, List, Tuple, Optional\nfrom collections import defaultdict\nimport time\nwarnings.filterwarnings('ignore')\n\n# ============================================================================\n# ENHANCED CONFIGURATION\n# ============================================================================\n\nclass Config:\n    # Paths\n    base_path = '/kaggle/input/biotrove-clustering'\n    metadata_path = os.path.join(base_path, 'metadata.csv')\n    images_dir = os.path.join(base_path, 'images/images')\n    submission_path = 'submission.csv'\n    features_cache_path = 'features_cache.pkl'\n    \n    # Image processing\n    img_size = 224\n    batch_size = 64 if torch.cuda.is_available() else 32\n    num_workers = 4\n    \n    # Feature extraction\n    use_ensemble = True  # Use multiple models\n    models_ensemble = ['resnet50', 'efficientnet_b3']  # Add 'vit_b_16' if you have memory\n    use_tta = False  # Disabled TTA for speed - set to True for better accuracy\n    tta_transforms = 3\n    use_cache = True  # Cache extracted features\n    \n    # Dimensionality reduction\n    use_pca = True\n    pca_components = 512\n    \n    # Clustering\n    clustering_method = 'ensemble'  # 'kmeans', 'hierarchical', 'gmm', 'ensemble'\n    n_families = None\n    n_genera = None  \n    n_species = None\n    family_ratio = 7  # genera per family\n    genus_ratio = 4  # species per genus\n    use_family_priors = True\n    refine_clusters = True\n    min_cluster_size = 5\n    max_cluster_size = 500\n    \n    # Advanced features\n    use_contrastive_learning = False  # Disabled for speed - set to True for better accuracy\n    contrastive_temperature = 0.5\n    contrastive_epochs = 3  # Reduced from 5 for speed\n    contrastive_batch_size = 512  # Increased for speed\n    \n    # System\n    seed = 42\n    device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\n    use_mixed_precision = True  # Use mixed precision for speed\n\nconfig = Config()\n\n# Set seeds for reproducibility\nnp.random.seed(config.seed)\ntorch.manual_seed(config.seed)\nif torch.cuda.is_available():\n    torch.cuda.manual_seed_all(config.seed)\n\nprint(f\"Using device: {config.device}\")\nif config.device.type == 'cuda':\n    print(f\"GPU: {torch.cuda.get_device_name(0)}\")\n\n# ============================================================================\n# ENHANCED DATASET WITH TTA\n# ============================================================================\n\nclass BioTroveDataset(Dataset):\n    def __init__(self, metadata_df, images_dir, transform=None, augmentation_transform=None, use_tta=False):\n        self.metadata_df = metadata_df\n        self.images_dir = images_dir\n        self.transform = transform\n        self.augmentation_transform = augmentation_transform\n        self.use_tta = use_tta\n        \n    def __len__(self):\n        return len(self.metadata_df)\n    \n    def __getitem__(self, idx):\n        row = self.metadata_df.iloc[idx]\n        hash_id = row['hash_id']\n        family = row['family']\n        \n        img_path = os.path.join(self.images_dir, f\"{hash_id}.jpg\")\n        try:\n            image = Image.open(img_path).convert('RGB')\n        except:\n            image = Image.new('RGB', (config.img_size, config.img_size), (128, 128, 128))\n        \n        if self.use_tta and self.augmentation_transform:\n            images = []\n            images.append(self.transform(image))\n            for _ in range(config.tta_transforms - 1):\n                images.append(self.augmentation_transform(image))\n            return torch.stack(images), hash_id, family\n        else:\n            if self.transform:\n                image = self.transform(image)\n            return image, hash_id, family\n\n# ============================================================================\n# CONTRASTIVE LEARNING MODULE (FIXED)\n# ============================================================================\n\nclass ContrastiveLearning(nn.Module):\n    def __init__(self, feature_dim, projection_dim=128):\n        super().__init__()\n        self.projection_head = nn.Sequential(\n            nn.Linear(feature_dim, feature_dim),\n            nn.ReLU(),\n            nn.Linear(feature_dim, projection_dim)\n        )\n        \n    def forward(self, features):\n        \"\"\"Pass features through projection head\"\"\"\n        return self.projection_head(features)\n        \n    def supervised_contrastive_loss(self, projected_features, labels, temperature=0.5):\n        \"\"\"Compute contrastive loss on projected features\"\"\"\n        device = projected_features.device\n        batch_size = projected_features.shape[0]\n        \n        # Normalize features\n        projected_features = F.normalize(projected_features, dim=1)\n        \n        # Compute similarity matrix\n        similarity = torch.matmul(projected_features, projected_features.T) / temperature\n        \n        # Create label mask\n        labels = labels.unsqueeze(1)\n        mask = (labels == labels.T).float().to(device)\n        mask.fill_diagonal_(0)\n        \n        # Numerical stability\n        logits_max, _ = torch.max(similarity, dim=1, keepdim=True)\n        logits = similarity - logits_max.detach()\n        \n        # Compute log probabilities\n        exp_logits = torch.exp(logits)\n        log_prob = logits - torch.log(exp_logits.sum(dim=1, keepdim=True))\n        \n        # Compute mean log likelihood for positive pairs\n        mask_sum = mask.sum(dim=1)\n        mask_sum = torch.clamp(mask_sum, min=1e-6)  # Avoid division by zero\n        mean_log_prob_pos = (mask * log_prob).sum(dim=1) / mask_sum\n        \n        # Loss is negative log likelihood\n        loss = -mean_log_prob_pos.mean()\n        \n        return loss\n\n# ============================================================================\n# FEATURE EXTRACTION WITH ENSEMBLE (OPTIMIZED)\n# ============================================================================\n\nclass FeatureExtractor:\n    def __init__(self, model_names, device):\n        self.device = device\n        self.models = self._load_models(model_names)\n        \n    def _load_models(self, model_names):\n        models_dict = {}\n        \n        for model_name in model_names:\n            print(f\"Loading {model_name}...\")\n            \n            if model_name == 'resnet50':\n                base_model = models.resnet50(pretrained=True)\n                model = nn.Sequential(*list(base_model.children())[:-1], nn.Flatten())\n                \n            elif model_name == 'efficientnet_b3':\n                base_model = models.efficientnet_b3(pretrained=True)\n                model = nn.Sequential(\n                    base_model.features,\n                    nn.AdaptiveAvgPool2d(1),\n                    nn.Flatten()\n                )\n                \n            elif model_name == 'vit_b_16':\n                base_model = models.vit_b_16(pretrained=True)\n                model = nn.Sequential(*list(base_model.children())[:-1])\n                \n            else:\n                continue\n            \n            model = model.to(self.device)\n            model.eval()\n            models_dict[model_name] = model\n            \n        return models_dict\n    \n    @torch.cuda.amp.autocast(enabled=True)  # Mixed precision for speed\n    def extract_features(self, dataloader, use_cache=True):\n        # Try to load from cache\n        if use_cache and os.path.exists(config.features_cache_path):\n            print(\"Loading features from cache...\")\n            with open(config.features_cache_path, 'rb') as f:\n                cache_data = pickle.load(f)\n            return cache_data['features'], cache_data['hash_ids'], cache_data['families']\n        \n        all_features = {name: [] for name in self.models.keys()}\n        hash_ids = []\n        families = []\n        \n        with torch.no_grad():\n            for batch_data in tqdm(dataloader, desc=\"Extracting features\"):\n                if config.use_tta and len(batch_data[0].shape) == 5:\n                    batch_imgs, batch_ids, batch_families = batch_data\n                    batch_size, n_aug = batch_imgs.shape[:2]\n                    batch_imgs = batch_imgs.view(-1, *batch_imgs.shape[2:])\n                else:\n                    batch_imgs, batch_ids, batch_families = batch_data\n                    batch_size = batch_imgs.shape[0]\n                    n_aug = 1\n                \n                batch_imgs = batch_imgs.to(self.device)\n                \n                for model_name, model in self.models.items():\n                    if config.use_mixed_precision:\n                        with torch.cuda.amp.autocast():\n                            features = model(batch_imgs)\n                    else:\n                        features = model(batch_imgs)\n                    \n                    if config.use_tta and n_aug > 1:\n                        features = features.view(batch_size, n_aug, -1).mean(dim=1)\n                    \n                    all_features[model_name].append(features.cpu().numpy())\n                \n                hash_ids.extend(batch_ids)\n                families.extend(batch_families)\n        \n        # Concatenate features from all models\n        ensemble_features = []\n        for model_name in self.models.keys():\n            model_features = np.vstack(all_features[model_name])\n            # Normalize each model's features\n            scaler = StandardScaler()\n            model_features = scaler.fit_transform(model_features)\n            ensemble_features.append(model_features)\n        \n        if len(ensemble_features) > 1:\n            features = np.hstack(ensemble_features)\n        else:\n            features = ensemble_features[0]\n        \n        # Cache the features\n        if use_cache:\n            print(\"Caching extracted features...\")\n            os.makedirs(os.path.dirname(config.features_cache_path) or '.', exist_ok=True)\n            with open(config.features_cache_path, 'wb') as f:\n                pickle.dump({\n                    'features': features,\n                    'hash_ids': hash_ids,\n                    'families': families\n                }, f)\n        \n        return features, hash_ids, families\n\n# ============================================================================\n# ADVANCED CLUSTERING (OPTIMIZED)\n# ============================================================================\n\nclass HierarchicalClusterer:\n    def __init__(self, n_families, n_genera, n_species):\n        self.n_families = n_families\n        self.n_genera = n_genera\n        self.n_species = n_species\n        \n    def adaptive_hierarchical_clustering(self, features, family_labels=None):\n        print(\"Computing hierarchical clustering...\")\n        \n        # Use faster distance computation for large datasets\n        if len(features) > 10000:\n            # Subsample for distance computation\n            sample_idx = np.random.choice(len(features), 10000, replace=False)\n            features_sample = features[sample_idx]\n            distances = pdist(features_sample, metric='cosine')\n            Z = linkage(distances, method='ward')\n            \n            # Predict full clusters using KMeans initialized with hierarchical results\n            species_clusters_sample = fcluster(Z, self.n_species, criterion='maxclust') - 1\n            genus_clusters_sample = fcluster(Z, self.n_genera, criterion='maxclust') - 1\n            family_clusters_sample = fcluster(Z, self.n_families, criterion='maxclust') - 1\n            \n            # Use sampled clusters to initialize KMeans for full dataset\n            species_clusters = self._kmeans_from_sample(features, features_sample, species_clusters_sample, self.n_species)\n            genus_clusters = self._kmeans_from_sample(features, features_sample, genus_clusters_sample, self.n_genera)\n            family_clusters = self._kmeans_from_sample(features, features_sample, family_clusters_sample, self.n_families)\n        else:\n            distances = pdist(features, metric='cosine')\n            Z = linkage(distances, method='ward')\n            \n            species_clusters = fcluster(Z, self.n_species, criterion='maxclust') - 1\n            genus_clusters = fcluster(Z, self.n_genera, criterion='maxclust') - 1\n            family_clusters = fcluster(Z, self.n_families, criterion='maxclust') - 1\n        \n        if family_labels is not None and config.use_family_priors:\n            family_clusters = family_labels\n        \n        return family_clusters, genus_clusters, species_clusters\n    \n    def _kmeans_from_sample(self, features, features_sample, clusters_sample, n_clusters):\n        \"\"\"Initialize KMeans with centers from sampled hierarchical clustering\"\"\"\n        centers = []\n        for i in range(n_clusters):\n            mask = clusters_sample == i\n            if mask.sum() > 0:\n                centers.append(features_sample[mask].mean(axis=0))\n            else:\n                centers.append(features_sample[np.random.choice(len(features_sample))])\n        \n        centers = np.array(centers)\n        kmeans = KMeans(n_clusters=n_clusters, init=centers, n_init=1, random_state=42)\n        return kmeans.fit_predict(features)\n    \n    def ensemble_clustering(self, features, family_labels=None):\n        print(\"Performing ensemble clustering...\")\n        \n        methods_predictions = []\n        \n        # Method 1: Hierarchical\n        hier_f, hier_g, hier_s = self.adaptive_hierarchical_clustering(features, family_labels)\n        methods_predictions.append((hier_f, hier_g, hier_s))\n        \n        # Method 2: KMeans (optimized with fewer iterations)\n        kmeans_s = KMeans(n_clusters=self.n_species, random_state=42, n_init=3, max_iter=100).fit_predict(features)\n        kmeans_g = KMeans(n_clusters=self.n_genera, random_state=42, n_init=3, max_iter=100).fit_predict(features)\n        if family_labels is not None and config.use_family_priors:\n            kmeans_f = family_labels\n        else:\n            kmeans_f = KMeans(n_clusters=self.n_families, random_state=42, n_init=3, max_iter=100).fit_predict(features)\n        methods_predictions.append((kmeans_f, kmeans_g, kmeans_s))\n        \n        # Skip GMM for speed (optional - uncomment if you want better accuracy)\n        # gmm_s = GaussianMixture(n_components=self.n_species, covariance_type='diag', random_state=42).fit_predict(features)\n        # gmm_g = GaussianMixture(n_components=self.n_genera, covariance_type='diag', random_state=42).fit_predict(features)\n        # if family_labels is not None and config.use_family_priors:\n        #     gmm_f = family_labels\n        # else:\n        #     gmm_f = GaussianMixture(n_components=self.n_families, covariance_type='diag', random_state=42).fit_predict(features)\n        # methods_predictions.append((gmm_f, gmm_g, gmm_s))\n        \n        # Consensus clustering\n        family_consensus = self._consensus_clustering([m[0] for m in methods_predictions], self.n_families)\n        genus_consensus = self._consensus_clustering([m[1] for m in methods_predictions], self.n_genera)\n        species_consensus = self._consensus_clustering([m[2] for m in methods_predictions], self.n_species)\n        \n        return family_consensus, genus_consensus, species_consensus\n    \n    def _consensus_clustering(self, predictions_list, n_clusters):\n        n_samples = len(predictions_list[0])\n        n_methods = len(predictions_list)\n        \n        coassoc = np.zeros((n_samples, n_samples))\n        \n        for pred in predictions_list:\n            for i in range(n_samples):\n                for j in range(i+1, n_samples):\n                    if pred[i] == pred[j]:\n                        coassoc[i, j] += 1\n                        coassoc[j, i] += 1\n        \n        coassoc /= n_methods\n        distance_matrix = 1 - coassoc\n        Z = linkage(squareform(distance_matrix), method='average')\n        consensus = fcluster(Z, n_clusters, criterion='maxclust') - 1\n        \n        return consensus\n    \n    def refine_clusters(self, family_clusters, genus_clusters, species_clusters):\n        print(\"Refining clusters for hierarchical consistency...\")\n        \n        n_samples = len(species_clusters)\n        \n        # Ensure species are nested within genera\n        species_to_genus = {}\n        for i in range(n_samples):\n            species_id = species_clusters[i]\n            genus_id = genus_clusters[i]\n            \n            if species_id not in species_to_genus:\n                species_to_genus[species_id] = genus_id\n            elif species_to_genus[species_id] != genus_id:\n                genus_counts = defaultdict(int)\n                for j in range(n_samples):\n                    if species_clusters[j] == species_id:\n                        genus_counts[genus_clusters[j]] += 1\n                most_common_genus = max(genus_counts, key=genus_counts.get)\n                species_to_genus[species_id] = most_common_genus\n        \n        for i in range(n_samples):\n            genus_clusters[i] = species_to_genus[species_clusters[i]]\n        \n        # Ensure genera are nested within families\n        genus_to_family = {}\n        for i in range(n_samples):\n            genus_id = genus_clusters[i]\n            family_id = family_clusters[i]\n            \n            if genus_id not in genus_to_family:\n                genus_to_family[genus_id] = family_id\n            elif genus_to_family[genus_id] != family_id:\n                family_counts = defaultdict(int)\n                for j in range(n_samples):\n                    if genus_clusters[j] == genus_id:\n                        family_counts[family_clusters[j]] += 1\n                most_common_family = max(family_counts, key=family_counts.get)\n                genus_to_family[genus_id] = most_common_family\n        \n        for i in range(n_samples):\n            family_clusters[i] = genus_to_family[genus_clusters[i]]\n        \n        return family_clusters, genus_clusters, species_clusters\n\n# ============================================================================\n# MAIN PIPELINE (FIXED AND OPTIMIZED)\n# ============================================================================\n\ndef main():\n    print(\"=\"*60)\n    print(\"BioTrove Enhanced Clustering Pipeline (Optimized)\")\n    print(\"=\"*60)\n    \n    # Load metadata\n    print(\"\\nLoading metadata...\")\n    metadata_df = pd.read_csv(config.metadata_path)\n    print(f\"Total images: {len(metadata_df)}\")\n    \n    # Estimate number of clusters\n    n_families = metadata_df['family'].nunique()\n    config.n_families = n_families\n    config.n_genera = int(n_families * config.family_ratio)\n    config.n_species = int(config.n_genera * config.genus_ratio)\n    \n    print(f\"Number of families: {config.n_families}\")\n    print(f\"Estimated genera: {config.n_genera}\")\n    print(f\"Estimated species: {config.n_species}\")\n    \n    # Define transforms\n    normalize = transforms.Normalize(mean=[0.485, 0.456, 0.406], std=[0.229, 0.224, 0.225])\n    \n    base_transform = transforms.Compose([\n        transforms.Resize((config.img_size, config.img_size)),\n        transforms.ToTensor(),\n        normalize\n    ])\n    \n    augmentation_transform = transforms.Compose([\n        transforms.Resize((config.img_size + 20, config.img_size + 20)),\n        transforms.RandomCrop(config.img_size),\n        transforms.RandomHorizontalFlip(p=0.5),\n        transforms.RandomRotation(degrees=10),\n        transforms.ColorJitter(brightness=0.1, contrast=0.1, saturation=0.1),\n        transforms.ToTensor(),\n        normalize\n    ])\n    \n    # Create dataset and dataloader\n    dataset = BioTroveDataset(\n        metadata_df, \n        config.images_dir, \n        transform=base_transform,\n        augmentation_transform=augmentation_transform if config.use_tta else None,\n        use_tta=config.use_tta\n    )\n    \n    dataloader = DataLoader(\n        dataset, \n        batch_size=config.batch_size,\n        shuffle=False, \n        num_workers=config.num_workers,\n        pin_memory=True if config.device.type == 'cuda' else False\n    )\n    \n    # Feature extraction\n    print(\"\\n\" + \"=\"*60)\n    print(\"Feature Extraction Phase\")\n    print(\"=\"*60)\n    \n    if config.use_ensemble:\n        model_names = config.models_ensemble\n    else:\n        model_names = ['resnet50']\n    \n    extractor = FeatureExtractor(model_names, config.device)\n    features, hash_ids, families = extractor.extract_features(dataloader, use_cache=config.use_cache)\n    print(f\"Features shape: {features.shape}\")\n    \n    # Clear GPU memory\n    if config.device.type == 'cuda':\n        torch.cuda.empty_cache()\n        gc.collect()\n    \n    # Apply contrastive learning (FIXED)\n    if config.use_contrastive_learning:\n        print(\"\\nApplying contrastive learning...\")\n        \n        feature_dim = features.shape[1]\n        contrastive_module = ContrastiveLearning(feature_dim, projection_dim=128)\n        contrastive_module = contrastive_module.to(config.device)\n        \n        # Convert to tensors\n        features_tensor = torch.FloatTensor(features).to(config.device)\n        family_to_id = {family: idx for idx, family in enumerate(sorted(set(families)))}\n        family_ids_tensor = torch.LongTensor([family_to_id[f] for f in families]).to(config.device)\n        \n        # Train contrastive projection head\n        optimizer = torch.optim.Adam(contrastive_module.projection_head.parameters(), lr=0.001)\n        \n        # Use larger batch size for speed\n        batch_size = config.contrastive_batch_size\n        \n        # Enable mixed precision for speed\n        scaler = torch.cuda.amp.GradScaler() if config.use_mixed_precision and config.device.type == 'cuda' else None\n        \n        for epoch in range(config.contrastive_epochs):\n            epoch_loss = 0\n            n_batches = 0\n            \n            indices = torch.randperm(len(features_tensor))\n            \n            for i in range(0, len(features_tensor), batch_size):\n                batch_indices = indices[i:i+batch_size]\n                batch_features = features_tensor[batch_indices]\n                batch_labels = family_ids_tensor[batch_indices]\n                \n                optimizer.zero_grad()\n                \n                if scaler:\n                    with torch.cuda.amp.autocast():\n                        # FIXED: Pass features through projection head first\n                        projected_features = contrastive_module(batch_features)\n                        loss = contrastive_module.supervised_contrastive_loss(\n                            projected_features, batch_labels, config.contrastive_temperature\n                        )\n                    \n                    scaler.scale(loss).backward()\n                    scaler.step(optimizer)\n                    scaler.update()\n                else:\n                    # FIXED: Pass features through projection head first\n                    projected_features = contrastive_module(batch_features)\n                    loss = contrastive_module.supervised_contrastive_loss(\n                        projected_features, batch_labels, config.contrastive_temperature\n                    )\n                    loss.backward()\n                    optimizer.step()\n                \n                epoch_loss += loss.item()\n                n_batches += 1\n            \n            avg_loss = epoch_loss / n_batches\n            print(f\"Epoch {epoch+1}/{config.contrastive_epochs}, Loss: {avg_loss:.4f}\")\n        \n        # Extract refined features\n        with torch.no_grad():\n            contrastive_module.eval()\n            refined_features = []\n            \n            for i in range(0, len(features_tensor), batch_size):\n                batch_features = features_tensor[i:i+batch_size]\n                batch_refined = contrastive_module(batch_features)\n                refined_features.append(batch_refined.cpu().numpy())\n            \n            features = np.vstack(refined_features)\n        \n        print(\"Contrastive learning completed\")\n        \n        # Clear GPU memory\n        if config.device.type == 'cuda':\n            torch.cuda.empty_cache()\n    \n    # Normalize and reduce dimensionality\n    print(\"\\nPreprocessing features...\")\n    scaler = StandardScaler()\n    features = scaler.fit_transform(features)\n    \n    if config.use_pca and features.shape[1] > config.pca_components:\n        print(f\"Applying PCA to reduce dimensions to {config.pca_components}...\")\n        pca = PCA(n_components=config.pca_components, random_state=42)\n        features = pca.fit_transform(features)\n        print(f\"Explained variance ratio: {pca.explained_variance_ratio_.sum():.3f}\")\n    \n    # Prepare family labels\n    family_to_id = {family: idx for idx, family in enumerate(sorted(set(families)))}\n    family_ids = np.array([family_to_id[f] for f in families])\n    \n    # Clustering\n    print(\"\\n\" + \"=\"*60)\n    print(\"Clustering Phase\")\n    print(\"=\"*60)\n    \n    clusterer = HierarchicalClusterer(config.n_families, config.n_genera, config.n_species)\n    \n    if config.clustering_method == 'ensemble':\n        family_clusters, genus_clusters, species_clusters = clusterer.ensemble_clustering(\n            features, family_ids if config.use_family_priors else None\n        )\n    elif config.clustering_method == 'hierarchical':\n        family_clusters, genus_clusters, species_clusters = clusterer.adaptive_hierarchical_clustering(\n            features, family_ids if config.use_family_priors else None\n        )\n    elif config.clustering_method == 'kmeans':\n        species_clusters = KMeans(n_clusters=config.n_species, random_state=42, n_init=3, max_iter=100).fit_predict(features)\n        genus_clusters = KMeans(n_clusters=config.n_genera, random_state=42, n_init=3, max_iter=100).fit_predict(features)\n        if config.use_family_priors:\n            family_clusters = family_ids\n        else:\n            family_clusters = KMeans(n_clusters=config.n_families, random_state=42, n_init=3, max_iter=100).fit_predict(features)\n    elif config.clustering_method == 'gmm':\n        species_clusters = GaussianMixture(n_components=config.n_species, covariance_type='diag', random_state=42).fit_predict(features)\n        genus_clusters = GaussianMixture(n_components=config.n_genera, covariance_type='diag', random_state=42).fit_predict(features)\n        if config.use_family_priors:\n            family_clusters = family_ids\n        else:\n            family_clusters = GaussianMixture(n_components=config.n_families, covariance_type='diag', random_state=42).fit_predict(features)\n    \n    # Refine clusters for consistency\n    if config.refine_clusters:\n        family_clusters, genus_clusters, species_clusters = clusterer.refine_clusters(\n            family_clusters, genus_clusters, species_clusters\n        )\n    \n    # Evaluate clustering (simplified for speed)\n    print(\"\\n\" + \"=\"*60)\n    print(\"Clustering Evaluation\")\n    print(\"=\"*60)\n    \n    try:\n        # Subsample for faster computation\n        sample_size = min(2000, len(features))  # Reduced from 5000 for speed\n        indices = np.random.choice(len(features), sample_size, replace=False)\n        features_sample = features[indices]\n        \n        # Compute only silhouette score for speed\n        sil_species = silhouette_score(features_sample, species_clusters[indices])\n        print(f\"Silhouette score for species: {sil_species:.3f}\")\n        \n    except Exception as e:\n        print(f\"Could not compute evaluation metrics: {e}\")\n    \n    # Create submission\n    print(\"\\n\" + \"=\"*60)\n    print(\"Creating Submission\")\n    print(\"=\"*60)\n    \n    submission_df = pd.DataFrame({\n        'hash_id': hash_ids,\n        'family_cluster': family_clusters,\n        'genus_cluster': genus_clusters,\n        'species_cluster': species_clusters\n    })\n    \n    # Handle missing samples\n    all_hash_ids = set(metadata_df['hash_id'].values)\n    submitted_hash_ids = set(submission_df['hash_id'].values)\n    \n    if all_hash_ids != submitted_hash_ids:\n        missing_ids = all_hash_ids - submitted_hash_ids\n        print(f\"Warning: {len(missing_ids)} hash_ids are missing. Adding with default clusters...\")\n        \n        for hash_id in missing_ids:\n            submission_df = pd.concat([\n                submission_df,\n                pd.DataFrame({\n                    'hash_id': [hash_id],\n                    'family_cluster': [0],\n                    'genus_cluster': [0],\n                    'species_cluster': [0]\n                })\n            ], ignore_index=True)\n    \n    # Sort and save\n    submission_df = submission_df.sort_values('hash_id').reset_index(drop=True)\n    submission_df.to_csv(config.submission_path, index=False)\n    print(f\"Submission saved to {config.submission_path}\")\n    \n    # Display statistics\n    print(\"\\n\" + \"=\"*60)\n    print(\"Final Statistics\")\n    print(\"=\"*60)\n    print(f\"Unique family clusters: {submission_df['family_cluster'].nunique()}\")\n    print(f\"Unique genus clusters: {submission_df['genus_cluster'].nunique()}\")\n    print(f\"Unique species clusters: {submission_df['species_cluster'].nunique()}\")\n    \n    # Cluster size statistics\n    print(\"\\nCluster size statistics:\")\n    for level, col in [('Family', 'family_cluster'), ('Genus', 'genus_cluster'), ('Species', 'species_cluster')]:\n        sizes = submission_df[col].value_counts()\n        print(f\"\\n{level}:\")\n        print(f\"  Mean size: {sizes.mean():.1f}\")\n        print(f\"  Median size: {sizes.median():.1f}\")\n        print(f\"  Min size: {sizes.min()}\")\n        print(f\"  Max size: {sizes.max()}\")\n    \n    # Show sample\n    print(\"\\nSample submission:\")\n    print(submission_df.head(10))\n    \n    print(\"\\n\" + \"=\"*60)\n    print(\"✅ Enhanced clustering pipeline completed!\")\n    print(\"=\"*60)\n    \n    return submission_df\n\nif __name__ == \"__main__\":\n    start_time = time.time()\n    submission = main()\n    elapsed_time = time.time() - start_time\n    print(f\"\\nTotal execution time: {elapsed_time/60:.2f} minutes\")","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}