QMatrix.R

##	
## Copyright 2009 Botond Sipos	
## See the package description for licensing information.	
##	
##########################################################################/** 
#
# @RdocClass QMatrix
# 
# @title "The QMatrix class"
# 
# \description{ 
#	This is the class representing rate matrices. The QMatrix objects store two rate matrices: one contains
#       the rates provided by the user (unscaled rate matrix), the other matrix (scaled rate matrix) is rescaled so the 
#       expected number of subsitutions per unit time is equal to one when the process is at equlibrium.
#	The scaled rate matrices, along with the site-process-specific rate multiplier
#	parameters define the rates of the Event objects generated by the associated GeneralSubstitution objects.          	
#
#	The QMatrix objects have a bidirectonal relationship with the GeneralSubstitution object as they store a reference 
#	to the associated process objects. QMatrix objects do not have a write-protection field, but they use the one
#	from the associated GeneralSubstitution object.
#
#	Usually there is no need to interact with the QMatrix objects directly, so this class
#	is mainly used internally.	
#	
#	@classhierarchy
# }
#	
# @synopsis
#	
# \arguments{
# 	\item{name}{The name of the QMatrix object (a character vector of length one).}
#	\item{alphabet}{An Alphabet object.}
#	\item{rate.list}{A list of unscaled event rates and the associated event names. It will be passed to the \code{setRateList.QMatrix} method.}
#	\item{process}{A GeneralSubstitution object to be associated with the QMatrix object.}
# 	\item{...}{Not used.}
#	}
# 
# \section{Fields and Methods}{ 
# 	@allmethods
# }
# 
# \examples{ 
#	# create a QMatrix object by providing an Alphabet object and rates
#	m<-QMatrix(name="Susie Q", alphabet=BinaryAlphabet(), rate.list=list("1->0"=2,"0->1"=3))
#	# get object summary
#	summary(m)
#       # create a GeneralSubstitution object by
#       # providing an Alphabet object and the rates
#       p<-GeneralSubstitution(alphabet=BinaryAlphabet(), rate.list=list("1->0"=1,"0->1"=3))
#       # get the QMatrix object from p
#       m<-p$QMatrix
#	# get various object properties
#	m
#	is.QMatrix(m)
#	m$name
#	m$id
#	m$alphabet
#	# get the associated process
#	m$process
#	# get the unscaled rate of "0->1"
#	getRate(m,"0->1")
#	# get the scaled rate of "0->1"
#	getEventRate(m,"0->1")
#	# get the list of unscaled rates	
#	m$rateList
#	# get unscaled rate matrix
#	m$matrix
#	# get scaled rate matrix
#	m$scaledMatrix
# }
# 
# @author
#
# \seealso{ 
# 	GeneralSubstitution Alphabet Process
# }
# 
#*/###########################################################################
setConstructorS3(
  "QMatrix",
  function( 
		name="Anonymous",
		alphabet=NA,
		rate.list=NA,
		process=NA,
		... 
		)	{

		this<-PSRoot();
   		this<-extend(
     			 this,
      			"QMatrix",
			.name=NA,
			.alphabet=NA,
			.rate.matrix=NA,
			.orig.matrix=NA,
		  	.norm.const=NA,
			.process=NA,
			.is.q.matrix=TRUE
    		);
		
		this$name<-name;

		if(!missing(process)){
			this$process<-process;
		}

		if(!missing(alphabet)){
			this$alphabet<-alphabet;
		}

		if(!missing(rate.list)){
			if(missing(alphabet)){
				throw("Cannot set rates because the alphabet is not specified!\n")
			}
			this$rateList<-rate.list;
		}

    return(this);
  },
  enforceRCC=TRUE
);

##	
## Method: .buildRateMatrix
##	
setMethodS3(
	".buildRateMatrix", 
	class="QMatrix", 
	function(
		this,
		...
	){

		size<-this$alphabet$size;
		symbols<-this$alphabet$symbols;

		# Setting the dimension names
		# for the original rates matrix:
		if(!isEmpty(this$.alphabet)){
			this$.orig.matrix<-matrix(nrow=size,ncol=size);
			colnames(this$.orig.matrix)<-symbols;
			rownames(this$.orig.matrix)<-symbols;
		}
		# Copy to the scaled matrix:
		this$.rate.matrix<-this$.orig.matrix;
	},
	private=TRUE,
	protected=FALSE,
	overwrite=FALSE,
	conflict="warning",
	validators=getOption("R.methodsS3:validators:setMethodS3")
);

