# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

#
# # Project 1 : Determining Probabilities of Handwriting Formations using PGMs
# ## UBIT Name : mgosi
# ## Person No   : 50290934
# ## Name : <NAME>

# + init_cell=true
import pandas as pd
import numpy as np
from pgmpy.estimators import K2Score
from pgmpy.models import BayesianModel
from IPython.display import display
from pgmpy.inference import BeliefPropagation
# -

# # Part 1 : Evaluate Correlations and Independencies
# We determine how the variables are related to eachother. This is tested by comparing the Joint Probability with the product of the probabilities of the corresponding variables.
# ## Retrieving the Probabilities into a Dataframe

# + init_cell=true
CondProbDist = {}
CondProbDist
for i in range(6):
    CondProbDist[i] = pd.read_csv('Table'+str(i+3)+'.csv')
    CondProbDist[i].replace(regex = True, inplace = True, to_replace = r'% *\(*[0-9]*\)*', value = r'')
MargProb = pd.read_csv('Table2.csv')
MargProb.replace(regex = True, inplace = True, to_replace = r'% *\(*[0-9]*\)*', value = r'')

# -

for i in range(6):
    print('='*76)
    print ('CPD '+str(i+1))
    print('='*76)
    display(CondProbDist[i])
    print('-'*76)

x_len = [0,4,5,3,4,4,5]
MargProb


# ### Here we define a function to calculate the product of probabilites of corresponding variables using the Marginal Probabilities of the Variables.

# + init_cell=true
def MulMargProb(x, y):
    x_val = MargProb.iloc[:,x].values
    y_val = MargProb.iloc[:,y].values
    prodXY = []
    for i in range(len(x_val)):
        pXY = []
        for j in range(len(y_val)):
            pXY.append(float(x_val[i])*float(y_val[j]))
        prodXY.append(pXY)
    return prodXY


# -

# ### Joint Probability of the Variables are calculated using the CPD Tables provided.

# + code_folding=[] init_cell=true
def jointProb(p, c):
    jP = []
    for i in range(len(c)):
        tempJp = []
        for j in range(p.shape[0]):
            tempJp.append(float(p[j][i]) * float(c[i]))
        jP.append(tempJp)
    jP = np.asarray(jP)
    return jP
        #print(np.asarray(jp).T)


# -

# ### Comparing both by subtracting

# + init_cell=true
def check_independence():
    ind_dic = {}
    ind_dic['12'] = abs(jointProb(p12,c1) - m12).sum()/10000
    ind_dic['14'] = abs(jointProb(p14,c1) - m14).sum()/10000
    ind_dic['16'] = abs(jointProb(p16,c1) - m16).sum()/10000
    
    ind_dic['23'] = abs(jointProb(p23,c2) - m23).sum()/10000
    ind_dic['25'] = abs(jointProb(p25,c2) - m25).sum()/10000
    
    ind_dic['32'] = abs(jointProb(p32,c3) - m32).sum()/10000
    ind_dic['35'] = abs(jointProb(p35,c3) - m35).sum()/10000
    ind_dic['36'] = abs(jointProb(p36,c3) - m36).sum()/10000
    
    ind_dic['41'] = abs(jointProb(p41,c4) - m41).sum()/10000
    ind_dic['42'] = abs(jointProb(p42,c4) - m42).sum()/10000
    ind_dic['46'] = abs(jointProb(p46,c4) - m46).sum()/10000
    
    ind_dic['52'] = abs(jointProb(p52,c5) - m52).sum()/10000
    ind_dic['53'] = abs(jointProb(p53,c5) - m53).sum()/10000

    ind_dic['61'] = abs(jointProb(p61,c6) - m61).sum()/10000
    ind_dic['62'] = abs(jointProb(p62,c6) - m62).sum()/10000
    ind_dic['63'] = abs(jointProb(p63,c6) - m63).sum()/10000
    ind_dic['64'] = abs(jointProb(p64,c6) - m64).sum()/10000
    
    return (ind_dic)

# -



# + code_folding=[] init_cell=true
probDist = CondProbDist[0]
p12 = probDist.iloc[1:6,1:].values
c1 = probDist.iloc[0,1:].values
p14 = probDist.iloc[6:10,1:].values
p16 = probDist.iloc[10:,1:].values

m12 = MulMargProb(1,2)[:4]
m14 = np.asarray(MulMargProb(1,4)[:4]).T[:4].T.tolist()
m16 = MulMargProb(1,6)[:4]



