text:
Back to the main Index

# Third (basic) lesson with Abinit and AbiPy

Crystalline silicon.

This lesson aims at showing you how to get the following physical properties, for an insulator:

• the total energy
• the lattice parameter
• the Kohn-Sham band structure

You will learn about the use of k-points, as well as the smearing of the plane-wave kinetic energy cut-off.

This tutorial is a complement to the standard ABINIT tutorial on silicon. Here, powerful flow and visualisation procedures will be demonstrated. Still, some basic understanding of the stand-alone working of ABINIT is a prerequisite. Also, in order to fully benefit from this Abipy tutorial, other more basic Abipy tutorials should have been followed, as suggested in the abitutorials index page.

In [1]:
# Use this at the beginning of your script so that your code will be compatible with python3
from __future__ import print_function, division, unicode_literals

import numpy as np

import warnings
warnings.filterwarnings("ignore")  # Ignore warnings

from abipy import abilab
abilab.enable_notebook() # This line tells AbiPy we are running inside a notebook

# This line configures matplotlib to show figures embedded in the notebook.
# Replace inline with notebook in classic notebook
%matplotlib inline

# Option available in jupyterlab. See https://github.com/matplotlib/jupyter-matplotlib
#%matplotlib widget

## Computing the total energy of silicon at fixed number of k-points¶

Our goal is to study the convergence of the total energy of silicon versus the number of k-points. So we start by defining a function that generates a Flow of SCF calculations by looping over a predefined list of ngkpt values. The crystalline structure is initialized from a CIF file while other parameters such as the cutoff energy are fixed:

In [2]:
from lesson_base3 import build_ngkpt_flow
abilab.print_source(build_ngkpt_flow)

Out[2]:
def build_ngkpt_flow(options):
"""
Crystalline silicon: computation of the total energy
Convergence with respect to the number of k points. Similar to tbase3_3.in

Args:
options: Command line options.

Return:
Abinit Flow object.
"""
# Definition of the different grids
ngkpt_list = [(2, 2, 2), (4, 4, 4), (6, 6, 6), (8, 8, 8)]

# These shifts will be the same for all grids
shiftk = [float(s) for s in "0.5 0.5 0.5 0.5 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.5".split()]

# Build MultiDataset object (container of ndtset inputs).
# Structure is initialized from CIF file.
multi = abilab.MultiDataset(structure=abidata.cif_file("si.cif"),
pseudos=abidata.pseudos("14si.pspnc"), ndtset=len(ngkpt_list))

# These variables are the same in each input.
multi.set_vars(ecut=8, toldfe=1e-6, diemac=12.0, iomode=3)

# Each input has its own value of ngkpt. shiftk is constant.
for i, ngkpt in enumerate(ngkpt_list):
multi[i].set_kmesh(ngkpt=ngkpt, shiftk=shiftk)

workdir = options.workdir if (options and options.workdir) else "flow_base3_ngkpt"

# Split the inputs by calling multi.datasets() and pass the list of inputs to Flow.from_inputs.
return flowtk.Flow.from_inputs(workdir, inputs=multi.split_datasets())

Let's call the function to build the flow:

In [3]:
flow = build_ngkpt_flow(options=None)

In total, we have four ScfTasks that will be executed in the flow_base3_ngkpt directory:

In [4]:
flow.get_graphviz()

Out[4]:

This is the input of the first task w0_t0:

In [5]:
flow[0][0].input

Out[5]:
############################################################################################
# SECTION: basic
############################################################################################
ecut 8
toldfe 1e-06
ngkpt 2 2 2
kptopt 1
nshiftk 4
shiftk
0.5 0.5 0.5
0.5 0.0 0.0
0.0 0.5 0.0
0.0 0.0 0.5
############################################################################################
# SECTION: dev
############################################################################################
iomode 3
############################################################################################
# SECTION: gstate
############################################################################################
diemac 12.0
############################################################################################
# STRUCTURE
############################################################################################
natom 2
ntypat 1
typat 1 1
znucl 14
xred
0.0000000000 0.0000000000 0.0000000000
0.2500000000 0.2500000000 0.2500000000
acell 1.0 1.0 1.0
rprim
6.3285005272 0.0000000000 3.6537614829
2.1095001757 5.9665675167 3.6537614829
0.0000000000 0.0000000000 7.3075229659

and these are the divisions of the k-mesh for the four different calculations:

In [6]:
for task in flow.iflat_tasks():

w0_t0 uses ngkpt: (2, 2, 2)
w0_t1 uses ngkpt: (4, 4, 4)
w0_t2 uses ngkpt: (6, 6, 6)
w0_t3 uses ngkpt: (8, 8, 8)

but we can achieve the same goal with:

In [7]:
flow.get_vars_dataframe("ngkpt", "ecut")

Out[7]:
ngkpt ecut class
w0_t0 (2, 2, 2) 8 ScfTask
w0_t1 (4, 4, 4) 8 ScfTask
w0_t2 (6, 6, 6) 8 ScfTask
w0_t3 (8, 8, 8) 8 ScfTask

At this point, we could run the flow in the notebook by just calling:

flow.make_scheduler().start()

or, alternatively, execute the lesson_base3.py script to build the directory with the flow and then use:

abirun.py flow_base3_ngkpt scheduler

inside the terminal.

## Analysis of the results¶

We could use the API provided by the flow to extract the total energies from the GSR files. Something like:

nkpt_list, ene_list = [], []
nkpt_list.append(gsr.nkpt)
ene_list.append(gsr.energy)

# Assuming values are already sorted wrt nkpt
import matplotlib.pyplot as plt
plt.plot(nkpt_list, ene_list, marker="o");

but it is much easier to create a GsrRobot that will do the work for us:

In [8]:
robot_enekpt = abilab.GsrRobot.from_dir("flow_base3_ngkpt")
robot_enekpt

Out[8]:
1. w0/t0/outdata/out_GSR.nc
2. w0/t1/outdata/out_GSR.nc
3. w0/t2/outdata/out_GSR.nc
4. w0/t3/outdata/out_GSR.nc

In the next lines, we are going to generate a pandas Dataframe with the most important results so that we can show how to use the pandas API to analyze the data:

In [9]:
ene_table = robot_enekpt.get_dataframe()
ene_table.keys()

Out[9]:
Index(['formula', 'natom', 'alpha', 'beta', 'gamma', 'a', 'b', 'c', 'volume',
'abispg_num', 'spglib_symb', 'spglib_num', 'spglib_lattice_type',
'energy', 'pressure', 'max_force', 'ecut', 'pawecutdg', 'tsmear',
'nkpt', 'nsppol', 'nspinor', 'nspden'],
dtype='object')

