# NOT RUN {
data("shpdata1")
retdata <- shpdata1$TS1$wrc
condata <- shpdata1$TS1$hcc
condata <- condata[!is.na(condata[,1]),]
weightmethod <- "range"
ivap <- NULL
set.itermax <- 1
LikModel <- "rss" # ALTERNATIVE OPTION: LikModel = "-2logLik"
ALG <- "DE" # ALTERNATIVE OPTION: ALG = "modMCMC"
parL<-list("p"=c("thr"=0.05,"ths"=0.45,"alf1"=0.01,"n"=2,"Ks"=100,"tau"=.5),
"psel" = c(1, 1, 1, 1, 1, 1),
"plo"= c(0.001 , 0.2 , 0.001 , 1.1, 1, -2),
"pup"= c(0.3 , 0.8 , .1, 11 , 1e4, 10))
out <- shypEstFun(shpmodel = "01110",
parL = parL,
retdata = retdata, condata = condata,
ivap = ivap,
hclip = FALSE,
weightmethod = weightmethod,
LikModel = LikModel,
ALG = ALG,
set.itermax = set.itermax,
lhs.query = FALSE)
\dontshow{
}
\donttest{
data("shpdata1")
retdata <- ret <- shpdata1$TS1$wrc
condata <- con <- shpdata1$TS1$hcc
condata <- condata[!is.na(condata[,1]),]
---
# 1 SET VARIABLES --------------------
# VARIABLES FOR PLOTTING
{pF <- seq(-3, 6.8, length = 201)
h <- 10^pF
ticksatmin <- -2
tcllen <- 0.4
ticksat <- seq(ticksatmin,5,1)
pow <- ticksatmin:6
# VARIABLES FOR THE FITTING ALGORITHM
weightmethod = "range"
ivap = NULL
set.itermax = 3e1
LikModel = "rss" # ALTERNATIVE OPTION: LikModel = "-2logLik"
ALG = "DE" # ALTERNATIVE OPTION: ALG = "modMCMC"
shpmodel.v <- c("01110", "01110FM")
plot.query = FALSE
no.shps <- length(shpmodel.v)
# initialising lists
out.L <- vector("list", no.shps)
gof.L <- vector("list", no.shps)
}
# Run comparison
for (i in 1:2) {
shpmodel = shpmodel.v[i]
# INITIAL PARAMETERS, BOUNDS, and SELECTED PARAMETERS FOR FITTING
switch(shpmodel,
"01110" = {
# van Genuchten-Mualem Model parameters
parL<-list("p"=c("thr"=0.05,"ths"=0.45,"alf1"=0.01,"n"=2,"Ks"=100,"tau"=.5),
"psel" = c(1, 1, 1, 1, 1, 1),
"plo"= c(0.001 , 0.2 , 0.001 , 1.1, 1, -2),
"pup"= c(0.3 , 0.8 , .1, 11 , 1e4, 10)
)
},
"01110FM" = {
# van Genuchten-Mualem Model parameters + BRUNSWICK MODEL
parL<-list("p"=c("thr"=0.05,"ths"=0.45,"alf1"=0.01,"n"=2,"Ksc"=100,
"tau"=.5,"Ksnc"=1e-4,"a"=1.5,"h0"=6.8),
"psel" = c( 1,1, 1 ,1 , 1,1,1, 0, 0),
"plo"= c(0.001 , 0.1 , 0.001 , 1.1, 1,0,1e-6 , 1, 6.5),
"pup"= c(0.35, 0.7 , .1, 11 , 1e4,10 ,1e0, 2, 6.9)
)
},
stop("Enter a meaningful shpmodel")
)
out <- shypEstFun(shpmodel = shpmodel,
parL = parL,
retdata = retdata, condata = condata,
ivap = ivap,
hclip = FALSE,
weightmethod = weightmethod,
LikModel = LikModel,
ALG = ALG,
set.itermax = set.itermax
,lhs.query = FALSE)
out$model <- shpmodel.v[i]
out.L[[i]] <- out
# Calculate the soil hydraulic properties for the plot
if(ALG == "DE"){
p <- out$out$optim$phattrans
}
if(ALG == "modMCMC"){
p <- out$out$bestpartrans
}
if(weightmethod =="est1"){
np <- length(p)
p <- p[-c(np-1, np)]
if(ALG =="modMCMC"){
parL$p[which(parL$psel==1)] <- p
p <- parL$p
}
}
if(plot.query==TRUE){
shyp.L<-shypFun(p,h,shpmodel=shpmodel.v[i],ivap.query=ivap)
if(shpmodel == c("01110")){
wrc<-shyp.L$theta
hcc<-log10(shyp.L$Kh)
# PLOT THE WATER RETENTION CURVE
par(mfrow=c(1,2),tcl=tcllen)
plot(retdata,ylim=c(.0,.50),xlim=c(0,6.8),ylab="",xlab="",
col="darkgrey",axes=FALSE,main="WaterRetentionCurve",cex=2)
lines(log10(abs(h)),wrc,col="darkblue",lwd=2)
legend("topright",c("observed","simulated"),pch=c(1,NA),
lty=c(NA,1),lwd=2,bty="n",cex=1.3,col=c("darkgrey","darkblue"))
axis(1,at=pow,labels=parse(text=paste('10^',(pow),sep="")),tcl=tcllen)
axis(2,tcl=tcllen)
axis(4,labels=NA)
axis(3,labels=NA)
mtext("pressurehead|h|[cm]",1,cex=1.2,line=2.8)
mtext("vol.watercontent[-]",2,cex=1.2,line=2.8)
box()
# PLOT THE MEASURED HYDRAULIC CONDUCTIVITY CURVE
plot(condata,ylim=c(-8,2),xlim=c(0,6.8),ylab="",xlab="",col="darkgrey",
axes=FALSE,main="HydraulicConductivityCurve",cex=2)
lines(log10(abs(h)),hcc,col="darkblue",lwd=2)
legend("topright",c("observed","simulated"),pch=c(1,NA),
lty=c(NA,1),lwd=2,bty="n",cex=1.3,col=c("darkgrey","darkblue"))
axis(1,at=pow,labels=parse(text=paste('10^',(pow),sep="")),tcl=tcllen)
axis(2)
axis(4,labels=NA)
axis(3,labels=NA)
mtext("log10K[cm/d]",2,cex=1.2,line=2.8)
mtext("pressurehead|h|[cm]",1,cex=1.2,line=2.8)
box()
par(mfrow=c(1,1))
}
if(shpmodel == "01110FM"){
wrc<-shyp.L$theta
wrccap<-shyp.L$thetacap
wrcnc<-shyp.L$thetanc
hcc<-log10(shyp.L$Kh)
hcccap<-log10(shyp.L$Kcap)
hccnc<-log10(shyp.L$Knc)
hcvap<-log10(shyp.L$Kvap)
par(mfrow=c(1,2),tcl=tcllen)
plot(retdata,ylim=c(.0,.50),xlim=c(0,6.8),ylab="",xlab="",
col="darkgrey",axes=FALSE,main="WaterRetentionCurve",cex=2)
lines(log10(h),wrccap,col="brown",lwd=2)
lines(log10(h),wrcnc,col="brown",lwd=2,lty=2)
lines(log10(h),wrc,col="darkblue",lwd=2)
legend("topright",c("observed","simulated"),pch=c(1,NA),
lty=c(NA,1),lwd=2,bty="n",cex=1.3,col=c("darkgrey","darkblue"))
axis(1,at=pow,labels=parse(text=paste('10^',(pow),sep="")),tcl=tcllen)
axis(2,tcl=tcllen)
axis(4,labels=NA)
axis(3,labels=NA)
mtext("pressurehead|h|[cm]",1,cex=1.2,line=2.8)
mtext("vol.watercontent[-]",2,cex=1.2,line=2.8)
box()
# PLOT THE HYDRAULIC CONDUCTIVITY CURVE
plot(condata,ylim=c(-8,max(max(condata[,1]),max(hcc)))
,xlim=c(0,6.8),ylab="",xlab="",col="darkgrey",
axes=FALSE,main="HydraulicConductivityCurve",cex=2)
lines(log10(h),hcccap,col="brown",lwd=2)
lines(log10(h),hccnc,col="brown",lwd=2,lty=2)
lines(log10(h),hcc,col="darkblue",lwd=2)
lines(log10(h),hcvap,col="darkblue",lwd=2)
legend("topright",c("observed","simulated"),pch=c(1,NA),
lty=c(NA,1),lwd=2,bty="n",cex=1.3,col=c("darkgrey","darkblue"))
axis(1,at=pow,labels=parse(text=paste('10^',(pow),sep="")),tcl=tcllen)
axis(2)
axis(4,labels=NA)
axis(3,labels=NA)
mtext("log10K[cm/d]",2,cex=1.2,line=2.8)
mtext("pressurehead|h|[cm]",1,cex=1.2,line=2.8)
box()
par(mfrow=c(1,1))
}
}
phattrans.m <- out$out$optim$phattrans
gof.L[[i]]<-gofFun(phattrans.m,shpmodel=shpmodel.v[i],retdata=retdata,condata=condata,
out.L[[i]]$settings$weight,parL$psel,ivap.query=NULL,hclip.query=FALSE)
}
statstab3 <- cbind("th_rmse" = unlist(lapply(lapply(gof.L, `[[`, "th"), '[[', "rmse")),
"log10Kh_rmse" = unlist(lapply(lapply(gof.L, `[[`, "log10Kh"), '[[', "rmse"))
)
}
# }
Run the code above in your browser using DataLab