1 /*
2
3 HyPhy - Hypothesis Testing Using Phylogenies.
4
5 Copyright (C) 1997-now
6 Core Developers:
7 Sergei L Kosakovsky Pond (sergeilkp@icloud.com)
8 Art FY Poon (apoon42@uwo.ca)
9 Steven Weaver (sweaver@temple.edu)
10
11 Module Developers:
12 Lance Hepler (nlhepler@gmail.com)
13 Martin Smith (martin.audacis@gmail.com)
14
15 Significant contributions from:
16 Spencer V Muse (muse@stat.ncsu.edu)
17 Simon DW Frost (sdf22@cam.ac.uk)
18
19 Permission is hereby granted, free of charge, to any person obtaining a
20 copy of this software and associated documentation files (the
21 "Software"), to deal in the Software without restriction, including
22 without limitation the rights to use, copy, modify, merge, publish,
23 distribute, sublicense, and/or sell copies of the Software, and to
24 permit persons to whom the Software is furnished to do so, subject to
25 the following conditions:
26
27 The above copyright notice and this permission notice shall be included
28 in all copies or substantial portions of the Software.
29
30 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
31 OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
32 MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
33 IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
34 CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
35 TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
36 SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
37
38 */
39
40 #include "baseobj.h"
41 #include "alignment.h"
42 #include "defines.h"
43 #include "batchlan.h"
44 #include "likefunc.h"
45 #include "bayesgraph.h"
46 #include "scfg.h"
47 #include "global_object_lists.h"
48 #include "mersenne_twister.h"
49 #include "global_things.h"
50 #include "hy_string_buffer.h"
51 #include "associative_list.h"
52 #include "tree_iterator.h"
53
54 #include "function_templates.h"
55
56 #include <signal.h>
57
58 #ifndef __HYPHY_NO_SQLITE__
59 #include "sqlite3.h"
60 #endif
61
62 using namespace hy_global;
63 using namespace hyphy_global_objects;
64
65 #include <ctype.h>
66 //____________________________________________________________________________________
67 /* various helper functions */
68
ExtractStatementAssignment(_String const & source,long & end_at,const bool validate,const bool exceptions,const long offset)69 const _String _ElementaryCommand::ExtractStatementAssignment (_String const& source, long& end_at, const bool validate, const bool exceptions, const long offset) {
70
71 _String id;
72
73 try {
74 long id_start = source.FirstNonSpaceFollowingSpace(offset);
75
76 end_at = id_start != kNotFound ? source.FindTerminator(id_start, '=') : kNotFound;
77
78 if (id_start == kNotFound || end_at == kNotFound) {
79 throw _String ("Missing an ID in 'Type <ID> = statement'");
80 }
81
82 id = source.Cut (id_start, end_at - 1);
83 if (validate) {
84 if (!id.IsValidIdentifier(fIDAllowCompound)) {
85 throw id.Enquote() & "is not a valid storage variable identifier";
86 }
87 }
88
89 end_at ++;
90
91 } catch (const _String& err) {
92 if (exceptions) {
93 throw err;
94 }
95 id.Clear();
96 end_at = kNotFound;
97 }
98
99 return id;
100
101 }
102
103 //____________________________________________________________________________________
104
ProcessProcedureCall(_String const & source,_String & procedure,_List & pieces)105 const _String _ElementaryCommand::ProcessProcedureCall (_String const& source,_String& procedure, _List& pieces) {
106 long op_start;
107 _String id = ExtractStatementAssignment (source, op_start, false);
108
109 long paren_start = op_start,
110 paren_end = source.ExtractEnclosedExpression(paren_start, '(', ')', fExtractRespectQuote | fExtractRespectEscape);
111
112 if (paren_end == kNotFound) {
113 throw _String ("Missing () enclosed argument list");
114 }
115
116 procedure = source.Cut (op_start,paren_start-1L);
117 pieces < new _String (id);
118 ExtractConditions (source,paren_start+1,pieces,',');
119
120 return id;
121 }
122
123 //____________________________________________________________________________________
124
_CheckExpressionForCorrectness(_Formula & parsed_expression,_String const & exp,_ExecutionList & program,long desired_type=HY_ANY_OBJECT)125 void _CheckExpressionForCorrectness (_Formula& parsed_expression, _String const& exp, _ExecutionList& program, long desired_type = HY_ANY_OBJECT) {
126 _String error_message;
127
128 long parse_result = parsed_expression.ParseFormula (exp,program.nameSpacePrefix, &error_message);
129
130 if (error_message.nonempty()) {
131 throw (_String ("Failed to parse ") & exp.Enquote () & " with the following error: " & error_message);
132 }
133 if (parse_result != HY_FORMULA_EXPRESSION) {
134 throw (exp.Enquote () & " did not parse to a simple expression");
135 }
136 if (parsed_expression.IsEmpty ()) {
137 throw (exp.Enquote () & " parsed to an empty expression");
138 }
139 if (!(desired_type == HY_ANY_OBJECT || parsed_expression.ObjectClass() & desired_type)) {
140 // TODO SLKP 20170704: ObjectClass will compute the expression with current values which may fail
141 throw (exp.Enquote () & " did not evaluate to a " & FetchObjectNameFromType (desired_type));
142 }
143 }
144
145 //____________________________________________________________________________________
146
_CheckForExistingVariableByType(_String const & name,_ExecutionList & program,long desired_type=HY_ANY_OBJECT)147 _Variable* _CheckForExistingVariableByType (_String const& name, _ExecutionList& program, long desired_type = HY_ANY_OBJECT) {
148 _String variable_id = AppendContainerName(name,program.nameSpacePrefix);
149 _Variable * target_variable = FetchVar(LocateVarByName (variable_id));
150
151 if (!target_variable) {
152 throw (variable_id.Enquote() & " is not an existing variable");
153 }
154
155 if (!(desired_type == HY_ANY_OBJECT || target_variable->ObjectClass() & desired_type)) {
156 throw (name.Enquote () & " is not of type " & FetchObjectNameFromType (desired_type));
157 }
158
159 return target_variable;
160 }
161
162 //____________________________________________________________________________________
163
_ProcessAnArgumentByType(_String const & expression,long desired_type,_ExecutionList & program,_List * reference_manager)164 HBLObjectRef _ProcessAnArgumentByType (_String const& expression, long desired_type, _ExecutionList& program, _List* reference_manager) {
165 // The return value needs to managed by the caller
166
167 /* first see if this is a simple expression of the form 'variable_id' */
168
169 HBLObjectRef simple_var = FetchObjectFromVariableByType (&AppendContainerName (expression, program.nameSpacePrefix), desired_type);
170 if (simple_var) {
171 if (reference_manager)
172 *reference_manager << simple_var;
173 else
174 simple_var->AddAReference();
175 return simple_var;
176 }
177
178 _Formula parsed_expression;
179 _CheckExpressionForCorrectness (parsed_expression, expression, program, desired_type);
180
181 HBLObjectRef expression_result = parsed_expression.Compute(0,program.nameSpacePrefix);
182 if (expression_result && (expression_result->ObjectClass() & desired_type)) {
183 if (reference_manager)
184 *reference_manager << expression_result;
185 else
186 expression_result->AddAReference();
187 return expression_result;
188 }
189
190 throw (expression.Enquote () & " did not evaluate to a " & FetchObjectNameFromType (desired_type));
191 return nil;
192 }
193
194 //____________________________________________________________________________________
195
_ProcessALiteralArgument(_String const & expression,_ExecutionList & program)196 const _String _ProcessALiteralArgument (_String const& expression, _ExecutionList& program) {
197 HBLObjectRef the_string = _ProcessAnArgumentByType (expression, STRING, program, nil);
198
199 _String result (((_FString*)the_string)->get_str());
200 DeleteObject (the_string);
201 return result;
202 }
203
204 //____________________________________________________________________________________
205
_GetHBLObjectByType(_String const & source_name,long & type,long * object_index,_ExecutionList * current_program)206 BaseRefConst _GetHBLObjectByType (_String const& source_name, long& type, long * object_index, _ExecutionList* current_program) {
207 long object_type = type;
208 BaseRefConst source_object = _HYRetrieveBLObjectByName (source_name, object_type,object_index,false, true, current_program);
209
210 if (source_object == nil) {
211 throw (source_name.Enquote('\'') & " is not a " & _HYHBLTypeToText(type));
212 }
213 type = object_type;
214 return source_object;
215 }
216
217 //____________________________________________________________________________________
218
_GetHBLObjectByTypeMutable(_String const & source_name,long & type,long * object_index=nil,bool do_literal_lookup=true)219 BaseRef _GetHBLObjectByTypeMutable (_String const& source_name, long& type, long * object_index = nil, bool do_literal_lookup = true) {
220 long object_type = type;
221 BaseRef source_object = _HYRetrieveBLObjectByNameMutable (source_name, object_type,object_index,false, do_literal_lookup);
222
223 if (source_object == nil) {
224 throw (source_name.Enquote('\'') & " is not a " & _HYHBLTypeToText(type));
225 }
226 type = object_type;
227 return source_object;
228 }
229
230
231 //____________________________________________________________________________________
232
_ValidateStorageVariable(_ExecutionList & program,unsigned long argument_index) const233 _Variable* _ElementaryCommand::_ValidateStorageVariable (_ExecutionList& program, unsigned long argument_index) const {
234 _String storage_id (program.AddNameSpaceToID(*GetIthParameter(argument_index)));
235 _Variable * receptacle = CheckReceptacleCommandIDException (&AppendContainerName(storage_id,program.nameSpacePrefix),get_code(), true, false, &program);
236 return receptacle;
237 }
238
239 //____________________________________________________________________________________
240
_DefaultExceptionHandler(_Variable * receptacle,_String const & error,_ExecutionList & current_program)241 bool _DefaultExceptionHandler (_Variable * receptacle, _String const& error, _ExecutionList& current_program) {
242 if (receptacle) { // if receptacle is nil, then we have already handled the error
243 receptacle->SetValue(new _MathObject, false, true, NULL);
244 }
245 current_program.ReportAnExecutionError (error);
246 return false;
247 }
248
249
250 //____________________________________________________________________________________
251
_EnsurePresenceOfKey(_AssociativeList * dict,_String const & key,long desired_type)252 HBLObjectRef _EnsurePresenceOfKey (_AssociativeList * dict, _String const& key, long desired_type) {
253 HBLObjectRef value = dict->GetByKey (key, desired_type);
254 if (!value) {
255 throw (key.Enquote() & " was not a key associated with a " & FetchObjectNameFromType (desired_type) & "-typed value");
256 }
257 return value;
258 }
259
260 //____________________________________________________________________________________
261
_NumericValueFromKey(_AssociativeList * dict,_String const & key,hyFloat default_value)262 hyFloat _NumericValueFromKey (_AssociativeList * dict, _String const& key, hyFloat default_value) {
263 HBLObjectRef value = dict->GetByKey (key, NUMBER);
264 if (!value) {
265 return default_value;
266 }
267 return value->Compute()->Value();
268 }
269
270 //____________________________________________________________________________________
271
HandleDifferentiate(_ExecutionList & current_program)272 bool _ElementaryCommand::HandleDifferentiate(_ExecutionList& current_program){
273 _Variable * receptacle = nil;
274 current_program.advance ();
275
276 try {
277 receptacle = _ValidateStorageVariable (current_program);
278 _String expression = *GetIthParameter(1);
279
280 _Formula parsed_expression;
281 _CheckExpressionForCorrectness (parsed_expression, expression, current_program);
282
283 long times = 1L;
284 if (parameter_count() >= 4UL) {
285 times = _ProcessNumericArgumentWithExceptions (*GetIthParameter(3),current_program.nameSpacePrefix);
286 if (times <= 0L) {
287 throw (GetIthParameter(3UL)->Enquote() & " (the number of times to differentiate) must be a non-negative integer");
288 }
289 }
290
291 _Variable * dx = _ValidateStorageVariable (current_program, 2);
292 _Formula * derivative = parsed_expression.Differentiate(*dx->GetName());
293
294 for (; times>1 && derivative; times--) {
295 _Formula * temp = derivative->Differentiate (*GetIthParameter(2));
296 delete derivative;
297 derivative = temp;
298 }
299 if (derivative) {
300 receptacle->SetFormula(*derivative);
301 delete derivative;
302 } else {
303 throw (_String ("Differentiation of ") & _String((_String*)GetIthParameter(1)->toStr()).Enquote() & " failed.");
304 }
305
306 } catch (const _String& error) {
307 return _DefaultExceptionHandler (receptacle, error, current_program);
308 }
309 return true;
310 }
311
312 //____________________________________________________________________________________
313
HandleFindRootOrIntegrate(_ExecutionList & currentProgram,bool do_integrate)314 bool _ElementaryCommand::HandleFindRootOrIntegrate (_ExecutionList& currentProgram, bool do_integrate){
315
316 _Variable * receptacle = nil;
317 currentProgram.advance ();
318
319 try {
320 receptacle = _ValidateStorageVariable (currentProgram);
321 _String expression = *GetIthParameter(1);
322
323 _Formula parsed_expression;
324 _CheckExpressionForCorrectness (parsed_expression, expression, currentProgram);
325 _ValidateStorageVariable (currentProgram, 2);
326 _Variable * target_variable = _CheckForExistingVariableByType (*GetIthParameter(2),currentProgram,NUMBER);
327
328 if (!parsed_expression.DependsOnVariable(target_variable->get_index()) && !do_integrate) {
329 throw (expression.Enquote() & " does not depend on the variable " & target_variable->GetName()->Enquote());
330 }
331
332
333 hyFloat lb = _ProcessNumericArgumentWithExceptions (*GetIthParameter(3),currentProgram.nameSpacePrefix),
334 ub = _ProcessNumericArgumentWithExceptions (*GetIthParameter(4),currentProgram.nameSpacePrefix);
335
336 if (ub<=lb) {
337 throw (_String ('[') & lb & ',' & ub & "] is not a valid interval");
338 return false;
339 }
340
341 _Formula * derivative = do_integrate ? nil : parsed_expression.Differentiate (*target_variable->GetName(),false);
342
343 if (!do_integrate) {
344 if (derivative) {
345 receptacle->SetValue (new _Constant (parsed_expression.Newton (*derivative,target_variable, 0.0, lb, ub)),false,true, NULL);
346 } else {
347 receptacle->SetValue (new _Constant (parsed_expression.Brent (target_variable, lb, ub)), false,true, NULL);
348 }
349 } else {
350 receptacle->SetValue (new _Constant (parsed_expression.Integral (target_variable, lb, ub, ub-lb>100)), false, true, NULL);
351 }
352
353 if (derivative) {
354 delete derivative;
355 }
356 } catch (const _String& error) {
357 return _DefaultExceptionHandler (receptacle, error, currentProgram);
358 }
359 return true;
360 }
361
362 //____________________________________________________________________________________
363
HandleExport(_ExecutionList & current_program)364 bool _ElementaryCommand::HandleExport(_ExecutionList& current_program){
365
366 _Variable * receptacle = nil;
367 current_program.advance();
368
369 try {
370 receptacle = _ValidateStorageVariable (current_program);
371
372 const _String source_name = AppendContainerName (*GetIthParameter(1), current_program.nameSpacePrefix);
373 long object_type = HY_BL_MODEL | HY_BL_LIKELIHOOD_FUNCTION | HY_BL_DATASET_FILTER | HY_BL_HBL_FUNCTION,
374 object_index;
375
376 BaseRef source_object;
377 try {
378 source_object = _GetHBLObjectByTypeMutable (source_name, object_type, &object_index);
379 } catch (const _String& ) {
380 receptacle->SetValue(new _MathObject, false, true, NULL);
381 }
382
383
384 switch (object_type) {
385 case HY_BL_LIKELIHOOD_FUNCTION: {
386 _StringBuffer * serialized_object = new _StringBuffer (8192L);
387 ((_LikelihoodFunction*)source_object)->SerializeLF (*serialized_object);
388 receptacle->SetValue(new _FString (serialized_object), false, true, NULL);
389 break;
390 }
391 case HY_BL_DATASET_FILTER: {
392 receptacle->SetValue(new _FString (new _String ((_String*)((_DataSetFilter*)source_object)->toStr())), false, true, NULL);
393 ReleaseDataFilterLock(object_index);
394 break;
395 }
396 case HY_BL_MODEL: {
397 _StringBuffer * serialized_object = new _StringBuffer (8192L);
398 SerializeModel (*serialized_object,object_index,nil,true);
399 receptacle->SetValue(new _FString (serialized_object), false, true, NULL);
400 break;
401 }
402 case HY_BL_HBL_FUNCTION: {
403 receptacle->SetValue(new _FString (new _String (ExportBFFunction (object_index))), false, true, NULL);
404 break;
405 }
406 }
407
408 } catch (const _String& error) {
409 return _DefaultExceptionHandler (receptacle, error, current_program);
410 }
411 return true;
412 }
413
414 //____________________________________________________________________________________
415
HandleGetDataInfo(_ExecutionList & current_program)416 bool _ElementaryCommand::HandleGetDataInfo (_ExecutionList& current_program) {
417 static const _String kPairwiseCountAmbiguitiesResolve ("RESOLVE_AMBIGUITIES"),
418 kPairwiseCountAmbiguitiesAverage ("AVERAGE_AMBIGUITIES"),
419 kPairwiseCountAmbiguitiesSkip ("SKIP_AMBIGUITIES"),
420 kCharacters ("CHARACTERS"),
421 kConsensus ("CONSENSUS"),
422 kParameters ("PARAMETERS"),
423 kPattern ("PATTERN"),
424 kSite ("SITE");
425
426
427 _Variable * receptacle = nil;
428 current_program.advance();
429
430 try {
431
432 receptacle = _ValidateStorageVariable (current_program);
433 const _String source_name = AppendContainerName (*GetIthParameter(1), current_program.nameSpacePrefix);
434
435 long object_type = HY_BL_DATASET|HY_BL_DATASET_FILTER;
436 BaseRefConst source_object = _GetHBLObjectByType(source_name, object_type, nil, ¤t_program);
437
438 _DataSetFilter const * filter_source = object_type == HY_BL_DATASET_FILTER ? (_DataSetFilter const *)source_object : nil;
439 _DataSet const * dataset_source = filter_source ? nil : (_DataSet const *)source_object;
440
441 switch (parameters.lLength) {
442 case 2UL: { // get site->pattern map
443 if (filter_source) {
444 receptacle->SetValue (new _Matrix (filter_source->duplicateMap),false,true, NULL);
445 } else {
446 receptacle->SetValue (new _Matrix (dataset_source->DuplicateMap()),false,true, NULL);
447 }
448 }
449 break;
450
451 case 3UL: { // data parameters, or sequence string
452
453 _String argument;
454 try {
455 argument = _ProcessALiteralArgument (*GetIthParameter(2),current_program);
456 } catch (const _String& err) {
457 //printf ("%s\n", err.get_str());
458 // not a string
459 }
460 if (argument.nonempty()) {
461 if (argument == kCharacters) {
462 _List characters;
463 if (filter_source) {
464 unsigned long character_count = filter_source->GetDimension(true),
465 fd = filter_source->GetUnitLength();
466
467 for (unsigned long idx = 0UL; idx < character_count; idx++) {
468 characters < new _String (filter_source->ConvertCodeToLetters (filter_source->CorrectCode(idx), fd));
469 }
470 } else {
471 _String alphabet_string = dataset_source->GetTT () ? dataset_source->GetTT ()->GetAlphabetString() : kEmptyString;
472 for (unsigned long idx = 0UL; idx < alphabet_string.length(); idx++) {
473 characters < new _String (alphabet_string (idx));
474 }
475 }
476 receptacle->SetValue (new _Matrix (characters), false, true, NULL);
477 } else if (argument == kParameters) {
478 if (filter_source) {
479 _AssociativeList * parameterInfo = new _AssociativeList;
480
481 (*parameterInfo) < (_associative_list_key_value){"ATOM_SIZE", new _Constant (filter_source->GetUnitLength())}
482 < (_associative_list_key_value){"EXCLUSIONS", new _FString (filter_source->GetExclusions())}
483 < (_associative_list_key_value){"SITES_STRING", new _FString ((_String*)filter_source->theOriginalOrder.ListToPartitionString())}
484 < (_associative_list_key_value){"SEQUENCES_STRING", new _FString ((_String*)filter_source->theNodeMap.ListToPartitionString())};
485
486 receptacle->SetValue (parameterInfo,false,true, NULL);
487
488 } else {
489 throw (argument.Enquote('\'') & " is only available for DataSetFilter objects");
490 }
491 } else if (argument == kConsensus) { // argument == _String("PARAMETERS")
492 if (filter_source) {
493 receptacle->SetValue (new _FString (new _String(filter_source->GenerateConsensusString())), false,true, NULL);
494 } else {
495 _DataSetFilter temp;
496 _SimpleList l1, l2;
497 temp.SetFilter (dataset_source, 1, l1, l2, false);
498 receptacle->SetValue (new _FString (new _String(temp.GenerateConsensusString())), false,true, NULL);
499 }
500 }
501 } else {
502 long seqID = _ProcessNumericArgumentWithExceptions (*GetIthParameter(2),current_program.nameSpacePrefix);
503
504 if (filter_source) {
505 if (seqID>=0 && seqID < filter_source->NumberSpecies()) {
506 receptacle->SetValue (new _FString (filter_source->GetSequenceCharacters(seqID)),false,true, NULL);
507 } else if (seqID >= -4 && seqID <= -1) {
508 _SimpleList indices, map, counts;
509 _hy_dataset_filter_unique_match match_mode = kUniqueMatchExact;
510 switch (seqID) {
511 case -2:
512 match_mode = kUniqueMatchExactOrGap;
513 break;
514 case -3:
515 match_mode = kUniqueMatchSuperset;
516 break;
517 case -4:
518 match_mode = kUniqueMatchPartialMatch;
519 break;
520 }
521 long uniqueSequences = filter_source->FindUniqueSequences(indices, map, counts, match_mode);
522 _AssociativeList * parameterInfo = new _AssociativeList;
523 parameterInfo->MStore ("UNIQUE_SEQUENCES", new _Constant (uniqueSequences), false);
524 parameterInfo->MStore ("UNIQUE_INDICES", new _Matrix (indices), false);
525 parameterInfo->MStore ("SEQUENCE_MAP", new _Matrix (map), false);
526 parameterInfo->MStore ("UNIQUE_COUNTS", new _Matrix (counts), false);
527 receptacle->SetValue (parameterInfo,false,true, NULL);
528 }
529 } else { // filter_source
530 if (seqID>=0 && seqID < dataset_source->NoOfSpecies()) {
531 receptacle->SetValue (new _FString (dataset_source->GetSequenceCharacters(seqID)),false,true, NULL);
532 }
533 }
534 } // else numeric cases
535 }
536 break;
537
538 case 4UL : {
539 if (filter_source) {
540 long seq = _ProcessNumericArgumentWithExceptions (*GetIthParameter(2),current_program.nameSpacePrefix),
541 site = _ProcessNumericArgumentWithExceptions (*GetIthParameter(3),current_program.nameSpacePrefix);
542
543 if (site >=0 && site<filter_source->GetPatternCount()) {
544 if ( seq>=0 && seq<filter_source->NumberSpecies()) {
545 bool count_gaps = hy_env::EnvVariableTrue(hy_env::harvest_frequencies_gap_options);
546 _Matrix * res = new _Matrix (filter_source->GetDimension (true), 1, false, true);
547
548 bool only_the_index = hy_env::EnvVariableTrue(hy_env::get_data_info_returns_only_the_index);
549
550 _String character (filter_source->RetrieveState(site, seq));
551
552
553 long theValue = filter_source->Translate2Frequencies (character, res->theData, count_gaps);
554
555 if (only_the_index) {
556 receptacle->SetValue (new _Constant (theValue),false,true, NULL);
557 DeleteObject (res);
558 } else {
559 receptacle->SetValue (res,false,true, NULL);
560 }
561 } else {
562 bool count_gaps = hy_env::EnvVariableTrue(hy_env::harvest_frequencies_gap_options);
563 long filter_dimension = filter_source->GetDimension (true);
564
565 _Matrix * accumulator = new _Matrix (filter_dimension, 1, false, true),
566 * storage = new _Matrix (filter_dimension, 1, false, true);
567
568 _String *buffer = filter_source->MakeSiteBuffer();
569
570 for (long species_index = filter_source->NumberSpecies()-1; species_index >= 0; species_index --) {
571 filter_source->RetrieveState(site,species_index,*buffer, false);
572 filter_source->Translate2Frequencies (*buffer, storage->theData, count_gaps);
573 *accumulator += *storage;
574 }
575 receptacle -> SetValue (accumulator, false,true, NULL);
576 BatchDelete(storage, buffer);
577
578 }
579 } else {
580 throw (_String ("Site index ") & site & " is invalid: must be in range " & "[0, " & (long)filter_source->GetPatternCount() & "]");
581 }
582 } else {
583 throw _String("This set of arguments is only supported for DataSetFilter objects");
584 }
585 }
586 break;
587
588 case 5UL: {
589 if (filter_source) {
590 long seq1 = _ProcessNumericArgumentWithExceptions (*GetIthParameter(2),current_program.nameSpacePrefix),
591 seq2 = _ProcessNumericArgumentWithExceptions (*GetIthParameter(3),current_program.nameSpacePrefix);
592
593 if ( seq1>=0 && seq2 >=0 && seq1< filter_source->NumberSpecies() && seq2 <filter_source->NumberSpecies()) {
594 _String const * res_flag = GetIthParameter(4);
595 _Matrix * res;
596
597 if (kPairwiseCountAmbiguitiesAverage == *res_flag) {
598 res = filter_source->ComputePairwiseDifferences (seq1,seq2,kAmbiguityHandlingAverageFrequencyAware);
599 } else if (kPairwiseCountAmbiguitiesResolve == *res_flag) {
600 res = filter_source->ComputePairwiseDifferences (seq1,seq2,kAmbiguityHandlingResolve);
601 } else if (kPairwiseCountAmbiguitiesSkip == *res_flag) {
602 res = filter_source->ComputePairwiseDifferences (seq1,seq2,kAmbiguityHandlingSkip);
603 } else {
604 res = filter_source->ComputePairwiseDifferences (seq1,seq2,kAmbiguityHandlingResolveFrequencyAware);
605 }
606
607 receptacle->SetValue (res,false,true, NULL);
608 } else {
609 throw (_String (seq1).Enquote() & "," & _String (seq2).Enquote() & " is an invalid sequence pair specification.");
610 }
611 } else {
612 throw _String("This set of options is not supported for DataSet arguments");
613 }
614 }
615 break;
616 // switch
617 }
618 } catch (const _String& error) {
619 return _DefaultExceptionHandler (receptacle, error, current_program);
620 }
621
622 return true;
623 }
624
625 //____________________________________________________________________________________
626
HandleGetInformation(_ExecutionList & current_program)627 bool _ElementaryCommand::HandleGetInformation (_ExecutionList& current_program) {
628 _Variable * receptacle = nil;
629 current_program.advance();
630 try {
631
632 _Matrix* result = nil;
633 receptacle = _ValidateStorageVariable (current_program);
634 const _String source_name = AppendContainerName (*GetIthParameter(1), current_program.nameSpacePrefix);
635
636 long object_type = HY_BL_LIKELIHOOD_FUNCTION | HY_BL_DATASET_FILTER | HY_BL_MODEL ,
637 object_index;
638 BaseRefConst source_object = _HYRetrieveBLObjectByName (source_name, object_type,&object_index,false);
639
640
641 if (source_object) {
642 switch (object_type) {
643 case HY_BL_LIKELIHOOD_FUNCTION: {
644 // list of ctagory variables
645 _LikelihoodFunction const * lf = (_LikelihoodFunction const *)source_object;
646
647 _List catVars;
648 for (unsigned long k=0UL; k<lf->GetCategoryVars().countitems(); k++) {
649 catVars << lf->GetIthCategoryVar (k)->GetName();
650 }
651 result = new _Matrix (catVars);
652 }
653 break;
654 case HY_BL_DATASET_FILTER : {
655 result = ((_DataSetFilter const *) source_object)->GetFilterCharacters();
656 }
657 break;
658 case HY_BL_MODEL: {
659
660 _List modelPNames;
661
662 PopulateAndSort ([&] (_AVLList & parameter_list) -> void {
663 if (IsModelOfExplicitForm (object_index)) {
664 ((_Formula const*)source_object)->ScanFForVariables(parameter_list,false);
665 } else {
666 ((_Variable const*)source_object)->ScanForVariables(parameter_list,false);
667 }
668 }).Each ([&] (long value, unsigned long) -> void {
669 modelPNames << LocateVar(value)->GetName();
670 });
671
672 result = new _Matrix (modelPNames);
673 }
674 break;
675 }
676 } else {
677 _Variable* source_object = FetchVar(LocateVarByName (source_name));
678
679 if (source_object && source_object->ObjectClass()==STRING) {
680 source_object = FetchVar (LocateVarByName (_String((_String*)source_object->Compute()->toStr())));
681 }
682 if (source_object) {
683 if (source_object->IsCategory()) {
684 _CategoryVariable * thisCV = (_CategoryVariable*)source_object;
685 thisCV->Refresh();
686
687 _Matrix *values = thisCV->GetValues(),
688 *weights = thisCV->GetWeights(!thisCV->IsUncorrelated());
689
690 long size = values->GetHDim()*values->GetVDim();
691 result = new _Matrix (2,size,false,true);
692
693 for (unsigned long k = 0UL; k< size ; k++) {
694 result->theData[k] = values->theData[k];
695 result->theData[size+k] = weights->theData[k];
696 }
697 } else {
698 if (source_object->ObjectClass()==TREE_NODE) {
699 _CalcNode* theNode = (_CalcNode*)source_object;
700 if (theNode->GetModelIndex() != HY_NO_MODEL) {
701 result = new _Matrix;
702 theNode->RecomputeMatrix (0,1,result);
703 }
704 } else {
705 if (source_object->ObjectClass() == TOPOLOGY || source_object->ObjectClass() == TREE) {
706
707 _List* map = ((_TreeTopology*)source_object)->MapNodesToModels ();
708 _AssociativeList* return_this = new _AssociativeList();
709
710 for (unsigned long i = 0; i < map->lLength; i++) {
711 _List * nodeInfo = (_List*) map->GetItem(i);
712 return_this->MStore(*(_String*)nodeInfo->GetItem(0), *(_String*)nodeInfo->GetItem (1));
713 }
714 result = (_Matrix*) return_this;
715 DeleteObject (map);
716 }
717 }
718
719 if ((!result)&& source_object->ObjectClass()==NUMBER) {
720 result = new _Matrix (1,3,false,true);
721 result->theData[0]=source_object->Compute()->Value();
722 result->theData[1]=source_object->GetLowerBound();
723 result->theData[2]=source_object->GetUpperBound();
724 }
725 }
726 } else {
727 // TODO : SLKP 20170702, check that this still works, eh?
728 _String reg_exp = GetStringFromFormula (&source_name,current_program.nameSpacePrefix);
729 if (reg_exp != *source_name) {
730 int errNo = 0;
731 regex_t* regex = PrepRegExp (reg_exp, errNo, true);
732 if (regex) {
733 _List matches;
734
735 for (AVLListXIteratorKeyValue variable_record : AVLListXIterator (&variableNames)) {
736 _String* vName = (_String*)variableNames.Retrieve (variable_record.get_index());
737 if (vName->RegExpMatch (regex).lLength) {
738 matches << vName;
739 }
740 }
741 if (matches.lLength) {
742 result = new _Matrix (matches);
743 }
744 _String::FlushRegExp (regex);
745 } else {
746 HandleApplicationError (_String::GetRegExpError (errNo));
747 }
748 }
749 }
750 }
751 if (!result) {
752 result = new _Matrix (0,0,false,false);
753 }
754 receptacle->SetValue(result, false,true, NULL);
755 } catch (const _String& error) {
756 return _DefaultExceptionHandler (receptacle, error, current_program);
757 }
758 return true;
759 }
760
761 //____________________________________________________________________________________
762
HandleConstructCategoryMatrix(_ExecutionList & current_program)763 bool _ElementaryCommand::HandleConstructCategoryMatrix (_ExecutionList& current_program) {
764
765 static _Trie kRunOptions (
766 _String ("COMPLETE"),_hyphyLFConstructCategoryMatrixConditionals,
767 _String ("WEIGHTS"),_hyphyLFConstructCategoryMatrixWeights,
768 _String ("SITE_LOG_LIKELIHOODS"), _hyphyLFConstructCategoryMatrixSiteProbabilities,
769 _String ("CLASSES"), _hyphyLFConstructCategoryMatrixClasses,
770 _String ("SHORT"), _hyphyLFConstructCategoryMatrixClasses
771 );
772
773
774 _Variable * receptacle = nil;
775 current_program.advance();
776
777 try {
778
779 _Matrix* result = nil;
780 receptacle = _ValidateStorageVariable (current_program);
781 const _String source_name = AppendContainerName (*GetIthParameter(1), current_program.nameSpacePrefix);
782 long object_type = HY_BL_LIKELIHOOD_FUNCTION | HY_BL_TREE,
783 object_index;
784 BaseRefConst source_object = _GetHBLObjectByType (source_name, object_type, &object_index, ¤t_program);
785
786 switch (object_type) {
787 case HY_BL_LIKELIHOOD_FUNCTION: {
788 _Matrix * partition_list = nil;
789 if (parameters.countitems () > 3) { // have a restricting partition
790 partition_list = (_Matrix*)_ProcessAnArgumentByType(*GetIthParameter(3), MATRIX, current_program, nil);
791 }
792 _SimpleList included_partitions;
793 _LikelihoodFunction * like_func = (_LikelihoodFunction * )source_object;
794
795 like_func->ProcessPartitionList(included_partitions, partition_list);
796 DeleteObject (partition_list);
797
798 long run_mode = _hyphyLFConstructCategoryMatrixConditionals;
799
800 if (parameters.countitems () > 2) {
801 long run_mode_long = kRunOptions.GetValueFromString(*GetIthParameter(2));
802 if (run_mode_long != kNotFound) {
803 run_mode = run_mode_long;
804 }
805 }
806
807 receptacle->SetValue(like_func->ConstructCategoryMatrix(included_partitions, run_mode ,true, receptacle->GetName()), false,true, NULL);
808 }
809 break;
810
811 case HY_BL_TREE: {
812 _TheTree *source_tree = (_TheTree*)source_object;
813
814
815 long which_partition = 0L,
816 linked_likelihood_id = source_tree->IsLinkedToALF (which_partition);
817
818 if (linked_likelihood_id >= 0) {
819 _LikelihoodFunction * linked_lf = (_LikelihoodFunction*) likeFuncList (linked_likelihood_id);
820 const _DataSetFilter * filter = linked_lf->GetIthFilter (which_partition);
821 linked_lf->PrepareToCompute();
822 linked_lf->Compute ();
823 long patterns = filter->GetPatternCount();
824
825 _Matrix *conditional_matrix = new _Matrix (2*patterns*(source_tree->GetLeafCount()
826 + source_tree->GetINodeCount()) * source_tree->categoryCount,
827 source_tree->GetCodeBase(),
828 false, true);
829
830 _List leaf_names,
831 internal_names;
832
833 _TreeIterator ti (source_tree, _HY_TREE_TRAVERSAL_POSTORDER);
834
835 while (_CalcNode * iterator = ti.Next()) {
836 if (ti.IsAtLeaf()) {
837 leaf_names.AppendNewInstance (new _String(iterator->ContextFreeName()));
838 } else {
839 internal_names.AppendNewInstance (new _String(iterator->ContextFreeName()));
840 }
841 }
842
843 leaf_names << internal_names;
844
845 for (unsigned long site = 0UL; site < patterns; site ++) {
846 source_tree->RecoverNodeSupportStates (filter,site,*conditional_matrix);
847 }
848
849 linked_lf->DoneComputing ();
850 receptacle->SetValue (
851 &(*(new _AssociativeList)
852 < _associative_list_key_value ({"Nodes", new _Matrix (leaf_names)})
853 < _associative_list_key_value ({"Values", conditional_matrix})),
854 false,true, NULL
855 );
856 }
857 }
858 break;
859 }
860 } catch (const _String& error) {
861 return _DefaultExceptionHandler (receptacle, error, current_program);
862 }
863
864 return true;
865
866 }
867 //____________________________________________________________________________________
868
HandleAlignSequences(_ExecutionList & current_program)869 bool _ElementaryCommand::HandleAlignSequences(_ExecutionList& current_program) {
870 _Variable * receptacle = nil;
871
872 static const _String kCharacterMap ("SEQ_ALIGN_CHARACTER_MAP"),
873 kScoreMatrix ("SEQ_ALIGN_SCORE_MATRIX"),
874 kGapChar ("SEQ_ALIGN_GAP_CHARACTER"),
875 kGapOpen ("SEQ_ALIGN_GAP_OPEN"),
876 kGapExtend ("SEQ_ALIGN_GAP_EXTEND"),
877 kGapOpen2 ("SEQ_ALIGN_GAP_OPEN2"),
878 kGapExtend2 ("SEQ_ALIGN_GAP_EXTEND2"),
879 kFrameShift ("SEQ_ALIGN_FRAMESHIFT"),
880 kGapLocal ("SEQ_ALIGN_NO_TP"),
881 kAffineGaps ("SEQ_ALIGN_AFFINE"),
882 kCodonAlign ("SEQ_ALIGN_CODON_ALIGN"),
883 kLinearSpace ("SEQ_ALIGN_LINEAR_SPACE"),
884 kScoreMatrixCodon3x1 ("SEQ_ALIGN_PARTIAL_3x1_SCORES"),
885 kScoreMatrixCodon3x2 ("SEQ_ALIGN_PARTIAL_3x2_SCORES"),
886 kScoreMatrixCodon3x4 ("SEQ_ALIGN_PARTIAL_3x4_SCORES"),
887 kScoreMatrixCodon3x5 ("SEQ_ALIGN_PARTIAL_3x5_SCORES"),
888 kLocalAlignment ("SEQ_ALIGN_LOCAL_ALIGNMENT");
889
890
891 current_program.advance();
892
893 _List dynamic_variable_cleanup;
894
895 try {
896 _Matrix * result = nil;
897 receptacle = _ValidateStorageVariable (current_program);
898 _Matrix * input_seqs = (_Matrix *)_ProcessAnArgumentByType(*GetIthParameter(1), MATRIX, current_program, &dynamic_variable_cleanup);
899
900 unsigned long input_seq_count = input_seqs->GetSize ();
901
902 auto string_validator = [] (long row, long col, _Formula* cell) -> bool {
903 if (cell) {
904 if (cell->ObjectClass() != STRING) {
905 throw (_String(" Matrix entry (") & row & "," & col & ") did not evaluate to a string");
906 }
907 return true;
908 }
909 throw (_String("Empty matrix entry (") & row & "," & col & ")");
910 return false;
911 };
912
913 if (! (input_seqs->IsAStringMatrix() && (input_seqs->is_row() || input_seqs->is_column()) && input_seq_count >= 2 && input_seqs->ValidateFormulaEntries (string_validator))) {
914 throw (GetIthParameter(1)->Enquote() & " did not evaluate to a dense string vector with ≥2 entries");
915 }
916
917
918 _AssociativeList* alignment_options = (_AssociativeList *)_ProcessAnArgumentByType(*GetIthParameter(2), ASSOCIATIVE_LIST, current_program, &dynamic_variable_cleanup);
919 _FString * char_vector = (_FString*) _EnsurePresenceOfKey (alignment_options, kCharacterMap, STRING);
920
921 unsigned long char_count = 0UL;
922 long character_map_to_integers [256];
923 InitializeArray(character_map_to_integers, 256, -1L);
924
925 for (unsigned long cc = 0UL; cc < char_vector->get_str().length(); cc++) {
926 unsigned char this_char = char_vector->get_str().get_uchar(cc);
927 if (character_map_to_integers [this_char]>=0) {
928 throw (_String ("Duplicate character ") & _String ((char)this_char).Enquote('\'') & " in " & kCharacterMap);
929 } else {
930 character_map_to_integers [this_char] = cc;
931 char_count ++;
932 }
933 }
934 if (char_count == 0) {
935 throw _String("Null alphabet supplied");
936 }
937
938 bool do_local = _NumericValueFromKey (alignment_options,kGapLocal,0.0) > 0.5,
939 do_affine = _NumericValueFromKey (alignment_options,kAffineGaps,0.0) > 0.5,
940 do_linear = _NumericValueFromKey (alignment_options,kLinearSpace,1.0) > 0.5,
941 do_codon = _NumericValueFromKey (alignment_options,kCodonAlign,0.0) > 0.5,
942 do_full_local = do_codon && _NumericValueFromKey (alignment_options,kLocalAlignment,0.0) > 0.5;
943
944
945 long codon_count = char_count * char_count * char_count,
946 expected_dimension = (do_codon ? codon_count : char_count) + 1UL;
947
948 _Matrix * score_matrix = (_Matrix*)_EnsurePresenceOfKey(alignment_options, kScoreMatrix, MATRIX);
949
950 if (!score_matrix->check_dimension (expected_dimension, expected_dimension)) {
951 throw (_String ("The dimension of the scoring matrix ") & kScoreMatrix.Enquote('(',')') & " was not the expected dimension: " & expected_dimension & 'x' & expected_dimension);
952 }
953
954 score_matrix = (_Matrix*)score_matrix->ComputeNumeric();
955 score_matrix->CheckIfSparseEnough(true); // force to not be sparse
956
957 _Matrix *codon3x5 = nil,
958 *codon3x4 = nil,
959 *codon3x2 = nil,
960 *codon3x1 = nil;
961
962 if (do_codon) {
963
964 unsigned long expected_columns [4] = {codon_count * 10UL,codon_count * 4UL,char_count * char_count * 3UL,char_count * 3UL};
965 _String const * keys [4] = {&kScoreMatrixCodon3x5, &kScoreMatrixCodon3x4, &kScoreMatrixCodon3x2, &kScoreMatrixCodon3x1};
966 _Matrix ** targets [4] = {&codon3x5, &codon3x4, &codon3x2, &codon3x1};
967
968 for (int i = 0; i < 4; i++) {
969 (*targets[i]) = (_Matrix*)_EnsurePresenceOfKey(alignment_options, *keys[i], MATRIX);
970 if (!(*targets[i])->check_dimension (expected_dimension, expected_columns[i])) {
971 throw (_String ("The dimension of the scoring matrix ") & keys[i]->Enquote('(',')') & " was not the expected dimension: " & expected_dimension & 'x' & (long)expected_columns[i]);
972 }
973 (*targets[i]) = (_Matrix*)(*targets[i])->ComputeNumeric();
974 (*targets[i]) -> CheckIfSparseEnough(true);
975 }
976
977 for (unsigned long i = 0UL; i < 256UL; i++) {
978 // this maps all undefined characters to '?' essentially
979 if (character_map_to_integers[i] < 0) {
980 character_map_to_integers[i] = -codon_count - 1;
981 }
982 }
983 }
984
985 char gap_character = '-';
986 if (_FString *gap_c = (_FString*)alignment_options->GetByKey (kGapChar, STRING)) {
987 if (gap_c->get_str().length () != 1UL) {
988 throw (_String ("Invalid gap character specification ") & gap_c->get_str());
989 }
990 gap_character = gap_c->get_str().char_at(0UL);
991 }
992
993 hyFloat gap_open = _NumericValueFromKey (alignment_options, kGapOpen,15.),
994 gap_open2 = _NumericValueFromKey (alignment_options, kGapOpen2,gap_open),
995 gap_extend = _NumericValueFromKey (alignment_options, kGapExtend,1.),
996 gap_extend2 = _NumericValueFromKey (alignment_options, kGapExtend2,gap_extend),
997 gap_frameshift = _NumericValueFromKey (alignment_options, kFrameShift,50.);
998
999 _StringBuffer settings_report (256L);
1000
1001 settings_report << "\n\tGap character : " << gap_character
1002 << "\n\tGap open cost [reference] : " << gap_open
1003 << "\n\tGap open cost [query] : " << gap_open2
1004 << "\n\tGap extend cost [reference] : " << gap_extend
1005 << "\n\tGap extend cost [query] : " << gap_extend2
1006 << "\n\tCodon frameshift cost : " << gap_frameshift
1007 << "\n\tIgnore terminal gaps : " << (do_local?"Yes":"No")
1008 << "\n\tPerform local alignment : " << (do_full_local?"Yes":"No");
1009
1010 if (do_codon) {
1011 settings_report << "\n\tUse codon alignment with frameshift routines";
1012 do_linear = false;
1013 }
1014
1015 _AssociativeList *aligned_strings = new _AssociativeList;
1016 _String const * reference_sequence = & ((_FString*)input_seqs->GetFormula(0,0)->Compute())->get_str();
1017
1018
1019 for (unsigned long index2 = 1UL; index2 < input_seq_count; index2++) {
1020 _String const * sequence2 = & ((_FString*)input_seqs->GetFormula(0,index2)->Compute())->get_str();
1021 _AssociativeList * pairwise_alignment = new _AssociativeList;
1022 hyFloat score = 0.0;
1023 if (do_linear) {
1024 unsigned long size_allocation = sequence2->length()+1UL;
1025
1026 _Matrix *buffers[6] = {nil};
1027
1028 ArrayForEach(buffers, 6, [=] (_Matrix* m, unsigned long) -> _Matrix* {
1029 return new _Matrix (size_allocation,1,false,true);
1030 });
1031
1032 char *alignment_route = new char[2*(size_allocation)] {0};
1033
1034 _SimpleList ops (reference_sequence->length()+2UL,-2,0);
1035 ops.list_data[reference_sequence->length()+1] = sequence2->length();
1036 ops.list_data[0] = -1;
1037
1038 score = LinearSpaceAlign(reference_sequence,sequence2,character_map_to_integers,score_matrix,
1039 gap_open,gap_extend,gap_open2,gap_extend2,
1040 do_local,do_affine,ops,score,0,
1041 reference_sequence->length(),0,sequence2->length(),buffers,0,alignment_route);
1042
1043 delete[] alignment_route;
1044
1045 _StringBuffer *result1 = new _StringBuffer (reference_sequence->length() + 1UL),
1046 *result2 = new _StringBuffer (size_allocation);
1047
1048 long last_column = ops.list_data[ops.lLength-1];
1049
1050 for (long position = (long)reference_sequence->length() - 1L; position>=0; position--) {
1051 long current_column = ops.list_data[position+1];
1052
1053 if (current_column<0) {
1054 if (current_column == -2 /*|| (current_column == -3 && last_column == string2->sLength)*/) {
1055 current_column = last_column;
1056 } else if (current_column == -3) {
1057 // find the next matched char or a -1
1058 long p = position, s2p;
1059 while ((ops.list_data[p+1]) < -1) {
1060 p--;
1061 }
1062
1063 s2p = ops.list_data[p+1];
1064
1065 for (long j = last_column-1; j>s2p;) {
1066 (*result1) << gap_character;
1067 (*result2) << sequence2->char_at(j--);
1068 }
1069
1070 last_column = s2p+1;
1071
1072 for (; position>p; position--) {
1073 (*result2) << gap_character;
1074 (*result1) << reference_sequence->char_at(position);
1075 }
1076 position ++;
1077 continue;
1078 } else {
1079 for (last_column--; last_column >=0L; last_column--) {
1080 (*result1) << gap_character;
1081 (*result2) << sequence2->char_at (last_column);
1082 }
1083 while (position>=0) {
1084 (*result1) << reference_sequence->char_at (position--);
1085 (*result2) << gap_character;
1086 }
1087 break;
1088 }
1089 }
1090
1091 if (current_column == last_column) { // insert in sequence 2
1092 (*result1) << reference_sequence->char_at (position);
1093 (*result2) << gap_character;
1094 } else {
1095 last_column--;
1096
1097 for (; last_column > current_column; last_column--) { // insert in column 1
1098 (*result2) << sequence2->char_at (last_column);
1099 (*result1) << gap_character;
1100 }
1101 (*result1) << reference_sequence->char_at (position);
1102 (*result2) << sequence2->char_at (current_column);
1103 }
1104 //printf ("%s\n%s\n", result1->sData, result2->sData);
1105 }
1106
1107 for (last_column--; last_column >=0; last_column--) {
1108 (*result1) << gap_character;
1109 (*result2) << sequence2->char_at(last_column);
1110 }
1111
1112 result1->Flip ();
1113 result2->Flip ();
1114 pairwise_alignment->MStore ("1", new _FString(result1), false);
1115 pairwise_alignment->MStore ("2", new _FString(result2), false);
1116 pairwise_alignment->MStore ("0", new _Constant (score), false);
1117 aligned_strings->MStore (_String((long)index2-1L), pairwise_alignment, false);
1118 } else { // not linear
1119 char * str1r = nil,
1120 * str2r = nil;
1121
1122
1123 score = AlignStrings (reference_sequence->get_str() ? reference_sequence->get_str() : "",
1124 sequence2->get_str() ? sequence2->get_str() : "",
1125 str1r,
1126 str2r,
1127 character_map_to_integers,
1128 score_matrix->fastIndex(),
1129 score_matrix->GetVDim(),
1130 gap_character,
1131 gap_open,
1132 gap_extend,
1133 gap_open2,
1134 gap_extend2,
1135 gap_frameshift,
1136 do_local,
1137 do_affine,
1138 do_codon,
1139 char_count,
1140 do_codon ? codon3x5->fastIndex() : nil,
1141 do_codon ? codon3x4->fastIndex() : nil,
1142 do_codon ? codon3x2->fastIndex() : nil,
1143 do_codon ? codon3x1->fastIndex() : nil,
1144 do_full_local);
1145
1146 if ( str1r && str2r ) {
1147 pairwise_alignment->MStore ("1", new _FString (new _String( str1r )), false);
1148 pairwise_alignment->MStore ("2", new _FString (new _String( str2r )), false);
1149 delete [] str1r;
1150 delete [] str2r;
1151 } else {
1152 throw _String( "Internal Error in AlignStrings" );
1153 }
1154 pairwise_alignment->MStore ("0", new _Constant (score), false);
1155 aligned_strings->MStore (_String((long)index2-1L), pairwise_alignment, false);
1156 }
1157 }
1158 receptacle->SetValue(aligned_strings, false,true, NULL);
1159
1160 } catch (const _String& error) {
1161 return _DefaultExceptionHandler (receptacle, error, current_program);
1162 }
1163
1164
1165 return true;
1166 }
1167 //____________________________________________________________________________________
1168
HandleHarvestFrequencies(_ExecutionList & current_program)1169 bool _ElementaryCommand::HandleHarvestFrequencies (_ExecutionList& current_program) {
1170
1171 _Variable * receptacle = nil;
1172 current_program.advance();
1173
1174 try {
1175 _Matrix * result = nil;
1176
1177 receptacle = _ValidateStorageVariable (current_program);
1178
1179 long object_type = HY_BL_DATASET|HY_BL_DATASET_FILTER;
1180 BaseRefConst source_object = _GetHBLObjectByType(*GetIthParameter(1), object_type, nil, ¤t_program);
1181
1182 long unit = _ProcessNumericArgumentWithExceptions(*GetIthParameter(2),current_program.nameSpacePrefix),
1183 atom = _ProcessNumericArgumentWithExceptions(*GetIthParameter(3),current_program.nameSpacePrefix);
1184
1185 bool position_specific = _ProcessNumericArgumentWithExceptions(*GetIthParameter(4),current_program.nameSpacePrefix) > 0.5,
1186 include_gaps = hy_env::EnvVariableTrue(hy_env::harvest_frequencies_gap_options);
1187
1188 switch (object_type) {
1189 case HY_BL_DATASET: {
1190 _String vertical_partition (parameters.countitems () > 5 ? *GetIthParameter(5) : kEmptyString),
1191 horizontal_partition (parameters.countitems () > 6 ? *GetIthParameter(6) : kEmptyString);
1192
1193 _DataSet const * dataset = (_DataSet const*)source_object;
1194 _SimpleList processed_sequence_partition, processed_site_partition;
1195 dataset->ProcessPartition (horizontal_partition,processed_sequence_partition,false, 1);
1196 dataset->ProcessPartition (vertical_partition,processed_site_partition,true, 1);
1197
1198 receptacle->SetValue (dataset->HarvestFrequencies(unit,atom,position_specific,processed_sequence_partition, processed_site_partition,include_gaps), false,true, NULL);
1199 }
1200 break;
1201 case HY_BL_DATASET_FILTER: {
1202 receptacle->SetValue (((_DataSetFilter const*)source_object)->HarvestFrequencies(unit,atom,position_specific,include_gaps), false,true, NULL);
1203 }
1204 break;
1205 }
1206
1207 } catch (const _String& error) {
1208 return _DefaultExceptionHandler (receptacle, error, current_program);
1209 }
1210
1211 return true;
1212 }
1213
1214 //____________________________________________________________________________________
1215
HandleOptimizeCovarianceMatrix(_ExecutionList & current_program,bool do_optimize)1216 bool _ElementaryCommand::HandleOptimizeCovarianceMatrix (_ExecutionList& current_program, bool do_optimize) {
1217 _Variable * receptacle = nil;
1218 current_program.advance();
1219
1220 try {
1221 receptacle = _ValidateStorageVariable (current_program);
1222
1223 long object_type = HY_BL_LIKELIHOOD_FUNCTION|HY_BL_SCFG|HY_BL_BGM;
1224 _String optimize_me = AppendContainerName (*GetIthParameter(1), current_program.nameSpacePrefix);
1225
1226 _LikelihoodFunction* source_object = (_LikelihoodFunction*)_HYRetrieveBLObjectByNameMutable (AppendContainerName (*GetIthParameter(1), current_program.nameSpacePrefix), object_type,nil,do_optimize==false);
1227
1228 if (!source_object) { // Custom function (expression based)
1229 source_object = new _CustomFunction (optimize_me, current_program.nameSpacePrefix);
1230 }
1231
1232 if (do_optimize) {
1233 if (parameters.countitems () > 2) { // have a restricting partition
1234 _List ref;
1235 receptacle -> SetValue(source_object->Optimize((_AssociativeList*)_ProcessAnArgumentByType(*GetIthParameter(2L), ASSOCIATIVE_LIST, current_program, &ref)),false,true, NULL);
1236 } else {
1237 receptacle -> SetValue(source_object->Optimize(),false,true, NULL);
1238 }
1239 } else {
1240 HBLObjectRef covariance_parameters = hy_env::EnvVariableGet(hy_env::covariance_parameter, ASSOCIATIVE_LIST|STRING);
1241 _SimpleList *restrictor = nil;
1242 switch (object_type) {
1243 case HY_BL_LIKELIHOOD_FUNCTION:
1244 case HY_BL_SCFG: {
1245 if (covariance_parameters) { // only consider some variables
1246 _SimpleList variable_ids;
1247 if (covariance_parameters->ObjectClass () == ASSOCIATIVE_LIST) {
1248 // a list of variables stored as keys in an associative array
1249 _List* restricted_variables = ((_AssociativeList*)covariance_parameters)->GetKeys();
1250 for (unsigned long iid = 0; iid < restricted_variables->lLength; iid++) {
1251 variable_ids << LocateVarByName (current_program.AddNameSpaceToID(*(_String*)(*restricted_variables)(iid)));
1252 }
1253 DeleteObject (restricted_variables);
1254 } else { // STRING
1255 variable_ids << LocateVarByName (current_program.AddNameSpaceToID(((_FString*)covariance_parameters)->get_str()));
1256 }
1257 if (!variable_ids.empty()) {
1258 restrictor = new _SimpleList ();
1259
1260 for (unsigned long var_index = 0UL; var_index < variable_ids.lLength; var_index++) {
1261 // TODO SLKP 20170706: this is a very inefficient linear search over and over again (in a large array)
1262 long vID = source_object->GetIndependentVars().Find(variable_ids (var_index));
1263 if (vID >= 0) (*restrictor) << vID;
1264 }
1265
1266 if (restrictor->empty()) {
1267 DeleteAndZeroObject (restrictor);
1268 }
1269 }
1270 }
1271 }
1272 break;
1273 case HY_BL_BGM: {
1274 _Matrix * bgm_cov = (_Matrix*)source_object->CovarianceMatrix(nil);
1275 if (bgm_cov) {
1276 receptacle->SetValue(bgm_cov,false,true, NULL);
1277 return true;
1278 } // TODO SLKP 20170706: handle the case when null is returned (why would that happen?); warn the user.
1279 }
1280 break;
1281 }
1282
1283 receptacle->SetValue(source_object->CovarianceMatrix(restrictor),false,true, NULL);
1284 DeleteObject (restrictor);
1285 }
1286
1287 if (object_type == HY_BL_NOT_DEFINED) {
1288 DeleteObject (source_object); // delete the custom function object
1289 }
1290
1291 } catch (const _String& error) {
1292 return _DefaultExceptionHandler (receptacle, error, current_program);
1293 }
1294
1295 return true;
1296 }
1297
1298 //____________________________________________________________________________________
1299
HandleReplicateConstraint(_ExecutionList & current_program)1300 bool _ElementaryCommand::HandleReplicateConstraint (_ExecutionList& current_program) {
1301 // TODO SLKP 20170706 this needs to be reimplemented; legacy code is ugly and buggy
1302
1303 static const _String kLastSetOfConstraints ("LAST_SET_OF_CONSTRAINTS");
1304
1305 current_program.advance();
1306
1307 _TreeIterator ** traversers = nil;
1308 unsigned long template_parameter_count = parameter_count() - 1L;
1309 deferSetFormula = new _SimpleList;
1310
1311 auto cleanup = [template_parameter_count] (_TreeIterator ** t) -> void {
1312 if (t) {
1313 for (long k = 0L; k < template_parameter_count; k++) {
1314 if (t[k]) {
1315 delete t[k];
1316 }
1317 }
1318 delete [] t;
1319 }
1320 delete deferSetFormula;
1321 };
1322
1323 auto resolve_name = [] (_String const& my_name, BaseRefConst parent_name) -> _String const {
1324 return my_name.Cut (((_String*)parent_name)->length() + 1, kStringEnd);
1325 };
1326
1327
1328 try {
1329
1330 /** check the plug-in first **/
1331
1332 _List templated_operations,
1333 parent_object_names;
1334
1335 traversers = new _TreeIterator* [parameter_count()-1] {NULL};
1336
1337 for (long k = 1L; k < parameter_count(); k++) {
1338 _CalcNode * this_object = (_CalcNode *)_CheckForExistingVariableByType (*GetIthParameter(k), current_program, TREE | TREE_NODE);
1339 if (this_object->ObjectClass () == TREE) {
1340 traversers[k-1] = new _TreeIterator ((_TheTree*)this_object, _HY_TREE_TRAVERSAL_POSTORDER);
1341 parent_object_names << this_object->GetName();
1342 } else {
1343 node<long> * cn = this_object->LocateMeInTree();
1344 if (!cn) {
1345 throw (*GetIthParameter(k) & " is not a part of a tree object");
1346 }
1347 parent_object_names << this_object->ParentTree()->GetName();
1348 traversers[k-1] = new _TreeIterator (this_object, cn, _HY_TREE_TRAVERSAL_POSTORDER);
1349 }
1350 templated_operations < new _List;
1351 }
1352
1353 /*
1354 _CalcNode * check = traversers[1]->Next();
1355 while (check) {
1356 printf ("%s\n", check->GetName()->get_str());
1357 check = traversers[1]->Next();
1358 }
1359
1360 throw (_String("BAIL"));
1361 */
1362
1363 _String const constraint_pattern = _ProcessALiteralArgument (*GetIthParameter(0), current_program);
1364
1365 /** parse the contraint pattern using the standard exprssion parser
1366 because of '?' in template names, and things like __ which need to be deferred,
1367 the standard parser will fail with standard argument calls
1368 */
1369
1370 _String parse_error, patttern_copy (constraint_pattern);
1371 _FormulaParsingContext fpc (&parse_error, current_program.nameSpacePrefix);
1372 fpc.allowTemplate() = '?';
1373
1374 _Formula pattern,
1375 lhs;
1376
1377 long return_value = Parse (&pattern, patttern_copy, fpc, &lhs);
1378
1379 if (return_value != HY_FORMULA_VARIABLE_FORMULA_ASSIGNMENT && return_value != HY_FORMULA_VARIABLE_VALUE_ASSIGNMENT && return_value != HY_FORMULA_EXPRESSION) {
1380 throw (_String ("The template for ReplicateConstraint must be an assignment, with an expression based left-hand side, or simple expression"));
1381 }
1382
1383
1384 /** for 0 -- (K-1), where thisK is the largest "this" argument,
1385 stores the operation objects that reference these variables
1386 **/
1387
1388 _Operation * lhs_op = nil;
1389
1390 if (return_value != HY_FORMULA_EXPRESSION) {
1391 lhs_op = new _Operation (*LocateVar(fpc.assignmentRefID()));
1392 lhs.PushTerm (lhs_op);
1393 lhs_op->RemoveAReference();
1394 }
1395
1396 _Formula * formulas [2] = {&lhs, &pattern};
1397
1398 _List substitution_variables (new _List, new _List);
1399 _List list_of_constraints;
1400
1401 auto add_a_constraint = [return_value, &list_of_constraints] (_String* lhs, _String *rhs) -> void {
1402 _StringBuffer * new_constraint = new _StringBuffer;
1403 switch (return_value) {
1404 case HY_FORMULA_EXPRESSION:
1405 new_constraint->AppendNewInstance(rhs);
1406 break;
1407 case HY_FORMULA_VARIABLE_FORMULA_ASSIGNMENT:
1408 (new_constraint->AppendNewInstance(lhs) << ":=").AppendNewInstance (rhs);
1409 break;
1410 case HY_FORMULA_VARIABLE_VALUE_ASSIGNMENT:
1411 (new_constraint->AppendNewInstance(lhs) << "=").AppendNewInstance (rhs);
1412 break;
1413
1414 }
1415
1416 list_of_constraints < new_constraint;
1417 };
1418
1419 unsigned long reference_argument = 0UL; // either the first argument
1420 // or the argument that appears on the LHS
1421
1422 _SimpleList substitution_variable_by_index;
1423 // this keeps track of which template parameter the i-th thisN.?.. argument comes from
1424
1425 for (_Formula * f : formulas) {
1426 bool is_lhs = f == &lhs;
1427 f->GetList().ForEach([&templated_operations, is_lhs, &reference_argument, &substitution_variables, &substitution_variable_by_index] (BaseRef op, unsigned long index) -> void {
1428 _Variable const * operation_var = ((_Operation*) op)->RetrieveVar();
1429 if (operation_var) {
1430 _SimpleList pattern_match (operation_var->GetName()->RegExpMatch(hy_replicate_constraint_regexp, 0));
1431 if (pattern_match.nonempty()) {
1432 unsigned long var_index = operation_var->GetName()->Cut (pattern_match (2), pattern_match (3)).to_long() - 1UL;
1433 if (var_index >= templated_operations.countitems()) {
1434 throw (operation_var->GetName()->Enquote() & " does not have a matched positional argument");
1435 }
1436 if (is_lhs) {
1437 reference_argument = var_index;
1438 }
1439 _List * term_record = new _List;
1440 term_record->AppendNewInstance (new _String (*operation_var->GetName(), pattern_match (4), pattern_match (5)));
1441 (*term_record) << op;
1442 //printf ("[term_record] %s\n", ((_String*) (term_record->toStr()))->get_str());
1443
1444 ((_List*)templated_operations.GetItem(var_index))->AppendNewInstance(term_record);
1445 *((_List*)(substitution_variables.GetItem(0))) << operation_var->GetName();
1446 substitution_variable_by_index << var_index;
1447 }
1448 }
1449 });
1450 }
1451
1452
1453 templated_operations.ForEach ([this] (BaseRef items, long index) -> void {
1454 if (((_List*)items)->empty()) {
1455 throw (GetIthParameter(index + 1L)->Enquote() & " (argument " & _String (index+1L) & ") did not appear in the contstraint expression");
1456 }
1457 });
1458
1459 /** perform a joint traversal of all the template objects
1460 the first template object will be used to create the list
1461 of local parameters that could be used to populate template contraints
1462 */
1463
1464
1465 auto incompatible_error = [this, reference_argument] (unsigned long i, const _String& name) -> void {
1466 throw (GetIthParameter(i + 1L)->Enquote() & " (argument " & _String ((long)i+1L) &
1467 ") is topologically incompatible with the reference argument " & GetIthParameter(reference_argument + 1UL)->Enquote()) & " at node " & name.Enquote();
1468 };
1469
1470 for (;;) {
1471 _CalcNode * reference_iteratee = traversers[reference_argument]->Next();
1472 node <long> * reference_node = traversers[reference_argument]->GetNode();
1473
1474
1475 if (reference_iteratee) {
1476 //printf ("%s\n", reference_iteratee->GetName()->get_str());
1477 _List iteratees;
1478 iteratees << reference_iteratee;
1479 for (unsigned long i = 0UL; i < template_parameter_count; i++) {
1480 if (reference_argument != i) {
1481 _CalcNode * current_iteratee = traversers[i]->Next();
1482 if (!current_iteratee) {
1483 incompatible_error (i, *reference_iteratee->GetName());
1484 }
1485 //printf ("%s\n", current_iteratee->GetName()->get_str());
1486 node <long> *current_node = traversers[i]->GetNode();
1487 if (current_node->get_num_nodes() != reference_node->get_num_nodes()) {
1488 incompatible_error (i,*reference_iteratee->GetName());
1489 }
1490 iteratees << current_iteratee;
1491 }
1492 }
1493
1494 auto handle_variable = [&parent_object_names,resolve_name] (_Variable* local_k, _List * conditions, _List * parameter_set, _List *matched_subexpressions, long condition_index) -> void {
1495 _String local_name = resolve_name (*local_k->GetName(), parent_object_names.GetItem(condition_index));
1496 _SimpleList matched_conditions;
1497
1498 if (conditions->Every ([&local_name,&matched_conditions] (long condition, unsigned long i) -> bool {
1499 _String const * match_pattern = (_String*)((_List*)condition)->GetItem (0);
1500 //printf ("[name] %s [pattern] %s\n",local_name.get_str(), match_pattern->get_str());
1501 if (i == 0) {
1502 return local_name.EqualWithWildChar (*match_pattern,'?', 0, 0, &matched_conditions);
1503
1504 } else {
1505 return local_name.EqualWithWildChar (*match_pattern,'?');
1506 }
1507 })) {
1508 //printf ("[%ld] %s (%s, %s)\n", condition_index,local_k->GetName()->get_str(), local_name.get_str(), ((_String*)parent_object_names.GetItem(condition_index))->get_str());
1509 *parameter_set << local_k;
1510 _List * subexpressions = new _List;
1511 for (long k = 0L; k < matched_conditions.countitems(); k+=2) {
1512 subexpressions->AppendNewInstance (new _String (local_name, matched_conditions(k), matched_conditions (k+1)));
1513 }
1514 //ObjectToConsole(subexpressions);
1515 *matched_subexpressions < subexpressions;
1516 }
1517 };
1518
1519 if (reference_iteratee->HasLocals() && !traversers[reference_argument]->IsAtRoot()) { // stuff to do
1520 _List parameter_sets,
1521 matched_subexpressions;
1522 for (unsigned long i = 0UL; i < template_parameter_count; i++) {
1523 parameter_sets < new _List;
1524 matched_subexpressions < new _List;
1525 }
1526
1527
1528 /** now see if there any applicable contstraints that can be generated
1529 if there is a LHS expression, then only independent variables will be pulled
1530 */
1531
1532 long independent_count = reference_iteratee->CountIndependents();
1533 _List * conditions = ((_List*)templated_operations.GetItem(reference_argument));
1534 _List * reference_parameters = (_List*)parameter_sets.GetItem(reference_argument);
1535 _List* matched_reference_subexpressions = (_List*)matched_subexpressions.GetItem(reference_argument);
1536 /*
1537 this will store, for each reference argument, the parts of parameter names
1538 that were matched by the '?' in the template expression.
1539 The convention is that the first set takes precedence, so that
1540 this1.?.? := this1.?.synRate will store the matches for the two '?' from the first expression
1541 */
1542
1543
1544 for (long k = 0L; k < independent_count; k++) {
1545 handle_variable ( reference_iteratee->GetIthIndependent(k), conditions, reference_parameters, matched_reference_subexpressions, reference_argument);
1546 }
1547
1548 if (!lhs_op) {
1549 for (long k = 0L; k < reference_iteratee->CountDependents(); k++) {
1550 handle_variable ( reference_iteratee->GetIthDependent(k), conditions, reference_parameters, matched_reference_subexpressions, reference_argument);
1551 }
1552 }
1553
1554 /** now check to see if subsequent template arguments fit reference template parameters
1555 If there are multiple parameter matches for any subsequent template, then
1556 */
1557
1558 if (reference_parameters->nonempty()) {
1559 for (unsigned long i = 0UL; i < template_parameter_count; i++) {
1560 if (i != reference_argument) {
1561 _CalcNode * template_iteratee = (_CalcNode*)iteratees.GetItem(i);
1562 _List * template_variables = (_List*)parameter_sets.GetItem(i);
1563 _List * matched_template_subexpressions = (_List*)matched_subexpressions.GetItem(i);
1564 _List * template_conditions = ((_List*)templated_operations.GetItem(i));
1565 long variable_count = template_iteratee->CountIndependents();
1566
1567
1568 // see which parameters match the template conditions
1569 for (long k = 0L; k < variable_count; k++) {
1570 handle_variable ( template_iteratee->GetIthIndependent(k), template_conditions, template_variables, matched_template_subexpressions, i);
1571 }
1572 variable_count = template_iteratee->CountDependents();
1573 for (long k = 0L; k < variable_count; k++) {
1574 handle_variable ( template_iteratee->GetIthDependent(k), template_conditions, template_variables, matched_template_subexpressions, i);
1575 }
1576 }
1577 }
1578 // now for each reference parameter, see if it is matched by the corresponding template paramters
1579
1580 reference_parameters->ForEach([&substitution_variables, &pattern, &lhs, &matched_subexpressions, reference_argument, ¶meter_sets, add_a_constraint, ¤t_program, &substitution_variable_by_index] (BaseRef ref, unsigned long i) -> void {
1581 _Variable * ref_var = (_Variable *)ref;
1582 _List * reference_subexpressions = (_List*)matched_subexpressions.GetItem(reference_argument,i);
1583 _List local_substitution_variables;
1584
1585 //ObjectToConsole(reference_subexpressions);
1586
1587 if (parameter_sets.Every ([reference_argument, ref_var, &local_substitution_variables, matched_subexpressions, reference_subexpressions] (long template_i, unsigned long i2) -> bool {
1588 if (i2 != reference_argument) {
1589 _List * template_parameter_set = (_List *)template_i;
1590 _List * matched_template_subexpressions = (_List*)matched_subexpressions.GetItem(i2);
1591
1592 return template_parameter_set->Any ([matched_template_subexpressions, reference_subexpressions, &local_substitution_variables] (long argument_ij, unsigned long i3) -> bool {
1593 _List * check_subs = (_List*) matched_template_subexpressions->GetItem (i3);
1594 for (long i = 0; i < Minimum(check_subs->countitems(), reference_subexpressions->countitems()); i++) {
1595 if (*(_String*)check_subs->GetItem(i) != *(_String*)reference_subexpressions->GetItem(i) ) {
1596 // printf ("[fail check] %s %s [%d]\n", ((_String*)check_subs->GetItem(i))->get_str(),((_String*)reference_subexpressions->GetItem(i))->get_str(), check_subs->countitems());
1597 return false;
1598 }
1599 }
1600 local_substitution_variables << (_Variable *)argument_ij;
1601 return true;
1602 });
1603 } else {
1604 local_substitution_variables << ref_var;
1605 return true;
1606 }
1607 })) { // this template argument is OK to go
1608 ((_List*)(substitution_variables.GetItem(1)))->Duplicate(substitution_variables.GetItem(0));
1609 //printf ("%s\n", ((_String*)local_substitution_variables.toStr())->get_str());
1610 substitution_variable_by_index.Each ([&substitution_variables,&local_substitution_variables] (long template_index, unsigned long idx) -> void {
1611 _String * sub_name = ((_Variable*) local_substitution_variables.GetItem(template_index))->GetName();
1612 ((_List*)(substitution_variables.GetItem(1)))->InsertElement (sub_name, idx, false, true);
1613 });
1614
1615
1616 /*local_substitution_variables.ForEach([&substitution_variables] (BaseRef v, unsigned long) -> void {
1617 *((_List*)(substitution_variables.GetItem(1))) << ((_Variable*) v)->GetName();
1618 });*/
1619
1620 add_a_constraint (lhs.IsEmpty() ? nil : ((_String*)lhs.toStr (kFormulaStringConversionNormal, &substitution_variables)),
1621 ((_String*)pattern.toStr (kFormulaStringConversionNormal, &substitution_variables)));
1622 // apply the constraint
1623 }
1624
1625 });
1626 }
1627
1628 }
1629 } else { // possible termination
1630 for (unsigned long i = 0UL; i < template_parameter_count; i++) {
1631 if (reference_argument != i) {
1632 if (traversers[i]->Next()) {
1633 incompatible_error (i, "Root node");
1634 }
1635 }
1636 }
1637 break;
1638 }
1639 }
1640
1641 list_of_constraints.ForEach ([¤t_program] (BaseRefConst constraint_i, unsigned long) -> void {
1642 _String constraint = *(_String*)constraint_i;
1643 _Formula rhs, lhs;
1644 _FormulaParsingContext fpc (nil, current_program.nameSpacePrefix);
1645 long result = Parse (&rhs,constraint,fpc,&lhs);
1646 ExecuteFormula(&rhs,&lhs,result,fpc.assignmentRefID(),current_program.nameSpacePrefix,fpc.assignmentRefType());
1647 });
1648
1649
1650 _StringBuffer * generated_constraints = new _StringBuffer (list_of_constraints.Join (";\n"));
1651 if (generated_constraints->nonempty()) {
1652 *generated_constraints << ";";
1653 }
1654
1655 hy_env::EnvVariableSet(kLastSetOfConstraints, new _FString (generated_constraints), false);
1656
1657 FinishDeferredSF ();
1658
1659 // first find the pattern
1660
1661
1662 } catch (const _String& error) {
1663 cleanup (traversers);
1664 return _DefaultExceptionHandler (nil, error, current_program);
1665 }
1666
1667 cleanup (traversers);
1668 return true;
1669 }
1670
1671 //____________________________________________________________________________________
1672
HandleComputeLFFunction(_ExecutionList & current_program)1673 bool _ElementaryCommand::HandleComputeLFFunction (_ExecutionList& current_program) {
1674
1675 const static _String kLFStartCompute ("LF_START_COMPUTE"),
1676 kLFDoneCompute ("LF_DONE_COMPUTE"),
1677 kLFTrackCache ("LF_TRACK_CACHE"),
1678 kLFAbandonCache ("LF_ABANDON_CACHE");
1679
1680 current_program.advance();
1681 _Variable * receptacle = nil;
1682
1683 try {
1684
1685
1686 _String const op_kind = * GetIthParameter(1UL);
1687
1688 long object_type = HY_BL_LIKELIHOOD_FUNCTION|HY_BL_SCFG|HY_BL_BGM;
1689 _LikelihoodFunction* source_object = (_LikelihoodFunction*)_GetHBLObjectByType(AppendContainerName (*GetIthParameter(0UL), current_program.nameSpacePrefix),object_type, nil,¤t_program);
1690
1691 if (op_kind == kLFStartCompute) {
1692 source_object->PrepareToCompute(true);
1693 } else if (op_kind == kLFDoneCompute) {
1694 source_object->FlushLocalUpdatePolicy();
1695 source_object->DoneComputing (true);
1696 } else {
1697 if (!source_object->HasBeenSetup()) {
1698 throw (_String("Please call LFCompute (") & *GetIthParameter (0UL)& ", " & kLFStartCompute & ") before evaluating the likelihood function");
1699 } else {
1700 if (op_kind == kLFTrackCache) {
1701 source_object->DetermineLocalUpdatePolicy();
1702 } else if (op_kind == kLFAbandonCache) {
1703 source_object->FlushLocalUpdatePolicy();
1704 } else {
1705 receptacle = _ValidateStorageVariable (current_program, 1UL);
1706 receptacle->SetValue (new _Constant (source_object->Compute()), false,true, NULL);
1707 }
1708 }
1709 }
1710
1711 } catch (const _String& error) {
1712 return _DefaultExceptionHandler (receptacle, error, current_program);
1713 }
1714
1715 return true;
1716 }
1717
1718 //____________________________________________________________________________________
1719
HandleUseModel(_ExecutionList & current_program)1720 bool _ElementaryCommand::HandleUseModel (_ExecutionList& current_program) {
1721 const static _String kUseNoModel ("USE_NO_MODEL");
1722 current_program.advance();
1723
1724 try {
1725
1726 _String raw_model_name = *GetIthParameter(0UL),
1727 source_name = AppendContainerName (raw_model_name, current_program.nameSpacePrefix);
1728
1729 long object_type_request = HY_BL_MODEL,
1730 object_type = object_type_request,
1731 model_index = HY_NO_MODEL;
1732
1733 _Variable * source_model = (_Variable*)_HYRetrieveBLObjectByNameMutable (source_name, object_type,&model_index,false);
1734
1735 if (!source_model && raw_model_name != kUseNoModel) {
1736 throw (source_name.Enquote() & " does not refer to a valid defined substitution model and is not " & kUseNoModel);
1737 }
1738
1739 lastMatrixDeclared = model_index;
1740
1741 } catch (const _String& error) {
1742 return _DefaultExceptionHandler (nil, error, current_program);
1743 }
1744
1745 return true;
1746 }
1747
1748
1749 //____________________________________________________________________________________
1750
HandleInitializeIterator(_ExecutionList & current_program)1751 bool _ElementaryCommand::HandleInitializeIterator (_ExecutionList& current_program) {
1752 /*
1753 parameters:
1754 [0] -> the string ID of the object to iterate over
1755 [1] -> the object to iteratre over
1756
1757
1758 simpleParameters:
1759 [0] -> command index for the corresponding advance command
1760 [1] -> type of object; MATRIX, etc
1761 [2] -> if dict, then this stores the pointer to AVLListXLIterator
1762 -> if tree, then this stores the pointer to a tree iterator
1763 */
1764 current_program.advance();
1765 try {
1766 // clean up previous iterator states
1767 if (parameters.countitems() > 1) {
1768 parameters.Delete (1L);
1769 if (simpleParameters.get (1) == ASSOCIATIVE_LIST) {
1770 delete ((AVLListXLIterator*)simpleParameters.get(2));
1771 } else {
1772 if (simpleParameters.get (1) == TREE || simpleParameters.get (1) == TOPOLOGY) {
1773 delete ((node_iterator<long>*)simpleParameters.get(2));
1774 }
1775 }
1776 }
1777
1778 if (simpleParameters.countitems() > 2) {
1779 simpleParameters.Delete (2);
1780 }
1781
1782
1783 HBLObjectRef iterator_substrate = _ProcessAnArgumentByType (*GetIthParameter(0UL), ASSOCIATIVE_LIST | MATRIX | TOPOLOGY | TREE, current_program, ¶meters);
1784
1785 _ElementaryCommand * advance_command = current_program.GetIthCommand(simpleParameters.get (0));
1786
1787 if (!advance_command || advance_command->code != HY_HBL_COMMAND_ADVANCE_ITERATOR) {
1788 throw _String ("Iterator init command is not linked with an advance command");
1789 }
1790
1791 advance_command->parameters.DeleteTail(advance_command->simpleParameters.get (1));
1792
1793 if (iterator_substrate->ObjectClass() == ASSOCIATIVE_LIST) {
1794 _AssociativeList* source_object = (_AssociativeList*)iterator_substrate;
1795 simpleParameters[1] = ASSOCIATIVE_LIST; // indicate that this is an ASSOCIATIVE_LIST
1796 simpleParameters << (long) new AVLListXLIterator (source_object->ListIterator ());
1797 } else { // MATRIX
1798 if (iterator_substrate->ObjectClass() == MATRIX) {
1799 simpleParameters[1] = MATRIX; // indicate that tjis is a MATRIX
1800 } else {
1801 simpleParameters[1] = iterator_substrate->ObjectClass();
1802 _TreeTopology * source_tree = (_TreeTopology*)iterator_substrate;
1803 simpleParameters << (long) new node_iterator<long> (&source_tree->GetRoot(), _HY_TREE_TRAVERSAL_POSTORDER);
1804 }
1805 }
1806 } catch (const _String& error) {
1807 return _DefaultExceptionHandler (nil, error, current_program);
1808 }
1809
1810 return true;
1811 }
1812
1813 //____________________________________________________________________________________
1814
HandleAdvanceIterator(_ExecutionList & current_program)1815 bool _ElementaryCommand::HandleAdvanceIterator(_ExecutionList& current_program) {
1816 /*
1817 parameters:
1818 [0-2] -> Id of the variable to store stuff in
1819 [+1] -> the source object
1820 [+2] -> _Variable to store the value in
1821 [+3] -> _Variable to store key/index in (could be nil)
1822 [+4] -> _Variable to store column index in (if the object is _Matrix)
1823
1824
1825
1826 simpleParameters:
1827 [0] -> index of the iterator init command in current_program
1828 [1] -> number of receptable variables
1829 [2] -> for _Matrix iterator, the row index
1830 [3] -> for _Matrix iterator, the column index
1831 */
1832 current_program.advance();
1833 try {
1834
1835 long reciever_count = simpleParameters.get (1);
1836
1837 _ElementaryCommand * init_command = current_program.GetIthCommand(simpleParameters.get (0));
1838 if (!init_command || init_command->code != HY_HBL_COMMAND_INIT_ITERATOR) {
1839 throw (_String ("Iterator advance command is not linked with an iterator initiaizer"));
1840 }
1841
1842 long source_object_class = init_command->simpleParameters.get (1);
1843
1844 bool first_in = false;
1845 if (source_object_class != MATRIX && reciever_count == 3) {
1846 throw (_String ("Iterators over dictionaries do not support the 3 argument form"));
1847 }
1848
1849 if (parameters.countitems () == reciever_count || parameters.GetItem(reciever_count) != init_command->parameters.GetItem(1)) {
1850 parameters.DeleteTail (reciever_count, true);
1851 parameters << init_command->parameters.GetItem(1);
1852 first_in = true;
1853 for (long k = 0L; k < reciever_count; k++) {
1854 parameters << _ValidateStorageVariable (current_program, k);
1855 }
1856 if (source_object_class == MATRIX) {
1857 if (simpleParameters.countitems() == 2) {
1858 simpleParameters << 0 << 0;
1859 } else {
1860 simpleParameters[2] = 0;
1861 simpleParameters[3] = 0;
1862 }
1863 }
1864 }
1865
1866 if (source_object_class == MATRIX) {
1867 _Matrix* source_object = (_Matrix*)parameters.GetItem (reciever_count);
1868 long row = simpleParameters.get (2),
1869 column = simpleParameters.get (3);
1870
1871 if (!first_in) {
1872 column ++;
1873 if (column >= source_object->GetVDim()) {
1874 column = 0;
1875 row ++;
1876 }
1877 }
1878
1879 if (row >= source_object->GetHDim() || source_object->GetVDim() == 0) {
1880 ((_Variable*)parameters.GetItem (reciever_count << 1))->SetValue (new _FString (kEndIteration), false, false, NULL); // end loop
1881 } else {
1882 ((_Variable*)parameters.GetItem (reciever_count << 1))->SetValue (source_object->GetMatrixCell(row, column), false, false, NULL);
1883 if (reciever_count == 2) {
1884 ((_Variable*)parameters.GetItem (reciever_count+1))->SetValue (new _Constant (row*source_object->GetVDim () + column), false, false, NULL);
1885 } else if (reciever_count == 3) {
1886 ((_Variable*)parameters.GetItem (reciever_count+1))->SetValue (new _Constant (row), false, false,NULL);
1887 ((_Variable*)parameters.GetItem (reciever_count+2))->SetValue (new _Constant (column ), false, false,NULL);
1888 }
1889 }
1890
1891 simpleParameters[2] = row;
1892 simpleParameters[3] = column;
1893
1894
1895
1896 } else {
1897 if (source_object_class == ASSOCIATIVE_LIST) {
1898 AVLListXLIterator * it = (AVLListXLIterator *)init_command->simpleParameters.get (2);
1899 if (first_in) {
1900 it->begin();
1901 } else {
1902 ++(*it);
1903 }
1904 if (!it->is_done()) { // iterator not finished
1905 AVLListXLIteratorKeyValue state = *(*it);
1906 state.get_object()->AddAReference();
1907 if (reciever_count > 1) {
1908 ((_Variable*)parameters.GetItem (reciever_count+2))->SetValue ((HBLObjectRef)state.get_object(), false, false, NULL);
1909 ((_Variable*)parameters.GetItem (reciever_count+1))->SetValue (new _FString (*state.get_key()), false, false,NULL);
1910 } else {
1911 ((_Variable*)parameters.GetItem (reciever_count+1))->SetValue ((HBLObjectRef)state.get_object(), false, false,NULL);
1912 }
1913 } else {
1914 ((_Variable*)parameters.GetItem (reciever_count << 1))->SetValue (new _FString (kEndIteration), false, false,NULL);
1915 // iterator done
1916 }
1917 } else {
1918 node_iterator<long> * tree_iterator = (node_iterator<long> *)init_command->simpleParameters.get (2);
1919 if (first_in) {
1920 if (simpleParameters.countitems() < 3) {
1921 simpleParameters << 0;
1922 } else {
1923 simpleParameters[2] = 0;
1924 }
1925 }
1926 node<long>*topTraverser = tree_iterator->Next();
1927 if (topTraverser && !topTraverser->is_root()) {
1928 _FString *node_name;
1929 if (source_object_class == TREE)
1930 node_name = new _FString (map_node_to_calcnode (topTraverser)->ContextFreeName());
1931 else
1932 node_name = new _FString (((_TreeTopology*)parameters.GetItem (1))->GetNodeName(topTraverser));
1933
1934 if (reciever_count > 1) {
1935 ((_Variable*)parameters.GetItem (reciever_count+2))->SetValue (node_name, false, false,NULL);
1936 ((_Variable*)parameters.GetItem (reciever_count+1))->SetValue (new _Constant (simpleParameters.get(2)), false, false,NULL);
1937 } else {
1938 ((_Variable*)parameters.GetItem (reciever_count+1))->SetValue (node_name, false, false,NULL);
1939 }
1940 simpleParameters[2]++;
1941 } else {
1942 ((_Variable*)parameters.GetItem (reciever_count << 1))->SetValue (new _FString (kEndIteration), false, false,NULL);
1943 }
1944
1945 }
1946 }
1947 } catch (const _String& error) {
1948 return _DefaultExceptionHandler (nil, error, current_program);
1949 }
1950
1951 return true;
1952 }
1953
1954
1955 //____________________________________________________________________________________
HandleRequireVersion(_ExecutionList & current_program)1956 bool _ElementaryCommand::HandleRequireVersion(_ExecutionList& current_program){
1957 current_program.advance();
1958
1959 try {
1960 _String requested_version = _ProcessALiteralArgument (*GetIthParameter(0UL),current_program);
1961
1962 auto throw_error = [&] (void) -> void {
1963 throw _String ("Current script requires at least version ")& requested_version &" of HyPhy. Had version " &kHyPhyVersion & ". Please download an updated version from http://www.hyphy.org or github.com/veg/hyphy and try again.";
1964 };
1965
1966 const _List local_version = kHyPhyVersion.Tokenize ("."),
1967 required_version = requested_version.Tokenize(".");
1968
1969 unsigned long const upper_bound = MIN (local_version.countitems (), required_version.countitems());
1970
1971
1972 for (unsigned long i = 0UL; i < upper_bound; i++) {
1973 hyFloat local_number = ((_String*)local_version.GetItem(i))->to_float(),
1974 required_number = ((_String*)required_version.GetItem(i))->to_float();
1975
1976 if (local_number > required_number) {
1977 return true;
1978 }
1979 if (local_number < required_number) {
1980 throw_error();
1981 }
1982 }
1983
1984 if (required_version.countitems() <= upper_bound) {
1985 return true;
1986 }
1987 throw_error ();
1988
1989 } catch (const _String& error) {
1990 return _DefaultExceptionHandler (nil, error, current_program);
1991 }
1992 return true;
1993 }
1994
1995
1996
1997 //____________________________________________________________________________________
1998
HandleDeleteObject(_ExecutionList & current_program)1999 bool _ElementaryCommand::HandleDeleteObject(_ExecutionList& current_program){
2000
2001 current_program.advance();
2002
2003 const static _String kShallow (":shallow");
2004
2005 bool do_shallow = parameter_count() >= 2 || *GetIthParameter(parameter_count()-1, false) == kShallow;
2006
2007 for (unsigned long i = 0UL; i < parameter_count(); i++) {
2008 long requested_type = HY_BL_LIKELIHOOD_FUNCTION,
2009 object_index = kNotFound;
2010 BaseRef source_object = _HYRetrieveBLObjectByNameMutable (AppendContainerName(*GetIthParameter(i),current_program.nameSpacePrefix), requested_type,&object_index,false);
2011
2012 if (source_object) {
2013 KillLFRecord (object_index,!do_shallow);
2014 } else {
2015 ReportWarning(GetIthParameter(i)->Enquote() & " is not a supported argument type for " & _HY_ValidHBLExpressions.RetrieveKeyByPayload(get_code()));
2016 }
2017 }
2018 return true;
2019 }
2020
2021 //____________________________________________________________________________________
2022
HandleClearConstraints(_ExecutionList & current_program)2023 bool _ElementaryCommand::HandleClearConstraints(_ExecutionList& current_program){
2024 current_program.advance();
2025
2026 for (unsigned long i = 0UL; i< parameter_count(); i++) {
2027
2028 _String source_name = AppendContainerName (*(_String*)parameters(i), current_program.nameSpacePrefix);
2029 _Variable *clear_me = (_Variable*) FetchVar (LocateVarByName (source_name));
2030 if (clear_me) { // variable exists
2031 clear_me->ClearConstraints();
2032 } else {
2033 ReportWarning(GetIthParameter(i)->Enquote() & " is not an existing variable in call to " & _HY_ValidHBLExpressions.RetrieveKeyByPayload(get_code()));
2034 }
2035 }
2036 return true;
2037 }
2038
2039 //____________________________________________________________________________________
2040
HandleGetURL(_ExecutionList & current_program)2041 bool _ElementaryCommand::HandleGetURL(_ExecutionList& current_program){
2042 const static _String save_to_file_action ("SAVE_TO_FILE");
2043
2044 current_program.advance();
2045 _Variable * receptacle = nil;
2046
2047 try {
2048 _String url = _ProcessALiteralArgument (*GetIthParameter(1UL),current_program),
2049 *action = GetIthParameter(2UL, false);
2050
2051 if (!action) { // store URL contents in a variable
2052 receptacle = _ValidateStorageVariable (current_program);
2053 if (Get_a_URL(url)) {
2054 receptacle->SetValue(new _FString (url,false),false,true,NULL);
2055 } else {
2056 throw (_String ("Could not fetch ") & url.Enquote());
2057 }
2058 } else {
2059 if (*action == save_to_file_action) {
2060 _String file_name = _ProcessALiteralArgument (*GetIthParameter(1UL),current_program);
2061 if (!ProcessFileName (file_name, true,true,(hyPointer)current_program.nameSpacePrefix,false,¤t_program)) {
2062 return false;
2063 }
2064 if (!Get_a_URL(url, &file_name)) {
2065 throw (_String ("Could not fetch ") & url.Enquote());
2066 }
2067 } else {
2068 throw (_String("Unknown action ") & action->Enquote());
2069 }
2070 }
2071
2072
2073 } catch (const _String& error) {
2074 return _DefaultExceptionHandler (nil, error, current_program);
2075 }
2076 return true;
2077 }
2078
2079 //____________________________________________________________________________________
2080
HandleAssert(_ExecutionList & current_program)2081 bool _ElementaryCommand::HandleAssert (_ExecutionList& current_program) {
2082 current_program.advance();
2083 static const _String kBreakpointTrap = ("__SIGTRAP__");
2084
2085 try {
2086 _Formula parsed_expression;
2087
2088 if (kBreakpointTrap == * GetIthParameter(0) ) {
2089 raise(SIGTRAP);
2090 return true;
2091 }
2092
2093 _CheckExpressionForCorrectness (parsed_expression, *GetIthParameter(0UL), current_program, NUMBER);
2094 if (CheckEqual (parsed_expression.Compute()->Value (), 0.0)) { // assertion failed
2095 bool soft_assertions = hy_env::EnvVariableTrue(hy_env::assertion_behavior);
2096 _String assertion_feedback;
2097 _String * custom_error_message = GetIthParameter(1UL,false);
2098 if (custom_error_message) {
2099 assertion_feedback = _ProcessALiteralArgument(*custom_error_message, current_program);
2100 } else {
2101 assertion_feedback = _String("Assertion ") & GetIthParameter(0UL)->Enquote() & " failed.";
2102 }
2103 if (soft_assertions) {
2104 StringToConsole (assertion_feedback);
2105 NLToConsole();
2106 current_program.GoToLastInstruction ();
2107 } else {
2108 current_program.ReportAnExecutionError (assertion_feedback);
2109 }
2110 }
2111 } catch (const _String& error) {
2112 return _DefaultExceptionHandler (nil, error, current_program);
2113 }
2114 return true;
2115 }
2116
2117
2118 //____________________________________________________________________________________
2119
HandleSelectTemplateModel(_ExecutionList & current_program)2120 bool _ElementaryCommand::HandleSelectTemplateModel (_ExecutionList& current_program) {
2121 static _String last_model_used, kPromptText ("Select a standard model");
2122
2123 current_program.advance();
2124 try {
2125 _String source_name = *GetIthParameter(0UL);
2126 if (source_name == hy_env::use_last_model) {
2127 if (last_model_used.nonempty()) {
2128 PushFilePath (last_model_used);
2129 } else {
2130 throw ( hy_env::use_last_model &" cannot be used before any models have been defined.");
2131 }
2132 } else {
2133 ReadModelList();
2134
2135 long object_type = HY_BL_DATASET|HY_BL_DATASET_FILTER;
2136 _DataSetFilter const * source_filter = (_DataSetFilter const*)_GetHBLObjectByType(source_name, object_type, nil, ¤t_program);
2137
2138
2139 _String data_type;
2140 unsigned long unit_length = source_filter->GetUnitLength();
2141
2142 _TranslationTable const* filter_table = source_filter->GetData()->GetTT();
2143
2144 if (unit_length==1UL) {
2145 if (filter_table->IsStandardNucleotide()) {
2146 data_type = "nucleotide";
2147 } else if (filter_table->IsStandardAA()) {
2148 data_type = "aminoacid";
2149 }
2150 } else {
2151 if (filter_table->IsStandardNucleotide()) {
2152 if (unit_length==3UL) {
2153 data_type = "codon";
2154 } else {
2155 if (unit_length==2UL)
2156 data_type = "dinucleotide";
2157 }
2158 }
2159 }
2160
2161 if (data_type.empty()) {
2162 throw (source_name.Enquote () & " contains non-standard data and template models can't be selected on it");
2163 }
2164
2165 _SimpleList matching_models;
2166
2167 for (unsigned long model_index = 0; model_index < templateModelList.lLength; model_index++) {
2168 _List *model_components = (_List*)templateModelList(model_index);
2169
2170 if (data_type == *(_String*)model_components->GetItem(3)) {
2171 _String * dim = (_String*)model_components->GetItem(2);
2172 if (*dim== _String("*")|| source_filter->GetDimension() == dim->to_long()) {
2173 matching_models << model_index;
2174 }
2175 }
2176 }
2177
2178 if (matching_models.empty()) {
2179 throw (source_name.Enquote () & " could not be matched with any template models");
2180 }
2181
2182 long model_id = kNotFound;
2183 bool need_to_prompt_user = true;
2184
2185 if (current_program.has_stdin_redirect() || current_program.has_keyword_arguments()) {
2186 _String option;
2187 bool has_input = true;
2188 try {
2189 _FString * redirect = (_FString*)hy_env::EnvVariableGet(hy_env::fprintf_redirect, STRING);
2190 option = current_program.FetchFromStdinRedirect (&kPromptText, false, !(redirect&&redirect->has_data()));
2191 } catch (const _String& e) {
2192 if (e != kNoKWMatch) {
2193 throw (e);
2194 }
2195 has_input = false;
2196 }
2197 if (has_input) {
2198 model_id = matching_models.FindOnCondition( [&] (long index, unsigned long) -> bool {
2199 return option == * (_String*) templateModelList.GetItem (index,0);
2200 });
2201
2202
2203 if (model_id == kNotFound) {
2204 throw (option.Enquote() & " is not a valid model (with input redirect)");
2205 }
2206
2207 need_to_prompt_user = false;
2208 }
2209 }
2210
2211 if (need_to_prompt_user) {
2212 #ifdef __HEADLESS__
2213 throw _String("Unhandled standard input interaction in SelectTemplateModel for headless HyPhy");
2214 #endif
2215
2216
2217
2218 for (int i = 0; i < kMaxDialogPrompts; i++) {
2219 printf ("\n\n +--------------------------+\n");
2220 printf ( " | %s. |\n", kPromptText.get_str());
2221 printf ( " +--------------------------+\n\n\n");
2222
2223 for (model_id = 0; model_id<matching_models.lLength; model_id++) {
2224 printf ("\n\t(%s):%s",((_String*)templateModelList.GetItem(matching_models(model_id),0))->get_str(),
2225 ((_String*)templateModelList.GetItem(matching_models(model_id),1))->get_str());
2226 }
2227 printf ("\n\n Please type in the abbreviation for the model you want to use:");
2228 _String const user_choice = StringFromConsole();
2229
2230 model_id = matching_models.FindOnCondition( [&] (long index, unsigned long) -> bool {
2231 return user_choice.EqualIgnoringCase(*(_String*) templateModelList.GetItem (index,0));
2232 });
2233
2234 if (model_id != kNotFound) {
2235 break;
2236 }
2237 }
2238
2239 if (model_id == kNotFound) {
2240 throw _String("Dialog did not return a valid choice after maximum allowed number of tries");
2241 return false;
2242 }
2243
2244 }
2245
2246 _String model_file = GetStandardDirectory (HY_HBL_DIRECTORY_TEMPLATE_MODELS) & *(_String*)templateModelList.GetItem(matching_models(model_id),4);
2247
2248 _ExecutionList std_model;
2249 PushFilePath (model_file, false);
2250 ReadBatchFile (model_file,std_model);
2251 PopFilePath ();
2252 last_model_used = model_file;
2253 std_model.Execute (¤t_program);
2254
2255 }
2256 } catch (const _String& error) {
2257 return _DefaultExceptionHandler (nil, error, current_program);
2258 }
2259
2260 return true;
2261
2262 }
2263
2264 //____________________________________________________________________________________
2265
HandleMolecularClock(_ExecutionList & current_program)2266 bool _ElementaryCommand::HandleMolecularClock(_ExecutionList& current_program){
2267 current_program.advance();
2268 try {
2269 _CalcNode * apply_clock_here = (_CalcNode *)_CheckForExistingVariableByType (*GetIthParameter(0), current_program, TREE | TREE_NODE);
2270 _TheTree * parent_tree;
2271
2272 _String clock_base;
2273
2274 if (apply_clock_here->ObjectClass() == TREE_NODE) {
2275 parent_tree = (_TheTree*)((_VariableContainer*)apply_clock_here)->GetTheParent();
2276 if (!parent_tree) {
2277 throw (_String("Internal error - orphaned tree node ") & apply_clock_here->GetName ()->Enquote());
2278 }
2279 clock_base = apply_clock_here->GetName() -> Cut (parent_tree->GetName()->length() + 1, kStringEnd);
2280 } else {
2281 parent_tree = (_TheTree*)apply_clock_here;
2282 }
2283
2284 parent_tree->MolecularClock(clock_base,parameters);
2285
2286 } catch (const _String& error) {
2287 return _DefaultExceptionHandler (nil, error, current_program);
2288 }
2289 return true;
2290 }
2291
2292 //____________________________________________________________________________________
2293
HandleMPIReceive(_ExecutionList & current_program)2294 bool _ElementaryCommand::HandleMPIReceive (_ExecutionList& current_program){
2295 current_program.advance();
2296 _Variable * receptacle = nil;
2297
2298 try {
2299
2300 #ifdef __HYPHYMPI__
2301
2302 receptacle = _ValidateStorageVariable (current_program, 2UL);
2303 _Variable* node_index_storage = _ValidateStorageVariable (current_program, 1UL);
2304
2305 long target_node = _ProcessNumericArgumentWithExceptions(*GetIthParameter(0UL), current_program.nameSpacePrefix),
2306 node_count = hy_env::EnvVariableGetNumber(hy_env::mpi_node_count);
2307
2308 if (target_node < -1L || target_node >= node_count) {
2309 throw (GetIthParameter(1UL)->Enquote () & " (=" & node_count & ") must be a valid MPI node index (or -1 to accept from any node");
2310 }
2311
2312 long received_from;
2313 receptacle->SetValue(new _FString (MPIRecvString (target_node,received_from)), false, true,NULL);
2314 node_index_storage->SetValue (new _Constant (received_from), false, true, NULL);
2315
2316 #else
2317 throw _String("Command not supported for non-MPI versions of HyPhy. HBL scripts need to check for MPI before calling MPI features");
2318 #endif
2319
2320 } catch (const _String& error) {
2321 return _DefaultExceptionHandler (receptacle, error, current_program);
2322 }
2323
2324 return true;
2325 }
2326
2327 //____________________________________________________________________________________
2328
HandleMPISend(_ExecutionList & current_program)2329 bool _ElementaryCommand::HandleMPISend (_ExecutionList& current_program){
2330 current_program.advance();
2331
2332 _List dynamic_variable_manager;
2333
2334 try {
2335
2336 #ifdef __HYPHYMPI__
2337 long target_node = _ProcessNumericArgumentWithExceptions(*GetIthParameter(0UL), current_program.nameSpacePrefix),
2338 node_count = hy_env::EnvVariableGetNumber(hy_env::mpi_node_count);
2339
2340 if (target_node < 0L || target_node >= node_count) {
2341 throw (GetIthParameter(1UL)->Enquote () & " (=" & node_count & ") is not a valid MPI node index; valud range is " & target_node & " to " & (node_count-1));
2342 }
2343
2344 _StringBuffer message_to_send (1024UL);
2345
2346 if (parameter_count() > 2) { // this is the case of MPISend (node, filepath, {input option}
2347 _AssociativeList * arguments = (_AssociativeList*)_ProcessAnArgumentByType(*GetIthParameter(2), ASSOCIATIVE_LIST, current_program, &dynamic_variable_manager);
2348 _String file_path = *GetIthParameter(1UL);
2349 if (! ProcessFileName(file_path,false,true,(hyPointer)current_program.nameSpacePrefix)) {
2350 throw (GetIthParameter(1UL)->Enquote() & " is an ivalid file path");
2351 }
2352 (message_to_send << _HY_ValidHBLExpressions.RetrieveKeyByPayload(HY_HBL_COMMAND_EXECUTE_A_FILE) << file_path.Enquote() << ',')
2353 .AppendNewInstance(arguments->Serialize(0UL)) << ");";
2354 } else {
2355 // is this a likelihood function?
2356 long type = HY_BL_LIKELIHOOD_FUNCTION;
2357 try {
2358 ((_LikelihoodFunction*) _GetHBLObjectByTypeMutable(AppendContainerName (*GetIthParameter(1), current_program.nameSpacePrefix), type))->SerializeLF(message_to_send, _hyphyLFSerializeModeOptimize);
2359 } catch (const _String &) {
2360 // catch literal cases here
2361 message_to_send = _ProcessALiteralArgument(*GetIthParameter(1UL), current_program);
2362 }
2363 }
2364
2365 if (message_to_send.nonempty()) {
2366 MPISendString(message_to_send, target_node);
2367 } else {
2368 throw (_String ("An invalid (empty) MPI message"));
2369 }
2370 #else
2371 throw _String("Command not supported for non-MPI versions of HyPhy. HBL scripts need to check for MPI before calling MPI features");
2372 #endif
2373
2374 } catch (const _String& error) {
2375 return _DefaultExceptionHandler (nil, error, current_program);
2376 }
2377
2378 return true;
2379 }
2380
2381
2382
2383 //____________________________________________________________________________________
2384
HandleSetParameter(_ExecutionList & current_program)2385 bool _ElementaryCommand::HandleSetParameter (_ExecutionList& current_program) {
2386
2387 static const _String kBGMNodeOrder ("BGM_NODE_ORDER"),
2388 kBGMGraph ("BGM_GRAPH_MATRIX"),
2389 kBGMScores ("BGM_SCORE_CACHE"),
2390 kBGMConstraintMx ("BGM_CONSTRAINT_MATRIX"),
2391 kBGMParameters ("BGM_NETWORK_PARAMETERS");
2392
2393 current_program.advance();
2394
2395 _List dynamic_variable_manager;
2396
2397 try {
2398 _String object_to_change = *GetIthParameter(0UL);
2399
2400 /* handle special cases */
2401 if (object_to_change == hy_env::random_seed) {
2402 hy_env::EnvVariableSet(hy_env::random_seed, new _Constant (hy_random_seed = _ProcessNumericArgumentWithExceptions (*GetIthParameter(1),current_program.nameSpacePrefix) ), false);
2403 init_genrand(hy_random_seed);
2404 return true;
2405 }
2406
2407 if (object_to_change == hy_env::defer_constrain_assignment) {
2408 bool defer_status = _ProcessNumericArgumentWithExceptions (*GetIthParameter(1),current_program.nameSpacePrefix);
2409 if (defer_status) {
2410 deferSetFormula = new _SimpleList;
2411 deferClearConstraint = new _AVLList (new _SimpleList);
2412 } else if (deferSetFormula || deferClearConstraint){
2413 FinishDeferredSF ();
2414 }
2415 return true;
2416 }
2417
2418 if (object_to_change == hy_env::execution_mode) {
2419 current_program.errorHandlingMode = _ProcessNumericArgumentWithExceptions (*GetIthParameter(1),current_program.nameSpacePrefix);
2420 return true;
2421 }
2422
2423 if (object_to_change == hy_env::status_bar_update_string) {
2424 SetStatusLineUser (_ProcessALiteralArgument (*GetIthParameter(1), current_program));
2425 return true;
2426 }
2427
2428
2429 const _String source_name = AppendContainerName (*GetIthParameter(0), current_program.nameSpacePrefix);
2430
2431 long object_type = HY_BL_ANY,
2432 object_index;
2433
2434 BaseRef source_object;
2435 _String set_this_attribute = *GetIthParameter(1UL);
2436
2437 try {
2438 source_object = _GetHBLObjectByTypeMutable (source_name, object_type, &object_index);
2439 } catch (const _String& error) { // handle cases when the source is not an HBL object
2440
2441 _CalcNode* tree_node = nil;
2442
2443
2444 if (source_name.IsValidIdentifier(fIDAllowFirstNumeric|fIDAllowCompound)) {
2445 tree_node = (_CalcNode*)FetchObjectFromVariableByType(&source_name, TREE_NODE);
2446 } else{
2447 _String converted = source_name.ConvertToAnIdent(fIDAllowFirstNumeric|fIDAllowCompound);
2448 tree_node = (_CalcNode*)FetchObjectFromVariableByType(&converted, TREE_NODE);
2449 }
2450
2451 if (tree_node) {
2452 if (set_this_attribute == _String("MODEL")) {
2453 _String model_name = AppendContainerName(*GetIthParameter(2UL),current_program.nameSpacePrefix);
2454 long model_type = HY_BL_MODEL, model_index;
2455 _Matrix* model_object = (_Matrix*)_GetHBLObjectByTypeMutable(model_name, model_type, &model_index);
2456 _TheTree * parent_tree = (_TheTree * )tree_node->ParentTree();
2457 if (!parent_tree) {
2458 throw (GetIthParameter(0UL)->Enquote() & " is an orphaned tree node (the parent tree has been deleted)");
2459 }
2460 tree_node->ReplaceModel (model_name, parent_tree);
2461 long partition_id, likelihood_function_id = ((_TheTree*)parent_tree->Compute())->IsLinkedToALF(partition_id);
2462 if (likelihood_function_id>=0){
2463 //throw(parent_tree->GetName()->Enquote() & " is linked to a likelihood function (" & *GetObjectNameByType (HY_BL_LIKELIHOOD_FUNCTION, likelihood_function_id) &") and cannot be modified ");
2464 //return false;
2465 ((_LikelihoodFunction*)likeFuncList (likelihood_function_id))->Rebuild(true);
2466 }
2467 } else {
2468 throw (set_this_attribute.Enquote() & " is not a supported parameter type for a tree node argument");
2469 }
2470 } else {
2471 throw (GetIthParameter(0UL)->Enquote() & " is not a supported object type");
2472 }
2473 return true;
2474 }
2475
2476
2477 switch (object_type) {
2478 case HY_BL_BGM: { // BGM Branch
2479
2480 _BayesianGraphicalModel * bgm = (_BayesianGraphicalModel *) source_object;
2481 long num_nodes = bgm->GetNumNodes();
2482
2483
2484 // set data matrix
2485 if (set_this_attribute == kBGMData) {
2486 _Matrix * data_mx = (_Matrix *) _ProcessAnArgumentByType(*GetIthParameter(2UL), MATRIX, current_program, &dynamic_variable_manager);
2487
2488 if (data_mx->GetVDim() == num_nodes) {
2489 bgm->SetDataMatrix (data_mx);
2490 } else {
2491 throw (_String("Data matrix columns (") & data_mx->GetVDim() & " ) does not match number of nodes in graph (" & num_nodes & ")");
2492 }
2493
2494 }
2495 else if (set_this_attribute == kBGMScores) {
2496 bgm->ImportCache((_AssociativeList *)_ProcessAnArgumentByType(*GetIthParameter(2UL), ASSOCIATIVE_LIST, current_program, &dynamic_variable_manager));
2497 } // set structure to user-specified adjacency matrix
2498 else if (set_this_attribute == kBGMGraph) {
2499 _Matrix * graphMx = (_Matrix *) _ProcessAnArgumentByType(*GetIthParameter(2UL), MATRIX, current_program, &dynamic_variable_manager);
2500
2501 if (graphMx->check_dimension(num_nodes, num_nodes)) {
2502 bgm->SetStructure ((_Matrix *) graphMx->makeDynamic());
2503 } else {
2504 throw _String("Dimension of graph does not match current graph");
2505 }
2506 } // set constraint matrix
2507 else if (set_this_attribute == kBGMConstraintMx) {
2508 _Matrix * constraint_mx = (_Matrix *) _ProcessAnArgumentByType(*GetIthParameter(2UL), MATRIX, current_program, &dynamic_variable_manager);
2509 if (constraint_mx->check_dimension(num_nodes, num_nodes)) {
2510 bgm->SetConstraints ((_Matrix *) constraint_mx->makeDynamic());
2511 } else {
2512 throw _String("Dimensions of constraint matrix do not match current graph");
2513 }
2514 } // set node order
2515 else if (set_this_attribute == kBGMNodeOrder) {
2516 _Matrix * order_mx = (_Matrix *) _ProcessAnArgumentByType(*GetIthParameter(2UL), MATRIX, current_program, &dynamic_variable_manager);
2517
2518 if (order_mx->check_dimension(1,num_nodes)) {
2519
2520 _SimpleList order_list;
2521 order_mx->ConvertToSimpleList(order_list);
2522
2523 bgm->SetNodeOrder ( &order_list );
2524 } else {
2525 throw _String("Order must be a row vector whose dimension matches the number of nodes in graph");
2526 }
2527 } else {
2528 throw (GetIthParameter (2UL)->Enquote() & " is not a valid parameter for BGM objects");
2529 }
2530 } // end BGM
2531 break;
2532
2533 case HY_BL_SCFG:
2534 case HY_BL_LIKELIHOOD_FUNCTION: {
2535
2536
2537 if (object_type == HY_BL_SCFG && set_this_attribute == hy_env::kSCFGCorpus) {
2538 HBLObjectRef corpus_source = _ProcessAnArgumentByType (*GetIthParameter(2UL), MATRIX|STRING, current_program, &dynamic_variable_manager);
2539 if (corpus_source->ObjectClass () == STRING) {
2540 _List single_string ( new _String (((_FString*)corpus_source)->get_str()));
2541 _Matrix wrapper (single_string);
2542 ((Scfg*)source_object)->SetStringCorpus (&wrapper);
2543 } else {
2544 _Matrix * matrix_corpus = (_Matrix*)corpus_source;
2545 if (matrix_corpus->IsAStringMatrix()) {
2546 ((Scfg*)source_object)->SetStringCorpus (matrix_corpus);
2547 } else {
2548 throw (set_this_attribute.Enquote() & " did not evaluate to a matrix of strings");
2549 }
2550 }
2551 } else {
2552 _LikelihoodFunction * lkf = (_LikelihoodFunction *) source_object;
2553 long parameter_index = _ProcessNumericArgumentWithExceptions (set_this_attribute ,current_program.nameSpacePrefix);
2554 if (lkf->GetIndependentVars().Map (parameter_index) < 0L) {
2555 throw (GetIthParameter(1)->Enquote() & " (=" & parameter_index & ") is not a valid parameter index");
2556 }
2557 lkf->SetIthIndependent (parameter_index,_ProcessNumericArgumentWithExceptions (*GetIthParameter(2),current_program.nameSpacePrefix));
2558 }
2559 } // end HY_BL_SCFG; HY_BL_LIKELIHOOD_FUNCTION
2560 break;
2561 case HY_BL_DATASET:
2562 case HY_BL_DATASET_FILTER: {
2563
2564 if (object_type == HY_BL_DATASET_FILTER) {
2565 ReleaseDataFilterLock (object_index);
2566 }
2567
2568 long sequence_index = _ProcessNumericArgumentWithExceptions (*GetIthParameter(1),current_program.nameSpacePrefix);
2569 _DataSet * ds = nil;
2570 if (object_type == HY_BL_DATASET) {
2571 ds = (_DataSet*) source_object;
2572 }
2573 else {
2574 _DataSetFilter *dsf = (_DataSetFilter*)source_object;
2575 ds = dsf->GetData ();
2576 sequence_index = dsf->theNodeMap.Map (sequence_index);
2577 }
2578
2579
2580 _String * sequence_name = new _String (_ProcessALiteralArgument (*GetIthParameter(2UL), current_program));
2581
2582 if (! ds->SetSequenceName (sequence_index, sequence_name)) {
2583 delete sequence_name;
2584 throw (GetIthParameter (1UL)->Enquote() & " (=" & sequence_index & ") is not a valid sequence index");
2585 }
2586
2587 } // end HY_BL_DATASET; HY_BL_DATASET_FILTER
2588 break;
2589
2590 } // end switch
2591
2592
2593 } catch (const _String& error) {
2594 return _DefaultExceptionHandler (nil, error, current_program);
2595 }
2596 return true;
2597 }
2598
2599
2600 //____________________________________________________________________________________
2601
HandleFprintf(_ExecutionList & current_program)2602 bool _ElementaryCommand::HandleFprintf (_ExecutionList& current_program) {
2603
2604 static const _String kFprintfStdout ("stdout"),
2605 kFprintfMessagesLog ("MESSAGE_LOG"),
2606 kFprintfClearFile ("CLEAR_FILE"),
2607 kFprintfKeepOpen ("KEEP_OPEN"),
2608 kFprintfCloseFile ("CLOSE_FILE"),
2609 kFprintfSystemVariableDump ("LIST_ALL_VARIABLES"),
2610 kFprintfSelfDump ("PRINT_SELF");
2611
2612
2613 static _List _open_file_handles_aux;
2614 static _AVLListX open_file_handles (&_open_file_handles_aux);
2615
2616
2617 current_program.advance();
2618
2619 bool do_close = true,
2620 print_to_stdout = false,
2621 skip_file_path_eval = false,
2622 success = true;
2623
2624
2625 FILE* destination_file = nil;
2626
2627 try {
2628
2629 _String destination = *GetIthParameter(0UL);
2630
2631 if (destination == kFprintfStdout) {
2632 _FString * redirect = (_FString*)hy_env::EnvVariableGet(hy_env::fprintf_redirect, STRING);
2633 if (redirect && redirect->has_data()) {
2634 destination = redirect->get_str();
2635 if (destination == hy_env::kDevNull) {
2636 return true; // "print" to /dev/null
2637 } else {
2638 skip_file_path_eval = true;
2639 }
2640 }
2641 else {
2642 print_to_stdout = true;
2643 }
2644 } else {
2645 if (destination == kPromptForFilePlaceholder) {
2646 skip_file_path_eval = true;
2647 }
2648 }
2649
2650 print_digit_specification = hy_env::EnvVariableGetDefaultNumber (hy_env::print_float_digits);
2651
2652 if (!print_to_stdout) {
2653 if (destination == kFprintfMessagesLog) {
2654 if ((destination_file = hy_message_log_file) == nil) {
2655 return true; // requested print to MESSAGE_LOG, but it does not exist
2656 // (e.g. remote MPI nodes, or running from a read only location
2657 }
2658 } else {
2659 if (skip_file_path_eval == false) {
2660 destination = _ProcessALiteralArgument (destination,current_program);
2661 }
2662
2663 if (!ProcessFileName(destination, true,false,(hyPointer)current_program.nameSpacePrefix, false, ¤t_program)) {
2664 return false;
2665 }
2666
2667 long open_handle = open_file_handles.Find (&destination);
2668
2669 do_close = open_handle < 0;
2670
2671 if (!do_close) {
2672 destination_file = (FILE*)open_file_handles.GetXtra (open_handle);
2673 } else {
2674 if ((destination_file = doFileOpen (destination.get_str(), "a")) == nil)
2675 throw (_String ("Could not create/open output file at path ") & destination.Enquote() & ".");
2676 }
2677 }
2678 }
2679
2680 for (unsigned long print_argument_idx = 1UL; print_argument_idx < parameter_count (); print_argument_idx++) {
2681 _String * current_argument = GetIthParameter(print_argument_idx);
2682 BaseRef managed_object_to_print = nil,
2683 dynamic_object_to_print = nil;
2684
2685 // handle special cases first
2686 if (*current_argument == kFprintfClearFile) {
2687 if (!print_to_stdout && destination_file) {
2688 fclose (destination_file);
2689 destination_file = doFileOpen (destination.get_str(), "w");
2690 if (!do_close) {
2691 _String* destination_copy = new _String (destination);
2692
2693 if (open_file_handles.UpdateValue(destination_copy, (long)destination_file, 1) >= 0) { // did not insert value
2694 DeleteObject (destination_copy);
2695 }
2696 }
2697
2698 }
2699 } else if (*current_argument == kFprintfKeepOpen) {
2700 if (!print_to_stdout) {
2701 open_file_handles.Insert (new _String (destination), (long)destination_file, false, true);
2702 do_close = false;
2703 }
2704 } else if (*current_argument == kFprintfCloseFile) {
2705 open_file_handles.Delete (&destination, true);
2706 do_close = true;
2707 } else if (*current_argument == kFprintfSystemVariableDump ) {
2708 managed_object_to_print = &variableNames;
2709 } else if (*current_argument == kFprintfSelfDump) {
2710 managed_object_to_print = ¤t_program;
2711 } else {
2712 _String namespaced_id = AppendContainerName (*current_argument, current_program.nameSpacePrefix);
2713
2714 // first check if the argument is a existing HBL object
2715 // TODO: this will go away in v3 when everything is in the same namespace
2716
2717
2718 if (!managed_object_to_print) { // not an existing variable
2719 long object_type = HY_BL_ANY, object_index;
2720 managed_object_to_print = //_GetHBLObjectByTypeMutable (namespaced_id, object_type, &object_index, false);
2721 _HYRetrieveBLObjectByNameMutable (namespaced_id, object_type,&object_index,false, false);
2722 if (managed_object_to_print) {
2723 if (object_type == HY_BL_DATASET_FILTER) {
2724 ReleaseDataFilterLock(object_index);
2725 }
2726 } else {
2727 dynamic_object_to_print = _ProcessAnArgumentByType (*current_argument, HY_ANY_OBJECT, current_program, nil);
2728 }
2729 }
2730 }
2731
2732 BaseRef printables [2] = {managed_object_to_print, dynamic_object_to_print};
2733 for (BaseRef obj : printables) {
2734 if (obj) {
2735 if (!print_to_stdout) {
2736 obj->toFileStr (destination_file);
2737 } else {
2738 StringToConsole (_String((_String*)obj->toStr()));
2739 }
2740 break;
2741 }
2742 }
2743 DeleteObject(dynamic_object_to_print);
2744
2745 } // end for
2746
2747 }
2748 catch (const _String& error) {
2749 if (hy_env::EnvVariableTrue(hy_env::soft_fileio_exceptions)) {
2750 hy_env::EnvVariableSet(hy_env::last_fileio_exception, new _FString (error,false), false);
2751 success = true;
2752 } else {
2753 success = _DefaultExceptionHandler (nil, error, current_program);
2754 }
2755 }
2756
2757 if (destination_file && destination_file != hy_message_log_file && do_close) {
2758 fclose (destination_file);
2759 }
2760
2761 return success;
2762 }
2763
2764 //____________________________________________________________________________________
2765
HandleExecuteCommandsCases(_ExecutionList & current_program,bool do_load_from_file,bool do_load_library)2766 bool _ElementaryCommand::HandleExecuteCommandsCases(_ExecutionList& current_program, bool do_load_from_file, bool do_load_library) {
2767 current_program.advance ();
2768 _StringBuffer * source_code = nil;
2769 bool pop_path = false;
2770
2771 auto cleanup = [&] () -> void {
2772 if (pop_path) {
2773 PopFilePath();
2774 }
2775 DeleteObject(source_code);
2776 };
2777
2778 _List dynamic_reference_manager;
2779
2780 try {
2781 bool has_redirected_input = false,
2782 has_user_kwargs = false;
2783
2784 _List _aux_argument_list;
2785 _AVLListXL argument_list (&_aux_argument_list);
2786 _AssociativeList * user_kwargs = nil;
2787
2788
2789 if (do_load_from_file) {
2790 _String file_path (*GetIthParameter(0UL));
2791
2792 if (file_path == kPromptForFilePlaceholder ){
2793 ProcessFileName (file_path, false, false, (hyPointer)current_program.nameSpacePrefix);
2794 } else {
2795 file_path = _ProcessALiteralArgument(file_path, current_program);
2796 }
2797 _String original_path (file_path);
2798
2799 FILE * source_file = nil;
2800
2801 bool reload = hy_env::EnvVariableTrue(hy_env::always_reload_libraries);
2802
2803
2804 if (do_load_library) {
2805 bool has_extension = file_path.FindBackwards (".",0,-1) != kNotFound;
2806
2807
2808 for (unsigned long p = 0; !source_file && p < _hy_standard_library_paths.countitems(); p++) {
2809 for (unsigned long e = 0; !source_file && e < _hy_standard_library_extensions.countitems(); e++) {
2810 _String try_path = *((_String*)_hy_standard_library_paths(p)) & file_path & *((_String*)_hy_standard_library_extensions(e));
2811
2812 ProcessFileName (try_path, false, false, (hyPointer)current_program.nameSpacePrefix, false, nil, false, true);
2813
2814 if (loadedLibraryPaths.Find(&try_path) >= 0 && parameter_count() == 2UL && !reload) {
2815 ReportWarning (_String("Already loaded ") & original_path.Enquote() & " from " & try_path);
2816 return true;
2817 }
2818 if ((source_file = doFileOpen (try_path.get_str (), "rb"))) {
2819 file_path = try_path;
2820 break;
2821 }
2822 if (has_extension) {
2823 break;
2824 }
2825 }
2826 }
2827 }
2828
2829 if (source_file == nil) {
2830 ProcessFileName (file_path, false,false,(hyPointer)current_program.nameSpacePrefix);
2831
2832 if (do_load_library && loadedLibraryPaths.Find(&file_path) >= 0 && parameter_count() == 2UL && !reload) {
2833 ReportWarning (_String("Already loaded ") & original_path.Enquote() & " from " & file_path);
2834 return true;
2835 }
2836
2837 if ((source_file = doFileOpen (file_path.get_str (), "rb")) == nil) {
2838 throw (_String("Could not read command file from '") &
2839 original_path & "' (expanded to '" & file_path & "')");
2840 }
2841 }
2842
2843 if (do_load_from_file && do_load_library && source_file) {
2844 ReportWarning (_String("Loaded ") & original_path.Enquote() & " from " & file_path.Enquote());
2845 loadedLibraryPaths.Insert (new _String (file_path),0,false,true);
2846 }
2847
2848 source_code = new _StringBuffer (source_file);
2849
2850 if (fclose (source_file) ) { // failed to fclose
2851 DeleteObject (source_code);
2852 throw (_String("Internal error: failed in a call to fclose ") & file_path.Enquote());
2853 }
2854 pop_path = true;
2855 PushFilePath (file_path);
2856 } else { // commands are not loaded from a file
2857 source_code = new _StringBuffer (_ProcessALiteralArgument(*GetIthParameter(0UL), current_program));
2858 }
2859
2860
2861 if (!source_code || source_code->empty()) {
2862 throw _String("Empty/missing source code string");
2863 }
2864
2865 if (!do_load_from_file && GetIthParameter(1UL)->nonempty()) {
2866 pop_path = true;
2867 PushFilePath (*GetIthParameter(1UL), false, false);
2868 }
2869
2870 _String * use_this_namespace = nil;
2871
2872 if (parameter_count() >= 3UL) { // stdin redirect (and/or name space prefix)
2873 _AssociativeList * input_arguments = nil;
2874 try {
2875 input_arguments = (_AssociativeList *)_ProcessAnArgumentByType(*GetIthParameter(2UL), ASSOCIATIVE_LIST, current_program, &dynamic_reference_manager);
2876 } catch (const _String& err) {
2877 if (parameter_count() == 3UL) {
2878 throw (err);
2879 }
2880 }
2881
2882
2883 if (input_arguments) {
2884
2885
2886
2887 _List *keys = input_arguments->GetKeys();
2888 dynamic_reference_manager < keys;
2889 keys->ForEach ([&] (BaseRef item, unsigned long) -> void {
2890 _String * key = (_String*) item;
2891 if (key) {
2892 HBLObjectRef payload = input_arguments->GetByKey (*key, STRING | ASSOCIATIVE_LIST);
2893 if (!payload) {
2894 throw ((_String("All entries in the associative array used as input redirect argument to ExecuteCommands/ExecuteAFile must be strings or associative lists (for keyword arguments). The following key was not: ") & key->Enquote()));
2895 }
2896 if (key->BeginsWith ("--") && key->length() > 2) {
2897 if (!user_kwargs) {
2898 dynamic_reference_manager < (user_kwargs = new _AssociativeList);
2899 }
2900 user_kwargs->MStore(key->Cut (2, kStringEnd), payload, true);
2901 has_user_kwargs = true;
2902 } else {
2903 argument_list.Insert (new _String (*key), (long)new _String (((_FString*)payload)->get_str()), false);
2904 if (payload->ObjectClass() != STRING) {
2905 throw ((_String("All entries in the associative array used as input redirect argument to ExecuteCommands/ExecuteAFile must be strings. The following key was not: ") & key->Enquote()));
2906 }
2907 has_redirected_input = true;
2908 }
2909 }
2910 });
2911
2912
2913
2914 if (parameter_count() > 3UL) {
2915 _String const namespace_for_code = _ProcessALiteralArgument (*GetIthParameter(3UL),current_program);
2916 if (namespace_for_code.nonempty()) {
2917 if (!namespace_for_code.IsValidIdentifier(fIDAllowCompound)) {
2918 throw (_String("Invalid namespace ID in call to ExecuteCommands/ExecuteAFile: ") & GetIthParameter(3UL)->Enquote());
2919 }
2920 dynamic_reference_manager < (use_this_namespace = new _String (namespace_for_code));
2921 }
2922 }
2923 }
2924 }
2925
2926 if (parameter_count () < 4UL && current_program.nameSpacePrefix) {
2927 dynamic_reference_manager < (use_this_namespace = new _String (*current_program.nameSpacePrefix->GetName()));
2928 }
2929
2930 if (source_code->BeginsWith ("#NEXUS")) {
2931 ReadDataSetFile (nil,1,source_code,nil,use_this_namespace);
2932 } else {
2933 bool result = false;
2934
2935 _ExecutionList code (*source_code, use_this_namespace, false, &result);
2936
2937 if (!result) {
2938 throw (_String("Encountered an error while parsing HBL"));
2939 } else {
2940
2941
2942 _AVLListXL * stash1 = nil;
2943 _List * stash2 = nil,
2944 * stash_kw_tags = nil;
2945
2946 _AssociativeList * stash_kw = nil;
2947
2948 bool update_kw = false;
2949
2950 if (has_redirected_input) {
2951 code.stdinRedirectAux = &_aux_argument_list;
2952 code.stdinRedirect = &argument_list;
2953 } else {
2954 if (current_program.has_stdin_redirect()) {
2955 stash1 = current_program.stdinRedirect;
2956 stash2 = current_program.stdinRedirectAux;
2957 current_program.stdinRedirect->AddAReference();
2958 current_program.stdinRedirectAux->AddAReference();
2959 }
2960 code.stdinRedirect = current_program.stdinRedirect;
2961 code.stdinRedirectAux = current_program.stdinRedirectAux;
2962 }
2963
2964
2965 bool ignore_ces_args = false;
2966
2967 if (has_user_kwargs) {
2968 code.SetKWArgs(user_kwargs);
2969 ignore_ces_args = true;
2970 } else {
2971 if (current_program.has_keyword_arguments()) {
2972 code.kwarg_tags = stash_kw_tags = current_program.kwarg_tags;
2973 code.kwargs = stash_kw = current_program.kwargs;
2974 if (stash_kw_tags) current_program.kwarg_tags->AddAReference();
2975 if (stash_kw) current_program.kwargs->AddAReference();
2976 code.currentKwarg = current_program.currentKwarg;
2977 update_kw = true;
2978 }
2979 }
2980
2981
2982
2983 if (!simpleParameters.empty() && code.TryToMakeSimple(true)) {
2984 ReportWarning (_String ("Successfully compiled an execution list (possibly partially).\n") & _String ((_String*)code.toStr()) );
2985 code.ExecuteSimple ();
2986 } else {
2987 code.Execute(nil, ignore_ces_args);
2988 }
2989
2990 if (stash1) {
2991 stash1->RemoveAReference();
2992 stash2->RemoveAReference();
2993 }
2994
2995 if (stash_kw_tags) stash_kw_tags->RemoveAReference();
2996 if (stash_kw) stash_kw->RemoveAReference();
2997
2998 code.stdinRedirectAux = nil;
2999 code.stdinRedirect = nil;
3000 if (update_kw) {
3001 code.kwarg_tags = nil;
3002 code.kwargs = nil;
3003 current_program.currentKwarg = code.currentKwarg;
3004 }
3005
3006 if (code.result) {
3007 DeleteObject (current_program.result);
3008 current_program.result = code.result;
3009 code.result = nil;
3010 }
3011 }
3012 }
3013 } catch (const _String& error) {
3014 cleanup ();
3015 return _DefaultExceptionHandler (nil, error, current_program);
3016 }
3017
3018 cleanup ();
3019 return true;
3020
3021 }
3022
3023 //____________________________________________________________________________________
3024
HandleDoSQL(_ExecutionList & current_program)3025 bool _ElementaryCommand::HandleDoSQL (_ExecutionList& current_program) {
3026
3027 static _SimpleList sql_databases;
3028
3029 static const _String kSQLOpen ("SQL_OPEN"),
3030 kSQLClose ("SQL_CLOSE");
3031
3032 static auto sql_handler = [] (void* exL,int cc, char** rd, char** cn) -> int
3033 {
3034 static const _String kSQLRowData ("SQL_ROW_DATA"),
3035 kSQLColumnNames ("SQL_COLUMN_NAMES");
3036
3037 _ExecutionList * caller = (_ExecutionList *)exL;
3038
3039 if (!terminate_execution)
3040 if (caller && cc && !caller->empty()) {
3041 _List row_data,
3042 column_names;
3043
3044 for (long cnt = 0; cnt < cc; cnt++) {
3045 if (rd[cnt]) {
3046 row_data.AppendNewInstance (new _String (rd[cnt]));
3047 } else {
3048 row_data.AppendNewInstance (new _String);
3049 }
3050
3051 if (cn[cnt]) {
3052 column_names.AppendNewInstance (new _String (cn[cnt]));
3053 } else {
3054 column_names.AppendNewInstance (new _String);
3055 }
3056 }
3057
3058 CheckReceptacleCommandIDException (&kSQLRowData,HY_HBL_COMMAND_DO_SQL,false)->SetValue (new _Matrix (row_data), false, true, NULL);
3059 CheckReceptacleCommandIDException (&kSQLColumnNames, HY_HBL_COMMAND_DO_SQL,false)->SetValue (new _Matrix (column_names), false,true, NULL);
3060
3061 caller->Execute();
3062
3063 }
3064 return 0;
3065 };
3066
3067
3068 _Variable * receptacle = nil;
3069
3070 current_program.advance ();
3071
3072 try {
3073 #ifdef __HYPHY_NO_SQLITE__
3074 throw _String("SQLite commands can not be used in a HyPhy version built with the __HYPHY_NO_SQLITE__ flag");
3075 #else
3076
3077 if (*GetIthParameter(0UL) == kSQLOpen) { // open a DB from file path
3078 receptacle = _ValidateStorageVariable (current_program, 2UL);
3079 _String file_name = _ProcessALiteralArgument (*GetIthParameter(1UL), current_program);
3080 if (!ProcessFileName(file_name, true, false, (hyPointer)current_program.nameSpacePrefix,false,¤t_program)) {
3081 return false;
3082 }
3083
3084 sqlite3 *db = nil;
3085 int err_code = sqlite3_open (file_name.get_str(),&db);
3086 if (err_code == SQLITE_OK) {
3087 err_code = sqlite3_exec(db, "SELECT COUNT(*) FROM sqlite_master", sql_handler, nil, nil);
3088 }
3089
3090 if (err_code != SQLITE_OK) {
3091 _String err_msg = sqlite3_errmsg(db);
3092 sqlite3_close(db);
3093 throw (err_msg);
3094 } else {
3095 long empty_slot = sql_databases.Find (0);
3096 if (empty_slot == kNotFound) {
3097 empty_slot = sql_databases.countitems();
3098 sql_databases << (long)db;
3099 } else {
3100 sql_databases.list_data[empty_slot] = (long)db;
3101 }
3102
3103 sqlite3_busy_timeout (db, 5000);
3104
3105 receptacle->SetValue (new _Constant (empty_slot), false,true, NULL);
3106 }
3107 } else {
3108 bool closing_db = *GetIthParameter(0UL) == kSQLClose;
3109 long db_index = _ProcessNumericArgumentWithExceptions(*GetIthParameter(closing_db ? 2UL : 0UL), current_program.nameSpacePrefix);
3110
3111 if (sql_databases.Map (db_index, 0L) == 0L) {
3112 throw (GetIthParameter(closing_db ? 2UL : 0UL)-> Enquote() & " is an invalid database index");
3113 }
3114
3115 if (closing_db) {
3116 sqlite3_close ((sqlite3*)sql_databases.get(db_index));
3117 sql_databases [db_index] = 0L;
3118 } else {
3119 _StringBuffer callback_code = _ProcessALiteralArgument(*GetIthParameter(2UL), current_program);
3120
3121 _ExecutionList callback_script (callback_code,current_program.nameSpacePrefix?(current_program.nameSpacePrefix->GetName()):nil);
3122
3123 if (!terminate_execution) {
3124 _String const sql_code = _ProcessALiteralArgument(*GetIthParameter(1UL), current_program);
3125
3126 if (sqlite3_exec((sqlite3*)sql_databases.get(db_index), sql_code.get_str(), sql_handler, (hyPointer)&callback_script, nil) != SQLITE_OK) {
3127 throw _String(sqlite3_errmsg((sqlite3*)sql_databases(db_index)));
3128 }
3129 }
3130 }
3131 }
3132 #endif
3133
3134 } catch (const _String& error) {
3135 return _DefaultExceptionHandler (receptacle, error, current_program);
3136 }
3137
3138 return true;
3139 }
3140
3141 //____________________________________________________________________________________
3142
HandleKeywordArgument(_ExecutionList & current_program)3143 bool _ElementaryCommand::HandleKeywordArgument (_ExecutionList& current_program) {
3144 current_program.advance ();
3145 #ifdef __HYPHYMPI__
3146 // ignore keyword options for MPI nodes
3147 if (hy_mpi_node_rank > 0L) {
3148 return true;
3149 }
3150 #endif
3151
3152 try {
3153 _String keyword = _ProcessALiteralArgument(*GetIthParameter(0UL), current_program),
3154 description = _ProcessALiteralArgument(*GetIthParameter(1UL), current_program),
3155 *default_value = nil;
3156
3157
3158
3159 _List reference_manager;
3160 if (parameter_count () > 2) {
3161 try {
3162 HBLObjectRef default_expression = _ProcessAnArgumentByType (*GetIthParameter(2UL), STRING|HY_UNDEFINED|NUMBER, current_program, nil);
3163 if (default_expression) {
3164 reference_manager < default_expression;
3165 if (default_expression->ObjectClass () == STRING) {
3166 default_value = new _String (((_FString*)default_expression)->get_str());
3167 reference_manager < default_value;
3168 } else {
3169 if (default_expression->ObjectClass () == NUMBER) {
3170 default_value = new _String (((_Constant*)default_expression)->Value());
3171 reference_manager < default_value;
3172 } else {
3173 // null values are handled separately
3174 }
3175 }
3176 }
3177 } catch (_String const& e){
3178 throw _String ("Not a valid type for the default expression value");
3179 }
3180 }
3181
3182 if (!current_program.kwarg_tags) {
3183 current_program.kwarg_tags = new _List;
3184 }
3185
3186 _List * new_tagged_keyword = new _List (new _String (keyword), new _String (description));
3187 if (default_value) {
3188 (*new_tagged_keyword) << default_value;
3189 } else {
3190 if (parameter_count () > 2) {
3191 (*new_tagged_keyword) < (_String*)nil;
3192 }
3193 }
3194
3195 if (parameter_count () > 3) {
3196 (*new_tagged_keyword) < new _String(_ProcessALiteralArgument(*GetIthParameter(3UL), current_program));
3197 }
3198
3199
3200 (*current_program.kwarg_tags) < new_tagged_keyword;
3201
3202 //printf ("%s\n", _String ((_String*)current_program.kwarg_tags->toStr()).get_str());
3203
3204 } catch (const _String& error) {
3205 return _DefaultExceptionHandler (nil, error, current_program);
3206 }
3207
3208 return true;
3209
3210 }
3211
3212 //____________________________________________________________________________________
3213
HandleGetString(_ExecutionList & current_program)3214 bool _ElementaryCommand::HandleGetString (_ExecutionList& current_program) {
3215
3216 auto make_fstring_pointer = [] (_String * s) -> _FString * {return new _FString (s);};
3217 auto make_fstring = [] (_String const s) -> _FString * {return new _FString (new _String (s));};
3218
3219 static const _String kVersionString ("HYPHY_VERSION"),
3220 kTimeStamp ("TIME_STAMP"),
3221 kListLoadedLibraries ("LIST_OF_LOADED_LIBRARIES");
3222
3223
3224 _Variable * receptacle = nil;
3225
3226 current_program.advance ();
3227
3228 try {
3229 receptacle = _ValidateStorageVariable (current_program);
3230
3231 long index1 = _ProcessNumericArgumentWithExceptions (*GetIthParameter(2),current_program.nameSpacePrefix),
3232 index2 = parameter_count() > 3 ? _ProcessNumericArgumentWithExceptions (*GetIthParameter(3),current_program.nameSpacePrefix) : -1L;
3233
3234
3235 // first, handle special, hardcoded cases
3236
3237 HBLObjectRef return_value = nil;
3238
3239 if (*GetIthParameter(1UL) == kVersionString) {
3240 if (index1 > 1.5) { // console header form
3241 return_value = make_fstring (_String ("HyPhy version ") & kHyPhyVersion);
3242 } else {
3243 if (index1 > 0.5) { // long form
3244 return_value = make_fstring(GetVersionString());
3245 } else {
3246 return_value = make_fstring (kHyPhyVersion);
3247 }
3248 }
3249 } else if (*GetIthParameter(1UL) == kTimeStamp) {
3250 return_value = make_fstring (GetTimeStamp (index1 < 0.5));
3251 } else if (*GetIthParameter(1UL) == kListLoadedLibraries) {
3252 return_value = new _Matrix (loadedLibraryPaths.Keys());
3253 }
3254
3255 if (!return_value) {
3256
3257 // next, handle cases like GetString (res, LikelihoodFunction, index)
3258
3259 long type_index = _HY_GetStringGlobalTypes.Find (GetIthParameter(1UL));
3260
3261 if (type_index != kNotFound) {
3262 type_index = _HY_GetStringGlobalTypes.GetXtra (type_index);
3263
3264
3265 if (type_index != HY_BL_TREE) {
3266 _String * object_name = (_String*)GetObjectNameByType (type_index, index1);
3267 if (!object_name) {
3268 throw (_String ("There is no ") & GetIthParameter(1UL)->Enquote() & " object with index " & index1);
3269 }
3270 if (type_index == HY_BL_HBL_FUNCTION) {
3271 return_value = &(
3272 (*new _AssociativeList)
3273 < (_associative_list_key_value){"ID", new _FString (*object_name) }
3274 < (_associative_list_key_value){"Arguments", new _Matrix(GetBFFunctionArgumentList(index1)) }
3275 );
3276 } else {
3277 return_value = make_fstring (*object_name);
3278 }
3279
3280 } else {
3281 _String * tree_name = FetchMathObjectNameOfTypeByIndex (TREE, index1);
3282 if (!tree_name) {
3283 throw (_String ("There is no ") & GetIthParameter(1UL)->Enquote() & " object with index " & index1);
3284 }
3285 return_value = make_fstring(*tree_name);
3286 }
3287
3288 receptacle->SetValue (return_value, false,true, NULL);
3289 return true;
3290 }
3291
3292 // next, handle lookup of HBL objects
3293
3294 const _String source_name = AppendContainerName (*GetIthParameter(1), current_program.nameSpacePrefix);
3295 long object_type = HY_BL_ANY,
3296 object_index;
3297
3298 BaseRefConst source_object = nil;
3299 try {
3300 source_object = _GetHBLObjectByTypeMutable (source_name, object_type, &object_index);
3301 }
3302 catch (_String const & err) { // lookup failed
3303 return_value = nil;
3304 object_type = HY_BL_NOT_DEFINED;
3305 }
3306
3307 switch (object_type) {
3308 case HY_BL_DATASET: {
3309 _DataSet const* data_set_object = (_DataSet const*)source_object;
3310 if (index1 >=0) { // get a specific sequence name
3311 if (index1 < data_set_object->NoOfSpecies()) {
3312 return_value = make_fstring (*(_String*)(data_set_object->GetNames().GetItem(index1)));
3313 } else {
3314 throw (_String (index1) & " exceeds the maximum index for the underlying DataSet object");
3315 }
3316 } else {
3317 return_value = new _Matrix (data_set_object->GetNames(), false);
3318 }
3319 break;
3320 }
3321 case HY_BL_DATASET_FILTER: {
3322 _DataSetFilter const* data_filter = (_DataSetFilter const*)source_object;
3323 ReleaseDataFilterLock (object_index);
3324 if (index1 >=0 ) { // get a specific sequence name
3325 if (index1 < data_filter->NumberSpecies()) {
3326 return_value = make_fstring (*(_String*)data_filter->GetData()->GetNames().GetItem(data_filter->theNodeMap.Element(index1)));
3327 } else {
3328 throw (_String (index1) & " exceeds the maximum index for the underlying DataSetFilter object");
3329 }
3330 } else {
3331 _List filter_seq_names;
3332 _List const * original_names = &data_filter->GetData()->GetNames();
3333 data_filter->theNodeMap.Each ([&] (long value, unsigned long) -> void {
3334 filter_seq_names << original_names->GetItem (value);
3335 });
3336 return_value = new _Matrix (filter_seq_names);
3337
3338 }
3339 break;
3340 }
3341
3342 case HY_BL_HBL_FUNCTION: {
3343 return_value = &(
3344 (*new _AssociativeList)
3345 < (_associative_list_key_value){"ID", new _FString (*GetObjectNameByType (HY_BL_HBL_FUNCTION, object_index, false)) }
3346 < (_associative_list_key_value){"Arguments", new _Matrix(GetBFFunctionArgumentList(object_index)) }
3347 < (_associative_list_key_value){"Body", new _FString (GetBFFunctionBody(object_index).sourceText,false) }
3348 );
3349 break;
3350 }
3351
3352 case HY_BL_LIKELIHOOD_FUNCTION:
3353 case HY_BL_SCFG: {
3354 _LikelihoodFunction const *lf = (_LikelihoodFunction const*) source_object;
3355 if (index1 >= 0L) {
3356 if (index1<lf->GetIndependentVars().countitems()) {
3357 return_value = make_fstring (*LocateVar(lf->GetIndependentVars().GetElement(index1))->GetName());
3358 } else {
3359 if (index1 < lf->GetIndependentVars().countitems()+lf->GetDependentVars().countitems()) {
3360 return_value = make_fstring (*LocateVar(lf->GetDependentVars().GetElement(index1-lf->GetIndependentVars().countitems()))->GetName());
3361 } else {
3362 throw (_String (index1) & " exceeds the maximum index for the underlying LikelihoodFunction/SCFG object");
3363 }
3364 }
3365 } else {
3366 return_value = lf->CollectLFAttributes ();
3367 if (object_type == HY_BL_SCFG) {
3368 ((Scfg* const)lf)->AddSCFGInfo ((_AssociativeList*)return_value);
3369 }
3370 }
3371 break;
3372 }
3373
3374 case HY_BL_MODEL: {
3375 if (index1 >= 0L) {
3376 if (index2 < 0L) { // get the name of index1 parameter
3377 long variable_index = PopulateAndSort ([&] (_AVLList & parameter_list) -> void {
3378 ScanModelForVariables (object_index, parameter_list, false, -1, false);
3379 }).Map (index1);
3380 if (variable_index >= 0) {
3381 return_value = make_fstring (*LocateVar(variable_index)->GetName());
3382 } else {
3383 throw (_String (index1) & " exceeds the maximum parameter index for the underlying Model object");
3384 }
3385 } else { // get the formula for cell (index1, index2)
3386 if (!IsModelOfExplicitForm (object_index)) {
3387 _Variable* rate_matrix = (_Variable*)source_object;
3388 _Formula * cell = ((_Matrix*)rate_matrix->GetValue())->GetFormula (index1,index2);
3389 if (cell) {
3390 return_value = make_fstring_pointer ((_String*)cell->toStr(kFormulaStringConversionNormal));
3391 } else {
3392 throw (_String("Invalid rate matrix cell index"));
3393 }
3394 } else {
3395 throw (_String("Direct indexing of rate matrix cells is not supported for expression based (e.g. mixture) substitution models"));
3396 }
3397 }
3398
3399 } else {
3400 _Variable * rates, * freqs;
3401 bool is_canonical;
3402 RetrieveModelComponents (object_index, rates, freqs, is_canonical);
3403
3404 if (rates) {
3405 if (index1 == -1) { // branch length expression
3406 return_value = make_fstring_pointer (((_Matrix*)rates->GetValue())->BranchLengthExpression((_Matrix*)freqs->GetValue(),is_canonical));
3407 } else
3408 /*
3409 returns an AVL with keys
3410 "RATE_MATRIX" - the ID of the rate matrix
3411 "EQ_FREQS" - the ID of eq. freq. vector
3412 "MULT_BY_FREQ" - a 0/1 flag to determine which format the matrix is in.
3413 */
3414 {
3415 return_value = &(
3416 (*new _AssociativeList)
3417 < (_associative_list_key_value){"RATE_MATRIX",new _FString(*rates->GetName())}
3418 < (_associative_list_key_value){"EQ_FREQS",new _FString(*freqs->GetName()) }
3419 < (_associative_list_key_value){"MULT_BY_FREQ",new _Constant (is_canonical) }
3420 );
3421 }
3422 } else {
3423 throw _String("Failed to retrieve model rate matrix");
3424 }
3425 }
3426 break;
3427 }
3428 case HY_BL_BGM: {
3429 //ReportWarning(_String("In HandleGetString() for case HY_BL_BGM"));
3430 _BayesianGraphicalModel * this_bgm = (_BayesianGraphicalModel *) source_object;
3431
3432 switch (index1) {
3433 case HY_HBL_GET_STRING_BGM_SCORE: { // return associative list containing node score cache
3434 _AssociativeList * export_alist = new _AssociativeList;
3435
3436 if (this_bgm -> ExportCache (export_alist)) {
3437 return_value = export_alist;
3438 } else {
3439 DeleteObject (export_alist);
3440 throw _String("Failed to export node score cache for BGM");
3441 }
3442
3443 break;
3444 }
3445 case HY_HBL_GET_STRING_BGM_SERIALIZE: { // return associative list with network structure and parameters
3446 _StringBuffer *serialized_bgm = new _StringBuffer (1024L);
3447 this_bgm -> SerializeBGM (*serialized_bgm);
3448 return_value = new _FString (serialized_bgm);
3449 break;
3450 }
3451 default: {
3452 throw _String ("Unrecognized index ") & index1 & " for a BGM object";
3453 }
3454 }
3455 break;
3456 }
3457 }
3458 }
3459
3460 // finally, try to look up a variable
3461 if (!return_value) {
3462 _Variable* var = FetchVar(LocateVarByName (AppendContainerName(*GetIthParameter(1UL), current_program.nameSpacePrefix)));
3463
3464 if (var) {
3465 if (var->IsIndependent() && index1 != -3) {
3466 if (!var->has_been_set()) {
3467 return_value = new _MathObject;
3468 } else {
3469 return_value = make_fstring_pointer ((_String*) var->toStr());
3470 }
3471 } else {
3472 if (index1 == -1){
3473 _SimpleList variable_list = PopulateAndSort ([&] (_AVLList & parameter_list) -> void {
3474 var->ScanForVariables (parameter_list, true);
3475 });
3476 _AssociativeList * var_list_by_kind = new _AssociativeList;
3477
3478 _List split_vars;
3479 SplitVariableIDsIntoLocalAndGlobal (variable_list, split_vars);
3480 InsertVarIDsInList (var_list_by_kind, "Global", *(_SimpleList*)split_vars(0));
3481 InsertVarIDsInList (var_list_by_kind, "Local", *(_SimpleList*)split_vars(1));
3482 return_value = var_list_by_kind;
3483 }
3484 else { // formula string
3485
3486 if (index1 == -3 || index1 == -4) {
3487 _StringBuffer local, global;
3488
3489 _SimpleList var_index;
3490 if (index1 == -3 || var->IsIndependent()) {
3491 var_index << var->get_index ();
3492 if (var->IsIndependent()) {
3493 //printf ("ExportIndVariables\n");
3494 ExportIndVariables (global, local, &var_index);
3495 } else {
3496 //printf ("ExportDepVariables\n");
3497 ExportDepVariables(global, local, &var_index);
3498 }
3499 } else {
3500 _AVLList vl (&var_index);
3501 var->ScanForVariables(vl, true);
3502 _SimpleList ind_vars = var_index.Filter([] (long index, unsigned long) -> bool {return LocateVar(index)->IsIndependent();}),
3503 dep_vars = var_index.Filter([] (long index, unsigned long) -> bool {return !LocateVar(index)->IsIndependent();});
3504 ExportIndVariables(global, local, &ind_vars);
3505 ExportDepVariables(global, local, &dep_vars);
3506 }
3507 return_value = make_fstring_pointer ( & ((*new _StringBuffer (128L)) << global << local << '\n'));
3508
3509 } else {
3510 _Matrix * formula_matrix = (index2 >= 0 && var->ObjectClass() == MATRIX) ? (_Matrix*)var->GetValue () : nil;
3511 if (formula_matrix) {
3512 _Formula* cell = formula_matrix->GetFormula(index1, index2);
3513 if (cell) {
3514 return_value = make_fstring_pointer ((_String*) cell->toStr(kFormulaStringConversionNormal));
3515 }
3516 } else {
3517 return_value = make_fstring_pointer ((_String*)var->GetFormulaString (kFormulaStringConversionNormal));
3518 }
3519 }
3520 }
3521 }
3522 }
3523 }
3524
3525 if (!return_value) {
3526 throw (_String("No viable object to obtain information from"));
3527 }
3528
3529 receptacle->SetValue (return_value, false,true, NULL);
3530
3531
3532 }
3533
3534 catch (const _String& error) {
3535 return _DefaultExceptionHandler (receptacle, error, current_program);
3536 }
3537
3538 return true;
3539
3540 }
3541
3542
3543 //____________________________________________________________________________________
3544
HandleFscanf(_ExecutionList & current_program,bool is_sscanf)3545 bool _ElementaryCommand::HandleFscanf (_ExecutionList& current_program, bool is_sscanf) {
3546
3547 const static _String kFscanfStdin ("stdin");
3548 static long last_call_stream_position = 0L;
3549
3550 current_program.advance();
3551
3552 _List dynamic_reference_manager;
3553
3554 try {
3555
3556 long current_stream_position = 0L,
3557 started_here_position = 0L;
3558
3559 unsigned long has_rewind = simpleParameters.Element(-1L) == -1 ? 1UL : 0UL,
3560 has_create = simpleParameters.Element(-1L) == -2 ? 1UL : 0UL;
3561
3562 _String const * input_data = nil;
3563
3564 _String source_name = *GetIthParameter(0UL);
3565
3566 unsigned long argument_index = 0UL,
3567 upper_bound = (has_rewind || has_create)? simpleParameters.countitems() - 1L : simpleParameters.countitems();
3568
3569 if (source_name == kFscanfStdin) {
3570 bool need_to_ask_user = true;
3571 if (current_program.has_stdin_redirect () || current_program.has_keyword_arguments()) {
3572 try {
3573 _FString * redirect = (_FString*)hy_env::EnvVariableGet(hy_env::fprintf_redirect, STRING);
3574 _String * redirected = current_program.FetchFromStdinRedirect(nil, false, !(redirect && redirect->has_data()));
3575 input_data = redirected;
3576 // echo the input if there is no fprintf redirect in effect
3577
3578
3579 dynamic_reference_manager < redirected;
3580 need_to_ask_user = false;
3581 } catch (const _String& e) {
3582 if (e != kNoKWMatch) {
3583 throw (e);
3584 }
3585 }
3586 }
3587
3588 if (need_to_ask_user){ // read from stdin
3589 if (hy_env::EnvVariableTrue(hy_env::end_of_file) == false && source_name == hy_scanf_last_file_path) {
3590 throw _String("Ran out of standard input");
3591 }
3592 _String * console_data = new _String (StringFromConsole());
3593 dynamic_reference_manager < console_data;
3594 input_data = console_data;
3595 }
3596 } else { // not stdin
3597
3598 if (is_sscanf) {
3599 _FString * source_string = (_FString*)_ProcessAnArgumentByType(source_name,STRING,current_program,&dynamic_reference_manager);
3600 input_data = & source_string->get_str();
3601 if (hy_env::EnvVariableTrue(hy_env::end_of_file)) { // reset path
3602 hy_scanf_last_file_path = kEmptyString;
3603 }
3604
3605 if (source_name != hy_scanf_last_file_path || has_rewind || has_create) { // new stream, or rewind on the current stream
3606 hy_scanf_last_file_path = source_name;
3607 current_stream_position = 0L;
3608 last_call_stream_position = 0L;
3609 } else {
3610 current_stream_position = last_call_stream_position;
3611 started_here_position = last_call_stream_position;
3612 if (last_call_stream_position >= input_data->length()) {
3613 // run out of input chars, set EOF and bail
3614 hy_env::EnvVariableSet(hy_env::end_of_file, new HY_CONSTANT_TRUE, false);
3615 return true;
3616 }
3617 }
3618 } else {
3619 _String file_path;
3620
3621
3622
3623 if (file_path == hy_env::kDevNull) {
3624 return true;
3625 }
3626
3627 if (source_name == kPromptForFilePlaceholder) {
3628 file_path = source_name;
3629 } else {
3630 file_path = ((_FString*)_ProcessAnArgumentByType(source_name,STRING,current_program,&dynamic_reference_manager))->get_str();
3631 }
3632
3633
3634 if (!ProcessFileName(file_path, true,false,(hyPointer)current_program.nameSpacePrefix, false, ¤t_program)) {
3635 return false;
3636 }
3637
3638 FILE * input_stream = doFileOpen (file_path.get_str(), "rb");
3639 hy_env::EnvVariableSet(hy_env::file_created, new HY_CONSTANT_FALSE, false);
3640 if (!input_stream) {
3641 if (has_create) {
3642 input_stream = doFileOpen (file_path.get_str(), "wb");
3643 if (input_stream){
3644 while (argument_index < upper_bound) {
3645 _Variable * store_here = _ValidateStorageVariable (current_program, argument_index + 1UL);
3646 store_here->SetValue(new _MathObject, false, true, nil);
3647 argument_index++;
3648 }
3649 hy_env::EnvVariableSet(hy_env::end_of_file, new HY_CONSTANT_TRUE, false);
3650 hy_env::EnvVariableSet(hy_env::file_created, new HY_CONSTANT_TRUE, false);
3651 return true;
3652 } else {
3653 throw (file_path.Enquote() & " could not be opened for reading or writing (CREATE mode) by fscanf. Path stack:\n\t" & GetPathStack("\n\t"));
3654 }
3655 }
3656
3657 throw (file_path.Enquote() & " could not be opened for reading by fscanf. Path stack:\n\t" & GetPathStack("\n\t"));
3658 }
3659 if (hy_env::EnvVariableTrue(hy_env::end_of_file)) { // reset path
3660 hy_scanf_last_file_path = kEmptyString;
3661 }
3662 if (source_name != hy_scanf_last_file_path || has_rewind || has_create) { // new stream, or rewind on the current stream
3663 hy_scanf_last_file_path = source_name;
3664 last_call_stream_position = 0L;
3665 }
3666
3667 fseek (input_stream,0,SEEK_END);
3668 current_stream_position = ftell (input_stream);
3669 current_stream_position -= last_call_stream_position;
3670
3671 if (current_stream_position<=0) {
3672 hy_env::EnvVariableSet(hy_env::end_of_file, new HY_CONSTANT_TRUE, false);
3673 fclose(input_stream);
3674 return true;
3675 }
3676
3677 rewind (input_stream);
3678 fseek (input_stream, last_call_stream_position, SEEK_SET);
3679 _String * file_data = new _String (input_stream, current_stream_position);
3680 fclose (input_stream);
3681 dynamic_reference_manager < file_data;
3682 input_data = file_data;
3683 current_stream_position = 0L;
3684 }
3685 }
3686
3687
3688 while (argument_index < upper_bound && current_stream_position < input_data->length()) {
3689 _Variable * store_here = _ValidateStorageVariable (current_program, argument_index + 1UL);
3690
3691 long lookahead = current_stream_position;
3692
3693 switch (simpleParameters.get(argument_index)) {
3694
3695 case 0L: { // number
3696 _SimpleList numerical_match (input_data->RegExpMatch(hy_float_regex, lookahead));
3697 if (numerical_match.empty()) {
3698 throw (_String("Failed to read a number from the input stream") & input_data->Cut (lookahead, kStringEnd));
3699 break;
3700 }
3701
3702 store_here->SetValue (new _Constant (input_data->Cut (numerical_match(0), numerical_match(1)).to_float ()), false,true, NULL);
3703 lookahead = input_data->FirstNonSpaceIndex(numerical_match (1) + 1, kStringEnd) ;
3704 }
3705 break;
3706
3707 case 3L: { // string
3708 lookahead = 0L;
3709 bool start_found=false;
3710 while (current_stream_position + lookahead < input_data->length ()) {
3711 char c = input_data->char_at (current_stream_position + lookahead);
3712 if (!start_found) {
3713 if (!isspace(c)) {
3714 current_stream_position += lookahead;
3715 start_found = true;
3716 lookahead = 0L;
3717 }
3718 } else if (c=='\n' || c=='\r' || c=='\t') {
3719 break;
3720 }
3721 lookahead++;
3722 }
3723
3724 if (start_found) {
3725 store_here->SetValue (new _FString (new _String(*input_data,current_stream_position,current_stream_position+lookahead-1)),false,true, NULL);
3726 } else {
3727 store_here->SetValue (new _FString, false,true, NULL);
3728 }
3729 lookahead = current_stream_position + lookahead - 1L;
3730 }
3731 break;
3732
3733 case 5L: { // Raw
3734 store_here->SetValue (new _FString (new _String (*input_data,current_stream_position,kStringEnd)), false,true, NULL);
3735 lookahead = input_data->length();
3736 }
3737 break;
3738
3739 case 6L: { // Lines
3740
3741
3742
3743 _String line_block (*input_data,current_stream_position,kStringEnd);
3744
3745 // break the line block by any of the three platform line breaks
3746 // \r, \n or \r\n
3747 _List lines;
3748 long last_break = 0L;
3749
3750 auto add_buffer = [&] (long s, long e) -> void {
3751 if (e > s) {
3752 lines.AppendNewInstance(new _String (line_block, s, e-1));
3753 } else {
3754 lines.AppendNewInstance(new _String (kEmptyString));
3755 }
3756 };
3757
3758 for (long i = 0UL; i < line_block.length(); i++) {
3759 char current_char = line_block.char_at(i);
3760 if (current_char == '\n') {
3761 add_buffer (last_break, i);
3762 last_break = i + 1L;
3763 } else {
3764 if (current_char == '\r') {
3765 add_buffer (last_break, i);
3766 if (line_block.char_at(i + 1L) == '\n') {
3767 i ++;
3768 }
3769 last_break = i + 1L;
3770 }
3771 }
3772 }
3773
3774 if (last_break < line_block.length()) {
3775 add_buffer (last_break, line_block.length ());
3776 }
3777
3778 store_here->SetValue (new _Matrix (lines, false), false,true,NULL);
3779 lookahead = input_data->length();
3780
3781 }
3782
3783 break;
3784
3785 default: {
3786 // TODO: 20170623 SLKP CHANGE: use expression extractor
3787
3788 lookahead = input_data->ExtractEnclosedExpression(current_stream_position, (simpleParameters.list_data[argument_index]==2)?'(':'{', (simpleParameters.list_data[argument_index]==2)?')':'}', fExtractRespectQuote | fExtractRespectEscape);
3789
3790 if (lookahead == kNotFound) {
3791 lookahead = input_data->length ();
3792 break;
3793 }
3794
3795 _String object_data (*input_data, current_stream_position, lookahead);
3796
3797 if (simpleParameters.list_data[argument_index] != 2) { // matrix
3798 _FormulaParsingContext def;
3799 store_here->SetValue (new _Matrix (object_data,simpleParameters.list_data[argument_index]==4, def), false,true,NULL);
3800 } else {
3801 _TheTree (*store_here->GetName(), object_data);
3802 }
3803
3804 }
3805 break;
3806
3807 } // end switch
3808
3809 current_stream_position = lookahead + 1L;
3810 argument_index ++;
3811 }
3812 if (argument_index + has_rewind + has_create <simpleParameters.countitems()) {
3813 hy_env::EnvVariableSet(hy_env::end_of_file, new HY_CONSTANT_TRUE, false);
3814 throw (_String("Could not read all the parameters requested."));
3815 } else {
3816 hy_env::EnvVariableSet(hy_env::end_of_file, new HY_CONSTANT_FALSE, false);
3817 }
3818
3819 last_call_stream_position += current_stream_position - started_here_position;
3820
3821 } catch (const _String& error) {
3822 if (hy_env::EnvVariableTrue(hy_env::soft_fileio_exceptions)) {
3823 hy_env::EnvVariableSet(hy_env::last_fileio_exception, new _FString (error,false), false);
3824 return true;
3825 }
3826 return _DefaultExceptionHandler (nil, error, current_program);
3827 }
3828
3829 return true;
3830
3831 }
3832
3833 //____________________________________________________________________________________
3834
HandleChoiceList(_ExecutionList & current_program)3835 bool _ElementaryCommand::HandleChoiceList (_ExecutionList& current_program) {
3836
3837 auto handle_exclusions = [] (long count, _SimpleList & excluded) -> const _SimpleList {
3838 _SimpleList lfids;
3839 lfids.Subtract(_SimpleList(count, 0L, 1L), excluded);
3840 return lfids;
3841 };
3842
3843 static const _String kSkipNone ("SKIP_NONE"),
3844 kNoSkip ("NO_SKIP"),
3845 kLikelihoodFunctions ("LikelihoodFunction");
3846
3847 static const unsigned long maximum_wrong_choices = 10UL;
3848
3849 current_program.advance();
3850
3851 _Variable * receptacle = nil;
3852
3853 _List local_dynamic_manager;
3854
3855
3856 try {
3857
3858 receptacle = _ValidateStorageVariable (current_program);
3859
3860 long number_of_choices = _ProcessNumericArgumentWithExceptions (*GetIthParameter(2UL),current_program.nameSpacePrefix);
3861 _String dialog_title = _ProcessALiteralArgument(*GetIthParameter(1UL), current_program),
3862 exclusions = *GetIthParameter(3UL);
3863
3864 _SimpleList selections,
3865 excluded;
3866
3867 bool variable_number = number_of_choices <= 0L,
3868 do_markdown = hy_env :: EnvVariableTrue(hy_env :: produce_markdown_output);
3869
3870
3871 if (exclusions != kSkipNone && exclusions != kNoSkip) {
3872 try {
3873 HBLObjectRef exlcusion_argument = _ProcessAnArgumentByType(exclusions, NUMBER | MATRIX, current_program, &local_dynamic_manager);
3874 if (exlcusion_argument->ObjectClass() == NUMBER) {
3875 excluded << exlcusion_argument->Compute ()->Value();
3876 } else {
3877 ((_Matrix*)exlcusion_argument)->ConvertToSimpleList (excluded);
3878 excluded.Sort();
3879 }
3880 } catch (_String const & e) {
3881 //printf ("%s\n", e.get_str());
3882 // no exclusions, so do nothing
3883 }
3884 }
3885
3886 _List * available_choices = nil;
3887
3888 if (simpleParameters.Element(0UL)) {
3889 // dynamically generated list of options
3890 _String const choices_parameter = *GetIthParameter(4UL);
3891 local_dynamic_manager < (available_choices = new _List);
3892
3893 if (choices_parameter == kLikelihoodFunctions) {
3894 // the list consists of all defined likelihood function objects
3895
3896 handle_exclusions (likeFuncNamesList.countitems(), excluded).Each([&] (long value, unsigned long) -> void {
3897 if (likeFuncList.GetItem(value)) {
3898 _String const * lf_name = (_String*) likeFuncNamesList (value);
3899 (*available_choices) < new _List (new _String (*lf_name), new _String ( _String ("Likelihood Function ") & *lf_name & "."));
3900 }
3901 });
3902 } else {
3903 // see if the argument is a reference to one of the standard HBL objects
3904 const _String source_name = AppendContainerName (choices_parameter, current_program.nameSpacePrefix);
3905
3906 long object_type = HY_BL_DATASET_FILTER | HY_BL_DATASET | HY_BL_MODEL,
3907 object_index;
3908
3909
3910 try {
3911 BaseRefConst source_object = _GetHBLObjectByType (source_name, object_type, &object_index, ¤t_program);
3912 // this wil also handle USE_LAST_MODEL
3913 switch (object_type) {
3914 case HY_BL_DATASET: {
3915 _DataSet const *ds = (_DataSet const*) source_object;
3916 handle_exclusions (ds->NoOfSpecies(), excluded).Each (
3917 [&] (long value, unsigned long) -> void {
3918 _String const * sequence_name = ds->GetSequenceName(value);
3919 (*available_choices) < new _List (new _String (*sequence_name), new _String ( _String ("Taxon ") & (value + 1L) & sequence_name->Enquote('(', ')') & "."));
3920 }
3921 );
3922 break;
3923 }
3924 case HY_BL_DATASET_FILTER: {
3925 _DataSetFilter const *df = (_DataSetFilter const*) source_object;
3926 handle_exclusions (df->NumberSpecies(), excluded).Each (
3927 [&] (long value, unsigned long) -> void {
3928 _String const * sequence_name = df->GetSequenceName(value);
3929 (*available_choices) < new _List (new _String (*sequence_name), new _String ( _String ("Taxon ") & (value + 1L) & sequence_name->Enquote('(', ')') & "."));
3930 }
3931 );
3932 break;
3933 }
3934 case HY_BL_MODEL: {
3935 (*available_choices) < new _List (new _String("All Parameters"), new _String("All local model parameters are constrained"));
3936
3937 _SimpleList model_indices = PopulateAndSort ([&] (_AVLList & parameter_list) -> void {
3938 if (IsModelOfExplicitForm (object_index)) {
3939 ((_Formula const*)source_object)->ScanFForVariables(parameter_list,false);
3940 } else {
3941 ((_Variable const*)source_object)->ScanForVariables(parameter_list,false);
3942 }
3943 });
3944 handle_exclusions (model_indices.countitems(), excluded).Each (
3945 [&] (long value, unsigned long) -> void {
3946 _String const * parameter_name = LocateVar(model_indices.get(value))->GetName();
3947 (*available_choices) < new _List (new _String (*parameter_name),
3948 new _String (_String ("Constrain parameter ") & *parameter_name & '.'));
3949 }
3950 );
3951 }
3952 break;
3953 }
3954 } catch (_String const & e) {
3955 // not an object
3956 try {
3957 _Matrix * target_variable = (_Matrix*)_ProcessAnArgumentByType(choices_parameter, MATRIX, current_program, &local_dynamic_manager);
3958 if (!target_variable->IsAStringMatrix()) {
3959 throw (choices_parameter.Enquote() & " is not a matrix of strings");
3960 }
3961 if (target_variable->GetVDim() != 2) {
3962 throw (choices_parameter.Enquote() & " is not a matrix with two columns");
3963 }
3964
3965 _List choices;
3966 target_variable->FillInList(choices, false);
3967 //choices.bumpNInst();
3968
3969 handle_exclusions (target_variable->GetHDim(), excluded).Each (
3970 [&] (long value, unsigned long ) -> void {
3971 _String const * parameter_name = LocateVar(value)->GetName();
3972 BaseRef key = choices.GetItem(value << 1),
3973 description = choices.GetItem(1L + (value << 1));
3974 key->AddAReference(); description->AddAReference();
3975 (*available_choices) < new _List (key,description);
3976 }
3977 );
3978
3979
3980 /*for (unsigned long k = 0UL; k < choices.countitems(); k+=2) {
3981 (*available_choices) < new _List (choices.GetItem(k), choices.GetItem(k+1));
3982 }*/
3983 } catch (_String const& e2) {
3984 throw (choices_parameter.Enquote() & " is not a supported object/literal for the list of choices");
3985 }
3986
3987 }
3988
3989 }
3990
3991 } else {
3992 available_choices = (_List*)parameters.GetItem(4UL);
3993 }
3994
3995 if (available_choices->empty()) {
3996 throw (_String("The list of choices is empty"));
3997 }
3998
3999 if ((long) available_choices->countitems() < number_of_choices) {
4000 throw (_String ("The list of choices has " ) & (long) available_choices->countitems() & " elements, but " & number_of_choices & " choices are required");
4001 }
4002
4003 auto validate_choice = [&] (_String const& user_choice) -> long {
4004 return available_choices->FindOnCondition([&] (BaseRefConst item, unsigned long) -> bool {
4005 return user_choice == *(_String*)((_List*)item)->GetItem (0);
4006 });
4007 };
4008
4009 long required = variable_number ? available_choices->countitems() : number_of_choices;
4010 bool need_to_prompt_user = true;
4011
4012 if (current_program.has_stdin_redirect() || current_program.has_keyword_arguments()) {
4013
4014 auto report_key_error = [&] (const _String & bad_key) -> void {
4015 _String choice_list_echo = available_choices->Map ([] (BaseRef object, unsigned long index) -> BaseRef {
4016 BaseRef title = ((_List*)object)->GetItem (0);
4017 title->AddAReference();
4018 return title;
4019 }, 0, MIN (32, available_choices->countitems())).Join(", ");
4020 if (available_choices->countitems() > 32) {
4021 choice_list_echo = choice_list_echo & " and " & _String(available_choices->countitems()-32) & " additional options";
4022 }
4023
4024 throw (bad_key.Enquote() & " is not a valid choice passed to '" & dialog_title & "' ChoiceList using redirected stdin input or keyword arguments. Valid choices are\n\t") & choice_list_echo & "\n";
4025 };
4026
4027
4028 while (selections.countitems() < required) {
4029 _String user_choice;
4030 try {
4031 _FString * redirect = (_FString*)hy_env::EnvVariableGet(hy_env::fprintf_redirect, STRING);
4032 user_choice = current_program.FetchFromStdinRedirect(&dialog_title, required > 1, !(redirect && redirect->has_data())); // handle multiple selections
4033 } catch (const _String& e) {
4034 if (e == kNoKWMatch) {
4035 break;
4036 } else {
4037 throw (e);
4038 }
4039 } catch (HBLObjectRef multiple_choice) {
4040 if (multiple_choice->ObjectClass() == ASSOCIATIVE_LIST) {
4041 _List option_list;
4042 ((_AssociativeList*)multiple_choice)->FillInList (option_list);
4043 option_list.ForEach([&] (BaseRef item, unsigned long) -> void {
4044 long match_found = validate_choice (_String ((_String*)item->toStr()));
4045 if (match_found == kNotFound) {
4046 report_key_error (_String ((_String*)item->toStr()));
4047 } else {
4048 selections << match_found;
4049 }
4050 });
4051 DeleteObject(multiple_choice);
4052 }
4053 }
4054
4055 if (variable_number && user_choice.empty()) {
4056 break;
4057 }
4058
4059 long match_found = validate_choice (user_choice);
4060 if (match_found == kNotFound) {
4061 report_key_error (user_choice);
4062 }
4063 selections << match_found;
4064 }
4065
4066 need_to_prompt_user = selections.countitems() != required && !variable_number || variable_number && selections.empty();
4067 }
4068
4069 if (need_to_prompt_user) {
4070 #ifdef __HEADLESS__
4071 throw (_String("Unhandled request for data from standard input (headless HyPhy)"));
4072 #endif
4073 if (do_markdown) {
4074 printf ("\n\n####%s\n", dialog_title.get_str());
4075 } else {
4076 const _String spacer (_String('-'), dialog_title.length());
4077 printf ("\n\n\t\t\t+%s+\n\t\t\t|%s|\n\t\t\t+%s+\n\n",
4078 spacer.get_str(),
4079 dialog_title.get_str(),
4080 spacer.get_str());
4081 }
4082
4083 unsigned long wrong_selections = 0UL;
4084
4085 _AVLList selection_tree (&selections);
4086
4087
4088 while (selections.countitems() < required) {
4089
4090 for (unsigned long option = 0UL; option < available_choices->countitems(); option ++) {
4091 if (selection_tree.FindLong (option) == kNotFound) {
4092 printf (do_markdown ? "\n%ld. [**%s**] %s" : "\n\t(%ld):[%s] %s", option+1UL, ((_String*)available_choices->GetItem (option, 0))->get_str(), ((_String*)available_choices->GetItem (option, 1))->get_str());
4093 }
4094
4095 }
4096
4097 if (variable_number) {
4098 printf ("\n\n%sPlease choose option %ld, enter d to complete selection, enter q to cancel:", do_markdown ? ">" : "",selections.countitems() + 1UL);
4099 } else {
4100 if (number_of_choices > 1L) {
4101 printf ("\n\n%sPlease choose option %ld of %ld (or press q to cancel selection):",do_markdown ? ">" : "", selections.countitems() + 1UL ,number_of_choices);
4102 } else {
4103 printf ("\n\n%sPlease choose an option (or press q to cancel selection):", do_markdown ? ">" : "");
4104 }
4105 }
4106
4107 _String user_choice (StringFromConsole());
4108 if (user_choice.length () == 1UL && toupper(user_choice.char_at(0L)) == 'Q') {
4109 selections.Clear();
4110 selections << -1L;
4111 break;
4112 } else {
4113 if (variable_number && user_choice.length () == 1UL && toupper(user_choice.char_at(0L)) == 'D') {
4114 if (selections.empty()) {
4115 selections << -1L;
4116 }
4117 break;
4118 }
4119 }
4120
4121 long integer_choice = user_choice.to_long();
4122
4123 if (integer_choice > 0L && integer_choice <= available_choices->countitems()) {
4124 selection_tree.Insert((BaseRef)(integer_choice-1));
4125 } else {
4126 ++ wrong_selections;
4127 }
4128
4129 if (wrong_selections > maximum_wrong_choices) {
4130 throw _String ("Too many invalid imputs");
4131 }
4132 }
4133
4134 selection_tree.ReorderList();
4135
4136
4137 }
4138
4139 if (selections.empty () || selections.GetElement(0UL) < 0) {
4140 // failed selection
4141 hy_env::EnvVariableSet(hy_env::selection_strings, new HY_NULL_RETURN, false);
4142 if (number_of_choices == 1L) {
4143 receptacle->SetValue (new _Constant (-1.), false,true,NULL);
4144 } else {
4145 receptacle->SetValue (new _Matrix (_SimpleList (1UL,-1L,0)), false,true,NULL);
4146 }
4147 terminate_execution = true;
4148 } else {
4149 if (variable_number || number_of_choices > 1L) {
4150 _SimpleList corrected_for_exclusions;
4151 _List chosen_strings;
4152 for (unsigned long i = 0UL; i < selections.countitems(); i++) {
4153 corrected_for_exclusions << excluded.SkipCorrect(selections.Element(i));
4154 chosen_strings < new _FString (*(_String*)available_choices->GetItem(selections.Element(i), 0),false);
4155 }
4156 receptacle->SetValue (new _Matrix (corrected_for_exclusions), false, true, NULL);
4157 } else {
4158 receptacle->SetValue (new _Constant (excluded.SkipCorrect (selections.Element(0UL))), false, true, NULL);
4159 hy_env::EnvVariableSet(hy_env::selection_strings, new _FString (*(_String*)available_choices->GetItem(selections.Element(0UL), 0),false), false);
4160 }
4161 }
4162
4163
4164 } catch (const _String& error) {
4165 return _DefaultExceptionHandler (receptacle, error, current_program);
4166 }
4167
4168 return true;
4169
4170 }
4171 //____________________________________________________________________________________
4172 // REQUIRES CODE REVIEW FROM THIS POINT ON
4173 //____________________________________________________________________________________
4174
ExecuteCase38(_ExecutionList & chain,bool sample)4175 void _ElementaryCommand::ExecuteCase38 (_ExecutionList& chain, bool sample) {
4176 chain.currentCommand++;
4177
4178 _List local_object_manager;
4179
4180 _String *likef = GetIthParameter(1);
4181
4182 _String name2lookup = AppendContainerName(*likef,chain.nameSpacePrefix);
4183 long objectID = FindLikeFuncName (name2lookup);
4184 try {
4185 if (objectID >= 0) {
4186 _DataSet * ds = new _DataSet;
4187 _String * dsName = new _String (AppendContainerName(*(_String*)parameters(0),chain.nameSpacePrefix));
4188 _LikelihoodFunction *lf = ((_LikelihoodFunction*)likeFuncList(objectID));
4189
4190 local_object_manager < ds < dsName;
4191
4192 _Matrix * partitionList = nil;
4193 if (parameters.lLength>2) {
4194 _String secondArg = *GetIthParameter(2);
4195 local_object_manager < (partitionList = (_Matrix*)ProcessAnArgumentByType (&secondArg, chain.nameSpacePrefix, MATRIX));
4196 }
4197 _SimpleList partsToDo;
4198
4199 if (lf->ProcessPartitionList(partsToDo, partitionList)) {
4200 lf->ReconstructAncestors(*ds, partsToDo, *dsName, sample, simpleParameters.Find(-1) >= 0, simpleParameters.Find(-2) >= 0 );
4201 }
4202
4203 ds->AddAReference();
4204 StoreADataSet (ds, dsName);
4205
4206 } else {
4207 objectID = FindSCFGName (name2lookup);
4208 if (objectID>=0)
4209 /* reconstruct best parse tree for corpus using SCFG */
4210 {
4211 CheckReceptacleAndStore (&AppendContainerName(*(_String*)parameters(0),chain.nameSpacePrefix)," ReconstructAncestors (SCFG)", true, new _FString( ((Scfg*)scfgList (objectID))->BestParseTree() ), false);
4212 } else {
4213 throw (_String("Likelihood Function/SCFG") & *likef & _String(" has not been initialized"));
4214 }
4215 }
4216 } catch (const _String& error) {
4217 _DefaultExceptionHandler (nil, error, chain);
4218 }
4219 }
4220
4221
4222 //____________________________________________________________________________________
4223
ExecuteCase31(_ExecutionList & chain)4224 void _ElementaryCommand::ExecuteCase31 (_ExecutionList& chain) {
4225 // 20100312 SLKP: added matrix-expression based model
4226 // definitions
4227 chain.advance();
4228
4229 try {
4230 // first check to see if matrix parameters here are valid
4231
4232 bool usingLastDefMatrix = false,
4233 doExpressionBased = false;
4234
4235 _Formula *isExpressionBased = nil;
4236
4237 _String* parameterName,
4238 errMsg,
4239 arg0 = chain.AddNameSpaceToID(*(_String*)parameters(0));
4240
4241 long f,
4242 f2=-1L,
4243 matrixDim,
4244 f3,
4245 multFreqs = 1;
4246
4247
4248
4249 if (parameters.lLength>3) {
4250 parameterName = (_String*)parameters.list_data[3];
4251 if (parameterName->Equal(explicitFormMExp)) {
4252 doExpressionBased = true;
4253 multFreqs = 0;
4254 } else {
4255 multFreqs = ProcessNumericArgument (parameterName,chain.nameSpacePrefix);
4256 }
4257 }
4258
4259 _Matrix* checkMatrix = nil;
4260
4261 parameterName = (_String*)parameters.list_data[1];
4262
4263 if (parameterName->Equal (useLastDefinedMatrix)) {
4264 if (lastMatrixDeclared<0) {
4265 throw _String ("First Call to Model. USE_LAST_DEFINED_MATRIX is meaningless.");
4266 }
4267 f3 = lastMatrixDeclared;
4268 f = modelMatrixIndices[f3];
4269 usingLastDefMatrix = true;
4270 }
4271
4272
4273 if (doExpressionBased) {
4274 _StringBuffer matrixExpression (ProcessLiteralArgument((_String*)parameters.list_data[1],chain.nameSpacePrefix)),
4275 defErrMsg = _String ("The expression for the explicit matrix exponential passed to Model must be a valid matrix-valued HyPhy formula that is not an assignment") & ':' & matrixExpression;
4276 // try to parse the expression, confirm that it is a square matrix,
4277 // and that it is a valid transition matrix
4278 isExpressionBased = new _Formula;
4279 _FormulaParsingContext fpc (nil, chain.nameSpacePrefix);
4280 _StringBuffer trimmed_expression;
4281 _ElementaryCommand::FindNextCommand (matrixExpression,trimmed_expression);
4282 matrixExpression = trimmed_expression;
4283 long parseCode = Parse(isExpressionBased,matrixExpression,fpc, nil);
4284 if (parseCode != HY_FORMULA_EXPRESSION || isExpressionBased->ObjectClass()!= MATRIX ) {
4285 throw (defErrMsg & " parse code = " & parseCode & " " & (parseCode == HY_FORMULA_EXPRESSION ? (_String(", object type code ") & _String((long) isExpressionBased->ObjectClass())) : kEmptyString ));
4286 }
4287
4288 //for (unsigned long k = 0; k < isExpressionBased
4289
4290 checkMatrix = (_Matrix*)isExpressionBased->Compute();
4291
4292
4293 } else {
4294 parameterName = (_String*)parameters.list_data[1];
4295
4296 _String augName (chain.AddNameSpaceToID(*parameterName));
4297 f = LocateVarByName (augName);
4298
4299 if (f<0) {
4300 throw (*parameterName & " has not been defined prior to the call to Model = ...");
4301 }
4302
4303 _Variable* checkVar = usingLastDefMatrix?LocateVar(f):FetchVar (f);
4304 if (checkVar->ObjectClass()!=MATRIX) {
4305 throw (*parameterName & " must refer to a matrix in the call to Model = ...");
4306 }
4307 checkMatrix = (_Matrix*)checkVar->GetValue();
4308 }
4309
4310
4311
4312 // so far so good
4313 matrixDim = checkMatrix->GetHDim();
4314 if ( matrixDim!=checkMatrix->GetVDim() || matrixDim<2 ) {
4315 throw (*parameterName & " must be a square matrix of dimension>=2 in the call to Model = ...");
4316 }
4317
4318 parameterName = (_String*)parameters.list_data[2]; // this is the frequency matrix (if there is one!)
4319 _String freqNameTag (chain.AddNameSpaceToID(*parameterName));
4320
4321 f2 = LocateVarByName (freqNameTag);
4322 if (f2<0) {
4323 throw(*parameterName & " has not been defined prior to the call to Model = ...");
4324 }
4325 _Variable * checkVar = FetchVar (f2);
4326 if (checkVar->ObjectClass()!=MATRIX) {
4327 throw (*parameterName & " must refer to a column/row vector in the call to Model = ...");
4328 }
4329
4330 checkMatrix = (_Matrix*)checkVar->GetValue();
4331
4332 if (checkMatrix->GetVDim()==1UL) {
4333 if (checkMatrix->GetHDim()!=matrixDim) {
4334 throw (*parameterName & " must be a column vector of the same dimension as the model matrix in the call to Model = ...");
4335 }
4336 } else if (checkMatrix->GetHDim()==1UL) {
4337 if (checkMatrix->GetVDim()!=matrixDim) {
4338 throw ( *parameterName & " must be a row vector of the same dimension as the model matrix in the call to Model = ...");
4339 }
4340 errMsg = *parameterName & " has been transposed to the default column vector setting ";
4341 checkMatrix->Transpose();
4342 ReportWarning (errMsg);
4343 } else {
4344 throw (*parameterName & " must refer to a column/row vector in the call to Model = ...");
4345 }
4346
4347 if (usingLastDefMatrix) {
4348 if (modelFrequenciesIndices[f3]<0) {
4349 f2 = -f2-1;
4350 }
4351 } else if (multFreqs == 0) { // optional flag present
4352 f2 = -f2-1;
4353 }
4354
4355 long existingIndex = modelNames.FindObject(&arg0);
4356
4357 if (existingIndex == -1) { // name not found
4358 lastMatrixDeclared = modelNames.FindObject (&kEmptyString);
4359
4360 if (lastMatrixDeclared>=0) {
4361 modelNames.Replace (lastMatrixDeclared,&arg0,true);
4362 modelTypeList.list_data[lastMatrixDeclared] = isExpressionBased?matrixDim:0;
4363 if (isExpressionBased) {
4364 modelMatrixIndices.list_data[lastMatrixDeclared] = (long)isExpressionBased;
4365 } else {
4366 modelMatrixIndices.list_data[lastMatrixDeclared] = (usingLastDefMatrix?f:variableNames.GetXtra(f));
4367 }
4368
4369 if (f2>=0) {
4370 modelFrequenciesIndices.list_data[lastMatrixDeclared] = variableNames.GetXtra(f2);
4371 } else {
4372 modelFrequenciesIndices.list_data[lastMatrixDeclared] = -variableNames.GetXtra(-f2-1)-1;
4373 }
4374 } else {
4375 modelNames && & arg0;
4376 modelTypeList << (isExpressionBased?matrixDim:0);
4377 if (isExpressionBased) {
4378 modelMatrixIndices << (long)isExpressionBased;
4379 } else {
4380 modelMatrixIndices << (usingLastDefMatrix?f:variableNames.GetXtra(f));
4381 }
4382 if (f2>=0) {
4383 modelFrequenciesIndices << variableNames.GetXtra(f2);
4384 } else {
4385 modelFrequenciesIndices << -variableNames.GetXtra(-f2-1)-1;
4386 }
4387 lastMatrixDeclared = modelNames.lLength-1;
4388 }
4389 } else {
4390 modelNames.Replace(existingIndex,&arg0,true);
4391 if (modelTypeList.list_data[existingIndex]) {
4392 delete ((_Formula*)modelMatrixIndices[existingIndex]);
4393 }
4394
4395 modelTypeList.list_data[existingIndex] = isExpressionBased?matrixDim:0;
4396 if (isExpressionBased) {
4397 modelMatrixIndices[existingIndex] = (long)isExpressionBased;
4398 } else {
4399 modelMatrixIndices[existingIndex] = usingLastDefMatrix?f:variableNames.GetXtra(f);
4400 }
4401
4402
4403 if (f2>=0) {
4404 modelFrequenciesIndices[existingIndex] = variableNames.GetXtra(f2);
4405 } else {
4406 modelFrequenciesIndices[existingIndex] = -variableNames.GetXtra(-f2-1)-1;
4407 }
4408
4409 lastMatrixDeclared = existingIndex;
4410 }
4411 } catch (const _String& error) {
4412 _DefaultExceptionHandler (nil, error, chain);
4413 }
4414 }
4415
4416
4417
4418
4419
4420