The dataframe contains several columns but we are mainly interested in the number of k-points nkpt and in the energy (given in eV). Let's massage a bit the data to facilitate the post-processing:

In [10]:
# We are gonna plot f(nkpt) so let's sort the rows first.
ene_table.sort_values(by="nkpt", inplace=True)

# Add a column with energies in Ha and another column with the difference wrt to the last point.
ene_table["energy_Ha"] = ene_table["energy"] * abilab.units.eV_to_Ha
ene_table["ediff_Ha"] = ene_table["energy_Ha"] - ene_table["energy_Ha"][-1]

before printing a subset of the columns with the syntax:

In [11]:
ene_table[["nkpt", "energy", "energy_Ha", "ediff_Ha"]]

Out[11]:
nkpt energy energy_Ha ediff_Ha
w0/t0/outdata/out_GSR.nc 2 -241.251546 -8.865831 0.006242
w0/t1/outdata/out_GSR.nc 10 -241.417959 -8.871946 0.000126
w0/t2/outdata/out_GSR.nc 28 -241.421158 -8.872064 0.000009
w0/t3/outdata/out_GSR.nc 60 -241.421391 -8.872073 0.000000

If you do not like tables and prefer figures, use:

In [12]:
ene_table.plot(x="nkpt", y=["energy_Ha", "ediff_Ha", "pressure"], style="-o", subplots=True);

The difference between dataset 3 and dataset 4 is rather small. Even the dataset 2 gives an accuracy of about 0.0001 Ha. So, our converged value for the total energy (at fixed acell and ecut) is -8.8726 Ha.

Now that we have learned a bit how to use pandas Dataframes, we can finally reveal that the AbiPy robots already provide methods to perform this kind of convergence studies so that we do not need to manipulate pandas dataframes explicitly. For example, we can perform the same analysis with a single line:

In [13]:
robot_enekpt.plot_gsr_convergence(sortby="nkpt");

We can also pass a function that will be called by the robot to compute the values along the x-axis and sort the results. The docstring of the function is used as label of the x-axis:

In [14]:
def inv_nkpt(abifile):
r"""$\dfrac{1}{nkpt}$"""
return 1 / abifile.nkpt

robot_enekpt.plot_gsr_convergence(sortby=inv_nkpt);

In [15]:
#robot_enekpt.plot_lattice_convergence(sortby="nkpt");

## Determination of the lattice parameters¶

At this point, the original Abinit tutorial proceeds with a convergence study for the optimized lattice parameters as function of the k-point sampling. In AbiPy, we only need to build different relaxation tasks with a slightly different input in which only ngkpt is changed.

In [16]:
from lesson_base3 import build_relax_flow
abilab.print_source(build_relax_flow)

Out[16]:
def build_relax_flow(options):
"""
Crystalline silicon: computation of the optimal lattice parameter.
Convergence with respect to the number of k points. Similar to tbase3_4.in
"""
# Structural relaxation for different k-point samplings.
ngkpt_list = [(2, 2, 2), (4, 4, 4)]

shiftk = [float(s) for s in "0.5 0.5 0.5 0.5 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.5".split()]

multi = abilab.MultiDataset(structure=abidata.cif_file("si.cif"),
pseudos=abidata.pseudos("14si.pspnc"), ndtset=len(ngkpt_list))

# Global variables
multi.set_vars(
ecut=8,
tolvrs=1e-9,
optcell=1,
ionmov=3,
ntime=10,
dilatmx=1.05,
ecutsm=0.5,
diemac=12,
iomode=3,
)

for i, ngkpt in enumerate(ngkpt_list):
multi[i].set_kmesh(ngkpt=ngkpt, shiftk=shiftk)

workdir = options.workdir if (options and options.workdir) else "flow_base3_relax"

In [17]:
relax_flow = build_relax_flow(options=None)
relax_flow.get_graphviz()

Out[17]:
If you want to run the flow from the shell, open lesson_base3.py and change the main function so that it calls build_relax_flow instead of build_ngkpt_flow.

This is our first structural relaxation with AbiPy and this gives us the opportunity to introduce the HIST.nc file. This file stores the history of the relaxation (energies, forces, stresses, lattice parameters and atomic positions at the different relaxation steps).

As usual, we use abiopen to generate an AbiPy object:

In [18]:
hist = abilab.abiopen("flow_base3_relax/w0/t1/outdata/out_HIST.nc")
print(hist)

================================= File Info =================================
Name: out_HIST.nc
Directory: /Users/gmatteo/git_repos/abitutorials/abitutorials/base3/flow_base3_relax/w0/t1/outdata
Size: 4.24 kb
Access Time: Mon Aug 13 16:49:10 2018
Modification Time: Wed Oct 11 00:28:35 2017
Change Time: Wed Oct 11 00:28:35 2017

============================= Initial Structure =============================
Full Formula (Si2)
Reduced Formula: Si
abc   :   3.866975   3.866975   3.866975
angles:  60.000000  60.000000  60.000000
Sites (2)
#  SP       a     b     c  cartesian_forces
---  ----  ----  ----  ----  -----------------------------------------------------------
0  Si    0     0     0     [ 4.10231606e-27 -5.80155101e-27 -3.31586260e-27] eV ang^-1
1  Si    0.25  0.25  0.25  [-4.10231606e-27  5.80155101e-27  3.31586260e-27] eV ang^-1

Number of relaxation steps performed: 4
============================== Final structure ==============================
Full Formula (Si2)
Reduced Formula: Si
abc   :   3.822962   3.822962   3.822962
angles:  60.000000  60.000000  60.000000
Sites (2)
#  SP       a      b     c  cartesian_forces
---  ----  ----  -----  ----  -----------------------------------------------------------
0  Si    0     -0     0     [5.53272686e-28 6.25956589e-27 9.58296402e-28] eV ang^-1
1  Si    0.25   0.25  0.25  [-5.53272686e-28 -6.25956589e-27 -9.58296402e-28] eV ang^-1

Volume change in percentage: -3.38%
Percentage lattice parameter changes:
a: -1.14%, b: -1.14%, c: -1%

Stress tensor (Cartesian coordinates in GPa):
[[8.75674654e-03 5.30998779e-11 7.50894518e-11]
[5.30998779e-11 8.75674654e-03 5.30995448e-11]
[7.50894518e-11 5.30995448e-11 8.75674652e-03]]
Pressure: -0.009 [GPa]

To plot the evolution of the most important physical quantities, use:

In [19]:
hist.plot();

In [20]:
hist_robot = abilab.HistRobot.from_dir("flow_base3_relax")

hist_table = hist_robot.get_dataframe()

There are several entries in the DataFrame:

In [21]:
hist_table.keys()

