RProjects/R/project/ecow_calculations.R

suppressPackageStartupMessages({ 
  library(roxygen2)
  library(dplyr)
  library(lubridate)
  library(purrr)
  library(stringr)
  library(DBI)
  library(RPostgres)
  library(here)
  library(httr)
  library(jsonlite)
})
source(here::here("R/common/ws_client.R"))
source(here::here("R/project/preprocessing.R"))
source(here::here("R/project/postprocessing.R"))

#' Fonction globale de mise à jour des indicateurs eCow pour vaches, taureaux et lignées
#' Stockage des données directement en base
#' step : 1 si calcul des vaches uniquement, 2 si vaches et taureaux, 3 global
#' @param cheptel character. Numéro du cheptel avec le FR devant
calcul_ecow_by_chep <- function(cheptel, step = 3) {
  t0 <- Sys.time()
  
  ######################################     Préparation des données       ##########################################

  # Récupère les paramètres de pondération
  params_ponderation <- get_params_ponderation(cheptel)
  # Récupère la liste des certificats zoo issue de Doli
  czhbc <- data.frame(ANIM = get_cztotaux())
  
  message("Début import des données")
  cheptel_ecow <- get_cheptel_ecow(cheptel)
  
  # Vaches actives ayant déjà eu une fin de gestation
  vaches <- cheptel_ecow$vaches
  
  if (length(vaches) == 0) {
    t1 <- Sys.time()
    message("Arrêt du traitement, aucunes vaches à traiter. Temps d'exécution : ", round(difftime(t1, t0, units = "secs"), 2), " sec")
    return(invisible(NULL))  
  }
  
  # Récupère les taureaux : tous les pères des vaches actives ou de tous les veaux des vaches actives
  taureaux <- cheptel_ecow$taureaux
  # Ajout du nom du père
  vaches <- vaches %>%
    dplyr::left_join(
      taureaux %>% select(anim, nom_pere = nom),
      by = c("pereGenetique" = "anim")
    )
  # Récupère les produits ajout des données manquantes
  produits_cheptel <- add_data_ecow(cheptel_ecow$produits, czhbc)
  
  
  ###################################################################################################################
  ####################################     Calcul ecow pour les vaches       ########################################
  ###################################################################################################################
  message("Début partie vaches : ", nrow(vaches), " vaches à traiter.")
  
  # Ajout info donneuses
  donneuses <- cheptel_ecow$donneuses
  vaches$donneuse <- vaches$anim %in% donneuses

  # Récupère les produits des vaches actives du cheptel et les embryons portés dans le cheptel
  produits_vaches <- produits_cheptel %>%
    filter(numeroMipg %in% vaches$anim)
  
  # Récupère les effets cheptel # TODO A REVOIR AVEC LAURENA
  effets <- get_effets_cheptel(produits_vaches)
  rapport_MF <- effets$rapport_MF

  prod_vaches_corr <- apply_effet_chep(produits_vaches, effets$effets_chep, effets$rapport_MF)
  
  # On récupère les coefficients de pondération pour les pointages au sevrage 
  pps <- params_ponderation$pointage$sevrage
  
  # Ajout à la table vache des informations synthétisées de leur veaux
  synth_brute_vaches <- get_synth_prod_parent(vaches, prod_vaches_corr, params_ponderation)
  
  # Normalisation et calcul des notes carrières
  synth_norm <- get_note_carriere(synth_brute_vaches, params_ponderation)
  synth_vaches <- synth_norm$synthese

  # calcul des stats, valeurs extremes et references pour la normalisation
  stats_chep <- synth_norm$stats_chep

  ####################################     Calcul des notes campagnes       ########################################
  ####### En réalité on travaille sur les rangs de velages, ce qui correspond dans 99% des cas aux campagnes #######
  
  # ==============================
  # 1. Aggrégation campagnes
  # ==============================
  synth_prod_vache <- prod_vaches_corr %>%
    group_by(numeroMipg, dateNaiss, rangVelageMipg) %>%
    summarise(
      ivv    = first(ivv1),  # ------------------------------------------------------------------------- TODO vraiment pas sure, à valider
      pn_c   = round(mean(pn_corr, na.rm = TRUE), 1),
      txvf   = round(mean(conditionNaiss %in% c('1','2')) * 100, 1),
      txm    = round(mean(sexe == '1') * 100, 1),
      ptgp   = round(mean(pps$devmus * dmSevrage + pps$devsqe * dsSevrage + pps$af * afSevrage, na.rm = TRUE), 1),
      p120_c = round(mean(pat120Corrige, na.rm = TRUE), 1),
      p210_c = round(mean(pat210Corrige, na.rm = TRUE), 1),
      prol   = n() * 100,  # ------------------------------- TODO a tester parce que je pense qu'il faudrait tous les produits d'une vache
      mort   = round(mean(mortsev == "O" | mortnat == "O") * 100, 1),
      pere   = first(pereGenetique),
      cheptel = first(cheptelNaiss),
      
      # noms des veaux pour simplifier affichage
      produits = list(
        pmap(
          list(nom = nom, anim = anim, sexe = sexe, mortnat = mortnat, mortsev = mortsev, NBPRODIPG = NBPRODIPG, embryon = embryon),
          \(nom, anim, sexe, mortnat, mortsev, NBPRODIPG, embryon) list(nom = nom, anim = anim, sexe = sexe, mortnat = mortnat, mortsev = mortsev, NBPRODIPG = NBPRODIPG, embryon = embryon)
        )
      ),
      .groups = "drop"
    ) 
  
  # Récupère le nom du père
  synth_prod_vache <- synth_prod_vache %>%
    dplyr::left_join(
      taureaux %>% select(anim, nom_pere = nom),
      by = c("pere" = "anim")
    )
  
  stats_camp <- get_stats_tbl(
    tab = synth_prod_vache,
    nom_tab = "synth_prod_vache",
    cols = c( 'ivv', 'prol', 'mort', 'txvf', 'txm', 'ptgp', 'pn_c', 'p120_c', 'p210_c')
  )
  
  stats_chep <- rbind(stats_chep, stats_camp)
  
  pond_camp_fin <- unlist(params_ponderation$campagne$final)
  pond_camp_ahp <- unlist(params_ponderation$campagne$AHPtech)
  
  # ==============================
  # 2. Normalisation complète
  # ==============================
  synth_prod_vache_n <- synth_prod_vache %>%
    mutate(
      # pn normalisé
      pn_n = case_when(
        is.na(pn_c) ~ NA_real_,
        pn_c >= 40 & pn_c <= 50 ~ 1,
        pn_c <= 22 | pn_c >= 68 ~ 0,
        pn_c > 22 & pn_c < 40 ~ round(0.056 * (pn_c - 22), 3),
        TRUE ~ round(1 - 0.056 * (pn_c - 50), 3)
      ),
      
      # 5 normalisations linéaires
      txvf_n = round(1 - abs(stats_chep$max[stats_chep$var=="synth_prod_vache$txvf"] - txvf) /
                       abs(diff(range(stats_chep[stats_chep$var=="synth_prod_vache$txvf",c("min","max")]))), 3),
      
      txm_n = round(1 - abs(stats_chep$max[stats_chep$var=="synth_prod_vache$txm"] - txm) /
                      abs(diff(range(stats_chep[stats_chep$var=="synth_prod_vache$txm",c("min","max")]))), 3),
      
      ptgp_n = round(1 - abs(stats_chep$max[stats_chep$var=="synth_prod_vache$ptgp"] - ptgp) /
                       abs(diff(range(stats_chep[stats_chep$var=="synth_prod_vache$ptgp",c("min","max")]))), 3),
      
      p120_n = round(1 - abs(stats_chep$max[stats_chep$var=="synth_prod_vache$p120_c"] - p120_c) /
                       abs(diff(range(stats_chep[stats_chep$var=="synth_prod_vache$p120_c",c("min","max")]))), 3),
      
      p210_n = round(1 - abs(stats_chep$max[stats_chep$var=="synth_prod_vache$p210_c"] - p210_c) /
                       abs(diff(range(stats_chep[stats_chep$var=="synth_prod_vache$p210_c",c("min","max")]))), 3),
      
      # prol normalisé
      prol_n = case_when(
        is.na(prol) ~ NA_real_,
        prol == 100 ~ 0.8,
        TRUE ~ 1
      ),
      
      # mortalité
      mort_n = round(exp(-0.031 * mort), 3),  # TODO sciender en mortsev et mortnat
      
      # IVV normalisé
      ivv_n = case_when(
        is.na(rangVelageMipg) | rangVelageMipg == 1 | is.na(ivv) ~ NA_real_,
        
        # ravelamere == 2
        rangVelageMipg == 2 & ivv > 460 ~ 0,
        rangVelageMipg == 2 & ivv < 390 ~ 1,
        rangVelageMipg == 2 ~ round(1 - abs(390 - ivv)/abs(390 - 460), 3),
        
        # autres ravelamere
        ivv > 435 ~ 0,
        ivv < 365 ~ 1,
        TRUE ~ round(1 - abs(365 - ivv)/abs(365 - 435), 3)
      )
    ) %>%
    
    # ==============================
    # 3. Score final ecowcamp
    # ==============================
    rowwise() %>%
      mutate(
        perf = list(c_across(c(
          ivv_n, mort_n, p120_n, p210_n, pn_n,
          prol_n, ptgp_n, txm_n, txvf_n
        ))),
        pond = sum(
          pond_camp_fin[!(is.na(perf) | is.nan(perf))]
        ),
        somme = sum(
          perf[!(is.na(perf) | is.nan(perf))] *
            pond_camp_ahp[!(is.na(perf) | is.nan(perf))]
        ),
        SOMME_tot = somme / pond * 10,
        
        ecowcamp = ifelse(
          is.na(ptgp_n) & is.na(p120_n) & is.na(p210_n), # Si pas de pointage, on réduit la note
          NA,
          round(SOMME_tot * 10, 0)
        ),
        produits = toJSON(produits, auto_unbox = TRUE)
      ) %>%
      ungroup()
  
  # remplissage de la table vaches avec les notes campagnes
  
  # 1) Moyenne ecowcamp par mère
  moy_camp <- synth_prod_vache_n %>%
    filter(ecowcamp > 10) %>%
    group_by(numeroMipg) %>%
    summarise(moyecowcamp = round(mean(ecowcamp, na.rm = TRUE), 1), .groups = "drop")
  
  # 3) Fusion + transformations
  v_camp <-synth_vaches %>%
    dplyr::left_join(moy_camp, by = c("anim" = "numeroMipg")) %>%
    mutate(
      rg_camp = as.integer(rank(1 / moyecowcamp, na.last='keep'))
    )
  
  message("Enregistrement données vaches")
  
  # Enregistrement des données des vaches en base
  v_camp <- v_camp %>%
    mutate(across(where(is.numeric), ~ trunc(.x * 100) / 100))
  
  save_data_vaches(v_camp)
  
  # Enregistrement des données de campagne en base
  synth_prod_vache_n <- synth_prod_vache_n %>%
    mutate(across(where(is.numeric), ~ trunc(.x * 100) / 100))
  save_data_campagne(synth_prod_vache_n)

  if (step == 1) {
    t1 <- Sys.time()
    message("Fin du traitement. Temps d'exécution : ", round(difftime(t1, t0, units = "secs"), 2), " sec")
    return(invisible(NULL))  
  }
  
  ###################################################################################################################
  ####################################     Calcul ecow pour les taureaux     ########################################
  ###################################################################################################################
  message("Début partie taureaux")
  # Recupère les produits des taureaux
  produits_taureaux <- produits_cheptel %>%
    filter(pereGenetique %in% taureaux$anim & embryon != 'O')
  
  # Récupères les filles des taureaux 
  filles_taureaux <- produits_taureaux %>%
    filter(sexe == 2)
  
  # Petits produits issus des filles des taureaux 
  pprod_filles_taureaux <- add_data_ecow(cheptel_ecow$petits_produits, czhbc)
  
  # Calcul de la synthèse par fille de leur produits
  synth_filles_taureaux <- get_synth_prod_parent(filles_taureaux, pprod_filles_taureaux, params_ponderation)
  
  # Calcul des stats par pere
  stats_prod_directe <- get_stats_parent(produits_taureaux, pereGenetique)

  stats_filles <- get_stats_filles(synth_filles_taureaux, "fillestot_", pereGenetique, actif = FALSE)

  stats_filles_act <- get_stats_filles(synth_filles_taureaux, "fillesact_", pereGenetique, actif = TRUE)
  
  grp_filles_et_act <- stats_filles %>%
    dplyr::left_join(stats_filles_act, by = "pereGenetique") %>%
    mutate(
      pctfillesact_avecprod = round( fillesact_nbavecprod / fillestot_nbavecprod * 100, 1 )
    )
  
  stats_filles_renouv <- produits_cheptel %>%
    filter(
      sexe == "2",
      nbFinGestation == 0 # TODO valider avec Lauréna que c'est bien ce champ là qu'on veut 
    ) %>%
    group_by(pereGenetique) %>%
    summarise(
      nbfilles_renouv = n(),
      .groups = "drop"
    )
      
  stats_taureaux <- stats_prod_directe %>%
    dplyr::left_join(grp_filles_et_act,    by = "pereGenetique") %>%
    dplyr::left_join(stats_filles_renouv,  by = "pereGenetique") %>%
    dplyr::left_join(taureaux %>% select(anim, nom, dateNaiss, nomCheptelNaiss),  by = c("pereGenetique" = "anim")) %>%
    mutate(
      nom = replace(nom, is.na(nom), ""),
      cheptel = cheptel,
      across(where(is.numeric), ~ trunc(.x * 100) / 100)
      )

  message("Enregistrement données taureaux")
  save_data_taureau(stats_taureaux)
  
  if (step == 2) {
    t1 <- Sys.time()
    message("Fin du traitement. Temps d'exécution : ", round(difftime(t1, t0, units = "secs"), 2), " sec")
    return(invisible(NULL))  
  }

  ###################################################################################################################
  ####################################     Remontee des lignees femelles     ########################################
  ###################################################################################################################
  
  fondatrices <- cheptel_ecow$fondatrices
  descendants <- add_data_ecow(cheptel_ecow$descendants, czhbc)

  vaches_lignees <- descendants %>% filter(anim %in% descendants$mereGenetique)
  produits_lignees <- descendants %>% filter(mereGenetique %in% vaches_lignees$anim)

  synth_vaches_lignees <- get_synth_prod_parent(vaches_lignees, produits_lignees, params_ponderation)

  # calcul des stats par fondatrice, ne gardant que celles ayant plus de 5 descendants dans le cheptel_________________

  stats_prod <- get_stats_parent(descendants, fondatrice) # produits lignées ou descendants ????
  
  stats_fem_tot <- get_stats_filles(synth_vaches_lignees, "femtot_", fondatrice, actif = FALSE)
  
  stats_fem_act <- get_stats_filles(synth_vaches_lignees, "femact_", fondatrice, actif = TRUE)

  stats_renouv <- descendants %>%
    filter(
      is.na(NBPRODIPG),# TODO valider avec Lauréna que c'est bien ce champ là qu'on veut 
      sexe == "2",
      is.na(dateSortDetenteur)
    ) %>%
    group_by(fondatrice) %>%
    summarise(
      nbfilles_renouv = n(),
      .groups = "drop"
    )
  
  stats_lignees <- stats_prod %>%
    dplyr::left_join(stats_fem_tot,  by = "fondatrice") %>%
    dplyr::left_join(stats_fem_act,  by = "fondatrice") %>%
    dplyr::left_join(stats_renouv,   by = "fondatrice") %>%
    dplyr::left_join(fondatrices %>% select(anim, nom, dateNaiss, nomCheptelNaiss),  by = c("fondatrice" = "anim")) %>%
    mutate(
      pctfem_avecprod =
        round(femtot_nbavecprod / nb_femelles* 100, 1),
      pctfemact_avecprod =
        round(femact_nbavecprod / nb_femelles * 100, 1),
      cheptel = cheptel
    )

  message("Enregistrement données lignées")
  save_data_lignees(stats_lignees)
  t1 <- Sys.time()
  message("Fin du traitement. Temps d'exécution : ", round(difftime(t1, t0, units = "secs"), 2), " sec")
}