require(tidyverse)
require(lpSolve)
nutrient_uptake <- function(S1=NA, S2=NA, d1=NA, a1=NA, d2=NA, a2=NA, r1=NA, r2=NA) {	
  # N, P and K uptakes based on QUEFTS
  if (S1 < r1 + ((S2 - r2) * a2 / d1)) {
    uptakeX_givenY = S1
  } else if (S1 > r1 + ((S2 - r2) * (2 * d2 / a1 - a2 / d1))) {
    uptakeX_givenY = r1 + (S2 - r2) * (d2 / a1)
  } else {
    uptakeX_givenY = S1 - 0.25 * (S1 - r1 - (S2 - r2) * (a2 / d1))^2 / ((S2 - r2) * (d2 / a1 - a2 / d1))
  }
  # Nutrient uptake given availability of other nutrient
  return(uptakeX_givenY)
}

yield_nutrients_combined <- function(U1=NA, d1=NA, a1=NA, Y2A=NA, Y2D=NA, Y3D=NA, r1=NA){	
  # Determine which nutrient limited to force yield being the lowest.
  YxD = min(Y2D, Y3D)	
  # If the uptake of one of the nutrients, and therefore the yield associated with that 
  # nutrient, is zero the overall yield is also zero.
  if (U1 == 0 || YxD == 0) {
    Y12 = 0
  }else{
    Y12 = Y2A + (2 * (YxD - Y2A) * (U1 - r1 - Y2A / d1)) / (YxD / a1 - Y2A / d1) -
      (YxD - Y2A) * (U1 - r1 - Y2A / d1)^2 / (YxD / a1 - Y2A / d1)^2
  }
  # Return the calculated yield based on the uptake of nutrients 1 and 2
  return(Y12)
}

water_dependent_nutrient_uptake <- function(S1=NA, WLY=NA, d1=NA, a1=NA, r1=NA) {	
  if (S1 < r1 + WLY / d1) {
    uptakeX_givenWater = S1    
  } else if (S1 > r1 + 2*WLY/a1 - WLY/d1) {
    uptakeX_givenWater = WLY / a1    
  } else {
    uptakeX_givenWater = S1 - 0.25 * (S1 - r1 - WLY/d1)^2 / (WLY / a1 - WLY / d1)    
  }
  
  return(uptakeX_givenWater)
}


NUE <- function(HI, CmaxNroots=2.5, CminNroots=0.9, CmaxNtops=3.2, CminNtops=1.4,
                CmaxProots=0.6, CminProots=0.1, CmaxPtops=1, CminPtops=0.13, 
                CmaxKroots=4.6, CminKroots=1.1, CmaxKtops=4.6, CminKtops=0.85 ){
  aN = round(1000 * HI/(HI * CmaxNroots + (1 - HI) * CmaxNtops), digits=0)
  dN = round(1000 * HI/(HI * CminNroots + (1 - HI) * CminNtops), digits=0)
  
  aP = round(1000 * HI/(HI * CmaxProots + (1 - HI) * CmaxPtops), digits=0)
  dP = round(1000 * HI/(HI * CminProots + (1 - HI) * CminPtops), digits=0)
  
  aK = round(1000 * HI/(HI * CmaxKroots + (1 - HI) * CmaxKtops), digits=0)
  dK = round(1000 * HI/(HI * CminKroots + (1 - HI) * CminKtops), digits=0)
  return(data.frame(aN=aN, dN=dN,aP=aP,dP=dP,aK=aK,dK=dK))
}

NUE(HI=0.75)


## Maize: aN=29	; dN=74	; aP=95	; dP=476; aK=38	; dK=143	; rN=5	; rP=0.4	, rK=2
## Cassava: aN=41, dN=96, rN=0, aP=233, dP=588, rP=0, aK=34, dK=161, rK=0)
##Yield_S <- function(SN, SP, SK, WLY, aN=41, dN=96, rN=0, aP=233, dP=588, rP=0, aK=34, dK=161, rK=0){ ## cassava
Yield_S <- function(SN, SP, SK, WLY, aN=175, dN=435, rN=0, aP=625, dP=4348, rP=0, aK=109, dK=513, rK=0){
  ## uptake of one nutrient conditioned by availablility of one other nutrient or water, rule of minimum
  UNP <- nutrient_uptake(S1 = SN, S2 = SP, d1 = dN, a1 = aN, d2 = dP, a2 = aP, r1 = rN, r2 = rP)
  UNK <- nutrient_uptake(S1 = SN, S2 = SK, d1 = dN, a1 = aN, d2 = dK, a2 = aK, r1 = rN, r2 = rK)
  UNW <- water_dependent_nutrient_uptake(S1 = SN, WLY = WLY, d1 = dN, a1 = aN, r1 = rN)
  UN <- min(UNP, UNK, UNW)
  
  UPN <- nutrient_uptake(S1 = SP, S2 = SN, d1 = dP, a1 = aP, d2 = dN, a2 = aN, r1 = rP, r2 = rN)
  UPK <- nutrient_uptake(S1 = SP, S2 = SK, d1 = dP, a1 = aP, d2 = dK, a2 = aK, r1 = rP, r2 = rK)
  UPW <- water_dependent_nutrient_uptake(S1 = SP, WLY = WLY, d1 = dP, a1 = aP, r1 = rP)
  UP <- min(UPN, UPK, UPW)				
  
  UKN <- nutrient_uptake(S1 = SK, S2 = SN, d1 = dK, a1 = aK, d2 = dN, a2 = aN, r1 = rK, r2 = rN)
  UKP <- nutrient_uptake(S1 = SK, S2 = SP, d1 = dK, a1 = aK, d2 = dP, a2 = aP, r1 = rK, r2 = rP)
  UKW <- water_dependent_nutrient_uptake(S1 = SK, WLY = WLY, d1 = dK, a1 = aK, r1 = rK)
  UK <- min(UKN, UKP, UKW)		
  
  ## yield base don uptake of Nitrogen at dilution or accumulation
  YNA <- max((UN - rN), 0) * aN
  YND <- max((UN - rN), 0) * dN
  YPA <- max((UP - rP), 0) * aP
  YPD <- max((UP - rP), 0) * dP
  YKA <- max((UK - rK), 0) * aK
  YKD <- max((UK - rK), 0) * dK	
  
  ## yield as defined by uptake of one nutrient conditioned to max and min yield based on secnd nutrient not to be constrined by a third nutrient
  YNP <- yield_nutrients_combined(U1 = UN, d1 = dN, a1 = aN, Y2A = YPA, Y2D = YPD, Y3D = YKD, r1 = rN)
  YNK <- yield_nutrients_combined(U1 = UN, d1 = dN, a1 = aN, Y2A = YKA, Y2D = YKD, Y3D = YPD, r1 = rN)
  YPN <- yield_nutrients_combined(U1 = UP, d1 = dP, a1 = aP, Y2A = YNA, Y2D = YND, Y3D = YKD, r1 = rP)
  YPK <- yield_nutrients_combined(U1 = UP, d1 = dP, a1 = aP, Y2A = YKA, Y2D = YKD, Y3D = YND, r1 = rP)
  YKN <- yield_nutrients_combined(U1 = UK, d1 = dK, a1 = aK, Y2A = YNA, Y2D = YND, Y3D = YPD, r1 = rK)
  YKP <- yield_nutrients_combined(U1 = UK, d1 = dK, a1 = aK, Y2A = YPA, Y2D = YPD, Y3D = YND, r1 = rK)
  
  # Make sure the nutrient limited yields do not exceed the maximum possible yield = WLY
  YNPc <- min(c(YNP, YND, YPD, YKD, WLY))
  YNKc <- min(c(YNK, YND, YPD, YKD, WLY))
  YPNc <- min(c(YPN, YND, YPD, YKD, WLY))
  YPKc <- min(c(YPK, YND, YPD, YKD, WLY))
  YKNc <- min(c(YKN, YND, YPD, YKD, WLY))
  YKPc <- min(c(YKP, YND, YPD, YKD, WLY))
  
  #Final estimate
  YEc <- mean(c(YNPc, YNKc, YPNc, YPKc, YKNc, YKPc))
  return(YEc)				
}