##	
## Method: checkConsistency
##	
###########################################################################/**
#
# @RdocMethod	checkConsistency
# 
# @title "Check object consistency"
# 
# \description{ 
#		@get "title".
# } 
# 
# @synopsis 
#
# \arguments{ 
#       \item{this}{An object.} 
#       \item{check.process}{Check the associated process object.} 
#       \item{...}{Not used.} 
# } 
# 
# 
# \value{ 
#		Returns an invisible TRUE if no inconsistencies found in the object, throws 
#		an error otherwise. 
# } 
# 
# @author 
# 
# \seealso{ 
# 	@seeclass 
# } 
# 
#*/###########################################################################
setMethodS3(
	"checkConsistency", 
	class="QMatrix", 
	function(
		this,
		check.process=TRUE,
		...
	){

			if(is.Process(this$.process)){

      		wp<-this$.process$.write.protected;
      		if (wp) {
			# Cannot use virtual field here because of 
			# some obscure errror:
        		this$.process$.write.protected<-FALSE;
      		}
			}
      may.fail<-function(this) {
					
				# Reassing name:
				this$name<-this$name;
				# Do not reassign alphabet as that will wipe out the rates!
				if(!is.na(this$.alphabet) & !is.Alphabet(this$.alphabet)){
					throw("QMatrix alphabet is invalid!\n");
				}
		
				# Check if the original rate matrix is a matrix or NA.
				if(!is.matrix(this$.orig.matrix) & !all(is.na(this$.orig.matrix))){
						throw("The original rates matrix is invalid!\n");
				}	
				# Check original rates matrix size:
				else if(!all(dim(this$.orig.matix) != c(this$.alphabet$.size, this$.alphabet$.size))) {
					throw("The original rates matrix is of wrong size!");
				}
				# Check if the rescaled rate matrix is a matrix or NA.
				if(!is.matrix(this$.rate.matrix) & !all(is.na(this$.rate.matrix))){
						throw("The rescaled rates matrix is invalid!\n");
				}	

				# Check rescaled rates matrix size:
				else if(!all(dim(this$.rates.matix) != c(this$.alphabet$.size, this$.alphabet$.size))) {
					throw("The original rates matrix is of wrong size!");
				}

				# Flag for not having NA-as in the matrices.
				COMPLETE<-TRUE;

			
				if(is.matrix(this$.orig.matrix) ){
					if(any(is.na(this$.orig.matrix))){
						warning("Some rates are undefined!\n");
						COMPLETE<-FALSE;
					}
					# Watch out for the operator precedence here!
					else if ( (!all(is.numeric(this$.orig.matrix))) & (!all(is.na(this$.orig.matrix))) ){
						print(this$.orig.matrix)
						throw("The original rates matrix has non-numeric elements!\n");
					}
				}
				
				if(is.matrix(this$.orig.matrix) ){
					if( any(is.na(this$.orig.matrix)) & COMPLETE ){
						COMPLETE<-FALSE;
						throw("The original rates matrix is complete, but the rescaled matrix has undefined elements!\n");
					}
					# Watch out for the operator precedence here!
					else if ( (!all(is.numeric(this$.orig.matrix)) & (!all(is.na(this$.orig.matrix)))) ){
						throw("The original rates matrix has non-numeric elements!\n");
					}
				}

				# Check the normalizing constant:
				if( length(this$.norm.const) != 1 | (!is.na(this$.norm.const) & !is.numeric(this$.norm.const)) ){
					 throw("Normalizing constant is invalid!\n");
				}

				# Check the normalization:
				if(is.matrix(this$.orig.matrix) & is.matrix(this$.rate.matrix) & COMPLETE ){
						if(!PSRoot$my.all.equal(this$.rate.matrix, (this$.norm.const * this$.orig.matrix)) ){
							throw("The scaled matrix is inconsistent with the original matrix and the normalizing constant!\n");
						}
				}

				# Check the process:
				if(check.process==TRUE){
					 if(is.Process(this$.process)){
							# Check for alphabet compatibility:
							if(this$.alphabet != this$.process$alphabet){
								throw("Process/QMatrix alphabet mismatch!\n");
							}
							# Check if the parent process QMatrix is this object:
							if(!equals(this$.process$.q.matrix, this) ){
								throw("Parent process QMatrix is not identical with self!\n");
							}
					 }
					 else if(!is.na(this$.process)){
							throw("Process entry is invalid!\n");
					 }
				}
				
      
			}
      tryCatch(may.fail(this),finally={if(is.Process(this$.process)){this$.process$.write.protected<-wp}});
			return(invisible(TRUE));

	},
	private=FALSE,
	protected=FALSE,
	overwrite=FALSE,
	conflict="warning",
	validators=getOption("R.methodsS3:validators:setMethodS3")
);

##	
## Method: .callRateRescaling
##	
setMethodS3(
	".callRateRescaling", 
	class="QMatrix", 
	function(
		this,
		guess.equ=TRUE,
		...
	){

		# Usually called when the rate matrix chenges.
		# If the Q matrix has a parent process with a valid equilibrium distribution:
		if(is.Process(this$.process) & !any(is.na(as.vector(this$.orig.matrix))) ){

		if(guess.equ){
			# Try to guess the new equlibrium distribution:
			if(!.setEquDistFromGuess(this$.process)){
						# Fill with NA-s if failed with guessing:
      	    .initEquDist(this$.process);
      	}
		}
			# Rescale if the equilibrium distribution was succesfully set:
			if(all(!is.na(this$.process$equDist))){
				rescaleQMatrix(this$.process);
			}

		}

	},
	private=TRUE,
	protected=FALSE,
	overwrite=FALSE,
	conflict="warning",
	validators=getOption("R.methodsS3:validators:setMethodS3")
);

##	
## Method: getName
##	
###########################################################################/**
#
# @RdocMethod getName
# 
# @title "Get the name of a QMatrix object" 
# 
# \description{ 
#	@get "title".
# } 
# 
# @synopsis 
# 
# \arguments{ 
# 	\item{this}{A QMatrix object.} 
# 	\item{...}{Not used.} 
# } 
# 
# \value{ 
# 	A charcter vector of length one.
# } 
# 
# \examples{
#	# create a QMatrix object
#	m<-QMatrix()
#	# set/get name
#	setName(m,"Susie Q")
#	getName(m)
#	# set/get name via virtual field
#	m$name<-"Q"
#	m$name
# } 
# 
# @author 
# 
# \seealso{ 
# 	@seeclass 
# } 
# 
#*/###########################################################################
setMethodS3(
	"getName", 
	class="QMatrix", 
	function(
		this,
		...
	){

		this$.name;

	},
	private=FALSE,
	protected=FALSE,
	overwrite=FALSE,
	conflict="warning",
	validators=getOption("R.methodsS3:validators:setMethodS3")
);

##	
## Method: setName
##	
###########################################################################/**
#
# @RdocMethod setName
# 
# @title "Set the name of a QMatrix object" 
# 
# \description{ 
#	@get "title".
# } 
# 
# @synopsis 
# 
# \arguments{ 
# 	\item{this}{A QMatrix object.} 
#	\item{value}{A character vector of length one.}
# 	\item{...}{Not used.} 
# } 
# 
# \value{ 
# 	The new object name.
# } 
# 
# \examples{
#	# create a QMatrix object
#	m<-QMatrix()
#	# set/get name
#	setName(m,"Susie Q")
#	getName(m)
#	# set/get name via virtual field
#	m$name<-"Q"
#	m$name
# } 
# 
# @author 
# 
# \seealso{ 
# 	@seeclass 
# } 
# 
#*/###########################################################################
setMethodS3(
	"setName", 
	class="QMatrix", 
	function(
		this,
		value,
		...
	){

	.checkWriteProtection(this);
	if(missing(value)){
		throw("No new value provided!\n");
	}
	else{
		value<-as.character(value);
		if(stringLength(value) == 0){
			throw("Cannot set empty name!");
		} else {
			this$.name<-value;
		}
	}
	return(this$.name);

	},
	private=FALSE,
	protected=FALSE,
	overwrite=FALSE,
	conflict="warning",
	validators=getOption("R.methodsS3:validators:setMethodS3")
);

##	
## Method: getProcess
##	
###########################################################################/**
#
# @RdocMethod getProcess
# 
# @title "Get the process object associated with a QMatrix object" 
# 
# \description{ 
#	@get "title".
# } 
# 
# @synopsis 
# 
# \arguments{ 
# 	\item{this}{A QMatrix object.} 
# 	\item{...}{Not used.} 
# } 
# 
# \value{ 
# 	A process object, most likely one which inherits from GeneralSubstitution.
# } 
# 
# \examples{
#	# Create a GeneralSubstitution object
#	p<-GeneralSubstitution(alphabet=BinaryAlphabet())
#	p
#	# get the associated QMatrix object from p
#	m<-p$qMatrix
#	summary(m)
#	# get the associated process from m
#	m$process
#	# clone p
#	pp<-clone(p)
#	# assotiate m with pp
#	pp$qMatrix<-m
#	# assotiate pp with m
#	m$process<-pp
#	m$process
# } 
# 
# @author 
# 
# \seealso{ 
# 	@seeclass 
# } 
# 
#*/###########################################################################
setMethodS3(
	"getProcess", 
	class="QMatrix", 
	function(
		this,
		...
	){

		this$.process;

	},
	private=FALSE,
	protected=FALSE,
	overwrite=FALSE,
	conflict="warning",
	validators=getOption("R.methodsS3:validators:setMethodS3")
);

