<jupyter_start><jupyter_text>* [Making packages in R using devtools](https://uoftcoders.github.io/studyGroup/lessons/r/packages/lesson/)    * `R CMD build bnlearn`* [Using devtools for Lazy People](http://rstudio-pubs-static.s3.amazonaws.com/2556_4e9f1c2af93b4683a19e2303a52bb2d5.html)* [Building and installing an R package](http://kbroman.org/pkg_primer/pages/build.html)* [devtools](https://www.rstudio.com/products/rpackages/devtools/)    * [github](https://github.com/r-lib/devtools)* [Can the code of a R package be changed? How?](https://www.biostars.org/p/146490/)    * `dev_mode`, `install_local`    * [load_all](https://rdrr.io/cran/devtools/man/load_all.html)    * `remove.packages("bnlearn")`<jupyter_code>%%R
ip <- as.data.frame(installed.packages()[,c(1,3:4)])
rownames(ip) <- NULL
ip <- ip[is.na(ip$Priority),1:2,drop=FALSE]
#print(ip, row.names=FALSE)
# /home/cs/local/install/Anaconda3-5.1.0-Linux-x86_64/lib/R/library/bnlearn
# %%R
# library("devtools")
# getwd()
# getLoadedDLLs()
# dyn.load("bnlearn/src/bnlearn.so")
# library.dynam("bnlearn", package = c("bnlearn_4.3"), lib.loc = c("bnlearn/src"))
# sessionInfo()
# devtools::build("/home/cs/workspaces/bnlearn/")
# %%R
# #devtools::install_local("../../bnlearn")
# load_all("/home/cs/workspaces/bnlearn/")
%load_ext autoreload
%autoreload 1
%aimport pybnl.bn
import locale
locale.setlocale(locale.LC_ALL, 'C')

import rpy2, rpy2.rinterface, rpy2.robjects, rpy2.robjects.packages, rpy2.robjects.lib, rpy2.robjects.lib.grid, \
    rpy2.robjects.lib.ggplot2, rpy2.robjects.pandas2ri, rpy2.interactive.process_revents, \
    rpy2.interactive, rpy2.robjects.lib.grdevices
# rpy2.interactive.process_revents.start()
rpy2.robjects.pandas2ri.activate()
%%R -o marks
library(bnlearn)
data(marks)
marks.head()<jupyter_output><empty_output><jupyter_text>* [networkx](https://github.com/networkx/networkx)* [tutorial](https://networkx.github.io/documentation/stable/tutorial.html)* [graphviz](https://github.com/xflr6/graphviz)* [graphviz.readthedocs](https://graphviz.readthedocs.io/en/stable/)* [examples/notebook.ipynb](https://github.com/xflr6/graphviz/blob/master/examples/notebook.ipynb)<jupyter_code>dg = nx.DiGraph()
# G.add_node(1)
dg.add_nodes_from(list(marks.columns))
dg.add_edges_from([
    ['STAT', 'ANL'],
    ['STAT', 'ALG'],
    ['ANL', 'ALG'],
    ['ALG', 'MECH'],
    ['ALG', 'VECT'],
    ['VECT', 'MECH'],
])
list(nx.connected_components(dg.to_undirected()))
nx.draw(dg, with_labels=True)
dg_dot = graphviz.Digraph(comment='disease-network')
for node in dg.nodes():
    dg_dot.node(node)

for edge in dg.edges():
    dg_dot.edge(edge[0],edge[1])

dg_dot
list(networkx.algorithms.dag.topological_sort(dg))
list(dg.nodes()),list(dg.edges())
nds = rpy2.robjects.StrVector(list(dg.nodes()))
connectivity_matrix = xr.DataArray(np.zeros((len(marks.columns), len(marks.columns))), 
             dims = ['from','to'], 
             coords= {'from': list(marks.columns), 'to': list(marks.columns)})
for edge in dg.edges():
    connectivity_matrix.loc[{'from': edge[0], 'to': edge[1]}] = 1
np_connectivity_matrix = np.array(connectivity_matrix.values)
np_connectivity_matrix
%%R -i np_connectivity_matrix -i nds
dag = empty.graph(nds)
amat(dag) = np_connectivity_matrix
dag
%%R
moral(dag)
%%R
node.ordering(dag)
%%R
nbr(dag, "ANL")
%%R
mb(dag, "ANL")
%%R
score(dag,data=marks,type="loglik-g")
%%R
vstructs(dag)
%%R
vstructs(dag, moral=TRUE)
%%R
cpdag(dag)<jupyter_output><empty_output><jupyter_text>[Principled Computational Methods for the Validation and Discovery of Genetic Regulatory Networks](http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.4.6170) by <NAME>, 2001<jupyter_code>%%R -o dmarks1 -o dmarks2
dmarks1 = discretize(marks, breaks=3, method="interval")
dmarks2 = discretize(marks, breaks=3, method="hartemink", ibreaks=5, idisc="quantile")
%%R
levels(dmarks1$MECH)
%%R
levels(dmarks2$MECH)
dmarks2.info()
display_side_by_side(dmarks1.head(), dmarks2.head())
fig=plt.figure(figsize=(15, 8), dpi= 80, facecolor='w', edgecolor='k')
ax  = plt.subplot(1,1,1)
dmarks1['MECH'].groupby(dmarks1['MECH']).agg(['count']).plot.bar(ax=ax);
import dsbasics.bin
x = marks['MECH'].values
X = x.reshape(-1,1)

bbt1 = dsbasics.bin.BayesianBlocksBinTransformer(p0=0.05)
bbt1.fit(X)

bbt2 = dsbasics.bin.DecisionTreeBinTransformer(max_leaf_nodes=3)
bbt2.fit(X)


fig=plt.figure(figsize=(15, 8), dpi= 80, facecolor='w', edgecolor='k')
ax  = plt.subplot(1,1,1)

ax.hist(x,bins=40, histtype='stepfilled', alpha=0.2, density=True)
ax.hist(x, bins=bbt1.bins_[0], histtype='step', color='black', density=True)
ax.hist(x, bins=bbt2.bins_[0], histtype='step', color='green', density=True)
ax.set_ylabel('P(MECH)')
# ax.set_xlim(0.0,400.0)
ax.set_xlabel('MECH')
x = dmarks1['MECH']
c = pd.Categorical(x)
c.categories
%%R
write.csv(marks, file = "marks.csv")
head(marks)
pd_marks = pd.read_csv('marks.csv', index_col=0).astype(np.float64)
pd_marks.head()
# %%R -i pd_marks
# dmarks2_ = discretize(pd_marks, breaks=3, method="hartemink", ibreaks=5, idisc="quantile")
# head(dmarks2_)
df = pybnl.bn.discretize(pd_marks)
df.head()
df.dtypes
df['MECH'].cat.categories
df['MECH'].groupby(df['MECH']).agg(['count']).plot.bar()
df['MECH'].head()
bn = pybnl.bn.NetAndDataDiscreteBayesNetwork(df, rnet=pybnl.bn.empty_graph(df.columns))
bn.fit()
bn = pybnl.bn.NetAndDataDiscreteBayesNetwork(df, dg=dg)
bn.fit()
bn.model_string
print(bn.rnet)
print(bn.rfit)
print(bn.grain)
df['STAT'].cat.categories[0]
evidence = dict(ALG=df['ALG'].cat.categories[0])
nodes_to_query = ['STAT']
bn.exact_query(evidence, nodes_to_query, only_python_result=True)
bn.write_net('marks.net')
# TODO
# Convert network from R to python xarray and pandas and allow for save and load
# Implement structure learning
# Ask question about discretization score available in bnlearn discretize?
# Plot discretization levels
# Implement error checking
#  Ensure that queries only use valid variables and values for these variables
# Implement unit tests
# Implement approximate queries
# Implement model averaging, cross validation, ...
# Implement ames house price example
#  Arc strengths
# Implement own problem with LLT dependencies
# Plot graphs as equivalence structures (cpdag) and use cextend for consistent extension
# Implement sampling from a network
rfit = bn.rfit
%%R -i rfit
# nodes(rfit)
# colnames, rownames, dimnames
# dimnames(rfit$STAT$prob)[1]
# dimnames(rfit$ANL$prob)[1]
dimnames(rfit$ANL$prob)
%%R -o anl_prob
anl_prob = rfit$ANL$prob
rfit.rx('ANL')
rpy2.robjects.pandas2ri.ri2py(anl_prob)
anl_prob
anl_prob.names.rx('ANL').names[0]
anl_prob.names
anl_prob.names
rfit.rx('ANL')[0].rx('prob')[0].names.rx('ANL')[0]
rpy2.robjects.pandas2ri.ri2py(rfit.rx('ANL')[0].rx('prob')[0])
rfit.rx('ANL')[0].rx('prob')[0].names.rx('ANL')[0]
rfit.rx('ANL')[0].rx('prob')[0].names
rfit.rx('STAT')[0].rx('prob')[0].names[0][0]
ds,lpd = pybnl.bn.convert_to_xarray_dataset(bn.rfit)
ds
lpd['cptALG']
ds['cptSTAT']
ds['cptANL']
ds['cptANL'].loc[{'STAT': '[9,29.4]'}]
ds['cptALG']
ds['cptALG'].loc[{'ANL': '(9,47]', 'STAT': '[9,29.4]'}]
ds['cptALG'].loc[{'ANL': '(47,60.2]', 'STAT': '(29.4,44]'}]
ds['cptALG'].values[:,1,1]
ds['cptALG'].values
ds['cptALG'].loc[{'ANL': '(60.2,70]', 'STAT': '[9,29.4]'}]
ds['cptALG'].transpose('ANL', 'STAT', 'ALG').stack(idx=['ANL', 'STAT', 'ALG']).to_pandas()#.reset_index()
ds.to_netcdf('marks.nc')<jupyter_output><empty_output><jupyter_text>* [NetCDF](https://www.unidata.ucar.edu/software/netcdf/software.html) viewers* use `ncdump` from the `netcdf-bin` package<jupyter_code>tmp = xr.open_dataset('marks.nc', autoclose=True)
tmp
tmp.equals(ds)
def tmpfn():
    for ar in tmp.data_vars:
        return ar
tmpfn()
bn.write_netcdf('marks.nc')
bn1 = pybnl.bn.bnnet_from_netcdf_file('marks.nc')
bn1.to_xrds().equals(bn.to_xrds())
bn.to_xrds()
tmp = xr.open_dataset('marks.nc',autoclose=True)
tmp.equals(bn.to_xrds())
tmp
bn1.to_xrds().equals(tmp)
bn1.to_dtcpm_dict()['cptALG'].sort_values(['ANL', 'STAT', 'ALG'])
p1 = bn1.to_dtcpm_dict()['cptALG'].sort_values(['ANL', 'STAT', 'ALG'])['p']
p1.head()
bn.to_dtcpm_dict()['cptALG'].sort_values(['ANL', 'STAT', 'ALG'])
p = bn.to_dtcpm_dict()['cptALG'].sort_values(['ANL', 'STAT', 'ALG'])['p']
p.head()
p.values == p1.values
bn.to_xrds()['cptALG']
np.all(bn1.to_xrds()['cptALG'] == bn.to_xrds()['cptALG'])
tmp = bn.rfit
%%R -i tmp
arcs(tmp)
%%R
arcs(cpdag(tmp))
%%R
amat(tmp)
%%R
amat(cpdag(tmp))
dg = nx.DiGraph()
# G.add_node(1)
dg.add_nodes_from(list(marks.columns))
dg.add_edges_from([
    ['STAT', 'ANL'],
    ['STAT', 'ALG'],
    ['ANL', 'ALG'],
    ['ALG', 'MECH'],
    ['ALG', 'VECT'],
    ['VECT', 'MECH'],
])
dg
dg_dot = graphviz.Digraph(comment='disease-network')
for node in dg.nodes():
    dg_dot.node(node)

for edge in dg.edges():
    dg_dot.edge(edge[0],edge[1])

print(dg_dot.source)
dg_dot = graphviz.Digraph(comment='disease-network')
for node in dg.nodes():
    dg_dot.node(node)

for edge in dg.edges():
    dg_dot.edge(edge[0],edge[1], dir='none')

dg_dot
# rm,_,_ = pybnl.bn.cpdag(bn)
# print(rm.dimnames[1])
# rm.rx(1,2)
pybnl.bn.cpdag(bn)
# for i in range(12):
#     print(rm.rx(i+1,1), rm.rx(i+1,2))
pybnl.bn.dag(bn)
pybnl.bn.score(bn, bn.df)
pybnl.bn.vstructs(bn.rnet)
tmp1 = bn.rfit
tmp2 = bn.rnet
%%R -i tmp1 -i tmp2
#print(tmp2)
names(tmp1)
print(pybnl.bn.rfit2rnet(bn.rfit))
pybnl.bn.dot(*pybnl.bn.dag(bn))
pybnl.bn.dot(*pybnl.bn.cpdag(bn))
print(pybnl.bn.drop_arc(bn.rnet, frm='STAT', to='ANL'))
pybnl.bn.rnet2cpdag(pybnl.bn.drop_arc(bn.rnet, frm='STAT', to='ANL'))
pybnl.bn.dot(*pybnl.bn.rnet2cpdag(pybnl.bn.drop_arc(bn.rnet, frm='STAT', to='ANL')),engine='dot')
tmp1 = bn.rfit
tmp2 = bn.rnet
%%R -i tmp1 -i tmp2
tmp2_da = drop.arc(tmp2, from="STAT", to="ANL")
# print(tmp2)
cpdag(tmp2_da)
%%R -o rdf_lt
data(learning.test)
rdf_lt = learning.test<jupyter_output><empty_output><jupyter_text>* [mapping-categorical-data-in-pandas](http://benalexkeen.com/mapping-categorical-data-in-pandas/)* [CategoricalDtype](https://pandas.pydata.org/pandas-docs/stable/categorical.htmlcategoricaldtype)<jupyter_code>#df_lt = rpy2.robjects.pandas2ri.ri2py(rdf_lt)
df_lt = rdf_lt

ct1 = pd.api.types.CategoricalDtype(['a', 'b', 'c'], ordered=True)
ct2 = pd.api.types.CategoricalDtype(['a', 'b'], ordered=True)

for c in 'ABCDE':
    df_lt[c] = df_lt[c].astype(ct1)
df_lt['F'] = df_lt['F'].astype(ct2)

df_lt.info()<jupyter_output><empty_output><jupyter_text>* [Fundamentals of Structure Learning](http://www.bnlearn.com/about/teaching/slides-bnshort.pdf) p. 157    * bnlearn implements several constraint-based algorithms, each with its own function: gs(), iamb(), mmpc(), si.hiton.pc(), etc.<jupyter_code>%%R -o cpdaghiton
cpdaghiton = si.hiton.pc(rdf_lt, test = "mc-mi", undirected = FALSE)
cpdaghiton<jupyter_output><empty_output><jupyter_text>[networks](http://www.bnlearn.com/documentation/networks/)<jupyter_code>pybnl.bn.dot(*pybnl.bn.rnet2dag(cpdaghiton),engine='dot')
%%R -o daghc
daghc = hc(rdf_lt, score = "bic", iss=1, restart = 10, perturb = 5, start = random.graph(names(rdf_lt)))
pybnl.bn.dot(*pybnl.bn.rnet2dag(daghc),engine='dot')<jupyter_output><empty_output><jupyter_text>* [Fundamentals of Structure Learning](http://www.bnlearn.com/about/teaching/slides-bnshort.pdf) p. 203    * bnlearn implements several constraint-based algorithms, each with its own function: gs(), iamb(), mmpc(), si.hiton.pc(), etc.<jupyter_code>%%R -o dagrsmax2
dagrsmax2 = rsmax2(rdf_lt, 
                   restrict = "si.hiton.pc", restrict.args = list(test = "x2", alpha = 0.01), 
                   maximize = "tabu", maximize.args = list(score = "bic", tabu = 10))
pybnl.bn.dot(*pybnl.bn.rnet2dag(dagrsmax2),engine='dot')
%%R -o dagmmhc
# dagmmhc = rsmax2(rdf_lt, restrict = "mmpc", maximize = "hc")
dagmmhc = mmhc(rdf_lt)
pybnl.bn.dot(*pybnl.bn.rnet2dag(dagmmhc),engine='dot')<jupyter_output><empty_output><jupyter_text>Model comparison: DAGs as Multivariate Trinomials[Fundamentals of Structure Learning](http://www.bnlearn.com/about/teaching/slides-bnshort.pdf) p. 214<jupyter_code># %%R
# bnlearn:::check.data(ldmarks, allow.levels = TRUE, allow.missing = TRUE, warn.if.no.missing = TRUE)#, 
# # ldmarks$LAT = as.numeric(ldmarks$LAT)
%%R
head(dmarks2)
%%R
ldmarks = data.frame(dmarks2, LAT = factor(rep(NA, nrow(dmarks2)), levels = c("A", "B")))
head(ldmarks)
%%R
imputed = ldmarks
imputed$LAT = sample(factor(c("A", "B")), nrow(dmarks2), replace = TRUE)
head(imputed)
%%R
head(ldmarks)
%%R -o fitted
fitted = bn.fit(empty.graph(names(ldmarks)), imputed)
fitted$LAT = array(c(0.5, 0.5), dim = 2, dimnames = list(c("A", "B")))
# three iterations of structural EM.
for (i in 1:15) {
    # expectation step.
    imputed = impute(fitted, ldmarks, method = "bayes-lw")
    # maximisation step (forcing LAT to be connected to the other nodes).
    dag = hc(imputed, whitelist = data.frame(from = "LAT", to = names(dmarks2)))
    fitted.new = bn.fit(dag, imputed, method = "bayes")


    if (isTRUE(all.equal(fitted, fitted.new))) {
      print(c("i: ", i))
      break
    } else {
      fitted = fitted.new
    }
}#FOR
print(score(dag,data=imputed,type="bic"))
modelstring(fitted)
%%R
table(imputed$LAT)
%%R
latent = factor(c(rep("A", 44), "B", rep("A", 7), rep("B", 36)))
print(table(latent))
# head(data.frame(dmarks2, LAT=latent))
%%R -o latent -o imputed
list(latent, imputed$LAT)
tmps1 = pd.Series(rpy2.robjects.pandas2ri.ri2py(latent))
tmps1.head()
imputed.LAT.head()
tmps2 = imputed.LAT.astype(pd.api.types.CategoricalDtype(['A', 'B']))
tmps2.head()
list(pd.Series(rpy2.robjects.pandas2ri.ri2py(latent)).astype(str).iloc[:10])
list(imputed.LAT.iloc[:10])
np.finfo(np.float64).eps
import sklearn, sklearn.metrics, sklearn.metrics.cluster.supervised
first_n = 80
tmps1.iloc[:first_n].values
tmps1.iloc[:first_n].value_counts()
tmps2.iloc[:first_n].value_counts()
sklearn.metrics.cluster.supervised.contingency_matrix(tmps1.iloc[:first_n], tmps2.iloc[:first_n], sparse=False)
tmpcm = pybnl.bn.contingency_matrix(tmps1.iloc[:first_n], tmps2.iloc[:first_n], sparse=False)
tmpcm
tmpcmp = tmpcm / tmpcm.sum() + np.finfo(np.float64).eps
tmpcmp
lxr_tmpcmp = xr.DataArray(tmpcmp, dims=['tmps1', 'tmps2'], coords={'tmps1': ['A', 'B'], 'tmps2': ['A', 'B'],})
lxr_tmpcmp
lxr_tmpcmp.loc[dict(tmps1='A',tmps2='A')]
lxr_tmpcmp.loc[dict(tmps1='A')]
lxr_tmpcmp.sum(dim=['tmps1'])
lxr_ps1 = lxr_tmpcmp.sum(dim=['tmps2'])
lxr_ps1
lxr_ps2 = lxr_tmpcmp.sum(dim=['tmps1'])
lxr_ps2
np.set_printoptions(suppress=True)
(lxr_tmpcmp/lxr_ps1/lxr_ps2)
float((lxr_tmpcmp*np.log(lxr_tmpcmp/lxr_ps1/lxr_ps2)).sum())
sklearn.metrics.mutual_info_score(tmps1.iloc[:first_n], tmps2.iloc[:first_n], contingency=tmpcm)
sklearn.metrics.mutual_info_score(tmps1.iloc[:first_n], tmps2.iloc[:first_n])
sklearn.metrics.adjusted_mutual_info_score(tmps1.iloc[:first_n], tmps2.iloc[:first_n])
sklearn.metrics.normalized_mutual_info_score(tmps1.iloc[:first_n], tmps2.iloc[:first_n])
'{0:.20f}'.format(pybnl.bn.relative_mutual_information(tmps1.iloc[:first_n], tmps2.iloc[:first_n]))
pybnl.bn.relative_mutual_information(tmps1, tmps2)
sklearn.metrics.normalized_mutual_info_score(tmps1, tmps2)
pybnl.bn.dot(*pybnl.bn.rnet2cpdag(pybnl.bn.rfit2rnet(fitted)))#,engine='dot'
%%R -o fitted
fitted = bn.fit(empty.graph(names(ldmarks)), imputed)
fitted$LAT = array(c(0.5, 0.5), dim = 2, dimnames = list(c("A", "B")))
# three iterations of structural EM.
for (i in 1:15) {
    # expectation step.
    imputed = impute(fitted, ldmarks, method = "bayes-lw")
    # maximisation step (forcing LAT to be connected to the other nodes).
    # dag = hc(imputed, whitelist = data.frame(from = "LAT", to = names(dmarks2)))
    dag = mmhc(imputed, whitelist = data.frame(from = "LAT", to = names(dmarks2)))
    fitted.new = bn.fit(dag, imputed, method = "bayes")


    if (isTRUE(all.equal(fitted, fitted.new))) {
      print(c("i: ", i))
      break
    } else {
      fitted = fitted.new
    }
}#FOR
print(score(dag,data=imputed,type="bic"))
modelstring(fitted)
%%R
table(imputed$LAT)
%%R
print(table(latent))
%%R
list(latent, imputed$LAT)
# mutual information between latent and imputed$LAT?
pybnl.bn.dot(*pybnl.bn.rnet2cpdag(pybnl.bn.rfit2rnet(fitted)))#,engine='dot'
# %%R
# # structural expectation-maximization.
# structural.em.2 = function(x, maximize = "hc", maximize.args = list(), fit = "mle", fit.args = list(), impute, impute.args = list(), return.all = FALSE, start = NULL, max.iter = 5, debug = FALSE) {

#   ntests = 0
#   print("###")

#   # check the data are there.
#   # data.info =  bnlearn:::check.data(x, allow.levels = TRUE, allow.missing = TRUE, warn.if.no.missing = TRUE, stop.if.all.missing = !is(start, "bn.fit"))
#   data.info =  bnlearn:::check.data(x, allow.levels = TRUE, allow.missing = TRUE, warn.if.no.missing = TRUE)
#   print(data.info)#$complete.nodes
#   print("###")

#   # check the max.iter parameter.
#   max.iter = bnlearn:::check.max.iter(max.iter)
#   print("###")
#   # check debug and return.data.
#   bnlearn:::check.logical(debug)
#   bnlearn:::check.logical(return.all)
#   print("###")

#   # check the arguments used for structure learning as in rsmax2().
#   bnlearn:::check.learning.algorithm(algorithm = maximize, class = "score")

#   print("###")
    
#   critical.arguments = c("x", "heuristic", "start", "debug")
#   bnlearn:::check.unused.args(intersect(critical.arguments, names(maximize.args)), character(0))
#   maximize.args[critical.arguments] = list(x = NULL, heuristic = maximize, start = NULL, debug = debug)

#   print("###")
    
#   # check the arguments used for parameter learning as in bn.cv().
#   bnlearn:::check.fitting.method(method = fit, data = x)
#   fit.args = bnlearn:::check.fitting.args(method = fit, network = NULL, data = x,  extra.args = fit.args)

#   print("###")
    
#   # check the arguments used for imputations.
#   impute = bnlearn:::check.imputation.method(impute, x)
#   impute.args = bnlearn:::check.imputation.extra.args(impute, impute.args)

#   # if there is no preseeded network, use an empty one.
#   print("###")
#   if (is.null(start)) {
#     print("is.null(start)")

#     dag = empty.graph(nodes = names(x))
#     fitted = bnlearn:::bn.fit.backend(dag, data = x, method = fit, extra.args = fit.args) # , data.info = data.info

#   }#THEN
#   else {
#     print("is.null(start): else")

#     print("---")
#     # check start's class.
#     bnlearn:::check.bn.or.fit(start)
#     print("---")
#     # check the preseeded network against the data set.
#     if (is(start, "bn")) {

#       print("--- is(start, bn)")
#       bnlearn:::check.bn.vs.data(start, x)
#       # check the preseeded network is valid for the model assumptions.
#       bnlearn:::check.arcs.against.assumptions(start$arcs, x, maximize.args$score)

#       dag = start
#       fitted = bnlearn:::bn.fit.backend(start, data = x, method = fit, extra.args = fit.args) # , data.info = data.info

#     }#THEN
#     else if (is(start, "bn.fit")) {
#       print("--- is(start, bn.fit)")

#       bnlearn:::check.fit.vs.data(start, x)
#       # fitted networks are necessarily valid for the model assumptions.

#       dag = bn.net(start)
#       fitted = start

#     }#THEN

#   }#ELSE

#   print("###.")
#   # initialize the algorithm.
#   if (debug) {

#     cat("* initializing the network to perform the first imputation.\n")
#     cat("* network structure:\n")
#     print(dag)
#     cat("* fitted parameters:\n")
#     print(fitted)

#   }#THEN

#   # from this point on the data are complete.
#   print(c("names(x) ", names(x)))
#   print(data.info)#$complete.nodes
#   print("###..")
#   #data.info$complete.nodes[names(x)] = TRUE

#   print("###...")
#   for (i in seq(max.iter)) {

#     if (debug) {

#       cat("----------------------------------------------------------------\n")
#       cat("* iteration", i, ", expectation step .\n")

#     }#THEN

#     # expectation step.
#     complete = bnlearn:::impute.backend(fitted = fitted, data = x, method = impute, extra.args = impute.args, debug = debug)

#     if (debug) {

#       cat("----------------------------------------------------------------\n")
#       cat("* iteration", i, ", maximization step .\n")

#     }#THEN

#     # maximization step, structure learning (starting from the previous network).
#     maximize.args$x = complete
#     maximize.args$start = dag
#     greedy.search.fn = get("greedy.search", asNamespace("bnlearn"))
#     dag = do.call(greedy.search.fn, maximize.args, quote=TRUE)

#     # maximization step, parameter learning.
#     fitted.new = bnlearn:::bn.fit.backend(dag, data = complete, method = fit, extra.args = fit.args) # , data.info = data.info

#     if (debug) {

#       cat("----------------------------------------------------------------\n")
#       cat("* fitted parameters:\n")
#       print(fitted)

#     }#THEN

#     # update the number of model comparisons.
#     ntests = ntests + dag$learning$ntests

#     # a rudimentary stopping rule (the fitted network must be the same as in
#     # the previous iteration).
#     if (isTRUE(all.equal(fitted, fitted.new)))
#       break
#     else
#       fitted = fitted.new

#   }#FOR

#   # set the metadata.
#   dag$learning$algo = "sem"
#   dag$learning$maximize = maximize
#   dag$learning$impute = impute
#   dag$learning$fit = fit
#   dag$learning$ntests = ntests

#   if (return.all)
#     invisible(list(dag = dag, imputed = complete, fitted = fitted))
#   else
#     invisible(dag)

# }#STRUCTURAL.EM
%%R
# ldmarks[1,"LAT"] = "A"
# ldmarks[2,"LAT"] = "B"
head(ldmarks)
%%R
imputed = ldmarks
imputed$LAT = sample(factor(c("A", "B")), nrow(dmarks2), replace = TRUE)

dag = empty.graph(nodes = names(ldmarks))
# fitted = bnlearn:::bn.fit.backend(dag, imputed)
fitted = bn.fit(dag, imputed)
fitted$LAT = array(c(0.5, 0.5), dim = 2, dimnames = list(c("A", "B")))
# complete = impute.backend(fitted = fitted, data = x, method = impute, extra.args = impute.args, debug = debug)
# %%R
# fitted<jupyter_output><empty_output><jupyter_text>* [bayesnetRtutorial](https://github.com/jacintoArias/bayesnetRtutorial/blob/master/index.Rmd) has a parametric em* [bnlearn-error-in-structural-em](https://stackoverflow.com/questions/47970302/bnlearn-error-in-structural-em)<jupyter_code># %%R
# bnlearn:::check.data(ldmarks, allow.levels = TRUE, allow.missing = TRUE, warn.if.no.missing = TRUE, stop.if.all.missing = !is(fitted, "bn.fit"))
# bnlearn:::check.data(ldmarks, allow.levels = TRUE, allow.missing = TRUE, warn.if.no.missing = TRUE, stop.if.all.missing = FALSE)#, 
# ldmarks$LAT = as.numeric(ldmarks$LAT)
# %%R
# trace(structural.em, edit="emacs")
#trace(bnlearn:::check.data)
%%R -o bn2
# structural.em(x, maximize = "hc", maximize.args = list(), fit = "mle", fit.args = list(), impute, impute.args = list(), return.all = FALSE, start = NULL, max.iter = 5, debug = FALSE)
# imputed2 = ldmarks
# imputed2$LAT = sample(factor(c("A", "B")), nrow(dmarks2), replace = TRUE)
# r = structural.em.2(ldmarks, fit = "bayes", impute="bayes-lw", start=fitted, return.all = TRUE)
# r = structural.em.2(ldmarks, fit = "bayes", impute="bayes-lw", return.all = TRUE)
# r = structural.em(ldmarks, fit = "bayes", impute="bayes-lw", return.all = TRUE)
r = structural.em(ldmarks, fit = "bayes", impute="bayes-lw", start=fitted, maximize.args=list(whitelist = data.frame(from = "LAT", to = names(dmarks2))), return.all = TRUE)
bn2 = r$dag
print(modelstring(r$dag))
%%R
table(r$imputed$LAT)
pybnl.bn.dot(*pybnl.bn.rnet2dag(bn2))#,engine='dot'
# %%R
# latent = factor(c(rep("A", 44), "B", rep("A", 7), rep("B", 36)))
# print(table(latent))
# head(data.frame(dmarks2, LAT=latent))
# %%R
# list(latent, r$imputed$LAT, imputed$LAT)
%%R -o rdf_lt
data(learning.test)
rdf_lt = learning.test
df_lt = rdf_lt

ct1 = pd.api.types.CategoricalDtype(['a', 'b', 'c'], ordered=True)
ct2 = pd.api.types.CategoricalDtype(['a', 'b'], ordered=True)

for c in 'ABCDE':
    df_lt[c] = df_lt[c].astype(ct1)
df_lt['F'] = df_lt['F'].astype(ct2)

df_lt.info()
cbnet = pybnl.bn.ConstraintBasedNetFromDataDiscreteBayesNetwork(df_lt)
cbnet.fit()<jupyter_output><empty_output><jupyter_text>[networks](http://www.bnlearn.com/documentation/networks/)<jupyter_code>cbnet.structure().dot()
cbnet.structure().cpdag().dot()
sbnet = pybnl.bn.ScoreBasedNetFromDataDiscreteBayesNetwork(df_lt)
sbnet.fit()
#sbnet.algorithm
sbnet.structure().cpdag().dot()
sbnet.structure().cpdag().vstructs()
#bn#.structure().cpdag().vstructs()
hnet1 = pybnl.bn.HybridScoreAndConstainedBasedNetFromDataDiscreteBayesNetwork(df_lt)
hnet1.fit()
hnet1.structure().cpdag().dot()
hnet2 = pybnl.bn.HybridScoreAndConstainedBasedNetFromDataDiscreteBayesNetwork(df_lt, algorithm='rxmax2_sihitonpc_tabu')
hnet2.fit()
hnet2.structure().cpdag().dot()<jupyter_output><empty_output><jupyter_text>[Online Density Estimates: A Probabilistic Condensed Representation of Data for Knowledge Discovery](https://publications.ub.uni-mainz.de/theses/volltexte/2017/100001664/pdf/100001664.pdf)<jupyter_code># mutual information for imputed data series as metric
# score method in BayesNetworkBase
# functionality for parametric.em and structural.em / imputing missing data
# naive bayes
# BayesNetworkStructure to/from model string
# parameter classes for constraint/score based algorithms to be included in rsmax2 hybrid / class -> dict
# example for ames house prices plus cross validation and comparison to random forest
# discretization sklearn transformer
# adapt CustomDiscreteBayesNetwork, NetAndDataDiscreteBayesNetwork to fit to new sklearn fit() methodology
# bayes factor: model comparison
# arc strength
# mixture models
lc = pd.Categorical([i%3 for i in range(9)], categories=[0,1,2], ordered=True)
lds1 = pd.Series(lc)

lds2 = pd.Series([i%3 for i in range(9)])

lds_ct = pd.api.types.CategoricalDtype([0,1,2], ordered=True)
lds2 = lds2.astype(lds_ct)

lds1.equals(lds2)
lds2[:] = np.nan
ldf = pd.DataFrame({'c1': lds1, 'c2': lds2})
pybnl.bn.identify_latent_variables(ldf)
ldf.c2.dtype.ordered
pybnl.bn.levels_of_latent_variable(ldf,'c2')
pd.Categorical([np.nan], categories=[0,1,2])
dmarks2.shape
dmarks2.info()
df = pybnl.bn.discretize(pd_marks)
df.info()
#print(pybnl.bn.pydf_to_factorrdf(df))
bn = pybnl.bn.NetAndDataDiscreteBayesNetwork(df, rnet=pybnl.bn.empty_graph(df.columns))
bn.fit()
ldmarks = df.copy()
pybnl.bn.augment_df_with_latent_variable(ldmarks, 'LAT', 3)
print(pybnl.bn.levels_of_latent_variable(ldmarks,'LAT'))
ldmarks.head()
ldmarks.info()
ldmarks['LAT'].dtype
%%R -o tmp
print(list(c("A", "B")))
tmp = array(c(0.5, 0.5), dim = 2, dimnames = list(c("A", "B")))
print(tmp)
tmp
count=3
dimnames = rpy2.robjects.r['list'](rpy2.robjects.StrVector(['l000','l001','l002']))
# tmp_ = rpy2.robjects.r['array'](np.full([count], float(1/count)), dim=rpy2.robjects.IntVector([count]), dimnames=dimnames)
tmp_ = rpy2.robjects.r['array'](np.full([count], float(1/count)), dim=count, dimnames=dimnames)
print(dimnames)
print(tmp_)
tmp_
sembn = pybnl.bn.StructuralEMNetFromDataDiscreteBayesNetwork(ldmarks)
tmp = sembn.fit()
#print(tmp.rfit)
tmp
tmp.structure().dot()
tmp.imputed_
# tmp.rfit[5][3] = tmp_
# tmp.rfit
ldmarks.head()
%%R
factor(c('a','b',NA))
rpy2.rinterface.NA_Integer
rpy2.robjects.DataFrame
ldmarks_ = ldmarks.copy()
ldmarks_.loc[:,'LAT'] = np.full([88],np.nan)
#ldmarks_.loc[:,'LAT'] = np.full([88],rpy2.rinterface.NA_Character)
ldmarks_['LAT'] = ldmarks_['LAT'].astype(ldmarks['LAT'].dtype)
tmp_ = pybnl.bn.pydf_to_factorrdf(ldmarks_)
# ldmarks_['LAT']
%%R -i tmp_
str(tmp_)
#unclass(tmp_$LAT)
# rpy2.robjects.r['unclass'](tmp_[5])
# pd.Series()
# rpy2.robjects.FloatVector(tmp_[5])
# rpy2.robjects.pandas2ri.ri2py(rpy2.robjects.r['unclass'](tmp_[5]))
rpy2.robjects.r('all')(rpy2.robjects.r('is.na')(tmp_[5]))[0]
pybnl.bn.factorrdf_to_pydf(pybnl.bn.pydf_to_factorrdf(ldmarks_)).info()
# https://stackoverflow.com/questions/30334385/display-svg-in-ipython-notebook-from-a-function
# http://nbviewer.jupyter.org/github/ipython/ipython/blob/3.x/examples/IPython%20Kernel/Rich%20Output.ipynb<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.19587075535835688, 'fraction_numerical': 0.01907181450602927, 'mean_word_length': 5.057217165149545, 'pattern_counts': {'":': 0, '<': 53, '<?xml version=': 0, '>': 51, 'https://': 16, 'lorem ipsum': 0, 'www.': 8, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '538432', 'license': 'Apache-2.0', 'n_tokens_mistral': 11721, 'n_tokens_neox': 11322, 'n_words': 1871, 'path': 'examples/010-2018-06-05-marks-scratchpad.ipynb', 'repo_name': 'cs224/pybnl', 'size': '27965'}