Out[21]:
Index(['formula', 'natom', 'alpha', 'beta', 'gamma', 'a', 'b', 'c', 'volume',
'abispg_num', 'spglib_symb', 'spglib_num', 'spglib_lattice_type',
'num_steps', 'final_energy', 'final_pressure', 'final_fmin',
'final_fmax', 'final_fmean', 'final_fstd', 'final_drift',
'initial_fmin', 'initial_fmax', 'initial_fmean', 'initial_fstd',
'initial_drift'],
dtype='object')

Let's select some of them with:

In [22]:
hist_table[["alpha", "a", "final_energy", "final_pressure", "num_steps"]]

Out[22]:
alpha a final_energy final_pressure num_steps
w0/t0/outdata/out_HIST.nc 60.0 3.829282 -241.255628 -0.013259 4
w0/t1/outdata/out_HIST.nc 60.0 3.822962 -241.425215 -0.008757 4

and print the evolution of important physical properties extracted from the two files:

In [23]:
hist_robot.gridplot(what_list=["energy", "abc", "pressure", "forces"]);

We can also compare the two structural relaxations with:

In [24]:
hist_robot.combiplot();

Unfortunately the HIST.nc file does not have enough metadata. In particular we would like to have information about the k-point sampling so that we can analyze the convergence of the optimized lattice parameters wrt nkpt. Fortunately the GSR.nc has all the information we need and it is just a matter of replacing the HistRobot with a GsrRobot:

In [25]:
with abilab.GsrRobot.from_dir("flow_base3_relax") as relkpt_robot:
relax_table = relkpt_robot.get_dataframe().sort_values(by="nkpt")
dfs = relkpt_robot.get_structure_dataframes()

In [26]:
relax_table[["energy", "a", "pressure", "max_force", "pressure"]]

Out[26]:
energy a pressure max_force pressure
w0/t0/outdata/out_GSR.nc -241.255628 3.829282 -0.013259 1.913429e-27 -0.013259
w0/t1/outdata/out_GSR.nc -241.425215 3.822962 -0.008757 6.356619e-27 -0.008757

Plotting the energy, the lattice parameter a in Bohr and the pressure in GPa vs nkpt is really a piece of cake!

In [27]:
relax_table.plot(x="nkpt", y=["energy", "a", "pressure"], subplots=True);

Alternatively, one can use the GsrRobot API:

In [28]:
relkpt_robot.plot_gsr_convergence(sortby="nkpt");

In [29]:
relkpt_robot.plot_lattice_convergence(what_list=["a"], sortby="nkpt");

We fix the parameters acell to the theoretical value of 3*10.216, and we fix also the grid of k points (the 4x4x4 FCC grid, equivalent to a 8x8x8 Monkhorst-pack grid). We will ask for 8 bands (4 valence and 4 conduction).

## Computing the band structure¶

A band structure can be computed by solving the Kohn-Sham equation for many different k points, along different lines of the Brillouin zone. The potential that enters the Kohn-Sham equation must be derived from a previous self-consistent calculation, and will not vary during the scan of different k-point lines.

This is our first Flow with dependencies in the sense that the band structure calculation must be connected to a previous SCF run. Fortunately AbiPy provides a factory function to generate this kind of workflow. We only need to focus on the definition of the two inputs:

In [30]:
from lesson_base3 import build_ebands_flow
abilab.print_source(build_ebands_flow)

Out[30]:
def build_ebands_flow(options):
"""
Band structure calculation.
First, a SCF density computation, then a non-SCF band structure calculation.
Similar to tbase3_5.in
"""
multi = abilab.MultiDataset(structure=abidata.cif_file("si.cif"),
pseudos=abidata.pseudos("14si.pspnc"), ndtset=2)
# Global variables
shiftk = [float(s) for s in "0.5 0.5 0.5 0.5 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.5".split()]
multi.set_vars(ecut=8, diemac=12, iomode=3)

# Dataset 1
multi[0].set_vars(tolvrs=1e-9)
multi[0].set_kmesh(ngkpt=[4, 4, 4], shiftk=shiftk)

# Dataset 2
multi[1].set_vars(tolwfr=1e-15)
multi[1].set_kpath(ndivsm=5)

scf_input, nscf_input = multi.split_datasets()

workdir = options.workdir if (options and options.workdir) else "flow_base3_ebands"

return flowtk.bandstructure_flow(workdir, scf_input=scf_input, nscf_input=nscf_input)

The Flow consists of a single Work with two Tasks (ScfTask with a k-mesh and a NscfTask performed on the k-path).

In [31]:
ebands_flow = build_ebands_flow(options=None)
ebands_flow.get_graphviz()

Out[31]:
If you want to run the flow from the shell, open lesson_base3.py and change the main function so that it calls build_ebands_flow.

Let's extract the band structure from the GSR.nc file produced by the NscfTask:

In [32]:
with abilab.abiopen("flow_base3_ebands/w0/t1/outdata/out_GSR.nc") as gsr:
ebands_kpath = gsr.ebands

and plot it with:

In [33]:
ebands_kpath.plot();

Visual inspection reveals that the width of the valence band is ~11.8 eV, the lowest unoccupied state at X is ~0.5 eV higher than the top of the valence band at $\Gamma$. Bulk silicon is described as an indirect band gap material (this is correct), with a band-gap of about 0.5 eV (this is quantitatively quite wrong: the experimental value is 1.17 eV at 25 degree Celsius, the famous DFT band-gap problem). The minimum of the conduction band is even slightly displaced with respect to X.

Unfortunately, it seems that AbiPy does not agree with us:

In [34]:
print(ebands_kpath)

================================= Structure =================================
Full Formula (Si2)
Reduced Formula: Si
abc   :   3.866975   3.866975   3.866975
angles:  60.000000  60.000000  60.000000
Sites (2)
#  SP       a     b     c
---  ----  ----  ----  ----
0  Si    0     0     0
1  Si    0.25  0.25  0.25

Abinit Spacegroup: spgid: 227, num_spatial_symmetries: 48, has_timerev: True, symmorphic: True

Number of electrons: 8.0, Fermi level: 5.280 (eV)
nsppol: 1, nkpt: 101, mband: 5, nspinor: 1, nspden: 1
smearing scheme: none, tsmear_eV: 0.272, occopt: 1
Direct gap:
Energy: 0.309 (eV)
Initial state: spin=0, kpt=[+0.050, +0.050, +0.100], weight: 0.000, band=2, eig=5.218, occ=2.000
Final state:   spin=0, kpt=[+0.050, +0.050, +0.100], weight: 0.000, band=3, eig=5.527, occ=2.000
Fundamental gap:
Energy: 0.021 (eV)
Initial state: spin=0, kpt=[+0.125, +0.125, +0.125], weight: 0.000, band=3, eig=5.280, occ=2.000
Final state:   spin=0, kpt=[+0.042, +0.042, +0.042], weight: 0.000, band=1, eig=5.301, occ=2.000
Bandwidth: 11.564 (eV)
Valence maximum located at:
spin=0, kpt=[+0.125, +0.125, +0.125], weight: 0.000, band=3, eig=5.280, occ=2.000
Conduction minimum located at:
spin=0, kpt=[+0.042, +0.042, +0.042], weight: 0.000, band=1, eig=5.301, occ=2.000

