PageRenderTime 169ms CodeModel.GetById 119ms app.highlight 31ms RepoModel.GetById 1ms app.codeStats 1ms

/PhyloSim.R

http://github.com/sbotond/phylosim
R | 6741 lines | 2944 code | 500 blank | 3297 comment | 491 complexity | 599a6b1f6c5fd46221d5d716d177dc82 MD5 | raw file

Large files files are truncated, but you can click here to view the full file

   1##	
   2## Copyright 2009 Botond Sipos	
   3## See the package description for licensing information.	
   4##	
   5##########################################################################/** 
   6#
   7# @RdocClass PhyloSim
   8# \alias{phylosim}
   9# 
  10# @title "The PhyloSim class"
  11# 
  12# \description{ 
  13#
  14#	PhyloSim is an extensible object-oriented framework for the Monte Carlo simulation 
  15#	of sequence evolution written in 100 percent \code{R}.
  16#	It is built on the top of the \code{\link{R.oo}} and \code{\link{ape}} packages and uses 
  17#	Gillespie's direct method to simulate substitutions, insertions and deletions.
  18#
  19#	Key features offered by the framework:	
  20#	\itemize{
  21#	\item  Simulation of the evolution of a set of discrete characters with arbitrary states evolving 
  22#		by a continuous-time Markov process with an arbitrary rate matrix.
  23#	\item Explicit implementations of the most popular substitution models (for nucleotides, amino acids and codons).
  24# 	\item Simulation under the popular models of among-sites rate variation, like the gamma (+G) and invariants plus gamma (+I+G) models.
  25#	\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.
  26#	\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).
  27#	\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.
  28#	\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.
  29#	\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.
  30#	\item The possibility to export the counts of various events ("branch statistics") as phylo objects (see \code{\link{exportStatTree.PhyloSim}}).
  31#	}
  32#
  33#	General notes:
  34#	\itemize{
  35#	\item The \code{Sequence} objects have no "immortal links". The simulation
  36#	is aborted if the sequence length shrinks to zero. It is up to the user 
  37#	to choose sensible indel rates and sequence lengths to prevent that.
  38#	\item The sites near the beginning and end of the sequences have less sites proposing
  39#	insertion and deletion events around the so the insertion and deletion processes 
  40#	have an "edge effect". The user can simulate
  41#	realistic flanking sequences to alleviate the edge effect in the simulation settings where
  42#	it may be an issue.
  43# }
  44#
  45#	Notes on performance: 
  46#	\itemize{
  47#	\item The pure \code{R} implementation offers felxibility, but also comes
  48#	with a slower simulation speed. If the \code{PSIM_FAST} object is present in the environment, a "fast & careless"
  49#	mode is enabled. In this mode most of the error checking is skipped, increasing the speed.
  50#	It is recomended that simulations are only run in fast mode if you are sure that the simulation
  51#	settings are free from errors. It is probably a good practice to set up the simulations in normal mode
  52#	with short sequences and enable fast mode when running the actual simulation with long sequences.
  53#	\item Please note, that no "branch statistics" are saved in fast mode.
  54#	\item Logging also has a negative impact on performance, so it's not a good idea to run
  55#	large simulations with the logging enabled.
  56#	\item The time needed to run a simulation depends not only on the number of the sites, 
  57#	but also on the length of the tree.
  58#	\item Constructing \code{Sequence} objects with large number of sites is expensive. Avoid doing
  59#	that inside a cycle.
  60#	\item In the case of \code{Sequence} objects with a large number of sites (more than 10 000) the 
  61#	amount of available memory can be limiting as well.
  62#  }
  63#
  64#	The examples below demonstrate only some more common simulation settings,
  65#	the framework offers much more flexibility. See the package
  66#	vignette (\code{vignette("PhyloSim",package="phylosim")}) and the
  67#	examples directory (\url{http://github.com/sbotond/phylosim/tree/master/examples/}) 
  68#	for additional examples.
  69#	
  70#	@classhierarchy
  71# }
  72#
  73#
  74# \references{
  75#	Gillespie, DT (1977) Exact stochastic simulation of coupled chemical reactions - 
  76#	J. Phys. Chem. 81 (25):2340-2361 \url{http://dx.doi.org/10.1021/j100540a008}
  77# }
  78#	
  79# @synopsis
  80#	
  81# \arguments{
  82# 	\item{phylo}{A rooted phylo object, constructed by the APE package.}
  83# 	\item{root.seq}{A valid Sequence object with Process objects attached. Used as the starting sequence during simulation.}
  84# 	\item{name}{The name of the object (a character vector of length one).}
  85# 	\item{log.file}{Name of the file used for logging.}
  86# 	\item{log.level}{An integer specifying the verbosity of logging (see \code{\link{setLogLevel.PhyloSim}}).}
  87# 	\item{...}{Not used.}
  88#	}
  89# 
  90# \section{Fields and Methods}{ 
  91# 	@allmethods
  92# }
  93# 
  94# \examples{ 
  95#   set.seed(1)
  96#	## The following examples demonstrate
  97#	## the typical use of the framework.
  98#	## See the package vignette and
  99#	## \url{http://github.com/sbotond/phylosim/tree/master/examples/}
 100#
 101#	## The ll() method gives information about the methods defined
 102#	## in the immediate class of an object.
 103#	## Useful when exploring the framework.
 104#
 105#	s<-Sequence()
 106#	ll(s)
 107#	ll(PhyloSim())
 108#	ll(GTR())
 109#
 110#	## Example 1 - A short simulation:
 111#	## simulate nucleotide seqeunces and display 
 112#	## the resulting alignment matrix.
 113#
 114#	Simulate(
 115#		PhyloSim(phy=rcoal(3),
 116#       root=NucleotideSequence(string="ATGC", proc=list(list(JC69())) ) )
 117#	)$alignment
 118#
 119#	# Construct a phylo object for the following
 120#	# simulations, scale total tree length to 1:
 121#
 122#	tmp<-PhyloSim(phylo=rcoal(3))
 123#	scaleTree(tmp,1/tmp$treeLength)
 124#	tmp$treeLength
 125#	t<-tmp$phylo
 126#
 127#	## Example 3 - simulating rate variation,
 128#	## insertions and deletions.
 129#	## See the examples/example_3_clean.R file
 130#	## in the phylosim GitHub repository.
 131#
 132#	# construct a GTR process object
 133#       gtr<-GTR(
 134#		name="MyGTR",
 135#               rate.params=list(
 136#                       "a"=1, "b"=2, "c"=3,
 137#                       "d"=1, "e"=2, "f"=3
 138#               ),
 139#               base.freqs=c(2,2,1,1)/6
 140#       )
 141#	# get object summary
 142#	summary(gtr)
 143#	# display a bubble plot
 144#	plot(gtr)
 145#
 146#	# construct root sequence object
 147#	s<-NucleotideSequence(length=20)
 148#	# attach process via virtual field
 149#	s$processes<-list(list(gtr))
 150#	# sample states from the equilibrium
 151#	# distribution of the attached processes
 152#	sampleStates(s)
 153#	# create among-site rate variation by sampling
 154#	# the "rate.multiplier" site-process-specific parameter
 155#	# from a discrete gamma distribution (GTR+G).
 156#	plusGamma(s,gtr,shape=0.1)
 157#	# make the range 11:12 invariable
 158#	setRateMultipliers(s,gtr,0,11:12)
 159#	# get the rate multipliers for s and gtr
 160#	getRateMultipliers(s,gtr)
 161#
 162	# construct a deletion process object
 163#	# proposing lengths in the range 1:3
 164#	d<-DiscreteDeletor(
 165#		rate=0.1,
 166#		name="MyDel",
 167#		sizes=c(1:3),
 168#		probs=c(3/6,2/6,1/6)
 169#	)
 170#	# get object 
 171#	summary(d)
 172#	# plot deletion length distribution
 173#	plot(d)
 174#	# attach deletion process d to sequence s
 175#	attachProcess(s,d)
 176# 	# create a region rejecting all deletions
 177#       setDeletionTolerance(s,d,0,11:12)
 178#
 179#	# construct an insertion process object
 180#	# proposing lengths in the range 1:3
 181#	i<-DiscreteInsertor(
 182#		rate=0.1,
 183#		name="MyDel",
 184#		sizes=c(1:2),
 185#		probs=c(1/2,1/2),
 186#		template.seq=NucleotideSequence(length=1,processes=list(list(JC69())))
 187#	) 
 188# 	# states will be sampled from the JC69 equilibrium distribution
 189#	# get object 
 190#	summary(i)
 191#	# plot insertion length distribution
 192#	plot(i)
 193#	# attach insertion process i to sequence s
 194#	attachProcess(s,i)
 195# 	# create a region rejecting all insertions
 196#       setInsertionTolerance(s,i,0,11:12)
 197#
 198#	# plot total site rates
 199#	plot(s)
 200#	# construct simulation object
 201#	sim<-PhyloSim(root.seq=s, phylo=t)
 202#	# get object summary
 203#	summary(sim)
 204#	# plot tree
 205#	plot(sim)
 206#	# run simulation
 207#	Simulate(sim)
 208# 	# get the list of recorded per-branch event counts
 209# 	getBranchEvents(sim)
 210#	# export the number of substitutions as a phylo object
 211#	subst<-exportStatTree(sim,"substitution")
 212#	# plot the exported phylo object
 213#	plot(subst)
 214#	# plot tree and alignment
 215#	plot(sim)
 216#	# save and display alingment
 217#	file<-paste("PhyloSim_dummy_fasta_",Sys.getpid(),".fas",sep="");
 218#	saveAlignment(sim,file=file);
 219#	# print out the Fasta file
 220#	cat(paste(scan(file=file,what=character(),sep="\n"),collapse="\n"));cat("\n");
 221#	# delete Fasta file
 222#	unlink(file);
 223#
 224#   # See \url{http://github.com/sbotond/phylosim/tree/master/examples/examples_class.R}
 225#   # for the full list of PhyloSim constructor examples.
 226#   
 227#	## See the package vignette and 
 228#	## the GitHub repository for even more examples.
 229# }
 230# 
 231# @author
 232#
 233# \seealso{ 
 234#	\code{\link{PSRoot} \link{Alphabet} \link{AminoAcidAlphabet} 
 235#	\link{AminoAcidSequence} \link{AminoAcidSubst}
 236#	\link{AnyAlphabet} \link{BinaryAlphabet} \link{BinarySequence} 
 237#	\link{BinarySubst} \link{BrownianInsertor} \link{CodonAlphabet}
 238#	\link{CodonSequence} \link{CodonUNREST} \link{ContinuousDeletor}
 239#	\link{ContinuousInsertor} \link{cpREV} \link{DiscreteDeletor}
 240#	\link{DiscreteInsertor} \link{Event} \link{F81} \link{F84}
 241#	\link{FastFieldDeletor} \link{GeneralDeletor}
 242#	\link{GeneralInDel} \link{GeneralInsertor} \link{GeneralSubstitution} 
 243#	\link{GTR} \link{GY94} \link{HKY} \link{JC69} \link{JTT} \link{JTT.dcmut}
 244#	\link{K80} \link{K81} \link{LG} \link{mtArt} \link{mtMam} \link{mtREV24}
 245#	\link{MtZoa} \link{NucleotideAlphabet} \link{NucleotideSequence} \link{PAM} 
 246#	\link{PAM.dcmut} \link{PhyloSim} \link{Process} \link{QMatrix} \link{Sequence}
 247#	\link{Site} \link{T92} \link{TN93} \link{UNREST} \link{WAG}}
 248# }
 249# 
 250#*/###########################################################################
 251setConstructorS3(
 252"PhyloSim",
 253  function(
 254		phylo=NA,
 255		root.seq=NA,
 256		name=NA,
 257		log.file=NA,
 258		log.level=-1, # no loggin is performed by default
 259		... 
 260		)	{
 261
 262		this<-PSRoot();	
 263		this<-extend(this,
 264			"PhyloSim",
 265			.name="Anonymous",
 266			.phylo=NA,
 267			.root.sequence=NA, 
 268			.sequences=list(),		# references to the sequence objects
 269			.node.hooks=list(),		# references to the node hook functions.
 270			.branch.stats=list(), # branch statistics.
 271			.alignment=NA,				# the resulting alignment in fasat format.
 272			.log.file=NA, 				# the name of the log file.
 273			.log.connection=NA,		# connection for the log file.
 274			.log.level=NA					# log level
 275		);
 276
 277		if(!all(is.na(phylo))){
 278			this$phylo<-phylo;
 279		}
 280
 281		if(!all(is.na(root.seq))){
 282			this$rootSeq<-root.seq;
 283		}
 284	
 285		if(!missing(name)){
 286			this$name<-name;
 287		}
 288
 289		if(!missing(log.file)){
 290			this$logFile<-log.file;
 291		} else {	
 292			# Setting default log file:
 293			tmp<-this$id;
 294			tmp<-gsub(":","_",tmp);
 295			this$logFile<-paste(tmp,".log",sep="");
 296		}
 297
 298		# Setting log level:
 299		this$logLevel<-log.level;
 300
 301		return(this);
 302
 303  },
 304  enforceRCC=TRUE
 305);
 306
 307##	
 308## Method: checkConsistency
 309##	
 310###########################################################################/**
 311#
 312# @RdocMethod	checkConsistency
 313# 
 314# @title "Check object consistency"
 315# 
 316# \description{ 
 317#		@get "title".
 318# } 
 319# 
 320# @synopsis 
 321#
 322# \arguments{ 
 323#       \item{this}{An object.} 
 324#       \item{...}{Not used.} 
 325# } 
 326# 
 327# 
 328# \value{ 
 329#		Returns an invisible TRUE if no inconsistencies found in the object, throws 
 330#		an error otherwise. 
 331# } 
 332# 
 333# @author 
 334# 
 335# \seealso{ 
 336#	@seeclass
 337# } 
 338# 
 339#*/###########################################################################
 340setMethodS3(
 341	"checkConsistency", 
 342	class="PhyloSim", 
 343	function(
 344		this,
 345		...
 346	){
 347
 348      may.fail<-function(this) {
 349					
 350				# Checking the name:	
 351				this$name<-this$name;
 352				# Checking the phylo object:
 353				if (!any(is.na(this$.phylo)) & !is.phylo(this$.phylo) ){
 354					throw("The phylo object is invalid!\n");
 355				}
 356				# Checking the log level:
 357				if(!is.numeric(this$.log.level) | (length(this$.log.level) != 1) ){
 358					throw("The log level must be numeric vector of length 1!\n");
 359				}
 360				# Checking lof file:
 361				if(!is.character(this$.log.file) | (length(this$.log.level) != 1) ){
 362					throw("The log file must be charcter vector of length 1!\n");
 363				}
 364				# Checking the sequences:
 365				for (seq in this$.sequences){
 366					if(is.Sequence(seq)){
 367						checkConsistency(seq);
 368					}
 369				}
 370				# Checking node hooks:
 371				for (hook in this$.node.hooks){
 372					if(!is.null(hook) & !is.function(hook)){
 373						throw("Invalid node hook found!\n");
 374					}
 375				}
 376				# Checking the alignment:
 377				if(!any(is.na(this$.alignment))){
 378					.checkAlignmentConsistency(this, this$.alignment);
 379				}
 380		
 381      }
 382      tryCatch(may.fail(this));
 383			return(invisible(TRUE));
 384
 385	},
 386	private=FALSE,
 387	protected=FALSE,
 388	overwrite=FALSE,
 389	conflict="warning",
 390	validators=getOption("R.methodsS3:validators:setMethodS3")
 391)
 392
 393##	
 394## Method: is.phylo.default
 395##	
 396###########################################################################/**
 397#
 398# @RdocDefault is.phylo
 399# 
 400# @title "Check if an object is an instance of the phylo class" 
 401# 
 402# \description{ 
 403#	@get "title".
 404#	Phylo objects are created by the \pkg{APE} package. This method just return the value of \code{inherits(this,"phylo")}.
 405# } 
 406# 
 407# @synopsis 
 408# 
 409# \arguments{ 
 410# 	\item{this}{An object.} 
 411# 	\item{...}{Not used.} 
 412# } 
 413# 
 414# \value{ 
 415# 	TRUE or FALSE.
 416# } 
 417# 
 418# \examples{
 419#	# load APE
 420#	library(ape);
 421#	# create some objects
 422#	o1<-Object();
 423#	o2<-rcoal(3);
 424#	# check if they are phylo objects
 425#	is.phylo(o1);
 426#	is.phylo(o2);			
 427#
 428# } 
 429# 
 430# @author 
 431# 
 432# \seealso{ 
 433# 	The \pkg{ape} package.
 434# } 
 435# 
 436#*/###########################################################################
 437setMethodS3(
 438	"is.phylo", 
 439	class="default", 
 440	function(
 441		this,
 442		...
 443	){
 444
 445		inherits(this,"phylo");
 446
 447	},
 448	private=FALSE,
 449	protected=FALSE,
 450	overwrite=FALSE,
 451	conflict="warning",
 452	validators=getOption("R.methodsS3:validators:setMethodS3")
 453)
 454
 455##	
 456## Method: setPhylo
 457##	
 458###########################################################################/**
 459#
 460# @RdocMethod setPhylo
 461# 
 462# @title "Set the phylo object for a PhyloSim object" 
 463# 
 464# \description{ 
 465#	@get "title".
 466#
 467#	The internal structure of the provided phylo object is reordered in a cladeweise fashion.
 468# } 
 469# 
 470# @synopsis 
 471# 
 472# \arguments{ 
 473# 	\item{this}{A PhyloSim object.} 
 474# 	\item{value}{A phylo object created by the \pkg{ape} package.} 
 475# 	\item{...}{Not used.} 
 476# } 
 477# 
 478# \value{ 
 479# 	A phylo object or FALSE.
 480# } 
 481# 
 482# \examples{
 483#	#create a PhyloSim object
 484#	sim<-PhyloSim();
 485#	# creat a phylo object
 486#	tree<-rcoal(3);
 487#	# get/set phylo object
 488#	setPhylo(sim,tree);
 489#	getPhylo(sim,tree);
 490#	# get/set phylo object via virtual field
 491#	sim$tree<-rcoal(5);
 492#	sim$tree;
 493# } 
 494# 
 495# @author 
 496# 
 497# \seealso{ 
 498# 	The PhyloSim class, the \pkg{ape} package.
 499# } 
 500# 
 501#*/###########################################################################
 502setMethodS3(
 503	"setPhylo", 
 504	class="PhyloSim", 
 505	function(
 506		this,
 507		value,
 508		...
 509	){
 510
 511		if(missing(value)){
 512			throw("No object provided!\n");
 513		}
 514		else if(!is.phylo(value)){
 515			throw("The new value must be a \"phylo\" object!\n");
 516		}
 517		else if(!is.rooted(value)){
 518			throw("The new value must be a rooted \"phylo\" object!\n");
 519		}
 520		else {
 521
 522			.checkTipLabels(value);
 523			this$.phylo<-value;
 524			this$.phylo<-reorder(this$.phylo, order="cladewise");
 525			for (i in this$nodes){
 526				this$.sequences[[i]]<-NA;
 527			}
 528			return(this$.phylo);
 529
 530		}
 531		return(FALSE);
 532
 533	},
 534	private=FALSE,
 535	protected=FALSE,
 536	overwrite=FALSE,
 537	conflict="warning",
 538	validators=getOption("R.methodsS3:validators:setMethodS3")
 539);
 540
 541##	
 542## Method: .checkTipLabels
 543##	
 544setMethodS3(
 545	".checkTipLabels", 
 546	class="phylo", 
 547	function(
 548		this,
 549		...
 550	){
 551
 552		for(label in this$tip.label){
 553			if(length(grep("^Node \\d+$",label,perl=TRUE,value=FALSE)) > 0){
 554					throw("Sorry, but the node labels matching \"Node \\d+\" are reserved for internal nodes! Blaming label: ",label,".\n");	
 555			}
 556			else if(length(grep("^Root node \\d+$",label,perl=TRUE,value=FALSE)) > 0){
 557					throw("Sorry, but the node labels matching \"Root node \\d+\" are reserved for the root node! Blaming label: ",label,".\n");	
 558			}
 559			
 560		}
 561
 562	},
 563	private=TRUE,
 564	protected=FALSE,
 565	overwrite=FALSE,
 566	conflict="warning",
 567	validators=getOption("R.methodsS3:validators:setMethodS3")
 568)
 569
 570##	
 571## Method: getPhylo
 572##	
 573###########################################################################/**
 574#
 575# @RdocMethod getPhylo
 576# 
 577# @title "Get the phylo object aggregated in a PhyloSim object" 
 578# 
 579# \description{ 
 580#	@get "title".
 581#
 582# } 
 583# 
 584# @synopsis 
 585# 
 586# \arguments{ 
 587# 	\item{this}{A PhyloSim object.} 
 588# 	\item{...}{Not used.} 
 589# } 
 590# 
 591# \value{ 
 592# 	A phylo object or NA.
 593# } 
 594# 
 595# \examples{
 596#	#create a PhyloSim object
 597#	sim<-PhyloSim();
 598#	# creat a phylo object
 599#	tree<-rcoal(3);
 600#	# get/set phylo object
 601#	setPhylo(sim,tree);
 602#	getPhylo(sim,tree);
 603#	# get/set phylo object via virtual field
 604#	sim$tree<-rcoal(5);
 605#	sim$tree;
 606# } 
 607# 
 608# @author 
 609# 
 610# \seealso{ 
 611# 	The PhyloSim class, the \pkg{ape} package.
 612# } 
 613# 
 614#*/###########################################################################
 615setMethodS3(
 616	"getPhylo", 
 617	class="PhyloSim", 
 618	function(
 619		this,
 620		...
 621	){
 622
 623		this$.phylo;
 624
 625	},
 626	private=FALSE,
 627	protected=FALSE,
 628	overwrite=FALSE,
 629	conflict="warning",
 630	validators=getOption("R.methodsS3:validators:setMethodS3")
 631)
 632
 633##	
 634## Method: setRootSeq
 635##	
 636###########################################################################/**
 637#
 638# @RdocMethod setRootSeq
 639# 
 640# @title "Set the root sequence for a PhyloSim object" 
 641# 
 642# \description{ 
 643#	@get "title".
 644#
 645#	The root sequence will be used as a starting point for the simulation. The phylo object must be set before
 646#	trying to set the root sequence object.
 647# } 
 648# 
 649# @synopsis 
 650# 
 651# \arguments{ 
 652# 	\item{this}{A PhyloSim object.} 
 653# 	\item{value}{A valid Sequence object.} 
 654# 	\item{...}{Not used.} 
 655# } 
 656# 
 657# \value{ 
 658# 	The root Sequence object if succesfull, FALSE otherwise.
 659# } 
 660# 
 661# \examples{
 662#	# create some objects
 663#	sim<-PhyloSim(phylo=rcoal(3));
 664#	seq<-NucleotideSequence(string="ATGCC");
 665#	# set/get root sequence
 666#	setRootSeq(sim, seq);
 667#	getRootSeq(sim, seq);
 668#	# set/get root sequence via virtual field
 669#	sim$rootSeq<-BinarySequence(string="111000111000");
 670#	sim$rootSeq;
 671#
 672# } 
 673# 
 674# @author 
 675# 
 676# \seealso{ 
 677# 	@seeclass Sequence Process
 678# } 
 679# 
 680#*/###########################################################################
 681setMethodS3(
 682	"setRootSeq", 
 683	class="PhyloSim", 
 684	function(
 685		this,
 686		value,
 687		...
 688	){
 689
 690		if(missing(value)){
 691			throw("No object provided!\n");
 692		}
 693		else if(!is.Sequence(value)){
 694			throw("The new value must be a sequence object!\n");
 695		}
 696		else {
 697
 698			this$.root.sequence<-clone(value);
 699			this$.root.sequence$name<-paste("Root node",this$rootNode);
 700        
 701            # Call garbage collection:
 702            gc();
 703            gc();
 704
 705			return(this$.root.sequence);
 706
 707		}
 708		return(FALSE);
 709
 710
 711	},
 712	private=FALSE,
 713	protected=FALSE,
 714	overwrite=FALSE,
 715	conflict="warning",
 716	validators=getOption("R.methodsS3:validators:setMethodS3")
 717)
 718
 719##	
 720## Method: getRootSeq
 721##	
 722###########################################################################/**
 723#
 724# @RdocMethod getRootSeq
 725# 
 726# @title "Get the root sequence aggregated by a PhyloSim object" 
 727# 
 728# \description{ 
 729#	@get "title".
 730#
 731# } 
 732# 
 733# @synopsis 
 734# 
 735# \arguments{ 
 736# 	\item{this}{A PhyloSim object.} 
 737# 	\item{...}{Not used.} 
 738# } 
 739# 
 740# \value{ 
 741# 	The root Sequence object or NA.
 742# } 
 743# 
 744# \examples{
 745#	# create some objects
 746#	sim<-PhyloSim(phylo=rcoal(3));
 747#	seq<-NucleotideSequence(string="ATGCC");
 748#	# set/get root sequence
 749#	setRootSeq(sim, seq);
 750#	getRootSeq(sim, seq);
 751#	# set/get root sequence via virtual field
 752#	sim$rootSeq<-BinarySequence(string="111000111000");
 753#	sim$rootSeq;
 754#
 755# } 
 756# 
 757# @author 
 758# 
 759# \seealso{ 
 760# 	@seeclass Sequence Process
 761# } 
 762# 
 763#*/###########################################################################
 764setMethodS3(
 765	"getRootSeq", 
 766	class="PhyloSim", 
 767	function(
 768		this,
 769		...
 770	){
 771
 772			this$.root.sequence;
 773
 774
 775	},
 776	private=FALSE,
 777	protected=FALSE,
 778	overwrite=FALSE,
 779	conflict="warning",
 780	validators=getOption("R.methodsS3:validators:setMethodS3")
 781)
 782
 783##	
 784## Method: as.character.PhyloSim
 785##	
 786###########################################################################/**
 787#
 788# @RdocMethod as.character
 789# 
 790# @title "Return the character representation of a PhyloSim object" 
 791# 
 792# \description{ 
 793#	@get "title".
 794#
 795#	The character representation is the identifier of the PhyloSim object as returned by the \code{getId} method.
 796# } 
 797# 
 798# @synopsis 
 799# 
 800# \arguments{ 
 801# 	\item{x}{A PhyloSim object.} 
 802# 	\item{...}{Not used.} 
 803# } 
 804# 
 805# \value{ 
 806# 	A character vector of length one.
 807# } 
 808# 
 809# \examples{
 810#	# create a PhyloSim object
 811#	o<-PhyloSim(name="MySim");
 812#	# get character representation
 813#	as.character(o)
 814# } 
 815# 
 816# @author 
 817# 
 818# \seealso{ 
 819# 	@seeclass 
 820# } 
 821# 
 822#*/###########################################################################
 823setMethodS3(
 824	"as.character", 
 825	class="PhyloSim", 
 826	function(
 827		x,
 828		...
 829	){
 830
 831		return(getId(x));
 832
 833	},
 834	private=FALSE,
 835	protected=FALSE,
 836	overwrite=FALSE,
 837	conflict="warning",
 838	validators=getOption("R.methodsS3:validators:setMethodS3")
 839)
 840
 841##
 842## Method: getName
 843##
 844###########################################################################/**
 845#
 846# @RdocMethod getName
 847# 
 848# @title "Get the name of a PhyloSim object" 
 849# 
 850# \description{ 
 851#	@get "title".
 852# } 
 853# 
 854# @synopsis 
 855# 
 856# \arguments{ 
 857# 	\item{this}{A PhyloSim object.} 
 858# 	\item{...}{Not used.} 
 859# } 
 860# 
 861# \value{ 
 862# 	A character vector of length one.
 863# } 
 864# 
 865# \examples{
 866#	# create a PhyloSim object
 867#	o<-PhyloSim();
 868#	# set/get name
 869#	setName(o,"MySim");
 870#	getName(o,"MySim");
 871#	# set/get name via virtual field
 872#	o$name<-"George";
 873#	o$name
 874# } 
 875# 
 876# @author 
 877# 
 878# \seealso{ 
 879# 	@seeclass 
 880# } 
 881# 
 882#*/###########################################################################
 883setMethodS3(
 884  "getName",
 885  class="PhyloSim",
 886  function(
 887    this,
 888    ...
 889  ){
 890
 891    this$.name;
 892  },
 893  private=FALSE,
 894  protected=FALSE,
 895  overwrite=FALSE,
 896  conflict="warning",
 897  validators=getOption("R.methodsS3:validators:setMethodS3")
 898);
 899
 900##
 901## Method: setName
 902##
 903###########################################################################/**
 904#
 905# @RdocMethod setName
 906# 
 907# @title "Set the name of a PhyloSim object" 
 908# 
 909# \description{ 
 910#	@get "title".
 911# } 
 912# 
 913# @synopsis 
 914# 
 915# \arguments{ 
 916# 	\item{this}{A PhyloSim object.} 
 917#	\item{new.name}{A character vector of length one.}
 918# 	\item{...}{Not used.} 
 919# } 
 920# 
 921# \value{ 
 922# 	The new name.
 923# } 
 924# 
 925# \examples{
 926#	# create a PhyloSim object
 927#	o<-PhyloSim();
 928#	# set/get name
 929#	setName(o,"MySim");
 930#	getName(o,"MySim");
 931#	# set/get name via virtual field
 932#	o$name<-"George";
 933#	o$name
 934# } 
 935# 
 936# @author 
 937# 
 938# \seealso{ 
 939# 	@seeclass 
 940# } 
 941# 
 942#*/###########################################################################
 943setMethodS3(
 944  "setName",
 945  class="PhyloSim",
 946  function(
 947    this,
 948    new.name,
 949    ...
 950  ){
 951
 952    this$.name<-as.character(new.name);
 953  },
 954  private=FALSE,
 955  protected=FALSE,
 956  overwrite=FALSE,
 957  conflict="warning",
 958  validators=getOption("R.methodsS3:validators:setMethodS3")
 959);
 960
 961##
 962## Method: getId
 963##
 964###########################################################################/**
 965#
 966# @RdocMethod getId
 967# 
 968# @title "Get the unique identifier of a PhyloSim object" 
 969# 
 970# \description{ 
 971#	@get "title".
 972#	The unique identifier is the concatenation of the class, the object name as returned by getName() and the object hash 
 973#       as returned by hashCode().
 974# } 
 975# 
 976# @synopsis 
 977# 
 978# \arguments{ 
 979# 	\item{this}{A PhyloSim object.}
 980# 	\item{...}{Not used.} 
 981# } 
 982# 
 983# \value{ 
 984# 	A character vector of length one.
 985# } 
 986# 
 987# \examples{
 988#	# create a PhyloSim object
 989#	o<-PhyloSim(name="MySim");
 990#	# get id
 991#	getId(o);
 992#	# get id via virtual field
 993#	o$id;
 994# } 
 995# 
 996# @author 
 997# 
 998# \seealso{ 
 999# 	@seeclass 
