##	
## Copyright 2009 Botond Sipos	
## See the package description for licensing information.	
##	
##########################################################################/** 
#
# @RdocClass PhyloSim
# \alias{phylosim}
# 
# @title "The PhyloSim class"
# 
# \description{ 
#
#	PhyloSim is an extensible object-oriented framework for the Monte Carlo simulation 
#	of sequence evolution written in 100 percent \code{R}.
#	It is built on the top of the \code{\link{R.oo}} and \code{\link{ape}} packages and uses 
#	Gillespie's direct method to simulate substitutions, insertions and deletions.
#
#	Key features offered by the framework:	
#	\itemize{
#	\item  Simulation of the evolution of a set of discrete characters with arbitrary states evolving 
#		by a continuous-time Markov process with an arbitrary rate matrix.
#	\item Explicit implementations of the most popular substitution models (for nucleotides, amino acids and codons).
# 	\item Simulation under the popular models of among-sites rate variation, like the gamma (+G) and invariants plus gamma (+I+G) models.
#	\item The possibility to simulate with arbitrarily complex patterns of among-sites rate variation by setting the site specific rates according to any \code{R} expression.
#	\item Simulation with one or more separate insertion and/or deletion processes acting on the sequences, which sample the insertion/deletion length from an arbitrary discrete distribution or an \code{R} expression (so all the probability distributions implemented in \code{R} are readily available for this purpose).
#	\item Simulation of the effects of variable functional constraints over the sites by site-process-specific insertion and deletion tolerance parameters, which determine the rejection probability of a proposed insertion/deletion.
#	\item The possibility of having a different set of processes and site-process-specific parameters for every site, which allow for an arbitrary number of partitions in the simulated data.
#	\item Simulation of heterotachy and other cases of non-homogeneous evolution by allowing the user to set "node hook" functions altering the site properties at internal nodes of the phylogeny.
#	\item The possibility to export the counts of various events ("branch statistics") as phylo objects (see \code{\link{exportStatTree.PhyloSim}}).
#	}
#
#	General notes:
#	\itemize{
#	\item The \code{Sequence} objects have no "immortal links". The simulation
#	is aborted if the sequence length shrinks to zero. It is up to the user 
#	to choose sensible indel rates and sequence lengths to prevent that.
#	\item The sites near the beginning and end of the sequences have less sites proposing
#	insertion and deletion events around the so the insertion and deletion processes 
#	have an "edge effect". The user can simulate
#	realistic flanking sequences to alleviate the edge effect in the simulation settings where
#	it may be an issue.
# }
#
#	Notes on performance: 
#	\itemize{
#	\item The pure \code{R} implementation offers felxibility, but also comes
#	with a slower simulation speed. If the \code{PSIM_FAST} object is present in the environment, a "fast & careless"
#	mode is enabled. In this mode most of the error checking is skipped, increasing the speed.
#	It is recomended that simulations are only run in fast mode if you are sure that the simulation
#	settings are free from errors. It is probably a good practice to set up the simulations in normal mode
#	with short sequences and enable fast mode when running the actual simulation with long sequences.
#	\item Please note, that no "branch statistics" are saved in fast mode.
#	\item Logging also has a negative impact on performance, so it's not a good idea to run
#	large simulations with the logging enabled.
#	\item The time needed to run a simulation depends not only on the number of the sites, 
#	but also on the length of the tree.
#	\item Constructing \code{Sequence} objects with large number of sites is expensive. Avoid doing
#	that inside a cycle.
#	\item In the case of \code{Sequence} objects with a large number of sites (more than 10 000) the 
#	amount of available memory can be limiting as well.
#  }
#
#	The examples below demonstrate only some more common simulation settings,
#	the framework offers much more flexibility. See the package
#	vignette (\code{vignette("PhyloSim",package="phylosim")}) and the
#	examples directory (\url{http://github.com/sbotond/phylosim/tree/master/examples/}) 
#	for additional examples.
#	
#	@classhierarchy
# }
#
#
# \references{
#	Gillespie, DT (1977) Exact stochastic simulation of coupled chemical reactions - 
#	J. Phys. Chem. 81 (25):2340-2361 \url{http://dx.doi.org/10.1021/j100540a008}
# }
#	
# @synopsis
#	
# \arguments{
# 	\item{phylo}{A rooted phylo object, constructed by the APE package.}
# 	\item{root.seq}{A valid Sequence object with Process objects attached. Used as the starting sequence during simulation.}
# 	\item{name}{The name of the object (a character vector of length one).}
# 	\item{log.file}{Name of the file used for logging.}
# 	\item{log.level}{An integer specifying the verbosity of logging (see \code{\link{setLogLevel.PhyloSim}}).}
# 	\item{...}{Not used.}
#	}
# 
# \section{Fields and Methods}{ 
# 	@allmethods
# }
# 
# \examples{ 
#   set.seed(1)
#	## The following examples demonstrate
#	## the typical use of the framework.
#	## See the package vignette and
#	## \url{http://github.com/sbotond/phylosim/tree/master/examples/}
#
#	## The ll() method gives information about the methods defined
#	## in the immediate class of an object.
#	## Useful when exploring the framework.
#
#	s<-Sequence()
#	ll(s)
#	ll(PhyloSim())
#	ll(GTR())
#
#	## Example 1 - A short simulation:
#	## simulate nucleotide seqeunces and display 
#	## the resulting alignment matrix.
#
#	Simulate(
#		PhyloSim(phy=rcoal(3),
#       root=NucleotideSequence(string="ATGC", proc=list(list(JC69())) ) )
#	)$alignment
#
#	# Construct a phylo object for the following
#	# simulations, scale total tree length to 1:
#
#	tmp<-PhyloSim(phylo=rcoal(3))
#	scaleTree(tmp,1/tmp$treeLength)
#	tmp$treeLength
#	t<-tmp$phylo
#
#	## Example 3 - simulating rate variation,
#	## insertions and deletions.
#	## See the examples/example_3_clean.R file
#	## in the phylosim GitHub repository.
#
#	# construct a GTR process object
#       gtr<-GTR(
#		name="MyGTR",
#               rate.params=list(
#                       "a"=1, "b"=2, "c"=3,
#                       "d"=1, "e"=2, "f"=3
#               ),
#               base.freqs=c(2,2,1,1)/6
#       )
#	# get object summary
#	summary(gtr)
#	# display a bubble plot
#	plot(gtr)
#
#	# construct root sequence object
#	s<-NucleotideSequence(length=20)
#	# attach process via virtual field
#	s$processes<-list(list(gtr))
#	# sample states from the equilibrium
#	# distribution of the attached processes
#	sampleStates(s)
#	# create among-site rate variation by sampling
#	# the "rate.multiplier" site-process-specific parameter
#	# from a discrete gamma distribution (GTR+G).
#	plusGamma(s,gtr,shape=0.1)
#	# make the range 11:12 invariable
#	setRateMultipliers(s,gtr,0,11:12)
#	# get the rate multipliers for s and gtr
#	getRateMultipliers(s,gtr)
#
	# construct a deletion process object
