################################################################################
## File hydroPSO-Rscript.R                                                     #
## Part of the hydroPSO R package:                                             #
##                             http://www.rforge.net/hydroPSO/ ;               #
##                             http://cran.r-project.org/web/packages/hydroPSO #
## Copyright 2011-2012 Mauricio Zambrano-Bigiarini & Rodrigo Rojas             #
## Distributed under GPL 2 or later                                            #
##                                                                             #
## R script created with he PEST2hydroPSO function                             #
##                                                                             #
## Created by Mauricio Zambrano-Bigiarini and Rodrigo Rojas. 08-Nov-2012       #
## Last saved: 08-Nov-2012                                                     #
################################################################################

###Loading required libraries
set.seed(2000)
library(hydroPSO)
library(hydroGOF)
library(hydroTSM)
source("read_output.R", echo=TRUE)

################################################################################
##################   User-defined variables - START  ###########################
################################################################################

###Definition of working directory: input, output and model files paths
model.drty <- "Phreeqc_hydroPSO"
model.drty <- getwd()

###Period of analysis. 
###Only required for out.FUN
###To define subperiod of analysis to compute the goodness-of-fit (GoF) measures
Sim.Ini <- "YYYY-MM-DD"
Sim.Fin <- "YYYY-MM-DD"
gof.Ini <- "YYYY-MM-DD"
gof.Fin <- "YYYY-MM-DD"


###Function for reading the simulated equivalents
out.FUN="read_output"   # name of user-defined function for reading model outputs
out.FUN.args=list( 							###START out.FUN.args 
             fname="U.prn",   # name of the output file
	      nobs=6)
             #,,,                  # additional arguments for 'out.FUN'
              ###END out.FUN.args
			  
###Goodness-of-fit function, either pre-defined from hydroGOF (e.g., ssq) or 
###customized
gof.FUN <- "ssq" # sum of squared residuals. PEST default. 
                 # any other model performance measure could be used
gof.FUN.args <- list()

### Parallel computations ?
  #parallel <- "none"            # No parallel computations
  #parallel <- "multicore"      # GNU/Linux, UNIX only. For multi-core machines only
  #parallel <- "parallel"       # GNU/Linux, UNIX only. For multi-core machines and network clusters
  parallel <- "parallelWin"    # Windows & other OS.   For multi-core machines and network clusters

  # Number of cores/nodes desired
  par.nnodes <- 5

  # OPTIONAL. Only required for parallel=="parallelWin"   
  # packages that have to be loaded in each core or node of the network cluster
  par.pkgs   <- c("hydroGOF", "hydroTSM") 
################################################################################
###################   User-defined variables - END  ############################
################################################################################

###Getting the OBSERVATIONS
obs.fname <- "PEST2hydroPSO_OBS.txt"
obs.fname <- paste(model.drty, "/PSO.in/", obs.fname, sep="")
obs 	<- read.table(obs.fname)[,1]

###MAIN model function
model.FUN.args=list(
   model.drty=model.drty, 
   param.files=paste(model.drty,"/PSO.in/ParamFiles.txt",sep=""), 
exe.fname= "p.bat",
   #stdout="",            
   out.FUN=out.FUN,  
   out.FUN.args=out.FUN.args,
   gof.FUN=gof.FUN,
   gof.FUN.args=gof.FUN.args,
   #gof.Ini=gof.Ini,        # un-comment if 'gof.Ini' is defined
   #gof.Fin=gof.Fin,        # un-comment if 'gof.Fin' is defined
   obs=obs
) ###END model.FUN.args

### MAIN PSO ALGORITHM
### For hydroPSO fine-tuning parameters, see Zambrano-Bigiarini and Rojas, 2012
hydroPSO(
   fn="hydromod",
   model.FUN="hydromod",        
   model.FUN.args=model.FUN.args,
   control=list(
      MinMax="min",                                                         ### Maximization of the objective function (Gaussian Likelihood).
      npart=10,                                                             ### Number of particles.
      maxit=200,                                                           ### Number of iterations.
      reltol=1e-30,                                                         ### Numeric tolerance to stop successive hydroPSO iterations. 
      c1=2.05,                                                              ### Acceleration coefficient c1.
      c2=2.05,                                                              ### Acceleration coefficient c2.
      use.IW=FALSE,                                                         ### Parameters for defining the inertia weight scheme employed.
	  use.CF=TRUE,                                                          ### Using the Clerc's contriction factor
      use.TVlambda=TRUE,TVlambda.type="linear",TVlambda.rng=c(1.0,0.5),TVlambda.exp=1,  ### Parameters for defining a time-varying velocity constraint lambda.
	  topology="random", K=5,                                              ### Parameters for defining a random topology with K informants
	  drty.out="SPSO2011_i31.parallel.out",                                          ### Directory that will store the results of hydroPSO.
      boundary.wall="reflecting",
      normalise=TRUE,
      REPORT=1,                                                             ### Interval to report messages to screen.
### control for multicore/Parallel options
      parallel="parallelWin",
      par.nnodes=5,
      par.pkgs=c("hydroGOF","hydroTSM")
###
    ) ##END control options
) ##END MAIN hydroPSO ALGORITHM
# Plotting the results
plot_results(drty.out="SPSO2011_i31.parallel.out", MinMax="min", do.png=TRUE, do.pairs=TRUE, ,gof.name="SSR")
