#Dan Blankenberg
#Contains file contents and helper methods for HYPHY configurations
import tempfile, os

def get_filled_temp_filename(contents):
    fh = tempfile.NamedTemporaryFile('w')
    filename = fh.name
    fh.close()
    fh = open(filename, 'w')
    fh.write(contents)
    fh.close()
    return filename

NJ_tree_shared_ibf = """
COUNT_GAPS_IN_FREQUENCIES = 0;
methodIndex                  = 1;

/*-----------------------------------------------------------------------------------------------------------------------------------------*/

function InferTreeTopology(verbFlag)
{
    distanceMatrix = {ds.species,ds.species};

    MESSAGE_LOGGING              = 0;
    ExecuteAFile (HYPHY_BASE_DIRECTORY+"TemplateBatchFiles"+DIRECTORY_SEPARATOR+"chooseDistanceFormula.def");
    InitializeDistances (0);
        
    for (i = 0; i<ds.species; i=i+1)
    {
        for (j = i+1; j<ds.species; j = j+1)
        {
            distanceMatrix[i][j] = ComputeDistanceFormula (i,j);
        }
    }

    MESSAGE_LOGGING              = 1;
    cladesMade                     = 1;
    

    if (ds.species == 2)
    {
        d1 = distanceMatrix[0][1]/2;
        treeNodes = {{0,1,d1__},
                     {1,1,d1__},
                     {2,0,0}};
                     
        cladesInfo = {{2,0}};
    }
    else
    {
        if (ds.species == 3)
        {
            /* generate least squares estimates here */
            
            d1 = (distanceMatrix[0][1]+distanceMatrix[0][2]-distanceMatrix[1][2])/2;
            d2 = (distanceMatrix[0][1]-distanceMatrix[0][2]+distanceMatrix[1][2])/2;
            d3 = (distanceMatrix[1][2]+distanceMatrix[0][2]-distanceMatrix[0][1])/2;
            
            treeNodes = {{0,1,d1__},
                         {1,1,d2__},
                         {2,1,d3__}
                         {3,0,0}};
                         
            cladesInfo = {{3,0}};        
        }
        else
        {    
            njm = (distanceMatrix > methodIndex)>=ds.species;
                
            treeNodes         = {2*(ds.species+1),3};
            cladesInfo        = {ds.species-1,2};
            
            for (i=Rows(treeNodes)-1; i>=0; i=i-1)
            {
                treeNodes[i][0] = njm[i][0];
                treeNodes[i][1] = njm[i][1];
                treeNodes[i][2] = njm[i][2];
            }

            for (i=Rows(cladesInfo)-1; i>=0; i=i-1)
            {
                cladesInfo[i][0] = njm[i][3];
                cladesInfo[i][1] = njm[i][4];
            }
            
            njm = 0;
        }
    }
    return 1.0;
}

/*-----------------------------------------------------------------------------------------------------------------------------------------*/

function TreeMatrix2TreeString (doLengths)
{
    treeString = "";
    p = 0;
    k = 0;
    m = treeNodes[0][1];
    n = treeNodes[0][0];
    treeString*(Rows(treeNodes)*25);

    while (m)
    {    
        if (m>p)
        {
            if (p)
            {
                treeString*",";
            }
            for (j=p;j<m;j=j+1)
            {
                treeString*"(";
            }
        }
        else
        {
            if (m<p)
            {
                for (j=m;j<p;j=j+1)
                {
                    treeString*")";
                }
            }    
            else
            {
                treeString*",";
            }    
        }
        if (n<ds.species)
        {
            GetString (nodeName, ds, n);
            if (doLengths != 1)
            {
                treeString*nodeName;            
            }
            else
            {    
                treeString*taxonNameMap[nodeName];
            }
        }
        if (doLengths>.5)
        {
            nodeName = ":"+treeNodes[k][2];
            treeString*nodeName;
        }
        k=k+1;
        p=m;
        n=treeNodes[k][0];
        m=treeNodes[k][1];
    }

    for (j=m;j<p;j=j+1)
    {
        treeString*")";
    }
    
    treeString*0;
    return treeString;
}
"""