# + code_folding=[] init_cell=true
probDist = CondProbDist[1]
#display(probDist)
c2 = probDist.iloc[0,1:].values
p23 = probDist.iloc[1:4,1:].values
p25 = probDist.iloc[4:,1:].values

m23 = np.asarray(MulMargProb(2,3)).T[:3].T.tolist()
m25 = np.asarray(MulMargProb(2,5)).T[:4].T.tolist()


# + init_cell=true
probDist = CondProbDist[2]
#display(probDist)
c3 = probDist.iloc[0,1:].values
p32 = probDist.iloc[1:6,1:].values
p35 = probDist.iloc[6:10,1:].values
p36 = probDist.iloc[10:,1:].values

m32 = MulMargProb(3,2)[:3]
m35 = np.asarray(MulMargProb(3,5)[:3]).T[:4].T.tolist()
m36 = MulMargProb(3,6)[:3]



# + init_cell=true
probDist = CondProbDist[3]
#display(probDist)
c4 = probDist.iloc[0,1:].values
p41 = probDist.iloc[1:5,1:].values
p42 = probDist.iloc[5:10,1:].values
p46 = probDist.iloc[10:,1:].values

m41 = np.asarray(MulMargProb(4,1)[:4]).T[:4].T.tolist()
m42 = MulMargProb(4,2)[:4]
m46 = MulMargProb(4,6)[:4]


# + init_cell=true
probDist = CondProbDist[4]
#display(probDist)
c5 = probDist.iloc[0,1:].values
p52 = probDist.iloc[1:6,1:].values
p53 = probDist.iloc[6:,1:].values

m52 = MulMargProb(5,2)[:4]
m53 = np.asarray(MulMargProb(5,3)[:4]).T[:3].T.tolist()


# + init_cell=true
probDist = CondProbDist[5]
#display(probDist)
c6 = probDist.iloc[0,1:].values
p61 = probDist.iloc[1:5,1:].values
p62 = probDist.iloc[5:10,1:].values
p63 = probDist.iloc[10:13,1:].values
p64 = probDist.iloc[13:,1:].values

m61 = np.asarray(MulMargProb(6,1)[:]).T[:4].T.tolist()
m62 = MulMargProb(6,2)[:]
m63 = np.asarray(MulMargProb(6,3)[:]).T[:3].T.tolist()
m64 = np.asarray(MulMargProb(6,4)[:]).T[:4].T.tolist()

# -

independence_level = check_independence()
independence_level

# # Part 2 - Constructing Bayesian Network
# We construct Bayseian Network by taking independence values which are high for a better model. 
# ## Storing Marginal and Conditional Probabilities in TabularCPD form for the model
#

# +
from pgmpy.models.BayesianModel import BayesianModel
from pgmpy.factors.discrete import TabularCPD
cpds = {}
cpds['x1'] = TabularCPD('x1', 4 , [MargProb.iloc[:4,1].values.astype(float)/100])
cpds['x2'] = TabularCPD('x2', 5 , [MargProb.iloc[:,2].values.astype(float)/100])
cpds['x3'] = TabularCPD('x3', 3 , [MargProb.iloc[:3,3].values.astype(float)/100])
cpds['x4'] = TabularCPD('x4', 4 , [MargProb.iloc[:4,4].values.astype(float)/100])
cpds['x5'] = TabularCPD('x5', 4 , [MargProb.iloc[:4,5].values.astype(float)/100])
cpds['x6'] = TabularCPD('x6', 5 , [MargProb.iloc[:,6].values.astype(float)/100])

cpds['x12'] = TabularCPD('x2', x_len[2],p12.astype(float)/100, ['x1'], [x_len[1]])
cpds['x14'] = TabularCPD('x4', x_len[4],p14.astype(float)/100, ['x1'], [x_len[1]])
cpds['x16'] = TabularCPD('x6', x_len[6],p16.astype(float)/100, ['x1'], [x_len[1]])

cpds['x23'] = TabularCPD('x3', x_len[3],p23.astype(float)/100, ['x2'], [x_len[2]])
cpds['x25'] = TabularCPD('x5', x_len[5],p25.astype(float)/100, ['x2'], [x_len[2]])

cpds['x32'] = TabularCPD('x2', x_len[2],p32.astype(float)/100, ['x3'], [x_len[3]])
cpds['x35'] = TabularCPD('x5', x_len[5],p35.astype(float)/100, ['x3'], [x_len[3]])
cpds['x36'] = TabularCPD('x6', x_len[6],p36.astype(float)/100, ['x3'], [x_len[3]])

