#!/usr/bin/env Rscript

##
## Example V2.7: Simulating many replicates in parallel
## See also the package vignette (vignette("PhyloSim",package="phylosim")).
##

# load PhyloSim
library("phylosim")

# The the speed of the above method for simulating replicates can be improved on a multicore machine by running
# many replicates in parallel by using the mclapply method from the multicore package. 

# Under default settings, the mclapply method launches one replication per core and this approach needs enough memory to run all
# of them in parallel.

# The following code illustrates how to simulate many replicates in parallel under the JC69+dG model.

# Construct the root sequence object and attach the substitution process:
p<-JC69();
root.seq<-NucleotideSequence(length=50)
attachProcess(root.seq,p)

# Read the required phylogeny from file (this will remain fixed in the simulated replicates):
tree<-read.tree("data/3taxa.nwk");

# Function to simulate a single replication: 
sim.replicate<-function(i){
        name<-paste("replication_",i,sep="")
        clearStates(root.seq)
        plusGamma(root.seq,p,0.25)
        sampleStates(root.seq)

        sim<-Simulate(PhyloSim(
                name=name,
                root.seq=root.seq,
                phylo=tree,
        ),
            quiet=TRUE
        )

        saveAlignment(sim,file=paste("aln_",i,".fas",sep=""))
        return(sim)
        # return(TRUE)
}

# Note that the states are cleared and resampled; the rate multipliers 
# are resampled as well. The resulting alignments are stored in files aln_1.fas, aln_2.fas, aln_3.fas.

# Memory can be saved by throwing away the objects generated by the replication by returning TRUE from the
# sim.replicate function.

# Load the multicore package:
library(multicore)

# Run replicates in parallel:
nr.replicates <-3
res.objects<-mclapply(1:nr.replicates, sim.replicate)

# Print the resulting PhyloSim objects:
print(res.objects)

# Plot one of the resulting alignments:
plot(res.objects[[1]])