def get_NJ_tree (filename):
    return """
DISTANCE_PROMPTS            = 1;
ExecuteAFile ("%s");

DataSet                 ds              = ReadDataFile (PROMPT_FOR_FILE);
DataSetFilter             filteredData = CreateFilter (ds,1);

/* do sequence to branch map */

taxonNameMap = {};

for (k=0; k<ds.species; k=k+1)
{
    GetString         (thisName, ds,k);
    shortName         = (thisName^{{"\\\\..+",""}})&&1;
    taxonNameMap[shortName] = thisName;
    SetParameter (ds,k,shortName);
}

DataSetFilter             filteredData = CreateFilter (ds,1);
InferTreeTopology         (0);
treeString                 = TreeMatrix2TreeString (1);

fprintf                    (PROMPT_FOR_FILE, CLEAR_FILE, treeString);
fscanf                    (stdin, "String", ps_file);

if (Abs(ps_file))
{
    treeString = TreeMatrix2TreeString (2);
    UseModel (USE_NO_MODEL);
    Tree givenTree = treeString;
    baseHeight         = TipCount (givenTree)*28;
    TREE_OUTPUT_OPTIONS = {};
    TREE_OUTPUT_OPTIONS["__FONT_SIZE__"] = 14;
    baseWidth = 0;
    treeAVL                = givenTree^0;
    drawLetter            = "/drawletter {"+TREE_OUTPUT_OPTIONS["__FONT_SIZE__"]$4+" -"+TREE_OUTPUT_OPTIONS["__FONT_SIZE__"]$2+ " show} def\\n";
    for (k3 = 1; k3 < Abs(treeAVL); k3=k3+1)
    {
        nodeName = (treeAVL[k3])["Name"];
        if(Abs((treeAVL[k3])["Children"]) == 0)
        {
            mySpecs = {};
            mySpecs ["TREE_OUTPUT_BRANCH_LABEL"] = "(" + taxonNameMap[nodeName] + ") drawLetter";
            baseWidth = Max (baseWidth, (treeAVL[k3])["Depth"]);
        }
    }
    baseWidth = 40*baseWidth;
    
    fprintf (ps_file,  CLEAR_FILE, drawLetter, PSTreeString (givenTree, "STRING_SUPPLIED_LENGTHS",{{baseWidth,baseHeight}}));
}
""" % (filename)

def get_NJ_treeMF (filename):
    return  """
ExecuteAFile ("%s");

VERBOSITY_LEVEL            = -1;
fscanf                  (PROMPT_FOR_FILE, "Lines", inLines);

_linesIn                        = Columns (inLines);
isomorphicTreesBySequenceCount  = {};

/*---------------------------------------------------------*/

_currentGene   = 1; 
_currentState = 0;
geneSeqs      = "";
geneSeqs      * 128;

fprintf (PROMPT_FOR_FILE, CLEAR_FILE, KEEP_OPEN);
treeOutFile = LAST_FILE_PATH;

fscanf  (stdin,"String", ps_file);
if (Abs(ps_file))
{
    fprintf (ps_file, CLEAR_FILE, KEEP_OPEN);
}

for (l=0; l<_linesIn; l=l+1)
{
    if (Abs(inLines[l]) == 0)
    {
        if (_currentState == 1)
        {
            geneSeqs      * 0;
            DataSet       ds            = ReadFromString (geneSeqs);
            _processAGene (_currentGene,treeOutFile,ps_file);
            geneSeqs * 128;
            _currentGene = _currentGene + 1;
        }
    }
    else
    {
        if (_currentState == 0)
        {
            _currentState = 1;
        }
        geneSeqs * inLines[l];
        geneSeqs * "\\n";
    }
}


if (_currentState == 1)
{
    geneSeqs      * 0;
    if (Abs(geneSeqs))
    {
        DataSet       ds            = ReadFromString (geneSeqs);
        _processAGene (_currentGene,treeOutFile,ps_file);
    }
}

fprintf (treeOutFile,CLOSE_FILE);
if (Abs(ps_file))
{
    fprintf (ps_file,CLOSE_FILE);
}
/*---------------------------------------------------------*/

function _processAGene (_geneID, nwk_file, ps_file)
{
    if (ds.species == 1)
    {
        fprintf                    (nwk_file, _geneID-1, "\\tNone \\tNone\\n");
        return 0;
        
    }
    
    DataSetFilter             filteredData = CreateFilter (ds,1);

    /* do sequence to branch map */

    taxonNameMap = {};

    for (k=0; k<ds.species; k=k+1)
    {
        GetString         (thisName, ds,k);
        shortName         = (thisName^{{"\\\\..+",""}});
        taxonNameMap[shortName] = thisName;
        SetParameter (ds,k,shortName);
    }

    DataSetFilter             filteredData = CreateFilter (ds,1);
    DISTANCE_PROMPTS        = (_geneID==1);

    InferTreeTopology         (0);
    baseTree                = TreeMatrix2TreeString (0);
    UseModel                 (USE_NO_MODEL);
    
    Tree             baseTop    = baseTree;
    
    /* standardize this top */
    
    for (k=0; k<Abs(isomorphicTreesBySequenceCount[filteredData.species]); k=k+1)
    {
        testString      = (isomorphicTreesBySequenceCount[filteredData.species])[k];
        Tree testTree = testString;
        if (testTree == baseTop)
        {
            baseTree = testString;
            break;
        }
    }
    if (k==Abs(isomorphicTreesBySequenceCount[filteredData.species]))
    {
        if (k==0)
        {
            isomorphicTreesBySequenceCount[filteredData.species] = {};
        }
        (isomorphicTreesBySequenceCount[filteredData.species])[k] = baseTree;
    }
    
    fprintf                    (nwk_file, _geneID-1, "\\t", baseTree, "\\t", TreeMatrix2TreeString (1), "\\n");
    if (Abs(ps_file))
    {
        treeString = TreeMatrix2TreeString (2);
        UseModel (USE_NO_MODEL);
        Tree givenTree = treeString;
        baseHeight         = TipCount (givenTree)*28;
        TREE_OUTPUT_OPTIONS = {};
        TREE_OUTPUT_OPTIONS["__FONT_SIZE__"] = 14;
        baseWidth = 0;
        treeAVL                = givenTree^0;
        drawLetter            = "/drawletter {"+TREE_OUTPUT_OPTIONS["__FONT_SIZE__"]$4+" -"+TREE_OUTPUT_OPTIONS["__FONT_SIZE__"]$2+ " show} def\\n";
        for (k3 = 1; k3 < Abs(treeAVL); k3=k3+1)
        {
            nodeName = (treeAVL[k3])["Name"];
            if(Abs((treeAVL[k3])["Children"]) == 0)
            {
                mySpecs = {};
                mySpecs ["TREE_OUTPUT_BRANCH_LABEL"] = "(" + taxonNameMap[nodeName] + ") drawLetter";
                baseWidth = Max (baseWidth, (treeAVL[k3])["Depth"]);
            }
        }
        baseWidth = 40*baseWidth;
        
        fprintf (stdout, _geneID, ":", givenTree,"\\n");
        fprintf (ps_file, PSTreeString (givenTree, "STRING_SUPPLIED_LENGTHS",{{baseWidth,baseHeight}}));
    }
    return 0;
}
""" % (filename)

