#!/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 MRBAYRS, PAUP ML BOOT and with parameteric bootstrapping.

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, MixtureSimulator
from PIPRes.service_impl.paup_wrap import *
from PIPRes.service_impl.mb_wrap import defaultMBWrap, defaultMBTwoExpWrap
from PIPRes.service_impl.param_boot import CIParamBootFreqSupplier
from PIPRes.util.crude_xml import xmlFilePathToPipres, toXMLGrouping, toPipresXML
from PIPRes.adhoc import *
import itertools
import sys, sets

_LOG = cipresGetLogger('bayesVsParaBoot')
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, trees, splits = readSimulationModelSpecFile(modelSpecFile)
    
    if len(simConditions.mixingProportion) != len(trees):
        sys.exit(modelSpecFile + ' does not hava simulation conditions mixing proportion equal in length to the number of trees')
    simModelProps = []
    for t, p in itertools.izip(trees, simConditions.mixingProportion):
        s = defaultSeqGenWrap()
        s.setTree(t)
        simModelProps.append([s, p])
    simulator = MixtureSimulator(simModelProps)
    simulator.nRealizations = simConditions.nRealizations or 1
    simulator.nSites = simConditions.nSites
    
    splitFreqSuppliers = []

    pBootInferer = defaultPaupWrap()
    pBootInferer.setCriterion(OptimalityCriterion.LIKELIHOOD)
    pBootInferer.execute(' lset nst = 1 basefreq = equal rates = equal ; ')
    pBootSimulator = defaultSeqGenWrap()
    pBootSimulator.nRealizations = 250
    pBootSimulator.nSites = simConditions.nSites
    paramBootstrapper = CIParamBootFreqSupplier(pBootSimulator, pBootInferer, pBootInferer, splits)
    paramBootstrapper.TMP_Name = 'ML parametric boot (PAUP)'
    paramBootstrapper.TMP_Type = 'PBP'
    splitFreqSuppliers.append(paramBootstrapper)
    
    stdMB = defaultMBWrap()
    stdMB.TMP_Name = 'MRBAYES'
    stdMB.TMP_Type = 'PP'
    splitFreqSuppliers.append(stdMB)
    
    npBootPaup = defaultPaupWrap()
    npBootPaup.TMP_Name = 'ML NPBoot (PAUP)'
    npBootPaup.TMP_Type = 'NPBP'
    npBootPaup.setCriterion(OptimalityCriterion.LIKELIHOOD)
    npBootPaup.execute(' set criterion = LIKELIHOOD; lset nst = 1 basefreq = equal rates = equal ; ')
    npBootPaup.nBootReps = 250
    splitFreqSuppliers.append(npBootPaup)
    
    simulator.debugging = debugging
    pBootSimulator.debugging = debugging
    pBootInferer.debugging = debugging
    for s in splitFreqSuppliers:
        s.debugging = debugging
    
    resultsList = simulateAndGetSplitFreqs(simulator, splitFreqSuppliers, splits, tempFile = '.results.xml')
    output.write(toPipresXML([resultsList]))
    print resultsList.__dict__