The reason is that ther Fermi energy in ebands_kpath is not completely consistent with the band structure. The Fermi energy, indeed, has been taken from the previous GS-SCF calculation performed on a shifted k-mesh, the $\Gamma$ point was not included and therefore the Fermi energy is underestimated.

To fix this problem we have to change manually the Fermi energy and set it to the maximum of the valence bands:

In [35]:
ebands_kpath.set_fermie_to_vbm()
print(ebands_kpath)

================================= Structure =================================
Full Formula (Si2)
Reduced Formula: Si
abc   :   3.866975   3.866975   3.866975
angles:  60.000000  60.000000  60.000000
Sites (2)
#  SP       a     b     c
---  ----  ----  ----  ----
0  Si    0     0     0
1  Si    0.25  0.25  0.25

Abinit Spacegroup: spgid: 227, num_spatial_symmetries: 48, has_timerev: True, symmorphic: True

Number of electrons: 8.0, Fermi level: 5.567 (eV)
nsppol: 1, nkpt: 101, mband: 5, nspinor: 1, nspden: 1
smearing scheme: none, tsmear_eV: 0.272, occopt: 1
Direct gap:
Energy: 2.512 (eV)
Initial state: spin=0, kpt=[+0.000, +0.000, +0.000], name: $\Gamma$, weight: 0.000, band=3, eig=5.567, occ=2.000
Final state:   spin=0, kpt=[+0.000, +0.000, +0.000], name: $\Gamma$, weight: 0.000, band=4, eig=8.080, occ=0.000
Fundamental gap:
Energy: 0.517 (eV)
Initial state: spin=0, kpt=[+0.000, +0.000, +0.000], name: $\Gamma$, weight: 0.000, band=3, eig=5.567, occ=2.000
Final state:   spin=0, kpt=[+0.429, +0.000, +0.429], weight: 0.000, band=4, eig=6.084, occ=0.000
Bandwidth: 11.852 (eV)
Valence maximum located at:
spin=0, kpt=[+0.000, +0.000, +0.000], name: $\Gamma$, weight: 0.000, band=3, eig=5.567, occ=2.000
Conduction minimum located at:
spin=0, kpt=[+0.429, +0.000, +0.429], weight: 0.000, band=4, eig=6.084, occ=0.000

Now the AbiPy results are consistent with our initial analysis:

In [36]:
ebands_kpath.plot(with_gaps=True);

In [37]:
# We can also plot the k-path in the Brillouin zone with:
#ebands_kpath.kpoints.plot();

The GSR file produced by the first task contains energies on a homogeneous k-mesh. We can therefore compute the DOS by invoking the get_edos method:

In [38]:
with abilab.abiopen("flow_base3_ebands/w0/t0/outdata/out_GSR.nc") as gsr:
ebands_kmesh = gsr.ebands

edos = ebands_kmesh.get_edos()

and plot the DOS with:

In [39]:
edos.plot();

where the zero of the energy axis is set to the Fermi level $\epsilon_F$ obtained by solving:

$$\int_{-\infty}^{\epsilon_F} g(\epsilon)\,d\epsilon = N$$

for $\epsilon_F$ with $N$ the number of electrons per unit cell. Note that the DOS is highly sensitive to the sampling of the IBZ and to the value of the broadening, especially in metallic systems.

In [40]:
edos.plot_dos_idos();

Want to make a nice picture of the band dispersion with a second panel for the DOS?

In [41]:
ebands_kpath.plot_with_edos(edos);

It is important to stress that each panel in the above figure is aligned with respect to its own Fermi energy and these values are not necessarly equal:

In [42]:
print(ebands_kpath.fermie, edos.fermie)

5.567493548211215 eV 5.838293265593386

We can always plot the bands and the DOS without setting their Fermi energy to zero by using:

In [43]:
ebands_kpath.plot_with_edos(edos, e0=0);

This figure shows that the bands and the DOS are not perfectly aligned. More specifically we would expect the DOS to be zero at the bottom/top of the conduction. This problems is essentially due to the use of a relatively large gaussian broadening. One should therefore compute the DOS with a much denser IBZ mesh and a much smaller broadening to solve this alignment issue

Back to the main Index

In [ ]:


text:
# Bayesian Networks¶

Bayesian networks are a probabilistic model that are especially good at inference given incomplete data. Much like a hidden Markov model, they consist of a directed graphical model (though Bayesian networks must also be acyclic) and a set of probability distributions. The edges encode dependency statements between the variables, where the lack of an edge between any pair of variables indicates a conditional independence. Each node encodes a probability distribution, where root nodes encode univariate probability distributions and inner/leaf nodes encode conditional probability distributions. Bayesian networks are exceptionally flexible when doing inference, as any subset of variables can be observed, and inference done over all other variables, without needing to define these groups in advance. In fact, the set of observed variables can change from one sample to the next without needing to modify the underlying algorithm at all.

Currently, pomegranate only supports discrete Bayesian networks, meaning that the values must be categories, i.e. ‘apples’ and ‘oranges’, or 1 and 2, where 1 and 2 refer to categories, not numbers, and so 2 is not explicitly ‘bigger’ than 1.

## Initialization¶

Bayesian networks can be initialized in two ways, depending on whether the underlying graphical structure is known or not: (1) the graphical structure can be built one node at a time with pre-initialized distributions set for each node, or (2) both the graphical structure and distributions can be learned directly from data. This mirrors the other models that are implemented in pomegranate. However, typically expectation maximization is used to fit the parameters of the distribution, and so initialization (such as through k-means) is typically fast whereas fitting is slow. For Bayesian networks, the opposite is the case. Fitting can be done quickly by just summing counts through the data, while initialization is hard as it requires an exponential time search through all possible DAGs to identify the optimal graph. More is discussed in the tutorials above and in the fitting section below.

Let’s take a look at initializing a Bayesian network in the first manner by quickly implementing the Monty Hall problem. The Monty Hall problem arose from the gameshow Let’s Make a Deal, where a guest had to choose which one of three doors had a prize behind it. The twist was that after the guest chose, the host, originally Monty Hall, would then open one of the doors the guest did not pick and ask if the guest wanted to switch which door they had picked. Initial inspection may lead you to believe that if there are only two doors left, there is a 50-50 chance of you picking the right one, and so there is no advantage one way or the other. However, it has been proven both through simulations and analytically that there is in fact a 66% chance of getting the prize if the guest switches their door, regardless of the door they initially went with.

