Nothing
###################################################################
#
# This function is part of WACSgen V1.0
# Copyright © 2013,2014,2015, D. Allard, BioSP,
# and Ronan Trépos MIA-T, INRA
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 2
# of the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details. http://www.gnu.org
#
###################################################################
wacs.simul_innercall=function(wacspar, from, to, first.day=NULL, REJECT=FALSE){
###################################################################
#
# WACSsimul_innercall : Main function for performing the simulations of the WACSgen simulator
#
# ARGUMENTS
# wacspar: Parameters, as created by WACSestim
# from : Starting date of the simulation (format: "yyyy-mm-dd")
# to : Ending date of the simuation (format: "yyyy-mm-dd")
# first.day : Conditioning values for first day (optional)
# REJECT : Boolean indicating wether a rejection technique is used to guarantee variables within bounds
# Default is FALSE. In this case, values outside bounds are forced the the bounds.
# VALUE
# A data frame with simulation results. Column 1-3 are year,month,day.
# Column 4 is rain; The, other variables
#
# CALLS
# map.wt
# rmcsnstar
# testvariables
# transform.rain
# rMarkov
# rmcsn.cond
# other functions from wacs.datefunctions.R
#
###################################################################
# Checking wacspar
if (class(wacspar) != "WACSpar") {
stop ("[WACSsimul] Parameters should be of class 'WACSpar', as generated by calling WACSestim")
}
###############################################
#
# Initialization
#
###############################################
##
Trange = wacspar$Trange;
seasons = wacspar$seasons;
bounds = wacspar$bounds;
Deviation = wacspar$Trend$Deviation;
Central = wacspar$Trend$Central;
RainModel = wacspar$Rain$RainModel;
RainPar = wacspar$Rain$RainPar
Nv = length(wacspar$varnames); #number of variables (including rain)
iyear =1; #indexes in global matrix
imonth =2;
iday =3;
iseason =4;
iWT =5;
irain =6;
iirain =1; #indexes in variable vector
iitmin =2;
iitrange=-1;
iitmax =-1;
if (Trange) {
iitrange = 3
} else {
iitmax=3;
}
#initialize output matrices
Nd = as.Date(to) - as.Date(from) +1; # number of days to simulate
Xsimul = matrix(0,Nd,(5+length(wacspar$varnames)))
Xresid = matrix(0,Nd,length(wacspar$varnames))
###############################################
#
# Simulating day 1
#
###############################################
j = 1; #day of simulation
currDate = as.Date(from);
if ((wacs.month(currDate) == 2) && (wacs.day(currDate) == 29)) {
stop ("[WACSsimul] cannot start at date 02-29")
}
currSeason = wacs.season(wacs.month(currDate),wacs.day(currDate),seasons);
doy = wacs.dayOfYear(currDate); #day of year
if (wacs.leapYear(wacs.year(currDate)) && doy >= 60) {
jd = doy -1;
} else {
jd = doy;
}
Xsimul[j,iyear] = wacs.year(currDate);
Xsimul[j,imonth] = wacs.month(currDate);
Xsimul[j,iday] = wacs.day(currDate);
Xsimul[j,iseason] = currSeason;
seasonPar = wacspar[[paste("Season",sep="_",currSeason)]];
Ndry = seasonPar$NumbWT[1];
Nwet = seasonPar$NumbWT[2];
limit.prob = abs(eigen(t(seasonPar$TransM))$vectors[,1]);
limit.prob = limit.prob/sum(limit.prob);
if (is.null(first.day)){
# Simulating weather type for day 1
wt = sample.int((Ndry+Nwet),size=1,prob=limit.prob);
Xsimul[j,iWT] = wt;
# Setting residual parameters for day 1
WET = (wt > Ndry)
vdry = 1-WET
Nvv = Nv - vdry
mu = seasonPar[[3 + wt]]$loc
sigma = seasonPar[[3 + wt]]$cov
skew = seasonPar[[3 + wt]]$skew
rho = seasonPar[[3 + wt]]$rho
MU = rep(mu,2)
SKEW = rep(skew,2)
SS = matrix(0,2*Nvv,2*Nvv)
SS[1:Nvv,1:Nvv] = sigma
SS[(Nvv+1):(2*Nvv),(Nvv+1):(2*Nvv)] = sigma
SS[1:Nvv,(Nvv+1):(2*Nvv)] = sqrtm(sigma)%*%diag(rho)%*%sqrtm(sigma)
SS[(Nvv+1):(2*Nvv),1:Nvv] = t(SS[1:Nvv,(Nvv+1):(2*Nvv)])
# Simulating residuals for day 1
COMP = FALSE
Xtest = rep(0,Nv-1)
while(!COMP){
#print (paste(" ! comp ", wacs.day(currDate) ));
Xresid[j,(1+vdry):Nv] = as.vector(rmcsnstar(1,MU,SS,SKEW)[1:Nvv])
for (v in 1:(Nv-1)){ #new
Xtest[v] = Xresid[j,v+1]*Deviation[jd,v] + Central[jd,v];
}
if (WET){
Xsimul[j,irain] = transform.rain(Xresid[iirain],rain.model=RainModel,
rain.par=RainPar[currSeason,])
}else{
Xsimul[j,irain]= 0.0
}
Xsimul[j,(irain+1):ncol(Xsimul)] = Xtest;
COMP = testvariables(Xsimul[j,irain:ncol(Xsimul)], bounds, iitmin, iitrange)
}
}
else{ # First day is given; therefore we compute the associated first
# Residual and set the first Simulated value
s = currSeason
if(RainModel=="Gamma"){
Xresid[1] = stats::qnorm(stats::pgamma(first.day[1],scale=RainPar[s,1],
shape=RainPar[s,2]))
}
if(RainModel=="None"){
Xresid[1] = first.day[1]
}
for (v in 2:Nv){
Xresid[v] = (first.day[v] -Central[jd,(v-1)]) / Deviation[jd,(v-1)];
}
Xsimul[1,iWT] = map.wt(Xresid[1,],seasonPar)
Xsimul[1,(iWT+1):(iWT+Nv)] = first.day
}
##################################
#
# Loop on all days:
#
# 1/ sample a new weather type (calls function rMarkov)
# 2/ conditionnaly to this weather type, sample a new vector of variables
# (calls function rvariables)
# 3/ test wether these variables are within vbounds (calls function testvariables);
# if yes OK; if not go to 2 (uses a while loop)
#
# newWT is a boolean; if TRUE must sample a weather type
# KRITER is a boolean; if TRUE the simulated variables are OK
#
#
##################################
Nhelp = 0
while (j < Nd){
j = j + 1;
if (j%%100 == 0) {
print(paste("day :",j,sep=""))
}
currDate = currDate + 1;
currSeason = wacs.season(wacs.month(currDate),wacs.day(currDate),seasons);
Xsimul[j,iyear] = wacs.year(currDate);
Xsimul[j,imonth] = wacs.month(currDate);
Xsimul[j,iday] = wacs.day(currDate);
Xsimul[j,iseason] = currSeason;
feb29 = ((Xsimul[j,imonth] == 2) && (Xsimul[j,iday] == 29));
if (!feb29) {
jd = jd+1;
if (jd == 366) {
jd = 1;
}
}
oldWT = Xsimul[j-1,iWT];
if (Xsimul[j-1,iseason] != currSeason) {
seasonPar = wacspar[[paste("Season",sep="_",currSeason)]];
Ndry = seasonPar$NumbWT[1];
Nwet = seasonPar$NumbWT[2];
oldWT = map.wt(Xresid[j-1,],seasonPar)
}
WET.old = (oldWT > Ndry)
vdry.old = 1-WET.old
mu.old = seasonPar[[3 + oldWT]]$loc
sigma.old = seasonPar[[3 + oldWT]]$cov
skew.old = seasonPar[[3 + oldWT]]$skew
rho.old = seasonPar[[3 + oldWT]]$rho
Nvv.old = Nv - vdry.old
COMP = FALSE;
Ntry = 0;
while(!COMP){
Ntry = Ntry+1;
# Simulating weather type for new day
wt = rMarkov(oldWT,seasonPar$TransM);
Xsimul[j,iWT] = wt
# new parameters
WET.new = (wt > Ndry);
vdry.new = 1-WET.new;
mu.new = seasonPar[[3 + wt]]$loc;
sigma.new = seasonPar[[3 + wt]]$cov;
skew.new = seasonPar[[3 + wt]]$skew;
rho.new = seasonPar[[3 + wt]]$rho;
Nvv.new = Nv - vdry.new;
# building the covariance matrix of (X_t, X_(t+1))
SS = matrix(0,(Nvv.old+Nvv.new),(Nvv.old+Nvv.new))
SS[1:Nvv.old,1:Nvv.old] = sigma.old
SS[(Nvv.old+1):(Nvv.old+Nvv.new),(Nvv.old+1):(Nvv.old+Nvv.new)] = sigma.new;
rho = NULL;
rho.mat = NULL;
if (WET.old == WET.new){
rho = apply(cbind(rho.old,rho.new),1,max);
rho.mat = diag(rho);
}
if (WET.old & !WET.new){
rho = apply(cbind(rho.old[-1],rho.new),1,max);
rho.mat = diag(rho);
rho.mat = rbind(rep(0,length(rho)),rho.mat);
}
if (!WET.old & WET.new){
rho = apply(cbind(rho.old,rho.new[-1]),1,max);
rho.mat = diag(rho);
rho.mat = cbind(rep(0,length(rho)),rho.mat);
}
SSdiag = sqrtm(sigma.old)%*%rho.mat%*%sqrtm(sigma.new);
SS[1:Nvv.old,(Nvv.old+1):(Nvv.old+Nvv.new)] = SSdiag;
SS[(Nvv.old+1):(Nvv.old+Nvv.new),1:Nvv.old] =
t(SS[1:Nvv.old,(Nvv.old+1):(Nvv.old+Nvv.new)]);
MU = c(mu.old,mu.new);
SKEW = c(skew.old,skew.new);
D = diag(SKEW)%*%solve(sqrtm(SS));
Delta = diag(Nvv.old+Nvv.new) - diag(SKEW^2);
# simulating residuals, conditional on the residuals of previous day
Xtilde.old = Xresid[j-1,];
if(!WET.old){
Xtilde.old = Xtilde.old[-1];
}
Xtilde.new = as.vector(rmcsn.cond(1,Xtilde.old,1:Nvv.old,MU,SS,D,
rep(0,Nvv.old+Nvv.new),Delta));
if (Ntry > 50){
Nhelp = Nhelp + 1;
Xtilde.new = as.vector(rmcsnstar(1,MU,SS,SKEW)[(Nvv.old+1):(Nvv.old+Nvv.new)]);
}
if (!WET.new){
Xtilde.new=c(-999,Xtilde.new);
}
#if (is.na(Xtilde.new[Nv])){
if(any(Xtilde.new[2:length(Xtilde.new)] > 3.5) | any(Xtilde.new[2:length(Xtilde.new)] < - 3.5) ){
COMP = FALSE;
}else{
Xresid[j,] = Xtilde.new;
if (WET.new){
Xsimul[j,irain] = transform.rain(Xtilde.new[iirain],
rain.model=RainModel,rain.par=RainPar[currSeason,])
if (Xtilde.new[iirain]>4.0) {
Xtilde.new[iirain]=4.0
}
}
Xtest = rep(0,Nv-1)
for (v in 1:(Nv-1)){
Xtest[v] = Xtilde.new[v+1]*Deviation[jd,v] + Central[jd,v]
}
Xsimul[j,(irain+1):ncol(Xsimul)] = Xtest;
if (REJECT){
COMP = testvariables(Xsimul[j,irain:ncol(Xsimul)],
bounds, iitmin, iitrange)
}else{
COMP=TRUE
}
}
} #closes loop on tries
if (!REJECT){
Xsimul[j,irain:ncol(Xsimul)] = boundsvariables(Xsimul[j,irain:ncol(Xsimul)],bounds,iitmin,iitrange)
}
} #closes loop on days
Xsimul = data.frame(Xsimul)
names(Xsimul) = c("year","month","day","season","WT",wacspar$varnames)
return(Xsimul)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.