##' @include hydroState.R
##' @export
hydroState.allModels <- setClass(
# Set the name for the class
"hydroState.allModels",
package='hydroState',
# Define the slots
slots = c(
siteID= 'character',
calib.reference.model.name = 'list',
calib.reference.criteria.met = 'logical',
models = 'list'
),
# Set the default values for the slots. (optional)
prototype=list(
siteID = '(not set)',
calib.reference.model.name= vector('list',1),
calib.reference.criteria.met= vector('logical',1),
models = vector('list',1)
)
)
# Valid object?
validObject <- function(object) {
TRUE
}
setValidity("hydroState.allModels", validObject)
# Initialise the object.
#setGeneric(name="initialize",def=function(.Object,input.data, Qhat.object, QhatModel.object, markov.model.object, ...){standardGeneric("initialize")})
setMethod(f="initialize",signature="hydroState.allModels",definition=function(.Object, siteID, input.data, allow.flickering=F, state.dependent.mean.trend=NA)
{
.Object@siteID <- siteID
# Define transition graphs
transition.graph.1State <- matrix(TRUE,1,1)
transition.graph.2State <- matrix(TRUE,2,2)
transition.graph.3State.unstructured <- matrix(TRUE,3,3)
transition.graph.3State <- matrix(TRUE,3,3)
transition.graph.3State[1,3] <- FALSE
transition.graph.3State[2,1] <- FALSE
transition.graph.3State[3,2] <- FALSE
# Build Qhat flow transform objects.
Qhat.boxcox = new('Qhat.boxcox', input.data=input.data)
Qhat.log = new('Qhat.log', input.data=input.data)
# Check if flickering between states is allows
model.extension=''
if (allow.flickering) {
model.extension.markov = '.flickering'
} else {
model.extension.markov = ''
}
# Build Qhat models
QhatModel.1State = new(paste('QhatModel.homo.normal.linear',model.extension,sep=''), input.data=input.data,state.dependent.mean.trend=state.dependent.mean.trend, transition.graph=transition.graph.1State)
QhatModel.1State.AR1 = new(paste('QhatModel.homo.normal.linear.AR1',model.extension,sep=''), input.data=input.data,state.dependent.mean.trend=state.dependent.mean.trend, transition.graph=transition.graph.1State)
QhatModel.1State.AR2 = new(paste('QhatModel.homo.normal.linear.AR2',model.extension,sep=''), input.data=input.data,state.dependent.mean.trend=state.dependent.mean.trend, transition.graph=transition.graph.1State)
QhatModel.1State.AR3 = new(paste('QhatModel.homo.normal.linear.AR3',model.extension,sep=''), input.data=input.data,state.dependent.mean.trend=state.dependent.mean.trend, transition.graph=transition.graph.1State)
QhatModel.1State.gamma = new(paste('QhatModel.homo.gamma.linear',model.extension,sep=''), input.data=input.data,state.dependent.mean.trend=state.dependent.mean.trend, transition.graph=transition.graph.1State)
QhatModel.1State.gamma.AR1 = new(paste('QhatModel.homo.gamma.linear.AR1',model.extension,sep=''), input.data=input.data,state.dependent.mean.trend=state.dependent.mean.trend, transition.graph=transition.graph.1State)
QhatModel.1State.gamma.AR2 = new(paste('QhatModel.homo.gamma.linear.AR2',model.extension,sep=''), input.data=input.data,state.dependent.mean.trend=state.dependent.mean.trend, transition.graph=transition.graph.1State)
QhatModel.1State.gamma.AR3 = new(paste('QhatModel.homo.gamma.linear.AR3',model.extension,sep=''), input.data=input.data,state.dependent.mean.trend=state.dependent.mean.trend, transition.graph=transition.graph.1State)
QhatModel.2State = new(paste('QhatModel.homo.normal.linear',model.extension,sep=''), input.data=input.data,state.dependent.mean.trend=state.dependent.mean.trend, transition.graph=transition.graph.2State)
QhatModel.2State.AR1 = new(paste('QhatModel.homo.normal.linear.AR1',model.extension,sep=''), input.data=input.data,state.dependent.mean.trend=state.dependent.mean.trend, transition.graph=transition.graph.2State)
QhatModel.2State.AR2 = new(paste('QhatModel.homo.normal.linear.AR2',model.extension,sep=''), input.data=input.data,state.dependent.mean.trend=state.dependent.mean.trend, transition.graph=transition.graph.2State)
QhatModel.2State.AR3 = new(paste('QhatModel.homo.normal.linear.AR3',model.extension,sep=''), input.data=input.data,state.dependent.mean.trend=state.dependent.mean.trend, transition.graph=transition.graph.2State)
QhatModel.2State.gamma = new(paste('QhatModel.homo.gamma.linear',model.extension,sep=''), input.data=input.data,state.dependent.mean.trend=state.dependent.mean.trend, transition.graph=transition.graph.2State)
QhatModel.2State.gamma.AR1 = new(paste('QhatModel.homo.gamma.linear.AR1',model.extension,sep=''), input.data=input.data,state.dependent.mean.trend=state.dependent.mean.trend, transition.graph=transition.graph.2State)
QhatModel.2State.gamma.AR2 = new(paste('QhatModel.homo.gamma.linear.AR2',model.extension,sep=''), input.data=input.data,state.dependent.mean.trend=state.dependent.mean.trend, transition.graph=transition.graph.2State)
QhatModel.2State.gamma.AR3 = new(paste('QhatModel.homo.gamma.linear.AR3',model.extension,sep=''), input.data=input.data,state.dependent.mean.trend=state.dependent.mean.trend, transition.graph=transition.graph.2State)
QhatModel.3State = new(paste('QhatModel.homo.normal.linear',model.extension,sep=''), input.data=input.data,state.dependent.mean.trend=state.dependent.mean.trend, transition.graph=transition.graph.3State)
QhatModel.3State.AR1 = new(paste('QhatModel.homo.normal.linear.AR1',model.extension,sep=''), input.data=input.data,state.dependent.mean.trend=state.dependent.mean.trend, transition.graph=transition.graph.3State)
QhatModel.3State.AR2 = new(paste('QhatModel.homo.normal.linear.AR2',model.extension,sep=''), input.data=input.data,state.dependent.mean.trend=state.dependent.mean.trend, transition.graph=transition.graph.3State)
QhatModel.3State.AR3 = new(paste('QhatModel.homo.normal.linear.AR3',model.extension,sep=''), input.data=input.data,state.dependent.mean.trend=state.dependent.mean.trend, transition.graph=transition.graph.3State)
QhatModel.3State.gamma = new(paste('QhatModel.homo.gamma.linear',model.extension,sep=''), input.data=input.data,state.dependent.mean.trend=state.dependent.mean.trend, transition.graph=transition.graph.3State)
QhatModel.3State.gamma.AR1 = new(paste('QhatModel.homo.gamma.linear.AR1',model.extension,sep=''), input.data=input.data,state.dependent.mean.trend=state.dependent.mean.trend, transition.graph=transition.graph.3State)
QhatModel.3State.gamma.AR2 = new(paste('QhatModel.homo.gamma.linear.AR2',model.extension,sep=''), input.data=input.data,state.dependent.mean.trend=state.dependent.mean.trend, transition.graph=transition.graph.3State)
QhatModel.3State.gamma.AR3 = new(paste('QhatModel.homo.gamma.linear.AR3',model.extension,sep=''), input.data=input.data,state.dependent.mean.trend=state.dependent.mean.trend, transition.graph=transition.graph.3State)
QhatModel.3State.US = new(paste('QhatModel.homo.normal.linear',model.extension,sep=''), input.data=input.data,state.dependent.mean.trend=state.dependent.mean.trend, transition.graph=transition.graph.3State.unstructured)
QhatModel.3State.US.AR1 = new(paste('QhatModel.homo.normal.linear.AR1',model.extension,sep=''), input.data=input.data,state.dependent.mean.trend=state.dependent.mean.trend, transition.graph=transition.graph.3State.unstructured)
QhatModel.3State.US.AR2 = new(paste('QhatModel.homo.normal.linear.AR2',model.extension,sep=''), input.data=input.data,state.dependent.mean.trend=state.dependent.mean.trend, transition.graph=transition.graph.3State.unstructured)
QhatModel.3State.US.AR3 = new(paste('QhatModel.homo.normal.linear.AR3',model.extension,sep=''), input.data=input.data,state.dependent.mean.trend=state.dependent.mean.trend, transition.graph=transition.graph.3State.unstructured)
QhatModel.3State.US.gamma = new(paste('QhatModel.homo.gamma.linear',model.extension,sep=''), input.data=input.data,state.dependent.mean.trend=state.dependent.mean.trend, transition.graph=transition.graph.3State.unstructured)
QhatModel.3State.US.gamma.AR1 = new(paste('QhatModel.homo.gamma.linear.AR1',model.extension,sep=''), input.data=input.data,state.dependent.mean.trend=state.dependent.mean.trend, transition.graph=transition.graph.3State.unstructured)
QhatModel.3State.US.gamma.AR2 = new(paste('QhatModel.homo.gamma.linear.AR2',model.extension,sep=''), input.data=input.data,state.dependent.mean.trend=state.dependent.mean.trend, transition.graph=transition.graph.3State.unstructured)
QhatModel.3State.US.gamma.AR3 = new(paste('QhatModel.homo.gamma.linear.AR3',model.extension,sep=''), input.data=input.data,state.dependent.mean.trend=state.dependent.mean.trend, transition.graph=transition.graph.3State.unstructured)
# Build Markov model object
markov.1State = new(paste('markov.annualHomogeneous',model.extension.markov,sep=''), transition.graph=transition.graph.1State)
markov.2State = new(paste('markov.annualHomogeneous',model.extension.markov,sep=''), transition.graph=transition.graph.2State)
markov.3State = new(paste('markov.annualHomogeneous',model.extension.markov,sep=''), transition.graph=transition.graph.3State)
markov.3State.US = new(paste('markov.annualHomogeneous',model.extension.markov,sep=''), transition.graph=transition.graph.3State.unstructured)
# Build hydrostate models
.Object@models <- list(
model.1State.log = new('hydroState', input.data=input.data, Qhat.object = Qhat.log, QhatModel.object = QhatModel.1State , markov.model.object=markov.1State),
model.1State.log.AR1 = new('hydroState',input.data=input.data, Qhat.object = Qhat.log, QhatModel.object = QhatModel.1State.AR1, markov.model.object=markov.1State),
model.1State.log.AR2 = new('hydroState',input.data=input.data, Qhat.object = Qhat.log, QhatModel.object = QhatModel.1State.AR2, markov.model.object=markov.1State),
model.1State.log.AR3 = new('hydroState',input.data=input.data, Qhat.object = Qhat.log, QhatModel.object = QhatModel.1State.AR3, markov.model.object=markov.1State),
model.1State.gamma.log = new('hydroState', input.data=input.data, Qhat.object = Qhat.log, QhatModel.object = QhatModel.1State.gamma , markov.model.object=markov.1State),
model.1State.gamma.log.AR1 = new('hydroState',input.data=input.data, Qhat.object = Qhat.log, QhatModel.object = QhatModel.1State.gamma.AR1, markov.model.object=markov.1State),
model.1State.gamma.log.AR2 = new('hydroState',input.data=input.data, Qhat.object = Qhat.log, QhatModel.object = QhatModel.1State.gamma.AR2, markov.model.object=markov.1State),
model.1State.gamma.log.AR3 = new('hydroState',input.data=input.data, Qhat.object = Qhat.log, QhatModel.object = QhatModel.1State.gamma.AR3, markov.model.object=markov.1State),
model.1State.BC = new('hydroState', input.data=input.data, Qhat.object = Qhat.boxcox, QhatModel.object = QhatModel.1State , markov.model.object=markov.1State),
model.1State.BC.AR1 = new('hydroState',input.data=input.data, Qhat.object = Qhat.boxcox, QhatModel.object = QhatModel.1State.AR1, markov.model.object=markov.1State),
model.1State.BC.AR2 = new('hydroState',input.data=input.data, Qhat.object = Qhat.boxcox, QhatModel.object = QhatModel.1State.AR2, markov.model.object=markov.1State),
model.1State.BC.AR3 = new('hydroState',input.data=input.data, Qhat.object = Qhat.boxcox, QhatModel.object = QhatModel.1State.AR3, markov.model.object=markov.1State),
model.1State.gamma.BC = new('hydroState', input.data=input.data, Qhat.object = Qhat.boxcox, QhatModel.object = QhatModel.1State.gamma , markov.model.object=markov.1State),
model.1State.gamma.BC.AR1 = new('hydroState',input.data=input.data, Qhat.object = Qhat.boxcox, QhatModel.object = QhatModel.1State.gamma.AR1, markov.model.object=markov.1State),
model.1State.gamma.BC.AR2 = new('hydroState',input.data=input.data, Qhat.object = Qhat.boxcox, QhatModel.object = QhatModel.1State.gamma.AR2, markov.model.object=markov.1State),
model.1State.gamma.BC.AR3 = new('hydroState',input.data=input.data, Qhat.object = Qhat.boxcox, QhatModel.object = QhatModel.1State.gamma.AR3, markov.model.object=markov.1State),
model.2State.log = new('hydroState', input.data=input.data, Qhat.object = Qhat.log, QhatModel.object = QhatModel.2State , markov.model.object=markov.2State),
model.2State.log.AR1 = new('hydroState',input.data=input.data, Qhat.object = Qhat.log, QhatModel.object = QhatModel.2State.AR1, markov.model.object=markov.2State),
model.2State.log.AR2 = new('hydroState',input.data=input.data, Qhat.object = Qhat.log, QhatModel.object = QhatModel.2State.AR2, markov.model.object=markov.2State),
model.2State.log.AR3 = new('hydroState',input.data=input.data, Qhat.object = Qhat.log, QhatModel.object = QhatModel.2State.AR3, markov.model.object=markov.2State),
model.2State.gamma.log = new('hydroState', input.data=input.data, Qhat.object = Qhat.log, QhatModel.object = QhatModel.2State.gamma , markov.model.object=markov.2State),
model.2State.gamma.log.AR1 = new('hydroState',input.data=input.data, Qhat.object = Qhat.log, QhatModel.object = QhatModel.2State.gamma.AR1, markov.model.object=markov.2State),
model.2State.gamma.log.AR2 = new('hydroState',input.data=input.data, Qhat.object = Qhat.log, QhatModel.object = QhatModel.2State.gamma.AR2, markov.model.object=markov.2State),
model.2State.gamma.log.AR3 = new('hydroState',input.data=input.data, Qhat.object = Qhat.log, QhatModel.object = QhatModel.2State.gamma.AR3, markov.model.object=markov.2State),
model.2State.BC = new('hydroState', input.data=input.data, Qhat.object = Qhat.boxcox, QhatModel.object = QhatModel.2State , markov.model.object=markov.2State),
model.2State.BC.AR1 = new('hydroState',input.data=input.data, Qhat.object = Qhat.boxcox, QhatModel.object = QhatModel.2State.AR1, markov.model.object=markov.2State),
model.2State.BC.AR2 = new('hydroState',input.data=input.data, Qhat.object = Qhat.boxcox, QhatModel.object = QhatModel.2State.AR2, markov.model.object=markov.2State),
model.2State.BC.AR3 = new('hydroState',input.data=input.data, Qhat.object = Qhat.boxcox, QhatModel.object = QhatModel.2State.AR3, markov.model.object=markov.2State),
model.2State.gamma.BC = new('hydroState', input.data=input.data, Qhat.object = Qhat.boxcox, QhatModel.object = QhatModel.2State.gamma , markov.model.object=markov.2State),
model.2State.gamma.BC.AR1 = new('hydroState',input.data=input.data, Qhat.object = Qhat.boxcox, QhatModel.object = QhatModel.2State.gamma.AR1, markov.model.object=markov.2State),
model.2State.gamma.BC.AR2 = new('hydroState',input.data=input.data, Qhat.object = Qhat.boxcox, QhatModel.object = QhatModel.2State.gamma.AR2, markov.model.object=markov.2State),
model.2State.gamma.BC.AR3 = new('hydroState',input.data=input.data, Qhat.object = Qhat.boxcox, QhatModel.object = QhatModel.2State.gamma.AR3, markov.model.object=markov.2State),
model.3State.log = new('hydroState', input.data=input.data, Qhat.object = Qhat.log, QhatModel.object = QhatModel.3State , markov.model.object=markov.3State),
model.3State.log.AR1 = new('hydroState',input.data=input.data, Qhat.object = Qhat.log, QhatModel.object = QhatModel.3State.AR1, markov.model.object=markov.3State),
model.3State.log.AR2 = new('hydroState',input.data=input.data, Qhat.object = Qhat.log, QhatModel.object = QhatModel.3State.AR2, markov.model.object=markov.3State),
model.3State.log.AR3 = new('hydroState',input.data=input.data, Qhat.object = Qhat.log, QhatModel.object = QhatModel.3State.AR3, markov.model.object=markov.3State),
#
model.3State.gamma.log = new('hydroState', input.data=input.data, Qhat.object = Qhat.log, QhatModel.object = QhatModel.3State.gamma , markov.model.object=markov.3State),
model.3State.gamma.log.AR1 = new('hydroState',input.data=input.data, Qhat.object = Qhat.log, QhatModel.object = QhatModel.3State.gamma.AR1, markov.model.object=markov.3State),
model.3State.gamma.log.AR2 = new('hydroState',input.data=input.data, Qhat.object = Qhat.log, QhatModel.object = QhatModel.3State.gamma.AR2, markov.model.object=markov.3State),
model.3State.gamma.log.AR3 = new('hydroState',input.data=input.data, Qhat.object = Qhat.log, QhatModel.object = QhatModel.3State.gamma.AR3, markov.model.object=markov.3State),
model.3State.BC = new('hydroState', input.data=input.data, Qhat.object = Qhat.boxcox, QhatModel.object = QhatModel.3State , markov.model.object=markov.3State),
model.3State.BC.AR1 = new('hydroState',input.data=input.data, Qhat.object = Qhat.boxcox, QhatModel.object = QhatModel.3State.AR1, markov.model.object=markov.3State),
model.3State.BC.AR2 = new('hydroState',input.data=input.data, Qhat.object = Qhat.boxcox, QhatModel.object = QhatModel.3State.AR2, markov.model.object=markov.3State),
model.3State.BC.AR3 = new('hydroState',input.data=input.data, Qhat.object = Qhat.boxcox, QhatModel.object = QhatModel.3State.AR3, markov.model.object=markov.3State),
#
model.3State.gamma.BC = new('hydroState', input.data=input.data, Qhat.object = Qhat.boxcox, QhatModel.object = QhatModel.3State.gamma , markov.model.object=markov.3State),
model.3State.gamma.BC.AR1 = new('hydroState',input.data=input.data, Qhat.object = Qhat.boxcox, QhatModel.object = QhatModel.3State.gamma.AR1, markov.model.object=markov.3State),
model.3State.gamma.BC.AR2 = new('hydroState',input.data=input.data, Qhat.object = Qhat.boxcox, QhatModel.object = QhatModel.3State.gamma.AR2, markov.model.object=markov.3State),
model.3State.gamma.BC.AR3 = new('hydroState',input.data=input.data, Qhat.object = Qhat.boxcox, QhatModel.object = QhatModel.3State.gamma.AR3, markov.model.object=markov.3State),
model.3StateUS.log = new('hydroState', input.data=input.data, Qhat.object = Qhat.log, QhatModel.object = QhatModel.3State , markov.model.object=markov.3State.US),
model.3StateUS.log.AR1 = new('hydroState',input.data=input.data, Qhat.object = Qhat.log, QhatModel.object = QhatModel.3State.AR1, markov.model.object=markov.3State.US),
model.3StateUS.log.AR2 = new('hydroState',input.data=input.data, Qhat.object = Qhat.log, QhatModel.object = QhatModel.3State.AR2, markov.model.object=markov.3State.US),
model.2StateUS.log.AR3 = new('hydroState',input.data=input.data, Qhat.object = Qhat.log, QhatModel.object = QhatModel.3State.AR3, markov.model.object=markov.3State.US),
model.3StateUS.gamma.log = new('hydroState', input.data=input.data, Qhat.object = Qhat.log, QhatModel.object = QhatModel.3State.gamma , markov.model.object=markov.3State.US),
model.3StateUS.gamma.log.AR1 = new('hydroState',input.data=input.data, Qhat.object = Qhat.log, QhatModel.object = QhatModel.3State.gamma.AR1, markov.model.object=markov.3State.US),
model.3StateUS.gamma.log.AR2 = new('hydroState',input.data=input.data, Qhat.object = Qhat.log, QhatModel.object = QhatModel.3State.gamma.AR2, markov.model.object=markov.3State.US),
model.3StateUS.gamma.log.AR3 = new('hydroState',input.data=input.data, Qhat.object = Qhat.log, QhatModel.object = QhatModel.3State.gamma.AR3, markov.model.object=markov.3State.US),
#
model.3StateUS.BC = new('hydroState', input.data=input.data, Qhat.object = Qhat.boxcox, QhatModel.object = QhatModel.3State , markov.model.object=markov.3State.US),
model.3StateUS.BC.AR1 = new('hydroState',input.data=input.data, Qhat.object = Qhat.boxcox, QhatModel.object = QhatModel.3State.AR1, markov.model.object=markov.3State.US),
model.3StateUS.BC.AR2 = new('hydroState',input.data=input.data, Qhat.object = Qhat.boxcox, QhatModel.object = QhatModel.3State.AR2, markov.model.object=markov.3State.US),
model.3StateUS.BC.AR3 = new('hydroState',input.data=input.data, Qhat.object = Qhat.boxcox, QhatModel.object = QhatModel.3State.AR3, markov.model.object=markov.3State.US),
#
model.3StateUS.gamma.BC = new('hydroState', input.data=input.data, Qhat.object = Qhat.boxcox, QhatModel.object = QhatModel.3State.gamma , markov.model.object=markov.3State.US),
model.3StateUS.gamma.BC.AR1 = new('hydroState',input.data=input.data, Qhat.object = Qhat.boxcox, QhatModel.object = QhatModel.3State.gamma.AR1, markov.model.object=markov.3State.US),
model.3StateUS.gamma.BC.AR2 = new('hydroState',input.data=input.data, Qhat.object = Qhat.boxcox, QhatModel.object = QhatModel.3State.gamma.AR2, markov.model.object=markov.3State.US),
model.3StateUS.gamma.BC.AR3 = new('hydroState',input.data=input.data, Qhat.object = Qhat.boxcox, QhatModel.object = QhatModel.3State.gamma.AR3, markov.model.object=markov.3State.US)
)
# Define Refernce models. That is, the calibration obecjive function that this models needs to meet or exceed.
.Object@calib.reference.model.name <- list(
model.1State.log = '',
model.1State.log.AR1 = 'model.1State.log',
model.1State.log.AR2 = 'model.1State.log.AR1',
model.1State.log.AR3 = 'model.1State.log.AR2',
model.1State.gamma.log = 'model.1State.log',
model.1State.gamma.log.AR1 = 'model.1State.gamma.log',
model.1State.gamma.log.AR2 = 'model.1State.gamma.log.AR1',
model.1State.gamma.log.AR3 = 'model.1State.gamma.log.AR2',
model.1State.BC = 'model.1State.log',
model.1State.BC = '',
model.1State.BC.AR1 = 'model.1State.BC',
model.1State.BC.AR2 = 'model.1State.BC.AR1',
model.1State.BC.AR3 = 'model.1State.BC.AR2',
model.1State.gamma.BC = 'model.1State.BC',
model.1State.gamma.BC.AR1 = 'model.1State.gamma.BC',
model.1State.gamma.BC.AR2 = 'model.1State.gamma.BC.AR1',
model.1State.gamma.BC.AR3 = 'model.1State.gamma.BC.AR2',
model.2State.log = 'model.1State.log',
model.2State.log.AR1 = 'model.2State.log',
model.2State.log.AR2 = 'model.2State.log.AR1',
model.2State.log.AR3 = 'model.2State.log.AR2',
model.2State.gamma.log = 'model.2State.log',
model.2State.gamma.log.AR1 = 'model.2State.gamma.log',
model.2State.gamma.log.AR2 = 'model.2State.gamma.log.AR1',
model.2State.gamma.log.AR3 = 'model.2State.gamma.log.AR2',
model.2State.BC = 'model.2State.log',
model.2State.BC = 'model.1State.BC',
model.2State.BC.AR1 = 'model.2State.BC',
model.2State.BC.AR2 = 'model.2State.BC.AR1',
model.2State.BC.AR3 = 'model.2State.BC.AR2',
model.2State.gamma.BC = 'model.1State.gamma.BC',
model.2State.gamma.BC = 'model.2State.BC',
model.2State.gamma.BC.AR1 = 'model.2State.gamma.BC',
model.2State.gamma.BC.AR2 = 'model.2State.gamma.BC.AR1',
model.2State.gamma.BC.AR3 = 'model.2State.gamma.BC.AR2',
model.3State.log = 'model.2State.log',
model.3State.log.AR1 = 'model.3State.log',
model.3State.log.AR2 = 'model.3State.log.AR1',
model.3State.log.AR3 = 'model.3State.log.AR2',
model.3State.gamma.log = 'model.3State.log',
model.3State.gamma.log.AR1 = 'model.3State.gamma.log',
model.3State.gamma.log.AR2 = 'model.3State.gamma.log.AR1',
model.3State.gamma.log.AR3 = 'model.3State.gamma.log.AR2',
model.3State.BC = 'model.3State.log',
model.3State.BC = 'model.2State.BC',
model.3State.BC.AR1 = 'model.3State.BC',
model.3State.BC.AR2 = 'model.3State.BC.AR1',
model.3State.BC.AR3 = 'model.3State.BC.AR2',
model.3State.gamma.BC = 'model.3State.BC',
model.3State.gamma.BC.AR1 = 'model.3State.gamma.BC',
model.3State.gamma.BC.AR2 = 'model.3State.gamma.BC.AR1',
model.3State.gamma.BC.AR3 = 'model.3State.gamma.BC.AR2',
model.3StateUS.log = 'model.3State.log',
model.3StateUS.log.AR1 = 'model.3StateUS.log',
model.3StateUS.log.AR2 = 'model.3StateUS.log.AR1',
model.2StateUS.log.AR3 = 'model.3StateUS.log.AR2',
model.3StateUS.gamma.log = 'model.3StateUS.log',
model.3StateUS.gamma.log.AR1 = 'model.3StateUS.gamma.log',
model.3StateUS.gamma.log.AR2 = 'model.3StateUS.gamma.log.AR1',
model.3StateUS.gamma.log.AR3 = 'model.3StateUS.gamma.log.AR2',
model.3StateUS.BC = 'model.3State.BC',
model.3StateUS.BC.AR1 = 'model.3StateUS.BC',
model.3StateUS.BC.AR2 = 'model.3StateUS.BC.AR1',
model.3StateUS.BC.AR3 = 'model.3StateUS.BC.AR2',
model.3StateUS.gamma.BC = 'model.3StateUS.BC',
model.3StateUS.gamma.BC.AR1 = 'model.3StateUS.gamma.BC',
model.3StateUS.gamma.BC.AR2 = 'model.3StateUS.gamma.BC.AR1',
model.3StateUS.gamma.BC.AR3 = 'model.3StateUS.gamma.BC.AR2'
)
model.names = names(.Object@models)
.Object@calib.reference.criteria.met = rep(F,length(model.names))
names(.Object@calib.reference.criteria.met) <- model.names
return(.Object)
}
)
setMethod(f = "fit",signature="hydroState.allModels",definition=function(.Object,
pop.size.perParameter=25,
max.generations=10000,
Domains,
reltol=1e-8,
steptol=50,
print.iterations = 25,
doParallel=F,
...)
{
# Set the min and max numbe rof calibrations per model
minTrialsPerModel=5
maxTrialsPerModel=20
model.names <- names(.Object@models)
nModels = length(.Object@models)
message(paste('Calibrating ',nModels,'Models.'))
message(paste('... Minimum number of calibration per model: ',minTrialsPerModel))
message(paste('... Maximum number of calibration per model: ',maxTrialsPerModel))
for (i in 1:nModels) {
# Get reference model calibration result.
current.model.name <- model.names[i]
ref.model.name <- .Object@calib.reference.model.name[[current.model.name]]
if (nchar(ref.model.name)==0){
max.calib.result=-Inf
} else {
ref.model.obj <- .Object@models[[ref.model.name]]
max.calib.result <- ref.model.obj@calibration.results$bestval
}
# Print user summary
message('************************************')
message(paste('Calibrating Model:',current.model.name))
message(paste('Reference model name for maximum acceptable calib. solution:',ref.model.name))
message(paste('Reference model maximum acceptable calib. solution:',max.calib.result))
message('')
current.model.err = Inf
j=1
while (j<minTrialsPerModel || (j<maxTrialsPerModel && current.model.err>max.calib.result)){
DEstrategy <- sample(c(2,3),1,replace=T)
message('')
message(paste('... Doing Calibration Trial ',j, ' using DEstrategy:',DEstrategy))
# Get number of params
nParams = length(getParameters.asVector(.Object@models[[i]]))
if (nParams>=doParallel || doParallel==T) {
model <- .Object@models[[i]]
assign("model", model, envir=globalenv())
model <- fit(model, DEstrategy=DEstrategy, pop.size.perParameter=pop.size.perParameter, max.generations=max.generations, reltol=reltol, steptol=steptol, print.iterations=print.iterations,
parallelType=1, packages = c('hydroState','truncnorm'),parVar=c('model'))
} else {
model <- fit(.Object@models[[i]], DEstrategy=DEstrategy, pop.size.perParameter=pop.size.perParameter, max.generations=max.generations, reltol=reltol, steptol=steptol, print.iterations=print.iterations)
}
# Store the best model to date int the object in case the reference maximum obj cannot be met.
if (model@calibration.results$bestval < current.model.err) {
.Object@models[[i]] <- model
current.model.err <- model@calibration.results$bestval
}
j <- j+1
}
if (current.model.err<=max.calib.result) {
.Object@calib.reference.criteria.met[[current.model.name]] = T
message('Best solution was <= reference model calibration solution!')
} else {
warning('Best solution was > reference model calibration solution!')
}
message('************************************')
}
return(.Object)
}
)
setMethod(f="getNegLogLikelihood",signature="hydroState.allModels",definition=function(.Object)
{
return(getNegLogLikelihood(.Object,plot.data=T))
}
)
setMethod(f="getNegLogLikelihood",signature="hydroState.allModels",definition=function(.Object, plot.data)
{
# Initialise outputs
model.names <- names(.Object@models)
nModels <- length(model.names)
nLL <- rep(Inf, nModels)
names(nLL) <- model.names
# Get AIC for each model
for (i in 1:nModels) {
nLL[i] <- getNegLogLikelihood(.Object@models[[i]])
}
# Plot AIC vs model name
if (missing(plot.data) || plot.data) {
par(mar=c(15,2,2,2))
plot(1:nModels, nLL, type='p',ylab = 'Neg. Log Liklihood')
axis(1, at=1:nModels, labels=model.names,las=2)
}
return(nLL)
}
)
# #
# # Get vector of AICc.
# setMethod(f="getAICc",signature="hydroState.allModels",definition=function(.Object)
# {
# return(getAICc(.Object,plot.data=T))
# }
# )
# # Get vector of AIC.
# setMethod(f="getAIC",signature="hydroState.allModels",definition=function(.Object)
# {
# return(getAIC(.Object,plot.data=T))
# }
# )
setMethod(f="getAIC",signature="hydroState.allModels",definition=function(.Object, plot.data, use.sampleSize.penality)
{
# Initialise outputs
model.names <- names(.Object@models)
nModels <- length(model.names)
AIC <- rep(Inf, nModels)
names(AIC) <- model.names
# Get AIC for each model
if (use.sampleSize.penality) {
for (i in 1:nModels) {
AIC[i] <- getAICc(.Object@models[[i]])
}
} else {
for (i in 1:nModels) {
AIC[i] <- getAIC(.Object@models[[i]])
}
}
# Build vector of model parameters.
nParameters <- rep(Inf, nModels)
names(nParameters) <- model.names
for (i in 1:nModels)
nParameters[i] <- length(getParameters.asVector(.Object@models[[i]]))
# Building a matrix of AIC and nParams.
AIC <- cbind(AIC, nParameters)
rownames(AIC) <- model.names
colnames(AIC) <- c('AIC','nParameters')
# Plot AIC vs model name
if (plot.data) {
par(mar=c(15,2,2,2))
plot(1:nModels, AIC[,1], type='p',ylab = 'AIC')
axis(1, at=1:nModels, labels=model.names,las=2)
}
return(AIC)
}
)
# Get best model (by mimimum AIC)
#' @exportMethod getAIC.bestModel
setGeneric(name="getAIC.bestModel",def=function(.Object, use.calib.reference.criteria=NA, use.sampleSize.penality=NA,min.obs.per.state=NA) {standardGeneric("getAIC.bestModel")})
setMethod(f="getAIC.bestModel",signature="hydroState.allModels",definition=function(.Object, use.calib.reference.criteria=T, use.sampleSize.penality=T, min.obs.per.state=3)
{
# Get AIC
model.names <- names(.Object@models)
nModels <- length(model.names)
AIC <- getAIC(.Object, plot.data=F, use.sampleSize.penality)
# Build filt to models that met calib. criteria
filt = rep(T, nModels)
if (use.calib.reference.criteria) {
for (i in 1:nModels) {
current.model.name <- model.names[i]
filt[i] <- F
if (.Object@calib.reference.criteria.met[[current.model.name]])
filt[i] <- T
}
}
AIC <- AIC[filt,]
# Filter out models that have a non-finite AIC.
nModels <- nrow(AIC)
filt = rep(F, nModels)
for (i in 1:nModels) {
if(is.finite(AIC[i,1]))
filt[i]=T
}
AIC <- AIC[filt,]
# Find those models where (by Viterbi decoding) it is on each state >= min.obs.per.state times.
nModels <- nrow(AIC)
model.names <- rownames(AIC)
filt = rep(T, nModels)
for (i in 1:nModels) {
current.model.name <- model.names[i]
# get Viterbi states.
viterbi.data <- viterbi(.Object@models[[current.model.name]], do.plot=F)
# Count frequency of each state.
if (any(table(viterbi.data[,'Viterbi State Number'])<min.obs.per.state)) {
filt[i] <- F
message(paste('... The following model has a state with only',min(table(viterbi.data[,'Viterbi State Number'])), ' time points: ',current.model.name ))
}
}
AIC <- AIC[filt,]
# Find best model by the lowest AIC.
# If there are >1 models with the best AIC, then take that with the lowest number of parameters.
AIC <- AIC[order(AIC[,1],AIC[,2]),]
bestModel.name <- rownames(AIC)[1]
# Get best 1 state model
ind = grep("model.1State",rownames(AIC))
oneStateModel.names = rownames(AIC)[ind]
ind = which(AIC[,1] == min(AIC[ind,1]))[1]
bestModel.1State.name <- rownames(AIC)[ind]
# Caclulate the delta AIC weights. From, Kenneth P. Burnham, David R. Anderson
# Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach
# Second Edition, section 2.9.
deltaAIC = AIC[,1] - AIC[1,1]
AIC.weights = exp(-deltaAIC/2)/sum(exp(-deltaAIC/2))
AIC = cbind(AIC,AIC.weights)
colnames(AIC) <- c('AIC','nParameters','AIC.weights')
# Calculate the Evidence Ratio between the best model and the
# best one-state model. The approach is from Burnham et al p78.
# To calculate this, first the best one-state model is found.
model.names <- rownames(AIC);
ind = grep('model.1State', model.names, value=F)
ind = ind[1]
if (ind>1) {
EF.oneState = AIC[1,3]/AIC[ind,3]
} else {
EF.oneState = 1;
}
# Get the best model
bestModel <- .Object@models[[bestModel.name]]
# Get the best 1 state model
bestModel.1State <- .Object@models[[bestModel.1State.name]]
return(list(model=bestModel, model.best1State = bestModel.1State, AIC=AIC, evidenceRatio.oneState=EF.oneState))
}
)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.