Our network will have three nodes, one for the guest, one for the prize, and one for the door Monty chooses to open. The door the guest initially chooses and the door the prize is behind are uniform random processes across the three doors, but the door which Monty opens is dependent on both the door the guest chooses (it cannot be the door the guest chooses), and the door the prize is behind (it cannot be the door with the prize behind it).

from pomegranate import *

guest = DiscreteDistribution({'A': 1./3, 'B': 1./3, 'C': 1./3})
prize = DiscreteDistribution({'A': 1./3, 'B': 1./3, 'C': 1./3})
monty = ConditionalProbabilityTable(
[['A', 'A', 'A', 0.0],
['A', 'A', 'B', 0.5],
['A', 'A', 'C', 0.5],
['A', 'B', 'A', 0.0],
['A', 'B', 'B', 0.0],
['A', 'B', 'C', 1.0],
['A', 'C', 'A', 0.0],
['A', 'C', 'B', 1.0],
['A', 'C', 'C', 0.0],
['B', 'A', 'A', 0.0],
['B', 'A', 'B', 0.0],
['B', 'A', 'C', 1.0],
['B', 'B', 'A', 0.5],
['B', 'B', 'B', 0.0],
['B', 'B', 'C', 0.5],
['B', 'C', 'A', 1.0],
['B', 'C', 'B', 0.0],
['B', 'C', 'C', 0.0],
['C', 'A', 'A', 0.0],
['C', 'A', 'B', 1.0],
['C', 'A', 'C', 0.0],
['C', 'B', 'A', 1.0],
['C', 'B', 'B', 0.0],
['C', 'B', 'C', 0.0],
['C', 'C', 'A', 0.5],
['C', 'C', 'B', 0.5],
['C', 'C', 'C', 0.0]], [guest, prize])

s1 = Node(guest, name="guest")
s2 = Node(prize, name="prize")
s3 = Node(monty, name="monty")

model = BayesianNetwork("Monty Hall Problem")
model.bake()

Note

The objects ‘state’ and ‘node’ are really the same thing and can be used interchangable. The only difference is the name, as hidden Markov models use ‘state’ in the literature frequently whereas Bayesian networks use ‘node’ frequently.

The conditional distribution must be explicitly spelled out in this example, followed by a list of the parents in the same order as the columns take in the tabble that is provided (e.g. the columns in the table correspond to guest, prize, monty, probability.)

However, one can also initialize a Bayesian network based completely on data. As mentioned before, the exact version of this algorithm takes exponential time with the number of variables and typically can’t be done on more than ~25 variables. This is because there are a super-exponential number of directed acyclic graphs that one could define over a set of variables, but fortunately one can use dynamic programming in order to reduce this complexity down to “simply exponential.” The implementation of the exact algorithm actually goes further than the original dynamic programing algorithm by implementing an A* search to somewhat reduce computational time but drastically reduce required memory, sometimes by an order of magnitude.

from pomegranate import *
import numpy

model = BayesianNetwork.from_samples(X, algorithm='exact')

The exact algorithm is not the default, though. The default is a novel greedy algorithm that greedily chooses a topological ordering of the variables, but optimally identifies the best parents for each variable given this ordering. It is significantly faster and more memory efficient than the exact algorithm and produces far better estimates than using a Chow-Liu tree. This is set to the default to avoid locking up the computers of users that unintentionally tell their computers to do a near-impossible task.

## Probability¶

You can calculate the probabiity of a sample under a Bayesian network as the product of the probability of each variable given its parents, if it has any. This can be expressed as $$P = \prod\limits_{i=1}^{d} P(D_{i}|Pa_{i})$$ for a sample with $d$ dimensions. For example, in the Monty Hal problem, the probability of a show is the probability of the guest choosing the respective door, times the probability of the prize being behind a given door, times the probability of Monty opening a given door given the previous two values. For example, using the manually initialized network above:

>>> print model.probability([['A', 'A', 'A'],
['A', 'A', 'B'],
['C', 'C', 'B']])
[ 0.          0.05555556  0.05555556]

## Prediction¶

Bayesian networks are frequently used to infer/impute the value of missing variables given the observed values. In other models, typically there is either a single or fixed set of missing variables, such as latent factors, that need to be imputed, and so returning a fixed vector or matrix as the predictions makes sense. However, in the case of Bayesian networks, we can make no such assumptions, and so when data is passed in for prediction it should be in the format as a matrix with None in the missing variables that need to be inferred. The return is thus a filled in matrix where the Nones have been replaced with the imputed values. For example:

>>> print model.predict([['A', 'B', None],
['A', 'C', None],
['C', 'B', None]])
[['A' 'B' 'C']
['A' 'C' 'B']
['C' 'B' 'A']]

In this example, the final column is the one that is always missing, but a more complex example is as follows:

>>> print model.predict([['A', 'B', None],
['A', None, 'C'],
[None, 'B', 'A']])
[['A' 'B' 'C']
['A' 'B' 'C']
['C' 'B' 'A']]

## Fitting¶

Fitting a Bayesian network to data is a fairly simple process. Essentially, for each variable, you need consider only that column of data and the columns corresponding to that variables parents. If it is a univariate distribution, then the maximum likelihood estimate is just the count of each symbol divided by the number of samples in the data. If it is a multivariate distribution, it ends up being the probability of each symbol in the variable of interest given the combination of symbols in the parents. For example, consider a binary dataset with two variables, X and Y, where X is a parent of Y. First, we would go through the dataset and calculate P(X=0) and P(X=1). Then, we would calculate P(Y=0|X=0), P(Y=1|X=0), P(Y=0|X=1), and P(Y=1|X=1). Those values encode all of the parameters of the Bayesian network.

## API Reference¶

class pomegranate.BayesianNetwork.BayesianNetwork

A Bayesian Network Model.

A Bayesian network is a directed graph where nodes represent variables, edges represent conditional dependencies of the children on their parents, and the lack of an edge represents a conditional independence.

Parameters: name : str, optional The name of the model. Default is None

Attributes

 states (list, shape (n_states,)) A list of all the state objects in the model graph (networkx.DiGraph) The underlying graph object.
add_edge()

Add a transition from state a to state b which indicates that B is dependent on A in ways specified by the distribution.

add_node()

Add a node to the graph.

add_nodes()

Add multiple states to the graph.

add_state()

Another name for a node.

add_states()

Another name for a node.

add_transition()

Transitions and edges are the same.

bake()

Finalize the topology of the model.

