1USE_MPI_CACHING = 1; 2PRINT_DIGITS = -1; 3 4 5function add_discrete_node (node_id, max_parents, sample_size, nlevels) 6{ 7 node = {}; 8 node["NodeID"] = node_id; 9 node["NodeType"] = 0; 10 node["MaxParents"] = max_parents; 11 node["PriorSize"] = sample_size; 12 node["NumLevels"] = nlevels; 13 return node; 14} 15 16function add_gaussian_node (node_id, max_parents, sample_size, mean, precision, scale) 17{ 18 node = {}; 19 node["NodeID"] = node_id; 20 node["NodeType"] = 1; 21 node["MaxParents"] = max_parents; 22 node["PriorSize"] = sample_size; 23 node["PriorMean"] = mean; 24 node["PriorPrecision"] = precision; 25 node["PriorScale"] = scale; 26 return node; 27} 28 29 30 31/* utility functions from ReadDelimitedFiles.bf */ 32function ReadCSVTable (fileName, haveHeader) 33{ 34 if (Abs(fileName) == 0) 35 { 36 fscanf (PROMPT_FOR_FILE, "Lines", inData); 37 } 38 else 39 { 40 fscanf (fileName, "Lines", inData); 41 } 42 if (haveHeader) 43 { 44 output = {}; 45 output[0] = splitOnRegExp (inData[0],"\\,"); 46 } 47 felMXString = ""; 48 felMXString * 256; 49 felMXString * "_tempMatrix={"; 50 for (lineID = haveHeader; lineID < Columns(inData); lineID = lineID + 1) 51 { 52 felMXString * ("{" + inData[lineID] + "}\n"); 53 } 54 felMXString * "}"; 55 felMXString * 0; 56 ExecuteCommands (felMXString); 57 felMXString = 0; 58 inData = 0; 59 if (haveHeader) 60 { 61 output[1] = _tempMatrix; 62 _tempMatrix = 0; 63 return output; 64 } 65 return _tempMatrix; 66} 67 68 69function splitOnRegExp (string, splitter) 70{ 71 matched = string || splitter; 72 splitBits = {}; 73 if (matched [0] < 0) 74 { 75 splitBits[0] = string; 76 } 77 else 78 { 79 mc = 0; 80 if (matched[0] == 0) 81 { 82 fromPos = matched[1]+1; 83 mc = 2; 84 } 85 else 86 { 87 fromPos = 0; 88 toPos = 0; 89 } 90 for (; mc < Rows (matched); mc = mc+2) 91 { 92 toPos = matched[mc]-1; 93 splitBits [Abs(splitBits)] = string[fromPos][toPos]; 94 fromPos = matched[mc+1]+1; 95 } 96 splitBits [Abs(splitBits)] = string[fromPos][Abs(string)-1]; 97 } 98 return splitBits; 99} 100 101 102/* a wrapper around ReadCSVTable */ 103function import_data (inData, hasHeader) 104{ 105 timer0 = Time(0); 106 file_input = ReadCSVTable (inData, hasHeader); 107 108 bgm_data_matrix = {{}}; 109 names = {{}}; 110 num_nodes = 0; 111 112 if (hasHeader) 113 { 114 names = file_input["0"]; 115 bgm_data_matrix = file_input["1"]; 116 117 fprintf (stdout, "Read ", Rows(bgm_data_matrix), " cases from file.\n"); 118 119 num_nodes = Columns(bgm_data_matrix); 120 121 if (Abs(file_input["0"]) != num_nodes) 122 { 123 fprintf (stdout, "ERROR! Number of items in header does not match the number of items in the data matrix."); 124 return 0; 125 } 126 127 fprintf (stdout, "Detected ", num_nodes, " variables.\n"); 128 } 129 else 130 { 131 bgm_data_matrix = file_input; 132 133 fprintf (stdout, "Read ", Rows(bgm_data_matrix), " cases from file.\n"); 134 135 num_nodes = Columns(bgm_data_matrix); 136 names = {num_nodes, 1}; 137 138 for (i = 0; i < num_nodes; i = i+1) 139 { 140 names[i] = i; 141 } 142 143 fprintf (stdout, "Detected ", num_nodes, " variables.\n"); 144 } 145 146 return bgm_data_matrix; 147} 148 149 150 151function import_cache (filename, cache_name) 152{ 153 fscanf (filename, "Raw", cacheStr); 154 ExecuteCommands(cache_name+" = "+cacheStr+";"); 155 return 0; 156} 157 158 159 160function attach_cache (_bgm, cache) 161{ 162 ExecuteCommands ("SetParameter("+_bgm+", BGM_SCORE_CACHE, cache);"); 163 return 0; 164} 165 166 167 168 169/* ____________________________________________________________ */ 170/* accessor functions */ 171function setStructure (_bgm, graph_matrix) 172{ 173 ExecuteCommands("SetParameter ("+_bgm+", BGM_GRAPH_MATRIX, graph_matrix);"); 174} 175 176function setOrder (_bgm, order_matrix) 177{ 178 if (Rows(order_matrix) > 1) 179 { 180 if (Columns(order_matrix) == 1) 181 { 182 t_order_matrix = Transpose(order_matrix); 183 ExecuteCommands("SetParameter ("+_bgm+", BGM_NODE_ORDER, t_order_matrix);"); 184 } 185 else 186 { 187 fprintf (stdout, "Warning: expecting row vector matrix, received non-vector matrix"); 188 fprintf (stdout, " with dimensions ", Rows(order_matrix), " x ", Columns(order_matrix), "\n"); 189 fprintf (stdout, "Node order not set!\n"); 190 } 191 } 192 else 193 { 194 ExecuteCommands ("SetParameter ("+_bgm+", BGM_NODE_ORDER, order_matrix);"); 195 } 196} 197 198 199function setConstraints (_bgm, constraint_matrix) 200{ 201 ExecuteCommands("SetParameter ("+_bgm+", BGM_CONSTRAINT_MATRIX, constraint_matrix);"); 202 return 0; 203} 204 205 206/* ____________________________________________________________ */ 207/* Assign data matrix to _BayesianGraphicalModel object */ 208function attach_data (_bgm, data, impute_max, impute_burn, impute_samp) 209{ 210 BGM_IMPUTE_MAXSTEPS = impute_max$1; 211 BGM_IMPUTE_BURNIN = impute_burn$1; 212 BGM_IMPUTE_SAMPLES = impute_samp$1; 213 214 ExecuteCommands("SetParameter ("+_bgm+", BGM_DATA_MATRIX, data);"); 215 return 0; 216} 217 218 219 220/* 221 Structural (graph) MCMC by Metropolis-Hastings 222 Returns matrix object containing chain trace, edge 223 marginal posterior probabilities, and best graph as 224 adjacency matrix. 225 226 rand_tolerance = maximum number of failed steps in graph randomization 227 to tolerate 228 229 prob_swap = probability of reversing an edge, instead of adding or deleting an edge 230 231 with_order = a vector containing node ordering to constrain graph MCMC 232 set to 0 to have unconstrained chain sample 233*/ 234 235BGM_MCMC_MAXFAILS = 100; 236BGM_MCMC_PROBSWAP = 0.1; 237 238function graph_MCMC (_bgm, duration, burnin, num_samples, with_order=0) 239{ 240 if (Rows(with_order) * Columns(with_order) > 0) 241 { 242 /* fixed node order */ 243 ExecuteCommands("setOrder ("+_bgm+", with_order);"); 244 BGM_OPTIMIZATION_METHOD = 2; 245 } 246 else 247 { 248 /* shuffle node order */ 249 BGM_OPTIMIZATION_METHOD = 3; 250 } 251 252 BGM_MCMC_MAXSTEPS = duration; 253 BGM_MCMC_BURNIN = burnin; 254 BGM_MCMC_SAMPLES = num_samples; 255 256 ExecuteCommands("Optimize(res, "+_bgm+");"); 257 258 return res; 259} 260 261 262/* 263 Order (node precedence permutation) MCMC by Metropolis-Hastings 264*/ 265function order_MCMC (_bgm, duration, burnin, num_samples) 266{ 267 BGM_OPTIMIZATION_METHOD = 4; 268 269 BGM_MCMC_MAXSTEPS = duration; 270 BGM_MCMC_BURNIN = burnin; 271 BGM_MCMC_SAMPLES = num_samples; 272 273 ExecuteCommands("Optimize(res, "+_bgm+");"); 274 275 return res; 276} 277 278 279 280 281 282function display_MCMC_chain (res) 283{ 284 if (Rows(res)*Columns(res) == 0) 285 { 286 fprintf (stdout, "ERROR: Cannot display MCMC chain for empty matrix\n"); 287 return 1; 288 } 289 290 pp_trace = res[-1][0]; 291 min_trace = pp_trace[0]; 292 max_trace = pp_trace[0]; 293 294 /* locate min/max and end of trace */ 295 for (k = 0; k < Rows(pp_trace); k = k+1) 296 { 297 if (pp_trace[k] == 0) 298 { 299 break; 300 } 301 if (pp_trace[k] < min_trace) 302 { 303 min_trace = pp_trace[k]; 304 } 305 if (pp_trace[k] > max_trace) 306 { 307 max_trace = pp_trace[k]; 308 } 309 } 310 k = k-1; 311 pp_trace = pp_trace[{{0,0}}][{{k-1,0}}]; 312 313 314 columnHeaders = {{"MCMC chain","sample;1;2;3;4;5;6;7;8;9"}}; 315 316 OpenWindow (CHARTWINDOW,{{"Posterior probability"} 317 {"columnHeaders"} 318 {"pp_trace"} 319 {"Step Plot"} 320 {"Index"} 321 {"MCMC chain"} 322 {"chain sample step"} 323 {"posterior prob."} 324 {""} 325 {"0"} 326 {""} 327 {"0;0"} 328 {"10;1.309;0.785398"} 329 {"Times:12:0;Times:10:0;Times:12:2"} 330 {"0;0;13816530;16777215;0;0;6579300;11842740;13158600;14474460;0;3947580;16777215;15670812;6845928;16771158;2984993;9199669;7018159;1460610;16748822;11184810;14173291"} 331 {"16,"+min_trace+","+max_trace} 332 }, 333 "405;462;105;100"); 334 335 return 0; 336} 337 338 339function get_MCMC_graph (res, num_nodes, mode) 340{ 341 /* mode = -1 : best_graph 342 mode = 0 : last_graph 343 0 < mode <= 1 : marginal posterior graph with threshold = mode (e.g. 0.9) 344 */ 345 graph = {num_nodes, num_nodes}; 346 347 if (mode > 0) 348 { 349 for (row = 0; row < num_nodes * num_nodes; row = row+1) 350 { 351 if (res[row][1] >= mode) 352 { 353 graph[row $ num_nodes][row % num_nodes] = 1; 354 } 355 } 356 } 357 else 358 { 359 for (row = 0; row < num_nodes; row = row+1) 360 { 361 for (col = 0; col < num_nodes; col = col+1) 362 { 363 graph[row][col] = res[row*num_nodes+col][mode+3]; 364 } 365 } 366 } 367 368 return graph; 369} 370 371 372function write_edgelist (filename,res,num_nodes,directed) 373{ 374 fprintf (filename, CLEAR_FILE, KEEP_OPEN); 375 if (directed) 376 { 377 for (row = 0; row < num_nodes; row = row+1) 378 { 379 for (col = 0; col < num_nodes; col = col+1) 380 { 381 fprintf (filename, names[row], ",", names[col], ",", res[row*num_nodes+col][1], "\n"); 382 } 383 } 384 } 385 else 386 { 387 for (row = 0; row < num_nodes-1; row = row+1) 388 { 389 for (col = row+1; col < num_nodes; col = col+1) 390 { 391 fprintf (filename, names[row], ",", names[col], ",", res[row*num_nodes+col][1] + res[col*num_nodes+row][1], "\n"); 392 } 393 } 394 } 395 fprintf (filename, CLOSE_FILE); 396 return 0; 397} 398 399 400function mcmc_graph_to_dotfile (filename, threshold, res, nodes) 401{ 402 fprintf (filename, CLEAR_FILE); 403 fprintf (filename, "digraph foo\n{\n"); 404 fprintf (filename, "\tnode [fontname=\"Helvetica\" style=\"filled\" fillcolor=\"white\"];\n"); 405 fprintf (filename, "\tedge [labelfontname=\"Helvetica\" labelangle=30 labeldistance=2];\n"); 406 407 for (_n = 0; _n < Abs(nodes); _n+=1) { 408 fprintf (filename, "\t", (nodes[_n])["NodeID"]); 409 if ((nodes[_n])["NodeType"]==0) { 410 fprintf (filename, " [shape=\"Msquare\"];\n"); 411 } else { 412 fprintf (filename, " [shape=\"circle\"];\n"); 413 } 414 } 415 416 417 // sum edge posteriors in both directions between nodes X and Y, 418 // and assign direction to the greater value 419 for (row = 0; row < num_nodes-1; row = row+1) { 420 for (col = row+1; col < num_nodes; col = col+1) { 421 xy = res[row*num_nodes+col][1]; 422 yx = res[col*num_nodes+row][1]; 423 if (xy+yx > threshold) { 424 /* 425 This is really annoying - order MCMC reports edge marginal matrix with rows = child 426 whereas graph MCMC reports rows = parent 427 */ 428 if ( xy > yx ) { 429 fprintf (filename, "\t", (nodes[row])["NodeID"], "->", (nodes[col])["NodeID"], ";\n"); 430 } else { 431 fprintf (filename, "\t", (nodes[col])["NodeID"], "->", (nodes[row])["NodeID"], ";\n"); 432 } 433 } 434 } 435 } 436 437 fprintf (filename, "}\n"); 438 return 0; 439} 440 441 442/* argument must be string identifier of BGM object */ 443function get_network_parameters (_bgm) 444{ 445 ExecuteCommands("GetString (res, "+_bgm+", 1);"); 446 ExecuteCommands(res); 447 /* returns string identifier to associative array */ 448 ExecuteCommands("params="+_bgm+"_export;"); 449 return params; 450} 451 452 453function get_node_score_cache (_bgm) 454{ 455 ExecuteCommands("GetString (res, "+_bgm+", 0);"); 456 return res; 457} 458 459 460/* 461function getStructure (_bgm) 462{ 463 ExecuteCommands("GetInformation (s, "+_bgm+", 0);"); 464 return s; 465} 466 467function getNodeOrder (_bgm) 468{ 469 ExecuteCommands("GetInformation (s, "+_bgm+", 1);"); 470 return s; 471} 472 473*/ 474 475 476 477 478/* 479 Simulation of data based on the inferred network 480 structure and parameters. 481 mode = 0 (local) : for each case, instantiate parameters de novo. 482 Better for assessing uncertainty. 483 mode = 1 (global) : instantiate all parameters once. 484 Assuming known network. 485*/ 486function instantiate_CPDFs (params) 487{ 488 node_names = Rows(params); 489 490 /* instantiate network parameters from conditional posterior distribution functions */ 491 for (i = 0; i < Abs(params); i = i + 1) { 492 /* stores instantiations */ 493 ExecuteCommands("(params[\""+node_names[i]+"\"])[\"Parameters\"] = {};"); 494 495 /* number of parent combinations */ 496 //ExecuteCommands("npac = Columns((params[\""+node_names[i]+"\"])[\"CPDFs\"]);"); 497 ExecuteCommands("npac = (params[\""+node_names[i]+"\"])[\"NParentCombs\"];"); // safe version 498 499 for (pa = 0; pa < npac; pa = pa+1) { 500 ExecuteCommands("_p = " + ((params[node_names[i]])["CPDFs"])[pa] + ";"); 501 ExecuteCommands("((params[\""+node_names[i]+"\"])[\"Parameters\"])[\""+pa+"\"] = "+_p+";"); 502 } 503 504 //ExecuteCommands("((params[\""+node_names[i]+"\"])[\"Levels\"] = Columns( ((params[\""+node_names[i]+"\"])[\"Parameters\"])[0] ));"); 505 } 506 return 0; 507} 508 509 510/* 511 Return a parameter vector for conditional Gaussian (CG) node given 512 hyperparameters passed as arguments. 513*/ 514function cg_params (mean_vec, rho, phi, tau) { 515 ExecuteCommands("sigma = Random({{"+phi+"}}, {\"PDF\":\"InverseWishart\", \"ARG0\":{{"+rho+"}} });"); 516 ExecuteCommands("em = Random("+mean_vec+", {\"PDF\":\"Gaussian\", \"ARG0\":(Inverse("+tau+") * "+sigma[0]+") } );"); 517 return ({"EM":em, "SIGMA":sigma}); 518} 519 520 521 522function simulate_data (params, num_cases) 523{ 524 // prepare matrix to store simulated data 525 result = {num_cases, Abs(params)}; 526 527 node_names = Rows(params); 528 if ( Columns(Rows((params[node_names[0]])["Parameters"])) == 0 ) 529 { 530 /* parameters have not been instantiated yet */ 531 instantiate_CPDFs(params); 532 } 533 534 535 // initialize State variables and generate root states 536 for (case = 0; case < num_cases; case = case+1) { 537 538 for (i = 0; i < Abs(params); i = i + 1) { 539 // set to String as a placeholder 540 (params[node_names[i]])["State"] = ""; 541 542 if ( Type((params[node_names[i]])["Parents"]) == "AssociativeList" ) { 543 // if condition is true then this is a root node (no parents) 544 if ( (params[node_names[i]])["NodeType"] == 0 ) { 545 // discrete node, parameters define conditional probability table 546 urn = Random(0,1); 547 cpt = ((params[node_names[i]])["Parameters"])[0]; 548 r_i = Columns(cpt); 549 for (k = 0; k < r_i; k = k+1) 550 { 551 if ( urn <= cpt[k] ) 552 { 553 (params[node_names[i]])["State"] = k; 554 break; 555 } 556 urn = urn - cpt[k]; 557 } 558 } else { 559 // conditional Gaussian node, parameter defines intercept 560 em = (((params[node_names[i]])["Parameters"])[0])["EM"]; 561 sigma = (((params[node_names[i]])["Parameters"])[0])["SIGMA"]; 562 (params[node_names[i]])["State"] = (Random(em, {"PDF":"Gaussian", "ARG0":sigma}))[0]; 563 } 564 } 565 } 566 567 while (1) 568 { 569 all_done = 1; 570 571 /* loop until parameters are instantiated for all nodes */ 572 for (i = 0; i < Abs(params); i = i+1) 573 { 574 if (Type(params[node_names[i]])["State"] == "String") 575 { 576 // Type String indicates no value - replace placeholder with NoneType when it becomes available 577 578 all_done = 0; 579 ok_to_go = 1; 580 581 parents = (params[node_names[i]])["Parents"]; 582 num_parent_combos = 1; 583 pa_index = 0; 584 585 for (p = 0; p < Abs(Rows(parents)); p = p+1) 586 { 587 pid = parents[p]; 588 if ( Type(params[pid])["State"] == "String" ) 589 { 590 // parents not resolved, skipping 591 ok_to_go = 0; 592 break; 593 } 594 595 // compute parental index for discrete parents 596 if ( (params[pid])["NodeType"] == 0 ) { 597 pa_index = pa_index + (params[pid])["State"] * num_parent_combos; 598 num_parent_combos = num_parent_combos * (params[pid])["Levels"]; 599 } 600 } 601 602 603 if (ok_to_go) 604 { 605 // instantiate this node's parameters 606 if ( (params[node_names[i]])["NodeType"] == 0 ) { 607 urn = Random(0,1); 608 cpt = ((params[node_names[i]])["Parameters"])[pa_index]; 609 r_i = Columns(cpt); 610 for (k = 0; k < r_i; k = k+1) { 611 if ( urn <= cpt[k] ) { 612 (params[node_names[i]])["State"] = k; 613 break; 614 } 615 urn = urn - cpt[k]; 616 } 617 } else { 618 em = ( ((params[node_names[i]])["Parameters"])[pa_index] )["EM"]; 619 sigma = ( ((params[node_names[i]])["Parameters"])[pa_index] )["SIGMA"]; 620 zvec = {Columns(em), 1}; 621 zvec[0] = 1; 622 623 // get states of continuous parents 624 cpar = 0; 625 for (p = 0; p < Abs(Rows(parents)); p += 1) { 626 pid = parents[p]; 627 if ( (params[pid])["NodeType"] == 1 ) { 628 zvec[cpar+1] = (params[pid])["State"]; 629 cpar += 1; 630 } 631 } 632 633 // conditional mean 634 cond_mean = em * zvec; 635 (params[node_names[i]])["State"] = (Random(cond_mean, {"PDF":"Gaussian", "ARG0":sigma}))[0]; 636 } 637 } 638 } 639 } 640 /* end for loop */ 641 642 if (all_done) break; 643 } 644 /* end while */ 645 646 /* add case to result */ 647 for (i = 0; i < Abs(params); i = i+1) { 648 result[case][i] = (params[node_names[i]])["State"]; 649 } 650 } 651 652 return result; 653} 654 655 656/* 657 Example: 658 import_xmlbif("/Users/apoon/svn/hyphy/HBL/art/BGM/alarm/alarm.xml", "Alarm"); 659*/ 660function import_xmlbif (filename, newname) 661{ 662 ExecuteCommands(newname+"={};"); 663 664 fscanf (filename, "Raw", input); 665 666 var_tags = input||"<VARIABLE"; 667 if (var_tags[0] < 0) 668 { 669 fprintf (stdout, "ERROR: <VARIABLE> tag absent from XML, exiting.."); 670 return 1; 671 } 672 673 ntags = Rows(var_tags)$2; 674 675 676 for (tag = 0; tag < ntags; tag = tag+1) 677 { 678 /* 679 search for <NAME> tag - note that we use an arbitrary character limit (1000) 680 for the last entry because if we use the rest of the XML file, it causes the 681 regular expression search to fail! - afyp, October 26, 2011 682 */ 683 start_char = var_tags[tag*2+1]; 684 if (tag == ntags-1) { end_char = start_char+1000; } 685 else { end_char = var_tags[(tag+1)*2]; } 686 substr = input[start_char][end_char]; 687 688 /* create node */ 689 name_tag = substr||"<NAME>.+</NAME>"; 690 node_name = substr[name_tag[0]+6][name_tag[1]-7]; 691 692 693 ExecuteCommands(newname+"[\""+node_name+"\"]= {};"); 694 695 outcome_tags = substr||"<OUTCOME>"; 696 ExecuteCommands("("+newname+"[\""+node_name+"\"])[\"Levels\"]= "+Rows(outcome_tags)$2+";"); 697 } 698 699 700 def_tags = input||"<DEFINITION>"; 701 if (def_tags[0] < 0) 702 { 703 fprintf (stdout, "ERROR: <DEFINITION> tag absent from XML, exiting.."); 704 return 1; 705 } 706 707 ntags = Rows(def_tags)$2; 708 for (tag = 0; tag < ntags; tag = tag+1) 709 { 710 /* parse definition tags */ 711 start_char = def_tags[tag*2+1]; 712 if (tag == ntags-1) { end_char = Abs(input); } 713 else { end_char = def_tags[(tag+1)*2]; } 714 substr = input[start_char][end_char]; 715 716 /* start a new node */ 717 for_tag = substr||"<FOR>.+</FOR>"; 718 node_name = substr[for_tag[0]+5][for_tag[1]-6]; 719 720 /* assign parents */ 721 exec_str = ""; 722 exec_str * 256; 723 exec_str * "("; 724 exec_str * newname; 725 exec_str * "[\""; 726 exec_str * node_name; 727 exec_str * "\"])[\"Parents\"]={"; 728 given_tags = substr||"<?GIVEN>"; 729 if (given_tags[0] >= 0) 730 { 731 for (gt = 1; gt < Rows(given_tags); gt = gt+4) 732 { 733 exec_str * "{\""; 734 exec_str * substr[given_tags[gt]+1][given_tags[gt+1]-3]; 735 exec_str * "\"}"; 736 if (gt < Rows(given_tags)-4) { exec_str * ","; } 737 } 738 } 739 exec_str * "};"; 740 exec_str * 0; 741 ExecuteCommands(exec_str); 742 743 744 /* assign conditional probability table - child state cycles fastest, then parents */ 745 table_tag = substr||"<TABLE>.+</TABLE>"; 746 table_str = substr[table_tag[0]+7][table_tag[1]-8]; 747 prob_tags = table_str||"[01]\.[0-9]+"; 748 749 n_parent_combos = 1; 750 ExecuteCommands("parents = ("+newname+"[\""+node_name+"\"])[\"Parents\"];"); 751 for (par = 0; par < Abs(Rows(parents)); par=par+1) 752 { 753 ExecuteCommands("n_parent_combos = n_parent_combos * ("+newname+"[\""+parents[par]+"\"])[\"Levels\"];"); 754 } 755 ExecuteCommands("n_levels = ("+newname+"[node_name])[\"Levels\"];"); 756 757 ExecuteCommands("("+newname+"[\""+node_name+"\"])[\"Parameters\"]= {};"); 758 759 for (pa = 0; pa < n_parent_combos; pa = pa+1) 760 { 761 ExecuteCommands("(("+newname+"[\""+node_name+"\"])[\"Parameters\"])[\""+pa+"\"]={1,n_levels};"); 762 for (lev = 0; lev < n_levels; lev=lev+1) 763 { 764 foo = lev * n_parent_combos + pa; 765 /* fprintf (stdout, lev, ",", pa, ",", table_str[prob_tags[foo*2]][prob_tags[foo*2+1]], "\n"); */ 766 ExecuteCommands("((("+newname+"[\""+node_name+"\"])[\"Parameters\"])[\""+pa+"\"])["+lev+"]="+table_str[prob_tags[foo*2]][prob_tags[foo*2+1]]+";"); 767 } 768 } 769 770 771 } 772 773 return 0; 774} 775 776 777 778function list2adjmat (alist) { 779/* 780 convert associative list returned by import_xmlbif into an adjacency matrix 781*/ 782 num_nodes = Abs(alist); 783 res = {num_nodes, num_nodes}; 784 node_names = Rows(alist); 785 name2index = {}; 786 787 // for indexing into adjacency matrix 788 for (node = 0; node < num_nodes; node += 1) { 789 name2index[node_names[node]] = node; 790 } 791 792 for (child = 0; child < num_nodes; child += 1) { 793 parents = (alist[node_names[child]])["Parents"]; 794 if (Type(parents) == "Matrix") { 795 for (par = 0; par < Rows(parents); par += 1) { 796 parent = name2index[parents[par]]; 797 res[parent][child] = 1; 798 } 799 } 800 } 801 802 return res; 803} 804 805 806function check_edgelist (results, adjmat, cutoff) { 807 // extract edge marginal posteriors vector from results matrix (in column 1) 808 edgep = results[-1][1]; 809 num_nodes = Rows(adjmat); 810 true_pos = 0; 811 false_pos = 0; 812 true_neg = 0; 813 false_neg = 0; 814 815 for (parent = 0; parent < (num_nodes-1); parent += 1) { 816 for (child = (parent+1); child < num_nodes; child += 1) { 817 x = edgep[parent * num_nodes + child] + edgep[child * num_nodes + parent]; 818 819 if (adjmat[parent][child] > 0 || adjmat[child][parent] > 0) { 820 if ( x > cutoff ) { 821 true_pos += 1; 822 } else { 823 false_neg += 1; 824 } 825 } else { 826 if ( x > cutoff ) { 827 false_pos += 1; 828 } else { 829 true_neg += 1; 830 } 831 } 832 } 833 } 834 835 result = {4,1}; /* TP, FN, FP, TN */ 836 result[0] = true_pos; 837 result[1] = false_neg; 838 result[2] = false_pos; 839 result[3] = true_neg; 840 841 return (result); 842} 843 844 845 846