##	
## Method: setProcess
##	
###########################################################################/**
#
# @RdocMethod setProcess
# 
# @title "Assotiate a process object with a QMatrix object" 
# 
# \description{ 
#	@get "title".
# } 
# 
# @synopsis 
# 
# \arguments{ 
# 	\item{this}{A QMatrix object.} 
#	\item{value}{A Process object.}
# 	\item{...}{Not used.} 
# } 
# 
# \value{ 
# 	The Process object.
# } 
# 
# \examples{
#	# Create a GeneralSubstitution object
#	p<-GeneralSubstitution(alphabet=BinaryAlphabet())
#	p
#	# get the associated QMatrix object from p
#	m<-p$qMatrix
#	summary(m)
#	# get the associated process from m
#	m$process
#	# clone p
#	pp<-clone(p)
#	# assotiate m with pp
#	pp$qMatrix<-m
#	# assotiate pp with m
#	m$process<-pp
#	m$process
# } 
# 
# @author 
# 
# \seealso{ 
# 	@seeclass 
# } 
# 
#*/###########################################################################
setMethodS3(
	"setProcess", 
	class="QMatrix", 
	function(
		this,
		value,
		...
	){

		.checkWriteProtection(this);
		if(missing(value)){
			throw("No new value provided!\n");
		}
		else if (!is.Process(value)){
			throw("Process object invalid!\n");
		} 
		else {
			this$.process<-value;
		}

	},
	private=FALSE,
	protected=FALSE,
	overwrite=FALSE,
	conflict="warning",
	validators=getOption("R.methodsS3:validators:setMethodS3")
);

##	
## Method: getAlphabet
##	
###########################################################################/**
#
# @RdocMethod getAlphabet
# 
# @title "Get the Alphabet object associated with a QMatrix object" 
# 
# \description{ 
#	@get "title".
# } 
# 
# @synopsis 
# 
# \arguments{ 
# 	\item{this}{A QMatrix object.} 
# 	\item{...}{Not used.} 
# } 
# 
# \value{ 
# 	An Alphabet object.
# } 
# 
# \examples{
#	# create a QMatrix object
#	m<-QMatrix()
#	# set the alphabet
#	setAlphabet(m,NucleotideAlphabet())
#	# get the alphabet
#	getAlphabet(m)
#	# set alphabet via virtual field
#	m$alphabet<-BinaryAlphabet()
#	summary(m)
# } 
# 
# @author 
# 
# \seealso{ 
# 	@seeclass 
# } 
# 
#*/###########################################################################
setMethodS3(
	"getAlphabet", 
	class="QMatrix", 
	function(
		this,
		...
	){

		this$.alphabet;

	},
	private=FALSE,
	protected=FALSE,
	overwrite=FALSE,
	conflict="warning",
	validators=getOption("R.methodsS3:validators:setMethodS3")
);

##	
## Method: setAlphabet
##	
###########################################################################/**
#
# @RdocMethod setAlphabet
# 
# @title "Set the Alphabet object for a QMatrix object" 
# 
# \description{ 
#	@get "title".
#
#	This method rebuilds the scaled and unscaled rate matrices and so sets all rates to NA.
# } 
# 
# @synopsis 
# 
# \arguments{ 
# 	\item{this}{A QMatrix object.} 
#	\item{value}{An Alphabet object.}
# 	\item{...}{Not used.} 
# } 
# 
# \value{ 
# 	The Alphabet object.
# } 
# 
# \examples{
#	# create a QMatrix object
#	m<-QMatrix()
#	# set the alphabet
#	setAlphabet(m,NucleotideAlphabet())
#	# get the alphabet
#	getAlphabet(m)
#	# set alphabet via virtual field
#	m$alphabet<-BinaryAlphabet()
#	summary(m)
# } 
# 
# @author 
# 
# \seealso{ 
# 	@seeclass 
# } 
# 
#*/###########################################################################
setMethodS3(
	"setAlphabet", 
	class="QMatrix", 
	function(
		this,
		value,
		...
	){

	.checkWriteProtection(this);

	if(!exists(x="PSIM_FAST")){

	if(missing(value)){
		throw("No new value provided!\n");
	}
	else if(!is.Alphabet(value)) {
		throw("Alphabet object invalid!\n");
	}
	else if(is.Process(this$.process)){
		if(value != this$.process$alphabet){
		throw("The new alphabet should match with the one from the subsitution process!\n");
		}
	}
	if(is.matrix(this$.rate.matrix)){
			warning("Be aware that setting a new alphabet wipes out completely the rate matrix!\n");
	}
	
	}
	
	this$.alphabet<-value;
	.buildRateMatrix(this);
	return(this$.alphabet);

	},
	private=FALSE,
	protected=FALSE,
	overwrite=FALSE,
	conflict="warning",
	validators=getOption("R.methodsS3:validators:setMethodS3")
);

##	
## Method: .nameToDim
##	
setMethodS3(
	".nameToDim", 
	class="QMatrix", 
	function(
		this,
		name,
		...
	){

		# split the name
		substitution<-rbind(strsplit(name,split="->",fixed=TRUE)[[1]]);
		if(length(substitution) != 2 ) {
			throw("Substitution event name was invalid!");
		}

		# Check if symbols are valid:
		if(!hasSymbols(this$.alphabet, substitution)){
				throw("All symbols must be in the alphabet!\n");
		}

		# Return a vector with the DIMNAMES:
		colnames(substitution)<-c("from","to");
		rownames(substitution)<-c("Substitution");
		substitution;

	},
	private=FALSE,
	protected=FALSE,
	overwrite=FALSE,
	conflict="warning",
	validators=getOption("R.methodsS3:validators:setMethodS3")
);

##	
## Method: .nameToDimFast
##	
setMethodS3(
	".nameToDimFast", 
	class="QMatrix", 
	function(
		this,
		name,
		...
	){

		substitution<-rbind(strsplit(name,split="->",fixed=TRUE)[[1]]);
		colnames(substitution)<-c("from","to");
		rownames(substitution)<-c("Substitution");
		substitution;

	},
	private=FALSE,
	protected=FALSE,
	overwrite=FALSE,
	conflict="warning",
	validators=getOption("R.methodsS3:validators:setMethodS3")
);

