[R-sig-dyn-mod] Qn about FME package

Thomas Petzoldt thomas.petzoldt at tu-dresden.de
Wed Jul 20 21:12:42 CEST 2016


I guess the error may be in:

Sfun=sensFun(gwcost,parms=pars,sensvar=c("h"))

See example of ?sensFun, where func is a wrapper around the model that 
returns the modelled (i.e. dependent or state) variables. In the 
statement above, it was called with a cost function.


Thomas P.





Am 20.07.2016 um 18:02 schrieb Vishal Mehta:
> #########################################
> #Jan 2016: Dr. Satkumar Tomer's 1-D model copied
> ####GW model 1-D GW model function
> gw_model <- function(pars, hmin.under.flow, pd.under.flow, rain){
> 	n <- length(rain)
> 	net_recharge <- rf*rain
> 	h <- array(0,n)
> 	discharge <- array(0,n)
> 	time = array(0,n)
> 	h[1] <- hini
> 	time[1] = 1
> 	for (i in 1:(n-1)){
> 	discharge[i] <- sy*(1-pd.under.flow)*(h[i]-hmin.under.flow+ net_recharge[i]/sy)
> 	if (discharge[i]<0){
> 		discharge[i] <- 0
>         }
> 	h[i+1] <- h[i] + (net_recharge[i]-discharge[i])/sy
>       time[i+1] = time[i] + 1
> }
> 	return(data.frame(time=time,h=h,discharge=discharge))
> }
> #synthetic forcing
> n = 12*20 # temporal length in months
> set.seed(123) # to have the repeatibility in results
> rain <- rgamma(n, shape=8, rate = 0.1) # synthetic rainfall
> #barplot(rain,xlab="Time",ylab="Rainfall (mm)",main="Monthly rainfall")
> rain=rain/1000#mm to m of rain
> #new
>
>
> #
> #input parameters
> rf <- 0.1 # recharge factor, fraction
> sy <- 0.005 #specific yield, dimensionless
> hini <- 150 # initial condition, metres
> hmin.under.flow <- 100
> pd.under.flow<-0.967
> pars=c("hini"=hini,"rf"=rf,"sy"=sy)##this is the subset of parameters to be SA'd, calibrated, and to inform UA
> ##
> pdf("output.pdf")
> output=gw_model(pars,hmin.under.flow,pd.under.flow,rain)
> plot(output$h,type='l',xlab='time',ylab='simulated gw level(m)')
> dev.off()
> ##########################
> #July 15, 2016
> #trial following FME tutorial
> #########
> library(FME)
> #Create dummy observations
> #n=240 length of simulated output
> #output$h contains the simulated head
> #below, observation dataframe is called DataT. It has to(?) contain a "time" and "h" variable like the simulation output
> DataT=data.frame(cbind("time"=seq(1,length(n)),"h"=output$h+rnorm(sd=0.5,mean=5,n=240),"sd"=4.5))
> DataT2=DataT[,c(1,2)]
> ##
> #Create a function to calculate fit
> gwcost=function(pars){
> 	out=gw_model(pars,hmin.under.flow, pd.under.flow, rain)
> 	#cost=modCost(model=out,obs=DataT,err="sd")
> 	cost=modCost(model=out,obs=DataT2)
> 	return(cost=cost)
> 	}
> ###Above seems to be working(maybe not quite?)
> gwcost(pars)
> gwcost(pars)$model
> #[1] 1075 with no weighting
> plot(gwcost(pars)$residuals$res)
> #
> ##Sensitivity using sensFun (here onwards results are not following the JSS article)
> ##local sensitivity does not seem to be working. does it work only with differential equations?
> Sfun=sensFun(gwcost,parms=pars,sensvar=c("h"))
>
> #parameter tinyy to add perturbation? not making any difference
> summary(Sfun)
> plot(Sfun)
> #########



More information about the R-sig-dynamic-models mailing list