Nothing
R/QMatrix.R
In phylosim: Flexible Simulations of Biological Sequence Evolution
##
## 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")
);
Try the phylosim package in your browser
Any scripts or data that you put into this service are public.
phylosim documentation built on Nov. 22, 2019, 1:07 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
##
## 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")
);
Try the phylosim package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.