##	
## Method: .dimToName
##	
setMethodS3(
	".dimToName", 
	class="QMatrix", 
	function(
		this,
		dim,
		...
	){
		paste(dim,collapse="->");
	},
	private=FALSE,
	protected=FALSE,
	overwrite=FALSE,
	conflict="warning",
	validators=getOption("R.methodsS3:validators:setMethodS3")
);

##	
## Method: getEventRate
##	
###########################################################################/**
#
# @RdocMethod getEventRate
# 
# @title "Get the unscaled rate of an event from a QMatrix object" 
# 
# \description{ 
#       @get "title".
#       
#       This method return the element corresponding to a given event from the scaled rate matrix stored in a QMatrix object.
#       The event can be specified by the inital and target states ("from" and "to" arguments), or by the
#       event name ("from->to"). The event name takes precedence over the "from" and "to" arguments. 
#
#	This method returns NA if the resacling of the rates was not performed.
#	The scaling is performed by the \code{rescaleQMatrix.GeneralSubstitution} method.
# } 
# 
# @synopsis 
# 
# \arguments{ 
#       \item{this}{A QMatrix object.} 
#       \item{name}{The name of the event.}
#       \item{from}{The initial state.}
#       \item{to}{Target state.}
#       \item{...}{Not used.} 
# } 
# 
# \value{ 
#       A Numeric vector of length one.
# } 
# 
# \examples{
#       # create a QMatrix object
#       # provide an Alphabet object and the rates
#       m<-QMatrix(alphabet=BinaryAlphabet(), rate.list=list("0->1"=1,"1->0"=1))
#       # get the unscaled rate of "0->1" by name
#       getEventRate(m,"0->1")  # retruns NA
#       # get the unscaled rate of "0->1" by states
#       getEventRate(m,from="0",to="1")	# returns NA
# } 
# 
# @author 
# 
# \seealso{ 
#       @seeclass 
# } 
# 
#*/###########################################################################
setMethodS3(
	"getEventRate", 
	class="QMatrix", 
	function(
		this,
		name=NA,
		from=NA,
		to=NA,
		...
	){

			# Event specified by name:
			if(!missing(name) & missing(from) & missing(to)){
				# convert to dimnames
				tmp<-.nameToDim(this, name);
				# Return the rate from the rescaled matrix: 
				return(this$.rate.matrix[tmp[1],tmp[2]]);
			}
			# Event specified by from= and to=
			else if(missing(name) & !missing(from) & !missing(to)){
				# Check symbols:
				if(!hasSymbols(this$.alphabet, c(from,to))){
					throw("All symbols must be in the alphabet!\n")
				}
				else{
					# Get the rate from the rescaled matrix:
					return(this$.rate.matrix[as.character(from),as.character(to)]);
				}
			}
			else {
				throw("The substitution should be specified by name or by the \"from\" and \"to\" arguments!\n");
			}

	},
	private=FALSE,
	protected=FALSE,
	overwrite=FALSE,
	conflict="warning",
	validators=getOption("R.methodsS3:validators:setMethodS3")
);

##	
## Method: .getEventRateFast
##	
setMethodS3(
	".getEventRateFast", 
	class="QMatrix", 
	function(
		this,
		name=NA,
		...
	){
				# convert to dimnames
				tmp<-.nameToDimFast(this, name);
				# Return the rate from the rescaled matrix: 
				return(this$.rate.matrix[tmp[1],tmp[2]]);
	},
	private=TRUE,
	protected=FALSE,
	overwrite=FALSE,
	conflict="warning",
	validators=getOption("R.methodsS3:validators:setMethodS3")
);

##	
## Method: getRate
##	
###########################################################################/**
#
# @RdocMethod getRate
# 
# @title "Get an unscaled rate of an event from a QMatrix object" 
# 
# \description{ 
#       @get "title".
#       
#       This method gets the element corresponding to a given event form the \emph{unscaled} rate matrix.
#       a given event. The event can be specified by the inital and target states ("from" and "to" arguments), or by the
#       event name ("from->to"). The event name takes precedence over the "from" and "to" arguments. 
#
# } 
# 
# @synopsis 
# 
# \arguments{ 
#       \item{this}{A QMatrix object.} 
#       \item{name}{The name of the event.}
#       \item{from}{The initial state.}
#       \item{to}{Target state.}
#       \item{...}{Not used.} 
# } 
# 
# \value{ 
#       A Numeric vector of length one.
# } 
# 
# \examples{
#       # create a QMatrix object
#       # provide an Alphabet object and the rates
#       m<-QMatrix(alphabet=BinaryAlphabet(), rate.list=list("1->0"=1,"0->1"=1))
#       # get the unscaled rate of "0->1" by name
#       getRate(m,"0->1")       
#       # get the unscaled rate of "0->1" by states
#       getRate(m,from="0",to="1")
# } 
# 
# 
# \seealso{ 
#       @seeclass 
# } 
# 
#*/###########################################################################
setMethodS3(
	"getRate", 
	class="QMatrix", 
	function(
		this,
		name=NA,
		from=NA,
		to=NA,
		...
	){

			# Event specified by name:
			if(!missing(name) & missing(from) & missing(to)){
				# Convert to dimnames:
				tmp<-.nameToDim(this, name);
				# return unscaled rate:
				return(this$.orig.matrix[tmp[1],tmp[2]]);
			}
			# Event specified by from= and to=:
			else if(missing(name) & !missing(from) & !missing(to)){
				# check symbols:
				if(!hasSymbols(this$.alphabet, c(from,to))){
					throw("All symbols must be in the alphabet!\n")
				}
				else{
				# return unscaled rate:
					return(this$.orig.matrix[as.character(from),as.character(to)]);
				}
			}
			else {
				throw("The substitution should be specified by name or by the \"from\" and \"to\" arguments!\n");
			}

	},
	private=FALSE,
	protected=FALSE,
	overwrite=FALSE,
	conflict="warning",
	validators=getOption("R.methodsS3:validators:setMethodS3")
);


