From fc8067dff8b55dcc13e62238b5bde9cecc64bba5 Mon Sep 17 00:00:00 2001 From: Segei L Kosakovsky Pond Date: Sat, 23 Jan 2021 10:59:58 -0500 Subject: [PATCH 1/4] Caching conditon fixes --- src/core/batchlan.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/core/batchlan.cpp b/src/core/batchlan.cpp index 34bd2be69..008346bf9 100644 --- a/src/core/batchlan.cpp +++ b/src/core/batchlan.cpp @@ -2592,10 +2592,9 @@ void _ElementaryCommand::ExecuteCase4 (_ExecutionList& chain) { HBLObjectRef result; if (expression) { //printf ("\n*** Interpreted condition\n"); - result = expression->Compute(); + result = expression->Compute(0, nil, nil, nil, HY_ANY_OBJECT, false); } else { - //printf ("\n*** Compiled condition\n"); - result = ((_Formula*)simpleParameters(2))->Compute(); + result = ((_Formula*)simpleParameters(2))->Compute(0, nil, nil, nil, HY_ANY_OBJECT, false); } // printf ("\n*** %s\n", ((_String*)result->toStr())->sData); From b4adfca8074830a73c3b93094c60069facd0c40a Mon Sep 17 00:00:00 2001 From: Sergei Date: Mon, 25 Jan 2021 10:01:31 -0500 Subject: [PATCH 2/4] Estimator utility functions --- res/TemplateBatchFiles/F_ST.bf | 3 +- .../SelectionAnalyses/FEL.bf | 95 +++++++++++++------ .../SelectionAnalyses/MEME.bf | 71 ++++++++------ .../libv3/tasks/estimators.bf | 30 ++++++ src/core/global_things.cpp | 3 +- 5 files changed, 141 insertions(+), 61 deletions(-) diff --git a/res/TemplateBatchFiles/F_ST.bf b/res/TemplateBatchFiles/F_ST.bf index 7e522f333..b8110e898 100644 --- a/res/TemplateBatchFiles/F_ST.bf +++ b/res/TemplateBatchFiles/F_ST.bf @@ -239,8 +239,7 @@ if (distanceChoice == 1) fprintf (stdout,"\nHYPHY is computing pairwise maximum likelihood distance estimates. A total of ", Format(ds.species*(ds.species-1)/2,0,0), " estimations will be performed.\n"); - MESSAGE_LOGGING = 0; - + for (i = 0; i Date: Mon, 25 Jan 2021 14:55:11 -0500 Subject: [PATCH 3/4] 2.5.27rc --- res/TemplateBatchFiles/Distances/JC69 | 12 +- res/TemplateBatchFiles/Distances/K2P | 2 +- res/TemplateBatchFiles/Distances/K2P_RV | 2 +- .../Distances/Modified_Nei_Gojobori | 2 +- res/TemplateBatchFiles/Distances/Nei_Gojobori | 2 +- res/TemplateBatchFiles/Distances/PC | 2 +- res/TemplateBatchFiles/Distances/PC_MH | 2 +- res/TemplateBatchFiles/Distances/PC_RV | 2 +- res/TemplateBatchFiles/Distances/T3P | 2 +- res/TemplateBatchFiles/Distances/TN84 | 2 +- res/TemplateBatchFiles/Distances/TN93 | 2 +- res/TemplateBatchFiles/Distances/TN93_RV | 2 +- res/TemplateBatchFiles/Distances/p_Distance | 2 +- .../Distances/p_Distance_aa | 2 +- .../Distances/p_Distance_binary | 2 +- .../Distances/p_Distance_codon | 2 +- res/TemplateBatchFiles/F_ST.bf | 147 +++++++++--------- res/TemplateBatchFiles/libv3/all-terms.bf | 4 +- .../libv3/tasks/estimators.bf | 3 + res/TemplateBatchFiles/libv3/tasks/trees.bf | 6 +- res/TemplateBatchFiles/partitionSequences.ibf | 1 + src/core/topology.cpp | 4 +- 22 files changed, 108 insertions(+), 99 deletions(-) diff --git a/res/TemplateBatchFiles/Distances/JC69 b/res/TemplateBatchFiles/Distances/JC69 index baed47642..33a9af4dc 100644 --- a/res/TemplateBatchFiles/Distances/JC69 +++ b/res/TemplateBatchFiles/Distances/JC69 @@ -1,13 +1,9 @@ -function InitializeDistances (dummy) -{ - summingVector = {4,1}["1"]; - return 0; +function InitializeDistances (dummy) { } -function ComputeDistanceFormula (s1,s2) -{ - GetDataInfo (siteDifferenceCount, filteredData, s1, s2, DIST_OPTION); - totalSitesCompared = (Transpose(summingVector)*(siteDifferenceCount*summingVector))[0]; +function ComputeDistanceFormula (s1,s2) { + GetDataInfo (siteDifferenceCount, filteredData, s1, s2, RESOLVE_AMBIGUITIES); + totalSitesCompared = +siteDifferenceCount; if (totalSitesCompared) { diff --git a/res/TemplateBatchFiles/Distances/K2P b/res/TemplateBatchFiles/Distances/K2P index ba2672aea..f75cd2702 100644 --- a/res/TemplateBatchFiles/Distances/K2P +++ b/res/TemplateBatchFiles/Distances/K2P @@ -6,7 +6,7 @@ function InitializeDistances (dummy) function ComputeDistanceFormula (s1,s2) { - GetDataInfo (siteDifferenceCount, filteredData, s1, s2, DISTANCE_AMBIG_OPTION); + GetDataInfo (siteDifferenceCount, filteredData, s1, s2, RESOLVE_AMBIGUITIES); totalSitesCompared = (Transpose(summingVector)*(siteDifferenceCount*summingVector))[0]; if (totalSitesCompared) diff --git a/res/TemplateBatchFiles/Distances/K2P_RV b/res/TemplateBatchFiles/Distances/K2P_RV index dc03ad8e1..0bede6ecb 100644 --- a/res/TemplateBatchFiles/Distances/K2P_RV +++ b/res/TemplateBatchFiles/Distances/K2P_RV @@ -15,7 +15,7 @@ function InitializeDistances (dummy) function ComputeDistanceFormula (s1,s2) { - GetDataInfo (siteDifferenceCount, filteredData, s1, s2, DISTANCE_AMBIG_OPTION); + GetDataInfo (siteDifferenceCount, filteredData, s1, s2, RESOLVE_AMBIGUITIES); totalSitesCompared = (Transpose(summingVector)*(siteDifferenceCount*summingVector))[0]; if (totalSitesCompared) diff --git a/res/TemplateBatchFiles/Distances/Modified_Nei_Gojobori b/res/TemplateBatchFiles/Distances/Modified_Nei_Gojobori index aa8310e52..d250d865c 100644 --- a/res/TemplateBatchFiles/Distances/Modified_Nei_Gojobori +++ b/res/TemplateBatchFiles/Distances/Modified_Nei_Gojobori @@ -46,7 +46,7 @@ function ComputeDistanceFormula (s1,s2) _SITE_OS_COUNT = {stateCharCount,stateCharCount}; _SITE_ON_COUNT = {stateCharCount,stateCharCount}; - GetDataInfo (siteDifferenceCount, filteredData, s1, s2, DISTANCE_AMBIG_OPTION); + GetDataInfo (siteDifferenceCount, filteredData, s1, s2, RESOLVE_AMBIGUITIES); _SITE_OS_COUNT = matrixTrick2*(_OBSERVED_S_$siteDifferenceCount)*Transpose(matrixTrick2); _SITE_ON_COUNT = matrixTrick2*(_OBSERVED_NS_$siteDifferenceCount)*Transpose(matrixTrick2); diff --git a/res/TemplateBatchFiles/Distances/Nei_Gojobori b/res/TemplateBatchFiles/Distances/Nei_Gojobori index a9dfd3435..ae6561c99 100644 --- a/res/TemplateBatchFiles/Distances/Nei_Gojobori +++ b/res/TemplateBatchFiles/Distances/Nei_Gojobori @@ -40,7 +40,7 @@ function ComputeDistanceFormula (s1,s2) _SITE_OS_COUNT = {stateCharCount,stateCharCount}; _SITE_ON_COUNT = {stateCharCount,stateCharCount}; - GetDataInfo (siteDifferenceCount, filteredData, s1, s2, DISTANCE_AMBIG_OPTION); + GetDataInfo (siteDifferenceCount, filteredData, s1, s2, RESOLVE_AMBIGUITIES); _SITE_OS_COUNT = matrixTrick2*(_OBSERVED_S_$siteDifferenceCount)*Transpose(matrixTrick2); _SITE_ON_COUNT = matrixTrick2*(_OBSERVED_NS_$siteDifferenceCount)*Transpose(matrixTrick2); diff --git a/res/TemplateBatchFiles/Distances/PC b/res/TemplateBatchFiles/Distances/PC index 98a1aee90..7691a9034 100644 --- a/res/TemplateBatchFiles/Distances/PC +++ b/res/TemplateBatchFiles/Distances/PC @@ -7,7 +7,7 @@ function InitializeDistances (dummy) function ComputeDistanceFormula (s1,s2) { - GetDataInfo (siteDifferenceCount, filteredData, s1, s2, DISTANCE_AMBIG_OPTION); + GetDataInfo (siteDifferenceCount, filteredData, s1, s2, RESOLVE_AMBIGUITIES); totalSitesCompared = (Transpose(summingVector)*(siteDifferenceCount*summingVector)); totalDifference = Transpose(summingVector)*(siteDifferenceCount$tracingVector*summingVector); diff --git a/res/TemplateBatchFiles/Distances/PC_MH b/res/TemplateBatchFiles/Distances/PC_MH index e66f9ba9b..45fc6cb73 100644 --- a/res/TemplateBatchFiles/Distances/PC_MH +++ b/res/TemplateBatchFiles/Distances/PC_MH @@ -7,7 +7,7 @@ function InitializeDistances (dummy) function ComputeDistanceFormula (s1,s2) { - GetDataInfo (siteDifferenceCount, filteredData, s1, s2, DISTANCE_AMBIG_OPTION); + GetDataInfo (siteDifferenceCount, filteredData, s1, s2, RESOLVE_AMBIGUITIES); totalDifference = Transpose(summingVector)*(siteDifferenceCount$tracingVector*summingVector); totalSitesCompared = Transpose(summingVector)*(siteDifferenceCount*summingVector); totalDifference = totalDifference[0]/totalSitesCompared[0]; diff --git a/res/TemplateBatchFiles/Distances/PC_RV b/res/TemplateBatchFiles/Distances/PC_RV index c2a79babb..029376024 100644 --- a/res/TemplateBatchFiles/Distances/PC_RV +++ b/res/TemplateBatchFiles/Distances/PC_RV @@ -16,7 +16,7 @@ function InitializeDistances (dummy) function ComputeDistanceFormula (s1,s2) { - GetDataInfo (siteDifferenceCount, filteredData, s1, s2, DISTANCE_AMBIG_OPTION); + GetDataInfo (siteDifferenceCount, filteredData, s1, s2, RESOLVE_AMBIGUITIES); totalDifference = Transpose(summingVector)*(siteDifferenceCount$tracingVector*summingVector); totalSitesCompared = Transpose(summingVector)*(siteDifferenceCount*summingVector); return _distance_alpha*((totalDifference[0]/totalSitesCompared[0])^(-1/_distance_alpha)-1); diff --git a/res/TemplateBatchFiles/Distances/T3P b/res/TemplateBatchFiles/Distances/T3P index e13845bab..71eeb63e9 100644 --- a/res/TemplateBatchFiles/Distances/T3P +++ b/res/TemplateBatchFiles/Distances/T3P @@ -11,7 +11,7 @@ function InitializeDistances (dummy) function ComputeDistanceFormula (s1,s2) { - GetDataInfo (siteDifferenceCount, filteredData, s1, s2, DISTANCE_AMBIG_OPTION); + GetDataInfo (siteDifferenceCount, filteredData, s1, s2, RESOLVE_AMBIGUITIES); totalSitesCompared = (Transpose(summingVector)*(siteDifferenceCount*summingVector))[0]; if (totalSitesCompared) diff --git a/res/TemplateBatchFiles/Distances/TN84 b/res/TemplateBatchFiles/Distances/TN84 index 4c7bbc531..3a644dc94 100644 --- a/res/TemplateBatchFiles/Distances/TN84 +++ b/res/TemplateBatchFiles/Distances/TN84 @@ -16,7 +16,7 @@ function InitializeDistances (dummy) function ComputeDistanceFormula (s1,s2) { - GetDataInfo (siteDifferenceCount, filteredData, s1, s2, DISTANCE_AMBIG_OPTION); + GetDataInfo (siteDifferenceCount, filteredData, s1, s2, RESOLVE_AMBIGUITIES); totalSitesCompared = (Transpose(summingVector)*(siteDifferenceCount*summingVector))[0]; if (totalSitesCompared) diff --git a/res/TemplateBatchFiles/Distances/TN93 b/res/TemplateBatchFiles/Distances/TN93 index 9b1c7200a..6ee5c3f2a 100644 --- a/res/TemplateBatchFiles/Distances/TN93 +++ b/res/TemplateBatchFiles/Distances/TN93 @@ -26,7 +26,7 @@ function InitializeDistances (void) function ComputeDistanceFormula (s1,s2) { - GetDataInfo (siteDifferenceCount, filteredData, s1, s2, DEFAULT); + GetDataInfo (siteDifferenceCount, filteredData, s1, s2, RESOLVE_AMBIGUITIES); totalSitesCompared = (Transpose(summingVector)*(siteDifferenceCount*summingVector))[0]; if (totalSitesCompared) diff --git a/res/TemplateBatchFiles/Distances/TN93_RV b/res/TemplateBatchFiles/Distances/TN93_RV index 437bb36df..985c5946e 100644 --- a/res/TemplateBatchFiles/Distances/TN93_RV +++ b/res/TemplateBatchFiles/Distances/TN93_RV @@ -25,7 +25,7 @@ function InitializeDistances (dummy) function ComputeDistanceFormula (s1,s2) { - GetDataInfo (siteDifferenceCount, filteredData, s1, s2, DISTANCE_AMBIG_OPTION); + GetDataInfo (siteDifferenceCount, filteredData, s1, s2, RESOLVE_AMBIGUITIES); totalSitesCompared = (Transpose(summingVector)*(siteDifferenceCount*summingVector))[0]; if (totalSitesCompared) diff --git a/res/TemplateBatchFiles/Distances/p_Distance b/res/TemplateBatchFiles/Distances/p_Distance index 2a4d80c12..313abec2d 100644 --- a/res/TemplateBatchFiles/Distances/p_Distance +++ b/res/TemplateBatchFiles/Distances/p_Distance @@ -6,7 +6,7 @@ function InitializeDistances (dummy) function ComputeDistanceFormula (s1,s2) { - GetDataInfo (siteDifferenceCount, filteredData, s1, s2, DISTANCE_AMBIG_OPTION); + GetDataInfo (siteDifferenceCount, filteredData, s1, s2, RESOLVE_AMBIGUITIES); totalSitesCompared = (Transpose(summingVector)*(siteDifferenceCount*summingVector))[0]; if (totalSitesCompared) { diff --git a/res/TemplateBatchFiles/Distances/p_Distance_aa b/res/TemplateBatchFiles/Distances/p_Distance_aa index 2e900be34..4335d19a2 100644 --- a/res/TemplateBatchFiles/Distances/p_Distance_aa +++ b/res/TemplateBatchFiles/Distances/p_Distance_aa @@ -12,7 +12,7 @@ function InitializeDistances (dummy) function ComputeDistanceFormula (s1,s2) { - GetDataInfo (siteDifferenceCount, filteredData, s1, s2, DISTANCE_AMBIG_OPTION); + GetDataInfo (siteDifferenceCount, filteredData, s1, s2, RESOLVE_AMBIGUITIES); totalDifference = Transpose(summingVector)*(siteDifferenceCount$tracingVector*summingVector); totalSitesCompared = Transpose(summingVector)*(siteDifferenceCount*summingVector); return 1-totalDifference[0]/totalSitesCompared[0]; diff --git a/res/TemplateBatchFiles/Distances/p_Distance_binary b/res/TemplateBatchFiles/Distances/p_Distance_binary index 761fa86d2..3b4adba09 100644 --- a/res/TemplateBatchFiles/Distances/p_Distance_binary +++ b/res/TemplateBatchFiles/Distances/p_Distance_binary @@ -6,7 +6,7 @@ function InitializeDistances (dummy) function ComputeDistanceFormula (s1,s2) { - GetDataInfo (siteDifferenceCount, filteredData, s1, s2, DISTANCE_AMBIG_OPTION); + GetDataInfo (siteDifferenceCount, filteredData, s1, s2, RESOLVE_AMBIGUITIES); totalSitesCompared = (Transpose(summingVector)*(siteDifferenceCount*summingVector))[0]; if (totalSitesCompared) { diff --git a/res/TemplateBatchFiles/Distances/p_Distance_codon b/res/TemplateBatchFiles/Distances/p_Distance_codon index ab8271136..129e38128 100644 --- a/res/TemplateBatchFiles/Distances/p_Distance_codon +++ b/res/TemplateBatchFiles/Distances/p_Distance_codon @@ -34,7 +34,7 @@ function ComputeDistanceFormula (s1,s2) { _SITE_ON_COUNT = {stateCharCount,stateCharCount}; - GetDataInfo (siteDifferenceCount, filteredData, s1, s2, DISTANCE_AMBIG_OPTION); + GetDataInfo (siteDifferenceCount, filteredData, s1, s2, RESOLVE_AMBIGUITIES); if (distChoice == 0) diff --git a/res/TemplateBatchFiles/F_ST.bf b/res/TemplateBatchFiles/F_ST.bf index b8110e898..0b2e30a0e 100644 --- a/res/TemplateBatchFiles/F_ST.bf +++ b/res/TemplateBatchFiles/F_ST.bf @@ -1,33 +1,28 @@ -RequireVersion ("2.2"); +RequireVersion ("2.5.25"); LoadFunctionLibrary ("GrabBag"); +LoadFunctionLibrary("libv3/UtilityFunctions.bf"); +LoadFunctionLibrary("libv3/IOFunctions.bf"); +LoadFunctionLibrary("libv3/all-terms.bf"); + /*------------------------------------------------------------------------------*/ -function doSNN (vpart1, vpart2, vsize1, vsize2) -{ - for (k=0; k" + namesMatrix [_i] + "\nA\n"); @@ -185,7 +180,11 @@ if (useSeqData) } promptFor2ndRegExp = 1; +KeywordArgument ("regexp1", "Regular expression for the first sequence partition"); +KeywordArgument ("regexp2", "Regular expression for the second sequence partition"); +skipConfirm = 1; ExecuteAFile ("partitionSequences.ibf"); +skipConfirm = 0; promptFor2ndRegExp = 0; bothCladesSize = clASize + clBSize; @@ -221,19 +220,18 @@ for (specIndex = 0; specIndex < ds.species; specIndex += 1) { } } -if (dataType) { +if (dataType == "Codon") { DataSetFilter filteredData = CreateFilter (ds,3,"","",GeneticCodeExclusions); } else { DataSetFilter filteredData = CreateFilter (ds,1,"",""); } -if (useSeqData) { +if (useSeqData != "Included") { distanceMatrix = {ds.species, ds.species}; } -if (distanceChoice == 1) -{ +if (distanceChoice == "Full likelihood") { SelectTemplateModel(filteredData); fprintf (stdout,"\nHYPHY is computing pairwise maximum likelihood distance estimates. A total of ", Format(ds.species*(ds.species-1)/2,0,0), @@ -296,9 +294,9 @@ if (distanceChoice == 1) } else { - if (distanceChoice == 2) + if (distanceChoice == "Load Matrix") { - if (useSeqData) { + if (useSeqData != "Included") { SetDialogPrompt ("Load the distance matrix"); fscanf (PROMPT_FOR_FILE,"NMatrix",distanceMatrix); } @@ -308,9 +306,11 @@ else return 0; } } - else - { + else { DISTANCE_PROMPTS = 1; + KeywordArgument ("formula", "Distance Formula", "JC69"); + + ExecuteAFile (HYPHY_LIB_DIRECTORY + "TemplateBatchFiles" + DIRECTORY_SEPARATOR + "chooseDistanceFormula.def"); InitializeDistances (0); @@ -325,8 +325,7 @@ else for (i = 0; i 0) -{ +if (resample > 0) { + KeywordArgument ("permutation-samples", "The number of permutations to perform", 100); sampleCount = prompt_for_a_value ("The number of permutations to perform", 100, 1, 100000, 1); diff --git a/res/TemplateBatchFiles/libv3/all-terms.bf b/res/TemplateBatchFiles/libv3/all-terms.bf index 98d164778..9d4382057 100644 --- a/res/TemplateBatchFiles/libv3/all-terms.bf +++ b/res/TemplateBatchFiles/libv3/all-terms.bf @@ -239,7 +239,7 @@ namespace terms{ NTP = "NTP"; SS = "SS"; NS = "NS"; - } + } /* Terms used in I/O */ namespace io { @@ -314,6 +314,8 @@ namespace terms{ version = "version"; convergence_failures = "convergence failures"; omega_ratio = "omega"; + dS = "dS"; + dN = "dN"; rate = "rate"; proportion = "proportion"; positive = "positive test results"; diff --git a/res/TemplateBatchFiles/libv3/tasks/estimators.bf b/res/TemplateBatchFiles/libv3/tasks/estimators.bf index 5f86ca40e..044292164 100644 --- a/res/TemplateBatchFiles/libv3/tasks/estimators.bf +++ b/res/TemplateBatchFiles/libv3/tasks/estimators.bf @@ -435,7 +435,10 @@ function estimators.ExtractMLEsOptions(likelihood_function_id, model_description GetInformation (estimators.ExtractMLEs.map, *_tree_name); + estimators.ExtractMLEs.branch_names = Rows(estimators.ExtractMLEs.map); + + (estimators.ExtractMLEs.results[terms.branch_length])[estimators.ExtractMLEs.i] = {}; for (estimators.ExtractMLEs.b = 0; estimators.ExtractMLEs.b < Abs(estimators.ExtractMLEs.map); estimators.ExtractMLEs.b += 1) { diff --git a/res/TemplateBatchFiles/libv3/tasks/trees.bf b/res/TemplateBatchFiles/libv3/tasks/trees.bf index 7bfe449b4..4ec46d2c8 100644 --- a/res/TemplateBatchFiles/libv3/tasks/trees.bf +++ b/res/TemplateBatchFiles/libv3/tasks/trees.bf @@ -447,8 +447,9 @@ lfunction trees.ExtractTreeInfoFromTopology(topology_object) { branch_lengths = BranchLength(^topology_object, -1); branch_names = BranchName(^topology_object, -1); - branch_count = utility.Array1D (branch_names) - 1; - + branch_count = Max (2,utility.Array1D (branch_names) - 1); + + bls = {}; for (k = 0; k < branch_count; k+=1) { @@ -456,6 +457,7 @@ lfunction trees.ExtractTreeInfoFromTopology(topology_object) { bls [branch_names[k]] = branch_lengths[k]; } } + GetInformation(modelMap, ^topology_object); diff --git a/res/TemplateBatchFiles/partitionSequences.ibf b/res/TemplateBatchFiles/partitionSequences.ibf index 9b956ea6b..df5736004 100644 --- a/res/TemplateBatchFiles/partitionSequences.ibf +++ b/res/TemplateBatchFiles/partitionSequences.ibf @@ -59,6 +59,7 @@ while (goOn) { fprintf (stdout, "\n\nClade 2 includes ", clBSize," sequences:\n"); cladeB ["__echo_clades"][""]; + if (skipConfirm) {break;} fprintf (stdout, "\nIs this partitioning correct (y/n)"); fscanf (stdin, "String", goOn); goOn = (goOn[0] == "n" || goOn[0] == "N"); diff --git a/src/core/topology.cpp b/src/core/topology.cpp index 5732a93d4..bd40651a0 100644 --- a/src/core/topology.cpp +++ b/src/core/topology.cpp @@ -2366,8 +2366,10 @@ HBLObjectRef _TreeTopology::BranchLength (HBLObjectRef p, HBLObjectRef cache) { // get ALL branch lengths _Vector * branch_lengths = new _Vector; + bool two = IsDegenerate(); + while (node* iterator = ni.Next()) { - if (!iterator->is_root()){ + if (!iterator->is_root() || two){ branch_lengths->Store(GetBranchLength (iterator)); } } From 462a9a6ec09da8f9f45d474bf5f8c50f4962515a Mon Sep 17 00:00:00 2001 From: Sergei Date: Sun, 7 Feb 2021 11:49:01 -0500 Subject: [PATCH 4/4] Bug fixes and MSS updates --- .../libv3/models/codon/MSS.bf | 124 ++++++++++++++++++ src/core/batchlanruntime.cpp | 2 +- src/core/dataset.cpp | 2 +- src/core/global_things.cpp | 2 +- src/core/likefunc.cpp | 2 +- src/core/string_buffer.cpp | 2 + 6 files changed, 130 insertions(+), 4 deletions(-) create mode 100644 res/TemplateBatchFiles/libv3/models/codon/MSS.bf diff --git a/res/TemplateBatchFiles/libv3/models/codon/MSS.bf b/res/TemplateBatchFiles/libv3/models/codon/MSS.bf new file mode 100644 index 000000000..2b63bf196 --- /dev/null +++ b/res/TemplateBatchFiles/libv3/models/codon/MSS.bf @@ -0,0 +1,124 @@ +LoadFunctionLibrary("MG_REV.bf"); + +/** @module models.codon.MG_REV */ + +/** + * @name models.codon.MSS.ModelDescription + * @param {String} type + * @param {String} code + * @param {Dict} codon_classes : codon => neutral class + */ + +terms.model.MSS.syn_rate_within = " within codon class "; +terms.model.MSS.syn_rate_between = " between codon classes "; + + +lfunction models.codon.MSS.ModelDescription(type, code, codon_classes) { + m = Call ("models.codon.MG_REV.ModelDescription", type, code); + if (^"fitter.frequency_type" == "F3x4") { + m[utility.getGlobalValue("terms.model.frequency_estimator")] = "frequencies.empirical.F3x4"; + } else { + if (^"fitter.frequency_type" == "F1x4") { + m[utility.getGlobalValue("terms.model.frequency_estimator")] = "frequencies.empirical.F1x4"; + } + } + m[utility.getGlobalValue("terms.description")] = "The Muse-Gaut 94 codon-substitution model coupled with the general time reversible (GTR) model of nucleotide substitution, which allows multiple classes of synonymous substitution rates"; + m[utility.getGlobalValue("terms.model.q_ij")] = "models.codon.MSS._GenerateRate"; + m[utility.getGlobalValue("terms.model.mss.codon_classes")] = codon_classes; + return m; +} + +//---------------------------------------------------------------------------------------------------------------- + +lfunction models.codon.MSS._GenerateRate (fromChar, toChar, namespace, model_type, model) { + + _GenerateRate.p = {}; + _GenerateRate.diff = models.codon.diff.complete(fromChar, toChar); + diff_count = utility.Array1D (_GenerateRate.diff); + + omega_term = utility.getGlobalValue ("terms.parameters.omega_ratio"); + alpha_term = utility.getGlobalValue ("terms.parameters.synonymous_rate"); + beta_term = utility.getGlobalValue ("terms.parameters.nonsynonymous_rate"); + omega = "omega"; + alpha = "alpha"; + beta = "beta"; + + _tt = model[utility.getGlobalValue("terms.translation_table")]; + + if (diff_count == 1) { + + _GenerateRate.p[model_type] = {}; + _GenerateRate.p[utility.getGlobalValue("terms.global")] = {}; + + nuc_rate = ""; + + for (i = 0; i < diff_count; i += 1) { + if ((_GenerateRate.diff[i])[utility.getGlobalValue("terms.diff.from")] > (_GenerateRate.diff[i])[utility.getGlobalValue("terms.diff.to")]) { + nuc_p = "theta_" + (_GenerateRate.diff[i])[utility.getGlobalValue("terms.diff.to")] + (_GenerateRate.diff[i])[utility.getGlobalValue("terms.diff.from")]; + } else { + nuc_p = "theta_" + (_GenerateRate.diff[i])[utility.getGlobalValue("terms.diff.from")] +(_GenerateRate.diff[i])[utility.getGlobalValue("terms.diff.to")]; + } + nuc_p = parameters.ApplyNameSpace(nuc_p, namespace); + (_GenerateRate.p[utility.getGlobalValue("terms.global")])[terms.nucleotideRateReversible((_GenerateRate.diff[i])[utility.getGlobalValue("terms.diff.from")], (_GenerateRate.diff[i])[utility.getGlobalValue("terms.diff.to")])] = nuc_p; + + nuc_rate = parameters.AppendMultiplicativeTerm (nuc_rate, nuc_p); + } + + + rate_entry = nuc_rate; + + if (_tt[fromChar] != _tt[toChar]) { + + if (model_type == utility.getGlobalValue("terms.global")) { + aa_rate = parameters.ApplyNameSpace(omega, namespace); + (_GenerateRate.p[model_type])[omega_term] = aa_rate; + } else { + aa_rate = beta; + (_GenerateRate.p[model_type])[beta_term] = aa_rate; + } + rate_entry += "*" + aa_rate; + } else { + + class_from = (model[^"terms.model.mss.codon_classes"])[fromChar]; + class_to = (model[^"terms.model.mss.codon_classes"])[toChar]; + + if (class_from == class_to) { + if (class_from == ^"mss.neutral_reference") { + if (model_type == utility.getGlobalValue("terms.local")) { + codon_rate = alpha + "_" + class_from; + (_GenerateRate.p[model_type])[alpha_term + ^"terms.model.MSS.syn_rate_within" + class_from] = codon_rate; + rate_entry += "*" + codon_rate; + } else { + rate_entry = nuc_rate; + } + } else { + if (model_type == utility.getGlobalValue("terms.local")) { + codon_rate = alpha + "_" + class_from; + } else { + codon_rate = parameters.ApplyNameSpace(alpha + "_" + class_from, namespace); + } + (_GenerateRate.p[model_type])[alpha_term + ^"terms.model.MSS.syn_rate_within" + class_from] = codon_rate; + rate_entry += "*" + codon_rate; + } + } else { + if (class_from > class_to) { + codon_rate = class_to; + class_to = class_from; + class_from = codon_rate; + } + if (model_type == utility.getGlobalValue("terms.local")) { + codon_rate = alpha + "_" + class_from + "_" + class_to; + } else { + codon_rate = parameters.ApplyNameSpace(alpha + "_" + class_from + "_" + class_to, namespace); + } + (_GenerateRate.p[model_type])[alpha_term + ^"terms.model.MSS.syn_rate_between" + class_from + " and " + class_to] = codon_rate; + rate_entry += "*" + codon_rate; + } + } + + _GenerateRate.p[utility.getGlobalValue("terms.model.rate_entry")] = rate_entry; + } + + return _GenerateRate.p; +} + diff --git a/src/core/batchlanruntime.cpp b/src/core/batchlanruntime.cpp index 5f781ec9d..069a7899c 100644 --- a/src/core/batchlanruntime.cpp +++ b/src/core/batchlanruntime.cpp @@ -1875,7 +1875,7 @@ bool _ElementaryCommand::HandleAdvanceIterator(_ExecutionList& current_prog } } - if (row >= source_object->GetHDim()) { + if (row >= source_object->GetHDim() || source_object->GetVDim() == 0) { ((_Variable*)parameters.GetItem (reciever_count << 1))->SetValue (new _FString (kEndIteration), false, false, NULL); // end loop } else { ((_Variable*)parameters.GetItem (reciever_count << 1))->SetValue (source_object->GetMatrixCell(row, column), false, false, NULL); diff --git a/src/core/dataset.cpp b/src/core/dataset.cpp index 860fa9d24..09c8b32a1 100644 --- a/src/core/dataset.cpp +++ b/src/core/dataset.cpp @@ -1614,7 +1614,7 @@ void FilterRawString (_String& s, FileState* fs, _DataSet & ds) { //_________________________________________________________________________________________________ -void ProcessTree (FileState *fState, FILE* f, _String& CurrentLine) { +void ProcessTree (FileState *fState, FILE* f, _StringBuffer& CurrentLine) { // TODO SLKP 20180921 this does extra work to read in the tree string multiple times; // the solution is to have a proper buffer wrapper, and to diff --git a/src/core/global_things.cpp b/src/core/global_things.cpp index 57db2e560..a2fb36a84 100644 --- a/src/core/global_things.cpp +++ b/src/core/global_things.cpp @@ -121,7 +121,7 @@ namespace hy_global { kErrorStringDatasetRefIndexError ("Dataset index reference out of range"), kErrorStringMatrixExportError ("Export matrix called with a non-polynomial matrix argument"), kErrorStringNullOperand ("Attempting to operate on an undefined value; this is probably the result of an earlier 'soft' error condition"), - kHyPhyVersion = _String ("2.5.27"), + kHyPhyVersion = _String ("2.5.28"), kNoneToken = "None", kNullToken = "null", diff --git a/src/core/likefunc.cpp b/src/core/likefunc.cpp index c67efc222..13e4674d9 100644 --- a/src/core/likefunc.cpp +++ b/src/core/likefunc.cpp @@ -7061,7 +7061,7 @@ void _LikelihoodFunction::LocateTheBump (long index,hyFloat gPrecision, hyFlo } if (CheckEqual(GetIthIndependent(index), bestVal) && fabs (middleValue-maxSoFar) > 1e-12*errorTolerance) { char buf[256]; - snprintf (buf, 256, " \n\tERROR: [_LikelihoodFunction::LocateTheBump (index %ld) current value %20.16g (parameter = %20.16g), best value %20.16g (parameter = %20.16g)); delta = %20.16g, tolerance = %20.16]\n\n", index, middleValue, GetIthIndependent(index), maxSoFar, bestVal, maxSoFar - middleValue, 1e-12*errorTolerance); + snprintf (buf, 256, " \n\tERROR: [_LikelihoodFunction::LocateTheBump (index %ld) current value %20.16g (parameter = %20.16g), best value %20.16g (parameter = %20.16g)); delta = %20.16g, tolerance = %20.16g]\n\n", index, middleValue, GetIthIndependent(index), maxSoFar, bestVal, maxSoFar - middleValue, 1e-12*errorTolerance); _TerminateAndDump (_String (buf) & "\n" & "\nParameter name " & *GetIthIndependentName(index)); } SetIthIndependent(index,bestVal); diff --git a/src/core/string_buffer.cpp b/src/core/string_buffer.cpp index 0e8453a55..26588c592 100644 --- a/src/core/string_buffer.cpp +++ b/src/core/string_buffer.cpp @@ -250,6 +250,8 @@ _StringBuffer& _StringBuffer::operator = (_StringBuffer && rhs) { return *this; } + + //=============================================================