#	# proposing lengths in the range 1:3
#	d<-DiscreteDeletor(
#		rate=0.1,
#		name="MyDel",
#		sizes=c(1:3),
#		probs=c(3/6,2/6,1/6)
#	)
#	# get object 
#	summary(d)
#	# plot deletion length distribution
#	plot(d)
#	# attach deletion process d to sequence s
#	attachProcess(s,d)
# 	# create a region rejecting all deletions
#       setDeletionTolerance(s,d,0,11:12)
#
#	# construct an insertion process object
#	# proposing lengths in the range 1:3
#	i<-DiscreteInsertor(
#		rate=0.1,
#		name="MyDel",
#		sizes=c(1:2),
#		probs=c(1/2,1/2),
#		template.seq=NucleotideSequence(length=1,processes=list(list(JC69())))
#	) 
# 	# states will be sampled from the JC69 equilibrium distribution
#	# get object 
#	summary(i)
#	# plot insertion length distribution
#	plot(i)
#	# attach insertion process i to sequence s
#	attachProcess(s,i)
# 	# create a region rejecting all insertions
#       setInsertionTolerance(s,i,0,11:12)
#
#	# plot total site rates
#	plot(s)
#	# construct simulation object
#	sim<-PhyloSim(root.seq=s, phylo=t)
#	# get object summary
#	summary(sim)
#	# plot tree
#	plot(sim)
#	# run simulation
#	Simulate(sim)
# 	# get the list of recorded per-branch event counts
# 	getBranchEvents(sim)
#	# export the number of substitutions as a phylo object
#	subst<-exportStatTree(sim,"substitution")
#	# plot the exported phylo object
#	plot(subst)
#	# plot tree and alignment
#	plot(sim)
#	# save and display alingment
#	file<-paste("PhyloSim_dummy_fasta_",Sys.getpid(),".fas",sep="");
#	saveAlignment(sim,file=file);
#	# print out the Fasta file
#	cat(paste(scan(file=file,what=character(),sep="\n"),collapse="\n"));cat("\n");
#	# delete Fasta file
#	unlink(file);
#
#   # See \url{http://github.com/sbotond/phylosim/tree/master/examples/examples_class.R}
#   # for the full list of PhyloSim constructor examples.
#   
#	## See the package vignette and 
#	## the GitHub repository for even more examples.
# }
# 
# @author
#
# \seealso{ 
#	\code{\link{PSRoot} \link{Alphabet} \link{AminoAcidAlphabet} 
#	\link{AminoAcidSequence} \link{AminoAcidSubst}
#	\link{AnyAlphabet} \link{BinaryAlphabet} \link{BinarySequence} 
#	\link{BinarySubst} \link{BrownianInsertor} \link{CodonAlphabet}
#	\link{CodonSequence} \link{CodonUNREST} \link{ContinuousDeletor}
#	\link{ContinuousInsertor} \link{cpREV} \link{DiscreteDeletor}
#	\link{DiscreteInsertor} \link{Event} \link{F81} \link{F84}
#	\link{FastFieldDeletor} \link{GeneralDeletor}
#	\link{GeneralInDel} \link{GeneralInsertor} \link{GeneralSubstitution} 
#	\link{GTR} \link{GY94} \link{HKY} \link{JC69} \link{JTT} \link{JTT.dcmut}
#	\link{K80} \link{K81} \link{LG} \link{mtArt} \link{mtMam} \link{mtREV24}
#	\link{MtZoa} \link{NucleotideAlphabet} \link{NucleotideSequence} \link{PAM} 
#	\link{PAM.dcmut} \link{PhyloSim} \link{Process} \link{QMatrix} \link{Sequence}
#	\link{Site} \link{T92} \link{TN93} \link{UNREST} \link{WAG}}
# }
# 
#*/###########################################################################
setConstructorS3(
"PhyloSim",
  function(
		phylo=NA,
		root.seq=NA,
		name=NA,
		log.file=NA,
		log.level=-1, # no loggin is performed by default
		... 
		)	{

		this<-PSRoot();	
		this<-extend(this,
			"PhyloSim",
			.name="Anonymous",
			.phylo=NA,
			.root.sequence=NA, 
			.sequences=list(),		# references to the sequence objects
			.node.hooks=list(),		# references to the node hook functions.
			.branch.stats=list(), # branch statistics.
			.alignment=NA,				# the resulting alignment in fasat format.
			.log.file=NA, 				# the name of the log file.
			.log.connection=NA,		# connection for the log file.
			.log.level=NA					# log level
		);

		if(!all(is.na(phylo))){
			this$phylo<-phylo;
		}

		if(!all(is.na(root.seq))){
			this$rootSeq<-root.seq;
		}
	
		if(!missing(name)){
			this$name<-name;
		}

		if(!missing(log.file)){
			this$logFile<-log.file;
		} else {	
			# Setting default log file:
			tmp<-this$id;
			tmp<-gsub(":","_",tmp);
			this$logFile<-paste(tmp,".log",sep="");
		}

		# Setting log level:
		this$logLevel<-log.level;

		return(this);

  },
  enforceRCC=TRUE
);