##	
## Method: setRate
##	
###########################################################################/**
#
# @RdocMethod setRate
# 
# @title "Set an unscaled rate in a QMatrix object" 
# 
# \description{ 
#       @get "title".
#       
#       This method sets the element corresponding to a given event in the unscaled rate matrix.
#       The event can be specified by the inital and target states ("from" and "to" arguments), or by the
#       event name ("from->to"). The event name takes precedence over the "from" and "to" arguments. 
#
# } 
# 
# @synopsis 
# 
# \arguments{ 
#       \item{this}{A QMatrix object.} 
#       \item{name}{The name of the event.}
#       \item{value}{The new value of the rate.}
#       \item{from}{The initial state.}
#       \item{to}{Target state.}
#       \item{scale}{Call rate rescaling.}
#       \item{diag}{Calculate diagonal elements.}
#       \item{guess.equ}{Guess equilibrium distribution.}
#       \item{...}{Not used.} 
# } 
# 
# \value{ 
#       A Numeric vector of length one.
# } 
# 
# \examples{
#       # create a QMatrix object
#       # provide an Alphabet object and the rates
#       m<-QMatrix(alphabet=BinaryAlphabet(), rate.list=list("1->0"=1,"0->1"=1))
#       # set the unscaled rate by event name
#       setRate(m,"0->1",2)
#       # get the unscaled rate of "0->1" by name
#       getRate(m,"0->1")       
#       # set the unscaled rate by states
#       setRate(m,"0->1",0.5)
#       # set the unscaled rate of "0->1" by states
#       setRate(m,"0->1",0.5)
#       # get the unscaled rate of "0->1" by states
#       getRate(m,from="0",to="1")
# } 
# 
# @author 
# 
# \seealso{ 
#       @seeclass 
# } 
# 
#*/###########################################################################
setMethodS3(
	"setRate", 
	class="QMatrix", 
	function(
		this,
		name=NA,
		value,
		from=NA,
		to=NA,
		scale=TRUE,
		diag=TRUE,
		guess.equ=TRUE,
		...
	){

		.checkWriteProtection(this);
	if(!exists(x="PSIM_FAST")){

		if (isEmpty(this$.alphabet)){
				throw("Alphabet is valid but empty, so no rates are defined!\n");
		}
		else if(missing(value)) {
			throw("No new value provided!\n");}
		else if(!is.numeric(value)) {
			throw("Rate must be numeric!\n");
		}
		else if (value < 0){
			throw("Cannot set negative rate!\n");
		}
	}
			
			.from<-character();		
			.to<-character();		
		
			# Event specified by name:		
			if(!missing(name) & missing(from) & missing(to)){
				# convert to dimnames:
				tmp<-.nameToDim(this, name);
				.from<-tmp[1];
				.to<-tmp[2];
				
			}
			# Event specified by from= and to=:
			else if(missing(name) & !missing(from) & !missing(to)){
				# check the symbols
				if(!hasSymbols(this$.alphabet, c(from,to))){
					throw("All symbols must be in the alphabet!\n")
				}
				else{
				.from<-as.character(from);
				.to<-as.character(to);
				}
			}
			else {
				throw("The substitution should be specified by name or by the \"from\" and \"to\" arguments!\n");
			}
				
				# Complain if tried to modify a diagonal element:
				if(.from == .to){
					throw("Modification of diagonal elements is not allowed!\n");
				}
				else {
				
					# Set the rate in the original rates matrix:	
			 		this$.orig.matrix[.from, .to]<-value;
			 		# Set the new diagonal element in the original rates matrix:
					if (diag == TRUE) {
			 			this$.orig.matrix[.from, .from]<-.calculateDiagonal(this, symbol=.from);
					}
					
					# Call rate rescaling, this will set the new values
					# in the rescaled rates matrix:
			 		if(scale==TRUE){
						.callRateRescaling(this,guess.equ);
			 		}
			}
		 return(invisible(value));


	},
	private=FALSE,
	protected=FALSE,
	overwrite=FALSE,
	conflict="warning",
	validators=getOption("R.methodsS3:validators:setMethodS3")
);

##	
## Method: getRateList
##	
###########################################################################/**
#
# @RdocMethod getRateList
# 
# @title "Get a list of events and their unscaled rates from a QMatrix object" 
# 
# \description{ 
#       @get "title".
#
#       This method returns the list of event rates from the \emph{unscaled} rate matrix. 
#       The returned list contains the rates associated with the corresponding event names.
# } 
# 
# @synopsis 
# 
# \arguments{ 
#       \item{this}{A GeneralSubstitution object.} 
#       \item{...}{Not used.} 
# } 
# 
# \value{ 
#       A list of event rates.
# } 
# 
# \examples{
#       # create a GeneralSubstitution object
#       # provide an Alphabet object and the rates
#       p<-GeneralSubstitution(alphabet=BinaryAlphabet(), rate.list=list("1->0"=1,"0->1"=3))
#	# get the QMatrix object from p
#	m<-p$QMatrix
#       # get the event rates from the unscaled Q matrix
#       getRateList(m)
#       # get rates from the unscaled rate matrix via virtual field
#       m$rateList
#       # set rates in the unscaled rate matrix
#       setRateList(m, list("0->1"=1,"1->0"=1))
#       m$rateList
#       # set rates in the unscaled rate matrix via virtual field
#       m$rateList<-list("0->1"=3,"1->0"=1);
#	m$rateList
# } 
# 
# @author 
# 
# \seealso{ 
#       @seeclass 
# } 
# 
#*/###########################################################################
setMethodS3(
	"getRateList", 
	class="QMatrix", 
	function(
		this,
		...
	){
		# Be gentle and return an empty list if the
		# alphabet is empty:
		if( isEmpty(this$.alphabet) ){
				return(list());	
		}
		else {
				# Fill in the rates list:
				rates<-list();
				for(i in this$.alphabet$symbols){
					for(j in this$.alphabet$symbols){
						if(i != j){
							name<-paste(i,j,sep="->");
							rate<-getRate(this, from=i, to=j);
							rates[[name]]<-rate;
						}
					}
				}
				return(rates);
		}

	},
	private=FALSE,
	protected=FALSE,
	overwrite=FALSE,
	conflict="warning",
	validators=getOption("R.methodsS3:validators:setMethodS3")
);

