1LoadFunctionLibrary("libv3/all-terms.bf"); 2LoadFunctionLibrary("../models/model_functions.bf"); 3LoadFunctionLibrary("../models/DNA/GTR.bf"); 4LoadFunctionLibrary("../convenience/regexp.bf"); 5LoadFunctionLibrary("mpi.bf"); 6LoadFunctionLibrary("libv3/convenience/math.bf"); 7LoadFunctionLibrary("libv3/IOFunctions.bf"); 8 9/** 10 * @name estimators.TakeLFStateSnapshot 11 * @param {String} lf_id 12 * @returns {Dict} parameter -> {"MLE" : value, "constraint" : string (if present)} 13 */ 14 15lfunction estimators.TakeLFStateSnapshot(lf_id) { 16 snapshot = {}; 17 GetString (info, ^lf_id,-1); 18 19 utility.ForEach (info[utility.getGlobalValue("terms.parameters.global_independent")], "_name_", 20 '`&snapshot`[_name_] = {terms.fit.MLE : Eval (_name_)};'); 21 utility.ForEach (info[utility.getGlobalValue("terms.parameters.local_independent")], "_name_", 22 '`&snapshot`[_name_] = {terms.fit.MLE : Eval (_name_)};'); 23 utility.ForEach (info[utility.getGlobalValue("terms.parameters.global_constrained")], "_name_", 24 '`&snapshot`[_name_] = {terms.fit.MLE : Eval (_name_), 25 terms.constraint : parameters.GetConstraint (_name_)};'); 26 utility.ForEach (info[utility.getGlobalValue("terms.parameters.local_constrained")], "_name_", 27 '`&snapshot`[_name_] = {terms.fit.MLE : Eval (_name_), 28 terms.constraint : parameters.GetConstraint (_name_)};'); 29 30 return snapshot; 31} 32 33lfunction estimators.RestoreLFStateFromSnapshot(lf_id, snapshot) { 34 p_names = utility.Keys (snapshot); 35 p_count = utility.Array1D (p_names); 36 37 for (k = 0; k < p_count; k += 1) { 38 _name_ = p_names [k]; 39 _info_ = snapshot [_name_]; 40 if (_info_ / ^"terms.constraint") { 41 parameters.SetConstraint (_name_, _info_ [^"terms.constraint"], ""); 42 } else { 43 parameters.SetValue (_name_, _info_ [^"terms.fit.MLE"]); 44 } 45 } 46} 47 48lfunction estimators.ConstrainAndRunLRT (lf_id, constraint) { 49 savedMLES = estimators.TakeLFStateSnapshot (lf_id); 50 currentLL = estimators.ComputeLF (lf_id); 51 52 df = Call (constraint, TRUE); 53 Optimize (res, ^lf_id); 54 lrt = math.DoLRT (res[1][0],currentLL,df); 55 56 estimators.RestoreLFStateFromSnapshot (lf_id, savedMLES); 57 58 Call (constraint, FALSE); 59 60 return lrt; 61} 62 63/** 64 * @name estimators.GetGlobalMLE 65 * @param {Dictionary} results 66 * @param {String} tag 67 * @returns None 68 */ 69lfunction estimators.GetGlobalMLE(results, tag) { 70 estimate = (results[ utility.getGlobalValue("terms.global")])[tag]; 71 72 if (Type(estimate) == "AssociativeList") { 73 return estimate[ utility.getGlobalValue("terms.fit.MLE")]; 74 } 75 return None; 76} 77 78/** 79 * @name estimators.GetBranchEstimates 80 * @param {Dictionary} results 81 * @param {Number} partiton_index 82 * @param {String} node_name 83 * @returns None 84 */ 85lfunction estimators.GetBranchEstimates (results, partition_index, node_name) { 86 87 estimate = ((results[ utility.getGlobalValue("terms.branch_length")])[partition_index])[node_name]; 88 89 if (Type(estimate) == "AssociativeList") { 90 return estimate; 91 } 92 return None; 93} 94 95 96/** 97 * Extract global scope parameter estimates that match a regular expression 98 * @name estimators.GetGlobalMLE_RegExp 99 * @param {Dictionary} results 100 * @param {String} regular expression to match 101 * @returns {Dict} parameter description : value (could be empty) 102 */ 103lfunction estimators.GetGlobalMLE_RegExp(results, re) { 104 105 names = utility.Filter (utility.Keys (results[ utility.getGlobalValue("terms.global")]), 106 "_parameter_description_", 107 "None != regexp.Find (_parameter_description_, `&re`)"); 108 109 result = {}; 110 count = utility.Array1D (names); 111 for (k = 0; k < count; k += 1) { 112 result [names[k]] = (results[ utility.getGlobalValue("terms.global")])[names[k]]; 113 } 114 115 return result; 116} 117 118/** 119 * @name estimators.copyGlobals2 120 * @private 121 * @param {String} key2 122 * @param {String} value2 123 * @returns nothing 124 */ 125function estimators.copyGlobals2(key2, value2) { 126 127 if (Type((estimators.ExtractMLEs.results[terms.global])[key2]) == "AssociativeList") { 128 key2 = "[`key`] `key2`"; 129 // this parameter has already been defined, need to prefix with model name 130 } 131 132 (estimators.ExtractMLEs.results[terms.global])[key2] = { 133 terms.id: value2, 134 terms.fit.MLE: Eval(value2) 135 }; 136 137 if (parameters.IsIndependent(value2) != TRUE) { 138 ((estimators.ExtractMLEs.results[terms.global])[key2])[terms.constraint] = parameters.GetConstraint(value2); 139 } 140} 141 142/** 143 * @name estimators.copyGlobals 144 * @private 145 * @param {String} key 146 * @param {String} value 147 * @returns nothing 148 */ 149function estimators.copyGlobals(key, value) { 150 ((value[terms.parameters])[terms.global])["estimators.copyGlobals2"][""]; 151} 152 153/** 154 * @name estimators.CopyFrequencies 155 * @private 156 * @param {String} key 157 * @param {Dictionary} value 158 * @returns nothing 159 */ 160function estimators.CopyFrequencies(model_name, model_decription) { 161 (estimators.ExtractMLEs.results[terms.efv_estimate])[model_name] = model_decription[terms.efv_estimate]; 162} 163 164 165function estimators.SetGlobals2(key2, value) { 166 167 168 if (Type(estimators.ApplyExistingEstimates.set_globals[key2]) == "AssociativeList") { 169 key3 = "[`key`] `key2`"; 170 } else { 171 key3 = key2; 172 } 173 174 175 176 177 178 estimators.ApplyExistingEstimates.set_globals[key3] = { 179 terms.id: key3, 180 terms.fit.MLE: value 181 }; 182 183 __init_value = (initial_values[terms.global])[key3]; 184 185 if (Type(__init_value) != "AssociativeList") { 186 __init_value = (initial_values[terms.global])[key2]; 187 } 188 189 190 if (Type(__init_value) == "AssociativeList") { 191 if (__init_value[terms.fix]) { 192 estimators.ApplyExistingEstimates.df_correction += parameters.IsIndependent(value); 193 ExecuteCommands("`value` := " + __init_value[terms.fit.MLE]); 194 //_did_set [value] = 1; 195 } else { 196 if (parameters.IsIndependent (value)) { 197 ExecuteCommands("`value` = " + __init_value[terms.fit.MLE]); 198 //_did_set [value] = 1; 199 } else { 200 messages.log (value + " was already constrained in estimators.SetGlobals2"); 201 } 202 } 203 } 204} 205 206/** 207 * @name estimators.SetGlobals 208 * @param {String} key 209 * @param {String} value 210 * @returns nothing 211 */ 212function estimators.SetGlobals(key, value) { 213 /*_did_set = {}; 214 for (i,v; in; ((value[terms.parameters])[terms.global])) { 215 _did_set [v] = 0; 216 }*/ 217 218 ((value[terms.parameters])[terms.global])["estimators.SetGlobals2"][""]; 219 220 //console.log (_did_set); 221} 222 223/** 224 * @name estimators.SetCategory 225 * @param {String} key 226 * @param {String} value 227 * @returns nothing 228 */ 229function estimators.SetCategory (key, value) { 230 ((value[terms.parameters])[terms.category])["estimators.SetCategory2"][""]; 231} 232 233function estimators.SetCategory2(key, value) { 234 ^key = 1; 235} 236 237/** 238 * @name estimators.RemoveBranchLengthConstraints 239 * @param {Dict} estimates 240 */ 241lfunction estimators.RemoveBranchLengthConstraints (estimates) { 242 for (part; in; estimates[^"terms.branch_length"]) { 243 for (b; in; part) { 244 for (p; in; b) { 245 if (Type (p) == "AssociativeList") { 246 p - ^"terms.constraint"; 247 } 248 } 249 250 } 251 } 252 return estimates; 253} 254 255/** 256 * @name estimators.RemoveGlobalConstraints 257 * @param {Dict} estimates 258 */ 259lfunction estimators.RemoveGlobalConstraints (estimates) { 260 for (pp; in; estimates[^"terms.global"]) { 261 if (Type (pp) == "AssociativeList") { 262 pp - ^"terms.constraint"; 263 } 264 } 265 return estimates; 266} 267 268/** 269 * @name estimators.ExtractBranchInformation.copy_local 270 * @param {String} key 271 * @param {String} value 272 */ 273function estimators.ExtractBranchInformation.copy_local(key, value) { 274 275 estimators.ExtractBranchInformation.copy_local.var_name = estimators.extractBranchLength.parameter_tag + "." + value; 276 277 estimators.extractBranchLength.result[key] = {}; 278 (estimators.extractBranchLength.result[key])[terms.id] = value; 279 (estimators.extractBranchLength.result[key])[terms.fit.MLE] = Eval(estimators.ExtractBranchInformation.copy_local.var_name); 280 /* terms.id: value, 281 terms.fit.MLE: Eval(estimators.ExtractBranchInformation.copy_local.var_name) 282 };*/ 283 284 if (parameters.IsIndependent(estimators.ExtractBranchInformation.copy_local.var_name) != TRUE) { 285 (estimators.extractBranchLength.result[key])[terms.constraint] = parameters.GetConstraint(estimators.ExtractBranchInformation.copy_local.var_name); 286 } 287 288} 289 290/** 291 * @name estimators.ExtractBranchInformation 292 * @param {String} type 293 * @param {String} node 294 * @param {String} model 295 * @returns {Dictionary} branch information 296 */ 297function estimators.ExtractBranchInformation(tree, node, model) { 298 estimators.extractBranchLength.result = {}; 299 300 if (Abs(model[terms.model.get_branch_length])) { 301 estimators.extractBranchLength.result[terms.fit.MLE] = Call(model[terms.model.get_branch_length], model, tree, node); 302 } else { 303 estimators.extractBranchLength.result[terms.fit.MLE] = Eval("BranchLength (`tree`, \"`node`\")"); 304 } 305 306 estimators.extractBranchLength.parameter_tag = tree + "." + node; 307 (model.parameters.Local(model))["estimators.ExtractBranchInformation.copy_local"][""]; 308 309 return estimators.extractBranchLength.result; 310} 311 312/** 313 * @name estimators.applyBranchLength 314 * @private 315 * @param {String} tree 316 * @param {String} node 317 * @param {String} model 318 * @param {String} length 319 */ 320function estimators.applyBranchLength(tree, node, model, length) { 321 return Call(model[terms.model.set_branch_length], model, length, tree + "." + node); 322} 323 324/** 325 * @name estimators.constrainBranchLength 326 * @private 327 * @param {String} tree 328 * @param {String} node 329 * @param {String} model 330 * @param {String} length 331 */ 332function estimators.constrainBranchLength(tree, node, model, length) { 333 return Call(model[terms.model.constrain_branch_length], model, length, tree + "." + node); 334} 335 336/** 337 * @name estimators.fixSubsetOfEstimates.helper 338 * @private 339 * @param {String} key 340 * @param {String} value 341 */ 342function estimators.fixSubsetOfEstimates.helper(key, value) { 343 if (value / terms.constraint == FALSE) { 344 estimators.fixSubsetOfEstimates.fixed + value[terms.id]; 345 value[terms.fix] = TRUE; 346 } 347} 348 349/** 350 * @name estimators.fixSubsetOfEstimates.helper_condition 351 * @private 352 * @param {String} key 353 */ 354function estimators.fixSubsetOfEstimates.helper_condition(key) { 355 return Type(variables[key]) != Unknown; 356} 357 358/** 359 * @name estimators.fixSubsetOfEstimates 360 * @private 361 * @param {String} estimates 362 * @param {String} variables 363 * @returns nothing 364 */ 365function estimators.fixSubsetOfEstimates(estimates, variables) { 366 estimators.fixSubsetOfEstimates.fixed = {}; 367 (estimates[terms.global])["estimators.fixSubsetOfEstimates.helper"]["estimators.fixSubsetOfEstimates.helper_condition"]; 368 return estimators.fixSubsetOfEstimates.fixed; 369} 370 371/** 372 * @name estimators.branch_lengths_in_string.map 373 * @private 374 * @param {String} id 375 * @param {String} value 376 */ 377function estimators.branch_lengths_in_string.map(id, value) { 378 estimators.branch_lengths_in_string.lookup[id] = value[terms.fit.MLE]; 379} 380 381/** 382 * @name estimators.branch_lengths_in_string 383 * @param {String} tree_id 384 * @param {String} lookup 385 * @returns {String} branch lenghts in tree string 386 */ 387function estimators.branch_lengths_in_string(tree_id, lookup) { 388 estimators.branch_lengths_in_string.lookup = {}; 389 lookup["estimators.branch_lengths_in_string.map"][""]; 390 utility.ToggleEnvVariable("BRANCH_LENGTH_STENCIL", estimators.branch_lengths_in_string.lookup); 391 estimators.branch_lengths_in_string.string = Eval("Format (`tree_id`,1,1)"); 392 utility.ToggleEnvVariable("BRANCH_LENGTH_STENCIL", None); 393 return estimators.branch_lengths_in_string.string; 394} 395 396/** 397 * @name estimators.ExtractMLEs 398 * @param {String} likelihood_function_id 399 * @param {String} model_descriptions 400 * @returns results 401 */ 402function estimators.ExtractMLEs(likelihood_function_id, model_descriptions) { 403 return estimators.ExtractMLEsOptions (likelihood_function_id, model_descriptions, {}); 404} 405 406/** 407 * @name estimators.ExtractMLEsOptions 408 * @param {String} likelihood_function_id 409 * @param {String} model_descriptions 410 * @param {Dict} options 411 * @returns results 412 */ 413function estimators.ExtractMLEsOptions(likelihood_function_id, model_descriptions, options) { 414 415 ExecuteCommands("GetString (estimators.ExtractMLEs.lfInfo, `likelihood_function_id`,-1)"); 416 417 estimators.ExtractMLEs.results = {}; 418 estimators.ExtractMLEs.partitions = utility.Array1D(estimators.ExtractMLEs.lfInfo[terms.fit.trees]); 419 420 // copy global variables first 421 422 estimators.ExtractMLEs.results[terms.global] = {}; 423 model_descriptions["estimators.copyGlobals"][""]; 424 estimators.ExtractMLEs.results[terms.efv_estimate] = {}; 425 model_descriptions["estimators.CopyFrequencies"][""]; 426 estimators.ExtractMLEs.results[terms.branch_length] = {}; 427 estimators.ExtractMLEs.results[terms.fit.trees] = estimators.ExtractMLEs.lfInfo[terms.fit.trees]; 428 429 if (options[utility.getGlobalValue("terms.globals_only")]) { 430 return estimators.ExtractMLEs.results; 431 } 432 433 for (estimators.ExtractMLEs.i = 0; estimators.ExtractMLEs.i < estimators.ExtractMLEs.partitions; estimators.ExtractMLEs.i += 1) { 434 _tree_name = (estimators.ExtractMLEs.lfInfo[terms.fit.trees])[estimators.ExtractMLEs.i]; 435 436 437 GetInformation (estimators.ExtractMLEs.map, *_tree_name); 438 439 estimators.ExtractMLEs.branch_names = Rows(estimators.ExtractMLEs.map); 440 441 442 (estimators.ExtractMLEs.results[terms.branch_length])[estimators.ExtractMLEs.i] = {}; 443 444 for (estimators.ExtractMLEs.b = 0; estimators.ExtractMLEs.b < Abs(estimators.ExtractMLEs.map); estimators.ExtractMLEs.b += 1) { 445 _branch_name = estimators.ExtractMLEs.branch_names[estimators.ExtractMLEs.b]; 446 447 ((estimators.ExtractMLEs.results[terms.branch_length])[estimators.ExtractMLEs.i])[_branch_name] = 448 estimators.ExtractBranchInformation(_tree_name, _branch_name, model_descriptions[estimators.ExtractMLEs.map[_branch_name]]); 449 450 } 451 452 453 (estimators.ExtractMLEs.results[terms.fit.trees])[estimators.ExtractMLEs.i] = 454 estimators.branch_lengths_in_string((estimators.ExtractMLEs.results[terms.fit.trees])[estimators.ExtractMLEs.i], (estimators.ExtractMLEs.results[terms.branch_length])[estimators.ExtractMLEs.i]); 455 456 } 457 458 return estimators.ExtractMLEs.results; 459} 460 461/** 462 * @name estimators.TraverseLocalParameters 463 * @param {String} likelihood_function_id 464 * @param {Dictionary} model_descriptions 465 * @param {String} callback (tree, node, parameter_list) 466 467 */ 468lfunction estimators.TraverseLocalParameters (likelihood_function_id, model_descriptions, callback) { 469 GetString(lf_info, ^ likelihood_function_id, -1); 470 partitions = utility.Array1D(lf_info[utility.getGlobalValue ("terms.fit.trees")]); 471 result = {}; 472 for (i = 0; i < partitions; i += 1) { 473 tree_name = (lf_info[utility.getGlobalValue ("terms.fit.trees")])[i]; 474 GetInformation (map, ^tree_name); 475 branch_names = Rows (map); 476 for (b = 0; b < Abs(map); b += 1) { 477 _branch_name = branch_names[b]; 478 result[_branch_name] = Call (callback, tree_name, _branch_name, (model_descriptions [map[_branch_name]])[utility.getGlobalValue ("terms.parameters")], i); 479 } 480 } 481 return result; 482} 483 484/** 485 * @name 486 * @param {String} tree_name 487 * @param {Dictionary} model_descriptions 488 * @param {Matrix} initial_values 489 * @param branch_length_conditions 490 * @returns number of constrained parameters; 491 */ 492function estimators.ApplyExistingEstimatesToTree (_tree_name, model_descriptions, initial_values, _application_type, keep_track_of_proportional_scalers) { 493 494 SetParameter (DEFER_CONSTRAINT_APPLICATION, 1, 0); 495 496 497 estimators.ApplyExistingEstimatesToTree.constraint_count = 0; 498 499 500 ExecuteCommands("GetInformation (estimators.ApplyExistingEstimatesToTree.map, `_tree_name`);"); 501 estimators.ApplyExistingEstimatesToTree.branch_names = Rows(estimators.ApplyExistingEstimatesToTree.map); 502 503 for (estimators.ApplyExistingEstimatesToTree.b = 0; estimators.ApplyExistingEstimatesToTree.b < Abs(estimators.ApplyExistingEstimatesToTree.map); estimators.ApplyExistingEstimatesToTree.b += 1) { 504 _branch_name = estimators.ApplyExistingEstimatesToTree.branch_names[estimators.ApplyExistingEstimatesToTree.b]; 505 506 if (initial_values / _branch_name) { // have an entry for this branch name 507 _existing_estimate = initial_values[_branch_name]; 508 509 if (Type(_existing_estimate) == "AssociativeList") { 510 _set_branch_length_to = (initial_values[_branch_name])[terms.fit.MLE]; 511 if (None != branch_length_conditions) { 512 if (None != _application_type) { 513 514 if (Type(_application_type) == "String") { 515 if (_application_type == terms.model.branch_length_constrain ) { 516 estimators.ApplyExistingEstimatesToTree.constraint_count += estimators.constrainBranchLength(_tree_name, _branch_name, model_descriptions[estimators.ApplyExistingEstimatesToTree.map[_branch_name]], _set_branch_length_to); 517 continue; 518 } 519 _set_branch_length_to = {}; 520 _set_branch_length_to[terms.branch_length] = _existing_estimate[terms.fit.MLE]; 521 _set_branch_length_to[terms.model.branch_length_scaler] = _application_type; 522 keep_track_of_proportional_scalers[_application_type] = 1; 523 524 } 525 } 526 } 527 528 estimators.ApplyExistingEstimatesToTree.constraint_count += estimators.applyBranchLength(_tree_name, _branch_name, model_descriptions[estimators.ApplyExistingEstimatesToTree.map[_branch_name]], _set_branch_length_to); 529 } else { 530 if (Type(_existing_estimate) != "Unknown") { 531 warning.log ("Incorrect type for the initial values object of for branch '" + _branch_name + "' : " + _existing_estimate); 532 } 533 } 534 } else { 535 //warning.log ("No initial branch length object of for branch '" + _branch_name); 536 } 537 } 538 539 SetParameter (DEFER_CONSTRAINT_APPLICATION, 0, 0); 540 541 542 //fprintf (stdout, Format (^_tree_name, 1,1), "\n"); 543 544 return estimators.ApplyExistingEstimatesToTree.constraint_count; 545} 546 547/** 548 * @name 549 * @param {String} likelihood_function_id 550 * @param {Dictionary} model_descriptions 551 * @param {Matrix} initial_values 552 * @param branch_length_conditions 553 * @returns estimators.ApplyExistingEstimates.df_correction - Abs(estimators.ApplyExistingEstimates.keep_track_of_proportional_scalers); 554 */ 555function estimators.ApplyExistingEstimates(likelihood_function_id, model_descriptions, initial_values, branch_length_conditions) { 556 //fprintf (stdout, model_descriptions, "\n", initial_values, "\n"); 557 558 /* set all category variable values to one */ 559 560 GetString(estimators.ApplyExistingEstimates.lfInfo, ^ likelihood_function_id, -1); 561 estimators.ApplyExistingEstimates.results = {}; 562 estimators.ApplyExistingEstimates.partitions = utility.Array1D(estimators.ApplyExistingEstimates.lfInfo[terms.fit.trees]); 563 564 565 estimators.ApplyExistingEstimates.df_correction = 0; 566 // copy global variables first 567 568 estimators.ApplyExistingEstimates.results[terms.global] = {}; 569 model_descriptions["estimators.SetCategory"][""]; 570 // the above line traverses all model descriptions and sets 571 // the _value_ of category variables to 1, so that we can 572 // compute branch lengths 573 574 estimators.ApplyExistingEstimates.set_globals = {}; 575 model_descriptions["estimators.SetGlobals"][""]; 576 577 578 if (Type(branch_length_conditions) == "String") { 579 if (branch_length_conditions == terms.globals_only) { 580 return estimators.ApplyExistingEstimates.df_correction; 581 } 582 assert("0", "Unsupported value for 'branch_length_conditions' in estimators.ApplyExistingEstimates"); 583 return 0; 584 } 585 586 estimators.ApplyExistingEstimates.keep_track_of_proportional_scalers = {}; 587 588 for (estimators.ApplyExistingEstimates.i = 0; estimators.ApplyExistingEstimates.i < estimators.ApplyExistingEstimates.partitions; estimators.ApplyExistingEstimates.i += 1) { 589 590 if (Type((initial_values[terms.branch_length])[estimators.ApplyExistingEstimates.i]) == "AssociativeList") { // have branch lengths for this partition 591 592 _application_type = None; 593 594 if (Type (branch_length_conditions) == "AssociativeList") { 595 if (Abs(branch_length_conditions) > estimators.ApplyExistingEstimates.i) { 596 _application_type = branch_length_conditions[estimators.ApplyExistingEstimates.i]; 597 } 598 } 599 600 estimators.ApplyExistingEstimates.df_correction += estimators.ApplyExistingEstimatesToTree ((estimators.ApplyExistingEstimates.lfInfo[terms.fit.trees])[estimators.ApplyExistingEstimates.i], 601 model_descriptions, 602 (initial_values[terms.branch_length])[estimators.ApplyExistingEstimates.i], 603 _application_type, 604 estimators.ApplyExistingEstimates.keep_track_of_proportional_scalers); 605 606 607 608 609 } else { 610 if (Type((initial_values[terms.branch_length])[estimators.ApplyExistingEstimates.i]) != "Unknown") { 611 warning.log ("Incorrect type for the initial values object for partition " + estimators.ApplyExistingEstimates.i 612 + ". " + (initial_values[terms.branch_length])[estimators.ApplyExistingEstimates.i]); 613 } 614 } 615 } // have branch lengths for this partition 616 617 618 return estimators.ApplyExistingEstimates.df_correction - Abs(estimators.ApplyExistingEstimates.keep_track_of_proportional_scalers); 619} 620 621/** 622 * @name estimators._aux.countEmpiricalParameters 623 * @private 624 * @param {String} id 625 * @param {Dictionary} model 626 * @returns nothing 627 */ 628function estimators._aux.countEmpiricalParameters(id, model) { 629 estimators._aux.parameter_counter += (model[terms.parameters])[terms.model.empirical]; 630} 631 632/** 633 * Fits a LikelihoodFunction 634 * @name estimators.FitExistingLF 635 * @param model_map 636 * @returns LF results 637 */ 638lfunction estimators.FitExistingLF (lf_id, model_objects) { 639 640 utility.ToggleEnvVariable("USE_LAST_RESULTS", TRUE); 641 Optimize (mles, ^lf_id); 642 utility.ToggleEnvVariable("USE_LAST_RESULTS", None); 643 644 results = estimators.ExtractMLEs( lf_id, model_objects); 645 results[utility.getGlobalValue ("terms.fit.log_likelihood")] = mles[1][0]; 646 results[utility.getGlobalValue ("terms.parameters")] = mles[1][1]; 647 648 return results; 649} 650 651/** 652 * Makes a likelihood function object with the desired parameters 653 * @name estimators.FitLF 654 * @param {Matrix} data_filters_list - a vector of {DataFilter}s 655 * @param {Matrix} tree_list - a vector of {Tree}s 656 * @param model_map 657 * @param initial_values 658 * @returns LF results 659 */ 660 661lfunction estimators.BuildLFObject (lf_id, data_filter, tree, model_map, initial_values, model_objects, run_options) { 662 663 if (Type(data_filter) == "String") { 664 return estimators.FitLF ({ 665 { 666 data_filter__ 667 } 668 }, { 669 "0": tree 670 }, 671 { 672 "0" : model_map 673 }, 674 initial_values, model_objects, run_options); 675 } 676 677 components = utility.Array1D(data_filter); 678 679 680 lf_components = { 681 2 * components, 682 1 683 }; 684 685 686 for (i = 0; i < components; i += 1) { 687 lf_components[2 * i] = data_filter[i]; 688 lf_components[2 * i + 1] = &tree_id + "_" + i; 689 model.ApplyModelToTree(lf_components[2*i + 1], tree[i], None, model_map[i]); 690 } 691 692 693 utility.ExecuteInGlobalNamespace ("LikelihoodFunction `lf_id` = (`&lf_components`)"); 694 695 696 697 df = 0; 698 699 if (Type(initial_values) == "AssociativeList") { 700 utility.ToggleEnvVariable("USE_LAST_RESULTS", 1); 701 df = estimators.ApplyExistingEstimates(lf_id, model_objects, initial_values, run_options[utility.getGlobalValue("terms.run_options.proportional_branch_length_scaler")]); 702 } 703 704 if (utility.Has (run_options,utility.getGlobalValue("terms.run_options.apply_user_constraints"),"String")) { 705 df += Call (run_options[utility.getGlobalValue("terms.run_options.apply_user_constraints")], lf_id, lf_components, data_filter, tree, model_map, initial_values, model_objects); 706 } 707 708 return estimators.ExtractMLEs( lf_id , model_objects); 709} 710 711/** 712 * Fits a LikelihoodFunction 713 * @name estimators.FitLF 714 * @param {Matrix} data_filters_list - a vector of {DataFilter}s 715 * @param {Matrix} tree_list - a vector of {Tree}s 716 * @param model_map 717 * @param initial_values 718 * @returns LF results 719 */ 720lfunction estimators.FitLF(data_filter, tree, model_map, initial_values, model_objects, run_options) { 721 722 723 if (Type(data_filter) == "String") { 724 return estimators.FitLF ({ 725 { 726 data_filter__ 727 } 728 }, { 729 "0": tree 730 }, 731 { 732 "0" : model_map 733 }, 734 initial_values, model_objects, run_options); 735 } 736 737 components = utility.Array1D(data_filter); 738 739 740 lf_components = { 741 2 * components, 742 1 743 }; 744 745 746 747 for (i = 0; i < components; i += 1) { 748 lf_components[2 * i] = data_filter[i]; 749 lf_components[2 * i + 1] = &tree_id + "_" + i; 750 model.ApplyModelToTree(lf_components[2*i + 1], tree[i], None, model_map[i]); 751 } 752 753 754 755 lf_id = &likelihoodFunction; 756 utility.ExecuteInGlobalNamespace ("LikelihoodFunction `lf_id` = (`&lf_components`)"); 757 758 759 df = 0; 760 761 if (utility.Has (run_options,utility.getGlobalValue("terms.run_options.apply_user_constraints"),"String")) { 762 df += Call (run_options[utility.getGlobalValue("terms.run_options.apply_user_constraints")], lf_id, lf_components, data_filter, tree, model_map, initial_values, model_objects); 763 } 764 765 if (Type(initial_values) == "AssociativeList") { 766 utility.ToggleEnvVariable("USE_LAST_RESULTS", 1); 767 //console.log (initial_values); 768 df = estimators.ApplyExistingEstimates("`&likelihoodFunction`", model_objects, initial_values, run_options[utility.getGlobalValue("terms.run_options.proportional_branch_length_scaler")]); 769 } 770 771 772 can_do_restarts = null; 773 774 775 776 if (utility.Has (run_options, utility.getGlobalValue("terms.search_grid"),"AssociativeList")) { 777 grid_results = mpi.ComputeOnGrid (&likelihoodFunction, run_options [utility.getGlobalValue("terms.search_grid")], "mpi.ComputeOnGrid.SimpleEvaluator", "mpi.ComputeOnGrid.ResultHandler"); 778 if (utility.Has (run_options, utility.getGlobalValue("terms.search_restarts"),"Number")) { 779 restarts = run_options[utility.getGlobalValue("terms.search_restarts")]; 780 if (restarts > 1) { 781 grid_results = utility.DictToSortedArray (grid_results); 782 can_do_restarts = {}; 783 for (i = 1; i <= restarts; i += 1) { 784 can_do_restarts + (run_options [utility.getGlobalValue("terms.search_grid")])[grid_results[Rows(grid_results)-i][1]]; 785 } 786 } 787 } 788 if (null == can_do_restarts) { 789 best_value = Max (grid_results, 1); 790 parameters.SetValues ((run_options [utility.getGlobalValue("terms.search_grid")])[best_value["key"]]); 791 } 792 //console.log (best_value); 793 //console.log ((run_options [utility.getGlobalValue("terms.search_grid")])[best_value["key"]]); 794 //assert (0); 795 } 796 797 798 799 if (Type (can_do_restarts) == "AssociativeList") { 800 io.ReportProgressBar("", "Working on crude initial optimizations"); 801 bestlog = -1e100; 802 for (i = 0; i < Abs (can_do_restarts); i += 1) { 803 parameters.SetValues (can_do_restarts[i]); 804 if (utility.Has (run_options,utility.getGlobalValue("terms.run_options.optimization_settings"),"AssociativeList")) { 805 Optimize (mles, likelihoodFunction, run_options[utility.getGlobalValue("terms.run_options.optimization_settings")]); 806 } else { 807 Optimize (mles, likelihoodFunction); 808 } 809 if (mles[1][0] > bestlog) { 810 811 //console.log ("\n\n**BEST LOG**\n\n"); 812 bestlog = mles[1][0]; 813 results = estimators.ExtractMLEs( & likelihoodFunction, model_objects); 814 results[utility.getGlobalValue ("terms.fit.log_likelihood")] = mles[1][0]; 815 } 816 io.ReportProgressBar("", "Starting point " + (i+1) + "/" + Abs (can_do_restarts) + ". Best LogL = " + bestlog); 817 818 } 819 io.ClearProgressBar(); 820 } else { 821 if (utility.Has (run_options,utility.getGlobalValue("terms.run_options.optimization_settings"),"AssociativeList")) { 822 Optimize (mles, likelihoodFunction, run_options[utility.getGlobalValue("terms.run_options.optimization_settings")]); 823 } else { 824 Optimize (mles, likelihoodFunction); 825 } 826 results = estimators.ExtractMLEs( & likelihoodFunction, model_objects); 827 results[utility.getGlobalValue ("terms.fit.log_likelihood")] = mles[1][0]; 828 } 829 830 831 832 if (Type(initial_values) == "AssociativeList") { 833 utility.ToggleEnvVariable("USE_LAST_RESULTS", None); 834 } 835 836 837 results[utility.getGlobalValue ("terms.parameters")] = mles[1][1] + df; 838 839 results[utility.getGlobalValue ("terms.fit.filters")] = { 840 1, 841 components 842 }; 843 844 for (i = 0; i < components; i += 1) { 845 (results[utility.getGlobalValue ("terms.fit.filters")])[i] = lf_components[2 * i]; 846 847 } 848 849 if (run_options[utility.getGlobalValue("terms.run_options.retain_lf_object")]) { 850 results[utility.getGlobalValue("terms.likelihood_function")] = & likelihoodFunction; 851 } else { 852 DeleteObject(likelihoodFunction); 853 } 854 855 return results; 856} 857 858lfunction estimators.CreateLFObject (context, data_filter, tree, model_template, initial_values, run_options, model_objects) { 859 860 if (Type(data_filter) == "String") { 861 return estimators.CreateLFObject (context, { 862 { 863 data_filter__ 864 } 865 }, { 866 "0": tree 867 }, model_template, initial_values, run_options, model_objects); 868 } 869 870 components = utility.Array1D(data_filter); 871 872 filters = utility.Map({ 873 components, 874 1 875 }["_MATRIX_ELEMENT_ROW_"], "_value_", "''+ '`context`.nuc_data_' + _value_"); 876 877 lf_components = { 878 2 * components, 879 1 880 }; 881 882 883 for (i = 0; i < components; i += 1) { 884 lf_components[2 * i] = filters[i]; 885 DataSetFilter ^ (filters[i]) = CreateFilter( ^ (data_filter[i]), 1); 886 } 887 888 889 user_model_id = context + ".user_model"; 890 utility.ExecuteInGlobalNamespace ("`user_model_id` = 0"); 891 892 893 ^(user_model_id) = model.generic.DefineModel(model_template, context + ".model", { 894 "0": "terms.global" 895 }, filters, None); 896 897 898 for (i = 0; i < components; i += 1) { 899 lf_components[2 * i + 1] = "`context`.tree_" + i; 900 model.ApplyModelToTree(lf_components[2 * i + 1], tree[i], { 901 "default": ^(user_model_id) 902 }, None); 903 } 904 905 906 lfid = context + ".likelihoodFunction"; 907 908 909 utility.ExecuteInGlobalNamespace ("LikelihoodFunction `lfid` = (`&lf_components`)"); 910 df = 0; 911 if (Type(initial_values) == "AssociativeList") { 912 if (None == model_objects) { 913 model_objects = { 914 (^user_model_id)[utility.getGlobalValue ("terms.id")]: ^(user_model_id) 915 }; 916 } 917 utility.ToggleEnvVariable("USE_LAST_RESULTS", 1); 918 df = estimators.ApplyExistingEstimates(lfid, model_objects, initial_values, run_options[utility.getGlobalValue("terms.run_options.proportional_branch_length_scaler")]); 919 } 920 921 return df; 922} 923 924/** 925 * @name estimators.FitSingleModel_Ext 926 * @param {DataFilter} data_filter 927 * @param {Tree} tree 928 * @param {Dict} model 929 * @param {Matrix} initial_values 930 * @param {Dict} run_options 931 * @returns results 932 */ 933 934lfunction estimators.FitSingleModel_Ext (data_filter, tree, model_template, initial_values, run_options) { 935 936 this_namespace = (&_); 937 this_namespace = this_namespace[0][Abs (this_namespace)-3]; 938 939 df = estimators.CreateLFObject (this_namespace, data_filter, tree, model_template, initial_values, run_options, None); 940 941 /* 942 partition parameters into groups 943 */ 944 945 pg = utility.getEnvVariable ("PARAMETER_GROUPING"); 946 947 if (Type (pg) != "AssociativeList") { 948 GetString (params, likelihoodFunction,-1); 949 pg = {"0" : params["Global Independent"]}; 950 utility.ToggleEnvVariable ("PARAMETER_GROUPING", pg); 951 } 952 953 if (utility.Has (run_options,utility.getGlobalValue("terms.run_options.optimization_settings"),"AssociativeList")) { 954 Optimize (mles, likelihoodFunction, run_options[utility.getGlobalValue("terms.run_options.optimization_settings")]); 955 } else { 956 Optimize (mles, likelihoodFunction); 957 } 958 959 if (Type (pg) == "AssociativeList") { 960 utility.ToggleEnvVariable ("PARAMETER_GROUPING", None); 961 } 962 963 if (Type(initial_values) == "AssociativeList") { 964 utility.ToggleEnvVariable("USE_LAST_RESULTS", None); 965 } 966 967 model_id_to_object = { 968 (this_namespace + ".model"): user_model 969 }; 970 971 972 results = estimators.ExtractMLEs( & likelihoodFunction, model_id_to_object); 973 974 975 results[utility.getGlobalValue("terms.fit.log_likelihood")] = mles[1][0]; 976 results[utility.getGlobalValue("terms.parameters")] = mles[1][1] + (user_model [utility.getGlobalValue("terms.parameters")]) [utility.getGlobalValue("terms.model.empirical")] + df; 977 978 979 if (option[utility.getGlobalValue("terms.run_options.retain_model_object")]) { 980 results[utility.getGlobalValue("terms.model")] = model_id_to_object; 981 } 982 983 if (run_options[utility.getGlobalValue("terms.run_options.retain_lf_object")]) { 984 results[utility.getGlobalValue("terms.likelihood_function")] = & likelihoodFunction; 985 } else { 986 DeleteObject(likelihoodFunction); 987 } 988 989 return results; 990} 991 992/** 993 * @name estimators.FitGTR_Ext 994 * @param {DataFilter} data_filter 995 * @param {Tree} tree 996 * @param {Matrix} initial_values 997 * @param {Dict} run_options 998 * @returns results 999 */ 1000 1001lfunction estimators.FitGTR_Ext (data_filter, tree, initial_values, run_options) { 1002 return estimators.FitSingleModel_Ext (data_filter, tree, "models.DNA.GTR.ModelDescription", initial_values, run_options); 1003} 1004 1005/** 1006 * @name estimators.FitGTR 1007 * @param {DataFilter} data_filter 1008 * @param {Tree} tree 1009 * @param {Matrix} initial_values 1010 * @returns results 1011 */ 1012lfunction estimators.FitGTR(data_filter, tree, initial_values) { 1013 return estimators.FitGTR_Ext(data_filter, tree, initial_values, {}); 1014} 1015 1016/** 1017 * @name estimators.FitMGREV.set_partition_omega 1018 * @private 1019 * @param {String} key 1020 * @param {String} value 1021 */ 1022function estimators.FitMGREV.set_partition_omega(key, value) { 1023 Eval("estimators.FitMGREV.tree.`key`.`estimators.FitMGREV.alpha` = 0.1"); 1024 ExecuteCommands("estimators.FitMGREV.tree.`key`.`estimators.FitMGREV.beta`:=estimators.FitMGREV.tree.`key`.`estimators.FitMGREV.alpha`*" + estimators.FitMGREV.partitioned_omega.parameters[value]); 1025} 1026 1027/** 1028 * @name estimators.FitMGREV.set_partition_omega 1029 * @private 1030 * @param {String} key 1031 * @param {String} value 1032 */ 1033lfunction estimators.FitMGREVExtractComponentBranchLengths(codon_data, fit_results) { 1034 1035 //extract fitted trees with branch lengths scaled on synonymous and non-synonymous 1036 //substitutions per site 1037 1038 stencils = genetic_code.ComputeBranchLengthStencils(codon_data[^"terms.code"]); 1039 1040 utility.SetEnvVariable ("BRANCH_LENGTH_STENCIL", stencils[^"terms.genetic_code.synonymous"]); 1041 fit_results[^"terms.fit.synonymous_trees"] = (estimators.ExtractMLEs(fit_results[^"terms.likelihood_function"], fit_results[^"terms.model"]))[^"terms.fit.trees"]; 1042 1043 utility.SetEnvVariable ("BRANCH_LENGTH_STENCIL", stencils[^"terms.genetic_code.nonsynonymous"]); 1044 fit_results[^"terms.fit.nonsynonymous_trees"] = (estimators.ExtractMLEs(fit_results[^"terms.likelihood_function"], fit_results[^"terms.model"]))[^"terms.fit.trees"]; 1045 1046 utility.SetEnvVariable ("BRANCH_LENGTH_STENCIL", None); 1047 1048 return fit_results; 1049} 1050 1051 1052/** 1053 * @name estimators.FitCodonModel 1054 * @param {DataFilter} codon_data 1055 * @param {Tree} tree 1056 * @param {String} genetic_code 1057 * @param {Dictionary} option 1058 * @param {Dictionary} initial_values 1059 * @returns MGREV results 1060 */ 1061lfunction estimators.FitCodonModel(codon_data, tree, generator, genetic_code, option, initial_values) { 1062 1063 1064 1065 //TODO: Where is data_filter being set? 1066 if (Type(data_filter) == "String") { 1067 return estimators.FitCodonModel({ 1068 { 1069 codon_data__ 1070 } 1071 }, { 1072 "0": tree 1073 }, 1074 generator, 1075 genetic_code, 1076 option, 1077 initial_values) 1078 } 1079 1080 components = utility.Array1D(codon_data); 1081 1082 lf_components = { 1083 2 * components, 1084 1 1085 }; 1086 1087 1088 1089 for (i = 0; i < components; i += 1) { 1090 GetDataInfo(fi, ^ (codon_data[i]), "PARAMETERS"); 1091 DataSetFilter * ("filter_" + i) = CreateFilter( ^ (codon_data[i]), 3, '', '', fi["EXCLUSIONS"]); 1092 // need to do this for global references 1093 lf_components[2 * i] = "filter_" + i; 1094 } 1095 1096 1097 name_space = & model_MGREV; 1098 1099 mg_rev = model.generic.DefineModel(generator, 1100 name_space, { 1101 "0": parameters.Quote(option[utility.getGlobalValue("terms.run_options.model_type")]), 1102 "1": genetic_code 1103 }, 1104 codon_data, 1105 None); 1106 1107 1108 //utility.ToggleEnvVariable("VERBOSITY_LEVEL", 10); 1109 1110 df = 0; 1111 model_assignment = { 1112 "default": mg_rev 1113 }; 1114 rules = None; 1115 model_id_to_object = { 1116 name_space: mg_rev 1117 }; 1118 1119 for (i = 0; i < components; i += 1) { 1120 lf_components[2 * i + 1] = "tree_" + i; 1121 model.ApplyModelToTree(Eval("&`lf_components[2*i + 1]`"), tree[i], model_assignment, None); 1122 } 1123 1124 1125 partition_omega = {}; 1126 1127 if (option[utility.getGlobalValue("terms.run_options.model_type")] == utility.getGlobalValue("terms.local") && Type(option[utility.getGlobalValue("terms.run_options.partitioned_omega")]) == "AssociativeList") { 1128 /** 1129 Assumes that option["partitioned-omega"] is a dictionary where each partition has 1130 an entry (0-index based), which itself is a dictionary of the form: "branch-name" : "branch-set" 1131 */ 1132 utility.ForEach(option[utility.getGlobalValue("terms.run_options.partitioned_omega")], "_value_", "utility.AddToSet(`&partition_omega`,utility.UniqueValues(_value_))"); 1133 } 1134 1135 1136 if (Abs(partition_omega)) { 1137 1138 /** 1139 declare the global ratios for each branch set 1140 and add them to the model parameter set 1141 */ 1142 1143 1144 1145 new_globals = {}; 1146 utility.ForEachPair(partition_omega, "_key_", "_value_", 1147 '`&new_globals` [_key_] = (`&name_space` + ".omega_" + Abs (`&new_globals`)); model.generic.AddGlobal (`&mg_rev`, `&new_globals` [_key_] , (utility.getGlobalValue("terms.parameters.omega_ratio")) + " for *" + _key_ + "*")'); 1148 parameters.DeclareGlobal(new_globals, None); 1149 1150 1151 /** 1152 now replicate the local constraint for individual branches 1153 */ 1154 1155 1156 alpha = model.generic.GetLocalParameter(mg_rev, utility.getGlobalValue("terms.parameters.synonymous_rate")); 1157 beta = model.generic.GetLocalParameter(mg_rev, utility.getGlobalValue("terms.parameters.nonsynonymous_rate")); 1158 io.CheckAssertion("None!=`&alpha` && None!=`&beta`", "Could not find expected local synonymous and non-synonymous rate parameters in \`estimators.FitMGREV\`"); 1159 1160 SetParameter (DEFER_CONSTRAINT_APPLICATION, 1, 0); 1161 1162 apply_constraint: = component_tree + "." + node_name + "." + beta + ":=" + component_tree + "." + node_name + "." + alpha + "*" + new_globals[branch_map[node_name]]; 1163 1164 for (i = 0; i < components; i += 1) { 1165 component_tree = lf_components[2 * i + 1]; 1166 ClearConstraints( * component_tree); 1167 branch_map = (option[utility.getGlobalValue("terms.run_options.partitioned_omega")])[i]; 1168 1169 component_branches = BranchName( * component_tree, -1); 1170 for (j = 0; j < Columns(component_branches) - 1; j += 1) { 1171 /** 1172 -1 in the upper bound because we don't want to count the root node 1173 */ 1174 1175 node_name = (component_branches[j]); 1176 ExecuteCommands(apply_constraint); 1177 } 1178 } 1179 SetParameter (DEFER_CONSTRAINT_APPLICATION, 0, 0); 1180 } else {} 1181 1182 1183 LikelihoodFunction likelihoodFunction = (lf_components); 1184 GetString (params, likelihoodFunction,-1); 1185 1186 utility.ToggleEnvVariable ("PARAMETER_GROUPING", {"0" : params["Global Independent"]}); 1187 1188 if (utility.Has (option,utility.getGlobalValue("terms.run_options.apply_user_constraints"),"String")) { 1189 df += Call (option[utility.getGlobalValue("terms.run_options.apply_user_constraints")], &likelihoodFunction, lf_components, codon_data, tree, model_map, initial_values, model_id_to_object); 1190 } 1191 1192 if (Type(initial_values) == "AssociativeList") { 1193 utility.ToggleEnvVariable("USE_LAST_RESULTS", 1); 1194 df += estimators.ApplyExistingEstimates("`&likelihoodFunction`", model_id_to_object, initial_values, option[utility.getGlobalValue("terms.run_options.proportional_branch_length_scaler")]); 1195 } 1196 1197 /*GetString (res, likelihoodFunction, -1); 1198 1199 utility.ForEach (res[utility.getGlobalValue ('terms.parameters.local_independent')], '_value_', ' 1200 parameters.SetRange (_value_,terms.range_clamp_locals); 1201 ');*/ 1202 1203 //io.SpoolLF ("likelihoodFunction", "/tmp/hyphy-spool", "cfit"); 1204 1205 //Export (lfe, likelihoodFunction); 1206 //console.log (lfe); 1207 1208 //utility.ToggleEnvVariable("VERBOSITY_LEVEL", 10); 1209 1210 if (utility.Has (option,utility.getGlobalValue("terms.run_options.optimization_settings"),"AssociativeList")) { 1211 Optimize(mles, likelihoodFunction, option[utility.getGlobalValue("terms.run_options.optimization_settings")]); 1212 } else { 1213 Optimize(mles, likelihoodFunction); 1214 } 1215 1216 if (Type(initial_values) == "AssociativeList") { 1217 utility.ToggleEnvVariable("USE_LAST_RESULTS", None); 1218 } 1219 1220 utility.ToggleEnvVariable ("PARAMETER_GROUPING", None); 1221 1222 results = estimators.ExtractMLEs( & likelihoodFunction, model_id_to_object); 1223 1224 1225 results[utility.getGlobalValue("terms.fit.log_likelihood")] = mles[1][0]; 1226 results[utility.getGlobalValue("terms.parameters")] = mles[1][1] + (mg_rev [utility.getGlobalValue("terms.parameters")]) [utility.getGlobalValue("terms.model.empirical")] + df; 1227 1228 1229 if (option[utility.getGlobalValue("terms.run_options.retain_model_object")]) { 1230 results[utility.getGlobalValue("terms.model")] = model_id_to_object; 1231 } 1232 1233 if (option[utility.getGlobalValue("terms.run_options.retain_lf_object")]) { 1234 results[utility.getGlobalValue("terms.likelihood_function")] = & likelihoodFunction; 1235 } else { 1236 DeleteObject(likelihoodFunction); 1237 } 1238 1239 if (option[utility.getGlobalValue("terms.run_options.retain_model_object")]) { 1240 results[utility.getGlobalValue("terms.model")] = model_id_to_object; 1241 } 1242 1243 return results; 1244} 1245 1246 1247 1248 1249 1250/** 1251 * @name estimators.FitMGREV 1252 * @param {DataFilter} codon_data 1253 * @param {Tree} tree 1254 * @param {String} genetic_code 1255 * @param {Dictionary} option 1256 * @param {Dictionary} initial_values 1257 * @returns MGREV results 1258 */ 1259lfunction estimators.FitMGREV(codon_data, tree, genetic_code, option, initial_values) { 1260 return estimators.FitCodonModel (codon_data, tree, "models.codon.MG_REV.ModelDescription", genetic_code, option, initial_values); 1261} 1262 1263/** 1264 * @name estimators.FitMGREV 1265 * @description compute the asymptotic (chi^2) p-value for the LRT 1266 * @param {Number} alternative log likelihood for the alternative (more general model) 1267 * @param {Number} Null log likelihood for the null 1268 * @param {Number} df degrees of freedom 1269 * @returns p-value 1270 */ 1271lfunction estimators.LRT (alternative, Null, df) { 1272 if (alternative > Null) { 1273 return 1-CChi2 (2*(alternative-Null), df); 1274 } 1275 return 1; 1276} 1277 1278/** 1279 * @name estimators.ComputeLF 1280 * @description compute the current value of the log likelihood function 1281 * @param {String} lfid name of the function 1282 * @returns log likelihood 1283 */ 1284lfunction estimators.ComputeLF (id) { 1285 LFCompute (^id,LF_START_COMPUTE); 1286 LFCompute (^id,logl); 1287 LFCompute (^id,LF_DONE_COMPUTE); 1288 return logl; 1289} 1290 1291/** 1292 * @name estimators.CreateInitialGrid 1293 * @description prepare a Dict object suitable for seeding initial LF values 1294 * @param {Dict} values : "parameter_id" -> {{initial values}} [row matrix], e.g. 1295 ... 1296 "busted.test.bsrel_mixture_aux_0": { 1297 {0.1, 0.25, 0.4, 0.55, 0.7, 0.85, 0.9} 1298 }, 1299 "busted.test.bsrel_mixture_aux_1": { 1300 {0.1, 0.25, 0.4, 0.55, 0.7, 0.85, 0.9} 1301 } 1302 ... 1303 1304 * @param {int} N how many points to sample 1305 * @param {Dict/null} init if not null, taken to be the initial template for variables 1306 i.e. random draws will be one index change from this vector 1307 * @returns {Dict} like in 1308 1309 { 1310 "0":{ 1311 "busted.test.bsrel_mixture_aux_0":{ 1312 "ID":"busted.test.bsrel_mixture_aux_0", 1313 "MLE":0.4 1314 }, 1315 "busted.test.bsrel_mixture_aux_1":{ 1316 "ID":"busted.test.bsrel_mixture_aux_1", 1317 "MLE":0.7 1318 }, 1319 "busted.test.omega1":{ 1320 "ID":"busted.test.omega1", 1321 "MLE":0.01 1322 }, 1323 "busted.test.omega2":{ 1324 "ID":"busted.test.omega2", 1325 "MLE":0.1 1326 }, 1327 "busted.test.omega3":{ 1328 "ID":"busted.test.omega3", 1329 "MLE":1.5 1330 } 1331 } 1332 .... 1333 */ 1334lfunction estimators.CreateInitialGrid (values, N, init) { 1335 result = {}; 1336 var_count = utility.Array1D (values); 1337 var_names = utility.Keys (values); 1338 var_dim = {var_count,1}; 1339 for (v = 0; v < var_count; v += 1) { 1340 var_dim [v] = utility.Array1D (values[var_names[v]]); 1341 } 1342 1343 if (null != init) { 1344 1345 toggle = Min (init[0], 0.5); 1346 1347 for (i = 0; i < N; i+=1) { 1348 entry = {}; 1349 for (v = 0; v < var_count; v += 1) { 1350 entry [var_names[v]] = {}; 1351 (entry [var_names[v]])[^"terms.id"] = var_names[v]; 1352 if (Random (0,1) < toggle) { 1353 (entry [var_names[v]])[^"terms.fit.MLE"] = (values[var_names[v]])[Random (0, var_dim[v])$1]; 1354 } else { 1355 (entry [var_names[v]])[^"terms.fit.MLE"] = (values[var_names[v]])[init[var_names[v]]]; 1356 } 1357 } 1358 result + entry; 1359 } 1360 } else { 1361 1362 1363 for (i = 0; i < N; i+=1) { 1364 entry = {}; 1365 for (v = 0; v < var_count; v += 1) { 1366 entry [var_names[v]] = {}; 1367 (entry [var_names[v]])[^"terms.id"] = var_names[v]; 1368 (entry [var_names[v]])[^"terms.fit.MLE"] = (values[var_names[v]])[Random (0, var_dim[v])$1]; 1369 } 1370 result + entry; 1371 } 1372 } 1373 1374 1375 return result; 1376} 1377 1378 1379/** 1380 * @name estimators.LHC 1381 * @description prepare a Dict object suitable for seeding initial LF values 1382 based on Latin Hypercube Sampling 1383 * @param {Dict} ranges : "parameter_id" -> range, i.e. { 1384 lower_bound: 0, 1385 upper_bound: 1 1386 }. 1387 1388 * @param {Number} samples : the # of samples to draw 1389 1390*/ 1391 1392lfunction estimators.LHC (ranges, samples) { 1393 1394 1395 result = {}; 1396 var_count = utility.Array1D (ranges); 1397 var_names = utility.Keys (ranges); 1398 var_def = {var_count,2}; 1399 var_samplers = {}; 1400 1401 for (v = 0; v < var_count; v += 1) { 1402 var_def [v][0] = (ranges[var_names[v]])[^"terms.lower_bound"]; 1403 var_def [v][1] = ((ranges[var_names[v]])[^"terms.upper_bound"] - var_def [v][0]) / (samples-1); 1404 var_samplers[v] = Random ({1,samples}["_MATRIX_ELEMENT_COLUMN_"], 0); 1405 } 1406 1407 1408 result = {}; 1409 1410 for (i = 0; i < samples; i+=1) { 1411 entry = {}; 1412 for (v = 0; v < var_count; v += 1) { 1413 entry [var_names[v]] = {}; 1414 (entry [var_names[v]])[^"terms.id"] = var_names[v]; 1415 (entry [var_names[v]])[^"terms.fit.MLE"] = var_def[v][0] + (var_samplers[v])[i] * var_def[v][1]; 1416 } 1417 result + entry; 1418 } 1419 1420 return result; 1421} 1422