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

[Vishal Mehta]

Thanks Karline,

I did attach the code; am new to this list, perhaps it does not accept attachments? I am attaching again, but also just copy-pasting it below: If you could take a look, I would greatly appreciate it.

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


-----Original Message-----
From: R-sig-dynamic-models [mailto:r-sig-dynamic-models-bounces at r-project.org] On Behalf Of Karline Soetaert
Sent: Tuesday, July 19, 2016 11:32 PM
To: Special Interest Group for Dynamic Simulation Models in R
Subject: Re: [R-sig-dyn-mod] Qn about FME package

Vishal, your code was not attached, so I can't have a look. 
However, this usually means that your sesnFun function does not return useful information, i.e. variables whose values are changed when parameters are altered.

There are several vignettes in the FME package that may help you. Try: 
vignette("FMEsteady")
vignette("FMEdyna")
vignette("FMEother")
vignette("FMEmcmc")

The first uses a 1-D steady-state model.

Hope it helps,

Karline

-----Original Message-----
From: R-sig-dynamic-models [mailto:r-sig-dynamic-models-bounces at r-project.org] On Behalf Of Vishal Mehta
Sent: maandag 18 juli 2016 23:28
To: r-sig-dynamic-models at r-project.org
Subject: [R-sig-dyn-mod] Qn about FME package

Dear list,

I've been trying to use FME for sensitivity analysis, parameter estimation and uncertainty analysis, using a simple 1d groundwater model in R (for a groundwater study in India).

I've been trying to follow along the JSS article ("Inverse modeling, sensitivity and monte carlo analysis in R using package FME"), using my model: but I have gotten stuck at the sensFun step: I am not able to get a sensitivity analysis (it comes out all zeros). I am guessing parameter values are not getting perturbed.

Does it work only for differential equations?

I am attaching the R code..if at all you can help, I would really appreciate it! Note, one of the codes is to calculate parameter sensitivity, the other was to try and calibrate (model fit)

Regards,


Vishal K. Mehta<http://sei-us.org/about/staff_person/19>, Ph.D.
Senior Scientist,
Stockholm Environment Institute-US
400 F St, Davis, CA 95616
+ 1 530 753 3035   x3
http://sei-us.org/



_______________________________________________
R-sig-dynamic-models mailing list
R-sig-dynamic-models at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models

_______________________________________________
R-sig-dynamic-models mailing list
R-sig-dynamic-models at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models