BranchLengthsMF = """
VERBOSITY_LEVEL            = -1;

fscanf                  (PROMPT_FOR_FILE, "Lines", inLines);



_linesIn = Columns (inLines);



/*---------------------------------------------------------*/



_currentGene   = 1; 

_currentState = 0;

geneSeqs      = "";

geneSeqs      * 128;



for (l=0; l<_linesIn; l=l+1)

{

    if (Abs(inLines[l]) == 0)

    {

        if (_currentState == 1)

        {

            geneSeqs      * 0;

            DataSet       ds            = ReadFromString (geneSeqs);

            _processAGene (_currentGene);

             geneSeqs * 128;

            _currentGene = _currentGene + 1;

        }

    }

    else

    {

        if (_currentState == 0)

        {

            _currentState = 1;

        }

        geneSeqs * inLines[l];

        geneSeqs * "\\n";

    }

}



if (_currentState == 1)

{

    geneSeqs      * 0;

    if (Abs(geneSeqs))

    {

        DataSet       ds            = ReadFromString (geneSeqs);

        _processAGene (_currentGene);

    }

}



fprintf (resultFile,CLOSE_FILE);



/*---------------------------------------------------------*/



function _processAGene (_geneID)

{

    DataSetFilter             filteredData = CreateFilter (ds,1);

    if (_currentGene == 1)

    {

        SelectTemplateModel        (filteredData);

        

        SetDialogPrompt           ("Tree file");

        fscanf                    (PROMPT_FOR_FILE, "Tree",  givenTree);

        fscanf                    (stdin, "String", resultFile);

        

        /* do sequence to branch map */

        

        validNames = {};

        taxonNameMap = {};

        

        for (k=0; k<TipCount(givenTree); k=k+1)

        {

            validNames[TipName(givenTree,k)&&1] = 1;

        }

        

        for (k=0; k<BranchCount(givenTree); k=k+1)

        {

            thisName = BranchName(givenTree,k);

            taxonNameMap[thisName&&1] = thisName;

        }

        

        storeValidNames = validNames;

        fprintf         (resultFile,CLEAR_FILE,KEEP_OPEN,"Block\\tBranch\\tLength\\tLowerBound\\tUpperBound\\n"); 

    }

    else

    {

        HarvestFrequencies (vectorOfFrequencies, filteredData, 1,1,1);

        validNames = storeValidNames;

    }

    

    for (k=0; k<ds.species; k=k+1)

    {

        GetString (thisName, ds,k);

        shortName = (thisName^{{"\\\\..+",""}})&&1;

        if (validNames[shortName])

        {

            taxonNameMap[shortName] = thisName;

            validNames - (shortName);

            SetParameter (ds,k,shortName);

        }

        else

        {

            fprintf         (resultFile,"ERROR:", thisName, " could not be matched to any of the leaves in tree ", givenTree,"\\n"); 

            return 0;

        }

    }

    

    /* */

    

    LikelihoodFunction lf = (filteredData,givenTree);

    Optimize                 (res,lf);

    

    timer = Time(0)-timer;

    

    branchNames   = BranchName   (givenTree,-1);

    branchLengths = BranchLength (givenTree,-1);

    

    

    for (k=0; k<Columns(branchNames)-1; k=k+1)

    {

        COVARIANCE_PARAMETER = "givenTree."+branchNames[k]+".t";

        COVARIANCE_PRECISION = 0.95;

        CovarianceMatrix (cmx,lf);

        if (k==0)

        {

            /* compute a scaling factor */

            ExecuteCommands ("givenTree."+branchNames[0]+".t=1");

            scaleFactor = BranchLength (givenTree,0);

            ExecuteCommands ("givenTree."+branchNames[0]+".t="+cmx[0][1]);

        }

        fprintf (resultFile,_geneID,"\\t",taxonNameMap[branchNames[k]&&1],"\\t",branchLengths[k],"\\t",scaleFactor*cmx[0][0],"\\t",scaleFactor*cmx[0][2],"\\n");

    }

    

    ttl = (branchLengths*(Transpose(branchLengths["1"])))[0];

    global treeScaler = 1;

    ReplicateConstraint ("this1.?.t:=treeScaler*this2.?.t__",givenTree,givenTree);

    COVARIANCE_PARAMETER = "treeScaler";

    COVARIANCE_PRECISION = 0.95;

    CovarianceMatrix (cmx,lf);

    fprintf (resultFile,_geneID,"\\tTotal Tree\\t",ttl,"\\t",ttl*cmx[0][0],"\\t",ttl*cmx[0][2],"\\n");

    ClearConstraints (givenTree);

    return 0;

}
"""