##	
## Method: setRateList
##	
###########################################################################/**
#
# @RdocMethod	setRateList
# 
# @title "Setting the unscaled rates stored in a QMatrix object" 
# 
# \description{ 
#       @get "title".
#       
#       This method set the rates in the \emph{unscaled} Q  matrix based on the provided list containing even names
#       and the associated rates. The rate must be specified for every event!
# } 
# 
# @synopsis 
# 
# \arguments{ 
#       \item{this}{A QMatrix object.} 
#       \item{value}{A list with the events names and the associated rates.}
#       \item{...}{Not used.} 
# } 
# 
# \value{ 
#       The QMatrix object (invisible).
# } 
# 
# \examples{
#       # create a GeneralSubstitution object
#       # provide an Alphabet object and the rates
#       p<-GeneralSubstitution(alphabet=BinaryAlphabet(), rate.list=list("1->0"=1,"0->1"=3))
#	# get the QMatrix object from p
#	m<-p$QMatrix
#       # get the event rates from the unscaled Q matrix
#       getRateList(m)
#       # get rates from the unscaled rate matrix via virtual field
#       m$rateList
#       # set rates in the unscaled rate matrix
#       setRateList(m, list("0->1"=1,"1->0"=1))
#       m$rateList
#       # set rates in the unscaled rate matrix via virtual field
#       m$rateList<-list("0->1"=3,"1->0"=1);
#	m$rateList
# } 
# 
# @author 
# 
# \seealso{ 
#       @seeclass 
# } 
# 
#*/###########################################################################
setMethodS3(
	"setRateList", 
	class="QMatrix", 
	function(
		this,
		value,
		...
	){

		.checkWriteProtection(this);
		if(!is.Alphabet(this$.alphabet)){
			throw("Cannot set rates because alphabet is undefined!\n");
		}
		if(missing(value)) {
			throw("No new value provided!\n");}
		else if(!is.list(value)) {
			throw("The new value must be a list!\n");
		} else {

			# Check if all of the rates are specified!
			expected<-.genExpected(this);
	
			# Check for missing rates:		
			if(length(tmp<-setdiff(expected,names(value))) != 0 ){
					throw("The rate matrix is not specified correctly, the following rates are missing: ",paste(tmp,collapse=" "),"!\n");
			}

			# Warn the user about the superfluous rates:
			if(length(tmp<-setdiff(names(value),expected)) != 0 ){
					warning("The following rates were not expected by this process: ",paste(tmp,collapse=" "),", so they were ignored!\n");
					# Getting rid of unexpected rates:
					value[tmp]<-NULL;
			}

				# set the rate matrix if all is OK!
				for (name in names(value)){
					setRate(this,name=name,value=(value[[name]]),scale=FALSE,diag=FALSE);
				}

				# Set diagonal elements:
				for (sym in this$alphabet$symbols){
					this$.orig.matrix[sym, sym]<-.calculateDiagonal(this, sym);
				}
	
				# Call the parent process object to guess the new equlibrium distribution and rescale
				# rates:
				.callRateRescaling(this);

		}
		return(invisible(this))


	},
	private=FALSE,
	protected=FALSE,
	overwrite=FALSE,
	conflict="warning",
	validators=getOption("R.methodsS3:validators:setMethodS3")
);

##	
## Method: .genExpected
##	
setMethodS3(
	".genExpected", 
	class="QMatrix", 
	function(
		this,
		...
	){

      expected<-list();
      sym<-this$alphabet$symbols;

      # Cretae the list of expected rates:
      for(i in sym){
        for(j in sym){
          if(i != j){
          expected<-c(expected,paste(i,j,sep="->"));
          }
        }
      }
			return(expected);
		

	},
	private=TRUE,
	protected=FALSE,
	overwrite=FALSE,
	conflict="warning",
	validators=getOption("R.methodsS3:validators:setMethodS3")
);

##	
## Method: .calculateDiagonal
##	
setMethodS3(
	".calculateDiagonal", 
	class="QMatrix", 
	function(
		this,
		symbol=NA,
		...
	){

		if(!missing(symbol)){
			# convert diname to dim:
			index<-.symToIndex(this, symbol=symbol);
		}
		else {
			throw("Symbol not specified!\n");
		}
			# Return -1 * sum of the off-diagonal elements
			# from the row specified by the index:
			return(-sum((this$.orig.matrix[symbol,])[-index] ));
	},
	private=FALSE,
	protected=FALSE,
	overwrite=FALSE,
	conflict="warning",
	validators=getOption("R.methodsS3:validators:setMethodS3")
);

##	
## Method: .symToIndex
##	
setMethodS3(
	".symToIndex", 
	class="QMatrix", 
	function(
		this,
		symbol=NA,
		...
	){

		if(exists(x="PSIM_FAST")){
			return(which(rownames(this$.orig.matrix) == symbol));
		}

		if(missing(symbol)){
			throw("No symbol specified");
		} else {
			index<-which(rownames(this$.orig.matrix) == symbol);
			if(length(index) == 0){
				print(symbol);
				throw("Symbol not in rate matrix!\n");
			}
			else if (length(index) > 1){
				throw("Rate matrix is inconsistent!\n");
			}
			return(index);
		}


	},
	private=FALSE,
	protected=FALSE,
	overwrite=FALSE,
	conflict="warning",
	validators=getOption("R.methodsS3:validators:setMethodS3")
);

##	
## Method: getMatrix
##	
###########################################################################/**
#
# @RdocMethod	getMatrix
# 
# @title "Get the unscaled rate matrix form a QMatrix object" 
# 
# \description{ 
#	@get "title".
# } 
# 
# @synopsis 
# 
# \arguments{ 
# 	\item{this}{A QMatrix object.} 
# 	\item{...}{Not used.} 
# } 
# 
# \value{ 
# 	A matrix object.
# } 
# 
# \examples{
#       # create a GeneralSubstitution object
#       # provide an Alphabet object and the rates
#       p<-GeneralSubstitution(alphabet=BinaryAlphabet(), rate.list=list("1->0"=1,"0->1"=3))
#	# get the QMatrix object from p
#	m<-p$QMatrix
#	# get the unscaled rate matrix from m
#	m$matrix
# } 
# 
# @author 
# 
# \seealso{ 
# 	@seeclass 
# } 
# 
#*/###########################################################################
setMethodS3(
	"getMatrix", 
	class="QMatrix", 
	function(
		this,
		...
	){

		this$.orig.matrix;

	},
	private=FALSE,
	protected=FALSE,
	overwrite=FALSE,
	conflict="warning",
	validators=getOption("R.methodsS3:validators:setMethodS3")
);

##	
## Method: setMatrix
##	
###########################################################################/**
#
# @RdocMethod setMatrix
#
# @title "Forbidden action: setting the unscaled rate matrix for a QMatrix object"
#
# \description{
#       @get "title".
# }
#
# @synopsis
#
# \arguments{
#       \item{this}{An object.}
#       \item{value}{Not used.}
#       \item{...}{Not used.}
# }
#
# \value{
#	Throws an error.
# }
#
# @author
#
# \seealso{
#       @seeclass
# }
#
#*/###########################################################################
setMethodS3(
	"setMatrix", 
	class="QMatrix", 
	function(
		this,
		value,
		...
	){

		virtualAssignmentForbidden(this);

	},
	private=FALSE,
	protected=FALSE,
	overwrite=FALSE,
	conflict="warning",
	validators=getOption("R.methodsS3:validators:setMethodS3")
);