## get yield estimates using QUEFTS and giving in SN, SP, SK and WLY.  get TSS for the estimates of yield and NOT blup. 
## SN, SP and SK are supposed to be the sum of what the soil originally has and added fertilizer. 
## added fertilizer is provided as it is defined in NOT trials and INS is the one to be defined by iterating on while controling the output
## by the squared yield difference between simulated and measured from NOT  (Blups of NOT)
optim_INS <- function(INS, addedFertilizer, YieldMatrix, WLY){	
  
  SoilN <- INS[1]
  SoilP <- INS[2]
  SoilK <- INS[3]
  
  #Yield_Control <- Yield_S(SN = addedFertilizer[1,1] + SoilN , SP = addedFertilizer[1,2] + SoilP, SK = addedFertilizer[1,3] + SoilK, WLY = WLY)
  Yield_pk <- Yield_S(SN = (addedFertilizer[2,1]*RecoveryFraction$rec_N) + SoilN , 
                      SP = (addedFertilizer[2,2]*RecoveryFraction$rec_P) + SoilP, 
                      SK = (addedFertilizer[2,3]*RecoveryFraction$rec_K) + SoilK, WLY = WLY)
  
  Yield_nk <- Yield_S(SN = (addedFertilizer[3,1]*RecoveryFraction$rec_N) + SoilN ,
                      SP = (addedFertilizer[3,2]*RecoveryFraction$rec_P) + SoilP, 
                      SK = (addedFertilizer[3,3]*RecoveryFraction$rec_K) + SoilK, WLY = WLY)
  
  Yield_np <- Yield_S(SN = (addedFertilizer[4,1]*RecoveryFraction$rec_N) + SoilN , 
                      SP = (addedFertilizer[4,2]*RecoveryFraction$rec_P) + SoilP, 
                      SK = (addedFertilizer[4,3]*RecoveryFraction$rec_K) + SoilK, WLY = WLY)
  #Yield_halfnpk <- Yield_S(SN = addedFertilizer[5,1] + SoilN , SP = addedFertilizer[5,2] + SoilP, SK = addedFertilizer[5,3] + SoilK, WLY = WLY)
  
  
  #control_SS <- (Yield_Control - YieldMatrix$NOT_yield[1])^2
  pk_SS <- (Yield_pk - YieldMatrix$NOT_yield[2])^2
  nk_SS <- (Yield_nk - YieldMatrix$NOT_yield[3])^2
  np_SS <- (Yield_np - YieldMatrix$NOT_yield[4])^2
  #halfnpk_SS <- (Yield_halfnpk - YieldMatrix$NOT_yield[5])^2
  
  TSS <- sum( pk_SS, nk_SS, np_SS)	
  return(TSS)	
}








#######################################################################
## FR & SP
#######################################################################
#'  @param fertilizers: data frame with type, N_cont, P_cont, K_cont, price. Price is per kg of fertilizer
#'  @param NG_CY_Fertdata: data frame with lat,long,fert_N,fert_P,fert_K, water_limited_yield, CurrentYield,location, pl_Date,zone, harvestDay, harvestmonth, daysOnField
#'  @param SoilData: data frame with lat,long,soilN,soilP,soilK,rec_N, rec_P, rec_K, rel_N,rel_P,rel_K
#'  @param rootUP: a price of 1 tonne of cassava in freshwt. It is used as freshwt price, after QUEFTS give drywt root yield (kg/ha) it is converted to freshwt in tonne/ha and this price is then used
#'  @param areaHa is area of land in ha
#'  @param WLY_FertRecom is the file with WLY, CY and fert recom for WLY for all lcoations. WLY and CY are in DM kg/ha
#'  @return a data frame with NPK elemental amount required for target yield with total cost, net revenu,and list of fertilizers with the 
#' amount required in kg for the user defined land size

# PD = "2018-03-01"; HD = "2019-05-31"; lat = 10.024; lon = 4.025; country = "NG"; cassUW = 1000; cassUP = 12000; maxInv = 72000;
# SoilData = SoilGridData_NG , userName = "acai cassava", userPhoneCC = 254, userPhoneNr = 702974480,userEmail = "acai.akilimo@gmail.com",
# cassPD = "roots"
#' 
#' @param dss 
#' @returnType 
#' @return 
#' 
#' @author Meklit
#' @export
getsupply <- function(dss){
  supply <- data.frame(lat=dss$lat, long=dss$long, rel_N=dss$rel_N, rel_P=dss$rel_P, rel_K=dss$rel_K, SN=dss$soilN, SP=dss$soilP, SK=dss$soilK, water_limited_yield = dss$water_limited_yield,
                       aN=dss$aN, dN=dss$dN, aP=dss$aP, dP=dss$dP, aK=dss$aK, dK=dss$dK, rN=dss$rN, rP=dss$rP, rK=dss$rK, max_yield=dss$max_yield,  tolerance=dss$tolerance,
                       WLY = dss$water_limited_yield)
}


#' The soil NPK as obtained from randdom forest model
#' @param zone, is used to define the varieties and HI to get NUE. Lake zone (the default) is proxy for Mkobozi and E & S zone is Kiroba  
#' @param Queft_Input_Data: per lat and long, crop param and soil param, water limited yield, fertlizer recommendation
#' @return 
#' 
#' @author Meklit
#' @export
QUEFTS_WLY_CY <- function(SoilData=SoilData, country=country, wlyd=wlyd){	
  #wly_plDate <- wly_data[wly_data$plantingDate == pl_Date, c("lat", "long", "wly_KgHa")]
  wly_plDate <- wlyd[,  c("lat", "long", "water_limited_yield")]
  
  # colnames(wly_plDate) <- c("lat", "long", "water_limited_yield")	
  Quefts_Input_Data_wly <- merge(SoilData, wly_plDate, by=c("lat", "long"))
  
  ## HI: Median for Nigeria=0.55 and Tanzania=0.52. Q3, Nigeria=0.63 and Tanzania=0.61
  if(country == "NG"){
    crop_param <- cbind(NUE(HI=0.55), data.frame(rN=0, rP=0, rK=0, max_yield=Quefts_Input_Data_wly$water_limited_yield, tolerance=0.01))
  }else{
    crop_param <- cbind(NUE(HI=0.55), data.frame(rN=0, rP=0, rK=0, max_yield=Quefts_Input_Data_wly$water_limited_yield, tolerance=0.01))
  }
  
  ## 1. get soil nutrient supply
  Queft_Input_Data_Var <- cbind(Quefts_Input_Data_wly, crop_param)	
  supply <- getsupply(Queft_Input_Data_Var) ## to get yield at zero input level
  
  
  ## 2. Current yield: 
  actualUptake <- merge(supply,ddply(supply,.(lat, long), actual_uptake_tool), by=c("lat","long"))
  minmax_Yield <-  merge(actualUptake, ddply(actualUptake,.(lat, long), max_min_yields_tools), by=c("lat","long"))
  Current_Yield <- ddply(minmax_Yield,.(lat, long), final_yield_tools)## yield at zero input
  colnames(Current_Yield) <- c("lat", "long", "CurrentYield")
  Yield_Fertilizer <- merge(wly_plDate, Current_Yield, by=c("lat", "long"))
  return(Yield_Fertilizer$CurrentYield)
}