BranchLengths = """
DataSet                 ds              = ReadDataFile (PROMPT_FOR_FILE);
DataSetFilter             filteredData = CreateFilter (ds,1);

SelectTemplateModel        (filteredData);

SetDialogPrompt         ("Tree file");
fscanf                    (PROMPT_FOR_FILE, "Tree",  givenTree);
fscanf                    (stdin, "String", resultFile);

/* do sequence to branch map */

validNames = {};
taxonNameMap = {};

for (k=0; k<TipCount(givenTree); k=k+1)
{
    validNames[TipName(givenTree,k)&&1] = 1;
}

for (k=0; k<BranchCount(givenTree); k=k+1)
{
    thisName = BranchName(givenTree,k);
    taxonNameMap[thisName&&1] = thisName;
}

for (k=0; k<ds.species; k=k+1)
{
    GetString (thisName, ds,k);
    shortName = (thisName^{{"\\\\..+",""}})&&1;
    if (validNames[shortName])
    {
        taxonNameMap[shortName] = thisName;
        validNames - (shortName);
        SetParameter (ds,k,shortName);
    }
    else
    {
        fprintf         (resultFile,CLEAR_FILE,"ERROR:", thisName, " could not be matched to any of the leaves in tree ", givenTree); 
        return 0;
    }
}

/* */

LikelihoodFunction lf = (filteredData,givenTree);

Optimize                 (res,lf);

timer = Time(0)-timer;

branchNames   = BranchName   (givenTree,-1);
branchLengths = BranchLength (givenTree,-1);

fprintf         (resultFile,CLEAR_FILE,KEEP_OPEN,"Branch\\tLength\\tLowerBound\\tUpperBound\\n"); 

for (k=0; k<Columns(branchNames)-1; k=k+1)
{
    COVARIANCE_PARAMETER = "givenTree."+branchNames[k]+".t";
    COVARIANCE_PRECISION = 0.95;
    CovarianceMatrix (cmx,lf);
    if (k==0)
    {
        /* compute a scaling factor */
        ExecuteCommands ("givenTree."+branchNames[0]+".t=1");
        scaleFactor = BranchLength (givenTree,0);
        ExecuteCommands ("givenTree."+branchNames[0]+".t="+cmx[0][1]);
    }
    fprintf (resultFile,taxonNameMap[branchNames[k]&&1],"\\t",branchLengths[k],"\\t",scaleFactor*cmx[0][0],"\\t",scaleFactor*cmx[0][2],"\\n");
}

ttl = (branchLengths*(Transpose(branchLengths["1"])))[0];
global treeScaler = 1;
ReplicateConstraint ("this1.?.t:=treeScaler*this2.?.t__",givenTree,givenTree);
COVARIANCE_PARAMETER = "treeScaler";
COVARIANCE_PRECISION = 0.95;
CovarianceMatrix (cmx,lf);
ClearConstraints (givenTree);
fprintf (resultFile,"Total Tree\\t",ttl,"\\t",ttl*cmx[0][0],"\\t",ttl*cmx[0][2],"\\n");
fprintf (resultFile,CLOSE_FILE);
"""

