## ----quefts_supply_1----------------------------------------------------------
library(Rquefts)
s <- nutSupply2(temp=25, pH=6, SOC=15, Kex=1, Polsen=10, Ptotal=500)
s


## ----quefts_supply_10, fig.with=8---------------------------------------------
ph <- seq(from=4.5, to=8.5, by=0.1)
s <- nutSupply2(temp=25, pH=ph, SOC=15, Kex=1, Polsen=10, Ptotal=500)
phs <- cbind(ph, s)
head(phs)
par(mfrow=c(1,3))
plot(phs[, c(1,2)])
plot(phs[, c(1,3)])
plot(phs[, c(1,4)])


## ----quefts_phsupply----------------------------------------------------------
s1 <- nutSupply2(25, ph, SOC=15, Kex=5,  Polsen=5, Ptotal=5*55)
s2 <- nutSupply2(25, ph, SOC=30, Kex=10, Polsen=10, Ptotal=10*55)
s3 <- nutSupply2(25, ph, SOC=45, Kex=20, Polsen=20, Ptotal=20*55)


## ----quefts_phsupply2, fig.width=10-------------------------------------------
par(mfrow=c(1,3))
sup <- c("N_base_supply", "P_base_supply", "K_base_supply")
for (i in 1:3) {
  yb <- paste(sup[i],"(kg/ha)")
  ym <- c(0, max(s3[, sup[i]]))
  plot(ph, s1[,sup[i]], xlab="pH", ylab=yb, type="l", ylim=ym, lty=1, las=2)
  lines(ph, s2[, sup[i]], lty=2)
  lines(ph, s3[, sup[i]], lty=3)
}
legend("bottomleft", paste("treatment", 1:3), lty=1:3)


## ----quefts_omsupply, fig.width=10--------------------------------------------
# unit is g/kg
om <- seq(0, 70, 5)
# three combinations
s1 <- nutSupply2(temp=25, pH=5.5, SOC=om, Kex=5, Polsen=5, Ptotal=5*55)
s2 <- nutSupply2(temp=25, pH=5.5, SOC=om, Kex=10, Polsen=10, Ptotal=10*55)
s3 <- nutSupply2(temp=25, pH=5.5, SOC=om, Kex=20, Polsen=20, Ptotal=20*55)

# and plot
par(mfrow=c(1,3))
for (supply in c("N_base_supply", "P_base_supply", "K_base_supply")) {
  plot(om, s3[,supply], xlab="Soil Carbon (g/kg)", ylab=paste(supply, "(kg/ha)"), type="l", ylim=c(0,max(s3[,supply])), lty=1)
  lines(om, s2[,supply], lty=2)
  lines(om, s1[,supply], lty=3)
}


## ----quefts_100---------------------------------------------------------------
library(Rquefts)
q <- quefts()


## ----quefts_110---------------------------------------------------------------
run(q)


## ----quefts_120soil-----------------------------------------------------------
soil <- quefts_soil()
class(soil)
str(soil)
soil$UptakeAdjust


## ----quefts_120crops, error=TRUE----------------------------------------------
quefts_crop(name="x")


## ----quefts_120bar, error=TRUE------------------------------------------------
crop <- quefts_crop(name="Barley")
str(crop)


## ----quefts_120fert-----------------------------------------------------------
fertilizer <- list(N=50, P=0, K=0)


## ----quefts_120fertbiom-------------------------------------------------------
biomass <- list(leaf_att=2200, stem_att=2700, store_att=4800, SeasonLength=110)


## ----quefts120mod-------------------------------------------------------------
q <- quefts(soil, crop, fertilizer, biomass)


## ----quefts120120-------------------------------------------------------------
str(q)


## ----quefts_120dollar---------------------------------------------------------
q$N
q$N <- 100
q$N


## ----quefts_120run------------------------------------------------------------
output <- run(q)
round(output, 1)


## ----quefts_130---------------------------------------------------------------
rice <- quefts_crop("Rice")
q <- quefts(crop=rice)
q$leaf_att <- 2651
q$stem_att <- 5053
q$store_att <- 8208
fert(q) <- list(N=0, P=0, K=0)
N <- seq(1, 200, 10)


## ----quefts_140---------------------------------------------------------------
results <- matrix(nrow=length(N), ncol=12)
colnames(results) <- names(run(q))
for (i in 1:length(N)) {
	q["N"] <- N[i]
	results[i,] <- run(q)
}
yield <- results[,"store_lim"]
yield


## ----quefts_150---------------------------------------------------------------
plot(N, yield, type="l", lwd=2)


## ----quefts_160---------------------------------------------------------------
f <- expand.grid(N=seq(0,400,10), P=seq(0,400,10), K=c(0,200))
x <- rep(NA, nrow(f))
for (i in 1:nrow(f)) {0
	q$N <- f$N[i] 
	q$P <- f$P[i] 
	q$K <- f$K[i] 
	x[i] <- run(q)["store_lim"]
}
x <- cbind(f, x)
head(x)


## ----quefts_plotraster, fig.width=10------------------------------------------
library(terra)
r0 <- rast(x[x$K==0, -3], type="xyz")
r200 <- rast(x[x$K==200, -3], type="xyz")

par(mfrow=c(1,2))
plot(r0, xlab="N fertilizer", ylab="P fertilizer", las=1, main="K=0")
contour(r0, add=TRUE)
plot(r200, xlab="N fertilizer", ylab="P fertilizer", las=1, main="K=200")
contour(r200, add=TRUE)


## ----quefts_180---------------------------------------------------------------
set.seed(99)
yldfun <- function(N, noise=TRUE) { 1000 +  500* log(N+1)/3 + noise * runif(length(N), -500, 500) }
N <- seq(0,300,25)
Y <- replicate(10, yldfun(N))
obs <- cbind(N, Y)


## ----quefts_190---------------------------------------------------------------
RMSE <- function(obs, pred) sqrt(mean((obs - pred)^2))


## ----quefts_200---------------------------------------------------------------
q <- quefts()
q$P <- 0
q$K <- 0

results <- matrix(nrow=length(N), ncol=12)
colnames(results) <- names(run(q))
for (i in 1:length(N)) {
	q$N <- N[i]
	results[i,] <- run(q)
}

yield <- results[,'store_lim']
RMSE(obs[,-1], yield)


## ----quefts_210---------------------------------------------------------------
vars <- c('soil$N_base_supply', 'soil$N_recovery', 'crop$NminStore', 'crop$NmaxStore')
f <- function(p) {
	if (any (p < 0)) return(Inf)
	if (p['crop$NminStore'] >= p['crop$NmaxStore']) return(Inf)
	if (p['soil$N_recovery'] > .6 | p['soil$N_recovery'] < .3) return(Inf)
	q[names(p)] <- p
	pred <- run(q)
	RMSE(obs[,-1], pred['store_lim'])
}

# try the function with some initial values
x <- c(50, 0.5,  0.011, 0.035)
names(x) <- vars
f(x)


## ----quefts_220---------------------------------------------------------------
par <- optim(x, f)
# RMSE
par$value
# optimal parameters
par$par