{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"nvidiaTeslaT4","dataSources":[{"sourceId":81000,"databundleVersionId":8812083,"sourceType":"competition"}],"dockerImageVersionId":30734,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"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\nimport torch\nfrom torch import nn\nfrom torch.utils.data import DataLoader, Dataset\nfrom tqdm.notebook import tqdm\nfrom sklearn.preprocessing import StandardScaler\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},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"DATA_DIR = r'/kaggle/input/the-future-crop-challenge'\ndevice = torch.device(\"cuda:0\" if torch.cuda.is_available() else \"cpu\")\nprint(f'Running on {device}')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class YieldNet(nn.Module):\n    \"\"\"\n    Defines yield neural network. A LSTM processes the daily climate feature timeseries, while a parallel branch \n    handles fixed soil properties per location and year in a fully connected network (FCN). \n    The resuls of both branches are merged to a single yield prediction (again, per location and year).\n    \"\"\"\n    def __init__(self, seq_input_size, seq_hidden_size, fcn_layers):\n        super().__init__()\n\n        self.seq_input_size = seq_input_size\n        self.seq_hidden_size = seq_hidden_size\n\n        self.lstm = nn.LSTM(input_size=seq_input_size, hidden_size=seq_hidden_size, num_layers=1, batch_first=True)\n        self.seq_fcn = nn.Linear(self.seq_hidden_size, 1)\n\n        self.fcn_layers = fcn_layers\n        self.activation = nn.Tanh()\n        self.linears = nn.ModuleList([nn.Linear(fcn_layers[i], fcn_layers[i+1]) for i in range(len(fcn_layers)-1)]) \n\n        for i in range(len(fcn_layers)-1):         \n            nn.init.xavier_normal_(self.linears[i].weight.data, gain=1.0)            \n            nn.init.zeros_(self.linears[i].bias.data)   \n    \n    def forward(self, x_seq, x_fcn):\n        h0 = torch.autograd.Variable(torch.zeros(1, x_seq.shape[0], self.seq_hidden_size)).to(x_seq.device)\n        c0 = torch.autograd.Variable(torch.zeros(1, x_seq.shape[0], self.seq_hidden_size)).to(x_seq.device)\n\n        out, _ = self.lstm(x_seq, (h0, c0))\n        out = self.seq_fcn(out[:, -1, :])\n\n        a = x_fcn\n        for i in range(len(self.fcn_layers) - 2):  \n            z = self.linears[i](a)\n            a = self.activation(z)\n        a = self.linears[-1](a)\n\n        return out.reshape(out.shape[0], out.shape[1]) * a","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class ClimateDataset(Dataset):\n    \"\"\"\n    The ClimateDataset class provides a convenient way for acessing, merging and scaling parquet data. \n    This is used in the main training loop to access features and target variables.\n    \"\"\"\n    def __init__(self, crop: str, mode: str, data_dir: str, scalers: list = None):\n        self.tasmax = pd.read_parquet(os.path.join(data_dir, f\"tasmax_{crop}_{mode}.parquet\"))\n        self.tasmin = pd.read_parquet(os.path.join(data_dir, f\"tasmin_{crop}_{mode}.parquet\"))\n        # self.tas = pd.read_parquet(os.path.join(data_dir, f\"tas_{crop}_{mode}.parquet\"))\n        self.pr = pd.read_parquet(os.path.join(data_dir, f\"pr_{crop}_{mode}.parquet\"))\n        self.rsds = pd.read_parquet(os.path.join(data_dir, f\"rsds_{crop}_{mode}.parquet\"))\n        self.soil_co2 = pd.read_parquet(os.path.join(data_dir, f\"soil_co2_{crop}_{mode}.parquet\"))\n        if mode == 'train':\n            self.yield_ = pd.read_parquet(os.path.join(data_dir, f\"{mode}_solutions_{crop}.parquet\"))\n        else:\n            self.yield_ = None\n\n        if scalers is None:\n            self._init_scalers()\n        else:\n            self.scaler_climate, self.scaler_soil, self.scaler_yield = scalers\n        self._check_data([self.tasmax, self.tasmin, self.pr, self.rsds, self.soil_co2], self.yield_)\n\n    def __getitem__(self, index):\n        # 240x4 climate matrix per location/year (features in last dimension by convention)\n        climate = np.vstack([\n            self.tasmax.iloc[index, 5:].astype(np.float32), \n            self.tasmin.iloc[index, 5:].astype(np.float32),\n            self.pr.iloc[index, 5:].astype(np.float32),\n            self.rsds.iloc[index, 5:].astype(np.float32),\n        ]).T\n\n        # Fixed soil properties per location/year\n        soil = self.soil_co2.iloc[index][['co2', 'nitrogen']].astype(np.float32)\n        id = soil.name\n        soil = soil.values\n\n        # Yield estimated by process model\n        if self.yield_ is not None:\n            yield_ = self.yield_.iloc[index].astype(np.float32).values\n            yield_ = self.scaler_yield.transform(yield_.reshape(1, -1)).reshape(-1)\n        else:\n            yield_ = None\n\n        climate = self.scaler_climate.transform(climate)\n        soil = self.scaler_soil.transform(soil.reshape(1, -1)).reshape(-1)\n\n        return torch.tensor(climate), torch.tensor(soil), torch.tensor(yield_ or []), id\n\n    def __len__(self):\n        return self.tasmax.shape[0]\n\n    def _init_scalers(self):\n        # Draw random sample from climate data to estimate distribution moments for scaler.\n        climate_sample = np.vstack([\n            self.tasmax.sample(1000).iloc[:, 5:].values.flatten(), \n            self.tasmin.sample(1000).iloc[:, 5:].values.flatten(),\n            self.pr.sample(1000).iloc[:, 5:].values.flatten(),\n            self.rsds.sample(1000).iloc[:, 5:].values.flatten(),\n        ]).T\n        self.scaler_climate = StandardScaler()\n        self.scaler_climate.fit(climate_sample)\n\n        # Scaler for fixed soil properties\n        self.scaler_soil = StandardScaler()\n        self.scaler_soil.fit(self.soil_co2[['co2', 'nitrogen']].values)\n\n        # Scaler for yield\n        self.scaler_yield = StandardScaler()\n        if self.yield_ is not None:\n            self.scaler_yield.fit(self.yield_.values)\n\n    def _check_data(self, climate: list, yield_: pd.DataFrame) -> bool:\n        # Check for matching year, lon, lat columns\n        for i in range(1, len(climate)):\n            assert np.all(climate[0][['year', 'lon', 'lat']] == climate[i][['year', 'lon', 'lat']])\n        # Check label for matching length\n        assert yield_ is None or climate[0].shape[0] == yield_.shape[0]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Intialize model, data loader, loss and optimizer\nds_train = ClimateDataset('maize', 'train', data_dir=DATA_DIR)\ntrain_loader = DataLoader(ds_train, batch_size=100, shuffle=True)\n\nmodel = YieldNet(seq_input_size=4, seq_hidden_size=128, fcn_layers=[2, 32, 32, 1]).to(device)\ncost_mse = torch.nn.MSELoss()\noptimizer = torch.optim.Adam(model.parameters(), lr=0.001)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Main training loop\n# In this example the model is only trained on maize data!\nloss = []\nnum_epochs = 1\npbar = tqdm(range(num_epochs))\n\nfor epoch in pbar:\n\n    for i, (climate, soil, yield_true, _), in enumerate(tqdm(train_loader, leave=False)):\n        optimizer.zero_grad()\n        yield_pred = model(climate.to(device), soil.to(device))\n        data_loss = cost_mse(yield_pred, yield_true.to(device))\n        data_loss.backward()\n        optimizer.step()\n\n        if i % 10 == 0:\n            loss_num = data_loss.item()\n            loss.append(loss_num)\n            pbar.set_description(f\"Loss: {loss_num:5.10f}\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Predict test set for submission\nds_test_maize = ClimateDataset('maize', 'test', data_dir=DATA_DIR, scalers=(ds_train.scaler_climate, ds_train.scaler_soil, ds_train.scaler_yield))\nds_test_wheat = ClimateDataset('wheat', 'test', data_dir=DATA_DIR, scalers=(ds_train.scaler_climate, ds_train.scaler_soil, ds_train.scaler_yield))\ntest_loader_maize = DataLoader(ds_test_maize, batch_size=100, shuffle=False)\ntest_loader_wheat = DataLoader(ds_test_wheat, batch_size=100, shuffle=False)\n\nids = []\nyields_pred = []\nfor climate, soil, _, id in tqdm(test_loader_maize):\n    ids.append(id.detach().numpy())\n    yields_pred.append(model(climate.to(device), soil.to(device)).detach().cpu().numpy())\n\nfor climate, soil, _, id in tqdm(test_loader_wheat):\n    ids.append(id.detach().numpy())\n    yields_pred.append(model(climate.to(device), soil.to(device)).detach().cpu().numpy())\n\nyields_pred = np.concatenate(yields_pred)\nyields_pred = ds_train.scaler_yield.inverse_transform(yields_pred)\nyields_pred = yields_pred.reshape(-1)\n\npredictions = pd.Series(yields_pred, index=np.concatenate(ids))\npredictions.index.name = 'ID'\npredictions.name = 'yield'\npredictions.to_csv('submission.csv')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}