SimpleLocalFitter = """
VERBOSITY_LEVEL               = -1;
COUNT_GAPS_IN_FREQUENCIES  = 0;

/*---------------------------------------------------------*/

function returnResultHeaders (dummy)
{
    _analysisHeaders = {};
    _analysisHeaders[0] = "BLOCK";
    _analysisHeaders[1] = "BP";
    _analysisHeaders[2] = "S_sites";
    _analysisHeaders[3] = "NS_sites";
    _analysisHeaders[4] = "Stop_codons";
    _analysisHeaders[5] = "LogL";
    _analysisHeaders[6] = "AC";
    _analysisHeaders[7] = "AT";
    _analysisHeaders[8] = "CG";
    _analysisHeaders[9] = "CT";
    _analysisHeaders[10] = "GT";
    _analysisHeaders[11] = "Tree";

    for (_biterator = 0; _biterator < treeBranchCount; _biterator = _biterator + 1)
    {
        branchName = treeBranchNames[_biterator];
        
        _analysisHeaders [Abs(_analysisHeaders)] = "length("+branchName+")";
        _analysisHeaders [Abs(_analysisHeaders)] = "dS("+branchName+")";
        _analysisHeaders [Abs(_analysisHeaders)] = "dN("+branchName+")";
        _analysisHeaders [Abs(_analysisHeaders)] = "omega("+branchName+")";
    }

    return _analysisHeaders;
}

/*---------------------------------------------------------*/

function runAGeneFit (myID)
{
    DataSetFilter filteredData = CreateFilter (ds,3,"","",GeneticCodeExclusions);

    if (_currentGene==1)
    {
        _MG94stdinOverload = {};
        _MG94stdinOverload ["0"] = "Local";
        _MG94stdinOverload ["1"] = modelSpecString;

        ExecuteAFile             (HYPHY_BASE_DIRECTORY+"TemplateBatchFiles"+DIRECTORY_SEPARATOR+"TemplateModels"+DIRECTORY_SEPARATOR+"MG94custom.mdl",
                                 _MG94stdinOverload);
        
        Tree    codonTree           = treeString;
    }
    else
    {
        HarvestFrequencies               (observedFreq,filteredData,3,1,1);
        MULTIPLY_BY_FREQS             = PopulateModelMatrix ("MG94custom", observedFreq);
        vectorOfFrequencies         = BuildCodonFrequencies (observedFreq);
        Model MG94customModel         = (MG94custom,vectorOfFrequencies,0);

        Tree    codonTree            = treeString;
    }

    LikelihoodFunction lf     = (filteredData,codonTree);

    Optimize                     (res,lf);

    _snsAVL       =                 _computeSNSSites ("filteredData", _Genetic_Code, vectorOfFrequencies, 0);
    _cL           =                  ReturnVectorsOfCodonLengths (ComputeScalingStencils (0), "codonTree");


    _returnMe = {};
    _returnMe ["BLOCK"]              = myID;
    _returnMe ["LogL"]              = res[1][0];
    _returnMe ["BP"]                 = _snsAVL ["Sites"];
    _returnMe ["S_sites"]             = _snsAVL ["SSites"];
    _returnMe ["NS_sites"]             = _snsAVL ["NSSites"];
    _returnMe ["AC"]                 = AC;
    _returnMe ["AT"]                 = AT;
    _returnMe ["CG"]                 = CG;
    _returnMe ["CT"]                 = CT;
    _returnMe ["GT"]                 = GT;
    _returnMe ["Tree"]                 = Format(codonTree,0,1);

    for (_biterator = 0; _biterator < treeBranchCount; _biterator = _biterator + 1)
    {
        branchName = treeBranchNames[_biterator];

        _returnMe ["length("+branchName+")"]         = (_cL["Total"])[_biterator];
        _returnMe ["dS("+branchName+")"]                 = (_cL["Syn"])[_biterator]*(_returnMe ["BP"]/_returnMe ["S_sites"]);
        _returnMe ["dN("+branchName+")"]                 = (_cL["NonSyn"])[_biterator]*(_returnMe ["BP"]/_returnMe ["NS_sites"]);
        
        ExecuteCommands ("_lom = _standardizeRatio(codonTree."+treeBranchNames[_biterator]+".nonSynRate,codonTree."+treeBranchNames[_biterator]+".synRate);");
        _returnMe ["omega("+branchName+")"]                 = _lom;
    }
    
    return _returnMe;
}

"""