cpds['x41'] = TabularCPD('x1', x_len[1],p41.astype(float)/100, ['x4'], [x_len[4]])
cpds['x42'] = TabularCPD('x2', x_len[2],p42.astype(float)/100, ['x4'], [x_len[4]])
cpds['x46'] = TabularCPD('x6', x_len[6],p46.astype(float)/100, ['x4'], [x_len[4]])

cpds['x52'] = TabularCPD('x2', x_len[2],p52.astype(float)/100, ['x5'], [x_len[5]])
cpds['x53'] = TabularCPD('x3', x_len[3],p53.astype(float)/100, ['x5'], [x_len[5]])

cpds['x61'] = TabularCPD('x1', x_len[1],p61.astype(float)/100, ['x6'], [x_len[6]])
cpds['x62'] = TabularCPD('x2', x_len[2],p62.astype(float)/100, ['x6'], [x_len[6]])
cpds['x63'] = TabularCPD('x3', x_len[3],p63.astype(float)/100, ['x6'], [x_len[6]])
cpds['x64'] = TabularCPD('x4', x_len[4],p64.astype(float)/100, ['x6'], [x_len[6]])

cpds
# -

cpds['x1'].values

# ## Setting up Different Bayesian Models to find a Good Model
# Add the links needed in the BayesianModel() and adding the corresponding CPD values with evidence (if needed) to the model. This allows the model to have the probabilities for each of the nodes considering its corresponding evidences.
#
# We are also converting the Bayesian Model into a Markov Model by 'BayesianModel.to_markov_model(). This makes a Markov Chain Model by utilizing the same nodes and relations betweent the features.

# +
from pgmpy.models.BayesianModel import BayesianModel
from pgmpy.factors.discrete import TabularCPD
from pgmpy.sampling import BayesianModelSampling
bModel = []
k2_B = []
k2_M = []
data = []
marModel = []

bModel.append(BayesianModel([('x1', 'x2'), ('x1', 'x6'), ('x2', 'x3'), ('x6', 'x4'), ('x2', 'x5')]))
bModel[0].add_cpds(cpds['x1'], cpds['x12'], cpds['x16'], cpds['x23'], cpds['x64'], cpds['x25'])
marModel.append(bModel[0].to_markov_model())
#print(gen_k2score(bModel[0]))

bModel.append(BayesianModel([('x1', 'x2'), ('x3', 'x6'), ('x2', 'x3'), ('x6','x4'), ('x2','x5')]))
bModel[1].add_cpds(cpds['x1'], cpds['x12'], cpds['x36'], cpds['x23'], cpds['x64'], cpds['x25'])
marModel.append(bModel[1].to_markov_model())
                                   
bModel.append(BayesianModel([('x5', 'x2'), ('x2', 'x3'), ('x3', 'x6'), ('x6', 'x1'), ('x1', 'x4')]))
bModel[2].add_cpds(cpds['x5'], cpds['x52'], cpds['x23'], cpds['x36'], cpds['x61'], cpds['x14'])
marModel.append(bModel[2].to_markov_model())

bModel.append(BayesianModel([('x4', 'x6'), ('x2', 'x3'), ('x3', 'x5'), ('x6', 'x2'), ('x1', 'x4')]))
bModel[3].add_cpds(cpds['x1'], cpds['x46'], cpds['x23'], cpds['x36'], cpds['x62'], cpds['x14'])
marModel.append(bModel[3].to_markov_model())

bModel.append(BayesianModel([('x5', 'x2'), ('x2', 'x3'), ('x3', 'x6'), ('x6', 'x1'), ('x1', 'x4')]))
bModel[4].add_cpds(cpds['x5'], cpds['x52'], cpds['x23'], cpds['x36'], cpds['x61'], cpds['x14'])
marModel.append(bModel[4].to_markov_model())

bModel.append(BayesianModel([('x6', 'x2'), ('x3', 'x5'), ('x6', 'x1'), ('x2', 'x3'), ('x1', 'x4')]))
bModel[5].add_cpds(cpds['x6'], cpds['x62'], cpds['x35'], cpds['x61'], cpds['x23'], cpds['x14'])
marModel.append(bModel[5].to_markov_model())


# -

# ### Evaluating the K2 Score of the Models
# Here we first need to gather data corresponding to the model passed. We do this by Ancestral Sampling (or Forward Sampling) to replicate the data based on the probabilities given for the model. We gather 1000 data samples for each of the models.
#
# We then find the K2 Score of the Model by using the the data sample generated.