QUEFTS_WLYCY_Potato <- function(SoilData_onePix){	
 
  
  # crop_param <- data.frame(aN=41, dN=96, aP=233, dP=588,  aK=34, dK=161, rN=0,rP=0, rK=0,
  #                          max_yield = SoilData_onePix$water_limited_yield, tolerance = 0.01 )
  # 
  
  crop_param <- data.frame(aN=175, dN=435, aP=625, dP=4348,  aK=109, dK=513, rN=0,rP=0, rK=0, 
                           max_yield = SoilData_onePix$water_limited_yield, tolerance = 0.01 )
  
  ## 1. get soil nutrient supply
  dss <- cbind(SoilData_onePix, crop_param)	
  
  supply <- data.frame(rel_N=dss$rel_N, rel_P=dss$rel_P, rel_K=dss$rel_K, SN=dss$soilN, SP=dss$soilP, SK=dss$soilK, water_limited_yield = dss$water_limited_yield,
             aN=dss$aN, dN=dss$dN, aP=dss$aP, dP=dss$dP, aK=dss$aK, dK=dss$dK, rN=dss$rN, rP=dss$rP, rK=dss$rK, max_yield=dss$max_yield,  tolerance=dss$tolerance,
             WLY = dss$water_limited_yield)
  
    ## 2. Current yield: 
  
  
  actualUptake <- cbind(supply, actual_uptake_tool(supply))
  minmax_Yield <-  cbind(actualUptake, max_min_yields_tools(actualUptake))
  getcy <- final_yield_tools(minmax_Yield)## yield at zero input
  
  return(cbind(SoilData_onePix, data.frame(Current_Yield= getcy)))
  
}




#' get optimized fertilizer rate, target yield for the recommended rate and net revenue given cost and investment
run_Optim_NG2 <- function(rootUP, QID, fertilizer, invest, plDate, WLYData, lat, lon, areaHa, HD, WLY, DCY, country){	
  
  
  ## input of CY and WLY are in dry wt in KG/ha
  QID$water_limited_yield <- WLY	
  initial <- rep(0, nrow(fertilizer))
  lowerST <- rep(0, nrow(fertilizer))
  
  ## both CY and TY should be changed to user land size in ton/ha and fresh wt
  CY_user <- ((getRFY(HD = HD, RDY =  DCY , country = country))/1000)  * areaHa
  WLY_user <- ((getRFY(HD = HD, RDY =  WLY , country = country))/1000)  * areaHa
  
  
  FR <- optim(par=initial, fn=optim_NR, lower=lowerST, method = "L-BFGS-B", control=list(fnscale=-1), rootUP=rootUP, 
              QID=QID, CY=DCY, fertilizer=fertilizer, invest=invest, HD=HD, country = country)$par
  
  if(all(FR == 0)){
    return(data.frame(lat=lat, lon=lon, plDate, N=0, P=0, K= 0,WLY=WLY_user, CurrentY = CY_user, TargetY = CY_user,  TC =0, NR=0))
  }else{
    
    fertilizer$FR <- FR	
    
    ## NPK rate for ha of land	
    N <- as.vector(FR %*% fertilizer$N_cont)
    P <- as.vector(FR %*% fertilizer$P_cont)
    K <- as.vector(FR %*% fertilizer$K_cont)
    rec <- c(N, P, K)	
    
    ## NPK rate for user land size
    NPK_user <- rec * areaHa
    
    ## TY for ha of land 	
    TY <- QUEFTS1_Pedotransfer(QID, rec)	# Yield possible at recommended NPK in kg/ha dry wt. 
    
    ## both CY and TY should be changed to user land size in ton/ha and fresh wt
    TY_user <- ((getRFY(HD = HD, RDY = TY, country = country))/1000) * areaHa
    
    # CY_user <- CY * areaHa * 0.003
    # TY_user <- TY * areaHa * 0.003
    # WLY_user <- QID$water_limited_yield * areaHa * 0.003
    
    ## total cost per ha
    TC <- (sum(FR * fertilizer$price))* areaHa
    
    ## net revenue on the users land size
    GR <- (TY_user - CY_user) * rootUP                      # Gross revenue given root up is for fresh wt ton/ha
    NR <- round(GR - TC, digits=0)    											# Net Revenue
    
    ## reporting the recommended fertilizers
    Recomfr <- fertilizer[fertilizer$FR > 0, ]
    Recomfr$FR <- Recomfr$FR * areaHa
    Recomfr_wide <- spread(Recomfr[, c('type', 'FR')], type, FR)
    
    d1 <- data.frame(lat=lat, lon=lon, plDate, N=NPK_user[1], P=NPK_user[2], K= NPK_user[3],
                     WLY=WLY_user, CurrentY = CY_user, TargetY = TY_user,  TC =TC, NR=NR)
    d2 <- cbind(d1, Recomfr_wide)
    row.names(d2) <- NULL
    if(d2$NR <=0 | d2$TargetY <= d2$CurrentY){
      fertinfo <- subset(d2, select = c(lat, lon, plDate, N, P, K, WLY, CurrentY, TargetY, TC, NR))
      fertinfo$N <- fertinfo$P <- fertinfo$K <- fertinfo$TC <- fertinfo$NR <- 0
      fertinfo$TargetY <- fertinfo$CurrentY
      d2 <- fertinfo
    }
    
    return(d2)
  }
  
}




#'  Optimize the UREA, TSP and MOP needed to maximize the NR. x1, x2, x3 = Urea, MOP and Nafaka kg/ha. 
optim_NR <- function(fertRate, rootUP, QID, CY, fertilizer, invest, HD, country){ 
  f_price <-fertilizer$price	
  TC <- sum(fertRate * f_price)					
  
  ## Kg of Urea, Kg of NPK151515, Kg of NPK201010, Kg of MOP
  
  N <- as.vector(fertRate %*% fertilizer$N_cont)
  P <- as.vector(fertRate %*% fertilizer$P_cont)
  K <- as.vector(fertRate %*% fertilizer$K_cont)		
  
  rec <- c(N,P,K)	
  
  TotalYield <- QUEFTS1_Pedotransfer(QID, rec)
  
  AdditionalYield <- (getRFY(HD = HD, RDY = (TotalYield - CY), country = country))/1000 ## DM is converted to FW and then from KG/ha to ton/ha
  #AdditionalYield <- (TotalYield - CY)*0.003
  PriceYield <- AdditionalYield * rootUP
  NetRev <- PriceYield - TC
  if (!is.na(invest) & TC > invest) {NetRev <- NetRev - (invest-TC)^2} #penalize NR if costs exceed investment cap	
  return(NetRev)		
}





optim_NR_potato <- function(fertRate, rootUP, WLY_CY, CY, fertilizer, invest){ 
  f_price <-fertilizer$price	
  TC <- sum(fertRate * f_price)					
  
  ## Kg of Urea, Kg of NPK151515, Kg of NPK201010, Kg of MOP
  
  # N <- as.vector(fertRate %*% fertilizer$N_cont)
  # P <- as.vector(fertRate %*% fertilizer$P_cont)
  # K <- as.vector(fertRate %*% fertilizer$K_cont)	
  
  N <- sum(fertRate * fertilizer$N_cont)
  P <- sum(fertRate * fertilizer$P_cont)
  K <- sum(fertRate * fertilizer$K_cont)
  
  rec <- c(N,P,K)	
  
  # rec <- rec * c(0.41, 0.25, 0.65)
  
  TotalYield <- QUEFTS1_Pedotransfer_potato(WLY_CY=WLY_CY, rec=rec)
  AdditionalYield <- (TotalYield - CY)/0.21

  PriceYield <- AdditionalYield * rootUP
  NetRev <- PriceYield - TC
  if (!is.na(invest) & TC > invest) {NetRev <- NetRev - (invest-TC)^2} #penalize NR if costs exceed investment cap	
  return(NetRev)		
}