SimpleGlobalFitter = """
VERBOSITY_LEVEL               = -1;
COUNT_GAPS_IN_FREQUENCIES  = 0;

/*---------------------------------------------------------*/

function returnResultHeaders (dummy)
{
    _analysisHeaders = {};
    _analysisHeaders[0]  = "BLOCK";
    _analysisHeaders[1]  = "BP";
    _analysisHeaders[2]  = "S_sites";
    _analysisHeaders[3]  = "NS_sites";
    _analysisHeaders[4]  = "Stop_codons";
    _analysisHeaders[5]  = "LogL";
    _analysisHeaders[6]  = "omega";
    _analysisHeaders[7]  = "omega_range";
    _analysisHeaders[8]  = "AC";
    _analysisHeaders[9]  = "AT";
    _analysisHeaders[10] = "CG";
    _analysisHeaders[11] = "CT";
    _analysisHeaders[12] = "GT";
    _analysisHeaders[13] = "Tree";

    return _analysisHeaders;
}

/*---------------------------------------------------------*/

function runAGeneFit (myID)
{
    fprintf (stdout, "[SimpleGlobalFitter.bf on GENE ", myID, "]\\n");
    taxonNameMap = {};
    
    for (k=0; k<ds.species; k=k+1)
    {
        GetString         (thisName, ds,k);
        shortName         = (thisName^{{"\\\\..+",""}})&&1;
        taxonNameMap[shortName] = thisName;
        SetParameter (ds,k,shortName);
    }

    DataSetFilter filteredData = CreateFilter (ds,1);
    _nucSites                    = filteredData.sites;
    
    if (Abs(treeString))
    {
        givenTreeString = treeString;
    }
    else
    {
        if (_currentGene==1)
        {
            ExecuteAFile             (HYPHY_BASE_DIRECTORY+"TemplateBatchFiles"+DIRECTORY_SEPARATOR+"Utility"+DIRECTORY_SEPARATOR+"NJ.bf");
        }
        givenTreeString = InferTreeTopology (0);
        treeString        = "";
    }

    DataSetFilter filteredData = CreateFilter (ds,3,"","",GeneticCodeExclusions);

    if (_currentGene==1)
    {
        _MG94stdinOverload = {};
        _MG94stdinOverload ["0"] = "Global";
        _MG94stdinOverload ["1"] = modelSpecString;

        ExecuteAFile             (HYPHY_BASE_DIRECTORY+"TemplateBatchFiles"+DIRECTORY_SEPARATOR+"TemplateModels"+DIRECTORY_SEPARATOR+"MG94custom.mdl",
                                 _MG94stdinOverload);
        
        Tree    codonTree           = givenTreeString;
    }
    else
    {
        HarvestFrequencies               (observedFreq,filteredData,3,1,1);
        MULTIPLY_BY_FREQS             = PopulateModelMatrix ("MG94custom", observedFreq);
        vectorOfFrequencies         = BuildCodonFrequencies (observedFreq);
        Model MG94customModel         = (MG94custom,vectorOfFrequencies,0);

        Tree    codonTree            = givenTreeString;
    }

    LikelihoodFunction lf     = (filteredData,codonTree);

    Optimize                     (res,lf);

    _snsAVL       =                 _computeSNSSites ("filteredData", _Genetic_Code, vectorOfFrequencies, 0);
    _cL           =                  ReturnVectorsOfCodonLengths (ComputeScalingStencils (0), "codonTree");


    _returnMe = {};
    _returnMe ["BLOCK"]              = myID;
    _returnMe ["LogL"]              = res[1][0];
    _returnMe ["BP"]                 = _snsAVL ["Sites"];
    _returnMe ["S_sites"]             = _snsAVL ["SSites"];
    _returnMe ["NS_sites"]             = _snsAVL ["NSSites"];
    _returnMe ["Stop_codons"]         = (_nucSites-filteredData.sites*3)$3;
    _returnMe ["AC"]                 = AC;
    _returnMe ["AT"]                 = AT;
    _returnMe ["CG"]                 = CG;
    _returnMe ["CT"]                 = CT;
    _returnMe ["GT"]                 = GT;
    _returnMe ["omega"]             = R;
    COVARIANCE_PARAMETER             = "R";
    COVARIANCE_PRECISION             = 0.95;
    CovarianceMatrix                 (cmx,lf);
    _returnMe ["omega_range"]         = ""+cmx[0]+"-"+cmx[2];
    _returnMe ["Tree"]                 = Format(codonTree,0,1);


    return _returnMe;
}
"""