Assign a numerical index to every state and create the underlying arrays corresponding to the states and edges between the states. This method must be called before any of the probability-calculating methods. This includes converting conditional probability tables into joint probability tables and creating a list of both marginal and table nodes.

Parameters: None None
clear_summaries()

Clear the summary statistics stored in the object.

copy()

Return a deep copy of this distribution object.

This object will not be tied to any other distribution or connected in any form.

Parameters: None distribution : Distribution A copy of the distribution with the same parameters.
dense_transition_matrix()

Returns the dense transition matrix. Useful if the transitions of somewhat small models need to be analyzed.

edge_count()

Returns the number of edges present in the model.

fit()

Fit the model to data using MLE estimates.

Fit the model to the data by updating each of the components of the model, which are univariate or multivariate distributions. This uses a simple MLE estimate to update the distributions according to their summarize or fit methods.

This is a wrapper for the summarize and from_summaries methods.

Parameters: items : array-like, shape (n_samples, n_nodes) The data to train on, where each row is a sample and each column corresponds to the associated variable. weights : array-like, shape (n_nodes), optional The weight of each sample as a positive double. Default is None. inertia : double, optional The inertia for updating the distributions, passed along to the distribution method. Default is 0.0. pseudocount : double, optional A pseudocount to add to the emission of each distribution. This effectively smoothes the states to prevent 0. probability symbols if they don’t happen to occur in the data. Only effects hidden Markov models defined over discrete distributions. Default is 0. None
freeze()

Freeze the distribution, preventing updates from occuring.

from_json()

Read in a serialized Bayesian Network and return the appropriate object.

Parameters: s : str A JSON formatted string containing the file. model : object A properly initialized and baked model.
from_samples()

Learn the structure of the network from data.

Find the structure of the network from data using a Bayesian structure learning score. This currently enumerates all the exponential number of structures and finds the best according to the score. This allows weights on the different samples as well. The score that is optimized is the minimum description length (MDL).

Parameters: X : array-like, shape (n_samples, n_nodes) The data to fit the structure too, where each row is a sample and each column corresponds to the associated variable. weights : array-like, shape (n_nodes), optional The weight of each sample as a positive double. Default is None. algorithm : str, one of ‘chow-liu’, ‘greedy’, ‘exact’, ‘exact-dp’ optional The algorithm to use for learning the Bayesian network. Default is ‘greedy’ that greedily attempts to find the best structure, and frequently can identify the optimal structure. ‘exact’ uses DP/A* to find the optimal Bayesian network, and ‘exact-dp’ tries to find the shortest path on the entire order lattice, which is more memory and computationally expensive. ‘exact’ and ‘exact-dp’ should give identical results, with ‘exact-dp’ remaining an option mostly for debugging reasons. ‘chow-liu’ will return the optimal tree-like structure for the Bayesian network, which is a very fast approximation but not always the best network. max_parents : int, optional The maximum number of parents a node can have. If used, this means using the k-learn procedure. Can drastically speed up algorithms. If -1, no max on parents. Default is -1. root : int, optional For algorithms which require a single root (‘chow-liu’), this is the root for which all edges point away from. User may specify which column to use as the root. Default is the first column. constraint_graph : networkx.DiGraph or None, optional A directed graph showing valid parent sets for each variable. Each node is a set of variables, and edges represent which variables can be valid parents of those variables. The naive structure learning task is just all variables in a single node with a self edge, meaning that you know nothing about pseudocount : double, optional A pseudocount to add to the emission of each distribution. This effectively smoothes the states to prevent 0. probability symbols if they don’t happen to occur in the data. Default is 0. state_names : array-like, shape (n_nodes), optional A list of meaningful names to be applied to nodes name : str, optional The name of the model. Default is None. n_jobs : int, optional The number of threads to use when learning the structure of the network. If a constraint graph is provided, this will parallelize the tasks as directed by the constraint graph. If one is not provided it will parallelize the building of the parent graphs. Both cases will provide large speed gains. model : BayesianNetwork The learned BayesianNetwork.
from_structure()

Return a Bayesian network from a predefined structure.

Pass in the structure of the network as a tuple of tuples and get a fit network in return. The tuple should contain n tuples, with one for each node in the graph. Each inner tuple should be of the parents for that node. For example, a three node graph where both node 0 and 1 have node 2 as a parent would be specified as ((2,), (2,), ()).

Parameters: X : array-like, shape (n_samples, n_nodes) The data to fit the structure too, where each row is a sample and each column corresponds to the associated variable. structure : tuple of tuples The parents for each node in the graph. If a node has no parents, then do not specify any parents. weights : array-like, shape (n_nodes), optional The weight of each sample as a positive double. Default is None. pseudocount : double, optional A pseudocount to add to the emission of each distribution. This effectively smoothes the states to prevent 0. probability symbols if they don’t happen to occur in the data. Default is 0. name : str, optional The name of the model. Default is None. state_names : array-like, shape (n_nodes), optional A list of meaningful names to be applied to nodes model : BayesianNetwoork A Bayesian network with the specified structure.
from_summaries()

Use MLE on the stored sufficient statistics to train the model.

This uses MLE estimates on the stored sufficient statistics to train the model.

Parameters: inertia : double, optional The inertia for updating the distributions, passed along to the distribution method. Default is 0.0. pseudocount : double, optional A pseudocount to add to the emission of each distribution. This effectively smoothes the states to prevent 0. probability symbols if they don’t happen to occur in the data. Default is 0. None
log_probability()

Return the log probability of samples under the Bayesian network.

The log probability is just the sum of the log probabilities under each of the components. The log probability of a sample under the graph A -> B is just P(A)*P(B|A). This will return a vector of log probabilities, one for each sample.

Parameters: X : array-like, shape (n_samples, n_dim) The sample is a vector of points where each dimension represents the same variable as added to the graph originally. It doesn’t matter what the connections between these variables are, just that they are all ordered the same. logp : double The log probability of that sample.
marginal()

Return the marginal probabilities of each variable in the graph.

This is equivalent to a pass of belief propogation on a graph where no data has been given. This will calculate the probability of each variable being in each possible emission when nothing is known.

Parameters: None marginals : array-like, shape (n_nodes) An array of univariate distribution objects showing the marginal probabilities of that variable.
node_count()

Returns the number of nodes/states in the model

plot()

Draw this model’s graph using NetworkX and matplotlib.

Note that this relies on networkx’s built-in graphing capabilities (and not Graphviz) and thus can’t draw self-loops.

See networkx.draw_networkx() for the keywords you can pass in.

Parameters: **kwargs : any The arguments to pass into networkx.draw_networkx() None
predict()

Predict missing values of a data matrix using MLE.

