#!/usr/bin/python
'''This script is an example of how to use PIPRes to perform a simulation study.
The script takes commmand line arguments:
    <infile> <outfile>
Data is simulated according to the model tree specified in infile.
Each realization is analyzed using a "default" MRBAYES run and a run using Yang and Rannala's TwoExpMRBAYES,
for comparison PAUP is used to generate parsimony and ML bootstrap proportions.

Split frequencies from the methods are written in outfile.
'''
from PIPRes.cipres_types import *
from PIPRes.basic import *
from PIPRes.tree import *
from PIPRes.service_impl.seq_gen_wrap import defaultSeqGenWrap
from PIPRes.service_impl.paup_wrap import defaultPaupWrap
from PIPRes.service_impl.mb_wrap import defaultMBWrap, defaultMBTwoExpWrap
from PIPRes.util.crude_xml import xmlFilePathToPipres, toXMLGrouping, toPipresXML
from PIPRes.adhoc import simulateAndGetSplitFreqs
import sys, sets

_LOG = cipresGetLogger('edgeLenSim')
class Statistic:pass

if __name__ == '__main__':
    if len(sys.argv) < 2:
        sys.exit('Expecting a model conditions file as an argument')
    modelSpecFile = sys.argv[1]
    if len(sys.argv) > 3:
        sys.exit('Expecting at most <infile> <outfile> as arguments')
    output = len(sys.argv) == 3 and open(sys.argv[2], 'a') or sys.stdout
    savingTemp = True
    debugging = True
    
    simConditions, tree, splits = readSimulationModelSpecFile(modelSpecFile)
    _LOG.debug(simConditions.__dict__)
    
    simulator = defaultSeqGenWrap()
    simulator.nRealizations = simConditions.nRealizations or 1
    simulator.nSites = simConditions.nSites
    #simulator.setSeed(1944846395) #############TEMPORARRY
    #_LOG.warn('HARD-CODED seed')
    simulator.setTree(tree)
    
    splitFreqSuppliers = []
    stdMB = defaultMBWrap()
    stdMB.TMP_Name = 'MRBAYES'
    stdMB.TMP_Type = 'PP'
    splitFreqSuppliers.append(stdMB)
    
    twoExpMB = defaultMBTwoExpWrap()
    twoExpMB.TMP_Name = 'TwoExpMRBAYES1e-6'
    twoExpMB.TMP_Type = 'PP'
    twoExpMB.internalMean = 1.0e-6
    twoExpMB.externalMean = 1.0
    twoExpMB.nGenerations = 1000000
    
    #_LOG.warn('very short mb run')
    #twoExpMB.nGenerations = 1000
    splitFreqSuppliers.append(twoExpMB)
   
    twoExpMBModerate = defaultMBTwoExpWrap()
    twoExpMBModerate.TMP_Name = 'TwoExpMRBAYES1e-3'
    twoExpMBModerate.TMP_Type = 'PP'
    twoExpMBModerate.internalMean = 1.0e-3
    twoExpMBModerate.externalMean = 1.0
    twoExpMBModerate.nGenerations = 1000000
    #_LOG.warn('very short mb run')
    #twoExpMBModerate.nGenerations = 1000
    splitFreqSuppliers.append(twoExpMBModerate)
   
    parsPaup = defaultPaupWrap()
    parsPaup.TMP_Name = 'Pars (PAUP)'
    parsPaup.TMP_Type = 'BP'
    splitFreqSuppliers.append(parsPaup)
    
    mlPaup = defaultPaupWrap()
    mlPaup.TMP_Name = 'ML (PAUP)'
    mlPaup.TMP_Type = 'BP'
    mlPaup.execute(' set criterion = likelihood; lset nst = 1 basefreq = equal rates = equal ; ')
    mlPaup.nBootReps = 250
    splitFreqSuppliers.append(mlPaup)
    
    nTaxa = 8
    if debugging:
        for s in splitFreqSuppliers:
            s.debugging = True
    
    resultsList = simulateAndGetSplitFreqs(simulator, splitFreqSuppliers, splits, tempFile = '.results.xml')
    output.write(toPipresXML([resultsList]))