Add module 5: Penalty in case of exceeding legal limit

DESCRIPTION

@@ -33,6 +33,7 @@ Imports:
     R6,
     torch
 Suggests: 
+    devtools,
     testthat (>= 3.0.0)
 Config/testthat/edition: 3
 Depends: 

NEWS.md

@@ -2,5 +2,5 @@
 
 ## Added
 * Adds `apus` object with the functions `addField`, `trainModel` and `optimizeFertilizerChoice`
-* Adds for the cost function module 1 and 2
+* Adds for the cost function module 1,5 and 6
 * Adds default table for `cultivations`, `fertilizers`, `parameters` and `fines`

R/apus.R

@@ -168,6 +168,7 @@ Apus <- R6::R6Class(
       dataset.train <- createApusDataset(farms = NULL,
                                          cultivation = self$cultivation,
                                          fertilizers = self$fertilizers,
+                                         fines = self$fines,
                                          fields_max = self$fields_max,
                                          device = device)
 
@@ -177,6 +178,7 @@ Apus <- R6::R6Class(
       dataset.valid <- createApusDataset(farms = farms.valid,
                                          cultivation = self$cultivation,
                                          fertilizers = self$fertilizers,
+                                         fines = self$fines,
                                          fields_max = self$fields_max,
                                          device = device)
 
@@ -217,7 +219,7 @@ Apus <- R6::R6Class(
         fields[is.na(fields)] <- 0
       }
 
-      dataset <- createApusDataset(farms = fields, cultivations = self$cultivations, fertilizers = self$fertilizers, fields_max = self$fields_max, device = self$device)
+      dataset <- createApusDataset(farms = fields, cultivations = self$cultivations, fertilizers = self$fertilizers, fines = self$fines, fields_max = self$fields_max, device = self$device)
       dl <- torch::dataloader(dataset, batch_size = 1)
 
 

R/dataset.R

@@ -6,6 +6,7 @@
 #' @param farms (data.table)
 #' @param cultivations (data.table)
 #' @param fertilizers (data.table)
+#' @param fines (data.table)
 #' @param fields_max (integer)
 #' @param device (character)
 #'
@@ -15,7 +16,7 @@
 #' @import torch
 #'
 #'@export
-createApusDataset <- function(farms = NULL, cultivations, fertilizers, fields_max, device) {
+createApusDataset <- function(farms = NULL, cultivations, fertilizers, fines, fields_max, device) {
 
   transformfieldsToTensor = createSyntheticfields = code = fields_count = self = NULL
   size = value_max = value_min = p_price = p_stored = b_id_farm = b_id_field = NULL
@@ -27,7 +28,7 @@ createApusDataset <- function(farms = NULL, cultivations, fertilizers, fields_ma
   apus_dataset <- torch::dataset(
     name = "apus_dataset",
 
-    initialize = function(farms = NULL, cultivations, fertilizers, fields_max, device) {
+    initialize = function(farms = NULL, cultivations, fertilizers, fines, fields_max, device) {
 
       # Check arguments -----------------------------------------------------
       # TODO
@@ -45,11 +46,14 @@ createApusDataset <- function(farms = NULL, cultivations, fertilizers, fields_ma
         self$farms_count <- 100
       }
 
+      fines <- dcast(fines, . ~ norm, value.var = 'fine')[, 2:4]
+      self$fines <- torch::torch_tensor(as.matrix(fines), device = device)
+
       # Set temporary
       fertilizers[, p_stored := 0]
       fertilizers[, p_price := 1]
 
-      fertilizers <- fertilizers[, c('p_stored', 'p_price', 'p_n_rt', 'p_n_wc', 'p_p_rt', 'p_k_rt')]
+      fertilizers <- fertilizers[, c('p_stored', 'p_price', 'p_n_rt', 'p_n_wc', 'p_p_rt', 'p_k_rt', 'p_type_manure', 'p_p_wcl')]
 
       self$fertilizers <- torch::torch_tensor(as.matrix(fertilizers), device = device)
     },
@@ -63,7 +67,7 @@ createApusDataset <- function(farms = NULL, cultivations, fertilizers, fields_ma
       }
       t.fields <- transformFieldsToTensor(farms, self$device)
 
-      return(list(fields = t.fields, fertilizers = self$fertilizers))
+      return(list(fields = t.fields, fertilizers = self$fertilizers, fines = self$fines))
     },
 
     .length = function() {
@@ -73,7 +77,7 @@ createApusDataset <- function(farms = NULL, cultivations, fertilizers, fields_ma
 
 
   # Create torch dataset for apus -------------------------------------------
-  dataset <- apus_dataset(farms = farms, cultivations = cultivations, fertilizers = fertilizers, fields_max = fields_max, device = device)
+  dataset <- apus_dataset(farms = farms, cultivations = cultivations, fertilizers = fertilizers, fines = fines, fields_max = fields_max, device = device)
 
   return(dataset)
 }

R/model.R

@@ -106,14 +106,14 @@ createApusModel <- function(dataset.train, dataset.valid, width = 12, layers = 1
 
       cli::cli_progress_bar(paste0('Training model [', epoch, '/', epochs, ']'), total = dl.train$.length())
 
-      # For testing
+      # For developing
       # b <- dl.train$.iter()
       # b <- b$.next()
 
       # Forward pass
       optimizer$zero_grad()
       doses <- model(b$fields, b$fertilizers)
-      cost <- calculateCost(doses, b$fields, b$fertilizers)
+      cost <- calculateCost(doses, b$fields, b$fertilizers, b$fines)
 
       # Backward pass
       cost$backward()
@@ -133,13 +133,13 @@ createApusModel <- function(dataset.train, dataset.valid, width = 12, layers = 1
 
       cli::cli_progress_bar(paste0('Validating model [', epoch, '/', epochs, ']'), total = dl.valid$.length())
 
-      # For testing
+      # For developing
       # b <- dl.valid$.iter()
       # b <- b$.next()
 
       # Forward pass
       doses <- model(b$fields, b$fertilizers)
-      cost <- calculateCost(doses, b$fields, b$fertilizers)
+      cost <- calculateCost(doses, b$fields, b$fertilizers, b$fines)
 
       losses.validation <- c(losses.validation, cost$item())
 
@@ -157,7 +157,7 @@ createApusModel <- function(dataset.train, dataset.valid, width = 12, layers = 1
   return(model)
 }
 
-calculateCost <- function(doses, fields, fertilizers, reduce_batches = TRUE) {
+calculateCost <- function(doses, fields, fertilizers, fines, reduce_batches = TRUE) {
 
 
   # Check arguments ---------------------------------------------------------
@@ -168,12 +168,16 @@ calculateCost <- function(doses, fields, fertilizers, reduce_batches = TRUE) {
   module1 <- calculateCostModule1(doses, fields, fertilizers)
 
 
-  # Module 4: Revenue from harvested crops ----------------------------------
-  module4 <- calculateRevenueModule4(doses, fields, fertilizers)
+  # Module 5: Revenue from harvested crops ----------------------------------
+  module5 <- calculateRevenueModule5(doses, fields, fertilizers)
+
+
+  # Module 6: Penalty for exceeding legal limit -----------------------------
+  module6 <- calculatePenaltyModule6(doses, fields, fertilizers, fines)
 
 
   # Combine the modules -----------------------------------------------------
-  cost <- torch::torch_zeros(dim(module1)) + module1 - module4
+  cost <- torch::torch_zeros(dim(module1)) + module1 - module5 + module6
 
 
   # Convert to € / ha -------------------------------------------------------
@@ -211,8 +215,8 @@ calculateCostModule1 <- function(doses, fields, fertilizers) {
   return(module1)
 }
 
-# Module 4: Revenue from harvested crops  ------------------------------------
-calculateRevenueModule4 <- function(doses, fields, fertilizers) {
+# Module 5: Revenue from harvested crops  ------------------------------------
+calculateRevenueModule5 <- function(doses, fields, fertilizers) {
 
   # Calculate N dose per fields
   fertilizers.p_n_rt <- fertilizers[,,3]
@@ -268,9 +272,71 @@ calculateRevenueModule4 <- function(doses, fields, fertilizers) {
   fields.b_area <- fields[,,1]
   fields.b_lu_yield <- fields[,,8]
   fields.b_lu_price <- fields[,,9]
-  module4 <- fields.b_area *  fields.b_lu_yield * fields.b_lu_price * fields.d_realized
-  module4 <- torch::torch_sum(module4, dim = 2L)
+  module5 <- fields.b_area *  fields.b_lu_yield * fields.b_lu_price * fields.d_realized
+  module5 <- torch::torch_sum(module5, dim = 2L)
 
-  return(module4)
+  return(module5)
 }
 
+# Module 6: Penalties in case of exceeding legal limits -----------------------
+calculatePenaltyModule6 <- function(doses, fields, fertilizers, fines) {
+
+  # Calculate d_n_norm_man per field
+  fertilizers.p_n_rt <- fertilizers[,,3]
+  fertilizers.p_type_manure <- fertilizers[,,7]
+  fertilizers.p_n_manure <- fertilizers.p_n_rt * fertilizers.p_type_manure
+  fertilizers.p_n_manure <- torch::torch_unsqueeze(fertilizers.p_n_manure, 2)
+  fertilizers.p_n_manure <- torch::torch_repeat_interleave(fertilizers.p_n_manure, repeats = dim(doses)[2], dim =2)
+  fields.fertilizers.dose.n_manure <- doses * fertilizers.p_n_manure
+  fields.dose.n_manure <- torch::torch_sum(fields.fertilizers.dose.n_manure, dim = 3)
+  farms.dose.n_manure <- torch::torch_sum(fields.dose.n_manure, dim = 2)
+
+  fields.d_n_norm_man <- fields[,,6]
+  fields.b_area <- fields[,,1]
+  fine.d_n_norm_man <- fines[,1,2]
+  farms.d_n_norm_man <- torch::torch_sum(fields.b_area * fields.d_n_norm_man, dim = 2L)
+  farms.exceeding.d_n_norm_man <- torch::torch_relu(farms.dose.n_manure - farms.d_n_norm_man)
+  farms.penalty.d_n_norm_man <- farms.exceeding.d_n_norm_man * fine.d_n_norm_man
+
+  # Calculate d_n_norm per field
+  fertilizers.p_n_rt <- fertilizers[,,3]
+  fertilizers.p_n_wc <- fertilizers[,,4] # TODO Replace with p_n_wcl
+  fertilizers.p_n_workable <- fertilizers.p_n_rt * fertilizers.p_n_wc
+  fertilizers.p_n_workable <- torch::torch_unsqueeze(fertilizers.p_n_workable, 2)
+  fertilizers.p_n_workable <- torch::torch_repeat_interleave(fertilizers.p_n_workable, repeats = dim(doses)[2], dim =2)
+  fields.fertilizers.dose.n_workable <- doses * fertilizers.p_n_workable
+  fields.dose.n_workable <- torch::torch_sum(fields.fertilizers.dose.n_workable, dim = 3)
+  farms.dose.n_workable <- torch::torch_sum(  fields.dose.n_workable , dim = 2)
+
+  fields.d_n_norm <- fields[,,5]
+  fields.b_area <- fields[,,1]
+  fine.d_n_norm <- fines[,1,1]
+  farms.d_n_norm <- torch::torch_sum(fields.b_area * fields.d_n_norm, dim = 2L)
+  farms.exceeding.d_n_norm <- torch::torch_relu(farms.dose.n_workable - farms.d_n_norm)
+  farms.penalty.d_n_norm <- farms.exceeding.d_n_norm * fine.d_n_norm
+
+  # Calculate d_p_norm per field
+  fertilizers.p_p_rt <- fertilizers[,,5]
+  fertilizers.p_p_wcl <- fertilizers[,,8]
+  fertilizers.p_p_legal <- fertilizers.p_p_rt * fertilizers.p_p_wcl
+  fertilizers.p_p_legal <- torch::torch_unsqueeze(fertilizers.p_p_legal, 2)
+  fertilizers.p_p_legal <- torch::torch_repeat_interleave(fertilizers.p_p_legal, repeats = dim(doses)[2], dim =2)
+  fields.fertilizers.dose.p_legal <- doses * fertilizers.p_p_legal
+  fields.dose.p_legal <- torch::torch_sum(fields.fertilizers.dose.p_legal, dim = 3)
+  farms.dose.p_legal <- torch::torch_sum(  fields.dose.p_legal, dim = 2)
+
+  fields.d_p_norm <- fields[,,7]
+  fields.b_area <- fields[,,1]
+  fine.d_p_norm <- fines[,1,3]
+  farms.d_p_norm <- torch::torch_sum(fields.b_area * fields.d_p_norm, dim = 2L)
+  farms.exceeding.d_p_norm <- torch::torch_relu(farms.dose.p_legal - farms.d_p_norm)
+  farms.penalty.d_p_norm <- farms.exceeding.d_p_norm * fine.d_p_norm
+
+
+  # Combine the penalties
+  module6 <- farms.penalty.d_n_norm_man + farms.penalty.d_n_norm + farms.penalty.d_p_norm
+
+  return(module6)
+}
+
+

README.md

@@ -73,13 +73,17 @@ For the v1 version of `apus` we plan to develop to following features:
 
 * [ ] Import and export trained models
 * [ ] Include a trained base model 
+* [x] Add function to train model
 * [ ] Enable fine-tuning of (base) models
-* [ ] Include cost function for module 2: Cost of storing fertilizers
-* [ ] Include cost function for module 3: Cost of applying fertilizers
-* [ ] Include cost function for module 5: Penalties in case of exceeding legal limits
-* [ ] Include cost function for module 6: Cost of greenhouse gas emissions
-* [ ] Include realistic cultivation response curves from module 4
-* [ ] Add other nutrients and organic matter to module 4
+* [x] Include cost function for module 1: Purchase of fertilizers
+* [ ] Include cost function for module 2: Disposal of manure
+* [ ] Include cost function for module 3: Cost of storing fertilizers
+* [ ] Include cost function for module 4: Cost of applying fertilizers
+* [x] Include cost function for module 5: Revenue of harvest
+* [x] Include cost function for module 6: Penalties in case of exceeding legal limits
+* [ ] Include cost function for module 7: Cost of greenhouse gas emissions
+* [ ] Include realistic cultivation response curves from module 5
+* [ ] Add other nutrients than NPK to module 5
 * [ ] Add custom fertilizers
 * [ ] Add custom cultivations
 * [ ] Add details of the optimization to the result

data-raw/fertilizers.R

@@ -13,5 +13,12 @@ token <- ''
 fertilizers <- fread(paste0('https://raw.githubusercontent.com/AgroCares/pandex/main/data-raw/b_fp/b_fp_srm.csv?token=', token))
 setnames(fertilizers, colnames(fertilizers), tolower(colnames(fertilizers)))
 
+fertilizers[, p_type_manure := fifelse(p_type_manure, 1, 0)]
+fertilizers[, p_p_wcl := 1]
+fertilizers[p_type_compost == TRUE, p_p_wcl := 0.25]
+fertilizers[p_name_nl == 'Champost', p_p_wcl := 0.75]
+fertilizers[p_name_nl == 'Rundvee vaste mest', p_p_wcl := 0.75]
+
+
 # Export table ------------------------------------------------------------
 usethis::use_data(fertilizers, overwrite = TRUE, version = 3, compress = 'xz')

tests/testthat/test-001-dataset.R

@@ -1,10 +1,10 @@
 test_that("Create training dataset", {
   fields_max <- 5
-  dataset <- apus::createApusDataset(farms = NULL, cultivations = apus::cultivations, fertilizers = apus::fertilizers, fields_max = fields_max, device = 'cpu')
+  dataset <- apus::createApusDataset(farms = NULL, cultivations = apus::cultivations, fertilizers = apus::fertilizers, fines = apus::fines, fields_max = fields_max, device = 'cpu')
   expect_contains(class(dataset), 'apus_dataset')
   expect_equal(dataset$.length(), dataset$farms_count)
   expect_contains(class(dataset$.getitem(1)$fields), 'torch_tensor')
-  expect_setequal(names(dataset$.getitem(1)), c('fields', 'fertilizers'))
+  expect_setequal(names(dataset$.getitem(1)), c('fields', 'fertilizers', 'fines'))
   expect_equal(dim(dataset$.getitem(1)$fields), c(dataset$fields_max, 9))
   expect_false(identical(dataset$.getitem(1), dataset$.getitem(2)))
 })
@@ -17,11 +17,11 @@ test_that("Create validation/test dataset", {
   expect_contains(class(farms), 'data.table')
   expect_equal(nrow(farms), farms_count * fields_max)
 
-  dataset.valid <- apus::createApusDataset(farms = farms, cultivations = apus::cultivations, fertilizers = apus::fertilizers, fields_max = fields_max, device = 'cpu')
+  dataset.valid <- apus::createApusDataset(farms = farms, cultivations = apus::cultivations, fertilizers = apus::fertilizers, fines = apus::fines, fields_max = fields_max, device = 'cpu')
   expect_contains(class(dataset.valid), 'apus_dataset')
   expect_equal(dataset.valid$.length(), dataset.valid$farms_count)
   expect_contains(class(dataset.valid$.getitem(1)$fields), 'torch_tensor')
-  expect_setequal(names(dataset.valid$.getitem(3)), c('fields', 'fertilizers'))
+  expect_setequal(names(dataset.valid$.getitem(3)), c('fields', 'fertilizers', 'fines'))
   expect_equal(dim(dataset.valid$.getitem(1)$fields), c(fields_max, 9))
 
   dl <- torch::dataloader(dataset.valid, batch_size = farms_count)
@@ -30,7 +30,7 @@ test_that("Create validation/test dataset", {
 
   expect_contains(class(batch$fields), 'torch_tensor')
   expect_equal(dim(batch$fields), c(farms_count, fields_max, 9))
-  expect_equal(dim(batch$fertilizers), c(farms_count, nrow(apus::fertilizers), 6))
+  expect_equal(dim(batch$fertilizers), c(farms_count, nrow(apus::fertilizers), 8))
 })