# + init_cell=true
def gen_k2score(model):
    inference = BayesianModelSampling(model)
    data = inference.forward_sample(size=1000)
    k2 = K2Score(data).score(model)
    return k2, data


# -

for i in range(len(bModel)):
    (k2, d) = (gen_k2score(bModel[0]))
    k2_B.append(k2)
    data.append(d)
    #k2_M.append(gen_k2score(marModel[0]))

# ### Finding Best Model
# To find the best model, we consider the model with the highest K2 Score. After finding the best model, we find the most probable 'th' image features for the model.
#
# This is carried out by evaluating the features of the data. We find the highest frequency of a combination of the features in the dataset. This will be the most probable features for the 'th' image. 
# The least probable images would be when the frequency of occurance in the data set is minimum.
#
# The K2 Scores of the models along with the best Model and the corresponding most and least probable image is shown below.
#
#

import networkx as nx
import pylab as plt
def drawModel(model):
    nx.draw(model, with_labels=True)
#plt.show()


# + init_cell=true
import collections
from collections import Counter
def most_prob(inp):
    counter = Counter()
    for a in inp:
        counter[tuple(a)]+=1
    return(counter.most_common(1)[0][0], counter.most_common()[:-4:-1])


# +
for key in range(len(k2_B)):
    print("Model " + str(key) +" Score : " + str(k2_B[key]))
best_model_thDS = [np.argmax(k2_B), max(k2_B)]
print("\nBest Model 'th' Dataset : " +  str(best_model_thDS[1]))
print ("Edges : " + str(bModel[best_model_thDS[0]].edges()))

feat = np.asarray(most_prob(data[best_model_thDS[0]].values.astype(int)))
print ("\nHigh Probable Features : "+str(feat[0]))
print ("Low Probable Features : " + str(feat[1][0][0]) +" , "+ str(feat[1][1][0]))
# -

print ("The best model for Bayesian Network is given as the following :")
drawModel(bModel[best_model_thDS[0]])

print ("The best model for Markov Model is given as the following :")
drawModel(marModel[best_model_thDS[0]])


# # Part 3 : Markov Models
#
# In the previous step, we have converted the models in to Markov Models. Using both Bayesian and Markov Models, we try to predict (or inference) the value of a feature assuming we know the rest of the feature values. This is accomplished by using Belief Propagation. 
#
# Belief Propagation (sum-product) is a message passing algorithm used for performing inferences of graphical models like Bayesian and Markov Models. To find the value of the missing node, it calculates the marginal distribution for each of the unknown nodes and conditional on the observed nodes. This enables us to determine the value for the unknown node.
#
# Below, the unknown variables are in 'variables' while the known variables and their values are stored in a dictionary as the 'evidence'. After querying for result, we get the predicted value of the 'variables'. This is executed for both the Bayesian and Markov Models for 'th' and 'and' Datasets.

# +

def cal_accuracy_th(data):
    m_result = []
    data = data[cols]
    for d in data.values:
        m_result.append(belief_propagation.map_query(variables=['x3'],evidence={'x1': d[0], 'x2':d[1],  'x4': d[3], 'x5': d[4], 'x6' : d[5]})['x3'])
    label = data['x3'].values
    acc = (m_result==label).sum()/len(label)*100
    return acc


# -

def cal_accuracy_and(data):
    m_result = []
    data = data[cols]
    for d in data.values:
        m_result.append(belief_propagation.map_query(variables=['f9'],evidence={'f1': d[0], 'f2':d[1], 'f3' : d[2], 'f4': d[3], 'f5': d[4], 'f6': d[5], 'f7':d[6], 'f8': d[7]})['f9'])
    label = data['f9'].values
    acc = (m_result==label).sum()/len(label)*100
    return acc


# +
# %%time
cols = ['x1', 'x2', 'x3', 'x4', 'x5', 'x6']

belief_propagation = BeliefPropagation(bModel[best_model_thDS[0]])
belief_propagation.calibrate()
acc_BN_th = cal_accuracy_th(data[best_model_thDS[0]])
belief_propagation = BeliefPropagation(marModel[best_model_thDS[0]])
belief_propagation.calibrate()
acc_MM_th = cal_accuracy_th(data[best_model_thDS[0]])
print ("Predicting node x3")
print ("Accuracy Bayesian Model : " + str(acc_BN_th))
print ("Accuracy Markov Model : " + str(acc_MM_th))
# -