1000# } 
1001# 
1002#*/###########################################################################
1003setMethodS3(
1004  "getId",
1005  class="PhyloSim",
1006  function(
1007    this,
1008    ...
1009  ){
1010
1011  	this.class<-class(this)[1];
1012	id<-paste(this.class,this$.name,hashCode(this),sep=":");
1013	return(id);
1014
1015  },
1016  private=FALSE,
1017  protected=FALSE,
1018  overwrite=FALSE,
1019  conflict="warning",
1020  validators=getOption("R.methodsS3:validators:setMethodS3")
1021);
1022
1023##
1024## Method: setId
1025##
1026###########################################################################/**
1027#
1028# @RdocMethod setId
1029#
1030# @title "Forbidden action: setting the unique identifier of a PhyloSim object"
1031#
1032# \description{
1033#       @get "title".
1034#
1035# }
1036#
1037# @synopsis
1038#
1039# \arguments{
1040#       \item{this}{An object.}
1041#       \item{value}{Not used.}
1042#       \item{...}{Not used.}
1043# }
1044#
1045# \value{
1046#	Throws an error.
1047# }
1048#
1049# @author
1050#
1051# \seealso{
1052#       @seeclass
1053# }
1054#
1055#*/###########################################################################
1056setMethodS3(
1057  "setId",
1058  class="PhyloSim",
1059  function(
1060    this,
1061    value,
1062    ...
1063  ){
1064
1065  throw("Id is generated automatically and it cannot be set!\n");
1066
1067  },
1068  private=FALSE,
1069  protected=FALSE,
1070  overwrite=FALSE,
1071  conflict="warning",
1072  validators=getOption("R.methodsS3:validators:setMethodS3")
1073);
1074
1075##
1076## Method: Simulate
1077##
1078###########################################################################/**
1079#
1080# @RdocMethod Simulate
1081# 
1082# @title "Run a simulation according to a PhyloSim object" 
1083# 
1084# \description{ 
1085#	@get "title".
1086#
1087#	The phylo object and the root sequence must be set before attempting to run a simulation.
1088#	Also the bigRate of the root sequence must not be NA or zero, so at least one sane
1089#	Process object must be attached to the root sequence object.
1090# } 
1091# 
1092# @synopsis 
1093# 
1094# \arguments{ 
1095# 	\item{this}{A PhyloSim object.} 
1096# 	\item{quiet}{TRUE or FALSE (default).} 
1097# 	\item{...}{Not used.} 
1098# } 
1099# 
1100# \value{ 
1101# 	The PhyloSim object (invisible).
1102# } 
1103# 
1104# \examples{
1105#	# Create a PhyloSim object.
1106#	# Provide the phylo object 
1107#	# and the root sequence.
1108#	sim<-PhyloSim(
1109#		name="TinySim",
1110#		phylo=rcoal(3),
1111#		root.seq=NucleotideSequence(string="ATGC",processes=list(list(JC69())))
1112#	);
1113#	# Run the simulation
1114#	Simulate(sim);
1115#	# Print the resulting sequences
1116#	sim$sequences
1117#	# Print the resulting alignment
1118#	sim$alignment
1119# } 
1120# 
1121# @author 
1122# 
1123# \seealso{ 
1124# 	@seeclass 
1125# } 
1126# 
1127#*/###########################################################################
1128setMethodS3(
1129  "Simulate",
1130  class="PhyloSim",
1131  function(
1132    		this,
1133		quiet=FALSE,
1134    ...
1135  ){
1136
1137
1138		if(!is.phylo(this$.phylo)){
1139			throw("Cannot simulate because the phylo object is not set or it is invalid!\n");
1140		}
1141		# Check for the root sequence:
1142		else if(!is.Sequence(this$.root.sequence)){
1143			throw("Cannot simulate because the root sequence is not set or it is invalid!\n");
1144		}
1145		# Check bigRate validity:
1146		else if(is.na(this$.root.sequence$bigRate)){
1147			throw("Cannot simulate because the bigRate of the root sequence is NA!\n");
1148		}
1149		else{
1150
1151			# Warn for zero bigRate:
1152			if(this$.root.sequence$bigRate == 0){
1153				warning("The bigRate of the root sequence is zero! You are running a pointless simulation!\n");
1154			}
1155
1156			if(exists(x="PSIM_FAST")){
1157				if(!quiet){
1158cat("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
1159cat("!! WARNING: fast & careless mode is on, most of the error checking is omitted!  !!\n");
1160cat("!!    Please note that this also disables the saving of branch statistics.      !!\n");
1161cat("!!       You can go back to normal mode by deleting the PSIM_FAST object.       !!\n");
1162cat("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
1163				}
1164			Log(this,"WARNING: fast & careless mode is on, most of the error checking is omitted!");
1165			}
1166		
1167			# Attach root sequence to root node:
1168			Log(this,paste("Attaching root sequence ",this$.root.sequence$id,sep=""));
1169			attachSeqToNode(this, node=getRootNode(this),seq=this$.root.sequence);
1170
1171			# Write protecting the root sequence:
1172			Log(this,paste("Write protecting root sequence ",this$.root.sequence$id,sep=""));
1173			this$.root.sequence$writeProtected<-TRUE;
1174
1175			# Traverse the tree and simulate:
1176			Log(this,paste("Starting simulation on the object",this$id));	
1177			edge.counter<-1;
1178			n.edges<-this$nedges;
1179			for(edge in 1:n.edges){
1180				if(!quiet){ cat("Simulating edge",edge,"of", n.edges,"\n");}
1181				Log(this,paste("Starting to simulate edge",edge,"of",n.edges));	
1182				.simulateEdge(this,number=edge);
1183				edge.counter<-edge.counter+1;
1184			}
1185		}
1186		Log(this, "Simulation finished, building alignment!\n");
1187		this$.alignment<-.recoverAlignment(this);
1188		# Flush the log connection:
1189		if(!is.na(this$.log.connection)){
1190				close(this$.log.connection);
1191		}
1192
1193        # Call the garbage collector:
1194        gc();
1195        gc();
1196
1197		return(invisible(this));
1198  },
1199  private=FALSE,
1200  protected=FALSE,
1201  overwrite=FALSE,
1202  conflict="warning",
1203  validators=getOption("R.methodsS3:validators:setMethodS3")
1204);
1205
1206##
1207## Method: .simulateEdge
1208##
1209setMethodS3(
1210  ".simulateEdge",
1211  class="PhyloSim",
1212  function(
1213    this,
1214		number=NA,
1215    ...
1216  ){
1217
1218		# Get edge:
1219		edge<-getEdge(this, number);
1220		# Get parent node:
1221		start.seq<-getSeqFromNode(this, edge[[1,"from"]]);
1222		# Evolve sequence:
1223		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);
1224		# Write protect the sequence:
1225		new.seq$writeProtected<-TRUE;
1226		# Attach sequence to children node:
1227		attachSeqToNode(this, node=edge[1,"to"], seq=new.seq);
1228
1229  },
1230  private=TRUE,
1231  protected=FALSE,
1232  overwrite=FALSE,
1233  conflict="warning",
1234  validators=getOption("R.methodsS3:validators:setMethodS3")
1235);
1236
1237##
1238## Method: .evolveBranch
1239##
1240setMethodS3(
1241  ".evolveBranch",
1242  class="PhyloSim",
1243  function(
1244    this,
1245    start.seq=NA,
1246    branch.length=NA,
1247		old.node=NA,
1248		new.node=NA,
1249		branch.number=NA,
1250    ...
1251  ){
1252
1253	if(!exists(x="PSIM_FAST")){
1254
1255		if(missing(start.seq)){
1256			throw("No starting sequence provided!\n");
1257		}
1258		else if(missing(branch.length)){
1259			throw("No branch length provided!\n");
1260		}
1261		else if(!is.numeric(branch.length)){
1262			throw("The branch length must be numeric!\n");
1263		}
1264	}
1265		if(.checkSeq(this, start.seq) ){
1266		
1267			# Cloning the starting sequence:
1268			seq<-clone(start.seq);
1269			
1270			# Set the name of the sequence object:
1271			if(is.tip(this, new.node)){
1272				seq$name<-this$tipLabels[[new.node]];
1273			}
1274			else {
1275				seq$name<-paste("Node",new.node);
1276			}
1277			.GillespieDirect(this, seq=seq, branch.length=branch.length, branch.number=branch.number);
1278			
1279			# Call the node hook if exists:
1280			hook<-this$.node.hooks[[as.character(new.node)]];
1281			if(!is.null(hook) & is.function(hook)){
1282				Log(this,paste("Calling node hook for node",new.node));
1283				seq<-hook(seq=seq);	
1284				if(!is.Sequence(seq)){
1285					throw("Node hook returned an invalid sequence object!\n");
1286				}
1287				else if(is.na(seq$bigRate)){
1288					throw("Node hook returned sequence with NA bigRate!\n");
1289				}
1290				else if(seq$bigRate == 0.0){
1291					throw("Node hook returned sequence with zero bigRate!\n");
1292				}
1293				else{
1294				 checkConsistency(seq, omit.sites=TRUE);
1295				}
1296			}
1297
1298			# Return the resulting sequence object:
1299			return(seq);
1300		
1301		}
1302  },
1303  private=FALSE,
1304  protected=FALSE,
1305  overwrite=FALSE,
1306  conflict="warning",
1307  validators=getOption("R.methodsS3:validators:setMethodS3")
1308);
1309
1310##
1311## Method: .GillespieDirect
1312##
1313setMethodS3(
1314  ".GillespieDirect",
1315  class="PhyloSim",
1316  function(
1317    this,
1318    seq=NA,
1319    branch.length=NA,
1320		branch.number=NA,
1321    ...
1322  ){
1323
1324		Debug(this, paste("Branch length is",branch.length));
1325		
1326		# Initialize time:
1327
1328		time<-0.0;
1329		
1330		# Sample the next waiting time until
1331		# the branch length is consumed:	
1332
1333		while( (time<-time + rexp(n=1, rate=(big.rate<-getBigRate(seq)))) <= branch.length){
1334
1335			# Generate a random number between zero and the bigRate:
1336
1337			E<-runif(n=1,min=0,max=big.rate);
1338
1339			# Identify the target site:
1340
1341			site.number<-which(.getCumulativeRatesFast(seq) >= E)[[1]];
1342
1343			# Get the events from the target site:
1344
1345			site<-seq$.sites[[site.number]];
1346			site$.position<-site.number;
1347			events<-getEvents(site);
1348			site$.position<-NULL;
1349			
1350			# Get the rates:
1351			rates<-double();
1352			for(e in events){
1353				rates<-c(rates,e$.rate);
1354			}
1355		
1356			# Calculate the corresponding cumulative rates:	
1357			if(site.number > 1){
1358				rates<-cumsum(c(seq$.cumulative.rates[[site.number - 1]], rates));
1359			}
1360			else {
1361				rates<-cumsum(c(0.0, rates));
1362			}
1363
1364			# Pick the event:
1365
1366			event.number<-which(rates >= E)[[1]] - 1;
1367			event<-events[[event.number]];
1368
1369			# Log the event:
1370
1371			Log(this,paste("Performing event [",event$.name,"] at position",event$.position,"generated by the process",event$.process$.id));
1372
1373			# Perform the event:
1374
1375			event.details<-Perform(event);
1376			Debug(this,paste("Remaining branch length is",(branch.length-time) ));
1377			
1378			# Log event details:
1379	
1380					# Log deletion event details:
1381						if(event$.name == "Deletion"){
1382						Log(this,paste("The process",event$.process,"proposed to delete range",paste(event.details$range,collapse="--"),". Accepted:",event.details$accepted));
1383					}
1384					# Log insertion event details:
1385					else if(event$.name == "Insertion"){
1386						message<-paste("The process ",event$.process," proposed insertion at position ",event.details$position,". Accepted: ",event.details$accepted,sep="");
1387						if(event.details$accepted == TRUE){
1388						message<-paste(message,"."," Insert length was ",event.details$length,sep="");
1389						}
1390						Log(this, message);
1391					}
1392
1393			# Update branch statistics if not in fast mode:
1394			if(!exists(x="PSIM_FAST")){
1395				.UpdateBranchStats(this,event,event.details, branch.number);
1396			}
1397
1398			# Abort if sequence length shrunk to zero:
1399
1400			if(seq$.length == 0){
1401				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");
1402				Log(this, message);
1403				throw(message);
1404			}	
1405
1406		} #/while
1407
1408        # Calling garbage collection:
1409        gc();
1410        gc();
1411		return(seq);
1412
1413  },
1414  private=FALSE,
1415  protected=FALSE,
1416  overwrite=FALSE,
1417  conflict="warning",
1418  validators=getOption("R.methodsS3:validators:setMethodS3")
1419);
1420
1421##
1422## Method: attachSeqToNode
1423##
1424###########################################################################/**
1425#
1426# @RdocMethod	attachSeqToNode
1427# 
1428# @title "Assotiate a Sequence object with a given node of a phylo object aggregated by a PhyloSim object" 
1429# 
1430# \description{ 
1431#	@get "title".
1432#
1433#	This method is mainly used internally.
1434#	
1435# } 
1436# 
1437# @synopsis 
1438# 
1439# \arguments{ 
1440# 	\item{this}{A PhyloSim object.} 
1441#	\item{node}{Node identifier.}
1442#	\item{seq}{A Sequence object.}
1443# 	\item{...}{Not used.} 
1444# } 
1445# 
1446# \value{ 
1447# 	The PhyloSim object (invisible).
1448# } 
1449# 
1450# 
1451# @author 
1452# 
1453# \seealso{ 
1454# 	@seeclass 
1455# } 
1456# 
1457#*/###########################################################################
1458setMethodS3(
1459  "attachSeqToNode",
1460  class="PhyloSim",
1461  function(
1462		this,
1463		node=NA,
1464		seq=NA,
1465    ...
1466  ){
1467
1468		if(!is.phylo(this$.phylo)){
1469			throw("The phylo object is not set, sequence to node is not possible!\n");
1470		}
1471		if(missing(node)){
1472			throw("No node specified!\n");
1473		}
1474		else if(missing(seq)){
1475			throw("No sequence object given");
1476		}
1477		else if(.checkNode(this,node) & .checkSeq(this, seq)){
1478			
1479			if(is.Sequence(this$.sequences[[node]])){
1480				throw("The node has already an attached sequence. Detach that before trying to attach a new one!\n");
1481			}
1482			else {
1483				this$.sequences[[as.numeric(node)]]<-seq;
1484				return(invisible(this));
1485			}
1486
1487		}
1488	
1489		return(invisible(this));
1490  },
1491  private=FALSE,
1492  protected=FALSE,
1493  overwrite=FALSE,
1494  conflict="warning",
1495  validators=getOption("R.methodsS3:validators:setMethodS3")
1496);
1497
1498##
1499## Method: attachHookToNode
1500##
1501###########################################################################/**
1502#
1503# @RdocMethod attachHookToNode
1504# 
1505# @title "Attach a callback function to a given node of a phylo object aggregated by a PhyloSim object" 
1506# 
1507# \description{ 
1508#	@get "title".
1509#
1510#	A "node hook" is a function which accepts a Sequence object through the named argument "seq" and returns a 
1511#	Sequence object. The node hook function must accept any object which inherits from the \code{Sequence} class!
1512#
1513#	After simulating the branch leading to the node, the resulting Sequence object is passed
1514#	to the node hook and the returned object is used to simulate the downstream branches.
1515#
1516#	By using node hooks the attached processes can be replaced during simulation, hence enabling the simulation of 
1517#	non-homogeneous sequence evolution.
1518# } 
1519# 
1520# @synopsis 
1521# 
1522# \arguments{ 
1523# 	\item{this}{A PhyloSim object.} 
1524# 	\item{node}{Node identifier.} 
1525# 	\item{fun}{A function (see above).} 
1526# 	\item{...}{Not used.} 
1527# } 
1528# 
1529# \value{ 
1530# 	The PhyloSim object (invisible).
1531# } 
1532# 
1533# \examples{
1534#	# Create a PhyloSim object.
1535#	# Provide the phylo object 
1536#	# and the root sequence.
1537#	sim<-PhyloSim(
1538#		name="TinySim",
1539#		phylo=rcoal(3),
1540#		root.seq=NucleotideSequence(string="ATGC",processes=list(list(JC69())))
1541#	);
1542#	# create a node hook function
1543#	hook<-function(seq=NA){
1544#		# replace the substitution process with F84
1545#		if(inherits(seq,"NucleotideSequence")){
1546#			cat("Replacing JC69 with F84.\n");
1547#			seq$processes<-list(list(F84(rate.params=list("Kappa" = 2)))); 
1548#		}
1549#		return(seq);
1550#	}
1551#	# attach hook function to node 5
1552#	attachHookToNode(sim,5,hook);
1553#	# Run the simulation
1554#	Simulate(sim);
1555#	# Check if the processes have been truly replaced
1556#	lapply(sim$sequences, getUniqueProcesses.Sequence)
1557#	# Print the resulting alignment
1558#	sim$alignment
1559# } 
1560# 
1561# @author 
1562# 
1563# \seealso{ 
1564# 	@seeclass 
1565# } 
1566# 
1567#*/###########################################################################
1568setMethodS3(
1569  "attachHookToNode",
1570  class="PhyloSim",
1571  function(
1572    		this,
1573		node=NA,
1574		fun=NA,
1575    ...
1576  ){
1577
1578		if(!is.phylo(this$.phylo)){
1579			throw("The phylo object is not set, attaching node hook is not possible!\n");
1580		}
1581		if(missing(node)){
1582			throw("No node specified!\n");
1583		}
1584		else if(missing(fun)){
1585			throw("No function given!");
1586		}
1587		else if(!is.function(fun)){
1588			throw("The argument \"fun\" must be a function!\n");
1589		}
1590		else if( length(intersect(names(formals(fun)), "seq")) == 0 ){
1591			throw("The function argument must have a an argument named \"seq\"");
1592		}
1593		else if(!is.Sequence(fun(Sequence(length=1)))){
1594      throw("The insert hook function must return a Sequence object!\n");
1595		}
1596		else if( .checkNode(this,node) ){
1597			if(is.function(this$.node.hooks[[as.character(node)]])){
1598				throw("The node has already an attached node hook. Detach that before trying to attach a new one!\n");
1599			}
1600			else {
1601				this$.node.hooks[[as.character(node)]]<-fun;
1602				return(invisible(this));
1603			}
1604
1605		}
1606	
1607		return(invisible(this));
1608  },
1609  private=FALSE,
1610  protected=FALSE,
1611  overwrite=FALSE,
1612  conflict="warning",
1613  validators=getOption("R.methodsS3:validators:setMethodS3")
1614);
1615
1616
1617##
1618## Method: .checkNode
1619##
1620setMethodS3(
1621  ".checkNode",
1622  class="PhyloSim",
1623  function(
1624    this,
1625		node=NA,
1626    ...
1627  ){
1628
1629		if(missing(node)){
1630			throw("No node specified!\n");
1631		} else if( length(intersect(node, getNodes(this))) != 1){
1632			throw("The specified node is invalid!\n");	
1633		}
1634		else {
1635			return(TRUE);
1636		}
1637  },
1638  private=FALSE,
1639  protected=FALSE,
1640  overwrite=FALSE,
1641  conflict="warning",
1642  validators=getOption("R.methodsS3:validators:setMethodS3")
1643);
1644
1645##
1646## Method: .checkSeq
1647##
1648setMethodS3(
1649  ".checkSeq",
1650  class="PhyloSim",
1651  function(
1652    this,
1653		seq=NA,
1654    ...
1655  ){
1656
1657		if(missing(seq)){
1658			throw("No sequence specified!\n");
1659		} else if(!is.Sequence(seq)){
1660			throw("The sequence object is invalid!\n");	
1661		}
1662		else {
1663			return(TRUE);
1664		}
1665  },
1666  private=FALSE,
1667  protected=FALSE,
1668  overwrite=FALSE,
1669  conflict="warning",
1670  validators=getOption("R.methodsS3:validators:setMethodS3")
1671);
1672
1673##
1674## Method: detachSeqFromNode
1675##
1676###########################################################################/**
1677#
1678# @RdocMethod	detachSeqFromNode
1679# 
1680# @title "Detach a Sequence object from a given node of a phylo object aggregated by a PhyloSim object" 
1681# 
1682# \description{ 
1683#	@get "title".
1684#
1685#	This method is mainly used internally.
1686#	
1687# } 
1688# 
1689# @synopsis 
1690# 
1691# \arguments{ 
1692# 	\item{this}{A PhyloSim object.} 
1693#	\item{node}{Node identifier.}
1694# 	\item{...}{Not used.} 
1695# } 
1696# 
1697# \value{ 
1698# 	The PhyloSim object (invisible).
1699# } 
1700# 
1701# @author 
1702# 
1703# \seealso{ 
1704# 	@seeclass 
1705# } 
1706# 
1707#*/###########################################################################
1708setMethodS3(
1709  "detachSeqFromNode",
1710  class="PhyloSim",
1711  function(
1712    		this,
1713		node=NA,
1714    		...
1715  ){
1716
1717		if(missing(node)){
1718			throw("No node specified!\n");
1719		}
1720		else if( .checkNode(this,node) ){
1721			
1722				this$.sequences[[as.numeric(node)]]<-NA;
1723		}
1724	
1725		return(invisible(this));
1726  },
1727  private=FALSE,
1728  protected=FALSE,
1729  overwrite=FALSE,
1730  conflict="warning",
1731  validators=getOption("R.methodsS3:validators:setMethodS3")
1732);
1733
1734##
1735## Method: detachHookFromNode
1736##
1737###########################################################################/**
1738#
1739# @RdocMethod	detachHookFromNode
1740# 
1741# @title "Detach a node hook function from a given node of a phylo object aggregated by a PhyloSim object" 
1742# 
1743# \description{ 
1744#	@get "title".
1745#
1746# } 
1747# 
1748# @synopsis 
1749# 
1750# \arguments{ 
1751# 	\item{this}{A PhyloSim object.} 
1752#	\item{node}{Node identifier.}
1753# 	\item{...}{Not used.} 
1754# } 
1755# 
1756# \value{ 
1757# 	The PhyloSim object (invisible).
1758# } 
1759# 
1760# \examples{
1761#	# Create a PhyloSim object.
1762#	# Provide the phylo object 
1763#	# and the root sequence.
1764#	sim<-PhyloSim(
1765#		name="TinySim",
1766#		phylo=rcoal(3),
1767#		root.seq=NucleotideSequence(string="ATGC",processes=list(list(JC69())))
1768#	);
1769#	# create a node hook function
1770#	hook<-function(seq=NA){
1771#		# replace the substitution process with F84
1772#		if(inherits(seq,"NucleotideSequence")){
1773#			cat("Replacing JC69 with F84.\n");
1774#			seq$processes<-list(list(F84(rate.params=list("Kappa" = 2)))); 
1775#		}
1776#		return(seq);
1777#	}
1778#	# attach hook function to node 5
1779#	attachHookToNode(sim,5,hook);
1780#	# detach hook from node 5
1781#	detachHookFromNode(sim,5);
1782#	# Run the simulation again
1783#	Simulate(sim);	# You should not see the message printed out by the "hook" function.
1784#
1785# } 
1786# 
1787# @author 
1788# 
1789# \seealso{ 
1790# 	attachHookToNode PhyloSim Simulate.PhyloSim
1791# } 
1792# 
1793#*/###########################################################################
1794setMethodS3(
1795  "detachHookFromNode",
1796  class="PhyloSim",
1797  function(
1798    this,
1799		node=NA,
1800    ...
1801  ){
1802
1803		if(missing(node)){
1804			throw("No node specified!\n");
1805		}
1806		else if( .checkNode(this,node) ){
1807			
1808				this$.node.hooks[[as.character(node)]]<-NA;
1809		}
1810	
1811		return(invisible(this));
1812  },
1813  private=FALSE,
1814  protected=FALSE,
1815  overwrite=FALSE,
1816  conflict="warning",
1817  validators=getOption("R.methodsS3:validators:setMethodS3")
1818);
1819
1820
1821##
1822## Method: getSeqFromNode
1823##
1824###########################################################################/**
1825#
1826# @RdocMethod getSeqFromNode
1827# 
1828# @title "Get the Sequence object associated with a given node of a phylo object aggregated by a PhyloSim object" 
1829# 
1830# \description{ 
1831#	@get "title".
1832# } 
1833# 
1834# @synopsis 
1835# 
1836# \arguments{ 
1837# 	\item{this}{A PhyloSim object.} 
1838# 	\item{node}{Node identifier.} 
1839# 	\item{...}{Not used.} 
1840# } 
1841# 
1842# \value{ 
1843# 	A Sequence object.
1844# } 
1845# 
1846# \examples{
1847#	# Create a PhyloSim object.
1848#	# Provide the phylo object 
1849#	# and the root sequence.
1850#	sim<-PhyloSim(
1851#		name="TinySim",
1852#		phylo=rcoal(3),
1853#		root.seq=NucleotideSequence(string="ATG",processes=list(list(JC69())))
1854#	);
1855#	# get the sequence associated with node 5
1856#	getSeqFromNode(sim,5)	# Should be NA
1857#	# Run the simulation
1858#	Simulate(sim)
1859#	# try again
1860#	getSeqFromNode(sim,5)
1861# } 
1862# 
1863# @author 
1864# 
1865# \seealso{ 
1866# 	@seeclass 
1867# } 
1868# 
1869#*/###########################################################################
1870setMethodS3(
1871  "getSeqFromNode",
1872  class="PhyloSim",
1873  function(
1874    		this,
1875		node=NA,
1876    		...
1877  ){
1878
1879		if(missing(node)){
1880			throw("No node specified!\n");
1881		}
1882		else if( .checkNode(this,node) ){
1883			
1884				return(this$.sequences[[as.numeric(node)]]);
1885		}
1886
1887  },
1888  private=FALSE,
1889  protected=FALSE,
1890  overwrite=FALSE,
1891  conflict="warning",
1892  validators=getOption("R.methodsS3:validators:setMethodS3")
1893);
1894
1895##
1896## Method: getSequences
1897##
1898###########################################################################/**
1899#
1900# @RdocMethod getSequences
1901# 
1902# @title "Gets all the Sequence objects associated with the nodes of a phylo object aggregated by a PhyloSim object" 
1903# 
1904# \description{ 
1905#	@get "title".
1906#
1907#	The order of the Sequence objects in the returned list reflects the identifiers of the associated nodes.
1908# } 
1909# 
1910# @synopsis 
1911# 
1912# \arguments{ 
1913# 	\item{this}{A PhyloSim object.} 
1914# 	\item{...}{Not used.} 
1915# } 
1916# 
1917# \value{ 
1918# 	A list of sequence objects.
1919# } 
1920# 
1921# \examples{
1922#	# Create a PhyloSim object.
1923#	# Provide the phylo object 
1924#	# and the root sequence.
1925#	sim<-PhyloSim(
1926#		name="TinySim",
1927#		phylo=rcoal(3),
1928#		root.seq=NucleotideSequence(string="ATG",processes=list(list(JC69())))
1929#	);
1930#	# run the simulation
1931#	Simulate(sim)
1932#	# get all the associated sequence objects
1933#	getSequences(sim)
1934#	# get the sequence associated with node 3
1935#	# via virtual field
1936#	sim$sequences[[3]]
1937# } 
1938# 
1939# @author 
1940# 
1941# \seealso{ 
1942# 	@seeclass 
1943# } 
1944# 
1945#*/###########################################################################
1946setMethodS3(
1947  "getSequences",
1948  class="PhyloSim",
1949  function(
1950    this,
1951    ...
1952  ){
1953
1954		slist<-list();
1955		for (node in getNodes(this)){
1956			slist[[node]]<-getSeqFromNode(this, node=node);
1957		}
1958		return(slist);
1959
1960  },
1961  private=FALSE,
1962  protected=FALSE,
1963  overwrite=FALSE,
1964  conflict="warning",
1965  validators=getOption("R.methodsS3:validators:setMethodS3")
1966);
1967
1968##
1969## Method: setSequences
1970##
1971###########################################################################/**
1972#
1973# @RdocMethod setSequences
1974#
1975# @title "Forbidden action: setting the Sequence objects associated with the nodes of a phylo object aggregated by a PhyloSim object"
1976#
1977# \description{
1978#       @get "title".
1979# }
1980#
1981# @synopsis
1982#
1983# \arguments{
1984#       \item{this}{An object.}
1985#       \item{value}{Not used.}
1986#       \item{...}{Not used.}
1987# }
1988#
1989# \value{
1990#	Throws an error.
1991# }
1992#
1993# @author
1994#
1995# \seealso{
1996#       @seeclass
1997# }
1998#
1999#*/###########################################################################
2000setMethodS3(
2001  "setSequences",
2002  class="PhyloSim",
2003  function(
2004    		this,
2005		value,
2006    		...
2007  ){
2008
2009		virtualAssignmentForbidden(this);
2010
2011  },
2012  private=FALSE,
2013  protected=FALSE,
2014  overwrite=FALSE,
2015  conflict="warning",
2016  validators=getOption("R.methodsS3:validators:setMethodS3")
2017);
2018
2019##
2020## Method: getAlignment
2021##
2022###########################################################################/**
2023#
2024# @RdocMethod	getAlignment
2025# 
2026# @title "Get the alignment stored in a PhyloSim object" 
2027# 
2028# \description{ 
2029#	@get "title".
2030# } 
2031# 
2032# @synopsis 
2033# 
2034# \arguments{ 
2035# 	\item{this}{A PhyloSim object.} 
2036# 	\item{...}{Not used.} 
2037# } 
2038# 
2039# \value{ 
2040# 	The alignment as a matrix. Gap are represented by strings composed of dashes.
2041# } 
2042# 
2043# \examples{
2044#	# Create a PhyloSim object.
2045#	# Provide the phylo object 
2046#	# and the root sequence.
2047#	sim<-PhyloSim(
2048#		name="TinySim",
2049#		phylo=rcoal(3),
2050#		root.seq=NucleotideSequence(string="ATG",processes=list(list(JC69())))
2051#	);
2052#	# run the simulation
2053#	Simulate(sim)
2054#	# get the resulting aligment
2055#	getAlignment(sim)
2056#	# via virtual fileld:
2057#	sim$alignment
2058# } 
2059# 
2060# @author 
2061# 
2062# \seealso{ 
2063# 	@seeclass 
2064# } 
2065# 
2066#*/###########################################################################
2067setMethodS3(
2068  "getAlignment",
2069  class="PhyloSim",
2070  function(
2071    this,
2072    ...
2073  ){
2074
2075		this$.alignment;
2076
2077  },
2078  private=FALSE,
2079  protected=FALSE,
2080  overwrite=FALSE,
2081  conflict="warning",
2082  validators=getOption("R.methodsS3:validators:setMethodS3")
2083);
2084
2085##
2086## Method: setAlignment
2087##
2088###########################################################################/**
2089#
2090# @RdocMethod setAlignment
2091#
2092# @title "Forbidden action: setting the alignment stored in a PhyloSim object"
2093#
2094# \description{
2095#       @get "title".
2096# }
2097#
2098# @synopsis
2099#
2100# \arguments{
2101#       \item{this}{An object.}
2102#       \item{value}{Not used.}
2103#       \item{...}{Not used.}
2104# }
2105#
2106# \value{
2107#	Throws an error.
2108# }
2109#
2110# @author
2111#
2112# \seealso{
2113#       @seeclass
2114# }
2115#
2116#*/###########################################################################
2117setMethodS3(
2118  "setAlignment",
2119  class="PhyloSim",
2120  function(
2121    this,
2122    value,
2123    ...
2124  ){
2125
2126               this$alignment <- value;
2127
2128  },
2129  private=FALSE,
2130  protected=FALSE,
2131  overwrite=FALSE,
2132  conflict="warning",
2133  validators=getOption("R.methodsS3:validators:setMethodS3")
2134);
2135
2136##
2137## Method: .recoverAlignment
2138##
2139setMethodS3(
2140  ".recoverAlignment",
2141  class="PhyloSim",
2142  function(
2143    		this,
2144		paranoid=FALSE,
2145    		...
2146  ){
2147
2148		# Refuse to build alignment if at least one of the sequences is NA:
2149		for (seq in this$.sequences){
2150			if(!is.Sequence(seq)){
2151				throw("Cannot build alignment because the simulation is incomplete!\n");
2152			}
2153		}
2154
2155		# The list holding all the partial alignment matrices:
2156		aln.mat<-list();
2157		
2158		# Assigning NA-s here to prevent creation of these variables in the global
2159		# environment.
2160		row.names<-NA;
2161		from.node<-NA;
2162		to.node<-NA;
2163		from.seq<-NA;
2164		to.seq<-NA;
2165		edge<-NA;
2166		from.name<-NA;
2167		to.name<-NA;
2168		from.mat<-NA;
2169		to.mat<-NA;
2170
2171		# Initialize the variables:
2172
2173		init.vars<-function(){
2174
2175			# Getting the edge:
2176			edge<<-getEdge(this, edge.number);
2177			
2178			# Getting the nodes:
2179			from.node<<-edge[[1,"from"]];
2180			to.node<<-edge[[1,"to"]];
2181
2182			# Getting the sequence objects:
2183			from.seq<<-getSeqFromNode(this, from.node)
2184			to.seq<<-getSeqFromNode(this, to.node)
2185	
2186			# Getting sequence names:	
2187			from.name<<-from.seq$name;
2188			to.name<<-to.seq$name;
2189
2190		}
2191
2192		# Initialize the aligment matrices:
2193		init.aln.mats

Large files files are truncated, but you can click here to view the full file