{"metadata":{"kernelspec":{"name":"ir","display_name":"R","language":"R"},"language_info":{"name":"R","codemirror_mode":"r","pygments_lexer":"r","mimetype":"text/x-r-source","file_extension":".r","version":"4.0.5"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":81000,"databundleVersionId":8812083,"sourceType":"competition"}],"dockerImageVersionId":30618,"isInternetEnabled":true,"language":"r","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# This R environment comes with many helpful analytics packages installed\n# It is defined by the kaggle/rstats Docker image: https://github.com/kaggle/docker-rstats\n# For example, here's a helpful package to load\n\n#library(tidyverse) # metapackage of all tidyverse packages\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\n#list.files(path = \"../input\")\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":"1479e5c0-ecdc-4bdb-816f-ea2d84d32ac0","_cell_guid":"fcccae82-43af-45c7-925d-ee0180c735f5","collapsed":false,"_execution_state":"idle","_kg_hide-input":true,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2024-07-06T12:18:47.891732Z","iopub.execute_input":"2024-07-06T12:18:47.924396Z","iopub.status.idle":"2024-07-06T12:18:47.936435Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Time-series aggregation\n\nThis example uses an approach to aggregate the time-series that takes into account what is mentioned in the description with regard to the variability of season length. We are going to aggregate the time-series using periods thresholded by what is called accumulated thermal time. Basically, the progress in a season for plant growth is connected to how warm are the growth conditions. If the temperatures are within an optimal range, this means there is no delay to growth. We are going to use a very simple approach to calculate these periods, which only refer to the temperature above which plants grow. There are more sophisticated approaches.\n\nThe functions below determine: 1. for each gridcell, in what day the accumulated thermal time threshold is reached, and 2. the summarized time-series for each gridcell.\n\nP.s.: This is the first script I posted for the competition and I'll leave it as is, but more recent scripts were optimized and allow for fractioning the season in more periods, which can't be done here because the notebook crashes.","metadata":{"_uuid":"793d56c6-5b39-48c6-a3a2-4381e117ae9e","_cell_guid":"206e0dd0-a850-4efa-904e-622d82f0bfb7","trusted":true}},{"cell_type":"code","source":"# Initialization and libraries\nlibraries <- c(\"here\", \"tidyverse\", \"data.table\", \"zoo\", \"arrow\",\n               # plots\n               \"ggplot2\",\n               # ML\n               \"xgboost\")\n\nlapply(libraries, require, character.only = TRUE)","metadata":{"_uuid":"68eecc28-460a-4d33-bb31-ceae8dd779fd","_cell_guid":"2a03b26e-b63a-48af-95db-cede2c1c4693","collapsed":false,"_kg_hide-output":true,"scrolled":true,"_kg_hide-input":true,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2024-07-06T12:18:47.939809Z","iopub.execute_input":"2024-07-06T12:18:47.941282Z","iopub.status.idle":"2024-07-06T12:18:49.778705Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def_tt <- function(x, tbase, thresholds){\n  \n    # Returns das for each predefined threshold of accumulated thermal time in the cycle\n    # x: ID and 30 to 209\n    # x_: ID, das (days after simulation started) and correspondent accumulated tt\n    \n    # tbase <- 8\n\n    # thresholds <- c(1000, 2000)\n    thresholds_name <- paste0(\"tt\", thresholds)\n    \n        \n    # Calculates thermal time for each day and accumulated thermal time for each ID\n    x_tt <- data.table::melt(x, id.vars = \"ID\", variable.name = \"das\",\n                           variable.factor = FALSE)\n    x_tt[, das := as.numeric(das)]\n    \n    x_tt[, value := ifelse(value - tbase < 0, 0, value - tbase)]\n    \n    x_tt <- x_tt[order(ID, das)][, tt := cumsum(value), \n                               keyby = c(\"ID\")][, value := NULL]\n    x_tt[, (thresholds_name) := as.list(thresholds\n                                        # 0.1*max(tt), 0.25*max(tt), 1/3)*max(tt), 0.5*max(tt), 800\n                                        ), \n       keyby = ID]\n    \n    # Identifies dates closest to the thresholds\n    x_tt <- data.table::melt(x_tt, \n                             id.vars = c(\"ID\", \"das\", \"tt\"), \n                             measure.vars = thresholds_name,\n                             variable.name = \"accum_tt\")\n    x_tt[, dif:= abs(tt-value)][, value := NULL]\n    x_tt[, flag := (dif == min(dif)), keyby = c(\"ID\", \"accum_tt\")]\n    \n    x_tt <- x_tt[flag == TRUE]\n    x_tt <- x_tt[x_tt[, .I[flag][1], by = .(ID, accum_tt)][, V1]]\n    x_tt <- x_tt[, .(ID, das, tt, accum_tt)]\n    x_tt <- unique(x_tt, by = c(\"das\", \"ID\"))\n    setkey(x_tt, NULL)\n    \n    return(x_tt)\n  \n}\n\nfeat_eng <- function(x, tt_das, var_){\n    \n    # Returns aggregated time-series for each threshold of tt of the variable\n    # x: ID and 30 to 239\n    # tt_das: ID, das and correspondent tt\n    # var_: variable to which dataset refers\n    \n    x_long <- data.table::melt(x, id.vars = \"ID\", variable.name = \"das\",\n                               variable.factor = FALSE)\n    x_long <- x_long[, das := as.numeric(das)][order(ID, das)]\n    \n    if(var_ == \"temp\"){\n        funcs <- list(\"avg\" = mean, \"min\" = min, \"max\" = max)\n    } else if (var_ == \"tmax\"){\n        funcs <- list(\"avg\" = mean, \"max\" = max)\n    }else if (var_ == \"tmin\"){\n        funcs <- list(\"avg\" = mean, \"min\" = min)\n    }else if (var_ == \"srad\"){\n        funcs <- list(\"avg\" = mean, \"sum\" = sum)\n    } else if (var_ == \"prec\"){\n        funcs <- list(\"sum\" = sum)\n    }\n    \n    if (!is.null(tt_das)){\n        \n        # Create groups by thresholded period\n        \n        x_mod <- merge(x_long, tt_das, all = TRUE)[, tt := NULL]\n        rm(x_long)\n        x_mod[, accum_tt := as.numeric(accum_tt)]\n        # We will ascribe the groups from last day to first, so we use 0 as a\n        # placeholder for group code to the last day if it exceeded the larger threshold,\n        # carry the value until the previous group is found code and\n        # end x_mod with ID, das, tt value and accum_tt group\n        x_mod[, accum_tt := ifelse(das == max(das) & is.na(accum_tt), 0, accum_tt),\n              keyby = \"ID\"]\n        x_mod[, accum_tt := zoo::na.locf(accum_tt, fromLast=TRUE)]\n        x_mod <- x_mod[accum_tt != 0]\n        \n        x_all <- setkey(data.table(unique(x_mod[, \"ID\"])), NULL)\n        \n        if(var_ == \"temp\"){\n            # We then add, just for analysis, which day of the cycle corresponds\n            # to the date of each threshold\n            tt_das <- data.table::dcast(tt_das, ID ~ accum_tt,\n                                        value.var = \"das\")[, ID := NULL]\n            x_all <- cbind(x_all, tt_das)\n            \n        }\n        \n        rm(tt_das)\n        \n        # Generates aggregation\n        for (func_name in names(funcs)) {\n            # Here we generate the aggregated series for each period and\n            # join them to the aggregation for the whole season\n            x_periods <- x_mod[, .(funcs[[func_name]](value)), by = c(\"ID\", \"accum_tt\")]\n            x_periods[, new_name := paste(var_, func_name,\n                                          accum_tt, sep = \"_\")][, accum_tt := NULL]\n            x_periods <- setkey(data.table::dcast(x_periods, \n                                                  ID ~ new_name, value.var = \"V1\"), NULL)\n            \n            x_full <- setkey(x_mod[, .(funcs[[func_name]](value)), by = c(\"ID\")], NULL)\n            x_full[, new_name := paste(var_, func_name, \"season\", sep = \"_\")]\n            x_full <- setkey(data.table::dcast(x_full, ID ~ new_name, value.var = \"V1\"), NULL)\n            \n            temp <- cbind(x_periods, x_full[, ID := NULL])\n            x_all <- setkey(merge(x_all, temp, all = TRUE), NULL)\n            \n        }\n        \n    } else {\n        \n        # Aggregation not performed using thermal time info (e.g., precipitation before growth)\n        x_mod <- x_long\n        rm(x_long)\n        x_all <- data.table(unique(x_mod[, \"ID\"]))\n        for (func_name in names(funcs)) {\n            x_full <- x_mod[, .(funcs[[func_name]](value)), by = c(\"ID\")]\n            x_full[, new_name := paste0(func_name, \"_season\")]\n            x_full <- data.table::dcast(x_full, ID ~ new_name, value.var = \"V1\")\n            \n            x_all <- setkey(merge(x_all, x_full, all = TRUE), NULL)\n            \n        }\n        \n    }\n    \n    return(x_all)\n    \n}","metadata":{"_uuid":"0ff0fc9c-15df-4831-aeec-e149302c53fa","_cell_guid":"18361303-1866-4feb-992a-d2c07c49cfa8","collapsed":false,"_kg_hide-input":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2024-06-30T14:23:34.120405Z","iopub.execute_input":"2024-06-30T14:23:34.122439Z","iopub.status.idle":"2024-06-30T14:23:34.145304Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Additional helper functions will be used to load and join the other datasets and perform initial preprocessing of soil data and missing values in generated in the feature engineering step.","metadata":{"_uuid":"a588411c-9b4e-41ce-8004-91edad494695","_cell_guid":"61d536b5-fad0-4a50-85e7-3f3c07e1d1dd","trusted":true}},{"cell_type":"code","source":"generate_set <- function(crop, set, max_rows, files_,\n                         tbase=8, thresholds=c(\"1000\", \"2000\")){\n    \n    cols_growth <- c(\"ID\", paste0(c(as.character(seq(30,239)))))\n    cols_pre <- c(\"ID\", paste0(c(as.character(seq(1,30)-1))))\n        \n    # Soil and IDs ------------------------------------------------------------\n    filename <- paste(\"soil_co2\", crop, set, sep = \"_\")\n    soil <- read_data(filename, files_, max_rows, cols_sel=NULL)\n    \n    soil <- soil[, texture_class := paste0(\"soil_\", letters[texture_class])] \n    soil <- soil[, crop := ifelse(crop == \"maize\", 0, 1)]\n    \n    soil <- data.table::dcast(soil, ID + crop + co2 + lat + lon ~ texture_class, \n                              fun.aggregate = length) %>%\n        complete_soil()\n\n    # Tmean -------------------------------------------------------------------\n    # Must be first to determine thresholds\n    filename <- paste(\"tas\", crop, set, sep = \"_\")\n    x <- read_data(filename, files_, max_rows, cols_sel=cols_growth)\n    \n    # Thermal time\n    tt_das <- def_tt(x, tbase, thresholds)\n    temp <- feat_eng(x, tt_das, \"temp\")\n\n    # TMax --------------------------------------------------------------------\n    filename <- paste(\"tasmax\", crop, set, sep = \"_\")\n    x <- read_data(filename, files_, max_rows, cols_sel=cols_growth)\n    \n    tmax <- feat_eng(x, tt_das, \"tmax\")\n\n    # TMin --------------------------------------------------------------------\n    filename <- paste(\"tasmin\", crop, set, sep = \"_\")\n    x <- read_data(filename, files_, max_rows, cols_sel=cols_growth)\n    \n    tmin <- feat_eng(x, tt_das, \"tmin\")\n\n    # SRad --------------------------------------------------------------------\n    filename <- paste(\"rsds\", crop, set, sep = \"_\")\n    x <- read_data(filename, files_, max_rows, cols_sel=cols_growth)\n    \n    srad <- feat_eng(x, tt_das, \"srad\")\n\n    # Prec --------------------------------------------------------------------\n    # Previous, to infer soil moisture at the start of the season\n    filename <- paste(\"pr\", crop, set, sep = \"_\")\n    \n    x <- read_data(filename, files_, max_rows, cols_sel=cols_pre)\n    prec_pre <- feat_eng(x, tt_das=NULL, \"prec\")\n    setnames(prec_pre, old = \"sum_season\", new = \"sum_pre\")\n    prec_pre[, max_season := NULL]\n    \n    # During growth\n    x <- read_data(filename, files_, max_rows, cols_sel=cols_growth)\n    prec <- feat_eng(x, tt_das, \"prec\")\n    \n    # Join data ---------------------------------------------------------------\n\n    if(!all((soil$ID == temp$ID),\n            (soil$ID == tmax$ID),\n            (soil$ID == tmin$ID),\n            (soil$ID == srad$ID),\n            (soil$ID == prec$ID),\n            (soil$ID == prec_pre$ID))){\n        return(\"dataset not aligned\")\n    }\n    \n    dataset <- cbind(soil, temp[, !c(\"ID\")]) %>%\n        cbind(tmax[, !c(\"ID\")]) %>%\n        cbind(tmin[, !c(\"ID\")]) %>%\n        cbind(srad[, !c(\"ID\")]) %>%\n        cbind(prec_pre[, !c(\"ID\")]) %>%\n        cbind(prec[, !c(\"ID\")])\n    \n    if(set==\"train\"){\n        \n        filename <- paste(set, \"solutions\", crop, sep = \"_\")\n        x <- read_data(filename, files_, max_rows, cols_sel=NULL)\n        \n        dataset <- cbind(dataset, x[, !c(\"ID\")]) \n        \n    }\n\n    return(dataset)\n    \n}","metadata":{"_uuid":"16a53748-8b31-41b3-b14e-5b3b9a66d582","_cell_guid":"3a7f5b3d-ae55-4446-8aab-237c185d5bf8","collapsed":false,"_kg_hide-input":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2024-06-30T14:23:34.149830Z","iopub.execute_input":"2024-06-30T14:23:34.151920Z","iopub.status.idle":"2024-06-30T14:23:34.174208Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"read_data <- function(filename, files_, max_rows, cols_sel=NULL){\n    \n    if(is.null(cols_sel)) {\n        \n        if(is.null(max_rows)){\n            \n            x <- read_parquet(files_[grepl(filename, files_)]) %>%\n                as.data.table()\n        } else {\n            \n            x <- open_dataset(files_[grepl(filename, files_)]) %>%\n                head(n=max_rows) %>% \n                compute() %>%\n                as.data.table()\n        }\n        \n    } else {\n     \n        if(is.null(max_rows)){\n            x <- read_parquet(files_[grepl(filename, files_)], \n                              col_select=cols_sel) %>%\n                as.data.table()\n        } else {\n            x <- open_dataset(files_[grepl(filename, files_)]) %>%\n                head(n=max_rows) %>% \n                compute() %>%\n                as.data.table()\n            \n            x <- x[, ..cols_sel]\n        }\n           \n    }\n    \n    return(x)\n}","metadata":{"_uuid":"6a25a266-94bb-47b2-8c82-364bba282d7f","_cell_guid":"ce1fff72-538d-4d48-8ce9-7a101ddbcac9","collapsed":false,"_kg_hide-input":true,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2024-06-30T14:23:34.178110Z","iopub.execute_input":"2024-06-30T14:23:34.181112Z","iopub.status.idle":"2024-06-30T14:23:34.199388Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"complete_soil <- function(x){\n    \n    all_soils <- paste0(\"soil_\", letters[(1:13)])\n    soils_missing <- all_soils[!all_soils %in% colnames(x)]\n    \n    if(length(soils_missing) > 0){\n        x[, (soils_missing) := NA]\n    }\n    \n    return(x)\n    \n}\n\nfill_na <- function(x){\n    \n    # Replace NAs in soil texture class variable\n    names_cols <- grep(\"soil\", colnames(x), value = TRUE)\n    x_ <- x[, (names_cols) := lapply(.SD, function(x) ifelse(is.na(x), 0, x)), .SDcols = names_cols]\n    \n    # Lazily replace NAs with constant value for the cases in which the season ended before the threshold\n    names_cols <- grep(\"temp|tmax|tmin|srad|prec|tt\", colnames(x), value = TRUE)\n    x_ <- x_[, (names_cols) := lapply(.SD, function(x) ifelse(is.na(x), -999, x)), .SDcols = names_cols]\n    \n    return(x_)\n                                      \n}","metadata":{"_uuid":"367a4b4c-0ba3-4523-8ff8-ec366d16fec1","_cell_guid":"c9a3a92a-5be9-45ec-bc0b-bdfcffc0bb46","collapsed":false,"_kg_hide-input":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2024-06-30T14:23:34.203693Z","iopub.execute_input":"2024-06-30T14:23:34.205722Z","iopub.status.idle":"2024-06-30T14:23:34.250852Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Obtaining the training set","metadata":{"_uuid":"739bc4d9-7e19-4c46-ac0b-b87747a93ec3","_cell_guid":"f02e678b-f971-4e77-aaca-e2f4d4d7f230","trusted":true}},{"cell_type":"code","source":"# Constants ---------------------------------------------------------------\nmax_rows_ <- NULL # Allows for testing a subset of the dataset; NULL reads the whole file\npath_inputs <- \"../input/the-future-crop-challenge/\"\nfiles_ <- list.files(path_inputs, full.names = TRUE, pattern = \".parquet\")\ntt_thresholds <- c(500, 1500)","metadata":{"_uuid":"95a7f2a1-7178-47b5-800e-97c4564fcdd9","_cell_guid":"7e540584-e913-4d61-ba79-a4af792639cd","collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2024-06-30T14:23:34.253957Z","iopub.execute_input":"2024-06-30T14:23:34.255573Z","iopub.status.idle":"2024-06-30T14:23:34.292012Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Train -------------------------------------------------------------------\ntr_maize <- generate_set(crop=\"maize\", set=\"train\", max_rows=max_rows_, files_,\n                         tbase=8, thresholds=tt_thresholds)\n\ntr_wheat <- generate_set(crop=\"wheat\", set=\"train\", max_rows=max_rows_, files_,\n                         tbase=0, thresholds=tt_thresholds)\n\ntrain <- rbind(tr_maize, tr_wheat, fill=TRUE) %>%\n    fill_na(.)\n\ntrain <- train[, !c(\"ID\")]\n\nrm(tr_maize, tr_wheat)\ngc()","metadata":{"_uuid":"4b84e272-ee6d-468d-aac2-0477143925ca","_cell_guid":"d99f82cb-97a8-411b-8db3-cdddad88c976","collapsed":false,"_kg_hide-output":true,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2024-06-30T14:23:34.296524Z","iopub.execute_input":"2024-06-30T14:23:34.298364Z","iopub.status.idle":"2024-06-30T14:39:22.748303Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"head(train)","metadata":{"_uuid":"7f7ba6ab-4fa9-4a10-9d29-8c6fd7c6b764","_cell_guid":"49e6ad3b-5e01-4534-bdb3-75455aa98f5c","collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2024-06-30T14:39:22.751602Z","iopub.execute_input":"2024-06-30T14:39:22.754151Z","iopub.status.idle":"2024-06-30T14:39:22.824232Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"colnames(train)","metadata":{"_uuid":"3b71945b-2d92-48ba-b979-3f0ec1a83edb","_cell_guid":"1126859d-c56c-46f5-a64e-2ad26745aee3","collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2024-06-30T14:39:22.827273Z","iopub.execute_input":"2024-06-30T14:39:22.828868Z","iopub.status.idle":"2024-06-30T14:39:22.847606Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dim(train)","metadata":{"_uuid":"3e8e341d-40ec-4932-a7de-d7b447666993","_cell_guid":"47b8035d-5df7-4673-93a8-6c16c434eeb4","collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2024-06-30T14:39:22.850558Z","iopub.execute_input":"2024-06-30T14:39:22.852149Z","iopub.status.idle":"2024-06-30T14:39:22.870648Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Season length analysis\nWe will inspect the features generated to assess to which day of the growing season correspond the thresholds assessed. We will divide them by crop because not only both crops may grow in different locations, but also because, as you may remember, we used different base temperatures. We will also remove the value that was included as a placeholder for missing values.\n\nIf 1500 is a good value to represent the end of the season, note how, for maize, it could be reached within 70 days (remember we start at day 30), but it could also take more than 200 days. There are, however, many details that are being overlooked in this analysis (such as which are good values to represent relevant periods in growth and season length, ways of calculating thermal time, better ways to deal with the missing values, etc) and by looking into them, it is possible that model performance could be improved.","metadata":{"_uuid":"0b98d88e-1832-4b32-b11b-058e8544287f","_cell_guid":"5ce05beb-df21-48a3-8c42-44983886eaee","trusted":true}},{"cell_type":"code","source":"dt_plot <- train[, .(crop, tt500, tt1500)][tt1500 != -999]\n\nggplot() +\n    facet_wrap(\"crop\", labeller = labeller(crop = c(\"0\" = \"maize\", \"1\" = \"wheat\"))) +\n    geom_histogram(data = dt_plot, aes(tt500), fill = \"lightblue\") +\n    geom_histogram(data = dt_plot, aes(tt1500), fill = \"lightgreen\") +\n    labs(x = \"Day of growth in which the threshold is reached\")","metadata":{"_uuid":"b10e51ee-09a6-4550-ade5-9bb028722d94","_cell_guid":"104523ec-596e-4361-ac19-9db2c6ed88f4","collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2024-06-30T14:40:07.820691Z","iopub.execute_input":"2024-06-30T14:40:07.822463Z","iopub.status.idle":"2024-06-30T14:40:10.601927Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Train model","metadata":{"_uuid":"03c7659a-d9ea-4a7d-844e-b2c006b8cf57","_cell_guid":"68da2797-bfad-4921-92dd-6b4efadef326","trusted":true}},{"cell_type":"code","source":"rem_cols <- grep(\"tt\", colnames(train), value = TRUE)\ntrain <- train[, !..rem_cols]\n\nmodel_str <- \"xgb\"\nsuffix <- \"002_stages500-1000\"\n\nif(model_str == \"avg\"){\nmodel_ <- round(mean(train$yield), 4)\n} else if (model_str == \"rpart\"){\nmodel_ <- rpart(yield~., data=train)\n} else if (model_str == \"rf\"){\nmodel_ <- randomForest(yield~., data=train, ntree=50)\n} else if (model_str == \"xgb\"){\nmodel_ <- xgboost(data=as.matrix(copy(train)[, yield:=NULL]),\n                  label=copy(train)[[\"yield\"]],\n                  nrounds = 20)\n}\n\nsolution <- train[, .(yield)]\n\n# Errors\nif(model_str != \"avg\"){\n  if(model_str == \"xgb\"){\n    yield_ <- predict(model_, as.matrix(copy(train)[, yield:=NULL])) \n  } else {\n    yield_ <- predict(model_, train) \n  }\n} else {\n  yield_ <- model_\n}\n\n# rm(train)\n\nsolution[, pred := round(as.numeric(yield_), 6)]\nsolution[, err := (pred - yield)]\nerror_est <- sqrt(mean(solution$err^2))\n\nprint(error_est)\nfilename <- paste0(\"errorest_\", model_str, suffix, \".csv\")\nwrite.csv(error_est, file=here(filename), row.names=FALSE, quote=FALSE)","metadata":{"_uuid":"e937c535-94ff-487a-8bf3-0c5fab8b602c","_cell_guid":"d8cf888a-2fa2-45ca-b000-e0fb30337965","collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2024-06-30T14:40:10.605090Z","iopub.execute_input":"2024-06-30T14:40:10.606824Z","iopub.status.idle":"2024-06-30T14:41:11.457540Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"xgb.importance(model=model_)","metadata":{"_uuid":"756d426c-4092-4bd8-96bb-1603f30a6784","_cell_guid":"911b4281-6c8e-4d49-aebd-fe85076675f9","collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2024-06-30T14:41:11.460694Z","iopub.execute_input":"2024-06-30T14:41:11.463113Z","iopub.status.idle":"2024-06-30T14:41:11.643547Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Feature generation for the test set","metadata":{"_uuid":"95abfd3f-c59d-4f66-bd3d-bde773627283","_cell_guid":"b83d2f75-873b-413b-a8fa-8de4129142a7","trusted":true}},{"cell_type":"code","source":"tst_maize <- generate_set(crop=\"maize\", set=\"test\", max_rows=max_rows_, files_,\n                         tbase=8, thresholds=tt_thresholds)\n\ntst_wheat <- generate_set(crop=\"wheat\", set=\"test\", max_rows=max_rows_, files_,\n                         tbase=0, thresholds=tt_thresholds)\n\ntest <- rbind(tst_maize, tst_wheat, fill=TRUE) %>%\n    fill_na(.)\n\nrem_cols <- c(grep(\"tt\", colnames(test), value = TRUE), \"ID\")\nids <- test[, c(\"ID\")]\ntest <- test[, !..rem_cols]\norder_cols <- colnames(train)[colnames(train) != \"yield\"]\ntest <- test[, ..order_cols]\n\nrm(tst_maize, tst_wheat)\ngc()","metadata":{"_uuid":"8e5d6194-53a6-4286-8dc7-49f902928da9","_cell_guid":"21787a99-43e7-4d76-8494-610c8d0b2fd5","collapsed":false,"_kg_hide-output":true,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2024-06-30T14:41:11.647246Z","iopub.execute_input":"2024-06-30T14:41:11.648916Z","iopub.status.idle":"2024-06-30T15:11:29.033993Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Prediction","metadata":{"_uuid":"5f554335-132b-4858-9f2d-444c74c21c2f","_cell_guid":"c4299553-614a-4922-bf52-88201e930efc","trusted":true}},{"cell_type":"code","source":"if(model_str != \"avg\"){\n  if(model_str == \"xgb\"){\n    yield_ <- predict(model_, as.matrix(copy(test)))\n  } else {\n    yield_ <- predict(model_, test) \n  }\n} else {\n  yield_ <- model_\n}\n\nsolution2 <- ids\nsolution2[, yield := round(as.numeric(yield_), 6)]\nfilename <- paste0(\"submission_\", model_str, suffix, \".csv\")\nwrite.csv(solution2, file=here(filename), row.names=FALSE, quote=FALSE)\n\nfilename <- \"submission.csv\"\nwrite.csv(solution2, file=here(filename), row.names=FALSE, quote=FALSE)\n\ngc()","metadata":{"_uuid":"06c06dfc-3843-4791-be69-b12f0fc6a24a","_cell_guid":"cd49769a-7ef2-45a0-b384-9c9eeef7ecc1","collapsed":false,"_kg_hide-output":true,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2024-06-30T15:11:29.038411Z","iopub.execute_input":"2024-06-30T15:11:29.040111Z","iopub.status.idle":"2024-06-30T15:11:36.805235Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"head(solution2)","metadata":{"_uuid":"f43f528e-b711-461b-917d-4445c5ad5dc9","_cell_guid":"05516036-bb45-4264-abef-4272bec93148","collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2024-06-30T15:11:36.809844Z","iopub.execute_input":"2024-06-30T15:11:36.811753Z","iopub.status.idle":"2024-06-30T15:11:36.837054Z"},"trusted":true},"execution_count":null,"outputs":[]}]}