# # Part 4 : Bayesian Network for 'and' Image Dataset
# Here we construct the Bayesian Network for the 'and' dataset. We have the image data with 9 features each. We first find the correlations and independencies of the features by using test_conditional_independence(). Only if the significance level is higher than a threshold, we consider those conditionals as true. 

# + init_cell=true
andDataset = pd.read_csv('AND-Features(1).csv')
from pgmpy.estimators import ConstraintBasedEstimator
est = ConstraintBasedEstimator(andDataset)
def is_independent(X, Y, Zs=[], significance_level=0.30):
    return est.test_conditional_independence(X, Y, Zs)[1] >= significance_level


# -

cols = andDataset.columns.values.tolist()[2:]
andIndependent = []
for i in range(len(cols)):
    for j in range(len(cols)):
        val = is_independent(cols[i], cols[j])
        #print(str(i+1) + " and " + str(j+1) +" : " + str(val))
        if val:
            andIndependent.append([cols[i], cols[j]])  
print ("Related Features of 'and' Dataset : \n" + str(andIndependent))
print ("No. of dependencies : " + str(len(andIndependent)))

# ## Best Models
# Here we use HillClimbSearch to estimate the best model for the given dataset. It performs a Heuristic Search which selects the best route out of possible routes.
#
# We also build various other models based on the independency levels shown previously. We estimate the CPDs for the different models using Maximum Likelihood Estimator and then find the corresponding K2 Scores. 

# +
from pgmpy.estimators import HillClimbSearch

hc = HillClimbSearch(andDataset.iloc[:,2:], scoring_method=K2Score(andDataset))
best_model_andDS = hc.estimate()

# + code_folding=[]
from pgmpy.estimators import BayesianEstimator
from pgmpy.estimators import MaximumLikelihoodEstimator
score = []
m = []
m.append([('f1', 'f3'), ('f2', 'f3'), ('f2', 'f8'), ('f3', 'f9'),('f4', 'f7'),('f8', 'f3'), ('f6', 'f1'), ('f5','f2')])
m.append([('f1', 'f3'), ('f2', 'f4'), ('f4', 'f5'), ('f2', 'f8'), ('f4','f7'), ('f9','f3'),('f3', 'f7'), ('f6','f1')])
m.append([('f1', 'f5'), ('f5', 'f2'), ('f2', 'f7'), ('f7', 'f3'),('f8', 'f2'),('f7','f4'), ('f6','f1'), ('f9','f3')])
m.append([('f1', 'f5'), ('f3', 'f2'), ('f8', 'f2'), ('f6', 'f1'),('f7', 'f2'), ('f7','f3'),('f3','f9'), ('f4','f7')])


for mod in m:
    model = BayesianModel(mod)
    estimator = MaximumLikelihoodEstimator(model, andDataset.iloc[:,2:])
    and_cpd = estimator.get_parameters()
    for cpd in and_cpd:
        model.add_cpds(cpd)
    score.append(K2Score(andDataset.iloc[:,2:]).score(model))

for i in range(len(score)):
    print ("Model "+str(i)+" : " + str(score[i]))
print ("\nBest Model 'and' Dataset " +  str(K2Score(andDataset.iloc[:,2:]).score(best_model_andDS)))
print("Best Model Edges : " + str(best_model_andDS.edges()))

# -

best_model_andDS.get_cpds()
estimator = MaximumLikelihoodEstimator(best_model_andDS, andDataset.iloc[:,2:])
and_cpd = estimator.get_parameters()
for cpd in and_cpd:
    best_model_andDS.add_cpds(cpd)
best_model_andDS.get_cpds()

# +
# %%time
cols = ['f1', 'f2', 'f3', 'f4', 'f5', 'f6', 'f7', 'f8', 'f9']

MM_bModel_andDS = best_model_andDS.to_markov_model()

belief_propagation = BeliefPropagation(best_model_andDS)
acc_BN_and = cal_accuracy_and(andDataset.iloc[:,2:])

belief_propagation = BeliefPropagation(MM_bModel_andDS)
acc_MM_and = cal_accuracy_and(andDataset.iloc[:,2:])

print ("'and' Dataset Inferences :")
print ("Predicting node f9 :")
print ("Bayesian Model Accuracy : " + str(acc_BN_and))
print ("Markov Model Accuracy   : " + str(acc_MM_and))

# -

print ("BN in 'and' Dataset")
drawModel(best_model_andDS)

# +

drawModel(MM_bModel_andDS)
# -



StarcoderdataJupyterScriptsDedupFiltered