##	
## Method: getScaledMatrix
##	
###########################################################################/**
#
# @RdocMethod	getScaledMatrix
# 
# @title "Get the scaled rate matrix form a QMatrix object" 
# 
# \description{ 
#	@get "title".
# } 
# 
# @synopsis 
# 
# \arguments{ 
# 	\item{this}{A QMatrix object.} 
# 	\item{...}{Not used.} 
# } 
# 
# \value{ 
# 	A matrix object.
# } 
# 
# \examples{
#       # create a GeneralSubstitution object
#       # provide an Alphabet object and the rates
#       p<-GeneralSubstitution(alphabet=BinaryAlphabet(), rate.list=list("1->0"=1,"0->1"=3))
#	# get the QMatrix object from p
#	m<-p$QMatrix
#	# get the scaled rate matrix from m
#	m$scaledMatrix
# } 
# 
# @author 
# 
# \seealso{ 
# 	@seeclass 
# } 
# 
#*/###########################################################################
setMethodS3(
	"getScaledMatrix", 
	class="QMatrix", 
	function(
		this,
		...
	){

		this$.rate.matrix;

	},
	private=FALSE,
	protected=FALSE,
	overwrite=FALSE,
	conflict="warning",
	validators=getOption("R.methodsS3:validators:setMethodS3")
);

##	
## Method: setScaledMatrix
##	
###########################################################################/**
#
# @RdocMethod setScaledMatrix
#
# @title "Forbidden action: setting the scaled rate matrix for a QMatrix object"
#
# \description{
#       @get "title".
# }
#
# @synopsis
#
# \arguments{
#       \item{this}{An object.}
#       \item{value}{Not used.}
#       \item{...}{Not used.}
# }
#
# \value{
#	Throws an error.
# }
#
# @author
#
# \seealso{
#       @seeclass
# }
#
#*/###########################################################################
setMethodS3(
	"setScaledMatrix", 
	class="QMatrix", 
	function(
		this,
		value,
		...
	){

		virtualAssignmentForbidden(this);

	},
	private=FALSE,
	protected=FALSE,
	overwrite=FALSE,
	conflict="warning",
	validators=getOption("R.methodsS3:validators:setMethodS3")
);

##	
## Method: Scale
##	
###########################################################################/**
#
# @RdocMethod Scale
# 
# @title "Scale the scaled rate matrix stored in a QMatrix object by the provided factor" 
# 
# \description{ 
#	@get "title".
#
#	This methods sets the scaled rate matrix to \code{unscaled_matrix * constant}.
# } 
# 
# @synopsis 
# 
# \arguments{ 
# 	\item{this}{A QMatrix object.} 
#	\item{constant}{The scaling factor (a numeric vector of length one).}
# 	\item{...}{Not used.} 
# } 
# 
# \value{ 
# 	The QMatrix object (invisible).
# } 
# 
# \examples{
#	# create a QMatrix object
#	# , provide Alphabet object and rates
#	m<-QMatrix(name="Susie Q", alphabet=BinaryAlphabet(), rate.list=list("1->0"=2,"0->1"=3))
#	# get object summary
#	summary(m)
#	# perform scaling
#	Scale(m, 1/0.666)
#	# get object summary
#	summary(m)
# } 
# 
# @author 
# 
# \seealso{ 
# 	@seeclass 
# } 
# 
#*/###########################################################################
setMethodS3(
	"Scale", 
	class="QMatrix", 
	function(
		this,
		constant=NA,
		...
	){

	if(!exists(x="PSIM_FAST")){

		if(missing(constant)){
			throw("No scaling constant specified!\n");
		}
		if(!is.numeric(constant)){
			throw("Scaling constant must be numeric!\n");
		}
	}
			 
			 # Set the rescaled matrix to the original matrix
			 # multiplied by the given constant:
			 this$.rate.matrix<-(this$.orig.matrix * constant);
			 # store the current rescaling constant:
			 this$.norm.const<-constant;
			 return(invisible(this));
	},
	private=FALSE,
	protected=FALSE,
	overwrite=FALSE,
	conflict="warning",
	validators=getOption("R.methodsS3:validators:setMethodS3")
);

##	
## Method: as.character
##	
###########################################################################/**
#
# @RdocMethod as.character
# 
# @title "Return the character representation of a QMatrix object" 
# 
# \description{ 
#	@get "title".
# } 
# 
# @synopsis 
# 
# \arguments{ 
# 	\item{x}{A QMatrix object.} 
# 	\item{...}{Not used.} 
# } 
# 
# \value{ 
# 	A character vector of length one.
# } 
# 
# \examples{
#	# create a QMatrix object
#	m<-QMatrix()
#	# get the character representation
#	as.character(m)
#	# the same, but implicitly
#	m
# } 
# 
# @author 
# 
# \seealso{ 
# 	@seeclass 
# } 
# 
#*/###########################################################################
setMethodS3(
	"as.character", 
	class="QMatrix", 
	function(
		x,
		...
	){
		
		x$id

	},
	private=FALSE,
	protected=FALSE,
	overwrite=FALSE,
	conflict="warning",
	validators=getOption("R.methodsS3:validators:setMethodS3")
);

##	
## Method: print
##	
###########################################################################/**
#
# @RdocMethod print
# 
# @title "Print the character representation of a QMatrix object" 
# 
# \description{ 
#	@get "title".
# } 
# 
# @synopsis 
# 
# \arguments{ 
# 	\item{x}{A QMatrix object.} 
# 	\item{...}{Not used.} 
# } 
# 
# \value{ 
# 	The character representation of the QMatrix object.
# } 
# 
# \examples{
#	# create a QMatrix object
#	m<-QMatrix()
#	# print the character representation
#	print(m)
#	# the same, but implicitly
#	m
# } 
# 
# @author 
# 
# \seealso{ 
# 	@seeclass 
# } 
# 
#*/###########################################################################
setMethodS3(
	"print", 
	class="QMatrix", 
	function(
		x,
		...
	){

		print.default(x$.orig.matrix);

	},
	private=FALSE,
	protected=FALSE,
	overwrite=FALSE,
	conflict="warning",
	validators=getOption("R.methodsS3:validators:setMethodS3")
);

##	
## Method: summary
##	
###########################################################################/**
#
# @RdocMethod summary
#
# @title "Summarize the properties of an object"
#
# \description{
#       @get "title".
# }
#
# @synopsis
#
# \arguments{
#       \item{object}{An object}
#       \item{...}{Not used.}
# }
#
# \value{
#  Returns a PSRootSummary object.
# }
#
# \examples{
#
#       # create an object
#       a<-QMatrix(alphabet=BinaryAlphabet(), rate.list=list("0->1"=1,"1->0"=3))
#       # get a summary
#       summary(a)
# }
#
# @author
#
# \seealso{
#       @seeclass
# }
#
#*/###########################################################################
setMethodS3(
	"summary", 
	class="QMatrix", 
	function(
		object,
		...
	){

		this<-object;
		this$.summary$"Name"<-this$name;
		this$.summary$"Id"<-this$id;
		this$.summary$"Attached to process"<-this$process;
	  this$.summary$"Unscaled rate matrix"<-paste( "\n\n\t",paste(capture.output(print(this$.orig.matrix)),collapse="\n\t"),"\n",sep="");
		this$.summary$"Scaling factor"<-this$.norm.const;
	  this$.summary$"Scaled rate matrix"<-paste( "\n\n\t",paste(capture.output(print(this$.rate.matrix)),collapse="\n\t"),"\n",sep="");
	
		
		NextMethod();

	},
	private=FALSE,
	protected=FALSE,
	overwrite=FALSE,
	conflict="warning",
	validators=getOption("R.methodsS3:validators:setMethodS3")
);