Impute the missing values of a data matrix using the maximally likely predictions according to the forward-backward algorithm. Run each sample through the algorithm (predict_proba) and replace missing values with the maximally likely predicted emission.

Parameters: items : array-like, shape (n_samples, n_nodes) Data matrix to impute. Missing values must be either None (if lists) or np.nan (if numpy.ndarray). Will fill in these values with the maximally likely ones. max_iterations : int, optional Number of iterations to run loopy belief propogation for. Default is 100. items : numpy.ndarray, shape (n_samples, n_nodes) This is the data matrix with the missing values imputed.
predict_proba()

Returns the probabilities of each variable in the graph given evidence.

This calculates the marginal probability distributions for each state given the evidence provided through loopy belief propogation. Loopy belief propogation is an approximate algorithm which is exact for certain graph structures.

Parameters: data : dict or array-like, shape <= n_nodes, optional The evidence supplied to the graph. This can either be a dictionary with keys being state names and values being the observed values (either the emissions or a distribution over the emissions) or an array with the values being ordered according to the nodes incorporation in the graph (the order fed into .add_states/add_nodes) and None for variables which are unknown. If nothing is fed in then calculate the marginal of the graph. Default is {}. max_iterations : int, optional The number of iterations with which to do loopy belief propogation. Usually requires only 1. Default is 100. check_input : bool, optional Check to make sure that the observed symbol is a valid symbol for that distribution to produce. Default is True. probabilities : array-like, shape (n_nodes) An array of univariate distribution objects showing the probabilities of each variable.
probability()

Return the probability of the given symbol under this distribution.

Parameters: symbol : object The symbol to calculate the probability of probability : double The probability of that point under the distribution.
sample()

Return a random item sampled from this distribution.

Parameters: n : int or None, optional The number of samples to return. Default is None, which is to generate a single sample. sample : double or object Returns a sample from the distribution of a type in the support of the distribution.
state_count()

Returns the number of states present in the model.

summarize()

Summarize a batch of data and store the sufficient statistics.

This will partition the dataset into columns which belong to their appropriate distribution. If the distribution has parents, then multiple columns are sent to the distribution. This relies mostly on the summarize function of the underlying distribution.

Parameters: items : array-like, shape (n_samples, n_nodes) The data to train on, where each row is a sample and each column corresponds to the associated variable. weights : array-like, shape (n_nodes), optional The weight of each sample as a positive double. Default is None. None
thaw()

Thaw the distribution, re-allowing updates to occur.

to_json()

Serialize the model to a JSON.

Parameters: separators : tuple, optional The two separaters to pass to the json.dumps function for formatting. indent : int, optional The indentation to use at each level. Passed to json.dumps for formatting. json : str A properly formatted JSON object.
class pomegranate.BayesianNetwork.ParentGraph

Generate a parent graph for a single variable over its parents.

This will generate the parent graph for a single parents given the data. A parent graph is the dynamically generated best parent set and respective score for each combination of parent variables. For example, if we are generating a parent graph for x1 over x2, x3, and x4, we may calculate that having x2 as a parent is better than x2,x3 and so store the value of x2 in the node for x2,x3.

Parameters: X : numpy.ndarray, shape=(n, d) The data to fit the structure too, where each row is a sample and each column corresponds to the associated variable. weights : numpy.ndarray, shape=(n,) The weight of each sample as a positive double. Default is None. key_count : numpy.ndarray, shape=(d,) The number of unique keys in each column. pseudocount : double A pseudocount to add to each possibility. max_parents : int The maximum number of parents a node can have. If used, this means using the k-learn procedure. Can drastically speed up algorithms. If -1, no max on parents. Default is -1. parent_set : tuple, default () The variables which are possible parents for this variable. If nothing is passed in then it defaults to all other variables, as one would expect in the naive case. This allows for cases where we want to build a parent graph over only a subset of the variables. structure : tuple, shape=(d,) The parents for each variable in this SCC
pomegranate.BayesianNetwork.discrete_exact_a_star()

Find the optimal graph over a set of variables with no other knowledge.

This is the naive dynamic programming structure learning task where the optimal graph is identified from a set of variables using an order graph and parent graphs. This can be used either when no constraint graph is provided or for a SCC which is made up of a node containing a self-loop. It uses DP/A* in order to find the optimal graph without considering all possible topological sorts. A greedy version of the algorithm can be used that massively reduces both the computational and memory cost while frequently producing the optimal graph.

Parameters: X : numpy.ndarray, shape=(n, d) The data to fit the structure too, where each row is a sample and each column corresponds to the associated variable. weights : numpy.ndarray, shape=(n,) The weight of each sample as a positive double. Default is None. key_count : numpy.ndarray, shape=(d,) The number of unique keys in each column. pseudocount : double A pseudocount to add to each possibility. max_parents : int The maximum number of parents a node can have. If used, this means using the k-learn procedure. Can drastically speed up algorithms. If -1, no max on parents. Default is -1. n_jobs : int The number of threads to use when learning the structure of the network. This parallelizes the creation of the parent graphs. structure : tuple, shape=(d,) The parents for each variable in this SCC
pomegranate.BayesianNetwork.discrete_exact_component()

Find the optimal graph over a multi-node component of the constaint graph.

The general algorithm in this case is to begin with each variable and add all possible single children for that entry recursively until completion. This will result in a far sparser order graph than before. In addition, one can eliminate entries from the parent graphs that contain invalid parents as they are a fast of computational time.

Parameters: X : numpy.ndarray, shape=(n, d) The data to fit the structure too, where each row is a sample and each column corresponds to the associated variable. weights : numpy.ndarray, shape=(n,) The weight of each sample as a positive double. Default is None. key_count : numpy.ndarray, shape=(d,) The number of unique keys in each column. pseudocount : double A pseudocount to add to each possibility. max_parents : int The maximum number of parents a node can have. If used, this means using the k-learn procedure. Can drastically speed up algorithms. If -1, no max on parents. Default is -1. n_jobs : int The number of threads to use when learning the structure of the network. This parallelizes the creation of the parent graphs. structure : tuple, shape=(d,) The parents for each variable in this SCC
pomegranate.BayesianNetwork.discrete_exact_dp()

Find the optimal graph over a set of variables with no other knowledge.

This is the naive dynamic programming structure learning task where the optimal graph is identified from a set of variables using an order graph and parent graphs. This can be used either when no constraint graph is provided or for a SCC which is made up of a node containing a self-loop. This is a reference implementation that uses the naive shortest path algorithm over the entire order graph. The ‘exact’ option uses the A* path in order to avoid considering the full order graph.