QUEFTS1_Pedotransfer_potato <- function(WLY_CY, rec){		
   WLY_CY$WLY <- WLY_CY$water_limited_yield
  
  # crop_param <- cbind(NUE(HI=0.75), data.frame(rN=0, rP=0, rK=0, max_yield=QID$WLY, tolerance=0.01))	## nutrient use efficiency of the crop
  crop_param <- data.frame(aN=41, dN=96, aP=233, dP=588,  aK=34, dK=161, rN=0,rP=0, rK=0, max_yield = WLY, tolerance = 0.01 )
 
  
  Queft_Input_Data_Var1 <- cbind(WLY_CY, crop_param)
 
  indata <- Queft_Input_Data_Var1[,c("WLY","aN", "dN", "aP", "dP","aK","dK", "rN", "rP", "rK", "soilN", "soilP", "soilK","max_yield", "tolerance")]
  
  N_rate <- rec[1]*0.41
  P_rate <- rec[2]*0.25
  K_rate <- rec[3]*0.65
  
  TargetYield_from_NPK <- NPK_TargetYield_potato(NutrUse_soilNPK=indata, N_rate, P_rate, K_rate)
  
  return(TargetYield_from_NPK)
}










NPK_TargetYield_potato <- function(NutrUse_soilNPK, N_rate, P_rate, K_rate){	
  NutrUse_soilNPK$N_rate <- N_rate
  NutrUse_soilNPK$P_rate <- P_rate
  NutrUse_soilNPK$K_rate <- K_rate
  
  ## Supply of nutrients to the crop
  NutrUse_soilNPK$SN <- NutrUse_soilNPK$N_rate + NutrUse_soilNPK$soilN
  NutrUse_soilNPK$SP <- NutrUse_soilNPK$P_rate + NutrUse_soilNPK$soilP
  NutrUse_soilNPK$SK <- NutrUse_soilNPK$K_rate + NutrUse_soilNPK$soilK
  
  ## Actual Uptake of nutrients: crop param + nutrient supply
  tmp <-  actual_uptake_tool(NutrUse_soilNPK)
  NutrUse_soilNPK <- cbind(NutrUse_soilNPK, tmp)
  
  ## max and min yield: actual uptake and crop param. min of N uptake constrianed by availability of P, K and water
  maxminY <- max_min_yields_tools(NutrUse_soilNPK)
  NutrUse_soilNPK <- cbind(NutrUse_soilNPK, maxminY)	
  
  ## final yield: min yield for combined uptake of 2 nutrients assuming the 3rd is not limiting, should be < WLY, and take meanof the six combinations 
  Target_Yield <- quefts_tools(NutrUse_soilNPK)	
  return(Target_Yield$FinalYield)
}


actual_uptake_tool <- function(ds_supply){	
  with(ds_supply,
       {
         UNP <- nutrient_uptake(S1 = SN, S2 = SP, d1 = dN, a1 = aN, d2 = dP, a2 = aP, r1 = rN, r2 = rP)
         UNK <- nutrient_uptake(S1 = SN, S2 = SK, d1 = dN, a1 = aN, d2 = dK, a2 = aK, r1 = rN, r2 = rK)
         UNW <- water_dependent_nutrient_uptake(S1 = SN, WLY = WLY, d1 = dN, a1 = aN, r1 = rN)
         UN <- min(UNP, UNK, UNW)
         
         
         UPN <- nutrient_uptake(S1 = SP, S2 = SN, d1 = dP, a1 = aP, d2 = dN, a2 = aN, r1 = rP, r2 = rN)
         UPK <- nutrient_uptake(S1 = SP, S2 = SK, d1 = dP, a1 = aP, d2 = dK, a2 = aK, r1 = rP, r2 = rK)
         UPW <- water_dependent_nutrient_uptake(S1 = SP, WLY = WLY, d1 = dP, a1 = aP, r1 = rP)
         UP <- min(UPN, UPK, UPW)
         
         
         UKN <- nutrient_uptake(S1 = SK, S2 = SN, d1 = dK, a1 = aK, d2 = dN, a2 = aN, r1 = rK, r2 = rN)
         UKP <- nutrient_uptake(S1 = SK, S2 = SP, d1 = dK, a1 = aK, d2 = dP, a2 = aP, r1 = rK, r2 = rP)
         UKW <- water_dependent_nutrient_uptake(S1 = SK, WLY = WLY, d1 = dK, a1 = aK, r1 = rK)
         UK <- min(UKN, UKP, UKW)
         
         
         return(data.frame(UN=UN, UP=UP, UK=UK))
       })
}



#' Nutrient uptake depends on the soil supply of the nutrient and the supply of other nutrients
nutrient_uptake <- function(S1=NA, S2=NA, d1=NA, a1=NA, d2=NA, a2=NA, r1=NA, r2=NA) {	
  # N, P and K uptakes based on QUEFTS
  if (S1 < r1 + ((S2 - r2) * a2 / d1)) {
    uptakeX_givenY = S1
  } else if (S1 > r1 + ((S2 - r2) * (2 * d2 / a1 - a2 / d1))) {
    uptakeX_givenY = r1 + (S2 - r2) * (d2 / a1)
  } else {
    uptakeX_givenY = S1 - 0.25 * (S1 - r1 - (S2 - r2) * (a2 / d1))^2 / ((S2 - r2) * (d2 / a1 - a2 / d1))
  }
  # Nutrient uptake given availability of other nutrient
  return(uptakeX_givenY)
}



water_dependent_nutrient_uptake <- function(S1=NA, WLY=NA, d1=NA, a1=NA, r1=NA) {	
  if (S1 < r1 + WLY / d1) {
    uptakeX_givenWater = S1    
  } else if (S1 > r1 + 2*WLY/a1 - WLY/d1) {
    uptakeX_givenWater = WLY / a1    
  } else {
    uptakeX_givenWater = S1 - 0.25 * (S1 - r1 - WLY/d1)^2 / (WLY / a1 - WLY / d1)    
  }
  
  return(uptakeX_givenWater)
}


max_min_yields_tools <- function(dss){	
  
  YNA <- max((dss$UN - dss$rN), 0) * dss$aN
  YND <- max((dss$UN - dss$rN), 0) * dss$dN
  YPA <- max((dss$UP - dss$rP), 0) * dss$aP
  YPD <- max((dss$UP - dss$rP), 0) * dss$dP
  YKA <- max((dss$UK - dss$rK), 0) * dss$aK
  YKD <- max((dss$UK - dss$rK), 0) * dss$dK
  
  
  return(data.frame(YNA=YNA, YND=YND, YPA=YPA, YPD=YPD, YKA=YKA, YKD=YKD))
  
}



quefts_tools <- function(supply_wly){	
  # Actual uptake of nutrients.
  tmp <- actual_uptake_tool(supply_wly)
  supply_wly$UN <- tmp[[1]]
  supply_wly$UP <- tmp[[2]]
  supply_wly$UK <- tmp[[3]]
  
  # Maximum and minimum yields, depending on maximum accumulation and dilution.
  yields <- max_min_yields_tools(supply_wly)
  supply_wly$YNA <- yields$YNA
  supply_wly$YND <- yields$YND
  supply_wly$YPA <- yields$YPA
  supply_wly$YPD <- yields$YPD
  supply_wly$YKA <- yields$YKA
  supply_wly$YKD <- yields$YKD
  
  # Final yield based on the combinations of nutrient uptake and minimum + maximum yields.
  supply_wly$FinalYield <- final_yield_tools(supply_wly)
  
  return(supply_wly)
}


