{"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":114201,"databundleVersionId":13622514,"sourceType":"competition"}],"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":"# ========================= Install Required Packages =========================\n!pip install -q rasterio albumentations\n\nimport os\nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\nimport numpy as np\nimport pandas as pd\nimport rasterio\nfrom torch.utils.data import Dataset, DataLoader\nfrom torch.optim.lr_scheduler import CosineAnnealingWarmRestarts\nfrom sklearn.metrics import f1_score\nfrom tqdm import tqdm\nimport warnings\nwarnings.filterwarnings('ignore')\n\n# ========================= Configuration =========================\nclass Config:\n    # Data paths\n    train_metadata_path = \"/kaggle/input/geoplant-at-paiss/GLC25_PA_metadata_train.csv\"\n    test_metadata_path = \"/kaggle/input/geoplant-at-paiss/GLC25_PA_metadata_test.csv\"\n    \n    # Modality paths\n    bioclim_train_path = \"/kaggle/input/geoplant-at-paiss/BioclimTimeSeries/cubes/PA-train\"\n    bioclim_test_path = \"/kaggle/input/geoplant-at-paiss/BioclimTimeSeries/cubes/PA-test\"\n    sentinel_train_path = \"/kaggle/input/geoplant-at-paiss/SatelitePatches/PA-train\"\n    sentinel_test_path = \"/kaggle/input/geoplant-at-paiss/SatelitePatches/PA-test\"\n    landsat_train_path = \"/kaggle/input/geoplant-at-paiss/SateliteTimeSeries-Landsat/cubes/PA-train\"\n    landsat_test_path = \"/kaggle/input/geoplant-at-paiss/SateliteTimeSeries-Landsat/cubes/PA-test\"\n    \n    # Training parameters\n    num_classes = 11255\n    batch_size = 32\n    num_epochs = 15\n    learning_rate = 3e-4\n    weight_decay = 1e-4\n    val_split = 0.1\n    \n    # Model parameters\n    use_bioclim = True\n    use_sentinel = True\n    use_landsat = True\n    hidden_dim = 768\n    dropout = 0.2\n    \n    # Advanced training\n    label_smoothing = 0.1\n    mixup_alpha = 0.2\n    gradient_clip = 1.0\n    warmup_epochs = 2\n    \n    # Device\n    device = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\n    seed = 42\n\n# ========================= Data Loading =========================\ndef set_seed(seed):\n    torch.manual_seed(seed)\n    np.random.seed(seed)\n    if torch.cuda.is_available():\n        torch.cuda.manual_seed_all(seed)\n        torch.backends.cudnn.deterministic = True\n        torch.backends.cudnn.benchmark = False\n\ndef construct_patch_path(data_path, survey_id):\n    \"\"\"Construct Sentinel-2 patch path\"\"\"\n    path = data_path\n    for d in (str(survey_id)[-2:], str(survey_id)[-4:-2]):\n        path = os.path.join(path, d)\n    return os.path.join(path, f\"{survey_id}.tiff\")\n\ndef quantile_normalize(band, low=2, high=98):\n    \"\"\"Normalize band using quantile normalization\"\"\"\n    sorted_band = np.sort(band.flatten())\n    quantiles = np.percentile(sorted_band, np.linspace(low, high, len(sorted_band)))\n    normalized_band = np.interp(band.flatten(), sorted_band, quantiles).reshape(band.shape)\n    \n    min_val, max_val = np.min(normalized_band), np.max(normalized_band)\n    if max_val == min_val:\n        return np.zeros_like(normalized_band, dtype=np.float32)\n    return ((normalized_band - min_val) / (max_val - min_val)).astype(np.float32)\n\nclass MultiModalDataset(Dataset):\n    def __init__(self, config, metadata, is_train=True):\n        self.config = config\n        self.is_train = is_train\n        self.metadata = metadata\n        \n        if is_train:\n            self.metadata = self.metadata.dropna(subset=\"speciesId\").reset_index(drop=True)\n            self.metadata['speciesId'] = self.metadata['speciesId'].astype(int)\n            self.label_dict = self.metadata.groupby('surveyId')['speciesId'].apply(list).to_dict()\n            self.metadata = self.metadata.drop_duplicates(subset=\"surveyId\").reset_index(drop=True)\n        \n        # Set paths based on train/test\n        if is_train:\n            self.bioclim_path = config.bioclim_train_path\n            self.sentinel_path = config.sentinel_train_path\n            self.landsat_path = config.landsat_train_path\n        else:\n            self.bioclim_path = config.bioclim_test_path\n            self.sentinel_path = config.sentinel_test_path\n            self.landsat_path = config.landsat_test_path\n    \n    def __len__(self):\n        return len(self.metadata)\n    \n    def load_bioclim(self, survey_id):\n        \"\"\"Load bioclimatic cube data\"\"\"\n        try:\n            path = os.path.join(self.bioclim_path, f\"GLC25-PA-{'train' if self.is_train else 'test'}-bioclimatic_monthly_{survey_id}_cube.pt\")\n            data = torch.load(path, weights_only=True)\n            # Normalize\n            data = (data - data.mean()) / (data.std() + 1e-6)\n            return data.float()\n        except:\n            return torch.zeros(4, 19, 12, dtype=torch.float32)\n    \n    def load_sentinel(self, survey_id):\n        \"\"\"Load Sentinel-2 patch data\"\"\"\n        try:\n            tiff_path = construct_patch_path(self.sentinel_path, survey_id)\n            with rasterio.open(tiff_path) as dataset:\n                image = dataset.read(out_dtype=np.float32)\n                image = np.array([quantile_normalize(band) for band in image])\n            return torch.from_numpy(image).float()\n        except:\n            return torch.zeros(4, 64, 64, dtype=torch.float32)\n    \n    def load_landsat(self, survey_id):\n        \"\"\"Load Landsat time series data\"\"\"\n        try:\n            path = os.path.join(self.landsat_path, f\"GLC25-PA-{'train' if self.is_train else 'test'}-landsat{'_time_series' if not self.is_train else '-time-series'}_{survey_id}_cube.pt\")\n            data = torch.nan_to_num(torch.load(path, weights_only=True))\n            # Normalize\n            data = (data - data.mean()) / (data.std() + 1e-6)\n            return data.float()\n        except:\n            return torch.zeros(6, 4, 21, dtype=torch.float32)\n    \n    def __getitem__(self, idx):\n        survey_id = self.metadata.surveyId.iloc[idx]\n        \n        # Load modalities\n        features = {}\n        if self.config.use_bioclim:\n            features['bioclim'] = self.load_bioclim(survey_id)\n        if self.config.use_sentinel:\n            features['sentinel'] = self.load_sentinel(survey_id)\n        if self.config.use_landsat:\n            features['landsat'] = self.load_landsat(survey_id)\n        \n        if self.is_train:\n            # Create multi-hot label\n            species_ids = self.label_dict.get(survey_id, [])\n            label = torch.zeros(self.config.num_classes)\n            for species_id in species_ids:\n                label[species_id] = 1\n            return features, label, survey_id\n        else:\n            return features, survey_id\n\n# ========================= Model Architecture =========================\nclass SpatialAttention(nn.Module):\n    \"\"\"Spatial attention module\"\"\"\n    def __init__(self, kernel_size=7):\n        super().__init__()\n        self.conv = nn.Conv2d(2, 1, kernel_size, padding=kernel_size//2)\n        self.sigmoid = nn.Sigmoid()\n    \n    def forward(self, x):\n        avg_out = torch.mean(x, dim=1, keepdim=True)\n        max_out, _ = torch.max(x, dim=1, keepdim=True)\n        x_cat = torch.cat([avg_out, max_out], dim=1)\n        x_out = self.conv(x_cat)\n        return self.sigmoid(x_out) * x\n\nclass ImprovedBlock(nn.Module):\n    \"\"\"Improved residual block with SE attention\"\"\"\n    def __init__(self, in_c, out_c, stride=1):\n        super().__init__()\n        self.conv1 = nn.Conv2d(in_c, out_c, 3, stride, 1, bias=False)\n        self.bn1 = nn.BatchNorm2d(out_c)\n        self.conv2 = nn.Conv2d(out_c, out_c, 3, 1, 1, bias=False)\n        self.bn2 = nn.BatchNorm2d(out_c)\n        \n        # SE module\n        self.se = nn.Sequential(\n            nn.AdaptiveAvgPool2d(1),\n            nn.Conv2d(out_c, out_c//4, 1),\n            nn.ReLU(inplace=True),\n            nn.Conv2d(out_c//4, out_c, 1),\n            nn.Sigmoid()\n        )\n        \n        self.downsample = None\n        if stride != 1 or in_c != out_c:\n            self.downsample = nn.Sequential(\n                nn.Conv2d(in_c, out_c, 1, stride, bias=False),\n                nn.BatchNorm2d(out_c)\n            )\n    \n    def forward(self, x):\n        identity = x\n        \n        out = F.relu(self.bn1(self.conv1(x)), inplace=True)\n        out = self.bn2(self.conv2(out))\n        \n        # Apply SE\n        se_weight = self.se(out)\n        out = out * se_weight\n        \n        if self.downsample:\n            identity = self.downsample(x)\n        \n        out += identity\n        return F.relu(out, inplace=True)\n\nclass ModalityEncoder(nn.Module):\n    \"\"\"Encoder for each modality\"\"\"\n    def __init__(self, in_channels, input_size, hidden_dim):\n        super().__init__()\n        \n        # Adaptive input normalization\n        self.norm = nn.LayerNorm([in_channels] + list(input_size))\n        \n        # Initial projection\n        self.stem = nn.Sequential(\n            nn.Conv2d(in_channels, 64, 3, 1, 1, bias=False),\n            nn.BatchNorm2d(64),\n            nn.ReLU(inplace=True)\n        )\n        \n        # Residual blocks with progressive channel expansion\n        self.blocks = nn.ModuleList([\n            ImprovedBlock(64, 128, stride=2 if min(input_size) > 4 else 1),\n            ImprovedBlock(128, 256, stride=2 if min(input_size) > 8 else 1),\n            ImprovedBlock(256, hidden_dim)\n        ])\n        \n        self.gap = nn.AdaptiveAvgPool2d(1)\n    \n    def forward(self, x):\n        x = self.norm(x)\n        x = self.stem(x)\n        \n        for block in self.blocks:\n            x = block(x)\n        \n        x = self.gap(x)\n        return x.flatten(1)\n\nclass MultiModalFusionNet(nn.Module):\n    \"\"\"Multi-modal fusion network with attention\"\"\"\n    def __init__(self, config):\n        super().__init__()\n        self.config = config\n        \n        # Modality-specific encoders\n        self.encoders = nn.ModuleDict()\n        if config.use_bioclim:\n            self.encoders['bioclim'] = ModalityEncoder(4, (19, 12), config.hidden_dim)\n        if config.use_sentinel:\n            self.encoders['sentinel'] = ModalityEncoder(4, (64, 64), config.hidden_dim)\n        if config.use_landsat:\n            self.encoders['landsat'] = ModalityEncoder(6, (4, 21), config.hidden_dim)\n        \n        # Cross-modal attention\n        num_modalities = len(self.encoders)\n        self.cross_attention = nn.MultiheadAttention(\n            config.hidden_dim, \n            num_heads=8, \n            dropout=config.dropout,\n            batch_first=True\n        )\n        \n        # Fusion and classification head\n        self.fusion = nn.Sequential(\n            nn.LayerNorm(config.hidden_dim * num_modalities),\n            nn.Linear(config.hidden_dim * num_modalities, config.hidden_dim),\n            nn.ReLU(inplace=True),\n            nn.Dropout(config.dropout),\n            nn.Linear(config.hidden_dim, config.hidden_dim // 2),\n            nn.ReLU(inplace=True),\n            nn.Dropout(config.dropout),\n            nn.Linear(config.hidden_dim // 2, config.num_classes)\n        )\n    \n    def forward(self, x):\n        # Encode each modality\n        encoded = []\n        for name, encoder in self.encoders.items():\n            if name in x:\n                feat = encoder(x[name])\n                encoded.append(feat)\n        \n        if len(encoded) > 1:\n            # Stack for attention [B, num_modalities, hidden_dim]\n            stacked = torch.stack(encoded, dim=1)\n            \n            # Apply cross-modal attention\n            attended, _ = self.cross_attention(stacked, stacked, stacked)\n            \n            # Concatenate all features\n            fused = torch.cat([attended[:, i] for i in range(attended.size(1))], dim=1)\n        else:\n            # Single modality\n            fused = encoded[0].repeat(1, len(self.encoders))  # Fix dimension mismatch\n        \n        return self.fusion(fused)\n\n# ========================= Loss Functions =========================\nclass FocalBCELoss(nn.Module):\n    \"\"\"Focal loss for addressing class imbalance\"\"\"\n    def __init__(self, alpha=0.25, gamma=2.0):\n        super().__init__()\n        self.alpha = alpha\n        self.gamma = gamma\n    \n    def forward(self, inputs, targets):\n        bce_loss = F.binary_cross_entropy_with_logits(inputs, targets, reduction='none')\n        pt = torch.exp(-bce_loss)\n        focal_loss = self.alpha * (1 - pt) ** self.gamma * bce_loss\n        return focal_loss.mean()\n\nclass AsymmetricLoss(nn.Module):\n    \"\"\"Asymmetric loss for multi-label classification\"\"\"\n    def __init__(self, gamma_neg=4, gamma_pos=1, clip=0.05, eps=1e-8):\n        super().__init__()\n        self.gamma_neg = gamma_neg\n        self.gamma_pos = gamma_pos\n        self.clip = clip\n        self.eps = eps\n    \n    def forward(self, x, y):\n        # Sigmoid activation\n        x_sigmoid = torch.sigmoid(x)\n        xs_pos = x_sigmoid\n        xs_neg = 1 - x_sigmoid\n        \n        # Asymmetric clipping\n        if self.clip is not None and self.clip > 0:\n            xs_neg = (xs_neg + self.clip).clamp(max=1)\n        \n        # Basic CE calculation\n        los_pos = y * torch.log(xs_pos.clamp(min=self.eps))\n        los_neg = (1 - y) * torch.log(xs_neg.clamp(min=self.eps))\n        \n        # Asymmetric focusing\n        pt0 = xs_pos * y\n        pt1 = xs_neg * (1 - y)\n        pt = pt0 + pt1\n        one_sided_gamma = self.gamma_pos * y + self.gamma_neg * (1 - y)\n        one_sided_w = torch.pow(1 - pt, one_sided_gamma)\n        \n        loss = -torch.sum(one_sided_w * (los_pos + los_neg))\n        return loss / x.shape[0]\n\n# ========================= Training Functions =========================\ndef mixup_data(x, y, alpha=1.0):\n    \"\"\"Mixup augmentation for multi-modal data\"\"\"\n    if alpha > 0:\n        lam = np.random.beta(alpha, alpha)\n    else:\n        lam = 1\n    \n    batch_size = list(x.values())[0].size(0)\n    index = torch.randperm(batch_size).to(list(x.values())[0].device)\n    \n    mixed_x = {}\n    for key in x:\n        mixed_x[key] = lam * x[key] + (1 - lam) * x[key][index]\n    \n    y_a, y_b = y, y[index]\n    return mixed_x, y_a, y_b, lam\n\ndef train_epoch(model, dataloader, criterion, optimizer, config, use_mixup=True):\n    \"\"\"Train for one epoch\"\"\"\n    model.train()\n    total_loss = 0\n    \n    pbar = tqdm(dataloader, desc=\"Training\")\n    for batch_idx, (data, targets, _) in enumerate(pbar):\n        # Move to device\n        for key in data:\n            data[key] = data[key].to(config.device)\n        targets = targets.to(config.device)\n        \n        # Apply mixup\n        if use_mixup and np.random.random() > 0.5:\n            data, targets_a, targets_b, lam = mixup_data(data, targets, config.mixup_alpha)\n            outputs = model(data)\n            loss = lam * criterion(outputs, targets_a) + (1 - lam) * criterion(outputs, targets_b)\n        else:\n            outputs = model(data)\n            loss = criterion(outputs, targets)\n        \n        # Backward pass\n        optimizer.zero_grad()\n        loss.backward()\n        \n        # Gradient clipping\n        torch.nn.utils.clip_grad_norm_(model.parameters(), config.gradient_clip)\n        \n        optimizer.step()\n        \n        total_loss += loss.item()\n        pbar.set_postfix({'loss': f'{loss.item():.4f}', 'avg_loss': f'{total_loss/(batch_idx+1):.4f}'})\n    \n    return total_loss / len(dataloader)\n\ndef validate(model, dataloader, criterion, config):\n    \"\"\"Validation function\"\"\"\n    model.eval()\n    total_loss = 0\n    all_preds = []\n    all_targets = []\n    \n    with torch.no_grad():\n        for data, targets, _ in tqdm(dataloader, desc=\"Validating\"):\n            for key in data:\n                data[key] = data[key].to(config.device)\n            targets = targets.to(config.device)\n            \n            outputs = model(data)\n            loss = criterion(outputs, targets)\n            \n            total_loss += loss.item()\n            \n            # Store predictions for F1 calculation\n            preds = torch.sigmoid(outputs)\n            all_preds.append(preds.cpu())\n            all_targets.append(targets.cpu())\n    \n    # Calculate F1 score\n    all_preds = torch.cat(all_preds)\n    all_targets = torch.cat(all_targets)\n    \n    # Get top-25 predictions for each sample\n    top_k = 25\n    pred_binary = torch.zeros_like(all_preds)\n    for i in range(all_preds.size(0)):\n        top_indices = torch.topk(all_preds[i], top_k).indices\n        pred_binary[i, top_indices] = 1\n    \n    f1 = f1_score(all_targets.numpy(), pred_binary.numpy(), average='samples', zero_division=0)\n    \n    return total_loss / len(dataloader), f1\n\ndef predict(model, dataloader, config, top_k=25):\n    \"\"\"Generate predictions\"\"\"\n    model.eval()\n    all_predictions = []\n    all_survey_ids = []\n    \n    with torch.no_grad():\n        for batch in tqdm(dataloader, desc=\"Predicting\"):\n            if len(batch) == 2:\n                data, survey_ids = batch\n            else:\n                data, _, survey_ids = batch\n            \n            for key in data:\n                data[key] = data[key].to(config.device)\n            \n            outputs = model(data)\n            probs = torch.sigmoid(outputs).cpu().numpy()\n            \n            # Get top-k predictions\n            top_k_indices = np.argsort(-probs, axis=1)[:, :top_k]\n            all_predictions.extend(top_k_indices)\n            all_survey_ids.extend(survey_ids.numpy())\n    \n    return all_predictions, all_survey_ids\n\n# ========================= Main Training Loop =========================\ndef main():\n    # Initialize configuration\n    config = Config()\n    set_seed(config.seed)\n    \n    print(f\"Device: {config.device}\")\n    print(f\"Using modalities: Bioclim={config.use_bioclim}, Sentinel={config.use_sentinel}, Landsat={config.use_landsat}\")\n    \n    # Load metadata\n    train_metadata = pd.read_csv(config.train_metadata_path)\n    test_metadata = pd.read_csv(config.test_metadata_path)\n    \n    # Create train/val split\n    n_train = int(len(train_metadata) * (1 - config.val_split))\n    train_meta = train_metadata.iloc[:n_train]\n    val_meta = train_metadata.iloc[n_train:]\n    \n    # Create datasets\n    train_dataset = MultiModalDataset(config, train_meta, is_train=True)\n    val_dataset = MultiModalDataset(config, val_meta, is_train=True)\n    test_dataset = MultiModalDataset(config, test_metadata, is_train=False)\n    \n    # Create dataloaders\n    train_loader = DataLoader(train_dataset, batch_size=config.batch_size, shuffle=True, num_workers=4, pin_memory=True)\n    val_loader = DataLoader(val_dataset, batch_size=config.batch_size*2, shuffle=False, num_workers=4, pin_memory=True)\n    test_loader = DataLoader(test_dataset, batch_size=config.batch_size*2, shuffle=False, num_workers=4, pin_memory=True)\n    \n    # Initialize model\n    model = MultiModalFusionNet(config).to(config.device)\n    \n    # Loss function (use AsymmetricLoss for better multi-label performance)\n    criterion = AsymmetricLoss()\n    \n    # Optimizer and scheduler\n    optimizer = torch.optim.AdamW(\n        model.parameters(), \n        lr=config.learning_rate, \n        weight_decay=config.weight_decay\n    )\n    \n    scheduler = CosineAnnealingWarmRestarts(\n        optimizer, \n        T_0=config.warmup_epochs, \n        T_mult=2, \n        eta_min=1e-6\n    )\n    \n    # Training loop\n    best_f1 = 0\n    for epoch in range(config.num_epochs):\n        print(f\"\\n{'='*50}\")\n        print(f\"Epoch {epoch+1}/{config.num_epochs}\")\n        print(f\"LR: {scheduler.get_last_lr()[0]:.6f}\")\n        \n        # Train\n        train_loss = train_epoch(model, train_loader, criterion, optimizer, config, use_mixup=epoch > config.warmup_epochs)\n        \n        # Validate\n        val_loss, val_f1 = validate(model, val_loader, criterion, config)\n        \n        print(f\"Train Loss: {train_loss:.4f}\")\n        print(f\"Val Loss: {val_loss:.4f}, Val F1: {val_f1:.4f}\")\n        \n        # Step scheduler\n        scheduler.step()\n        \n        # Save best model\n        if val_f1 > best_f1:\n            best_f1 = val_f1\n            torch.save(model.state_dict(), 'best_model.pth')\n            print(f\"New best model saved! F1: {best_f1:.4f}\")\n    \n    # Load best model for prediction\n    model.load_state_dict(torch.load('best_model.pth'))\n    \n    # Generate predictions\n    predictions, survey_ids = predict(model, test_loader, config)\n    \n    # Create submission\n    submission_data = []\n    for survey_id, preds in zip(survey_ids, predictions):\n        pred_str = ' '.join(map(str, preds))\n        submission_data.append({'surveyId': survey_id, 'predictions': pred_str})\n    \n    submission_df = pd.DataFrame(submission_data)\n    submission_df.to_csv('submission.csv', index=False)\n    print(\"\\nSubmission saved to submission.csv\")\n    print(f\"Best validation F1: {best_f1:.4f}\")\n\nif __name__ == \"__main__\":\n    main()","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}