## ----setup, include=FALSE------------------------------------------------


## ----package, message=FALSE----------------------------------------------
if (!require(Rquefts)) { install.packages("Rquefts") }
library(Rquefts)  


## ----supply--------------------------------------------------------------
library(Rquefts)
# range of pH values
ph <- seq(4, 8, 1)
s <- nutSupply(ph, SOC=1.5, Kex=1, Polsen=1)
cbind(ph, s)


## ----phsupply------------------------------------------------------------
ph <- seq(4, 8, 0.25)
# three combinations
s1 <- nutSupply(ph, SOC=15, Kex=5, Polsen=5)
s2 <- nutSupply(ph, SOC=30, Kex=10, Polsen=10)
s3 <- nutSupply(ph, SOC=45, Kex=20, Polsen=20)


## ----phsupply2, fig.width=10---------------------------------------------
par(mfrow=c(1,3))

for (supply in c("N_base_supply", "P_base_supply", "K_base_supply")) {
  plot(ph, s3[,supply], xlab="pH", ylab=paste(supply, "(kg/ha)"), type="l", ylim=c(0,max(s3[,supply])))
  lines(ph, s2[,supply])
  lines(ph, s1[,supply])
}


## ----omsupply, width=10--------------------------------------------------
# unit is g/kg
om <- seq(0, 70, 5)
# three combinations
s1 <- nutSupply(pH=5.5, SOC=om, Kex=5, Polsen=5)
s2 <- nutSupply(pH=5.5, SOC=om, Kex=10, Polsen=10)
s3 <- nutSupply(pH=5.5, SOC=om, Kex=20, Polsen=20)

# 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])))
  lines(om, s2[,supply])
  lines(om, s1[,supply])
}


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


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


## ------------------------------------------------------------------------
soiltype <- quefts_soil()
barley <- quefts_crop('Barley')
fertilizer <- list(N=50, P=0, K=0)
biomass <- list(leaf_att=2200, stem_att=2700, store_att=4800, SeasonLength=110)
q <- quefts(soiltype, barley, fertilizer, biomass)


## ------------------------------------------------------------------------

library(Rquefts)

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)

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"]

par(mfrow=c(1,2))
plot(N, yield, type="l", lwd=2)
Efficiency = yield / N
plot(N, Efficiency, type="l", lwd=2)


## ------------------------------------------------------------------------

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)

# display results
library(raster  )
r0 <- rasterFromXYZ(x[x$K==0, -3])
r200 <- rasterFromXYZ(x[x$K==200, -3])

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)
#plot(stack(r0, r200))


## ------------------------------------------------------------------------
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)


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


## ------------------------------------------------------------------------
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)


## ------------------------------------------------------------------------

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)


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