#' after setting fertilizer recommendation <25 kg/ha Urea, MOP or Nafaka, target yield with the remaining recommended fertilizer is  re-estimated  and 
#'  total cost, gross and net revenue are re calcuated.
#' @param rootUP cassava root price 
#' @param zone 
#' @param wdd has dry wt
#' @param rdd has fresh wt
#' @param fertilizer 
#' @author Meklit
#' @export

Rerun_25kgKa_potato <- function(rootUP, rdd, fertilizer, WLY_CY, onlyFert){
  
  WLY <- WLY_CY$water_limited_yield
  DCY <- WLY_CY$Current_Yield
  
  fertilizer <- merge(fertilizer, onlyFert, by='type')
  
  TC <- sum(fertilizer$price * fertilizer$amount) 
  
  # N  <- as.vector(fertilizer$amount %*% fertilizer$N_cont)
  # P  <- as.vector(fertilizer$amount %*% fertilizer$P_cont)
  # K  <- as.vector(fertilizer$amount %*% fertilizer$K_cont)
  
  
  N  <- sum(fertilizer$amount * fertilizer$N_cont)
  P  <- sum(fertilizer$amount * fertilizer$P_cont)
  K  <- sum(fertilizer$amount * fertilizer$K_cont)
  
  
  
  rec <- c(N, P, K)
  TY  <- QUEFTS1_Pedotransfer_potato(WLY_CY=WLY_CY, rec=rec)					#dry wt yield in kg/ha
  
  TY_user  <- TY / 0.21
  CY_user  <- DCY / 0.21
  
  
  rdd$CurrentY <- CY_user
  rdd$TargetY <- TY_user
  rdd$TC <- TC
  rdd$NR <- ((rdd$TargetY - rdd$CurrentY)*rootUP) - rdd$TC
  rdd$N <- N
  rdd$P <- P
  rdd$K <- K
  
  if(rdd$TargetY <= rdd$CurrentY | rdd$NR <= 0 ){
    rdd$N <- rdd$P <- rdd$K <- rdd$TC <- rdd$NR <- 0
    rdd$TargetY <- rdd$CurrentY
  }
  
  return(rdd)
  #	
  #	return(data.frame(lat=recalc_data$lat, lon=recalc_data$lon, rateUrea,  rateNPK151515, rateNPK201010,rateMOP, currentY= CY, 
  #					targetY=TY, WLY=WLY, netRev=NR,totalCost = TC, N=N, P=P, K=K, plDate=recalc_data$plDate))
}


NRabove18Cost_potato <- function(ds) {
  if (ds$NR < ds$TC * 0.2) {
    fertRecom <- subset(ds, select = c(N, P, K, WLY, CurrentY, TargetY, TC, NR,soilN,soilP,soilK ))
    fertRecom$N <- fertRecom$P <- fertRecom$K <- fertRecom$TC <- fertRecom$NR <- 0
    fertRecom$TargetY <- fertRecom$CurrentY
    
    onlyFert <- subset(ds, select = -c(N, P, K, WLY, CurrentY, TargetY, TC, NR,soilN,soilP,soilK ))
    row.names(onlyFert) <- NULL
    onlyFert[, 1:ncol(onlyFert)] <- 0
    
    fertRecom <- cbind(fertRecom, onlyFert)
    ds <- fertRecom
  }
  row.names(ds) <- NULL
  return(ds)
}



#' Yield calculated based on the combined uptake of 2 nutrients, while taking into account the availability of the third nutrient.  
yield_nutrients_combined <- function(U1=NA, d1=NA, a1=NA, Y2A=NA, Y2D=NA, Y3D=NA, r1=NA){	
  # Determine which nutrient limited yield is lowest.
  YxD = min(Y2D, Y3D)	
  # If the uptake of one of the nutrients, and therefore the yield associated with that 
  # nutrient, is zero the overall yield is also zero.
  if (U1 == 0 || YxD == 0) {
    Y12 = 0
  }else{
    Y12 = Y2A + (2 * (YxD - Y2A) * (U1 - r1 - Y2A / d1)) / (YxD / a1 - Y2A / d1) -
      (YxD - Y2A) * (U1 - r1 - Y2A / d1)^2 / (YxD / a1 - Y2A / d1)^2
  }
  # Return the calculated yield based on the uptake of nutrients 1 and 2
  return(Y12)
}


final_yield_tools <- function(Uptake_Yield){	
  with(Uptake_Yield, 
       {
         YNP <- yield_nutrients_combined(U1 = UN, d1 = dN, a1 = aN, Y2A = YPA, Y2D = YPD, Y3D = YKD, r1 = rN)
         YNK <- yield_nutrients_combined(U1 = UN, d1 = dN, a1 = aN, Y2A = YKA, Y2D = YKD, Y3D = YPD, r1 = rN)
         YPN <- yield_nutrients_combined(U1 = UP, d1 = dP, a1 = aP, Y2A = YNA, Y2D = YND, Y3D = YKD, r1 = rP)
         YPK <- yield_nutrients_combined(U1 = UP, d1 = dP, a1 = aP, Y2A = YKA, Y2D = YKD, Y3D = YND, r1 = rP)
         YKN <- yield_nutrients_combined(U1 = UK, d1 = dK, a1 = aK, Y2A = YNA, Y2D = YND, Y3D = YPD, r1 = rK)
         YKP <- yield_nutrients_combined(U1 = UK, d1 = dK, a1 = aK, Y2A = YPA, Y2D = YPD, Y3D = YND, r1 = rK)
         
         # Make sure the nutrient limited yields do not exceed the maximum possible yield = WLY
         YNPc <- min(c(YNP, YND, YPD, YKD, WLY))
         YNKc <- min(c(YNK, YND, YPD, YKD, WLY))
         YPNc <- min(c(YPN, YND, YPD, YKD, WLY))
         YPKc <- min(c(YPK, YND, YPD, YKD, WLY))
         YKNc <- min(c(YKN, YND, YPD, YKD, WLY))
         YKPc <- min(c(YKP, YND, YPD, YKD, WLY))
         
         #Final estimate
         YEc <- mean(c(YNPc, YNKc, YPNc, YPKc, YKNc, YKPc))
         
         if(YEc > WLY){
           YEc == WLY
         }
         
         return(YEc)         
       })
}