FastaReader = """
fscanf            (stdin, "String", _coreAnalysis);
fscanf            (stdin, "String", _outputDriver);

ExecuteAFile             (HYPHY_BASE_DIRECTORY+"TemplateBatchFiles"+DIRECTORY_SEPARATOR+"TemplateModels"+DIRECTORY_SEPARATOR+"chooseGeneticCode.def");
ExecuteAFile             (HYPHY_BASE_DIRECTORY+"TemplateBatchFiles"+DIRECTORY_SEPARATOR+"dSdNTreeTools.ibf");
ExecuteAFile             (HYPHY_BASE_DIRECTORY+"TemplateBatchFiles"+DIRECTORY_SEPARATOR+"Utility"+DIRECTORY_SEPARATOR+"CodonTools.bf");
ExecuteAFile             (HYPHY_BASE_DIRECTORY+"TemplateBatchFiles"+DIRECTORY_SEPARATOR+"Utility"+DIRECTORY_SEPARATOR+"GrabBag.bf");

SetDialogPrompt ("Tree file");
fscanf            (PROMPT_FOR_FILE, "Tree",  givenTree);

treeBranchNames               = BranchName (givenTree,-1);
treeBranchCount               = Columns    (treeBranchNames)-1;
treeString                    = Format (givenTree,1,1);

SetDialogPrompt ("Multiple gene FASTA file");
fscanf          (PROMPT_FOR_FILE, "Lines", inLines);
fscanf            (stdin, "String",  modelSpecString);
fscanf            (stdin, "String", _outPath);
 
ExecuteAFile    (_outputDriver);
ExecuteAFile    (_coreAnalysis);
 
/*---------------------------------------------------------*/

_linesIn     = Columns (inLines);
_currentGene   = 1; 
 _currentState = 0;
/* 0 - waiting for a non-empty line */
/* 1 - reading files */

geneSeqs       = "";
geneSeqs        * 0;

_prepareFileOutput (_outPath);

for (l=0; l<_linesIn; l=l+1)
{
    if (Abs(inLines[l]) == 0)
    {
        if (_currentState == 1)
        {
            geneSeqs      * 0;
            DataSet       ds            = ReadFromString (geneSeqs);
            _processAGene (ds.species == treeBranchCount,_currentGene);
            geneSeqs * 128;
            _currentGene = _currentGene + 1;
        }
    }
    else
    {
        if (_currentState == 0)
        {
            _currentState = 1;
        }
        geneSeqs * inLines[l];
        geneSeqs * "\\n";
    }
}

if (_currentState == 1)
{
    geneSeqs      * 0;
    DataSet       ds            = ReadFromString (geneSeqs);
    _processAGene (ds.species == treeBranchCount,_currentGene);
}

_finishFileOutput (0);
"""

TabWriter = """
/*---------------------------------------------------------*/
function _prepareFileOutput (_outPath)
{
    _outputFilePath = _outPath;
    
    _returnHeaders     = returnResultHeaders(0);
    
    fprintf (_outputFilePath, CLEAR_FILE, KEEP_OPEN, _returnHeaders[0]);
    for (_biterator = 1; _biterator < Abs(_returnHeaders); _biterator = _biterator + 1)
    {
        fprintf (_outputFilePath,"\\t",_returnHeaders[_biterator]);
    }


    
    fprintf (_outputFilePath,"\\n");
    return 0;
}

/*---------------------------------------------------------*/

function _processAGene (valid, _geneID)
{
    if (valid)
    {
        returnValue = runAGeneFit (_geneID);
        fprintf (_outputFilePath, returnValue[_returnHeaders[0]]);
        for (_biterator = 1; _biterator < Abs(_returnHeaders); _biterator = _biterator + 1)
        {
            fprintf (_outputFilePath,"\\t",returnValue[_returnHeaders[_biterator]]);
        }
        fprintf (_outputFilePath, "\\n");
    }
    /*
    else
    {
        fprintf (_outputFilePath, 
                _geneID, ", Incorrect number of sequences\\n");
    }
    */
    _currentState = 0;
    return 0;
}

/*---------------------------------------------------------*/
function _finishFileOutput (dummy)
{
    return 0;
}
"""

