R/wacs.simul_innercall.R

Defines functions wacs.simul_innercall

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

Try the WACS package in your browser

Any scripts or data that you put into this service are public.

WACS documentation built on July 1, 2020, 5:22 p.m.