Parameters: X : numpy.ndarray, shape=(n, d) The data to fit the structure too, where each row is a sample and each column corresponds to the associated variable. weights : numpy.ndarray, shape=(n,) The weight of each sample as a positive double. Default is None. key_count : numpy.ndarray, shape=(d,) The number of unique keys in each column. pseudocount : double A pseudocount to add to each possibility. max_parents : int The maximum number of parents a node can have. If used, this means using the k-learn procedure. Can drastically speed up algorithms. If -1, no max on parents. Default is -1. n_jobs : int The number of threads to use when learning the structure of the network. This parallelizes the creation of the parent graphs. structure : tuple, shape=(d,) The parents for each variable in this SCC
pomegranate.BayesianNetwork.discrete_exact_slap()

Find the optimal graph in a node with a Self Loop And Parents (SLAP).

Instead of just performing exact BNSL over the set of all parents and removing the offending edges there are efficiencies that can be gained by considering the structure. In particular, parents not coming from the main node do not need to be considered in the order graph but simply added to each entry after creation of the order graph. This is because those variables occur earlier in the topological ordering but it doesn’t matter how they occur otherwise. Parent graphs must be defined over all variables however.

Parameters: X : numpy.ndarray, shape=(n, d) The data to fit the structure too, where each row is a sample and each column corresponds to the associated variable. weights : numpy.ndarray, shape=(n,) The weight of each sample as a positive double. Default is None. key_count : numpy.ndarray, shape=(d,) The number of unique keys in each column. pseudocount : double A pseudocount to add to each possibility. max_parents : int The maximum number of parents a node can have. If used, this means using the k-learn procedure. Can drastically speed up algorithms. If -1, no max on parents. Default is -1. n_jobs : int The number of threads to use when learning the structure of the network. This parallelizes the creation of the parent graphs. structure : tuple, shape=(d,) The parents for each variable in this SCC
pomegranate.BayesianNetwork.discrete_exact_with_constraints()

This returns the optimal Bayesian network given a set of constraints.

This function controls the process of learning the Bayesian network by taking in a constraint graph, identifying the strongly connected components (SCCs) and solving each one using the appropriate algorithm. This is mostly an internal function.

Parameters: X : numpy.ndarray, shape=(n, d) The data to fit the structure too, where each row is a sample and each column corresponds to the associated variable. weights : numpy.ndarray, shape=(n,) The weight of each sample as a positive double. Default is None. key_count : numpy.ndarray, shape=(d,) The number of unique keys in each column. pseudocount : double A pseudocount to add to each possibility. max_parents : int The maximum number of parents a node can have. If used, this means using the k-learn procedure. Can drastically speed up algorithms. If -1, no max on parents. Default is -1. constraint_graph : networkx.DiGraph A directed graph showing valid parent sets for each variable. Each node is a set of variables, and edges represent which variables can be valid parents of those variables. The naive structure learning task is just all variables in a single node with a self edge, meaning that you know nothing about n_jobs : int The number of threads to use when learning the structure of the network. This parallelized both the creation of the parent graphs for each variable and the solving of the SCCs. structure : tuple, shape=(d,) The parents for each variable in the network.
pomegranate.BayesianNetwork.discrete_exact_with_constraints_task()

This is a wrapper for the function to be parallelzied by joblib.

This function takes in a single task as an id and a set of parents and children and calls the appropriate function. This is mostly a wrapper for joblib to parallelize.

Parameters: X : numpy.ndarray, shape=(n, d) The data to fit the structure too, where each row is a sample and each column corresponds to the associated variable. weights : numpy.ndarray, shape=(n,) The weight of each sample as a positive double. Default is None. key_count : numpy.ndarray, shape=(d,) The number of unique keys in each column. pseudocount : double A pseudocount to add to each possibility. max_parents : int The maximum number of parents a node can have. If used, this means using the k-learn procedure. Can drastically speed up algorithms. If -1, no max on parents. Default is -1. task : tuple A 3-tuple containing the id, the set of parents and the set of children to learn a component of the Bayesian network over. The cases represent a SCC of the following: 0 - Self loop and no parents 1 - Self loop and parents 2 - Parents and no self loop 3 - Multiple nodes n_jobs : int The number of threads to use when learning the structure of the network. This parallelizes the creation of the parent graphs for each task or the finding of best parents in case 2. structure : tuple, shape=(d,) The parents for each variable in this SCC
pomegranate.BayesianNetwork.discrete_greedy()

Find the optimal graph over a set of variables with no other knowledge.

This is the naive dynamic programming structure learning task where the optimal graph is identified from a set of variables using an order graph and parent graphs. This can be used either when no constraint graph is provided or for a SCC which is made up of a node containing a self-loop. It uses DP/A* in order to find the optimal graph without considering all possible topological sorts. A greedy version of the algorithm can be used that massively reduces both the computational and memory cost while frequently producing the optimal graph.

Parameters: X : numpy.ndarray, shape=(n, d) The data to fit the structure too, where each row is a sample and each column corresponds to the associated variable. weights : numpy.ndarray, shape=(n,) The weight of each sample as a positive double. Default is None. key_count : numpy.ndarray, shape=(d,) The number of unique keys in each column. pseudocount : double A pseudocount to add to each possibility. max_parents : int The maximum number of parents a node can have. If used, this means using the k-learn procedure. Can drastically speed up algorithms. If -1, no max on parents. Default is -1. greedy : bool, default is True Whether the use a heuristic in order to massive reduce computation and memory time, but without the guarantee of finding the best network. n_jobs : int The number of threads to use when learning the structure of the network. This parallelizes the creation of the parent graphs. structure : tuple, shape=(d,) The parents for each variable in this SCC
pomegranate.BayesianNetwork.generate_parent_graph()

Generate a parent graph for a single variable over its parents.

This will generate the parent graph for a single parents given the data. A parent graph is the dynamically generated best parent set and respective score for each combination of parent variables. For example, if we are generating a parent graph for x1 over x2, x3, and x4, we may calculate that having x2 as a parent is better than x2,x3 and so store the value of x2 in the node for x2,x3.

Parameters: X : numpy.ndarray, shape=(n, d) The data to fit the structure too, where each row is a sample and each column corresponds to the associated variable. weights : numpy.ndarray, shape=(n,) The weight of each sample as a positive double. Default is None. key_count : numpy.ndarray, shape=(d,) The number of unique keys in each column. pseudocount : double A pseudocount to add to each possibility. max_parents : int The maximum number of parents a node can have. If used, this means using the k-learn procedure. Can drastically speed up algorithms. If -1, no max on parents. Default is -1. parent_set : tuple, default () The variables which are possible parents for this variable. If nothing is passed in then it defaults to all other variables, as one would expect in the naive case. This allows for cases where we want to build a parent graph over only a subset of the variables. structure : tuple, shape=(d,) The parents for each variable in this SCC

format: plain