def get_dnds_config_filename(Fitter_filename, TabWriter_filename, genetic_code, tree_filename, input_filename, nuc_model, output_filename, FastaReader_filename ):
    contents = """
_genomeScreenOptions = {};

/* all paths are either absolute or relative 
to the DATA READER */

_genomeScreenOptions ["0"] = "%s"; 
    /* which analysis to run on each gene; */    
_genomeScreenOptions ["1"] = "%s"; 
    /* what output to produce; */    
_genomeScreenOptions ["2"] = "%s"; 
    /* genetic code */    
_genomeScreenOptions ["3"] = "%s"; 
    /* tree file */
_genomeScreenOptions ["4"] = "%s"; 
    /* alignment file */
_genomeScreenOptions ["5"] = "%s"; 
    /* nucleotide bias string; can define any of the 203 models */
_genomeScreenOptions ["6"] = "%s"; 
    /* output csv file */
    
ExecuteAFile ("%s", _genomeScreenOptions);    
""" % (Fitter_filename, TabWriter_filename, genetic_code, tree_filename, input_filename, nuc_model, output_filename, FastaReader_filename )
    return get_filled_temp_filename(contents)
    
    
def get_branch_lengths_config_filename(input_filename, nuc_model, model_options, base_freq, tree_filename, output_filename, BranchLengths_filename):
    contents = """
_genomeScreenOptions = {};

/* all paths are either absolute or relative 
to the NucDataBranchLengths.bf */

_genomeScreenOptions ["0"] = "%s"; 
    /* the file to analyze; */    
_genomeScreenOptions ["1"] = "CUSTOM"; 
    /* use an arbitrary nucleotide model */    
_genomeScreenOptions ["2"] = "%s";
    /* which model to use */ 
_genomeScreenOptions ["3"] = "%s"; 
    /* model options */
_genomeScreenOptions ["4"] = "Estimated"; 
    /* rate parameters */
_genomeScreenOptions ["5"] = "%s"; 
    /* base frequencies */
_genomeScreenOptions ["6"] = "%s"; 
    /* the tree to use; */    
_genomeScreenOptions ["7"] = "%s"; 
    /* write .csv output to; */    
    
ExecuteAFile ("%s", _genomeScreenOptions);    
""" % (input_filename, nuc_model, model_options, base_freq, tree_filename, output_filename, BranchLengths_filename)
    return get_filled_temp_filename(contents)
    
    
def get_nj_tree_config_filename(input_filename, distance_metric, output_filename1, output_filename2, NJ_tree_filename):
    contents = """
_genomeScreenOptions = {};

/* all paths are either absolute or relative 
to the BuildNJTree.bf */

_genomeScreenOptions ["0"] = "%s"; 
    /* the file to analyze; */    
_genomeScreenOptions ["1"] = "%s"; 
    /* pick which distance metric to use; TN93 is a good default */    
_genomeScreenOptions ["2"] = "%s"; 
    /* write Newick tree output to; */    
_genomeScreenOptions ["3"] = "%s"; 
    /* write a postscript tree file to this file; leave blank to not write a tree */    
    
ExecuteAFile ("%s", _genomeScreenOptions);    
""" % (input_filename, distance_metric, output_filename1, output_filename2, NJ_tree_filename)
    return get_filled_temp_filename(contents)
    
    
def get_nj_treeMF_config_filename(input_filename, output_filename1, output_filename2, distance_metric, NJ_tree_filename):
    contents = """
_genomeScreenOptions = {};

/* all paths are either absolute or relative 
to the BuildNJTreeMF.bf */

_genomeScreenOptions ["0"] = "%s"; 
    /* the multiple alignment file to analyze; */    
_genomeScreenOptions ["1"] = "%s"; 
    /* write Newick tree output to; */    
_genomeScreenOptions ["2"] = "%s"; 
    /* write a postscript tree file to this file; leave blank to not write a tree */    
_genomeScreenOptions ["3"] = "%s"; 
    /* pick which distance metric to use; TN93 is a good default */    
    
ExecuteAFile ("%s", _genomeScreenOptions);        
""" % (input_filename, output_filename1, output_filename2, distance_metric, NJ_tree_filename)
    return get_filled_temp_filename(contents)