##	
## Method: checkConsistency
##	
###########################################################################/**
#
# @RdocMethod	checkConsistency
# 
# @title "Check object consistency"
# 
# \description{ 
#		@get "title".
# } 
# 
# @synopsis 
#
# \arguments{ 
#       \item{this}{An 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="PhyloSim", 
	function(
		this,
		...
	){

      may.fail<-function(this) {
					
				# Checking the name:	
				this$name<-this$name;
				# Checking the phylo object:
				if (!any(is.na(this$.phylo)) & !is.phylo(this$.phylo) ){
					throw("The phylo object is invalid!\n");
				}
				# Checking the log level:
				if(!is.numeric(this$.log.level) | (length(this$.log.level) != 1) ){
					throw("The log level must be numeric vector of length 1!\n");
				}
				# Checking lof file:
				if(!is.character(this$.log.file) | (length(this$.log.level) != 1) ){
					throw("The log file must be charcter vector of length 1!\n");
				}
				# Checking the sequences:
				for (seq in this$.sequences){
					if(is.Sequence(seq)){
						checkConsistency(seq);
					}
				}
				# Checking node hooks:
				for (hook in this$.node.hooks){
					if(!is.null(hook) & !is.function(hook)){
						throw("Invalid node hook found!\n");
					}
				}
				# Checking the alignment:
				if(!any(is.na(this$.alignment))){
					.checkAlignmentConsistency(this, this$.alignment);
				}
		
      }
      tryCatch(may.fail(this));
			return(invisible(TRUE));

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

##	
## Method: is.phylo.default
##	
###########################################################################/**
#
# @RdocDefault is.phylo
# 
# @title "Check if an object is an instance of the phylo class" 
# 
# \description{ 
#	@get "title".
#	Phylo objects are created by the \pkg{APE} package. This method just return the value of \code{inherits(this,"phylo")}.
# } 
# 
# @synopsis 
# 
# \arguments{ 
# 	\item{this}{An object.} 
# 	\item{...}{Not used.} 
# } 
# 
# \value{ 
# 	TRUE or FALSE.
# } 
# 
# \examples{
#	# load APE
#	library(ape);
#	# create some objects
#	o1<-Object();
#	o2<-rcoal(3);
#	# check if they are phylo objects
#	is.phylo(o1);
#	is.phylo(o2);			
#
# } 
# 
# @author 
# 
# \seealso{ 
# 	The \pkg{ape} package.
# } 
# 
#*/###########################################################################
setMethodS3(
	"is.phylo", 
	class="default", 
	function(
		this,
		...
	){

		inherits(this,"phylo");

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

##	
## Method: setPhylo
##	
###########################################################################/**
#
# @RdocMethod setPhylo
# 
# @title "Set the phylo object for a PhyloSim object" 
# 
# \description{ 
#	@get "title".
#
#	The internal structure of the provided phylo object is reordered in a cladeweise fashion.
# } 
# 
# @synopsis 
# 
# \arguments{ 
# 	\item{this}{A PhyloSim object.} 
# 	\item{value}{A phylo object created by the \pkg{ape} package.} 
# 	\item{...}{Not used.} 
# } 
# 
# \value{ 
# 	A phylo object or FALSE.
# } 
# 
# \examples{
#	#create a PhyloSim object
#	sim<-PhyloSim();
#	# creat a phylo object
#	tree<-rcoal(3);
#	# get/set phylo object
#	setPhylo(sim,tree);
#	getPhylo(sim,tree);
#	# get/set phylo object via virtual field
#	sim$tree<-rcoal(5);
#	sim$tree;
# } 
# 
# @author 
# 
# \seealso{ 
# 	The PhyloSim class, the \pkg{ape} package.
# } 
# 
#*/###########################################################################
setMethodS3(
	"setPhylo", 
	class="PhyloSim", 
	function(
		this,
		value,
		...
	){

		if(missing(value)){
			throw("No object provided!\n");
		}
		else if(!is.phylo(value)){
			throw("The new value must be a \"phylo\" object!\n");
		}
		else if(!is.rooted(value)){
			throw("The new value must be a rooted \"phylo\" object!\n");
		}
		else {

			.checkTipLabels(value);
			this$.phylo<-value;
			this$.phylo<-reorder(this$.phylo, order="cladewise");
			for (i in this$nodes){
				this$.sequences[[i]]<-NA;
			}
			return(this$.phylo);

		}
		return(FALSE);

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

##	
## Method: .checkTipLabels
##	
setMethodS3(
	".checkTipLabels", 
	class="phylo", 
	function(
		this,
		...
	){

		for(label in this$tip.label){
			if(length(grep("^Node \\d+$",label,perl=TRUE,value=FALSE)) > 0){
					throw("Sorry, but the node labels matching \"Node \\d+\" are reserved for internal nodes! Blaming label: ",label,".\n");	
			}
			else if(length(grep("^Root node \\d+$",label,perl=TRUE,value=FALSE)) > 0){
					throw("Sorry, but the node labels matching \"Root node \\d+\" are reserved for the root node! Blaming label: ",label,".\n");	
			}
			
		}

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

##	
## Method: getPhylo
##	
###########################################################################/**
#
# @RdocMethod getPhylo
# 
# @title "Get the phylo object aggregated in a PhyloSim object" 
# 
# \description{ 
#	@get "title".
#
# } 
# 
# @synopsis 
# 
# \arguments{ 
# 	\item{this}{A PhyloSim object.} 
# 	\item{...}{Not used.} 
# } 
# 
# \value{ 
# 	A phylo object or NA.
# } 
# 
# \examples{
#	#create a PhyloSim object
#	sim<-PhyloSim();
#	# creat a phylo object
#	tree<-rcoal(3);
#	# get/set phylo object
#	setPhylo(sim,tree);
#	getPhylo(sim,tree);
#	# get/set phylo object via virtual field
#	sim$tree<-rcoal(5);
#	sim$tree;
# } 
# 
# @author 
# 
# \seealso{ 
# 	The PhyloSim class, the \pkg{ape} package.
# } 
# 
#*/###########################################################################
setMethodS3(
	"getPhylo", 
	class="PhyloSim", 
	function(
		this,
		...
	){

		this$.phylo;

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

##	
## Method: setRootSeq
##	
###########################################################################/**
#
# @RdocMethod setRootSeq
# 
# @title "Set the root sequence for a PhyloSim object" 
# 
# \description{ 
#	@get "title".
#
#	The root sequence will be used as a starting point for the simulation. The phylo object must be set before
#	trying to set the root sequence object.
# } 
# 
# @synopsis 
# 
# \arguments{ 
# 	\item{this}{A PhyloSim object.} 
# 	\item{value}{A valid Sequence object.} 
# 	\item{...}{Not used.} 
# } 
# 
# \value{ 
# 	The root Sequence object if succesfull, FALSE otherwise.
# } 
# 
# \examples{
#	# create some objects
#	sim<-PhyloSim(phylo=rcoal(3));
#	seq<-NucleotideSequence(string="ATGCC");
#	# set/get root sequence
#	setRootSeq(sim, seq);
#	getRootSeq(sim, seq);
#	# set/get root sequence via virtual field
#	sim$rootSeq<-BinarySequence(string="111000111000");
#	sim$rootSeq;
#
# } 
# 
# @author 
# 
# \seealso{ 
# 	@seeclass Sequence Process
# } 
# 
#*/###########################################################################
setMethodS3(
	"setRootSeq", 
	class="PhyloSim", 
	function(
		this,
		value,
		...
	){

		if(missing(value)){
			throw("No object provided!\n");
		}
		else if(!is.Sequence(value)){
			throw("The new value must be a sequence object!\n");
		}
		else {

			this$.root.sequence<-clone(value);
			this$.root.sequence$name<-paste("Root node",this$rootNode);
        
            # Call garbage collection:
            gc();
            gc();

			return(this$.root.sequence);

		}
		return(FALSE);


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

##	
## Method: getRootSeq
##	
###########################################################################/**
#
# @RdocMethod getRootSeq
# 
# @title "Get the root sequence aggregated by a PhyloSim object" 
# 
# \description{ 
#	@get "title".
#
# } 
# 
# @synopsis 
# 
# \arguments{ 
# 	\item{this}{A PhyloSim object.} 
# 	\item{...}{Not used.} 
# } 
# 
# \value{ 
# 	The root Sequence object or NA.
# } 
# 
# \examples{
#	# create some objects
#	sim<-PhyloSim(phylo=rcoal(3));
#	seq<-NucleotideSequence(string="ATGCC");
#	# set/get root sequence
#	setRootSeq(sim, seq);
#	getRootSeq(sim, seq);
#	# set/get root sequence via virtual field
#	sim$rootSeq<-BinarySequence(string="111000111000");
#	sim$rootSeq;
#
# } 
# 
# @author 
# 
# \seealso{ 
# 	@seeclass Sequence Process
# } 
# 
#*/###########################################################################
setMethodS3(
	"getRootSeq", 
	class="PhyloSim", 
	function(
		this,
		...
	){

			this$.root.sequence;


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

##	
## Method: as.character.PhyloSim
##	
###########################################################################/**
#
# @RdocMethod as.character
# 
# @title "Return the character representation of a PhyloSim object" 
# 
# \description{ 
#	@get "title".
#
#	The character representation is the identifier of the PhyloSim object as returned by the \code{getId} method.
# } 
# 
# @synopsis 
# 
# \arguments{ 
# 	\item{x}{A PhyloSim object.} 
# 	\item{...}{Not used.} 
# } 
# 
# \value{ 
# 	A character vector of length one.
# } 
# 
# \examples{
#	# create a PhyloSim object
#	o<-PhyloSim(name="MySim");
#	# get character representation
#	as.character(o)
# } 
# 
# @author 
# 
# \seealso{ 
# 	@seeclass 
# } 
# 
#*/###########################################################################
setMethodS3(
	"as.character", 
	class="PhyloSim", 
	function(
		x,
		...
	){

		return(getId(x));

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

##
## Method: getName
##
###########################################################################/**
#
# @RdocMethod getName
# 
# @title "Get the name of a PhyloSim object" 
# 
# \description{ 
#	@get "title".
# } 
# 
# @synopsis 
# 
# \arguments{ 
# 	\item{this}{A PhyloSim object.} 
# 	\item{...}{Not used.} 
# } 
# 
# \value{ 
# 	A character vector of length one.
# } 
# 
# \examples{
#	# create a PhyloSim object
#	o<-PhyloSim();
#	# set/get name
#	setName(o,"MySim");
#	getName(o,"MySim");
#	# set/get name via virtual field
#	o$name<-"George";
#	o$name
# } 
# 
# @author 
# 
# \seealso{ 
# 	@seeclass 
# } 
# 
#*/###########################################################################
setMethodS3(
  "getName",
  class="PhyloSim",
  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 PhyloSim object" 
# 
# \description{ 
#	@get "title".
# } 
# 
# @synopsis 
# 
# \arguments{ 
# 	\item{this}{A PhyloSim object.} 
#	\item{new.name}{A character vector of length one.}
# 	\item{...}{Not used.} 
# } 
# 
# \value{ 
# 	The new name.
# } 
# 
# \examples{
#	# create a PhyloSim object
#	o<-PhyloSim();
#	# set/get name
#	setName(o,"MySim");
#	getName(o,"MySim");
#	# set/get name via virtual field
#	o$name<-"George";
#	o$name
# } 
# 
# @author 
# 
# \seealso{ 
# 	@seeclass 
# } 
# 
#*/###########################################################################
setMethodS3(
  "setName",
  class="PhyloSim",
  function(
    this,
    new.name,
    ...
  ){

    this$.name<-as.character(new.name);
  },
  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 PhyloSim object" 
# 
# \description{ 
#	@get "title".
#	The unique identifier is the concatenation of the class, the object name as returned by getName() and the object hash 
#       as returned by hashCode().
# } 
# 
# @synopsis 
# 
# \arguments{ 
# 	\item{this}{A PhyloSim object.}
# 	\item{...}{Not used.} 
# } 
# 
# \value{ 
# 	A character vector of length one.
# } 
# 
# \examples{
#	# create a PhyloSim object
#	o<-PhyloSim(name="MySim");
#	# get id
#	getId(o);
#	# get id via virtual field
#	o$id;
# } 
# 
# @author 
# 
# \seealso{ 
# 	@seeclass 
# } 
# 
#*/###########################################################################
setMethodS3(
  "getId",
  class="PhyloSim",
  function(
    this,
    ...
  ){

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

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

##
## Method: setId
##
###########################################################################/**
#
# @RdocMethod setId
#
# @title "Forbidden action: setting the unique identifier of a PhyloSim 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="PhyloSim",
  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: Simulate
##
###########################################################################/**
#
# @RdocMethod Simulate
# 
# @title "Run a simulation according to a PhyloSim object" 
# 
# \description{ 
#	@get "title".
#
#	The phylo object and the root sequence must be set before attempting to run a simulation.
#	Also the bigRate of the root sequence must not be NA or zero, so at least one sane
#	Process object must be attached to the root sequence object.
# } 
# 
# @synopsis 
# 
# \arguments{ 
# 	\item{this}{A PhyloSim object.} 
# 	\item{quiet}{TRUE or FALSE (default).} 
# 	\item{...}{Not used.} 
# } 
# 
# \value{ 
# 	The PhyloSim object (invisible).
# } 
# 
# \examples{
#	# Create a PhyloSim object.
#	# Provide the phylo object 
#	# and the root sequence.
#	sim<-PhyloSim(
#		name="TinySim",
#		phylo=rcoal(3),
#		root.seq=NucleotideSequence(string="ATGC",processes=list(list(JC69())))
#	);
#	# Run the simulation
#	Simulate(sim);
#	# Print the resulting sequences
#	sim$sequences
#	# Print the resulting alignment
#	sim$alignment
# } 
# 
# @author 
# 
# \seealso{ 
# 	@seeclass 
# } 
# 
#*/###########################################################################
setMethodS3(
  "Simulate",
  class="PhyloSim",
  function(
    		this,
		quiet=FALSE,
    ...
  ){


		if(!is.phylo(this$.phylo)){
			throw("Cannot simulate because the phylo object is not set or it is invalid!\n");
		}
		# Check for the root sequence:
		else if(!is.Sequence(this$.root.sequence)){
			throw("Cannot simulate because the root sequence is not set or it is invalid!\n");
		}
		# Check bigRate validity:
		else if(is.na(this$.root.sequence$bigRate)){
			throw("Cannot simulate because the bigRate of the root sequence is NA!\n");
		}
		else{

			# Warn for zero bigRate:
			if(this$.root.sequence$bigRate == 0){
				warning("The bigRate of the root sequence is zero! You are running a pointless simulation!\n");
			}

			if(exists(x="PSIM_FAST")){
				if(!quiet){
cat("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
cat("!! WARNING: fast & careless mode is on, most of the error checking is omitted!  !!\n");
cat("!!    Please note that this also disables the saving of branch statistics.      !!\n");
cat("!!       You can go back to normal mode by deleting the PSIM_FAST object.       !!\n");
cat("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
				}
			Log(this,"WARNING: fast & careless mode is on, most of the error checking is omitted!");
			}
		
			# Attach root sequence to root node:
			Log(this,paste("Attaching root sequence ",this$.root.sequence$id,sep=""));
			attachSeqToNode(this, node=getRootNode(this),seq=this$.root.sequence);

			# Write protecting the root sequence:
			Log(this,paste("Write protecting root sequence ",this$.root.sequence$id,sep=""));
			this$.root.sequence$writeProtected<-TRUE;

			# Traverse the tree and simulate:
			Log(this,paste("Starting simulation on the object",this$id));	
			edge.counter<-1;
			n.edges<-this$nedges;
			for(edge in 1:n.edges){
				if(!quiet){ cat("Simulating edge",edge,"of", n.edges,"\n");}
				Log(this,paste("Starting to simulate edge",edge,"of",n.edges));	
				.simulateEdge(this,number=edge);
				edge.counter<-edge.counter+1;
			}
		}
		Log(this, "Simulation finished, building alignment!\n");
		this$.alignment<-.recoverAlignment(this);
		# Flush the log connection:
		if(!is.na(this$.log.connection)){
				close(this$.log.connection);
		}

        # Call the garbage collector:
        gc();
        gc();

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

##
## Method: .simulateEdge
##
setMethodS3(
  ".simulateEdge",
  class="PhyloSim",
  function(
    this,
		number=NA,
    ...
  ){

		# Get edge:
		edge<-getEdge(this, number);
		# Get parent node:
		start.seq<-getSeqFromNode(this, edge[[1,"from"]]);
		# Evolve sequence:
		new.seq<-.evolveBranch(this, start.seq=start.seq, branch.length=edge[1,"length"], old.node=edge[[1,"from"]],new.node=edge[[1,"to"]], branch.number=number);
		# Write protect the sequence:
		new.seq$writeProtected<-TRUE;
		# Attach sequence to children node:
		attachSeqToNode(this, node=edge[1,"to"], seq=new.seq);

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

##
## Method: .evolveBranch
##
setMethodS3(
  ".evolveBranch",
  class="PhyloSim",
  function(
    this,
    start.seq=NA,
    branch.length=NA,
		old.node=NA,
		new.node=NA,
		branch.number=NA,
    ...
  ){

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

		if(missing(start.seq)){
			throw("No starting sequence provided!\n");
		}
		else if(missing(branch.length)){
			throw("No branch length provided!\n");
		}
		else if(!is.numeric(branch.length)){
			throw("The branch length must be numeric!\n");
		}
	}
		if(.checkSeq(this, start.seq) ){
		
			# Cloning the starting sequence:
			seq<-clone(start.seq);
			
			# Set the name of the sequence object:
			if(is.tip(this, new.node)){
				seq$name<-this$tipLabels[[new.node]];
			}
			else {
				seq$name<-paste("Node",new.node);
			}
			.GillespieDirect(this, seq=seq, branch.length=branch.length, branch.number=branch.number);
			
			# Call the node hook if exists:
			hook<-this$.node.hooks[[as.character(new.node)]];
			if(!is.null(hook) & is.function(hook)){
				Log(this,paste("Calling node hook for node",new.node));
				seq<-hook(seq=seq);	
				if(!is.Sequence(seq)){
					throw("Node hook returned an invalid sequence object!\n");
				}
				else if(is.na(seq$bigRate)){
					throw("Node hook returned sequence with NA bigRate!\n");
				}
				else if(seq$bigRate == 0.0){
					throw("Node hook returned sequence with zero bigRate!\n");
				}
				else{
				 checkConsistency(seq, omit.sites=TRUE);
				}
			}

			# Return the resulting sequence object:
			return(seq);
		
		}
  },
  private=FALSE,
  protected=FALSE,
  overwrite=FALSE,
  conflict="warning",
  validators=getOption("R.methodsS3:validators:setMethodS3")
);

##
## Method: .GillespieDirect
##
setMethodS3(
  ".GillespieDirect",
  class="PhyloSim",
  function(
    this,
    seq=NA,
    branch.length=NA,
		branch.number=NA,
    ...
  ){

		Debug(this, paste("Branch length is",branch.length));
		
		# Initialize time:

		time<-0.0;
		
		# Sample the next waiting time until
		# the branch length is consumed:	

		while( (time<-time + rexp(n=1, rate=(big.rate<-getBigRate(seq)))) <= branch.length){

			# Generate a random number between zero and the bigRate:

			E<-runif(n=1,min=0,max=big.rate);

			# Identify the target site:

			site.number<-which(.getCumulativeRatesFast(seq) >= E)[[1]];

			# Get the events from the target site:

			site<-seq$.sites[[site.number]];
			site$.position<-site.number;
			events<-getEvents(site);
			site$.position<-NULL;
			
			# Get the rates:
			rates<-double();
			for(e in events){
				rates<-c(rates,e$.rate);
			}
		
			# Calculate the corresponding cumulative rates:	
			if(site.number > 1){
				rates<-cumsum(c(seq$.cumulative.rates[[site.number - 1]], rates));
			}
			else {
				rates<-cumsum(c(0.0, rates));
			}

			# Pick the event:

			event.number<-which(rates >= E)[[1]] - 1;
			event<-events[[event.number]];

			# Log the event:

			Log(this,paste("Performing event [",event$.name,"] at position",event$.position,"generated by the process",event$.process$.id));

			# Perform the event:

			event.details<-Perform(event);
			Debug(this,paste("Remaining branch length is",(branch.length-time) ));
			
			# Log event details:
	
					# Log deletion event details:
						if(event$.name == "Deletion"){
						Log(this,paste("The process",event$.process,"proposed to delete range",paste(event.details$range,collapse="--"),". Accepted:",event.details$accepted));
					}
					# Log insertion event details:
					else if(event$.name == "Insertion"){
						message<-paste("The process ",event$.process," proposed insertion at position ",event.details$position,". Accepted: ",event.details$accepted,sep="");
						if(event.details$accepted == TRUE){
						message<-paste(message,"."," Insert length was ",event.details$length,sep="");
						}
						Log(this, message);
					}

			# Update branch statistics if not in fast mode:
			if(!exists(x="PSIM_FAST")){
				.UpdateBranchStats(this,event,event.details, branch.number);
			}

			# Abort if sequence length shrunk to zero:

			if(seq$.length == 0){
				message<-paste("Terminating the simulation because the length of the sequence ",seq$name," shrunk to zero! Please be more careful when tuning the indel rates!\n");
				Log(this, message);
				throw(message);
			}	

		} #/while

        # Calling garbage collection:
        gc();
        gc();
		return(seq);

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

##
## Method: attachSeqToNode
##
###########################################################################/**
#
# @RdocMethod	attachSeqToNode
# 
# @title "Assotiate a Sequence object with a given node of a phylo object aggregated by a PhyloSim object" 
# 
# \description{ 
#	@get "title".
#
#	This method is mainly used internally.
#	
# } 
# 
# @synopsis 
# 
# \arguments{ 
# 	\item{this}{A PhyloSim object.} 
#	\item{node}{Node identifier.}
#	\item{seq}{A Sequence object.}
# 	\item{...}{Not used.} 
# } 
# 
# \value{ 
# 	The PhyloSim object (invisible).
# } 
# 
# 
# @author 
# 
# \seealso{ 
# 	@seeclass 
# } 
# 
#*/###########################################################################
setMethodS3(
  "attachSeqToNode",
  class="PhyloSim",
  function(
		this,
		node=NA,
		seq=NA,
    ...
  ){

		if(!is.phylo(this$.phylo)){
			throw("The phylo object is not set, sequence to node is not possible!\n");
		}
		if(missing(node)){
			throw("No node specified!\n");
		}
		else if(missing(seq)){
			throw("No sequence object given");
		}
		else if(.checkNode(this,node) & .checkSeq(this, seq)){
			
			if(is.Sequence(this$.sequences[[node]])){
				throw("The node has already an attached sequence. Detach that before trying to attach a new one!\n");
			}
			else {
				this$.sequences[[as.numeric(node)]]<-seq;
				return(invisible(this));
			}

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

##
## Method: attachHookToNode
##
###########################################################################/**
#
# @RdocMethod attachHookToNode
# 
# @title "Attach a callback function to a given node of a phylo object aggregated by a PhyloSim object" 
# 
# \description{ 
#	@get "title".
#
#	A "node hook" is a function which accepts a Sequence object through the named argument "seq" and returns a 
#	Sequence object. The node hook function must accept any object which inherits from the \code{Sequence} class!
#
#	After simulating the branch leading to the node, the resulting Sequence object is passed
#	to the node hook and the returned object is used to simulate the downstream branches.
#
#	By using node hooks the attached processes can be replaced during simulation, hence enabling the simulation of 
#	non-homogeneous sequence evolution.
# } 
# 
# @synopsis 
# 
# \arguments{ 
# 	\item{this}{A PhyloSim object.} 
# 	\item{node}{Node identifier.} 
# 	\item{fun}{A function (see above).} 
# 	\item{...}{Not used.} 
# } 
# 
# \value{ 
# 	The PhyloSim object (invisible).
# } 
# 
# \examples{
#	# Create a PhyloSim object.
#	# Provide the phylo object 
#	# and the root sequence.
#	sim<-PhyloSim(
#		name="TinySim",
#		phylo=rcoal(3),
#		root.seq=NucleotideSequence(string="ATGC",processes=list(list(JC69())))
#	);
#	# create a node hook function
#	hook<-function(seq=NA){
#		# replace the substitution process with F84
#		if(inherits(seq,"NucleotideSequence")){
#			cat("Replacing JC69 with F84.\n");
#			seq$processes<-list(list(F84(rate.params=list("Kappa" = 2)))); 
#		}
#		return(seq);
#	}
#	# attach hook function to node 5
#	attachHookToNode(sim,5,hook);
#	# Run the simulation
#	Simulate(sim);
#	# Check if the processes have been truly replaced
#	lapply(sim$sequences, getUniqueProcesses.Sequence)
#	# Print the resulting alignment
#	sim$alignment
# } 
# 
# @author 
# 
# \seealso{ 
# 	@seeclass 
# } 
# 
#*/###########################################################################
setMethodS3(
  "attachHookToNode",
  class="PhyloSim",
  function(
    		this,
		node=NA,
		fun=NA,
    ...
  ){

		if(!is.phylo(this$.phylo)){
			throw("The phylo object is not set, attaching node hook is not possible!\n");
		}
		if(missing(node)){
			throw("No node specified!\n");
		}
		else if(missing(fun)){
			throw("No function given!");
		}
		else if(!is.function(fun)){
			throw("The argument \"fun\" must be a function!\n");
		}
		else if( length(intersect(names(formals(fun)), "seq")) == 0 ){
			throw("The function argument must have a an argument named \"seq\"");
		}
		else if(!is.Sequence(fun(Sequence(length=1)))){
      throw("The insert hook function must return a Sequence object!\n");
		}
		else if( .checkNode(this,node) ){
			if(is.function(this$.node.hooks[[as.character(node)]])){
				throw("The node has already an attached node hook. Detach that before trying to attach a new one!\n");
			}
			else {
				this$.node.hooks[[as.character(node)]]<-fun;
				return(invisible(this));
			}

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


##
## Method: .checkNode
##
setMethodS3(
  ".checkNode",
  class="PhyloSim",
  function(
    this,
		node=NA,
    ...
  ){

		if(missing(node)){
			throw("No node specified!\n");
		} else if( length(intersect(node, getNodes(this))) != 1){
			throw("The specified node is invalid!\n");	
		}
		else {
			return(TRUE);
		}
  },
  private=FALSE,
  protected=FALSE,
  overwrite=FALSE,
  conflict="warning",
  validators=getOption("R.methodsS3:validators:setMethodS3")
);

##
## Method: .checkSeq
##
setMethodS3(
  ".checkSeq",
  class="PhyloSim",
  function(
    this,
		seq=NA,
    ...
  ){

		if(missing(seq)){
			throw("No sequence specified!\n");
		} else if(!is.Sequence(seq)){
			throw("The sequence object is invalid!\n");	
		}
		else {
			return(TRUE);
		}
  },
  private=FALSE,
  protected=FALSE,
  overwrite=FALSE,
  conflict="warning",
  validators=getOption("R.methodsS3:validators:setMethodS3")
);

##
## Method: detachSeqFromNode
##
###########################################################################/**
#
# @RdocMethod	detachSeqFromNode
# 
# @title "Detach a Sequence object from a given node of a phylo object aggregated by a PhyloSim object" 
# 
# \description{ 
#	@get "title".
#
#	This method is mainly used internally.
#	
# } 
# 
# @synopsis 
# 
# \arguments{ 
# 	\item{this}{A PhyloSim object.} 
#	\item{node}{Node identifier.}
# 	\item{...}{Not used.} 
# } 
# 
# \value{ 
# 	The PhyloSim object (invisible).
# } 
# 
# @author 
# 
# \seealso{ 
# 	@seeclass 
# } 
# 
#*/###########################################################################
setMethodS3(
  "detachSeqFromNode",
  class="PhyloSim",
  function(
    		this,
		node=NA,
    		...
  ){

		if(missing(node)){
			throw("No node specified!\n");
		}
		else if( .checkNode(this,node) ){
			
				this$.sequences[[as.numeric(node)]]<-NA;
		}
	
		return(invisible(this));
  },
  private=FALSE,
  protected=FALSE,
  overwrite=FALSE,
  conflict="warning",
  validators=getOption("R.methodsS3:validators:setMethodS3")
);

##
## Method: detachHookFromNode
##
###########################################################################/**
#
# @RdocMethod	detachHookFromNode
# 
# @title "Detach a node hook function from a given node of a phylo object aggregated by a PhyloSim object" 
# 
# \description{ 
#	@get "title".
#
# } 
# 
# @synopsis 
# 
# \arguments{ 
# 	\item{this}{A PhyloSim object.} 
#	\item{node}{Node identifier.}
# 	\item{...}{Not used.} 
# } 
# 
# \value{ 
# 	The PhyloSim object (invisible).
# } 
# 
# \examples{
#	# Create a PhyloSim object.
#	# Provide the phylo object 
#	# and the root sequence.
#	sim<-PhyloSim(
#		name="TinySim",
#		phylo=rcoal(3),
#		root.seq=NucleotideSequence(string="ATGC",processes=list(list(JC69())))
#	);
#	# create a node hook function
#	hook<-function(seq=NA){
#		# replace the substitution process with F84
#		if(inherits(seq,"NucleotideSequence")){
#			cat("Replacing JC69 with F84.\n");
#			seq$processes<-list(list(F84(rate.params=list("Kappa" = 2)))); 
#		}
#		return(seq);
#	}
#	# attach hook function to node 5
#	attachHookToNode(sim,5,hook);
#	# detach hook from node 5
#	detachHookFromNode(sim,5);
#	# Run the simulation again
#	Simulate(sim);	# You should not see the message printed out by the "hook" function.
#
# } 
# 
# @author 
# 
# \seealso{ 
# 	attachHookToNode PhyloSim Simulate.PhyloSim
# } 
# 
#*/###########################################################################
setMethodS3(
  "detachHookFromNode",
  class="PhyloSim",
  function(
    this,
		node=NA,
    ...
  ){

		if(missing(node)){
			throw("No node specified!\n");
		}
		else if( .checkNode(this,node) ){
			
				this$.node.hooks[[as.character(node)]]<-NA;
		}
	
		return(invisible(this));
  },
  private=FALSE,
  protected=FALSE,
  overwrite=FALSE,
  conflict="warning",
  validators=getOption("R.methodsS3:validators:setMethodS3")
);


##
## Method: getSeqFromNode
##
###########################################################################/**
#
# @RdocMethod getSeqFromNode
# 
# @title "Get the Sequence object associated with a given node of a phylo object aggregated by a PhyloSim object" 
# 
# \description{ 
#	@get "title".
# } 
# 
# @synopsis 
# 
# \arguments{ 
# 	\item{this}{A PhyloSim object.} 
# 	\item{node}{Node identifier.} 
# 	\item{...}{Not used.} 
# } 
# 
# \value{ 
# 	A Sequence object.
# } 
# 
# \examples{
#	# Create a PhyloSim object.
#	# Provide the phylo object 
#	# and the root sequence.
#	sim<-PhyloSim(
#		name="TinySim",
#		phylo=rcoal(3),
#		root.seq=NucleotideSequence(string="ATG",processes=list(list(JC69())))
#	);
#	# get the sequence associated with node 5
#	getSeqFromNode(sim,5)	# Should be NA
#	# Run the simulation
#	Simulate(sim)
#	# try again
#	getSeqFromNode(sim,5)
# } 
# 
# @author 
# 
# \seealso{ 
# 	@seeclass 
# } 
# 
#*/###########################################################################
setMethodS3(
  "getSeqFromNode",
  class="PhyloSim",
  function(
    		this,
		node=NA,
    		...
  ){

		if(missing(node)){
			throw("No node specified!\n");
		}
		else if( .checkNode(this,node) ){
			
				return(this$.sequences[[as.numeric(node)]]);
		}

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

##
## Method: getSequences
##
###########################################################################/**
#
# @RdocMethod getSequences
# 
# @title "Gets all the Sequence objects associated with the nodes of a phylo object aggregated by a PhyloSim object" 
# 
# \description{ 
#	@get "title".
#
#	The order of the Sequence objects in the returned list reflects the identifiers of the associated nodes.
# } 
# 
# @synopsis 
# 
# \arguments{ 
# 	\item{this}{A PhyloSim object.} 
# 	\item{...}{Not used.} 
# } 
# 
# \value{ 
# 	A list of sequence objects.
# } 
# 
# \examples{
#	# Create a PhyloSim object.
#	# Provide the phylo object 
#	# and the root sequence.
#	sim<-PhyloSim(
#		name="TinySim",
#		phylo=rcoal(3),
#		root.seq=NucleotideSequence(string="ATG",processes=list(list(JC69())))
#	);
#	# run the simulation
#	Simulate(sim)
#	# get all the associated sequence objects
#	getSequences(sim)
#	# get the sequence associated with node 3
#	# via virtual field
#	sim$sequences[[3]]
# } 
# 
# @author 
# 
# \seealso{ 
# 	@seeclass 
# } 
# 
#*/###########################################################################
setMethodS3(
  "getSequences",
  class="PhyloSim",
  function(
    this,
    ...
  ){

		slist<-list();
		for (node in getNodes(this)){
			slist[[node]]<-getSeqFromNode(this, node=node);
		}
		return(slist);

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

##
## Method: setSequences
##
###########################################################################/**
#
# @RdocMethod setSequences
#
# @title "Forbidden action: setting the Sequence objects associated with the nodes of a phylo object aggregated by a PhyloSim 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(
  "setSequences",
  class="PhyloSim",
  function(
    		this,
		value,
    		...
  ){

		virtualAssignmentForbidden(this);

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

##
## Method: getAlignment
##
###########################################################################/**
#
# @RdocMethod	getAlignment
# 
# @title "Get the alignment stored in a PhyloSim object" 
# 
# \description{ 
#	@get "title".
# } 
# 
# @synopsis 
# 
# \arguments{ 
# 	\item{this}{A PhyloSim object.} 
# 	\item{...}{Not used.} 
# } 
# 
# \value{ 
# 	The alignment as a matrix. Gap are represented by strings composed of dashes.
# } 
# 
# \examples{
#	# Create a PhyloSim object.
#	# Provide the phylo object 
#	# and the root sequence.
#	sim<-PhyloSim(
#		name="TinySim",
#		phylo=rcoal(3),
#		root.seq=NucleotideSequence(string="ATG",processes=list(list(JC69())))
#	);
#	# run the simulation
#	Simulate(sim)
#	# get the resulting aligment
#	getAlignment(sim)
#	# via virtual fileld:
#	sim$alignment
# } 
# 
# @author 
# 
# \seealso{ 
# 	@seeclass 
# } 
# 
#*/###########################################################################
setMethodS3(
  "getAlignment",
  class="PhyloSim",
  function(
    this,
    ...
  ){

		this$.alignment;

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

##
## Method: setAlignment
##
###########################################################################/**
#
# @RdocMethod setAlignment
#
# @title "Forbidden action: setting the alignment stored in a PhyloSim 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(
  "setAlignment",
  class="PhyloSim",
  function(
    this,
    value,
    ...
  ){

               this$alignment <- value;

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

##
## Method: .recoverAlignment
##
setMethodS3(
  ".recoverAlignment",
  class="PhyloSim",
  function(
    		this,
		paranoid=FALSE,
    		...
  ){

		# Refuse to build alignment if at least one of the sequences is NA:
		for (seq in this$.sequences){
			if(!is.Sequence(seq)){
				throw("Cannot build alignment because the simulation is incomplete!\n");
			}
		}

		# The list holding all the partial alignment matrices:
		aln.mat<-list();
		
		# Assigning NA-s here to prevent creation of these variables in the global
		# environment.
		row.names<-NA;
		from.node<-NA;
		to.node<-NA;
		from.seq<-NA;
		to.seq<-NA;
		edge<-NA;
		from.name<-NA;
		to.name<-NA;
		from.mat<-NA;
		to.mat<-NA;

		# Initialize the variables:

		init.vars<-function(){

			# Getting the edge:
			edge<<-getEdge(this, edge.number);
			
			# Getting the nodes:
			from.node<<-edge[[1,"from"]];
			to.node<<-edge[[1,"to"]];

			# Getting the sequence objects:
			from.seq<<-getSeqFromNode(this, from.node)
			to.seq<<-getSeqFromNode(this, to.node)
	
			# Getting sequence names:	
			from.name<<-from.seq$name;
			to.name<<-to.seq$name;

		}

		# Initialize the aligment matrices:
		init.aln.mats…