rm(list=ls())

library(tidyverse)

# working directory must be defined (setwd()) 

###### importing data ------

eco<-read.table("data/20221107-1015-TG3-ecoTable.csv",header=T,sep=";")%>%
  mutate(scenario = factor(scenario_case,levels = c(0,1,2,4),
                           labels = c("Unprotected","Bark-scraping 400", "Bark-scraping all","Fence")))%>%
  filter(label != "latex")

stand<-read.table("data/20221107-1015-TG3-standTable.csv",
                  header=T,
                  sep=";")%>%
  mutate(scenario_stand = factor(scenario_case, levels = c(0,1,2,4),
                                 labels = c("Unprotected","Bark-scraping 400", "Bark-scraping all","Fence")))%>%
  group_by(siteIndex,bpr,scenario_case,lofflerModel,repetition)%>%
  mutate(isIntervenerResult = date == lag(date),
         isIntervenerResult = ifelse(is.na(isIntervenerResult),F,isIntervenerResult))

###### some functions -----
npv<-function(year,net,r,final_year){
  n<-final_year
  v = sum(net/(1+r)^year) * (1+r)^n/((1+r)^n-1)
  return(v)
}
getbillbook2<-function(stand,billbook){
  stumpage <- stand%>%
    group_by(date)%>%
    summarise(stumpage = max(stumpage))
  
  billboook2<-billbook%>%
    group_by(date)%>%
    summarise(income = sum(income),
              expense = sum(expense))%>%
    left_join(stumpage, by = "date")
  return(billboook2)
}
getnpvs<-function(billbook2){
  
  npv_results<-tibble()
  
  rotations<-unique(billbook2$date[billbook2$date>17])
  
  for(rotation in rotations){
    
    nets <- billbook2%>%
      mutate(net = ifelse(date<rotation,income-expense,ifelse(date==rotation,stumpage-expense,0)))
    
    npv_results<-npv_results%>%
      bind_rows(tibble("rotation"=rotation,"discountRate"=0.01,"npv"=npv(nets$date,nets$net,r = 0.01,final_year = rotation)))%>%
      bind_rows(tibble("rotation"=rotation,"discountRate"=0.02,"npv"=npv(nets$date,nets$net,r = 0.02,final_year = rotation)))%>%
      bind_rows(tibble("rotation"=rotation,"discountRate"=0.03,"npv"=npv(nets$date,nets$net,r = 0.03,final_year = rotation)))%>%
      bind_rows(tibble("rotation"=rotation,"discountRate"=0.04,"npv"=npv(nets$date,nets$net,r = 0.04,final_year = rotation)))
    
  }
  
  return(npv_results)
}

#### parameter ----
trituPrice <- 5 ### the debased volume was considered to have no value

#### manipulating data -----

## calculer le volume des purges
tmp<-stand%>%
  mutate(vol_debased = vha - vha_eco)%>%
  mutate(value_debased = vol_debased * trituPrice)%>%
  mutate(stumpage2 = stumpage + value_debased)

eco_purge<-tmp%>%
  group_by(siteIndex, bpr, scenario_case, scenario_stand, lofflerModel, repetition, date)%>%
  mutate(intervention_year = any(isIntervenerResult))%>%
  filter(intervention_year)%>%
  summarise(vha_ecl_grume = first(vha_eco) - last(vha_eco),
            vha_ecl_total = first(vha) - last(vha),
            vha_ecl_tritu = vha_ecl_total - vha_ecl_grume,
            value_ecl_tritu = vha_ecl_tritu * trituPrice)%>%
  mutate(label = "thinning_debased",
         income= value_ecl_tritu,
         volume = vha_ecl_tritu,
         expense = 0)%>%
  select(siteIndex, bpr, scenario_case, lofflerModel, repetition, date, label, income, expense, volume, scenario = scenario_stand)%>%
  filter(volume > 0.001)

### combine existing simulation
nestedStand<-tmp%>%
  group_by(siteIndex, bpr, scenario_case, lofflerModel, repetition)%>%
  nest()%>%
  rename(stand = data)

nestedEco <- eco%>%
  bind_rows(eco_purge)%>%
  ungroup()%>%
  group_by(siteIndex, bpr, scenario_case, scenario, lofflerModel, repetition)%>%
  nest()%>%
  rename(billbook = data)

#### !!! Takes time !!!!
results <- nestedEco %>%
 left_join(nestedStand, by = c("siteIndex", "bpr", "scenario_case", "lofflerModel", "repetition"))%>%
  mutate(billbook2 = map2(stand, billbook, getbillbook2))%>%
  mutate(npvs = map(billbook2, getnpvs))

results_npv<-results%>%
  unnest(cols = c(npvs))

results_maxnpv<-results_npv%>%
  group_by(siteIndex, bpr, scenario_case, scenario, lofflerModel, repetition, discountRate)%>%
  summarise(maxnpv = max(npv),
            rotation = rotation[npv == maxnpv])

getStandRotation<-function(stand,rotation){
  
  stand2<-stand%>%
    filter(date<=rotation)
  
  lastline<-slice(stand2,n())%>%
    mutate_at(vars(dmoy:peeled_prop), ~0)
  
  stand2<-stand2%>%
    bind_rows(lastline)
  
}

results2<-results%>%
  left_join(filter(results_maxnpv,discountRate == 0.03), 
            by = c("siteIndex", "bpr", "scenario_case", "repetition", "scenario"))%>%
  mutate(stand2 = map2(stand,rotation, getStandRotation)) #slow 


#### export data  ----
save(eco,stand,results,results2,results_maxnpv, file="tmp/prepared_data.Rdata")

