#TODO: price of fertilizers for tanzania is different so from GPS we need to define the zone and take the correct price accordingly. default is at the end of the script
fertilizerFunc <- function(ureaavailable=TRUE, ureaCostperBag=NA,ureaBagWt=50,
                           MOPavailable=TRUE, MOPCostperBag=NA, MOPBagWt=50,
                           DAPavailable=TRUE, DAPCostperBag=NA, DAPBagWt=50, 
                           NPK201010available=TRUE, NPK201010CostperBag=NA, NPK201010BagWt=50, 
                           NPK151515available=TRUE, NPK151515CostperBag=NA, NPK151515BagWt=50,
                           TSPavailable=TRUE, TSPCostperBag=NA, TSPBagWt=50, 
                           NPK171717available=FALSE, NPK171717CostperBag=NA, NPK171717BagWt=50, 
                           Nafakaavailable=TRUE, NafakaCostperBag=NA, NafakaBagWt=50,
                           CANavailable=TRUE, CANCostperBag=NA, CANBagWt=50, 
                           SSPavailable=TRUE, SSPCostperBag=NA, SSPBagWt=50,
                           YaraMila_UNIKavailable=TRUE, YaraMila_UNIKCostperBag=NA, YaraMila_UNIKBagWt=50,country){
  if(country == "NG"){
    if(ureaCostperBag == 0) {ureaCostperBag <- 7500} else {ureaCostperBag <- as.numeric(ureaCostperBag)}
    if(MOPCostperBag == 0) {MOPCostperBag <- 13500}else {MOPCostperBag <- as.numeric(MOPCostperBag)}
    if(DAPCostperBag == 0) {DAPCostperBag <- 13250}else {DAPCostperBag <- as.numeric(DAPCostperBag)}
    if(NPK201010CostperBag == 0) {NPK201010CostperBag <- 7200}else {NPK201010CostperBag <- as.numeric(NPK201010CostperBag)}
    if(NPK151515CostperBag == 0) {NPK151515CostperBag <- 8500}else {NPK151515CostperBag <- as.numeric(NPK151515CostperBag)}
    if(TSPCostperBag == 0) {TSPCostperBag <- 13250}else {TSPCostperBag <- as.numeric(TSPCostperBag)}
    if(NPK171717CostperBag == 0) {NPK171717CostperBag <- 7800}else {NPK171717CostperBag <- as.numeric(NPK171717CostperBag)}
    if(CANCostperBag == 0) {CANCostperBag <- 7500}else {CANCostperBag <- as.numeric(CANCostperBag)}
    if(SSPCostperBag == 0) {SSPCostperBag <- 22364}else {SSPCostperBag <- as.numeric(SSPCostperBag)}
  }else if (country == "TZ"){
    if(ureaCostperBag == 0) {ureaCostperBag <- 58000}else {ureaCostperBag <- as.numeric(ureaCostperBag)}##65000
    if(MOPCostperBag == 0) {MOPCostperBag <- 67000}else {MOPCostperBag <- as.numeric(MOPCostperBag)}##120000
    if(DAPCostperBag == 0)  {DAPCostperBag <- 60000}else {DAPCostperBag <- as.numeric(DAPCostperBag)}##85000
    #if(levels(as.factor(NafakaCostperBag)) == "NA") {NafakaCostperBag = NA}
    if(NPK201010CostperBag == 0) {NPK201010CostperBag <- 68000}else {NPK201010CostperBag <- as.numeric(NPK201010CostperBag)}
    # if(NPK151515CostperBag == 0) {NPK151515CostperBag <- 60000}else {NPK151515CostperBag <- as.numeric(NPK151515CostperBag)}
    if(TSPCostperBag == 0) {TSPCostperBag <- 64000}else {TSPCostperBag <- as.numeric(TSPCostperBag)}#80000
    if(NPK171717CostperBag == 0) {NPK171717CostperBag <- 63000}else {NPK171717CostperBag <- as.numeric(NPK171717CostperBag)}##61000						
    if(CANCostperBag == 0) {CANCostperBag <- 53000}else {CANCostperBag <- as.numeric(CANCostperBag)}#65000
    # if(SSPCostperBag == 0) {SSPCostperBag <- 80000}else {SSPCostperBag <- as.numeric(SSPCostperBag)}		
    if(YaraMila_UNIKCostperBag == 0) {YaraMila_UNIKCostperBag <- 60000}	else {YaraMila_UNIKCostperBag <- as.numeric(YaraMila_UNIKCostperBag)}##61000
  }	else if( (country == "RW")){
    if(ureaCostperBag == 0) {ureaCostperBag <- 22350}else {ureaCostperBag <- as.numeric(ureaCostperBag)}##65000
    if(MOPCostperBag == 0) {MOPCostperBag <- 24050}else {MOPCostperBag <- as.numeric(MOPCostperBag)}##120000
    if(DAPCostperBag == 0)  {DAPCostperBag <- 24000}else {DAPCostperBag <- as.numeric(DAPCostperBag)}##85000
    if(TSPCostperBag == 0) {TSPCostperBag <- 26150}else {TSPCostperBag <- as.numeric(TSPCostperBag)}#80000
    if(NPK171717CostperBag == 0) {NPK171717CostperBag <- 30150}else {NPK171717CostperBag <- as.numeric(NPK171717CostperBag)}##61000						
    
  }					
  
  ureaFert <- data.frame(type = 'Urea', available =ureaavailable,  N_cont = 0.46, P_cont = 0, K_cont=0, costPerBag = ureaCostperBag, bagWeight=ureaBagWt ) 
  MOPFert <- data.frame(type = 'MOP', available = MOPavailable, N_cont = 0.00, P_cont = 0.00, K_cont=0.60, costPerBag = MOPCostperBag, bagWeight=MOPBagWt)
  DAPFert <- data.frame(type = 'DAP',available = DAPavailable,  N_cont = 0.18, P_cont = 0.20, K_cont=0.0, costPerBag = DAPCostperBag, bagWeight=DAPBagWt)
  NPK201010Fert <- data.frame(type = 'NPK20_10_10',available = NPK201010available,  N_cont = 0.20, P_cont = 0.044, K_cont=0.083, costPerBag = NPK201010CostperBag, bagWeight=NPK201010BagWt)
  NPK151515Fert <- data.frame(type = 'NPK15_15_15',available = NPK151515available,  N_cont = 0.15, P_cont = 0.07, K_cont=0.125, costPerBag = NPK151515CostperBag, bagWeight=NPK151515BagWt)
  TSPFert <- data.frame(type = 'TSP',available = TSPavailable,  N_cont = 0.0, P_cont = 0.2, K_cont=0.0, costPerBag = TSPCostperBag, bagWeight=TSPBagWt)
  NPK171717Fert <- data.frame(type = 'NPK17_17_17',available = NPK171717available,  N_cont = 0.17, P_cont = 0.083, K_cont=0.15, costPerBag = NPK171717CostperBag, bagWeight=NPK171717BagWt)# TODO get price
  # Minjingu_Nafaka <- data.frame(type = 'Minjingu_Nafaka+',available = Nafakaavailable,  N_cont = 0.09, P_cont = 0.07, K_cont=0.06, costPerBag = NafakaCostperBag, bagWeight=NafakaBagWt)
  CANFert <- data.frame(type = 'CAN',available = CANavailable,  N_cont = 0.01, P_cont = 0.01, K_cont=0.01,costPerBag = CANCostperBag, bagWeight=CANBagWt)## not correct value TODO check
  SSPFert <- data.frame(type = 'SSP', available = SSPavailable, N_cont = 0.01, P_cont = 0.01, K_cont=0.01, costPerBag = SSPCostperBag, bagWeight=SSPBagWt)## not correct value TODO check
  YaraMila_UNIKFert <- data.frame(type = 'YaraMila_UNIK', available = YaraMila_UNIKavailable, N_cont = 0.17, P_cont = 0.17, K_cont=0.17, costPerBag = YaraMila_UNIKCostperBag, bagWeight=YaraMila_UNIKBagWt)## 
  
  
  if(country == "NG"){
    fd_cont <- rbind(ureaFert,MOPFert , DAPFert,CANFert,NPK171717Fert,NPK151515Fert, NPK201010Fert,TSPFert,SSPFert)
  }else if(country == "TZ"){
    fd_cont <- rbind(ureaFert,MOPFert , DAPFert, CANFert, NPK171717Fert, NPK151515Fert, NPK201010Fert,TSPFert,SSPFert,YaraMila_UNIKFert)
  }else  if(country == "RW"){
    fd_cont <- rbind(ureaFert,MOPFert , DAPFert, NPK171717Fert, TSPFert)
  }
  
  
  fd_cont <- droplevels(fd_cont[fd_cont$available == TRUE, ])
  fd_cont$costPerBag <- as.numeric(fd_cont$costPerBag)
  fd_cont$price <- fd_cont$costPerBag / fd_cont$bagWeight	
  fd_cont <- subset(fd_cont, select=-c(available))
  
 
  return(fd_cont)
}