##
## setId
##
###########################################################################/**
#
# @RdocMethod setId
#
# @title "Forbidden action: setting the unique identifier for a QMatrix object"
#
# \description{
#       @get "title".
# }
#
# @synopsis
#
# \arguments{
#       \item{this}{An object.}
#       \item{value}{Not used.}
#       \item{...}{Not used.}
# }
#
# \value{
#	Throws an error.
# }
#
# @author
#
# \seealso{
#       @seeclass
# }
#
#*/###########################################################################
setMethodS3(
  "setId",
  class="QMatrix",
  function(
    this,
    value,
    ...
  ){

  throw("Id is generated automatically and it cannot be set!\n");

  },
  private=FALSE,
  protected=FALSE,
  overwrite=FALSE,
  conflict="warning",
  validators=getOption("R.methodsS3:validators:setMethodS3")
);

##
## Method: getId
##
###########################################################################/**
#
# @RdocMethod	getId
# 
# @title "Get the unique identifier of a QMatrix object" 
# 
# \description{ 
#	@get "title".
# } 
# 
# @synopsis 
# 
# \arguments{ 
# 	\item{this}{A QMatrix object.} 
# 	\item{...}{Not used.} 
# } 
# 
# \value{ 
# 	A character vector of length one.
# } 
# 
# \examples{
#	# create a QMatrix object
#	m<-QMatrix()
#	# get object id
#	getId(m)
#	# get object id via virtual field
#	m$id
# } 
# 
# @author 
# 
# \seealso{ 
# 	@seeclass 
# } 
# 
#*/###########################################################################
setMethodS3(
  "getId",
  class="QMatrix",
  function(
    this,
    ...
  ){

  this.class<-class(this)[1];
  id<-paste(this.class,this$.name,hashCode(this),sep=":");

  },
  private=FALSE,
  protected=FALSE,
  overwrite=FALSE,
  conflict="warning",
  validators=getOption("R.methodsS3:validators:setMethodS3")
);


##
## Method: getWriteProtected
##
###########################################################################/**
#
# @RdocMethod getWriteProtected
#  
# @title "Check if the object is write protected" 
# 
# \description{ 
#       @get "title".
#
#	QMatrix object do not have a write protection flag of their own, but they use the one from the 
#	associated Process object.
#       Write protected objects cannot be modified through get/set methods and virtual fields.
# } 
# 
# @synopsis 
# 
# \arguments{ 
#       \item{this}{A QMatrix object.} 
#       \item{...}{Not used.} 
# } 
# 
# \value{ 
#       TRUE or FALSE
# } 
# 
# 
# @author 
# 
# \seealso{ 
#       getWriteProtected.Process
# } 
# 
#*/###########################################################################
setMethodS3(
  "getWriteProtected",
  class="QMatrix",
  function(
    this,
    ...
  ){

	# return false if no process is attached:
	if(!is.Process(this$.process)) {
			return(FALSE);
	}
	else {
			# The flag from the parent process is used!
			return(getWriteProtected(this$.process));
	}

  },
  private=FALSE,
  protected=FALSE,
  overwrite=FALSE,
  conflict="warning",
  validators=getOption("R.methodsS3:validators:setMethodS3")
);

##
## Method: is.QMatrix
##
###########################################################################/**
#
# @RdocDefault is.QMatrix
# 
# @title "Check if an object is an instance of the QMatrix class" 
# 
# \description{ 
#       @get "title".
# } 
# 
# @synopsis 
# 
# \arguments{ 
#       \item{this}{An object.} 
#       \item{...}{Not used.} 
# } 
# 
# \value{ 
#       TRUE or FALSE.
# } 
# 
# \examples{
#       # create some objects
#       m<-QMatrix()
#       p<-Process()
#       # chek if they inherit from QMatrix
#       is.GeneralSubstitution(m)
#       is.GeneralSubstitution(p)
# } 
# 
# @author 
# 
# \seealso{ 
#       @seeclass
# } 
# 
#*/###########################################################################
setMethodS3(
  "is.QMatrix",
  class="default",
  function(
    this,
    ...
  ){

    if(!is.PSRoot(this)) {return(FALSE)}
    if(!is.null(this$.is.process)){return(TRUE)}
    if ( inherits(this, "QMatrix")) {
      this$.is.q.matrix<-TRUE;
      return(TRUE);
    } else {
      return(FALSE)
    }

  },
  private=FALSE,
  protected=FALSE,
  overwrite=FALSE,
  conflict="warning",
  validators=getOption("R.methodsS3:validators:setMethodS3")
);


##
## Method: setWriteProtected
##
###########################################################################/**
#
# @RdocMethod setWriteProtected
#  
# @title "Set the write protection field for a QMatrix object" 
# 
# \description{ 
#       @get "title".
#
#	QMatrix object do not have a write protection flag of their own, but they use the one from the 
#	associated Process object.
#       Write protected objects cannot be modified through get/set methods and virtual fields.
# } 
# 
# @synopsis 
# 
# \arguments{ 
#       \item{this}{A QMatrix object.} 
#	\item{value}{A logical vector of size one.}
#       \item{...}{Not used.} 
# } 
# 
# \value{ 
#       TRUE or FALSE
# } 
# 
# 
# @author 
# 
# \seealso{ 
#       setWriteProtected.Process
# } 
# 
#*/###########################################################################
setMethodS3(
  "setWriteProtected",
  class="QMatrix",
  function(
    this,
		value,
    ...
  ){

		throw("The QMatrix objects use the write protection flags of the enclosing substitution process, modify that (if exists)!\n");

  },
  private=FALSE,
  protected=FALSE,
  overwrite=FALSE,
  conflict="warning",
  validators=getOption("R.methodsS3:validators:setMethodS3")
);

##
## Method: .checkWriteProtection
##
setMethodS3(
  ".checkWriteProtection",
  class="QMatrix",
  function(
    this,
    ...
  ){

   if(!exists(x="PSIM_FAST")){
	return(FALSE);
   }
   
 	if(this$writeProtected) {throw("Cannot set value because the object is write protected!\n")}
    	else {return(FALSE)}

  },
  private=FALSE,
  protected=FALSE,
  overwrite=FALSE,
  conflict="warning",
  validators=getOption("R.methodsS3:validators:setMethodS3")
);