{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.14","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"gpu","dataSources":[{"sourceId":81000,"databundleVersionId":8812083,"sourceType":"competition"}],"dockerImageVersionId":30761,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport pyarrow\nimport os\nfrom os import listdir\nfrom os.path import isfile, join\nfrom time import time\nimport re\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.metrics import r2_score\n\nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\nimport torch.optim as optim\n# from torch.utils.data import Dataset, DataLoader, TensorDataset\nfrom torch.utils.data import TensorDataset, DataLoader, random_split,DataLoader, Dataset\nfrom torch.nn.utils.rnn import pack_padded_sequence, pad_packed_sequence","metadata":{"execution":{"iopub.status.busy":"2024-09-07T14:27:53.217285Z","iopub.execute_input":"2024-09-07T14:27:53.217747Z","iopub.status.idle":"2024-09-07T14:27:56.147465Z","shell.execute_reply.started":"2024-09-07T14:27:53.217695Z","shell.execute_reply":"2024-09-07T14:27:56.146157Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mypath = '/kaggle/input/the-future-crop-challenge/'","metadata":{"execution":{"iopub.status.busy":"2024-09-07T14:28:03.644816Z","iopub.execute_input":"2024-09-07T14:28:03.645688Z","iopub.status.idle":"2024-09-07T14:28:03.650415Z","shell.execute_reply.started":"2024-09-07T14:28:03.645643Z","shell.execute_reply":"2024-09-07T14:28:03.649240Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"USE_CUDA = torch.cuda.is_available()\n# DEVICE = torch.device('cuda:0')\nDEVICE = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")","metadata":{"execution":{"iopub.status.busy":"2024-09-07T14:28:28.105284Z","iopub.execute_input":"2024-09-07T14:28:28.105836Z","iopub.status.idle":"2024-09-07T14:28:28.114184Z","shell.execute_reply.started":"2024-09-07T14:28:28.105781Z","shell.execute_reply":"2024-09-07T14:28:28.112845Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load train","metadata":{}},{"cell_type":"code","source":"def get_region(row):\n    if -125 <= row['lon'] <= -66 and 24 <= row['lat'] <= 50:\n        return 'USA'\n    elif -75 <= row['lon'] <= -35 and -40 <= row['lat'] <= -10:\n        return 'South America'\n    elif -10 <= row['lon'] <= 50 and 35 <= row['lat'] <= 60:\n        return 'Europe'\n    else:\n        return 'Other'","metadata":{"execution":{"iopub.status.busy":"2024-09-07T11:02:46.873774Z","iopub.execute_input":"2024-09-07T11:02:46.874146Z","iopub.status.idle":"2024-09-07T11:02:46.884115Z","shell.execute_reply.started":"2024-09-07T11:02:46.874100Z","shell.execute_reply":"2024-09-07T11:02:46.882910Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def load(crop,mode,data_dir):\n    tasmax = pd.read_parquet(os.path.join(data_dir, \"tasmax_{}_{}.parquet\".format(crop,mode)))\n    tasmin = pd.read_parquet(os.path.join(data_dir, \"tasmin_{}_{}.parquet\".format(crop,mode)))\n    tas = pd.read_parquet(os.path.join(data_dir, \"tas_{}_{}.parquet\".format(crop,mode)))\n    pr = pd.read_parquet(os.path.join(data_dir, \"pr_{}_{}.parquet\".format(crop,mode)))\n    rsds = pd.read_parquet(os.path.join(data_dir, \"rsds_{}_{}.parquet\".format(crop,mode)))\n    soil_co2 = pd.read_parquet(os.path.join(data_dir, \"soil_co2_{}_{}.parquet\".format(crop,mode)))\n\n    if mode == 'train':\n        yield_ = pd.read_parquet(os.path.join(data_dir, \"{}_solutions_{}.parquet\".format(mode,crop)))\n        yield_ = yield_.values.astype(np.float32)\n\n    if mode == 'test':\n        yield_ = None\n        \n    climate = np.concatenate([\n        tas.values[:, 5:,np.newaxis].astype(np.float32),\n        tasmax.values[:, 5:,np.newaxis].astype(np.float32), \n        tasmin.values[:, 5:,np.newaxis].astype(np.float32),\n        pr.values[:, 5:,np.newaxis].astype(np.float32),\n        rsds.values[:, 5:,np.newaxis].astype(np.float32),\n    ], axis=2)    \n    return(climate,yield_,soil_co2)","metadata":{"execution":{"iopub.status.busy":"2024-09-07T11:02:46.885662Z","iopub.execute_input":"2024-09-07T11:02:46.886133Z","iopub.status.idle":"2024-09-07T11:02:46.899922Z","shell.execute_reply.started":"2024-09-07T11:02:46.886082Z","shell.execute_reply":"2024-09-07T11:02:46.898767Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_training_dataloader(crop,mode,mypath,region,static_var,detrend=True):\n    (climate, yield_label, soil) = load(crop,mode,mypath)\n\n    # Reshape and calculate mean and sd along the grouped axis\n    climate_mean = climate.reshape(climate.shape[0], 8, 30, 5).mean(axis=2)\n    climate_sd = climate.reshape(climate.shape[0], 8, 30, 5).std(axis=2)\n\n    GPP = np.maximum(climate.reshape(climate.shape[0], 8, 30, 5)[:,:,:,3] - 10, 0).sum(axis=2,keepdims=True).cumsum(axis=1) # cumulative GDD by phases\n    heat_stress_day =  (climate.reshape(climate.shape[0], 8, 30, 5)[:,:,:,1] > 30).sum(axis=2,keepdims=True).cumsum(axis=1)\n    frost_days = (climate.reshape(climate.shape[0], 8, 30, 5)[:,:,:,2] < 0).sum(axis=2,keepdims=True).cumsum(axis=1)\n\n    climate_stats = np.concatenate([climate_sd,climate_mean,GPP,heat_stress_day,frost_days],axis=2)\n    \n    #annual mean\n    annual_df = pd.DataFrame(climate.mean(axis=1),columns=['tas', 'tasmax', 'tasmin', 'pr', 'rsds'])\n    soil = soil.reset_index()\n    soil = pd.concat([annual_df, soil], axis=1)\n    soil['region'] = soil.apply(get_region, axis=1)\n    soil['co2'] = soil['co2']/1500\n    \n    #detrend yield\n    soil['yield'] = yield_label\n    grouped_mean = soil.groupby(['lon', 'lat'])['yield'].mean().reset_index()\n    grouped_mean = grouped_mean.rename(columns={'yield': 'yield_mean'})\n    soil = pd.merge(soil, grouped_mean, on=['lon', 'lat'])\n    soil['de_yield'] = soil['yield']-soil['yield_mean']\n    \n    #scale input features\n    climate_arr = climate_stats.reshape(-1, climate_stats.shape[2])\n    scaler = StandardScaler()\n    scaler.fit(climate_arr)\n    scaled_climate = scaler.transform(climate_arr).reshape(climate_stats.shape)\n    \n    #scale static features\n    scaler_static = StandardScaler()\n    static_arr = soil[static_var].values.astype(np.float32)\n    scaler_static.fit(static_arr)\n    scaled_static = scaler_static.transform(static_arr)\n    #add co2\n    scaled_static = np.concatenate([scaled_static,soil['co2'].values.reshape(-1, 1)],axis = 1)\n\n    X_train_to_loader = torch.from_numpy(scaled_climate)\n    X_train_to_loader= X_train_to_loader.type(torch.FloatTensor)\n    X_static_to_loader = torch.from_numpy(scaled_static).type(torch.FloatTensor)\n    if detrend:\n        y_train_to_loader = torch.from_numpy(soil['de_yield'].values.reshape(-1, 1))\n    else:\n        y_train_to_loader = torch.from_numpy(soil['yield'].values.reshape(-1, 1))\n        \n    dataset = TensorDataset(X_train_to_loader, X_static_to_loader, y_train_to_loader)\n\n    # subset for REGION\n    Region_index = soil[soil['region'].isin(region)].index.values\n    dataset = torch.utils.data.Subset(dataset, Region_index)\n    print(len(dataset))\n    \n    # Define the sizes of train and validation sets\n    train_size = int(0.8 * len(dataset))  # 80% of the data for training\n    val_size = int(0.1 *len(dataset))  # Remaining 10% for validation\n    test_size = len(dataset)- train_size -val_size# Remaining 10% for test\n    \n    # Split dataset into train and validation sets\n    train_dataset, val_dataset,test_dataset = random_split(dataset, [train_size, val_size,test_size])\n    \n    train_loader = torch.utils.data.DataLoader(train_dataset, 3000, shuffle=True)   \n    val_loader = torch.utils.data.DataLoader(val_dataset, 3000, shuffle=True)  \n    test_loader = torch.utils.data.DataLoader(test_dataset, 3000, shuffle=True)  \n    return train_loader,val_loader,test_loader","metadata":{"execution":{"iopub.status.busy":"2024-09-07T11:02:46.901339Z","iopub.execute_input":"2024-09-07T11:02:46.901761Z","iopub.status.idle":"2024-09-07T11:02:47.057609Z","shell.execute_reply.started":"2024-09-07T11:02:46.901721Z","shell.execute_reply":"2024-09-07T11:02:47.056100Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Train","metadata":{}},{"cell_type":"code","source":"class LSTMForecaster_static(nn.Module):\n  def __init__(self, n_features, n_hidden, n_outputs,n_static, n_static_hidden, sequence_len,DEVICE, n_lstm_layers=1, n_deep_layers=10, use_cuda=True, dropout=0.2):\n    '''\n    n_features: number of input features (1 for univariate forecasting)\n    n_hidden: number of neurons in each hidden layer\n    n_outputs: number of outputs to predict for each training example\n    n_static: number of static features\n    n_deep_layers: number of hidden dense layers after the lstm layer\n    sequence_len: number of steps to look back at for prediction\n    dropout: float (0 < dropout < 1) dropout ratio between dense layers\n    '''\n    super().__init__()\n    self.n_lstm_layers = n_lstm_layers\n    self.nhid = n_hidden\n    self.use_cuda = use_cuda # set option for device selection\n    self.DEVICE = DEVICE # set option for device selection\n    \n    # LSTM Layer\n    self.lstm = nn.LSTM(n_features,\n                        n_hidden,\n                        num_layers=n_lstm_layers,\n                        batch_first=True) # As we have transformed our data in this way\n\n      \n    # First dense layer to expand static variables\n    self.fc1 = nn.Linear(n_static,n_static_hidden)\n    self.relu1 = nn.ReLU()\n      \n    # Dropout layer \n    self.dropout = nn.Dropout(p=dropout)\n    # Create fully connected layers (n_hidden x n_deep_layers)\n    dnn_layers = []\n      \n    # the first layer to concatenate hidden and static properties, 2 additional elements\n    dnn_layers.append(nn.ReLU())\n    dnn_layers.append(nn.Linear(n_hidden*n_lstm_layers+n_static_hidden, n_hidden)) \n      \n    for i in range(n_deep_layers):\n      # Last layer (n_hidden x n_outputs)\n      if i == n_deep_layers - 1:\n        dnn_layers.append(nn.ReLU())\n        dnn_layers.append(nn.Linear(n_hidden, n_outputs))\n      # All other layers (n_hidden x n_hidden) with dropout option\n      else:\n        dnn_layers.append(nn.ReLU())\n        dnn_layers.append(nn.Linear(n_hidden, n_hidden))\n        if dropout:\n          dnn_layers.append(nn.Dropout(p=dropout))\n    # compile DNN layers\n    self.dnn = nn.Sequential(*dnn_layers)\n\n  def forward(self, x,x_static):\n\n    # Initialize hidden state\n    hidden_state = torch.zeros(self.n_lstm_layers, x.shape[0], self.nhid)\n    cell_state = torch.zeros(self.n_lstm_layers, x.shape[0], self.nhid)\n\n    # move hidden state to device\n    if self.use_cuda:\n      hidden_state = hidden_state.to(self.DEVICE)\n      cell_state = cell_state.to(self.DEVICE)\n        \n    self.hidden = (hidden_state, cell_state)\n    # Forward Pass\n    output, (h,c) = self.lstm(x, self.hidden) # LSTM\n    \n    # Flatten hidden state of the last step \n    x_ = self.dropout(h.contiguous().view(x.shape[0], -1)) \n    x_static = self.relu1(self.fc1(x_static))\n    # Pass forward hidden state and static features through fully connected DNN.  \n    return self.dnn(torch.cat((x_,x_static),axis=1)).squeeze() ","metadata":{"execution":{"iopub.status.busy":"2024-09-07T11:07:42.411833Z","iopub.execute_input":"2024-09-07T11:07:42.412290Z","iopub.status.idle":"2024-09-07T11:07:42.429466Z","shell.execute_reply.started":"2024-09-07T11:07:42.412246Z","shell.execute_reply":"2024-09-07T11:07:42.428205Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def train_model(model, train_loader, val_loader, test_loader,static_var,DEVICE,lr, weight_decay,step_size,sch_gamma,num_epochs=100):\n    # Initialize the loss function and optimizer\n    criterion = nn.MSELoss().to(DEVICE)\n    optimizer = torch.optim.AdamW(model.parameters(), lr=lr, weight_decay=weight_decay)\n    scheduler = torch.optim.lr_scheduler.StepLR(optimizer, step_size=step_size, gamma=sch_gamma)\n    #static & scaled\n    # num_epochs = 100\n    train_losses = [] #save performance metrics\n    train_r2s = []\n    val_losses = []\n    val_r2s = []\n    model.train()\n    for epoch in range(num_epochs):\n        start = time()\n        model.train()\n        for i, (x, x_static, y_true), in enumerate(train_loader):\n            # model.train()\n            # print(i)\n            xt,x_st, yt = x.to(DEVICE),x_static.to(DEVICE),y_true.to(DEVICE).squeeze()\n            \n            y_pred = model(xt,x_st)\n            loss_batch = criterion(y_pred, yt)\n            optimizer.zero_grad()\n            loss_batch.backward()\n            optimizer.step()\n        scheduler.step()\n    \n        if epoch % 10 != 0:\n            continue\n        train_loss,train_r2 = test_mse_r2_static(model, criterion, train_loader,DEVICE)\n        val_loss,val_r2 = test_mse_r2_static(model, criterion, val_loader,DEVICE)\n        test_loss,test_r2 = test_mse_r2_static(model, criterion, test_loader,DEVICE)\n        \n        train_losses.append(train_loss)\n        train_r2s.append(train_r2)\n        val_losses.append(val_loss)\n        val_r2s.append(val_r2)\n        end = time()\n        \n        print(\"Epoch %d: TRAIN MSE loss: %.6f, TRAIN r2 loss: %.6f\" % (epoch, train_loss,train_r2))\n        print(\"Epoch %d: VAL MSE loss: %.6f, VAL r2: %.6f\" % (epoch, val_loss, val_r2))\n        print(\"Epoch %d: TEST MSE loss: %.6f, TEST r2: %.6f\" % (epoch, test_loss, test_r2))\n        print(\"Finish with for one epoch:{} second\".format(end - start))\n    return model","metadata":{"execution":{"iopub.status.busy":"2024-09-07T11:07:43.250421Z","iopub.execute_input":"2024-09-07T11:07:43.251570Z","iopub.status.idle":"2024-09-07T11:07:43.265794Z","shell.execute_reply.started":"2024-09-07T11:07:43.251500Z","shell.execute_reply":"2024-09-07T11:07:43.264527Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def r2_loss(output, target):\n    target_mean = torch.mean(target)\n    ss_tot = torch.sum((target - target_mean) ** 2)\n    ss_res = torch.sum((target - output) ** 2)\n    r2 = 1 - ss_res / ss_tot\n    return r2","metadata":{"execution":{"iopub.status.busy":"2024-09-07T11:07:43.986280Z","iopub.execute_input":"2024-09-07T11:07:43.986743Z","iopub.status.idle":"2024-09-07T11:07:43.993280Z","shell.execute_reply.started":"2024-09-07T11:07:43.986703Z","shell.execute_reply":"2024-09-07T11:07:43.991981Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def test_mse_r2_static(model,loss_fn, dataloader,DEVICE):\n    model.eval()\n    loss_test_cum = 0\n    y_test_pred = []\n    y_test_label = []\n    with torch.no_grad():\n        for i, (x, z, y) in enumerate(dataloader):        \n            y_pred = model(x.to(DEVICE),z.to(DEVICE))\n            # predicted = model(input)\n            loss   = loss_fn(y_pred, y.to(DEVICE).squeeze())\n            loss_test_cum += loss\n            y_test_pred.append(y_pred)\n            y_test_label.append(y.to(DEVICE).squeeze())\n        r2 = r2_loss(torch.cat(y_test_pred),torch.cat(y_test_label))\n    return loss_test_cum/(i+1),r2","metadata":{"execution":{"iopub.status.busy":"2024-09-07T11:07:44.284277Z","iopub.execute_input":"2024-09-07T11:07:44.284729Z","iopub.status.idle":"2024-09-07T11:07:44.293543Z","shell.execute_reply.started":"2024-09-07T11:07:44.284689Z","shell.execute_reply":"2024-09-07T11:07:44.292167Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# (climate, yield_label, soil) = load('maize','train',mypath)\nstatic_var=['nitrogen','texture_class','lon','lat','yield_mean','tas', 'tasmax', 'tasmin', 'pr', 'rsds']","metadata":{"execution":{"iopub.status.busy":"2024-09-07T11:07:45.515436Z","iopub.execute_input":"2024-09-07T11:07:45.516484Z","iopub.status.idle":"2024-09-07T11:07:45.521693Z","shell.execute_reply.started":"2024-09-07T11:07:45.516435Z","shell.execute_reply":"2024-09-07T11:07:45.520466Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Maize","metadata":{}},{"cell_type":"code","source":"## make maize train_loader\ntrain_loader,val_loader,test_loader = get_training_dataloader('maize','train',mypath,region=['USA', 'Other', 'South America', 'Europe'],static_var=static_var)","metadata":{"execution":{"iopub.status.busy":"2024-09-07T11:02:47.138558Z","iopub.execute_input":"2024-09-07T11:02:47.138931Z","iopub.status.idle":"2024-09-07T11:05:09.661876Z","shell.execute_reply.started":"2024-09-07T11:02:47.138880Z","shell.execute_reply":"2024-09-07T11:05:09.659996Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"(x,z,y) = next(iter(train_loader))","metadata":{"execution":{"iopub.status.busy":"2024-09-07T11:05:09.664640Z","iopub.execute_input":"2024-09-07T11:05:09.665173Z","iopub.status.idle":"2024-09-07T11:05:09.802610Z","shell.execute_reply.started":"2024-09-07T11:05:09.665117Z","shell.execute_reply":"2024-09-07T11:05:09.801343Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# # global model training\n# model = LSTMForecaster_static(n_features = x.shape[2], n_hidden = 200, n_outputs = 1, \n#                               n_static = len(static_var)+1,n_static_hidden=200,sequence_len = x.shape[1], \n#                               DEVICE = DEVICE, n_deep_layers=4, use_cuda=USE_CUDA, dropout=0.4).to(DEVICE)\n\n# maize_unscale_model_GLOB2 = train_model(model, train_loader,val_loader,test_loader,static_var,DEVICE,\n#                                   lr=4e-4, weight_decay=0.1,step_size=100,sch_gamma=0.8,num_epochs=200)","metadata":{"execution":{"iopub.status.busy":"2024-09-07T11:05:09.803935Z","iopub.execute_input":"2024-09-07T11:05:09.804308Z","iopub.status.idle":"2024-09-07T11:05:09.809758Z","shell.execute_reply.started":"2024-09-07T11:05:09.804270Z","shell.execute_reply":"2024-09-07T11:05:09.808624Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# global model training\nmodel = LSTMForecaster_static(n_features = x.shape[2], n_hidden = 200, n_outputs = 1, \n                              n_static = len(static_var)+1,n_static_hidden=200,sequence_len = x.shape[1], \n                              DEVICE = DEVICE, n_deep_layers=4, use_cuda=USE_CUDA, dropout=0.4).to(DEVICE)\n\nmaize_unscale_model_GLOB = train_model(model, train_loader,val_loader,test_loader,static_var,DEVICE,\n                                  lr=4e-4, weight_decay=0.1,step_size=100,sch_gamma=0.8,num_epochs=200)","metadata":{"execution":{"iopub.status.busy":"2024-09-07T12:07:41.175818Z","iopub.execute_input":"2024-09-07T12:07:41.176333Z","iopub.status.idle":"2024-09-07T12:07:41.183090Z","shell.execute_reply.started":"2024-09-07T12:07:41.176282Z","shell.execute_reply":"2024-09-07T12:07:41.181749Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Wheat","metadata":{}},{"cell_type":"code","source":"## make wheat train_loader\ntrain_loader,val_loader,test_loader = get_training_dataloader('wheat','train',mypath,\n                                                              region=['USA', 'Other', 'South America', 'Europe'],static_var=static_var)","metadata":{"execution":{"iopub.status.busy":"2024-09-07T11:05:10.603251Z","iopub.status.idle":"2024-09-07T11:05:10.603750Z","shell.execute_reply.started":"2024-09-07T11:05:10.603507Z","shell.execute_reply":"2024-09-07T11:05:10.603529Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"(x,z,y) = next(iter(train_loader))","metadata":{"execution":{"iopub.status.busy":"2024-09-07T11:05:10.605543Z","iopub.status.idle":"2024-09-07T11:05:10.606043Z","shell.execute_reply.started":"2024-09-07T11:05:10.605822Z","shell.execute_reply":"2024-09-07T11:05:10.605845Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# # global model training\n# model = LSTMForecaster_static(n_features = x.shape[2], n_hidden = 200, n_outputs = 1, \n#                               n_static = len(static_var)+1,n_static_hidden=200,sequence_len = x.shape[1], \n#                               DEVICE = DEVICE, n_deep_layers=4, use_cuda=USE_CUDA, dropout=0.4).to(DEVICE)\n\n# wheat_unscale_model_GLOB = train_model(model, train_loader,val_loader,test_loader,static_var,DEVICE,\n#                                   lr=4e-4, weight_decay=0.01,step_size=50,sch_gamma=0.8,num_epochs=400)","metadata":{"execution":{"iopub.status.busy":"2024-09-07T11:05:10.608282Z","iopub.status.idle":"2024-09-07T11:05:10.608949Z","shell.execute_reply.started":"2024-09-07T11:05:10.608624Z","shell.execute_reply":"2024-09-07T11:05:10.608660Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# global model training\nmodel = LSTMForecaster_static(n_features = x.shape[2], n_hidden = 200, n_outputs = 1, \n                              n_static = len(static_var)+1,n_static_hidden=200,sequence_len = x.shape[1], \n                              DEVICE = DEVICE, n_deep_layers=4, use_cuda=USE_CUDA, dropout=0.4).to(DEVICE)\n\nwheat_unscale_model_GLOB2 = train_model(model, train_loader,val_loader,test_loader,static_var,DEVICE,\n                                  lr=4e-4, weight_decay=0.1,step_size=100,sch_gamma=0.8,num_epochs=200)","metadata":{"execution":{"iopub.status.busy":"2024-09-07T11:05:10.611395Z","iopub.status.idle":"2024-09-07T11:05:10.611858Z","shell.execute_reply.started":"2024-09-07T11:05:10.611645Z","shell.execute_reply":"2024-09-07T11:05:10.611667Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Evaluation","metadata":{}},{"cell_type":"markdown","source":"## Test climate stats","metadata":{}},{"cell_type":"code","source":"def get_test_dataloader(climate,climate_test,soil,soil_test,region,static_var):\n    # Reshape and calculate mean and sd along the grouped axis\n    climate_test_mean = climate_test.reshape(climate_test.shape[0], 8, 30, 5).mean(axis=2)\n    climate_test_sd = climate_test.reshape(climate_test.shape[0], 8, 30, 5).std(axis=2)\n    GPP_test = np.maximum(climate_test.reshape(climate_test.shape[0], 8, 30, 5)[:,:,:,3] - 10, 0).sum(axis=2,keepdims=True).cumsum(axis=1) # cumulative GDD by phases\n    heat_stress_day_test =  (climate_test.reshape(climate_test.shape[0], 8, 30, 5)[:,:,:,1] > 30).sum(axis=2,keepdims=True).cumsum(axis=1)\n    frost_days_test = (climate_test.reshape(climate_test.shape[0], 8, 30, 5)[:,:,:,2] < 0).sum(axis=2,keepdims=True).cumsum(axis=1)\n    climate_test_stats = np.concatenate([climate_test_sd,climate_test_mean,GPP_test,heat_stress_day_test,frost_days_test],axis=2)\n\n    climate_mean = climate.reshape(climate.shape[0], 8, 30, 5).mean(axis=2)\n    climate_sd = climate.reshape(climate.shape[0], 8, 30, 5).std(axis=2)\n    GPP = np.maximum(climate.reshape(climate.shape[0], 8, 30, 5)[:,:,:,3] - 10, 0).sum(axis=2,keepdims=True).cumsum(axis=1) # cumulative GDD by phases\n    heat_stress_day =  (climate.reshape(climate.shape[0], 8, 30, 5)[:,:,:,1] > 30).sum(axis=2,keepdims=True).cumsum(axis=1)\n    frost_days = (climate.reshape(climate.shape[0], 8, 30, 5)[:,:,:,2] < 0).sum(axis=2,keepdims=True).cumsum(axis=1)\n    climate_stats = np.concatenate([climate_sd,climate_mean,GPP,heat_stress_day,frost_days],axis=2)\n    \n    #annual mean for test\n    annual_test_df = pd.DataFrame(climate_test.mean(axis=1),columns=['tas', 'tasmax', 'tasmin', 'pr', 'rsds'])\n    soil_test = pd.concat([annual_test_df, soil_test], axis=1)\n    \n    #annual mean for train\n    annual_df = pd.DataFrame(climate.mean(axis=1),columns=['tas', 'tasmax', 'tasmin', 'pr', 'rsds'])\n    soil = pd.concat([annual_df, soil], axis=1)\n                                                 \n    #scale input features\n    climate_arr = climate_stats.reshape(-1, climate_stats.shape[2])\n    scaler = StandardScaler()\n    scaler.fit(climate_arr)\n    # scaled_climate = scaler.transform(climate_arr).reshape(climate_stats.shape)\n    \n    #scale static features\n    scaler_static = StandardScaler()\n    static_arr = soil[static_var].values.astype(np.float32)\n    scaler_static.fit(static_arr)\n    # scaler_static.fit(static_arr)\n\n    climate_test_arr = climate_test_stats.reshape(-1, climate_test_stats.shape[2])\n    climate_test_scale = scaler.transform(climate_test_arr).reshape(climate_test_stats.shape)\n\n    static_test_arr = soil_test[static_var].values.astype(np.float32)\n    static_test = scaler_static.transform(static_test_arr)\n    #add co2\n    soil_test['co2'] = soil_test['co2']/1500\n    static_test = np.concatenate([static_test,soil_test['co2'].values.reshape(-1, 1)],axis = 1)\n    \n    X_test_to_loader = torch.from_numpy(climate_test_scale).type(torch.FloatTensor)\n    static_test_loader = torch.from_numpy(static_test).type(torch.FloatTensor)\n    \n    dataset = TensorDataset(X_test_to_loader, static_test_loader)\n\n    # subset by REGION\n    Region_index = soil_test[soil_test['region'].isin(region)].index.values\n    dataset = torch.utils.data.Subset(dataset, Region_index)\n    print(len(dataset))\n    \n    test_loader_stats = torch.utils.data.DataLoader(dataset, 2000, shuffle=False) \n    \n    return soil_test, Region_index, test_loader_stats","metadata":{"execution":{"iopub.status.busy":"2024-09-07T11:05:10.613209Z","iopub.status.idle":"2024-09-07T11:05:10.613698Z","shell.execute_reply.started":"2024-09-07T11:05:10.613443Z","shell.execute_reply":"2024-09-07T11:05:10.613466Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Maize","metadata":{}},{"cell_type":"code","source":"## make maize test_loader\ncrop = 'maize'\n(climate, yield_label, soil) = load(crop,'train', mypath)\n(climate_test, yield_, soil_test_maize) = load(crop,'test', mypath)","metadata":{"execution":{"iopub.status.busy":"2024-09-07T11:05:10.614917Z","iopub.status.idle":"2024-09-07T11:05:10.615329Z","shell.execute_reply.started":"2024-09-07T11:05:10.615131Z","shell.execute_reply":"2024-09-07T11:05:10.615151Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#detrend yield\nsoil = soil.reset_index()\nsoil['yield'] = yield_label\ngrouped_mean = soil.groupby(['lon', 'lat'])['yield'].mean().reset_index()\ngrouped_mean = grouped_mean.rename(columns={'yield': 'yield_mean'})\nsoil = pd.merge(soil, grouped_mean, on=['lon', 'lat'])\n\nsoil_test_maize = soil_test_maize.reset_index()\nsoil_test_maize = pd.merge(soil_test_maize, grouped_mean, on=['lon', 'lat'])\n\nsoil['de_yield'] = soil['yield']-soil['yield_mean']\nsoil['region'] = soil.apply(get_region, axis=1)\nsoil_test_maize['region'] = soil_test_maize.apply(get_region, axis=1)","metadata":{"execution":{"iopub.status.busy":"2024-09-07T11:05:10.616998Z","iopub.status.idle":"2024-09-07T11:05:10.617409Z","shell.execute_reply.started":"2024-09-07T11:05:10.617208Z","shell.execute_reply":"2024-09-07T11:05:10.617227Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_soil, Region_index_US, test_loader_US = get_test_dataloader(climate,climate_test,soil,soil_test_maize,region=['USA', 'Other', 'South America', 'Europe'],static_var=static_var)","metadata":{"execution":{"iopub.status.busy":"2024-09-07T11:05:10.618573Z","iopub.status.idle":"2024-09-07T11:05:10.619030Z","shell.execute_reply.started":"2024-09-07T11:05:10.618802Z","shell.execute_reply":"2024-09-07T11:05:10.618826Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"maize_pred_stats_test = []\nfor i, (x,z) in enumerate(test_loader_US):\n    with torch.no_grad():\n        # y_pred = maize_unscale_model_US(x.to(DEVICE),z.to(DEVICE))\n        y_pred = maize_unscale_model_GLOB(x.to(DEVICE),z.to(DEVICE))\n        maize_pred_stats_test.append(y_pred.detach().cpu().numpy())","metadata":{"execution":{"iopub.status.busy":"2024-09-07T11:05:10.620385Z","iopub.status.idle":"2024-09-07T11:05:10.620827Z","shell.execute_reply.started":"2024-09-07T11:05:10.620614Z","shell.execute_reply":"2024-09-07T11:05:10.620636Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Wheat","metadata":{}},{"cell_type":"code","source":"## make wheat test_loader\ncrop = 'wheat'\n(climate, yield_label, soil) = load(crop,'train', mypath)\n(climate_test, yield_, soil_test_wheat) = load(crop,'test', mypath)","metadata":{"execution":{"iopub.status.busy":"2024-09-07T11:05:10.622456Z","iopub.status.idle":"2024-09-07T11:05:10.622917Z","shell.execute_reply.started":"2024-09-07T11:05:10.622706Z","shell.execute_reply":"2024-09-07T11:05:10.622728Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#detrend yield\nsoil = soil.reset_index()\nsoil['yield'] = yield_label\ngrouped_mean = soil.groupby(['lon', 'lat'])['yield'].mean().reset_index()\ngrouped_mean = grouped_mean.rename(columns={'yield': 'yield_mean'})\nsoil = pd.merge(soil, grouped_mean, on=['lon', 'lat'])\n\nsoil_test_wheat = soil_test_wheat.reset_index()\nsoil_test_wheat = pd.merge(soil_test_wheat, grouped_mean, on=['lon', 'lat'])\n\nsoil['de_yield'] = soil['yield']-soil['yield_mean']\nsoil['region'] = soil.apply(get_region, axis=1)\nsoil_test_wheat['region'] = soil_test_wheat.apply(get_region, axis=1)","metadata":{"execution":{"iopub.status.busy":"2024-09-07T11:05:10.624447Z","iopub.status.idle":"2024-09-07T11:05:10.624899Z","shell.execute_reply.started":"2024-09-07T11:05:10.624688Z","shell.execute_reply":"2024-09-07T11:05:10.624711Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_soil, Region_index_US, test_loader_US = get_test_dataloader(climate,climate_test,soil,soil_test_wheat,region=['USA', 'Other', 'South America', 'Europe'],static_var=static_var)","metadata":{"execution":{"iopub.status.busy":"2024-09-07T11:05:10.626930Z","iopub.status.idle":"2024-09-07T11:05:10.627518Z","shell.execute_reply.started":"2024-09-07T11:05:10.627205Z","shell.execute_reply":"2024-09-07T11:05:10.627236Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"wheat_pred_stats_test = []\nfor i, (x,z) in enumerate(test_loader_US):\n    with torch.no_grad():\n        # y_pred = maize_unscale_model_US(x.to(DEVICE),z.to(DEVICE))\n        y_pred = wheat_unscale_model_GLOB2(x.to(DEVICE),z.to(DEVICE))\n        wheat_pred_stats_test.append(y_pred.detach().cpu().numpy())","metadata":{"execution":{"iopub.status.busy":"2024-09-07T11:05:10.629973Z","iopub.status.idle":"2024-09-07T11:05:10.630619Z","shell.execute_reply.started":"2024-09-07T11:05:10.630281Z","shell.execute_reply":"2024-09-07T11:05:10.630321Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Create submission","metadata":{}},{"cell_type":"code","source":"sub_v5 = pd.read_csv('/kaggle/input/the-future-crop-challenge/sample_submission.csv')","metadata":{"execution":{"iopub.status.busy":"2024-09-07T11:05:10.631978Z","iopub.status.idle":"2024-09-07T11:05:10.632411Z","shell.execute_reply.started":"2024-09-07T11:05:10.632204Z","shell.execute_reply":"2024-09-07T11:05:10.632225Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"soil_test_maize['preds_stats'] = np.concatenate(maize_pred_stats_test)\nsoil_test_maize['preds_yield']  = soil_test_maize['preds_stats'] +soil_test_maize['yield_mean']\nsoil_test_maize['preds_yield'] =  soil_test_maize['preds_yield'].where(soil_test_maize['preds_yield'] >=0, soil_test_maize['yield_mean'] )\nsoil_test_maize.describe()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"soil_test_wheat['preds_stats'] = np.concatenate(wheat_pred_stats_test)\nsoil_test_wheat['preds_yield']  = soil_test_wheat['preds_stats'] +soil_test_wheat['yield_mean']\nsoil_test_wheat['preds_yield'] =  soil_test_wheat['preds_yield'].where(soil_test_wheat['preds_yield'] >=0, soil_test_wheat['yield_mean'] )\nsoil_test_wheat.describe()","metadata":{"execution":{"iopub.status.busy":"2024-09-07T11:05:10.634146Z","iopub.status.idle":"2024-09-07T11:05:10.634633Z","shell.execute_reply.started":"2024-09-07T11:05:10.634375Z","shell.execute_reply":"2024-09-07T11:05:10.634397Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sub_v5.loc[sub_v5['ID'].isin(soil_test_maize['ID']),'yield']  = soil_test_maize['preds_yield'].values\nsub_v5.loc[sub_v5['ID'].isin(soil_test_wheat['ID']),'yield']  = soil_test_wheat['preds_yield'].values","metadata":{"execution":{"iopub.status.busy":"2024-09-07T11:05:10.636243Z","iopub.status.idle":"2024-09-07T11:05:10.636721Z","shell.execute_reply.started":"2024-09-07T11:05:10.636468Z","shell.execute_reply":"2024-09-07T11:05:10.636490Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sub_v5.to_csv('submission.csv',index=0)","metadata":{"execution":{"iopub.status.busy":"2024-09-07T11:05:10.639006Z","iopub.status.idle":"2024-09-07T11:05:10.639445Z","shell.execute_reply.started":"2024-09-07T11:05:10.639228Z","shell.execute_reply":"2024-09-07T11:05:10.639250Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}