getFRrecommendations_potato <- function(maxInv, fertilizers, rootUP, areaHa,WLY_CY){
  
  
  if(WLY_CY$water_limited_yield == WLY_CY$Current_Yield) {
     fertinfo <- subset(WLY_CY, select = c(lat, lon, CONclass,NAME_1, NAME_2, WLY, CY))
     fertinfo$TargetY <- fertinfo$CY
     fertinfo$N <- 0
     fertinfo$P <- 0
     fertinfo$K <- 0
     fertinfo$TC <- 0
     fertinfo$NR <- 0
     fertinfo$Urea <- 0
     fertinfo$NPK17_17_17 <- 0
     fertinfo$DAP <- 0
     fertinfo <- fertinfo[, c("lat","lon","N","P","K","CONclass","NAME_1","NAME_2","WLY",
                              "CY","TargetY","TC","NR","Urea","NPK17_17_17", "DAP")]
     return(list(rec=fertinfo, fertilizer_rates=NULL))
             
  }else{
    #############################
    ## 1. get WLY, CY, fert recom and soil data
    WLY <- WLY_CY$water_limited_yield ## DM in kg/ha
    DCY <- WLY_CY$Current_Yield## DM in kg/ha
    
    
    ## 2. change investment from given areaHa to 1ha
    invest <- (maxInv / areaHa)	
    
    
    
    ## 3. optimize the fertilizer recommendation for maxInv in local currency and provide expected target yield in kg
    fert_optim <- run_Optim_Potato(rootUP=rootUP, fertilizer=fertilizers, invest=invest, WLY_CY=WLY_CY, areaHa=areaHa)
    fert_optim <- subset(fert_optim, select=-c(location,soilN,soilP,soilK, CON, 
                                               water_limited_yield, Current_Yield))
    
    if(fert_optim$NR == 0){ ## no fertilizer recommendation
      fertilizer_rates <- NULL
      return(list(rec=fert_optim, fertilizer_rates=fertilizer_rates))
    }else{
      fertinfo <- subset(fert_optim, select = c(lat, lon, N, P, K, CONclass,NAME_1, NAME_2, WLY, CY, TargetY, TC, NR))
      onlyFert <- subset(fert_optim, select = -c(lat, lon, N, P, K, CONclass,NAME_1, NAME_2, WLY, CY, TargetY, TC, NR))
      
      ## 4. remove ferilizer application < 25 kg/ha and re run the TY and NR calculation
      RecomperHa <- onlyFert/areaHa
      RecomperHa2 <- gather(RecomperHa, type, rate)
      onlyFert2 <- droplevels(RecomperHa2[RecomperHa2$rate > 25, ])
      
      if(nrow(onlyFert2) == 0 ){ ## if all fertilizer recom < 25 kg/ha all will be set to 0
        fertinfo$N <- fertinfo$P <- fertinfo$K <- fertinfo$NR <- fertinfo$TC <- 0
        fertinfo$TargetY <- fertinfo$CY
        fertilizer_rates <- NULL
        return(list(rec=fertinfo, fertilizer_rates=fertilizer_rates))
      }else if (ncol(onlyFert) == nrow(onlyFert2)){ ## if all fertilizer recom are >= 25 kg/ha they will be kept and only checked for NR >= 18% of invest
        Reset_fert_Cont <- fert_optim
        GPS_fertRecom <- NRabove18Cost(ds=Reset_fert_Cont)
        rec <- subset(GPS_fertRecom, select = c(lat, lon, N, P, K, CONclass,NAME_1, NAME_2, WLY, CY, TargetY, TC, NR))
        frates <- subset(GPS_fertRecom, select = -c(lat, lon, N, P, K, CONclass,NAME_1, NAME_2, WLY, CY, TargetY, TC, NR))
        frates2 <- gather(frates, type, rate)
        return(list(rec=rec, fertilizer_rates=frates2))
        
      }else{
        fert25 <- spread(onlyFert2, type, rate) ## when some fertilizer recom are dropped b/c < 25 kg/ha, ty and NR should be recalculated
        fert_optim2 <- cbind(fertinfo, fert25)	
        fertilizer <- fertilizers[fertilizers$type %in% onlyFert2$type, ]
        Reset_fert_Cont <- Rerun_25kgKa_try(rootUP=rootUP, rdd=fert_optim2, fertilizer=fertilizer, WLY_CY=WLY_CY, onlyFert=onlyFert2, areaHa=areaHa)
        if(Reset_fert_Cont$NR <= 0){ ## after rerunning after avoiding <25KG/ha fertilizers, if NR <=0
          fertilizer_rates <- NULL
          return(list(rec=Reset_fert_Cont, fertilizer_rates=fertilizer_rates))
        }else{
          GPS_fertRecom <- NRabove18Cost(ds=Reset_fert_Cont)
          rec <- subset(GPS_fertRecom, select = c(lat, lon, N, P, K, CONclass,NAME_1, NAME_2, WLY, CY, TargetY, TC, NR))
          frates <- subset(GPS_fertRecom, select = -c(lat, lon, N, P, K, CONclass,NAME_1, NAME_2, WLY, CY, TargetY, TC, NR))
          frates2 <- gather(frates, type, rate)
          return(list(rec=rec, fertilizer_rates=frates2))
          
        }
      }
    }
  }
}



Rerun_25kgKa_try <- function(rootUP, rdd, fertilizer, WLY_CY, onlyFert, areaHa=areaHa){
  
  fertilizer <- merge(fertilizer, onlyFert, by='type')
  TC <- (sum(fertilizer$price %*% fertilizer$rate) ) * areaHa
  N  <- as.vector(fertilizer$rate %*% fertilizer$N_cont)
  P  <- as.vector(fertilizer$rate %*% fertilizer$P_cont)
  K  <- as.vector(fertilizer$rate %*% fertilizer$K_cont)
  rec <- c(N, P, K)
  
  ## NPK rate for user land size
  NPK_user <- rec * areaHa
  
  TY <- QUEFTS1_Pedotransfer_potato(WLY_CY = WLY_CY, rec = rec)
  #TY_user  <- ((getRFY(HD = as.Date(HD), RDY = TY, country = country))/1000) * areaHa
  TY_user  <- TY/0.21 * areaHa
  CY_user  <- rdd$CY
  
  
  rdd$CurrentY <- CY_user
  rdd$TargetY <- TY_user
  rdd$TC <- TC
  rdd$NR <- ((rdd$TargetY - rdd$CurrentY)*rootUP) - rdd$TC
  rdd$N <- NPK_user[1]
  rdd$P <- NPK_user[2]
  rdd$K <- NPK_user[3]
  
  if(rdd$TargetY <= rdd$CurrentY){
    rdd$N <- rdd$P <- rdd$K <- rdd$TC <- rdd$NR <- 0
    rdd$TargetY <- CY_user
  }
  
  if(rdd$NR <=0 | rdd$TargetY <= rdd$CurrentY){
    fertinfo <- subset(rdd, select = c(lat, lon, N, P, K, CONclass,NAME_1, NAME_2, WLY, CY, TargetY, TC, NR))
    fertinfo$N <- fertinfo$P <- fertinfo$K <- fertinfo$TC <- fertinfo$NR <- 0
    fertinfo$TargetY <- fertinfo$CY
    rdd <- fertinfo
  }
  
  return(rdd)
}



