We start, as normal, by creating a Gobnilp
object to work with:
from pygobnilp.gobnilp import Gobnilp
m = Gobnilp()
Using license file /home/james/gurobi.lic Academic license - for non-commercial use only Changed value of parameter PreCrush to 1 Prev: 0 Min: 0 Max: 1 Default: 0 Changed value of parameter CutPasses to 100000 Prev: -1 Min: -1 Max: 2000000000 Default: -1 Changed value of parameter GomoryPasses to 100000 Prev: -1 Min: -1 Max: 2000000000 Default: -1 Changed value of parameter MIPFocus to 2 Prev: 0 Min: 0 Max: 3 Default: 0 Changed value of parameter ZeroHalfCuts to 2 Prev: -1 Min: -1 Max: 2 Default: -1 Changed value of parameter MIPGap to 0.0 Prev: 0.0001 Min: 0.0 Max: inf Default: 0.0001 Changed value of parameter MIPGapAbs to 0.0 Prev: 1e-10 Min: 0.0 Max: inf Default: 1e-10
Now we will create and store 'local scores' from the built-in alarm_10000.dat
dataset:
m.learn('alarm_10000.dat',end='local scores')
Note that, since we have stated that learning should end
after local scores have been computed and stored, a BN is not learned. The computation of these local scores is not immediate, but should not take too long. Since no local score function was specified and no data type specified either, the data will be assumed to be discrete and a BDeu score will be used.
The local scores are saved in a dictionary which is the local_scores
attribute of the Gobnilp
object. The keys of this dictionary are the BN variables (the names of these variables are given on the first line of alarm_10000.dat
).
print(m.local_scores.keys())
dict_keys(['ANAPHYLAXIS', 'ARTCO2', 'BP', 'CATECHOL', 'CO', 'CVP', 'DISCONNECT', 'ERRCAUTER', 'ERRLOWOUTPUT', 'EXPCO2', 'FIO2', 'HISTORY', 'HR', 'HRBP', 'HREKG', 'HRSAT', 'HYPOVOLEMIA', 'INSUFFANESTH', 'INTUBATION', 'KINKEDTUBE', 'LVEDVOLUME', 'LVFAILURE', 'MINVOL', 'MINVOLSET', 'PAP', 'PCWP', 'PRESS', 'PULMEMBOLUS', 'PVSAT', 'SAO2', 'SHUNT', 'STROKEVOLUME', 'TPR', 'VENTALV', 'VENTLUNG', 'VENTMACH', 'VENTTUBE'])
The value associated with each BN variable is a dictionary mapping potential parent sets for that variable to a local score - which states, intuitively, how good it would be for that parent set to be selected. We find that there is variation in the number of parent sets stored:
for child, scored_parentsets in m.local_scores.items():
print(child, len(scored_parentsets))
ANAPHYLAXIS 5 ARTCO2 473 BP 266 CATECHOL 383 CO 469 CVP 32 DISCONNECT 127 ERRCAUTER 22 ERRLOWOUTPUT 14 EXPCO2 176 FIO2 18 HISTORY 36 HR 154 HRBP 135 HREKG 122 HRSAT 122 HYPOVOLEMIA 57 INSUFFANESTH 1 INTUBATION 211 KINKEDTUBE 41 LVEDVOLUME 64 LVFAILURE 84 MINVOL 371 MINVOLSET 96 PAP 3 PCWP 43 PRESS 212 PULMEMBOLUS 10 PVSAT 440 SAO2 471 SHUNT 162 STROKEVOLUME 109 TPR 173 VENTALV 472 VENTLUNG 448 VENTMACH 134 VENTTUBE 317
ANAPHYLAXIS has a mere 5 parent sets whereas ARTCO2 has 473! This is because Gobnilp
is using pruning. If a variable $X$ has potential parent set $S$ which is a superset of another potential parent set $T$ which has a higher score then $S$ is 'pruned' away. This is because, assuming we do not have any special constraints on allowed BNs, we can never have an optimal BN where $S$ is a parent set, since we could always replace $S$ with $T$ and get a higher scoring BN without creating a cycle, since we would remove arrows.
It just so happens that more parent sets were pruned for ANAPHYLAXIS than ARTCO2. Let's look at the scored parent sets for ANAPHYLAXIS:
for parentset, score in m.local_scores['ANAPHYLAXIS'].items():
print(parentset, score)
frozenset() -578.587070487134 frozenset({'BP'}) -523.5582005530596 frozenset({'CATECHOL'}) -577.4536694322014 frozenset({'TPR'}) -471.48556937259855 frozenset({'HYPOVOLEMIA', 'BP'}) -518.1711042499956
The empty set has the worst score of the 5 - if any of the non-empty parent set had a worse score it would be pruned.
Each local score will be eventually be associated with a Gurobi MIP variable called a family variable (since it is associated with a 'child' variable and a parent set). The fewer of these family variables the better since the MIP will be smaller and thus easier for Gurobi to solve. If we have prior knowledge then this can
To provide an example of how prior knowledge can help, let's see how often variables occur in the potential parent sets for ARTCO2:
from collections import Counter
d = Counter()
for parent_set in m.local_scores['ARTCO2']:
for v in parent_set:
d[v] += 1
print(d)
Counter({'EXPCO2': 101, 'SHUNT': 98, 'INTUBATION': 93, 'DISCONNECT': 91, 'CATECHOL': 85, 'KINKEDTUBE': 75, 'MINVOLSET': 72, 'VENTMACH': 62, 'PRESS': 60, 'SAO2': 56, 'VENTLUNG': 50, 'VENTTUBE': 50, 'HR': 49, 'HREKG': 42, 'HRSAT': 42, 'HRBP': 41, 'TPR': 40, 'PVSAT': 37, 'MINVOL': 31, 'CO': 25, 'FIO2': 12, 'BP': 7, 'VENTALV': 5, 'STROKEVOLUME': 1})
So, for example KINKEDTUBE occurs in 75 parent sets. Now, suppose we knew that ARTCO2 and KINKEDTUNE had to be independent. Then we should provide this information to Gobnilp. We do this as follows:
m.add_obligatory_independence(['ARTCO2'],['KINKEDTUBE'])
We can recompute the local scores with this new information available:
m.learn(start='data',end='local scores')
And now we can check that there should be 75 fewer parents for ARTCO2, and also some reduction in the parent sets for KINKEDTUBE.
print(len(m.local_scores['ARTCO2']))
print(len(m.local_scores['KINKEDTUBE']))
398 32
As expected there are $398 = 473 - 75$ fewere parent sets for ARTCO2, and also some reduction in the number of parent sets for KINKEDTUBE. This is because we asserted an independence constraint, which we can get to via the obligatroy_conditional_indepedences
attribute:
m.obligatory_conditional_independences
{(frozenset({'ARTCO2'}), frozenset({'KINKEDTUBE'}), frozenset())}
Note that the independence constraint is stored as a conditional independence constraint with an empty conditioning set (the 3rd set in the tuple above). Also note that Gobnilp
is taking advantage of this constraint before a MIP problem has been created - to help reduce the number of local scores computed. Naturally this constraint is also taken into account when learning a BN (which is done by constructing and solving a MIP model) to ensure that any learned BN satisfies the constraint.
When an independence (or conditional independence) constraint is added (by the user) Gobnilp
deduces that certain other simpler constraints follow. In our current example, it is deduced that
m.forbidden_arrows
{('ARTCO2', 'KINKEDTUBE'), ('KINKEDTUBE', 'ARTCO2')}
m.forbidden_adjacencies
{frozenset({'ARTCO2', 'KINKEDTUBE'})}
Next, let's learn an optimal BN and check that it satisfies the given independence constraint. Since this is a reasonably hard BN learning problem, we'll get the output from Gurobi printed out so we can at least see the progress. We do this by setting the appropriate Gurobi attribute.
m.learn(start='local scores',gurobi_output=True) # construct a MIP model and get Gurobi to solve it
Changed value of parameter OutputFlag to 1 Prev: 0 Min: 0 Max: 1 Default: 1 Changed value of parameter PoolSolutions to 1 Prev: 10 Min: 1 Max: 2000000000 Default: 10 Changed value of parameter LazyConstraints to 1 Prev: 0 Min: 0 Max: 1 Default: 0 Discarded solution information Gurobi Optimizer version 9.0.0 build v9.0.0rc2 (linux64) Optimize a model with 2288 rows, 7349 columns and 158765 nonzeros Model fingerprint: 0x05469d81 Variable types: 0 continuous, 7349 integer (7349 binary) Coefficient statistics: Matrix range [1e+00, 1e+00] Objective range [1e+02, 1e+04] Bounds range [1e+00, 1e+00] RHS range [1e+00, 1e+00] Presolve removed 663 rows and 383 columns Presolve time: 0.59s Presolved: 1625 rows, 6966 columns, 85959 nonzeros Variable types: 0 continuous, 6966 integer (6966 binary) Found heuristic solution: objective -205904.1188 Presolve removed 382 rows and 0 columns Presolved: 1245 rows, 6966 columns, 29441 nonzeros Root relaxation: objective -8.362617e+04, 205 iterations, 0.07 seconds Nodes | Current Node | Objective Bounds | Work Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time 0 0 -83626.171 0 45 -205904.12 -83626.171 59.4% - 0s 0 0 -85269.604 0 80 -205904.12 -85269.604 58.6% - 1s 0 0 -85480.549 0 27 -205904.12 -85480.549 58.5% - 1s 0 0 -86385.401 0 150 -205904.12 -86385.401 58.0% - 1s 0 0 -86613.668 0 38 -205904.12 -86613.668 57.9% - 2s 0 0 -87014.774 0 55 -205904.12 -87014.774 57.7% - 2s 0 0 -87143.296 0 144 -205904.12 -87143.296 57.7% - 2s 0 0 -88898.576 0 213 -205904.12 -88898.576 56.8% - 2s 0 0 -89005.873 0 217 -205904.12 -89005.873 56.8% - 2s 0 0 -89824.417 0 88 -205904.12 -89824.417 56.4% - 2s 0 0 -89896.780 0 246 -205904.12 -89896.780 56.3% - 3s 0 0 -91142.411 0 193 -205904.12 -91142.411 55.7% - 3s 0 0 -92995.992 0 186 -205904.12 -92995.992 54.8% - 3s 0 0 -93483.431 0 180 -205904.12 -93483.431 54.6% - 3s 0 0 -93641.013 0 213 -205904.12 -93641.013 54.5% - 3s 0 0 -94145.102 0 170 -205904.12 -94145.102 54.3% - 3s 0 0 -94604.223 0 202 -205904.12 -94604.223 54.1% - 3s 0 0 -94694.544 0 222 -205904.12 -94694.544 54.0% - 3s 0 0 -97376.832 0 310 -205904.12 -97376.832 52.7% - 4s 0 0 -98992.648 0 256 -205904.12 -98992.648 51.9% - 4s 0 0 -99201.329 0 325 -205904.12 -99201.329 51.8% - 4s 0 0 -101251.98 0 319 -205904.12 -101251.98 50.8% - 4s 0 0 -101572.82 0 328 -205904.12 -101572.82 50.7% - 4s 0 0 -101934.25 0 319 -205904.12 -101934.25 50.5% - 4s 0 0 -102294.43 0 337 -205904.12 -102294.43 50.3% - 5s 0 0 -102508.31 0 337 -205904.12 -102508.31 50.2% - 5s 0 0 -102800.77 0 329 -205904.12 -102800.77 50.1% - 5s 0 0 -103069.80 0 346 -205904.12 -103069.80 49.9% - 5s 0 0 -103214.26 0 341 -205904.12 -103214.26 49.9% - 6s 0 0 -103325.21 0 321 -205904.12 -103325.21 49.8% - 6s 0 0 -103379.05 0 341 -205904.12 -103379.05 49.8% - 6s 0 0 -103587.29 0 376 -205904.12 -103587.29 49.7% - 6s 0 0 -103737.57 0 362 -205904.12 -103737.57 49.6% - 7s 0 0 -103816.43 0 401 -205904.12 -103816.43 49.6% - 7s 0 0 -103978.24 0 356 -205904.12 -103978.24 49.5% - 7s 0 0 -104047.70 0 331 -205904.12 -104047.70 49.5% - 8s 0 0 -104077.72 0 350 -205904.12 -104077.72 49.5% - 8s 0 0 -104085.90 0 275 -205904.12 -104085.90 49.4% - 8s 0 0 -104114.21 0 263 -205904.12 -104114.21 49.4% - 8s 0 0 -104114.48 0 76 -205904.12 -104114.48 49.4% - 8s 0 0 -104193.19 0 324 -205904.12 -104193.19 49.4% - 9s 0 0 -104193.19 0 324 -205904.12 -104193.19 49.4% - 10s 0 0 -104243.53 0 336 -205904.12 -104243.53 49.4% - 14s 0 0 -104243.53 0 336 -205904.12 -104243.53 49.4% - 15s 0 0 -104257.79 0 368 -205904.12 -104257.79 49.4% - 15s 0 0 -104277.94 0 335 -205904.12 -104277.94 49.4% - 15s 0 0 -104284.42 0 331 -205904.12 -104284.42 49.4% - 16s 0 0 -104286.60 0 219 -205904.12 -104286.60 49.4% - 16s 0 0 -104288.40 0 328 -205904.12 -104288.40 49.4% - 16s 0 0 -104404.44 0 211 -205904.12 -104404.44 49.3% - 17s 0 0 -104419.82 0 223 -205904.12 -104419.82 49.3% - 18s 0 0 -104458.37 0 231 -205904.12 -104458.37 49.3% - 19s 0 0 -104489.15 0 376 -205904.12 -104489.15 49.3% - 20s 0 0 -104498.32 0 252 -205904.12 -104498.32 49.2% - 20s 0 0 -104545.94 0 362 -205904.12 -104545.94 49.2% - 21s 0 0 -104550.03 0 380 -205904.12 -104550.03 49.2% - 21s 0 0 -104554.28 0 376 -205904.12 -104554.28 49.2% - 21s 0 0 -104555.92 0 404 -205904.12 -104555.92 49.2% - 22s 0 0 -104589.75 0 361 -205904.12 -104589.75 49.2% - 23s 0 0 -104591.08 0 340 -205904.12 -104591.08 49.2% - 23s 0 0 -104597.12 0 331 -205904.12 -104597.12 49.2% - 24s 0 0 -104599.14 0 383 -205904.12 -104599.14 49.2% - 24s 0 0 -104600.47 0 374 -205904.12 -104600.47 49.2% - 25s 0 0 -104602.19 0 428 -205904.12 -104602.19 49.2% - 25s 0 0 -104603.66 0 386 -205904.12 -104603.66 49.2% - 26s 0 0 -104604.35 0 425 -205904.12 -104604.35 49.2% - 27s 0 0 -104604.35 0 422 -205904.12 -104604.35 49.2% - 29s 0 0 -104604.35 0 422 -205904.12 -104604.35 49.2% - 30s 0 2 -104619.66 0 422 -205904.12 -104619.66 49.2% - 34s 3 6 -104702.79 2 287 -205904.12 -104619.66 49.2% 281 35s 68 71 -105210.03 18 36 -205904.12 -104636.05 49.2% 213 41s 156 149 -105793.19 59 - -205904.12 -104701.70 49.2% 165 46s 245 235 -105077.65 21 146 -205904.12 -104701.96 49.2% 142 50s 320 314 -105305.76 38 137 -205904.12 -104701.96 49.2% 145 55s 465 418 -105067.93 10 303 -205904.12 -104703.32 49.1% 131 60s 557 496 -105318.16 62 79 -205904.12 -104730.35 49.1% 123 65s 680 604 -105519.33 86 24 -205904.12 -104755.42 49.1% 116 70s 804 685 -105256.96 35 422 -205904.12 -104775.45 49.1% 106 78s 806 686 -105415.26 77 133 -205904.12 -104775.45 49.1% 106 80s 822 697 -105357.94 7 388 -205904.12 -104775.45 49.1% 104 85s H 822 661 -105380.9245 -104775.45 0.57% 104 86s 829 666 -105380.92 34 420 -105380.92 -104775.45 0.57% 103 90s H 829 633 -105335.8191 -104775.45 0.53% 103 90s 834 636 -105335.82 61 402 -105335.82 -104775.45 0.53% 102 100s 842 641 -105335.82 53 440 -105335.82 -104795.70 0.51% 102 107s 845 646 -104844.79 11 363 -105335.82 -104806.09 0.50% 109 111s 855 653 -105181.25 14 219 -105335.82 -104892.57 0.42% 112 115s 907 670 -105106.96 26 188 -105335.82 -104892.57 0.42% 118 120s 1000 672 -105028.09 16 233 -105335.82 -104991.77 0.33% 119 125s 1139 700 -105222.53 25 72 -105335.82 -105011.81 0.31% 116 130s 1252 723 -105246.52 35 274 -105335.82 -105028.84 0.29% 114 135s 1383 715 -105295.50 28 250 -105335.82 -105053.05 0.27% 116 140s 1517 715 -105277.77 23 4 -105335.82 -105073.71 0.25% 114 146s 1587 729 -105260.87 22 240 -105335.82 -105102.16 0.22% 113 150s H 1679 668 -105331.8157 -105112.55 0.21% 111 154s 1681 661 -105328.57 57 5 -105331.82 -105118.18 0.20% 112 155s 1772 665 -105240.74 45 108 -105331.82 -105137.80 0.18% 113 161s 1855 649 cutoff 21 -105331.82 -105148.83 0.17% 113 165s 2038 574 -105298.95 27 187 -105331.82 -105168.83 0.15% 110 171s 2182 489 -105315.14 31 - -105331.82 -105193.02 0.13% 110 177s 2267 469 -105315.29 32 - -105331.82 -105207.31 0.12% 110 181s 2380 421 -105317.13 31 - -105331.82 -105216.09 0.11% 108 186s 2608 242 -105331.06 50 6 -105331.82 -105250.11 0.08% 104 193s 2722 245 -105326.39 51 - -105331.82 -105263.63 0.06% 102 195s 2832 233 cutoff 51 -105331.82 -105263.63 0.06% 99.5 200s 3069 249 cutoff 51 -105331.82 -105263.63 0.06% 93.1 208s 3221 316 -105309.24 62 14 -105331.82 -105263.63 0.06% 89.7 212s 3391 381 cutoff 44 -105331.82 -105265.65 0.06% 86.0 217s 3531 436 cutoff 66 -105331.82 -105269.75 0.06% 83.5 222s 3653 430 -105309.61 34 - -105331.82 -105269.75 0.06% 81.5 225s 3906 452 -105305.75 42 10 -105331.82 -105281.87 0.05% 77.3 234s 4256 534 -105304.72 43 19 -105331.82 -105291.81 0.04% 72.3 246s 4615 552 cutoff 44 -105331.82 -105294.08 0.04% 67.9 260s 5048 555 infeasible 67 -105331.82 -105298.01 0.03% 63.8 273s 5513 602 -105331.59 66 - -105331.82 -105298.41 0.03% 60.1 285s 5698 528 cutoff 66 -105331.82 -105304.00 0.03% 58.7 293s 6226 462 cutoff 45 -105331.82 -105308.49 0.02% 55.4 302s 6461 499 cutoff 46 -105331.82 -105309.80 0.02% 54.0 305s 6711 350 -105330.38 48 - -105331.82 -105310.64 0.02% 52.6 311s 6943 401 -105328.10 44 - -105331.82 -105310.64 0.02% 51.3 315s 7154 229 cutoff 48 -105331.82 -105310.64 0.02% 50.2 321s 7343 272 cutoff 44 -105331.82 -105316.75 0.01% 49.2 325s 7538 96 -105321.73 26 5 -105331.82 -105316.75 0.01% 48.4 330s 7906 0 -105329.92 50 11 -105331.82 -105317.78 0.01% 46.7 335s Cutting planes: User: 792 Gomory: 3 Cover: 1 MIR: 5 Flow cover: 14 Zero half: 109 Lazy constraints: 1110 Explored 8178 nodes (381475 simplex iterations) in 338.88 seconds Thread count was 4 (of 4 available processors) Solution count 1: -105332 No other solutions better than -105332 Optimal solution found (tolerance 0.00e+00) Best objective -1.053318157194e+05, best bound -1.053318157194e+05, gap 0.0000% User-callback calls 50750, time in user-callback 149.28 sec ********** BN has score -105331.81571940292 ********** ANAPHYLAXIS<-TPR -471.48556937259855 TPR<- -10901.914716929197 ARTCO2<-VENTALV -1816.018107669981 VENTALV<-INTUBATION -7506.456092395529 BP<-CO,TPR -4927.80460674889 CO<-HR,STROKEVOLUME -2678.9117199075117 CATECHOL<-ARTCO2,TPR -1521.301873422577 HR<-CATECHOL -3744.599767304855 STROKEVOLUME<-HYPOVOLEMIA,LVFAILURE -4437.288938131547 CVP<-LVEDVOLUME -3006.6153021322098 LVEDVOLUME<-HYPOVOLEMIA,LVFAILURE -3613.5057310384145 DISCONNECT<-VENTTUBE -1567.6065855929337 VENTTUBE<-INTUBATION,KINKEDTUBE,VENTLUNG -4542.004684191925 ERRCAUTER<- -3222.576822231058 ERRLOWOUTPUT<- -1972.2788189759885 EXPCO2<-ARTCO2,VENTLUNG -1772.659127224848 VENTLUNG<-INTUBATION,VENTALV -1274.187317113814 FIO2<- -1888.6286705336533 HISTORY<-LVFAILURE -637.4903136947833 LVFAILURE<- -2010.543210255084 HRBP<-ERRLOWOUTPUT,HR -1331.6218853034952 HREKG<-ERRCAUTER,HR -1507.9736736634077 HRSAT<-ERRCAUTER,HR -1452.085102920697 HYPOVOLEMIA<- -4990.780506105162 INSUFFANESTH<- -3303.7325344726705 INTUBATION<- -3428.3602719444025 KINKEDTUBE<- -1558.2831144287484 MINVOL<-INTUBATION,VENTLUNG -2061.9057414226554 MINVOLSET<-VENTMACH -946.556491622774 VENTMACH<-DISCONNECT,VENTTUBE -1896.468405037318 PAP<-PULMEMBOLUS -4008.254809145714 PULMEMBOLUS<- -583.1474353126687 PCWP<-LVEDVOLUME -2233.6875379606645 PRESS<-INTUBATION,KINKEDTUBE,VENTTUBE -8478.131859909474 PVSAT<-FIO2,VENTALV -846.3768118197622 SAO2<-PVSAT,SHUNT -1137.5909747223777 SHUNT<-INTUBATION,PULMEMBOLUS -2052.9806218132726 ********** bnlearn modelstring = [ANAPHYLAXIS|TPR][TPR][ARTCO2|VENTALV][VENTALV|INTUBATION][BP|TPR:CO][CO|HR:STROKEVOLUME][CATECHOL|TPR:ARTCO2][HR|CATECHOL][STROKEVOLUME|LVFAILURE:HYPOVOLEMIA][CVP|LVEDVOLUME][LVEDVOLUME|LVFAILURE:HYPOVOLEMIA][DISCONNECT|VENTTUBE][VENTTUBE|VENTLUNG:INTUBATION:KINKEDTUBE][ERRCAUTER][ERRLOWOUTPUT][EXPCO2|VENTLUNG:ARTCO2][VENTLUNG|INTUBATION:VENTALV][FIO2][HISTORY|LVFAILURE][LVFAILURE][HRBP|HR:ERRLOWOUTPUT][HREKG|HR:ERRCAUTER][HRSAT|HR:ERRCAUTER][HYPOVOLEMIA][INSUFFANESTH][INTUBATION][KINKEDTUBE][MINVOL|VENTLUNG:INTUBATION][MINVOLSET|VENTMACH][VENTMACH|VENTTUBE:DISCONNECT][PAP|PULMEMBOLUS][PULMEMBOLUS][PCWP|LVEDVOLUME][PRESS|INTUBATION:KINKEDTUBE:VENTTUBE][PVSAT|FIO2:VENTALV][SAO2|PVSAT:SHUNT][SHUNT|PULMEMBOLUS:INTUBATION] ********** CPDAG: Vertices: ANAPHYLAXIS,TPR,ARTCO2,VENTALV,BP,CO,CATECHOL,HR,STROKEVOLUME,CVP,LVEDVOLUME,DISCONNECT,VENTTUBE,ERRCAUTER,ERRLOWOUTPUT,EXPCO2,VENTLUNG,FIO2,HISTORY,LVFAILURE,HRBP,HREKG,HRSAT,HYPOVOLEMIA,INSUFFANESTH,INTUBATION,KINKEDTUBE,MINVOL,MINVOLSET,VENTMACH,PAP,PULMEMBOLUS,PCWP,PRESS,PVSAT,SAO2,SHUNT TPR-ANAPHYLAXIS TPR->BP TPR->CATECHOL ARTCO2->CATECHOL ARTCO2->EXPCO2 VENTALV-ARTCO2 VENTALV->PVSAT VENTALV-VENTLUNG CO->BP CATECHOL->HR HR->CO HR->HRBP HR->HREKG HR->HRSAT STROKEVOLUME->CO LVEDVOLUME->CVP LVEDVOLUME->PCWP DISCONNECT-VENTMACH VENTTUBE->DISCONNECT VENTTUBE->PRESS VENTTUBE->VENTMACH ERRCAUTER->HREKG ERRCAUTER->HRSAT ERRLOWOUTPUT->HRBP VENTLUNG->EXPCO2 VENTLUNG-MINVOL VENTLUNG->VENTTUBE FIO2->PVSAT LVFAILURE-HISTORY LVFAILURE->LVEDVOLUME LVFAILURE->STROKEVOLUME HYPOVOLEMIA->LVEDVOLUME HYPOVOLEMIA->STROKEVOLUME INTUBATION-MINVOL INTUBATION->PRESS INTUBATION->SHUNT INTUBATION-VENTALV INTUBATION-VENTLUNG INTUBATION->VENTTUBE KINKEDTUBE->PRESS KINKEDTUBE->VENTTUBE VENTMACH->MINVOLSET PULMEMBOLUS-PAP PULMEMBOLUS->SHUNT PVSAT->SAO2 SHUNT->SAO2
Gobnilp/Gurobi took quite a while to find a provably optimal BN in this case (and recall that we are running with the default parent set limit of 3 - were we to raise that solving would take longer still). As is normal Gobnilp spent some time adding cutting planes to reduce the upper bound on the score of an optimal BN from an original upper bound of -83626.171 down to -104619.66. Then Gurobi builds a search tree by branching on the values of MIP variables, finally finding a BN with a reasonable score at node 822 of this tree. Although at this point we know that the incumbent is close to the best possible BN (gap = 0.57%) it takes quite a while longer and thousands more search nodes before we can be sure we have an optimal BN.
It is interesting to compare this to what happens when we learn from this dataset with no extra constraints:
m2 = Gobnilp()
m2.learn('alarm_10000.dat',gurobi_output=True)
Changed value of parameter PreCrush to 1 Prev: 0 Min: 0 Max: 1 Default: 0 Changed value of parameter CutPasses to 100000 Prev: -1 Min: -1 Max: 2000000000 Default: -1 Changed value of parameter GomoryPasses to 100000 Prev: -1 Min: -1 Max: 2000000000 Default: -1 Changed value of parameter MIPFocus to 2 Prev: 0 Min: 0 Max: 3 Default: 0 Changed value of parameter ZeroHalfCuts to 2 Prev: -1 Min: -1 Max: 2 Default: -1 Changed value of parameter MIPGap to 0.0 Prev: 0.0001 Min: 0.0 Max: inf Default: 0.0001 Changed value of parameter MIPGapAbs to 0.0 Prev: 1e-10 Min: 0.0 Max: inf Default: 1e-10 Parameter OutputFlag unchanged Value: 1 Min: 0 Max: 1 Default: 1 Changed value of parameter PoolSolutions to 1 Prev: 10 Min: 1 Max: 2000000000 Default: 10 Changed value of parameter LazyConstraints to 1 Prev: 0 Min: 0 Max: 1 Default: 0 Discarded solution information Gurobi Optimizer version 9.0.0 build v9.0.0rc2 (linux64) Optimize a model with 2304 rows, 7437 columns and 161370 nonzeros Model fingerprint: 0xd35b1531 Variable types: 0 continuous, 7437 integer (7437 binary) Coefficient statistics: Matrix range [1e+00, 1e+00] Objective range [1e+02, 1e+04] Bounds range [1e+00, 1e+00] RHS range [1e+00, 1e+00] Presolve removed 663 rows and 382 columns Presolve time: 0.61s Presolved: 1641 rows, 7055 columns, 87103 nonzeros Variable types: 0 continuous, 7055 integer (7055 binary) Found heuristic solution: objective -205904.1188 Presolve removed 384 rows and 0 columns Presolved: 1259 rows, 7055 columns, 30194 nonzeros Root relaxation: objective -8.362617e+04, 209 iterations, 0.08 seconds Nodes | Current Node | Objective Bounds | Work Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time 0 0 -83626.171 0 45 -205904.12 -83626.171 59.4% - 0s 0 0 -85269.604 0 80 -205904.12 -85269.604 58.6% - 1s 0 0 -86855.165 0 119 -205904.12 -86855.165 57.8% - 1s 0 0 -89330.686 0 124 -205904.12 -89330.686 56.6% - 2s 0 0 -90135.061 0 28 -205904.12 -90135.061 56.2% - 2s 0 0 -90607.947 0 35 -205904.12 -90607.947 56.0% - 2s 0 0 -90928.085 0 95 -205904.12 -90928.085 55.8% - 2s 0 0 -91351.872 0 97 -205904.12 -91351.872 55.6% - 3s 0 0 -93564.967 0 187 -205904.12 -93564.967 54.6% - 4s 0 0 -93855.817 0 168 -205904.12 -93855.817 54.4% - 4s 0 0 -96373.054 0 280 -205904.12 -96373.054 53.2% - 4s 0 0 -96925.946 0 272 -205904.12 -96925.946 52.9% - 4s 0 0 -98423.047 0 293 -205904.12 -98423.047 52.2% - 4s 0 0 -98926.822 0 265 -205904.12 -98926.822 52.0% - 5s 0 0 -99592.098 0 308 -205904.12 -99592.098 51.6% - 5s 0 0 -99999.530 0 251 -205904.12 -99999.530 51.4% - 5s 0 0 -101214.96 0 293 -205904.12 -101214.96 50.8% - 5s 0 0 -102461.94 0 303 -205904.12 -102461.94 50.2% - 5s 0 0 -102943.08 0 279 -205904.12 -102943.08 50.0% - 6s 0 0 -103233.78 0 275 -205904.12 -103233.78 49.9% - 6s 0 0 -103346.74 0 238 -205904.12 -103346.74 49.8% - 6s 0 0 -103455.02 0 298 -205904.12 -103455.02 49.8% - 6s 0 0 -103474.22 0 300 -205904.12 -103474.22 49.7% - 6s 0 0 -103678.28 0 300 -205904.12 -103678.28 49.6% - 7s 0 0 -103814.01 0 294 -205904.12 -103814.01 49.6% - 7s 0 0 -103880.24 0 268 -205904.12 -103880.24 49.5% - 7s 0 0 -104035.50 0 277 -205904.12 -104035.50 49.5% - 8s 0 0 -104097.82 0 356 -205904.12 -104097.82 49.4% - 9s 0 0 -104122.78 0 307 -205904.12 -104122.78 49.4% - 9s 0 0 -104315.35 0 272 -205904.12 -104315.35 49.3% - 10s H 0 0 -107316.9764 -104315.35 2.80% - 10s 0 0 -104361.56 0 103 -107316.98 -104361.56 2.75% - 11s 0 0 -104371.37 0 316 -107316.98 -104371.37 2.74% - 11s 0 0 -104482.11 0 271 -107316.98 -104482.11 2.64% - 11s 0 0 -104492.67 0 330 -107316.98 -104492.67 2.63% - 13s 0 0 -104493.91 0 288 -107316.98 -104493.91 2.63% - 13s 0 0 -104525.72 0 354 -107316.98 -104525.72 2.60% - 14s 0 0 -104533.65 0 310 -107316.98 -104533.65 2.59% - 14s 0 0 -104559.47 0 281 -107316.98 -104559.47 2.57% - 14s 0 0 -104559.76 0 279 -107316.98 -104559.76 2.57% - 15s 0 0 -104578.83 0 181 -107316.98 -104578.83 2.55% - 16s 0 0 -104579.25 0 163 -107316.98 -104579.25 2.55% - 16s 0 0 -104590.93 0 362 -107316.98 -104590.93 2.54% - 17s 0 0 -104591.84 0 335 -107316.98 -104591.84 2.54% - 17s 0 0 -104595.89 0 305 -107316.98 -104595.89 2.54% - 18s H 0 0 -106678.6779 -104595.89 1.95% - 19s 0 0 -104597.72 0 331 -106678.68 -104597.72 1.95% - 19s 0 0 -104602.74 0 328 -106678.68 -104602.74 1.95% - 20s 0 0 -104603.30 0 294 -106678.68 -104603.30 1.95% - 20s 0 0 -104614.14 0 322 -106678.68 -104614.14 1.94% - 21s 0 0 -104614.14 0 322 -106678.68 -104614.14 1.94% - 26s 0 0 -104616.57 0 354 -106678.68 -104616.57 1.93% - 29s 0 0 -104620.54 0 387 -106678.68 -104620.54 1.93% - 30s 0 0 -104621.47 0 405 -106678.68 -104621.47 1.93% - 30s 0 0 -104625.45 0 352 -106678.68 -104625.45 1.92% - 31s 0 0 -104626.82 0 382 -106678.68 -104626.82 1.92% - 31s 0 0 -104627.84 0 395 -106678.68 -104627.84 1.92% - 32s 0 0 -104628.36 0 399 -106678.68 -104628.36 1.92% - 32s 0 0 -104641.34 0 361 -106678.68 -104641.34 1.91% - 33s 0 0 -104642.64 0 349 -106678.68 -104642.64 1.91% - 34s 0 0 -104653.53 0 321 -106678.68 -104653.53 1.90% - 34s 0 0 -104655.18 0 364 -106678.68 -104655.18 1.90% - 35s 0 0 -104662.42 0 356 -106678.68 -104662.42 1.89% - 36s 0 0 -104664.81 0 359 -106678.68 -104664.81 1.89% - 36s 0 0 -104666.09 0 358 -106678.68 -104666.09 1.89% - 36s 0 0 -104669.01 0 323 -106678.68 -104669.01 1.88% - 37s 0 0 -104669.27 0 347 -106678.68 -104669.27 1.88% - 37s 0 0 -104678.47 0 378 -106678.68 -104678.47 1.87% - 38s 0 0 -104679.25 0 368 -106678.68 -104679.25 1.87% - 39s 0 0 -104680.62 0 383 -106678.68 -104680.62 1.87% - 40s 0 0 -104680.66 0 383 -106678.68 -104680.66 1.87% - 40s 0 0 -104682.86 0 404 -106678.68 -104682.86 1.87% - 41s 0 0 -104682.99 0 404 -106678.68 -104682.99 1.87% - 41s 0 0 -104684.53 0 411 -106678.68 -104684.53 1.87% - 42s 0 0 -104684.65 0 423 -106678.68 -104684.65 1.87% - 42s 0 0 -104686.65 0 392 -106678.68 -104686.65 1.87% - 43s 0 0 -104687.03 0 389 -106678.68 -104687.03 1.87% - 43s 0 0 -104701.25 0 379 -106678.68 -104701.25 1.85% - 44s 0 0 -104703.02 0 426 -106678.68 -104703.02 1.85% - 45s 0 0 -104703.15 0 415 -106678.68 -104703.15 1.85% - 46s H 0 0 -105226.5119 -104703.15 0.50% - 47s 0 0 -104722.17 0 85 -105226.51 -104722.17 0.48% - 49s 0 0 -104722.17 0 85 -105226.51 -104722.17 0.48% - 50s 0 0 -104722.17 0 236 -105226.51 -104722.17 0.48% - 51s 0 0 -104722.17 0 133 -105226.51 -104722.17 0.48% - 51s 0 0 -104722.17 0 321 -105226.51 -104722.17 0.48% - 51s 0 0 -104722.17 0 404 -105226.51 -104722.17 0.48% - 51s 0 0 -104722.17 0 389 -105226.51 -104722.17 0.48% - 52s 0 0 -104722.17 0 425 -105226.51 -104722.17 0.48% - 52s 0 0 -104722.17 0 388 -105226.51 -104722.17 0.48% - 52s 0 0 -104722.17 0 385 -105226.51 -104722.17 0.48% - 53s 0 0 -104722.17 0 419 -105226.51 -104722.17 0.48% - 54s 0 0 -104722.17 0 317 -105226.51 -104722.17 0.48% - 55s 0 0 -104722.17 0 355 -105226.51 -104722.17 0.48% - 55s 0 0 -104722.17 0 414 -105226.51 -104722.17 0.48% - 56s 0 0 -104722.17 0 412 -105226.51 -104722.17 0.48% - 57s 0 2 -104722.17 0 412 -105226.51 -104722.17 0.48% - 60s 124 83 -105141.98 40 155 -105226.51 -104756.15 0.45% 155 65s 229 111 cutoff 32 -105226.51 -104826.92 0.38% 148 70s 353 170 -104990.43 14 353 -105226.51 -104845.79 0.36% 141 75s 424 200 cutoff 26 -105226.51 -104878.16 0.33% 146 80s 558 228 cutoff 19 -105226.51 -104947.94 0.26% 140 85s 742 242 -105193.27 25 209 -105226.51 -105000.91 0.21% 136 90s 881 232 -105097.21 11 169 -105226.51 -105031.30 0.19% 133 95s 1128 246 cutoff 24 -105226.51 -105061.14 0.16% 124 100s 1356 203 cutoff 24 -105226.51 -105099.49 0.12% 119 105s 1571 78 -105167.11 33 53 -105226.51 -105151.51 0.07% 113 110s Cutting planes: User: 403 Gomory: 5 MIR: 1 Flow cover: 4 Zero half: 73 Lazy constraints: 6 Explored 1945 nodes (204343 simplex iterations) in 112.85 seconds Thread count was 4 (of 4 available processors) Solution count 1: -105227 No other solutions better than -105227 Optimal solution found (tolerance 0.00e+00) Best objective -1.052265119168e+05, best bound -1.052265119168e+05, gap 0.0000% User-callback calls 21479, time in user-callback 24.94 sec ********** BN has score -105226.51191678 ********** ANAPHYLAXIS<-TPR -471.48556937259855 TPR<- -10901.914716929197 ARTCO2<-VENTALV -1816.018107669981 VENTALV<-INTUBATION,VENTLUNG -1785.1568955858238 BP<-CO,TPR -4927.80460674889 CO<-HR,STROKEVOLUME -2678.9117199075117 CATECHOL<-ARTCO2,TPR -1521.301873422577 HR<-CATECHOL -3744.599767304855 STROKEVOLUME<-HYPOVOLEMIA,LVFAILURE -4437.288938131547 CVP<-LVEDVOLUME -3006.6153021322098 LVEDVOLUME<-HYPOVOLEMIA,LVFAILURE -3613.5057310384145 DISCONNECT<- -3295.0317430806754 ERRCAUTER<- -3222.576822231058 ERRLOWOUTPUT<- -1972.2788189759885 EXPCO2<-ARTCO2,VENTLUNG -1772.659127224848 VENTLUNG<-INTUBATION,KINKEDTUBE,VENTTUBE -3618.0085247885872 FIO2<- -1888.6286705336533 HISTORY<-LVFAILURE -637.4903136947833 LVFAILURE<- -2010.543210255084 HRBP<-ERRLOWOUTPUT,HR -1331.6218853034952 HREKG<-ERRCAUTER,HR -1507.9736736634077 HRSAT<-ERRCAUTER,HR -1452.085102920697 HYPOVOLEMIA<- -4990.780506105162 INSUFFANESTH<- -3303.7325344726705 INTUBATION<- -3428.3602719444025 KINKEDTUBE<- -1558.2831144287484 MINVOL<-INTUBATION,VENTLUNG -2061.9057414226554 MINVOLSET<- -3952.777366733397 PAP<- -4194.791057672541 PCWP<-LVEDVOLUME -2233.6875379606645 PRESS<-INTUBATION,KINKEDTUBE,VENTTUBE -8478.131859909474 VENTTUBE<-DISCONNECT,VENTMACH -1814.1386675404356 PULMEMBOLUS<-PAP -396.6111867858417 PVSAT<-FIO2,VENTALV -846.3768118197622 SAO2<-PVSAT,SHUNT -1137.5909747223777 SHUNT<-INTUBATION,PULMEMBOLUS -2052.9806218132726 VENTMACH<-MINVOLSET -3162.862542532719 ********** bnlearn modelstring = [ANAPHYLAXIS|TPR][TPR][ARTCO2|VENTALV][VENTALV|VENTLUNG:INTUBATION][BP|TPR:CO][CO|HR:STROKEVOLUME][CATECHOL|TPR:ARTCO2][HR|CATECHOL][STROKEVOLUME|LVFAILURE:HYPOVOLEMIA][CVP|LVEDVOLUME][LVEDVOLUME|LVFAILURE:HYPOVOLEMIA][DISCONNECT][ERRCAUTER][ERRLOWOUTPUT][EXPCO2|VENTLUNG:ARTCO2][VENTLUNG|INTUBATION:KINKEDTUBE:VENTTUBE][FIO2][HISTORY|LVFAILURE][LVFAILURE][HRBP|HR:ERRLOWOUTPUT][HREKG|HR:ERRCAUTER][HRSAT|HR:ERRCAUTER][HYPOVOLEMIA][INSUFFANESTH][INTUBATION][KINKEDTUBE][MINVOL|VENTLUNG:INTUBATION][MINVOLSET][PAP][PCWP|LVEDVOLUME][PRESS|INTUBATION:KINKEDTUBE:VENTTUBE][VENTTUBE|VENTMACH:DISCONNECT][PULMEMBOLUS|PAP][PVSAT|FIO2:VENTALV][SAO2|PVSAT:SHUNT][SHUNT|PULMEMBOLUS:INTUBATION][VENTMACH|MINVOLSET] ********** CPDAG: Vertices: ANAPHYLAXIS,TPR,ARTCO2,VENTALV,BP,CO,CATECHOL,HR,STROKEVOLUME,CVP,LVEDVOLUME,DISCONNECT,ERRCAUTER,ERRLOWOUTPUT,EXPCO2,VENTLUNG,FIO2,HISTORY,LVFAILURE,HRBP,HREKG,HRSAT,HYPOVOLEMIA,INSUFFANESTH,INTUBATION,KINKEDTUBE,MINVOL,MINVOLSET,PAP,PCWP,PRESS,VENTTUBE,PULMEMBOLUS,PVSAT,SAO2,SHUNT,VENTMACH TPR-ANAPHYLAXIS TPR->BP TPR->CATECHOL ARTCO2->CATECHOL ARTCO2->EXPCO2 VENTALV->ARTCO2 VENTALV->PVSAT CO->BP CATECHOL->HR HR->CO HR->HRBP HR->HREKG HR->HRSAT STROKEVOLUME->CO LVEDVOLUME->CVP LVEDVOLUME->PCWP DISCONNECT->VENTTUBE ERRCAUTER->HREKG ERRCAUTER->HRSAT ERRLOWOUTPUT->HRBP VENTLUNG->EXPCO2 VENTLUNG->MINVOL VENTLUNG->VENTALV FIO2->PVSAT LVFAILURE-HISTORY LVFAILURE->LVEDVOLUME LVFAILURE->STROKEVOLUME HYPOVOLEMIA->LVEDVOLUME HYPOVOLEMIA->STROKEVOLUME INTUBATION-MINVOL INTUBATION->PRESS INTUBATION->SHUNT INTUBATION-VENTALV INTUBATION->VENTLUNG KINKEDTUBE->PRESS KINKEDTUBE->VENTLUNG MINVOLSET-VENTMACH PAP-PULMEMBOLUS VENTTUBE->PRESS VENTTUBE->VENTLUNG PULMEMBOLUS->SHUNT PVSAT->SAO2 SHUNT->SAO2 VENTMACH->VENTTUBE
Without the independence constraint solving is a lot faster, even though the set of allowed BNs is greater without the constraint. This is basically because Gobnilp is optimised for the case where the only constraint is that the BN is acyclic.
The number of "Lazy constraints" posted in the two runs is worth examining: without the constraint there are only 6 but with the independence constraint there are 1100. This is because Gobnilp is adding constraints ruling out BNs (strictly speaking ancestral subgraphs of BNs) which do not satisfy the independence constraint. This is a crude method of enforcing (conditional) independence constraints and this adds to the slowness of solving.