#' get optimized fertilizer rate, target yield for the recommended rate and net revenue given cost and investment
run_Optim_Potato <- function(rootUP=rootUP, fertilizer=fertilizers, invest=invest, WLY_CY=WLY_CY, areaHa=areaHa){	
  
  
  ## input of CY and WLY are in dry wt in KG/ha
  initial <- rep(0, nrow(fertilizer))
  lowerST <- rep(0, nrow(fertilizer))
  
  ## both CY and TY should be changed to user land size in ton/ha and fresh wt
  CY_fresh <- WLY_CY$CY
  WLY_fresh <- WLY_CY$WLY
  
  
  FR <- optim(par=initial, fn=optim_NR_potato, lower=lowerST, method = "L-BFGS-B", control=list(fnscale=-1),
              rootUP=rootUP, WLY_CY=WLY_CY, fertilizer=fertilizer, invest=invest)$par
  
  # if(all(FR == 0)){}
  
  fertilizer$FR <- FR	
  
  # fertilizer <- fertilizer[fertilizer$FR >=25, ] ## at least 25 kg/ha fertilizer
  # FR <- FR[FR>=25]
  
  ## NPK rate for ha of land	
  N <- as.vector(FR %*% fertilizer$N_cont)
  P <- as.vector(FR %*% fertilizer$P_cont)
  K <- as.vector(FR %*% fertilizer$K_cont)
  rec <- c(N, P, K)	
  
  ## NPK rate for user land size
  NPK_user <- rec * areaHa
  
  ## TY for ha of land 	
  TY <- QUEFTS1_Pedotransfer_potato(WLY_CY=WLY_CY, rec=rec)	# Yield possible at recommended NPK in kg/ha dry wt. 
  
  
  ## both CY and TY should be changed to user land size in ton/ha and fresh wt
  TY_user <- (TY/0.21) * areaHa
  
  if(TY_user > WLY_fresh){
    TY_user <- WLY_fresh
  }
  
  ## total cost per ha
  TC <- (sum(FR * fertilizer$price))* areaHa
  
  ## net revenue on the users land size
  GR <- (TY_user - CY_fresh) * rootUP                      # Gross revenue given root up is for fresh wt ton/ha
  NR <- round(GR - TC, digits=0)    											# Net Revenue
  
  ## reporting the recommended fertilizers
  Recomfr <- fertilizer[fertilizer$FR > 0, ]
  Recomfr$FR <- Recomfr$FR * areaHa
  Recomfr_wide <- spread(Recomfr[, c('type', 'FR')], type, FR)
  
  d1 <- data.frame( N=NPK_user[1], P=NPK_user[2], K= NPK_user[3], TargetY = TY_user,  TC =TC, NR=NR)
  
  if(all(FR == 0)){
    Recomfr_wide <- data.frame(urea = 0, NPK171717 = 0, DAP = 0)
    d1$NR <- 0
    d1$TC <- 0
    d1$TargetY <- WLY_CY$Current_Yield
  }
  
  d2 <- cbind(WLY_CY, d1, Recomfr_wide)
  row.names(d2) <- NULL
  if(d2$NR <=0 | d2$TargetY <= d2$CY){
    fertinfo <- subset(d2, select = c(N, P, K, TargetY, TC, NR))
    fertinfo$N <- fertinfo$P <- fertinfo$K <- fertinfo$TC <- fertinfo$NR <- 0
    fertinfo$TargetY <- CY_fresh
    Recomfr_wide[, 1:ncol(Recomfr_wide)]  <- 0
    d2 <- cbind(WLY_CY, fertinfo, Recomfr_wide)
  }
  return(d2)
  
}


#'  Optimize the UREA, TSP and MOP needed to maximize the NR. x1, x2, x3 = Urea, MOP and Nafaka kg/ha. 
optim_NR_potato <- function(fertRate, rootUP, WLY_CY, fertilizer, invest){ 
  f_price <-fertilizer$price	
  TC <- sum(fertRate * f_price)					
  
  ## Kg of Urea, Kg of NPK151515, Kg of NPK201010, Kg of MOP
  
  N <- as.vector(fertRate %*% fertilizer$N_cont)
  P <- as.vector(fertRate %*% fertilizer$P_cont)
  K <- as.vector(fertRate %*% fertilizer$K_cont)		
  
  rec <- c(N,P,K)	
  
  TotalYield <- QUEFTS1_Pedotransfer_potato(WLY_CY=WLY_CY, rec=rec)
  
  AdditionalYield <- (TotalYield - WLY_CY$Current_Yield)/0.21 ## DM is converted to FW and then from KG/ha to ton/ha
  #AdditionalYield <- (TotalYield - CY)*0.003
  PriceYield <- AdditionalYield * rootUP
  NetRev <- PriceYield - TC
  if (!is.na(invest) & TC > invest) {NetRev <- NetRev - (invest-TC)^2} #penalize NR if costs exceed investment cap	
  return(NetRev)		
}



# QUEFTS1_Pedotransfer_potato <- function(WLY_CY, rec){	
#   # crop_param <- data.frame(aN=41, dN=96, aP=233, dP=588,  aK=34, dK=161, rN=0,rP=0, rK=0,
#   #                         max_yield = WLY_CY$water_limited_yield, tolerance = 0.01 )
#   crop_param <- data.frame(aN=175, dN=435, aP=625, dP=4348,  aK=109, dK=513, rN=0.41,rP=0.25, rK=0.65,
#                            max_yield = WLY, tolerance = 0.01 )
#   # 
#   dss <- cbind(WLY_CY, crop_param)
#   
#   supply <- data.frame(rel_N=dss$rel_N, rel_P=dss$rel_P, rel_K=dss$rel_K, SN=dss$soilN+rec[1], SP=dss$soilP+rec[2], SK=dss$soilK+rec[3], water_limited_yield = dss$water_limited_yield,
#                        aN=dss$aN, dN=dss$dN, aP=dss$aP, dP=dss$dP, aK=dss$aK, dK=dss$dK, rN=dss$rN, rP=dss$rP, rK=dss$rK, max_yield=dss$max_yield,  tolerance=dss$tolerance,
#                        WLY = dss$water_limited_yield)
#   
#   actualUptake <- cbind(supply, actual_uptake_tool(supply))
#   minmax_Yield <-  cbind(actualUptake, max_min_yields_tools(actualUptake))
#   getcy <- final_yield_tools(minmax_Yield)
#   
#   return(getcy)
# }

# QUEFTS1_Pedotransfer_potato <- function(WLY_CY, rec){	
#   crop_param <- data.frame(aN=41, dN=96, aP=233, dP=588,  aK=34, dK=161, rN=0,rP=0, rK=0,
#                            max_yield = WLY_CY$water_limited_yield, tolerance = 0.01 )
#   indata <- cbind(WLY_CY, crop_param)
#   indata <- indata[,c("lat","lon" ,"WLY","aN", "dN", "aP", "dP","aK","dK", "rN", "rP", "rK", "soilN", "soilP", "soilK","max_yield", "tolerance")]
#   
#   N_rate <- rec[1]
#   P_rate <- rec[2]
#   K_rate <- rec[3]
#   
#   TargetYield_from_NPK <- NPK_TargetYield_potato(NutrUse_soilNPK=indata, N_rate, P_rate, K_rate)
#   
#   return(TargetYield_from_NPK$TargetYield)
# }





### see if profit is > (0.18 * total cost) + total cost
NRabove18Cost <- function(ds){
  
  if(ds$NR < ds$TC * 0.18){
    
    fertRecom <- subset(ds, select = c(lat, lon, N, P, K, CONclass,NAME_1, NAME_2, WLY, CY, TargetY, TC, NR))
    fertRecom$N <- fertRecom$P <- fertRecom$K <- fertRecom$TC <- fertRecom$NR <- 0
    fertRecom$TargetY <- fertRecom$CY
    
    onlyFert <- subset(ds, select = -c(lat, lon, N, P, K, CONclass,NAME_1, NAME_2, WLY, CY, TargetY, TC, NR))
    row.names(onlyFert) <- NULL
    for(j in 1:ncol(onlyFert)){
      onlyFert[,j] <-0
    }
    fertRecom <- cbind(fertRecom, onlyFert)
    ds <- fertRecom
  }	
  row.names(ds) <- NULL
  return(ds)
}