diff --git a/examples/decisionExample.py b/examples/decisionExample.py deleted file mode 100644 index a86c34e..0000000 --- a/examples/decisionExample.py +++ /dev/null @@ -1,178 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Thu Dec 7 17:59:32 2017 - -@author: jpoeppel -""" - -from primo2.networks import DecisionNetwork -from primo2.nodes import DiscreteNode, DecisionNode, UtilityNode - -from primo2.inference.decision import VariableElimination - -""" -PHD example 7.3 from Bayesian Reasoning and Machine Learning - Barber -""" -net = DecisionNetwork() - -education = DecisionNode("education", decisions=["do Phd", "no Phd"]) #E - -income = DiscreteNode("income", values=["low", "average", "high"]) #I -nobel = DiscreteNode("nobel", values=["prize", "no prize"]) #P - -costs = UtilityNode("costs") #UC -gains = UtilityNode("gains") #UB - - -#Add nodes to network. They can be treated the same -net.add_node(education) -net.add_node(income) -net.add_node(nobel) - -net.add_node(costs) -net.add_node(gains) - - -#Add edges. Edges can either be acutal dependencies or information links. -#The type is figured out by the nodes themsevles -net.add_edge(education, costs) -net.add_edge(education, nobel) -net.add_edge(education, income) - -net.add_edge(nobel, income) -net.add_edge(income, gains) - -#Define CPTs: (Needs to be done AFTER the structure is defined as that) -#determines the table structure for the different nodes - -income.set_probability("low", 0.1, parentValues={"education":"do Phd", "nobel":"no prize"}) -income.set_probability("low", 0.2, parentValues={"education":"no Phd", "nobel":"no prize"}) -income.set_probability("low", 0.01, parentValues={"education":"do Phd", "nobel":"prize"}) -income.set_probability("low", 0.01, parentValues={"education":"no Phd", "nobel":"prize"}) - -income.set_probability("average", 0.5, parentValues={"education":"do Phd", "nobel":"no prize"}) -income.set_probability("average", 0.6, parentValues={"education":"no Phd", "nobel":"no prize"}) -income.set_probability("average", 0.04, parentValues={"education":"do Phd", "nobel":"prize"}) -income.set_probability("average", 0.04, parentValues={"education":"no Phd", "nobel":"prize"}) - -income.set_probability("high", 0.4, parentValues={"education":"do Phd", "nobel":"no prize"}) -income.set_probability("high", 0.2, parentValues={"education":"no Phd", "nobel":"no prize"}) -income.set_probability("high", 0.95, parentValues={"education":"do Phd", "nobel":"prize"}) -income.set_probability("high", 0.95, parentValues={"education":"no Phd", "nobel":"prize"}) - - -nobel.set_probability("prize", 0.0000001, parentValues={"education":"no Phd"}) -nobel.set_probability("prize", 0.001, parentValues={"education":"do Phd"}) - -nobel.set_probability("no prize", 0.9999999, parentValues={"education":"no Phd"}) -nobel.set_probability("no prize", 0.999, parentValues={"education":"do Phd"}) - - -#Define utilities - -costs.set_utility(-50000, parentValues={"education":"do Phd"}) -costs.set_utility(0, parentValues={"education":"no Phd"}) - -gains.set_utility(100000, parentValues={"income":"low"}) -gains.set_utility(200000, parentValues={"income":"average"}) -gains.set_utility(500000, parentValues={"income":"high"}) - - -ve = VariableElimination(net) - -print("Expected Utility for doing a Phd: {}".format(ve.expected_utility(decisions={"education":"do Phd"}))) -print("Expected Utility for not doing a Phd: {}".format(ve.expected_utility(decisions={"education":"no Phd"}))) - -print("Optimal deciosn: ", ve.get_optimal_decisions(["education"])) - - - - -""" -PHD + Startup example 7.4 from Bayesian Reasoning and Machine Learning - Barber -""" - -net = DecisionNetwork() - -education = DecisionNode("education", decisions=["do Phd", "no Phd"]) #E -startup = DecisionNode("startup", decisions=["start up", "no start up"]) # S - -income = DiscreteNode("income", values=["low", "average", "high"]) #I -nobel = DiscreteNode("nobel", values=["prize", "no prize"]) #P - -costsEducation = UtilityNode("costsE") #UC -costsStartUp = UtilityNode("costsS") #US -gains = UtilityNode("gains") #UB - - - -#Add nodes to network. They can be treated the same -net.add_node(education) -net.add_node(startup) -net.add_node(income) -net.add_node(nobel) - -net.add_node(costsEducation) -net.add_node(costsStartUp) -net.add_node(gains) - - -#Add edges. Edges can either be acutal dependencies or information links. -#The type is figured out by the nodes themsevles -net.add_edge(education, costsEducation) -net.add_edge(education, nobel) - -net.add_edge(startup, income) -net.add_edge(startup, costsStartUp) - -net.add_edge(nobel, income) -net.add_edge(income, gains) - -#Define CPTs: (Needs to be done AFTER the structure is defined as that) -#determines the table structure for the different nodes - -income.set_probability("low", 0.1, parentValues={"startup":"start up", "nobel":"no prize"}) -income.set_probability("low", 0.2, parentValues={"startup":"no start up", "nobel":"no prize"}) -income.set_probability("low", 0.005, parentValues={"startup":"start up", "nobel":"prize"}) -income.set_probability("low", 0.05, parentValues={"startup":"no start up", "nobel":"prize"}) - -income.set_probability("average", 0.5, parentValues={"startup":"start up", "nobel":"no prize"}) -income.set_probability("average", 0.6, parentValues={"startup":"no start up", "nobel":"no prize"}) -income.set_probability("average", 0.005, parentValues={"startup":"start up", "nobel":"prize"}) -income.set_probability("average", 0.15, parentValues={"startup":"no start up", "nobel":"prize"}) - -income.set_probability("high", 0.4, parentValues={"startup":"start up", "nobel":"no prize"}) -income.set_probability("high", 0.2, parentValues={"startup":"no start up", "nobel":"no prize"}) -income.set_probability("high", 0.99, parentValues={"startup":"start up", "nobel":"prize"}) -income.set_probability("high", 0.8, parentValues={"startup":"no start up", "nobel":"prize"}) - - -nobel.set_probability("prize", 0.0000001, parentValues={"education":"no Phd"}) -nobel.set_probability("prize", 0.001, parentValues={"education":"do Phd"}) - -nobel.set_probability("no prize", 0.9999999, parentValues={"education":"no Phd"}) -nobel.set_probability("no prize", 0.999, parentValues={"education":"do Phd"}) - - -#Define utilities - -costsEducation.set_utility(-50000, parentValues={"education":"do Phd"}) -costsEducation.set_utility(0, parentValues={"education":"no Phd"}) - -costsStartUp.set_utility(-200000, parentValues={"startup":"start up"}) -costsStartUp.set_utility(0, parentValues={"startup":"no start up"}) - -gains.set_utility(100000, parentValues={"income":"low"}) -gains.set_utility(200000, parentValues={"income":"average"}) -gains.set_utility(500000, parentValues={"income":"high"}) - -ve = VariableElimination(net) - -print("Expected Utility for doing a Phd + startup: {}".format(ve.expected_utility(decisions={"education":"do Phd", "startup": "start up"}))) -print("Expected Utility for doing a Phd + no startup: {}".format(ve.expected_utility(decisions={"education":"do Phd", "startup": "no start up"}))) -print("Expected Utility for not doing a Phd + startup: {}".format(ve.expected_utility(decisions={"education":"no Phd", "startup": "start up"}))) -print("Expected Utility for not doing a Phd + no startup: {}".format(ve.expected_utility(decisions={"education":"no Phd", "startup": "no start up"}))) - - -print("Optimal deciosn: ", ve.get_optimal_decisions(["startup", "education"])) \ No newline at end of file diff --git a/examples/decision_example.py b/examples/decision_example.py new file mode 100644 index 0000000..9fda6d7 --- /dev/null +++ b/examples/decision_example.py @@ -0,0 +1,410 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Thu Dec 7 17:59:32 2017 + +@author: jpoeppel +""" + +from primo2.networks import DecisionNetwork, DynamicDecisionNetwork, d0_net, Two_TDN +from primo2.nodes import DiscreteNode, DecisionNode, UtilityNode + +from primo2.inference.decision import VariableElimination +from math import exp, sqrt, pi +import matplotlib.pyplot as plt +import numpy as np +import matplotlib.mlab as mlab +import pandas as pd + +print ("") +print ("PHD example 7.3 from Bayesian Reasoning and Machine Learning - Barber") +print ("") + +net = DecisionNetwork() + +education = DecisionNode("education", decisions=["do Phd", "no Phd"]) # E + +income = DiscreteNode("income", values=["low", "average", "high"]) # I +nobel = DiscreteNode("nobel", values=["prize", "no prize"]) # P + +costs = UtilityNode("costs") # UC +gains = UtilityNode("gains") # UB + +# Add nodes to network. They can be treated the same +net.add_node(education) +net.add_node(income) +net.add_node(nobel) + +net.add_node(costs) +net.add_node(gains) + +# Add edges. Edges can either be acutal dependencies or information links. +# The type is figured out by the nodes themsevles +net.add_edge(education, costs) +net.add_edge(education, nobel) +net.add_edge(education, income) + +net.add_edge(nobel, income) +net.add_edge(income, gains) + +# Define CPTs: (Needs to be done AFTER the structure is defined as that) +# determines the table structure for the different nodes + +income.set_probability("low", 0.1, parentValues={"education": "do Phd", "nobel": "no prize"}) +income.set_probability("low", 0.2, parentValues={"education": "no Phd", "nobel": "no prize"}) +income.set_probability("low", 0.01, parentValues={"education": "do Phd", "nobel": "prize"}) +income.set_probability("low", 0.01, parentValues={"education": "no Phd", "nobel": "prize"}) + +income.set_probability("average", 0.5, parentValues={"education": "do Phd", "nobel": "no prize"}) +income.set_probability("average", 0.6, parentValues={"education": "no Phd", "nobel": "no prize"}) +income.set_probability("average", 0.04, parentValues={"education": "do Phd", "nobel": "prize"}) +income.set_probability("average", 0.04, parentValues={"education": "no Phd", "nobel": "prize"}) + +income.set_probability("high", 0.4, parentValues={"education": "do Phd", "nobel": "no prize"}) +income.set_probability("high", 0.2, parentValues={"education": "no Phd", "nobel": "no prize"}) +income.set_probability("high", 0.95, parentValues={"education": "do Phd", "nobel": "prize"}) +income.set_probability("high", 0.95, parentValues={"education": "no Phd", "nobel": "prize"}) + +nobel.set_probability("prize", 0.0000001, parentValues={"education": "no Phd"}) +nobel.set_probability("prize", 0.001, parentValues={"education": "do Phd"}) + +nobel.set_probability("no prize", 0.9999999, parentValues={"education": "no Phd"}) +nobel.set_probability("no prize", 0.999, parentValues={"education": "do Phd"}) + +# Define utilities + +costs.set_utility(-50000, parentValues={"education": "do Phd"}) +costs.set_utility(0, parentValues={"education": "no Phd"}) + +gains.set_utility(100000, parentValues={"income": "low"}) +gains.set_utility(200000, parentValues={"income": "average"}) +gains.set_utility(500000, parentValues={"income": "high"}) + +net.set_partial_ordering([education, [income, nobel]]) +ve = VariableElimination(net) + +print("Expected Utility for doing a Phd: {}".format(ve.expected_utility(decisions={"education": "do Phd"}))) +print("Expected Utility for not doing a Phd: {}".format(ve.expected_utility(decisions={"education": "no Phd"}))) +print "Optimal decision using max_sum: ", ve.max_sum("education") +print "Get optimal decision using classic algorithm: ", ve.get_optimal_decisions(["education"]) + +print ("") +print ("PHD example 7.4 from Bayesian Reasoning and Machine Learning - Barber") +print ("") + +net = DecisionNetwork() + +education = DecisionNode("education", decisions=["do Phd", "no Phd"]) # E +startup = DecisionNode("startup", decisions=["start up", "no start up"]) # S + +income = DiscreteNode("income", values=["low", "average", "high"]) # I +nobel = DiscreteNode("nobel", values=["prize", "no prize"]) # P + +costsEducation = UtilityNode("costsE") # UC +costsStartUp = UtilityNode("costsS") # US +gains = UtilityNode("gains") # UB + +# Add nodes to network. They can be treated the same +net.add_node(education) +net.add_node(startup) +net.add_node(income) +net.add_node(nobel) + +net.add_node(costsEducation) +net.add_node(costsStartUp) +net.add_node(gains) + +# Add edges. Edges can either be acutal dependencies or information links. +# The type is figured out by the nodes themsevles +net.add_edge(education, costsEducation) +net.add_edge(education, nobel) + +net.add_edge(startup, income) +net.add_edge(startup, costsStartUp) + +net.add_edge(nobel, income) +net.add_edge(income, gains) + +# Define CPTs: (Needs to be done AFTER the structure is defined as that) +# determines the table structure for the different nodes + +income.set_probability("low", 0.1, parentValues={"startup": "start up", "nobel": "no prize"}) +income.set_probability("low", 0.2, parentValues={"startup": "no start up", "nobel": "no prize"}) +income.set_probability("low", 0.005, parentValues={"startup": "start up", "nobel": "prize"}) +income.set_probability("low", 0.05, parentValues={"startup": "no start up", "nobel": "prize"}) + +income.set_probability("average", 0.5, parentValues={"startup": "start up", "nobel": "no prize"}) +income.set_probability("average", 0.6, parentValues={"startup": "no start up", "nobel": "no prize"}) +income.set_probability("average", 0.005, parentValues={"startup": "start up", "nobel": "prize"}) +income.set_probability("average", 0.15, parentValues={"startup": "no start up", "nobel": "prize"}) + +income.set_probability("high", 0.4, parentValues={"startup": "start up", "nobel": "no prize"}) +income.set_probability("high", 0.2, parentValues={"startup": "no start up", "nobel": "no prize"}) +income.set_probability("high", 0.99, parentValues={"startup": "start up", "nobel": "prize"}) +income.set_probability("high", 0.8, parentValues={"startup": "no start up", "nobel": "prize"}) + +nobel.set_probability("prize", 0.0000001, parentValues={"education": "no Phd"}) +nobel.set_probability("prize", 0.001, parentValues={"education": "do Phd"}) + +nobel.set_probability("no prize", 0.9999999, parentValues={"education": "no Phd"}) +nobel.set_probability("no prize", 0.999, parentValues={"education": "do Phd"}) + +# Define utilities + +costsEducation.set_utility(-50000, parentValues={"education": "do Phd"}) +costsEducation.set_utility(0, parentValues={"education": "no Phd"}) + +costsStartUp.set_utility(-200000, parentValues={"startup": "start up"}) +costsStartUp.set_utility(0, parentValues={"startup": "no start up"}) + +gains.set_utility(100000, parentValues={"income": "low"}) +gains.set_utility(200000, parentValues={"income": "average"}) +gains.set_utility(500000, parentValues={"income": "high"}) + +net.set_partial_ordering([education, nobel, startup, income]) +ve = VariableElimination(net) + +print("Expected Utility for doing a Phd + startup: {}".format( + ve.expected_utility(decisions={"education": "do Phd", "startup": "start up"}))) +print("Expected Utility for doing a Phd + no startup: {}".format( + ve.expected_utility(decisions={"education": "do Phd", "startup": "no start up"}))) +print("Expected Utility for not doing a Phd + startup: {}".format( + ve.expected_utility(decisions={"education": "no Phd", "startup": "start up"}))) +print("Expected Utility for not doing a Phd + no startup: {}".format( + ve.expected_utility(decisions={"education": "no Phd", "startup": "no start up"}))) +print "Get optimal decision for education using max_sum Algorithm: {}".format(ve.max_sum("education")) +print "Get optimal decision for startup using max_sum Algorithm: {}".format(ve.max_sum("startup")) + +print "" +print "Fever Example 4.3.5 from Bayesian Artificial Intelligence - Korb Nicholson" +print "" + +d0 = d0_net() + +flu_0 = DiscreteNode("flu_0", values=["True", "False"]) +fever_0 = DiscreteNode("fever_0", values=["True", "False"]) +therm_0 = DiscreteNode("therm_0", values=["True", "False"]) +take_aspirin_0 = DecisionNode("take_aspirin_0", decisions=["Yes", "No"]) +fever_later_0 = DiscreteNode("fever_later_0", values=["True", "False"]) +reaction_0 = DiscreteNode("reaction_0", values=["Yes", "No"]) +utility_0 = UtilityNode("utility_0") + +d0.add_nodes([flu_0, fever_0, therm_0, take_aspirin_0, fever_later_0, reaction_0, utility_0]) + +d0.add_edge(flu_0, fever_0) +d0.add_edge(fever_0, therm_0) +d0.add_edge(fever_0, fever_later_0) +d0.add_edge(fever_later_0, utility_0) +d0.add_edge(take_aspirin_0, reaction_0) +d0.add_edge(take_aspirin_0, fever_later_0) +d0.add_edge(reaction_0, utility_0) + +flu_0.set_probability("True", 0.05) +flu_0.set_probability("False", 0.95) + +fever_0.set_probability("True", 0.95, parentValues={"flu_0": "True"}) +fever_0.set_probability("True", 0.02, parentValues={"flu_0": "False"}) +fever_0.set_probability("False", 0.05, parentValues={"flu_0": "True"}) +fever_0.set_probability("False", 0.98, parentValues={"flu_0": "False"}) + +therm_0.set_probability("True", 0.90, parentValues={"fever_0": "True"}) +therm_0.set_probability("True", 0.05, parentValues={"fever_0": "False"}) +therm_0.set_probability("False", 0.10, parentValues={"fever_0": "True"}) +therm_0.set_probability("False", 0.95, parentValues={"fever_0": "False"}) + +fever_later_0.set_probability("True", 0.05, parentValues={"fever_0": "True", "take_aspirin_0": "Yes"}) +fever_later_0.set_probability("True", 0.90, parentValues={"fever_0": "True", "take_aspirin_0": "No"}) +fever_later_0.set_probability("True", 0.01, parentValues={"fever_0": "False", "take_aspirin_0": "Yes"}) +fever_later_0.set_probability("True", 0.02, parentValues={"fever_0": "False", "take_aspirin_0": "No"}) +fever_later_0.set_probability("False", 0.95, parentValues={"fever_0": "True", "take_aspirin_0": "Yes"}) +fever_later_0.set_probability("False", 0.10, parentValues={"fever_0": "True", "take_aspirin_0": "No"}) +fever_later_0.set_probability("False", 0.99, parentValues={"fever_0": "False", "take_aspirin_0": "Yes"}) +fever_later_0.set_probability("False", 0.98, parentValues={"fever_0": "False", "take_aspirin_0": "No"}) + +reaction_0.set_probability("Yes", 0.05, parentValues={"take_aspirin_0": "Yes"}) +reaction_0.set_probability("Yes", 0.00, parentValues={"take_aspirin_0": "No"}) +reaction_0.set_probability("No", 0.95, parentValues={"take_aspirin_0": "Yes"}) +reaction_0.set_probability("No", 1.0, parentValues={"take_aspirin_0": "No"}) + +utility_0.set_utility(-50, {"fever_later_0": "True", "reaction_0": "Yes"}) +utility_0.set_utility(-10, {"fever_later_0": "True", "reaction_0": "No"}) +utility_0.set_utility(-30, {"fever_later_0": "False", "reaction_0": "Yes"}) +utility_0.set_utility(50, {"fever_later_0": "False", "reaction_0": "No"}) + +d0.set_partial_ordering([["flu_0", "fever_0", "therm_0"], "take_aspirin_0", ["fever_later_0", "reaction_0"]]) + +ve = VariableElimination(d0) + +print("Expected Utility for Taking Aspirin: {}".format(ve.expected_utility(decisions={"take_aspirin_0": "Yes"}))) +print("Expected Utility for not Taking Aspirin: {}".format(ve.expected_utility(decisions={"take_aspirin_0": "No"}))) +print "Get optimal decision using classic algorithm: ", ve.get_optimal_decisions(["take_aspirin_0"]) +print "Get optimal decision for education using max_sum Algorithm: {}".format(ve.max_sum("take_aspirin_0")) + +print "" +print "L2Tor example" +print "" + +d0 = d0_net() + +skill_0 = DiscreteNode("skill_0", values=[0, 1, 2, 3, 4, 5]) +action_0 = DecisionNode("action_0", decisions=["1", "2", "3", "4"]) +observation_0 = DiscreteNode("observation_0", values=["+O", "-O"]) + +d0.add_nodes([skill_0, action_0, observation_0]) +d0.add_edge(skill_0, observation_0) +d0.add_edge(skill_0, action_0) +d0.add_edge(action_0, observation_0) + +d0.set_partial_ordering(["skill_0", "action_0", "observation_0"]) + +skill_0.set_cpd([1.0 / 6, 1.0 / 6, 1.0 / 6, 1.0 / 6, 1.0 / 6, 1.0 / 6]) + +action_0.set_cpd([[0.40, 0.30, 0.20, 0.15, 0.10, 0.05], + [0.30, 0.35, 0.40, 0.25, 0.25, 0.25], + [0.25, 0.25, 0.25, 0.40, 0.35, 0.30], + [0.05, 0.10, 0.15, 0.20, 0.30, 0.40]]) + +action_0.set_action_norm({'1': {'loc': .15, 'scale': .05}, + '2': {'loc': .40, 'scale': .05}, + '3': {'loc': .625, 'scale': .05}, + '4': {'loc': .85, 'scale': .05}}) + +observation_0.set_cpd([[[0.50, 0.3, 0.25, 0.15], + [0.55, 0.40, 0.30, 0.20], + [0.65, 0.55, 0.40, 0.30], + [0.75, 0.65, 0.50, 0.40], + [0.85, 0.75, 0.60, 0.50], + [0.95, 0.85, 0.70, 0.60]], + [[0.50, 0.67, 0.75, 0.85], + [0.45, 0.60, 0.70, 0.80], + [0.35, 0.45, 0.60, 0.70], + [0.25, 0.35, 0.50, 0.60], + [0.15, 0.25, 0.40, 0.50], + [0.05, 0.15, 0.30, 0.40]]]) + +two_tdn = Two_TDN() + +skill_t = DiscreteNode("skill_t", values=[0, 1, 2, 3, 4, 5]) + +action_t = DecisionNode("action_t", decisions=["1", "2", "3", "4"]) +observation_t = DiscreteNode("observation_t", values=["+O", "-O"]) + +skill_t_plus_one = DiscreteNode("skill_t_plus_one", values=[0, 1, 2, 3, 4, 5]) + +action_t_plus_one = DecisionNode("action_t_plus_one", decisions=["1", "2", "3", "4"]) +observation_t_plus_one = DiscreteNode("observation_t_plus_one", values=["+O", "-O"]) + +transition_dict = {'+O': {"1": {0: {0: 0.60, 1: 0.00, 2: 0.00, 3: 0.00, 4: 0.00, 5: 0.00}, + 1: {0: 0.25, 1: 0.70, 2: 0.00, 3: 0.00, 4: 0.00, 5: 0.00}, + 2: {0: 0.15, 1: 0.20, 2: 0.80, 3: 0.00, 4: 0.00, 5: 0.00}, + 3: {0: 0.00, 1: 0.10, 2: 0.15, 3: 0.85, 4: 0.00, 5: 0.00}, + 4: {0: 0.00, 1: 0.00, 2: 0.05, 3: 0.10, 4: 0.90, 5: 0.00}, + 5: {0: 0.00, 1: 0.00, 2: 0.00, 3: 0.05, 4: 0.10, 5: 1.00} + }, + "2": {0: {0: 0.55, 1: 0.00, 2: 0.00, 3: 0.00, 4: 0.00, 5: 0.00}, + 1: {0: 0.35, 1: 0.60, 2: 0.00, 3: 0.00, 4: 0.00, 5: 0.00}, + 2: {0: 0.10, 1: 0.30, 2: 0.70, 3: 0.00, 4: 0.00, 5: 0.00}, + 3: {0: 0.00, 1: 0.10, 2: 0.20, 3: 0.70, 4: 0.00, 5: 0.00}, + 4: {0: 0.00, 1: 0.00, 2: 0.10, 3: 0.25, 4: 0.75, 5: 0.00}, + 5: {0: 0.00, 1: 0.00, 2: 0.00, 3: 0.05, 4: 0.25, 5: 1.00} + }, + "3": {0: {0: 0.40, 1: 0.00, 2: 0.00, 3: 0.00, 4: 0.00, 5: 0.00}, + 1: {0: 0.40, 1: 0.40, 2: 0.00, 3: 0.00, 4: 0.00, 5: 0.00}, + 2: {0: 0.20, 1: 0.40, 2: 0.40, 3: 0.00, 4: 0.00, 5: 0.00}, + 3: {0: 0.00, 1: 0.20, 2: 0.40, 3: 0.50, 4: 0.00, 5: 0.00}, + 4: {0: 0.00, 1: 0.00, 2: 0.20, 3: 0.40, 4: 0.60, 5: 0.00}, + 5: {0: 0.00, 1: 0.00, 2: 0.00, 3: 0.10, 4: 0.40, 5: 1.00} + }, + "4": {0: {0: 0.30, 1: 0.00, 2: 0.00, 3: 0.00, 4: 0.00, 5: 0.00}, + 1: {0: 0.45, 1: 0.30, 2: 0.00, 3: 0.00, 4: 0.00, 5: 0.00}, + 2: {0: 0.20, 1: 0.45, 2: 0.30, 3: 0.00, 4: 0.00, 5: 0.00}, + 3: {0: 0.05, 1: 0.20, 2: 0.45, 3: 0.30, 4: 0.00, 5: 0.00}, + 4: {0: 0.00, 1: 0.05, 2: 0.20, 3: 0.45, 4: 0.40, 5: 0.00}, + 5: {0: 0.00, 1: 0.00, 2: 0.05, 3: 0.25, 4: 0.60, 5: 1.00} + } + }, + '-O': {"1": {0: {0: 1.00, 1: 0.10, 2: 0.05, 3: 0.00, 4: 0.00, 5: 0.00}, + 1: {0: 0.00, 1: 0.90, 2: 0.10, 3: 0.15, 4: 0.00, 5: 0.00}, + 2: {0: 0.00, 1: 0.00, 2: 0.85, 3: 0.25, 4: 0.20, 5: 0.00}, + 3: {0: 0.00, 1: 0.00, 2: 0.00, 3: 0.60, 4: 0.30, 5: 0.25}, + 4: {0: 0.00, 1: 0.00, 2: 0.00, 3: 0.00, 4: 0.50, 5: 0.35}, + 5: {0: 0.00, 1: 0.00, 2: 0.00, 3: 0.00, 4: 0.00, 5: 0.40} + }, + "2": {0: {0: 1.00, 1: 0.25, 2: 0.05, 3: 0.00, 4: 0.00, 5: 0.00}, + 1: {0: 0.00, 1: 0.75, 2: 0.25, 3: 0.10, 4: 0.00, 5: 0.00}, + 2: {0: 0.00, 1: 0.00, 2: 0.70, 3: 0.25, 4: 0.10, 5: 0.00}, + 3: {0: 0.00, 1: 0.00, 2: 0.00, 3: 0.65, 4: 0.30, 5: 0.10}, + 4: {0: 0.00, 1: 0.00, 2: 0.00, 3: 0.00, 4: 0.60, 5: 0.35}, + 5: {0: 0.00, 1: 0.00, 2: 0.00, 3: 0.00, 4: 0.00, 5: 0.55} + }, + "3": {0: {0: 1.00, 1: 0.20, 2: 0.20, 3: 0.00, 4: 0.00, 5: 0.00}, + 1: {0: 0.00, 1: 0.80, 2: 0.40, 3: 0.20, 4: 0.00, 5: 0.00}, + 2: {0: 0.00, 1: 0.00, 2: 0.40, 3: 0.40, 4: 0.20, 5: 0.00}, + 3: {0: 0.00, 1: 0.00, 2: 0.00, 3: 0.40, 4: 0.40, 5: 0.20}, + 4: {0: 0.00, 1: 0.00, 2: 0.00, 3: 0.00, 4: 0.40, 5: 0.40}, + 5: {0: 0.00, 1: 0.00, 2: 0.00, 3: 0.00, 4: 0.00, 5: 0.40} + }, + "4": {0: {0: 1.00, 1: 0.40, 2: 0.15, 3: 0.00, 4: 0.00, 5: 0.00}, + 1: {0: 0.00, 1: 0.60, 2: 0.35, 3: 0.15, 4: 0.00, 5: 0.00}, + 2: {0: 0.00, 1: 0.00, 2: 0.50, 3: 0.35, 4: 0.15, 5: 0.00}, + 3: {0: 0.00, 1: 0.00, 2: 0.00, 3: 0.50, 4: 0.35, 5: 0.15}, + 4: {0: 0.00, 1: 0.00, 2: 0.00, 3: 0.00, 4: 0.50, 5: 0.35}, + 5: {0: 0.00, 1: 0.00, 2: 0.00, 3: 0.00, 4: 0.00, 5: 0.50} + } + } + } + +two_tdn.add_nodes([skill_t, action_t, observation_t, skill_t_plus_one, action_t_plus_one, observation_t_plus_one]) + +two_tdn.add_intra_edge(skill_t, action_t) +two_tdn.add_intra_edge(skill_t, observation_t) +two_tdn.add_intra_edge(action_t, observation_t) + +two_tdn.add_inter_edge(skill_t, skill_t_plus_one) +two_tdn.add_inter_edge(action_t, skill_t_plus_one) +two_tdn.add_inter_edge(observation_t, skill_t_plus_one) + + +# This is just a helper function to translate the given dictionary to an array (for this specific example. +# Should not be applied to different examples. + +def convert_to_array(transition_): + if isinstance(transition_, dict): + cpd = [] + + def get_table(transition, answer, action): + return pd.DataFrame(transition[answer][action]) + + for i in transition_.keys(): + temp = [] + for j in transition_[i]: + temp.append(get_table(transition_, i, j).as_matrix()) + + cpd.append(temp) + + # return np.transpose(np.array(cpd).T, (0, 1, 3, 2)) + return np.array(cpd).T + + +print(convert_to_array(transition_dict).shape) +two_tdn.add_transition(skill_t_plus_one, convert_to_array(transition_dict)) + +DDN = DynamicDecisionNetwork(d0, two_tdn) +new_net = DDN.unroll(2) + +for i in new_net.get_all_nodes(): + print "Current One:" + print i.name + print(i.cpd) + +# ve = VariableElimination(new_net) +# print("Optimal decision: ", ve.get_optimal_decisions(["action_0"])) +# print("Get optimal decision for action_0 using max_sum Algorithm: {}".format(ve.max_sum("action_0"))) + +# mu = [.15, .40, .625, .85] +# variance = [.05, .05, .05, .05] +# sigma = [sqrt(i) for i in variance] +# x = [np.linspace(mu[i] - 3 * sigma[i], mu[i] + 3 * sigma[i], 100) for i, _ in enumerate(mu)] +# for i, v in enumerate(x): +# plt.plot(v, mlab.normpdf(v, mu[i], sigma[i])) +# plt.show() diff --git a/primo2/inference/decision.py b/primo2/inference/decision.py index 5a5d5fb..c821c7f 100644 --- a/primo2/inference/decision.py +++ b/primo2/inference/decision.py @@ -6,25 +6,27 @@ @author: jpoeppel """ -from __future__ import division +from __future__ import division import functools from ..nodes import UtilityNode from .factor import Factor +from primo2.nodes import DiscreteNode, DecisionNode, UtilityNode +import numpy as np + class VariableElimination(object): - + def __init__(self, decisionNetwork): self.net = decisionNetwork - + def get_decision(self, decisionNode, otherDecisions=None): """ """ pass - def _combine_factors(self, factors): """ Helper function to combine multiple joint factors. @@ -40,11 +42,12 @@ def _combine_factors(self, factors): A joint factor representing the combination of all the given factors. """ + def _combine_two_factors(f1, f2): - return (f1[0]*f2[0], f1[1]+f2[1]) - + return (f1[0] * f2[0], f1[1] + f2[1]) + return functools.reduce(_combine_two_factors, factors) - + def _marginalize_joint_factor(self, factor, variable): """ Helper function to marginalize a variable from a joint factor, @@ -64,7 +67,7 @@ def _marginalize_joint_factor(self, factor, variable): A joint factor with the given variable marginalized out. """ tmp = factor[0].marginalize(variable) - return (tmp, (factor[0]*factor[1]).marginalize(variable)/tmp) + return (tmp, (factor[0] * factor[1]).marginalize(variable) / tmp) def generalized_VE(self, joint_factors, elimination_variables): """ @@ -87,20 +90,20 @@ def generalized_VE(self, joint_factors, elimination_variables): utility factor over all the variables from the initial joint_factors which have not been eliminated. """ - + working_factors = set(joint_factors) for var in elimination_variables: relevant_factors = set([f for f in working_factors if var in f[0].variableOrder]) tmp_factor = self._combine_factors(relevant_factors) marginalized_factor = self._marginalize_joint_factor(tmp_factor, var) - + working_factors = (working_factors - relevant_factors) working_factors.add(marginalized_factor) - + res_factor = self._combine_factors(working_factors) - + return res_factor - + def expected_utility(self, decisions=None): """ Computes the expected utility fo the decision network, given the @@ -120,21 +123,20 @@ def expected_utility(self, decisions=None): """ if decisions is None: decisions = {} - #Set given decisions + # Set given decisions for decision in decisions: decisionNode = self.net.node_lookup[decision] -# decisionNode.clear() + # decisionNode.clear() decisionNode.set_decision(decisions[decision]) - - #Create joint factors + + # Create joint factors factors = set([]) for node in self.net.node_lookup.values(): factors.add(Factor.joint_factor(node)) - + eliminations = [node for node in self.net.node_lookup.values() if not isinstance(node, UtilityNode)] return self.generalized_VE(factors, eliminations)[1].potentials - - + def _optimize_locally(self, decisionNodeName): """ Helper function to choose the optimal decision for the given @@ -150,19 +152,18 @@ def _optimize_locally(self, decisionNodeName): string The optimal decision for this node """ - nodes = [node for node in self.net.node_lookup.keys() - if node != decisionNodeName and - not isinstance(self.net.node_lookup[node], UtilityNode) and - node not in self.net.node_lookup[decisionNodeName].parents] + nodes = [node for node in self.net.node_lookup.keys() + if node != decisionNodeName and + not isinstance(self.net.node_lookup[node], UtilityNode) and + node not in self.net.node_lookup[decisionNodeName].parents] decisionNode = self.net.node_lookup[decisionNodeName] - #Create joint factors + # Create joint factors factors = set([]) for node in self.net.node_lookup.values(): factors.add(Factor.joint_factor(node)) - + reducedFactor = self.generalized_VE(factors, nodes) - - + argmax = None maxVal = None for decision in decisionNode.values: @@ -170,10 +171,9 @@ def _optimize_locally(self, decisionNodeName): if maxVal is None or eu > maxVal: maxVal = eu argmax = decision - + return argmax - - + def get_optimal_decisions(self, decisionOrder, fixedDecisions=None): """ Iteratively computes the optimal decisions for all given decision @@ -200,20 +200,126 @@ def get_optimal_decisions(self, decisionOrder, fixedDecisions=None): """ if fixedDecisions is None: fixedDecisions = {} - - #Initialize random fully mixed strategy + + # Initialize random fully mixed strategy for decisionNode in decisionOrder: if not decisionNode in fixedDecisions: self.net.node_lookup[decisionNode].fully_mixed() else: self.net.node_lookup[decisionNode].set_decision(fixedDecisions[decisionNode]) - + solution = {} for decisionNode in decisionOrder: localDecision = self._optimize_locally(decisionNode) solution[decisionNode] = localDecision - + return solution + + @staticmethod + def inner_product(factors, decision_factors=None, utilities=None): + """ + Helper Function to multiply all factors in the given list. + This function is guided by the Algorithm 7.3.2 at page 112 in "Bayesian Reasoning and Machine Learning" + from David Barber. + Here the inner_product is calculated before the summing/maxing out using factors + + + Parameters + ---------- + factors: list + A list of factors to be multiplied with + + decision_factors: list + A list of decisions to be multiplied with the factors. Not none if the decisions + in the net have a decision rule (CPD) + + utilities: list + A list of utilities to be multiplied with the factors + + + + Returns + ------- + factor + The resulting (combined) factor + """ + prob_product = factors[0] + + for i, v in enumerate(factors): + if i != len(factors) - 1: + prob_product = prob_product * factors[i + 1] + combined_factors = prob_product + + if decision_factors: + decision_product = decision_factors[0] + + for i, v in enumerate(decision_factors): + if i != len(decision_factors) - 1: + decision_product = decision_product * decision_factors[i + 1] + combined_factors = combined_factors * decision_product + + if utilities: + utility_factors = [Factor.from_utility_node(i) for i in utilities] + utility_sum = utility_factors[0] + + for i, v in enumerate(utility_factors): + if i != len(utility_factors) - 1: + utility_sum = utility_sum + utility_factors[i + 1] + + combined_factors = combined_factors * utility_sum + + return combined_factors + + def max_sum(self, decisionNode): + """ + Max Sum Algorithm taken from the Algorithm 7.3.2 at page 112 in + "Bayesian Reasoning and Machine Learning" from David Barber. + + NEW: Now works with decision rules (That means the decisions have a CPD) + - - \ No newline at end of file + Parameters + ---------- + decisionNode: String + The optimal decision + + Returns + ------- + List + A list containing the optimal variable and corresponding optimal utility value of the resulting max sum algorithm + """ + + partialOrder = self.net.get_partial_ordering() + reverseOrder = partialOrder[::-1] + randomVariables = self.net.get_random_nodes() + utilities = self.net.get_utility_nodes() + decisions = self.net.get_decision_nodes() + + factors = [] + for node in randomVariables: + factors.append(Factor.from_node(node)) + + decision_factors = [] + + 'If All decision Nodes have a decision Rule' + if all(val.check_decision_rule() for val in decisions): + for node in decisions: + decision_factors.append(Factor.from_node(node)) + else: + decision_factors = None + + current = self.inner_product(factors, decision_factors, utilities) + for i in reverseOrder: + if isinstance(i, list): + if all(isinstance(self.net.node_lookup[val], DiscreteNode) for val in i) \ + or (all(isinstance(self.net.node_lookup[val], DecisionNode) for val in i) + and all([val.check_decision_rule() for val in i])): + current = current.marginalize(i) + else: + raise Exception("Marginalizing failed: Not all elements in the list are Discrete Nodes or " + "Decision Nodes with a specified decision rule") + elif isinstance(self.net.node_lookup[i], DiscreteNode): + current = current.marginalize(i) + elif i != decisionNode: + current = current.maximize(i) + return [current.values.values()[0][np.argmax(current.potentials)], max(current.potentials)] diff --git a/primo2/inference/factor.py b/primo2/inference/factor.py index e58c8e3..f274f73 100644 --- a/primo2/inference/factor.py +++ b/primo2/inference/factor.py @@ -508,7 +508,6 @@ def marginalize(self, variables): Factor A new factor where the given variables has been summed out. """ - if not isinstance(variables, (list,set)): variables = [variables] @@ -519,6 +518,37 @@ def marginalize(self, variables): res.variableOrder.remove(v) return res + + def maximize(self, variables): + """ + Function to create a new factor with the given variable maxed out + of this factor. + + Parameter + --------- + variables: String, RandomNode, [String,], [RandomNode,], set(String,) or set(RandomNode) + Either a single variable or a list of variables that are to + be removed. + Variables can either be addressed by their names or by the + nodes themselves. + + Returns + ------ + Factor + A new factor where the given variables has been summed out. + """ + + if not isinstance(variables, (list,set)): + variables = [variables] + + res = self.copy() + for v in variables: + res.potentials = np.amax(res.potentials, axis=res.variableOrder.index(v)) + del res.values[v] + res.variableOrder.remove(v) + + return res + def get_potential(self, variables=None): """ diff --git a/primo2/networks.py b/primo2/networks.py index 686b9a8..3282f01 100644 --- a/primo2/networks.py +++ b/primo2/networks.py @@ -23,6 +23,11 @@ from . import exceptions from . import nodes +from primo2.nodes import RandomNode, DiscreteNode, DecisionNode, UtilityNode +from math import exp, pi, sqrt +import pandas as pd +import numpy as np +from primo2.inference.factor import Factor class BayesianNetwork(object): @@ -34,7 +39,7 @@ def __init__(self): self.meta = [] # Used to be compatible with XMLBIF, stores properties def add_node(self, node): - if isinstance(node, nodes.RandomNode): + if isinstance(node, RandomNode): if node.name in self.node_lookup: raise ValueError("The network already contains a node " \ "called '{}'.".format(node.name)) @@ -42,7 +47,7 @@ def add_node(self, node): self.graph.add_node(node) else: raise TypeError("Only subclasses of RandomNode are valid nodes.") - + def remove_node(self, node): if node in self.graph: # Go over all children of this node @@ -50,7 +55,7 @@ def remove_node(self, node): child.remove_parent(self.node_lookup[node]) self.graph.remove_node(node) del self.node_lookup[node] - + def remove_edge(self, from_name, to_name): if from_name in self.graph and to_name in self.graph: self.node_lookup[to_name].remove_parent(self.node_lookup[from_name]) @@ -71,7 +76,7 @@ def get_node(self, node_name): except KeyError: raise Exception("There is no node with name {} in the " \ "BayesianNetwork".format(node_name)) - + def change_node_values(self, node, new_values): """ Updates the values of the given node. This will invalidate the node @@ -92,7 +97,7 @@ def change_node_values(self, node, new_values): else: raise Exception("There is no node with name {} in " \ "the network.".format(node)) - + def change_node_name(self, old_name, new_name): """ Renames the given node to the new name. @@ -107,27 +112,27 @@ def change_node_name(self, old_name, new_name): child.parents[new_name] = n idx = child.parentOrder.index(old_name) child.parentOrder[idx] = new_name - + # Fix nx graph self.graph.remove_node(n) - + n.name = new_name del self.node_lookup[old_name] self.node_lookup[new_name] = n - + self.graph.add_node(n) for child in children: self.graph.add_edge(n, self.node_lookup[child]) for parent in parents: self.graph.add_edge(self.node_lookup[parent], n) - + else: raise Exception("There is no node with name {} in the " \ "network.".format(old_name)) def get_all_nodes(self): return self.graph.nodes() - + def get_all_node_names(self): return self.node_lookup.keys() @@ -139,7 +144,7 @@ def get_nodes(self, node_names=None): for node_name in node_names: nodes.append(self.get_node(node_name)) return nodes - + def get_children(self, node_name): """ Returns a list of all the children of the given node. @@ -156,17 +161,17 @@ def get_children(self, node_name): node as parent. """ return self.graph.succ[node_name] - + def get_sample(self, evidence): sample = {} for e in evidence: sample[e] = evidence[e] # Make sure evidence contains RandomNodes and not only names for n in nx.topological_sort(self.graph): - + if n not in evidence: sample[n.name] = n.sample_value(sample, self.get_children(n), forward=True) - + return sample def clear(self): @@ -182,82 +187,476 @@ def number_of_nodes(self): def __len__(self): """Return the number of nodes in the graph.""" return len(self.graph) - -class DynamicBayesianNetwork(object): - """Class representing a the structure of a dynamic Bayesian network. - This temporal relationship is modelled as a 2-time-slice Bayesian - network (2-TBN; Koller & Friedman, 2009, §6.2.2). The Bayesian network - B_0 represents the initial distribution. The Bayesian network B_{->}, - a 2-TBN, represents the process. +# class DynamicBayesianNetwork(object): +# """Class representing a the structure of a dynamic Bayesian network. +# +# This temporal relationship is modelled as a 2-time-slice Bayesian +# network (2-TBN; Koller & Friedman, 2009, §6.2.2). The Bayesian network +# B_0 represents the initial distribution. The Bayesian network B_{->}, +# a 2-TBN, represents the process. +# """ +# +# def __init__(self, b0=None, two_tbn=None, transitions=None): +# """Create a dynamic Bayesian network. +# +# Parameters +# ---------- +# b0 : BayesianNetwork +# The network representing the initial distribution. +# two_tbn : BayesianNetwork +# The two-time-slice network representing the process. +# transitions : [(node, node_p),] +# A list of pairs, each of which represents one transition. +# See add_transition for more information. +# """ +# super(DynamicBayesianNetwork, self).__init__() +# self._b0 = BayesianNetwork() if b0 is None else b0 +# self._two_tbn = BayesianNetwork() if two_tbn is None else two_tbn +# self._transitions = [] +# if transitions is not None: +# self.add_transitions(transitions) +# +# @property +# def b0(self): +# """Get the Bayesian network B_0. +# +# Returns +# ------- +# BayesianNetwork +# The network representing the initial distribution. +# """ +# return self._b0 +# +# @b0.setter +# def b0(self, bn): +# """Set the Bayesian network B_0. +# +# Parameters +# ---------- +# bn : BayesianNetwork +# The network representing the initial distribution. +# """ +# self._b0 = bn +# +# @property +# def two_tbn(self): +# """Get the 2-TBN B_{->}. +# +# Returns +# ------- +# BayesianNetwork +# The two-time-slice network representing the process. +# """ +# return self._two_tbn +# +# @two_tbn.setter +# def two_tbn(self, bn): +# """Set the 2-TBN B_{->}. +# +# Parameters +# ---------- +# bn : BayesianNetwork +# The two-time-slice network representing the process. +# """ +# self._two_tbn = bn +# +# def add_transition(self, node, node_p): +# """Add a transition connecting nodes when unrolling the network. +# +# The transition is a directed edge from node X_i (`node`) to node X_i' +# (`node_p`). +# +# Parameters +# ---------- +# node : RandomNode, String +# The node X_i in the next time-slice. +# +# node_p: RandomNode, String +# The corresponding node X_i' in the current time-slice. +# +# Raises +# ------ +# primo2.exceptions.StructureError +# If a node specified in the transition model cannot be found in the +# corresponding networks. +# """ +# if node_p not in self._b0.get_all_nodes(): +# raise exceptions.StructureError( +# 'Node "{}" is not found in B_0 network "{}".'.format( +# node_p, self._b0.name)) +# if node_p not in self._two_tbn.get_all_nodes(): +# raise exceptions.StructureError( +# 'Node "{}" is not found in B_{->} network "{}".'.format( +# node_p, self._two_tbn.name)) +# if node not in self._two_tbn.get_all_nodes(): +# raise exceptions.StructureError( +# 'Node "{}" is not found in B_{{->}} network "{}".'.format( +# node, self._two_tbn.name)) +# self._transitions.append((node, node_p)) +# +# def add_transitions(self, transitions): +# """Add multiple transitions connecting nodes when unrolling the network. +# +# Parameters +# ---------- +# transitions : [(node, node_p),] +# A list of pairs, each of which represents one transition. +# See add_transition for more information. +# """ +# for transition in transitions: +# self.add_transition(transition[0], transition[1]) +# +# @property +# def transitions(self): +# """Get the transition model. +# +# Returns +# ------- +# [(node, node_p),] +# A list of pairs, each of which represents one transition. +# See add_transition for more information. +# """ +# return self._transitions +# +# +# class DecisionNetwork(object): +# +# def __init__(self): +# self.node_lookup = {} +# self.partialOrdering = [] +# self.random_nodes = [] +# self.decision_nodes = [] +# self.utility_nodes = [] +# +# def add_node(self, node): +# if isinstance(node, primo2.nodes.RandomNode): +# if node.name in self.node_lookup.keys(): +# raise Exception("Node name already exists in Bayesnet: " + node.name) +# if isinstance(node, primo2.nodes.DiscreteNode): +# self.random_nodes.append(node) +# elif isinstance(node, primo2.nodes.UtilityNode): +# self.utility_nodes.append(node) +# elif isinstance(node, primo2.nodes.DecisionNode): +# self.decision_nodes.append(node) +# else: +# raise Exception("Tried to add a node which the Bayesian Decision Network can not work with") +# self.node_lookup[node.name] = node +# else: +# raise Exception("Can only add 'Node' and its subclasses as nodes into the BayesianNetwork") +# +# def add_edge(self, node_from, node_to): +# +# self.node_lookup[node_to].add_parent(self.node_lookup[node_from]) +# +# def get_PartialOrdering(self): +# return self.partialOrdering +# +# def set_PartialOrdering(self, partialOrder): +# self.partialOrdering = partialOrder +# +# def get_all_nodes(self): +# '''Returns all RandomNodes''' +# return self.random_nodes +# +# def get_all_decision_nodes(self): +# return self.decision_nodes +# +# def get_all_utility_nodes(self): +# return self.utility_nodes + + +class DynamicDecisionNetwork(object): + """Class representing a the structure of a dynamic Decision network. + This temporal relationship is modelled as a 2-time-slice Decision + network. The Decision network + d0 represents the initial distribution. The Decision network D_{->}, + a 2-TDN, represents the process. """ - def __init__(self, b0=None, two_tbn=None, transitions=None): + def __init__(self, d0=None, two_tdn=None): """Create a dynamic Bayesian network. - Parameters ---------- - b0 : BayesianNetwork + d0 : d0_net The network representing the initial distribution. - two_tbn : BayesianNetwork + two_tdn : Two_TDN The two-time-slice network representing the process. - transitions : [(node, node_p),] - A list of pairs, each of which represents one transition. - See add_transition for more information. """ - super(DynamicBayesianNetwork, self).__init__() - self._b0 = BayesianNetwork() if b0 is None else b0 - self._two_tbn = BayesianNetwork() if two_tbn is None else two_tbn - self._transitions = [] - if transitions is not None: - self.add_transitions(transitions) - - @property - def b0(self): - """Get the Bayesian network B_0. - + super(DynamicDecisionNetwork, self).__init__() + + self._d0 = d0_net() if d0 is None else d0 + self._two_tdn = Two_TDN() if two_tdn is None else two_tdn + + @staticmethod + def update(node): + print ("In Update Method \n") + factor_parents_all_ = [] + factor_parents_persist = [] + node_factor = Factor.from_node(node) + print "Ordering of the parents \n" + for j in node.parents.values(): + factor_parents_all_.append(Factor.from_node(j)) + if node.get_unindexed_node() == j.get_unindexed_node(): + factor_parents_persist.append(Factor.from_node(node.parents[j])) + + prob_product = node_factor + print "Initial Shape" + print prob_product.potentials.shape + for i, _ in enumerate(factor_parents_all_): + print "We are multiplying with:", node.parents.keys()[i] + + if i != len(factor_parents_all_): + print "Shape after Multiplication" + prob_product = prob_product * factor_parents_all_[i] + print prob_product.potentials.shape + + print "Final shape" + print(prob_product.get_potential().shape) + prob_product = prob_product.marginalize(node.parentOrder) + print prob_product.potentials.shape + return prob_product.get_potential() + + def unroll(self, length): + + """Unrolling the network over specified length. Returns ------- - BayesianNetwork - The network representing the initial distribution. + DecisionNetwork + An unrolled DDN as a DecisionNetwork """ - return self._b0 - @b0.setter - def b0(self, bn): - """Set the Bayesian network B_0. + unrolled_net = DecisionNetwork() + + if length > 0: + for i in range(0, length): + if i == 0: + unrolled_net.add_nodes(self._d0.get_all_nodes()) + elif i == 1: + children = [] + """copy nodes from the two_tdn Network into unrolled_net""" + unrolled_net.copy_nodes_indexed(self._two_tdn.get_future_nodes(), i) + for inter_pair in self._two_tdn.get_inter_edges(): + parent_node_name = inter_pair[0].get_indexed_node(i - 1) + child_node_name = inter_pair[1].get_indexed_node(i) + """Connect copied nodes via inter edges like described in the two tdn""" + unrolled_net.add_edge(unrolled_net.node_lookup[parent_node_name], + unrolled_net.node_lookup[child_node_name]) + "Copy the CPD's describing the transition probability (defined by the two_tdn) into the new net" + for current_transition in self._two_tdn.transitions.keys(): + node_name = current_transition.get_indexed_node(i) + unrolled_net.node_lookup[node_name].set_cpd( + self._two_tdn.get_transition_probability(current_transition)) + + potential = self.update(unrolled_net.node_lookup[node_name]) + + for inter_pair in self._two_tdn.get_inter_edges(): + parent_node_name = inter_pair[0].get_indexed_node(i - 1) + child_node_name = inter_pair[1].get_indexed_node(i) + unrolled_net.remove_edge(unrolled_net.node_lookup[parent_node_name], + unrolled_net.node_lookup[child_node_name]) + unrolled_net.node_lookup[node_name].set_cpd(potential) + + for intra_pair in self._two_tdn.get_intra_edges(): + parent_node_name = intra_pair[0].get_indexed_node(i) + child_node_name = intra_pair[1].get_indexed_node(i) + """Connect copied nodes via intra edges like described in the two tdn""" + children.append(child_node_name) + unrolled_net.add_edge(unrolled_net.node_lookup[parent_node_name], + unrolled_net.node_lookup[child_node_name]) + "Copy the CPD's of the intra connections to the next timeslice" + for intra_child in children: + unrolled_net.node_lookup[intra_child].set_cpd( + self._d0.node_lookup[intra_child.split("_", 1)[0] + "_" + str(0)].cpd) + + else: + children = [] + unrolled_net.copy_nodes_indexed(self._two_tdn.get_future_nodes(), i) + # in case we dont want the decisions to be copied for every new timeslice but only for every + # two. Write this and leave the decisions node out of the timeslices in both d0 and two_tbn network: + # unrolled_net.copy_nodes_indexed(self._two_tdn.decision_nodes, i - 1) + + for inter_pair in self._two_tdn.get_inter_edges(): + parent_node_name = inter_pair[0].get_indexed_node(i - 1) + child_node_name = inter_pair[1].get_indexed_node(i) + unrolled_net.add_edge(unrolled_net.node_lookup[parent_node_name], + unrolled_net.node_lookup[child_node_name]) + + for current_transition in self._two_tdn.transitions.keys(): + node_name = current_transition.get_indexed_node(i) + + unrolled_net.node_lookup[node_name].set_cpd( + self._two_tdn.get_transition_probability(current_transition)) + + potential = self.update(unrolled_net.node_lookup[node_name]) + + for inter_pair in self._two_tdn.get_inter_edges(): + parent_node_name = inter_pair[0].get_indexed_node(i - 1) + child_node_name = inter_pair[1].get_indexed_node(i) + unrolled_net.remove_edge(unrolled_net.node_lookup[parent_node_name], + unrolled_net.node_lookup[child_node_name]) + unrolled_net.node_lookup[node_name].set_cpd(potential) + + for intra_pair in self._two_tdn.get_intra_edges(): + parent_node_name = intra_pair[0].get_indexed_node(i) + child_node_name = intra_pair[1].get_indexed_node(i) + children.append(child_node_name) + unrolled_net.add_edge(unrolled_net.node_lookup[parent_node_name], + unrolled_net.node_lookup[child_node_name]) + + for intra_child in children: + unrolled_net.node_lookup[intra_child].set_cpd( + self._d0.node_lookup[intra_child.split("_", 1)[0] + "_" + str(0)].cpd) + "Since we are unrolling over a specified length, the partial ordering for the DN has to be extended" + unrolled_net.extend_partial_ordering(self._d0.get_partial_ordering(), length) + else: + raise Exception("Can't unroll over length: 0") + return unrolled_net + + +class d0_net(object): + + def __init__(self): + self.utility_nodes = [] + self.decision_nodes = [] + self.random_nodes = [] + self.node_lookup = {} + self.partialOrdering = [] + + def add_nodes(self, nodes): + for i in nodes: + if isinstance(i, RandomNode): + if i.name in self.node_lookup.keys(): + pass + if isinstance(i, DiscreteNode): + self.random_nodes.append(i) + elif isinstance(i, UtilityNode): + self.utility_nodes.append(i) + elif isinstance(i, DecisionNode): + self.decision_nodes.append(i) + else: + raise Exception("Tried to add a node which the D0 Network can not work with") + self.node_lookup[i.name] = i + else: + raise Exception("Can only add 'Node' and its subclasses as nodes into the D0 Network") + + def add_edge(self, node_from, node_to): + self.node_lookup[node_to].add_parent(self.node_lookup[node_from]) + + def remove_edge(self, node_from, node_to): + self.node_lookup[node_to].remove_parent(self.node_lookup[node_from]) + + def get_partial_ordering(self): + return self.partialOrdering + + def set_partial_ordering(self, partial_order): + self.partialOrdering = partial_order + + def get_random_nodes(self): + return self.random_nodes + + def get_decision_nodes(self): + return self.decision_nodes + + def get_utility_nodes(self): + return self.utility_nodes + + def get_all_nodes(self): + return self.random_nodes + self.decision_nodes + self.utility_nodes + + def get_all_node_names(self): + return self.node_lookup.keys() + + +class Two_TDN(object): + def __init__(self, inter_edges=None, intra_edges=None): + + self.utility_nodes = [] + self.decision_nodes = [] + self.random_nodes = [] + self.node_lookup = {} + self._inter_edges = [] + self._intra_edges = [] + self.transitions = {} + if inter_edges is not None: + self.add_inter_edges(inter_edges) + if intra_edges is not None: + self.add_intra_edges(intra_edges) + + def add_nodes(self, nodes): + for i in nodes: + if isinstance(i, RandomNode): + if i.name in self.node_lookup.keys(): + pass + if isinstance(i, DiscreteNode): + self.random_nodes.append(i) + elif isinstance(i, UtilityNode): + self.utility_nodes.append(i) + elif isinstance(i, DecisionNode): + self.decision_nodes.append(i) + else: + raise Exception("Tried to add a node which the D0 Network can not work with") + self.node_lookup[i.name] = i + else: + raise Exception("Can only add 'Node' and its subclasses as nodes into the D0 Network") + + def get_random_nodes(self): + return self.random_nodes + + def get_decision_nodes(self): + return self.decision_nodes + + def get_utility_nodes(self): + return self.utility_nodes + + def get_all_nodes(self): + return self.random_nodes + self.decision_nodes + self.utility_nodes + + def get_all_node_names(self): + return self.node_lookup.keys() + + def add_inter_edge(self, node, node_p): + """Add a transition connecting nodes when unrolling the network. + + The transition is a directed edge from node X_i (`node`) to node X_i' + (`node_p`). Parameters ---------- - bn : BayesianNetwork - The network representing the initial distribution. - """ - self._b0 = bn + node : RandomNode, String + The node X_i in the next time-slice. - @property - def two_tbn(self): - """Get the 2-TBN B_{->}. + node_p: RandomNode, String + The corresponding node X_i' in the current time-slice. - Returns - ------- - BayesianNetwork - The two-time-slice network representing the process. + Raises + ------ + primo2.exceptions.StructureError + If a node specified in the transition model cannot be found in the + corresponding networks. """ - return self._two_tbn + if node_p not in self.get_all_nodes(): + raise exceptions.StructureError( + 'Node "{}" is not found in D_{->} network "{}".') + if node not in self.get_all_nodes(): + raise exceptions.StructureError( + 'Node "{}" is not found in D_{{->}} network "{}".') + self._inter_edges.append((node, node_p)) - @two_tbn.setter - def two_tbn(self, bn): - """Set the 2-TBN B_{->}. + def add_inter_edges(self, inter_edges): + """Add multiple transitions connecting nodes when unrolling the network. Parameters ---------- - bn : BayesianNetwork - The two-time-slice network representing the process. + inter_edges : [(node, node_p),] + A list of pairs, each of which represents one transition. + See add_transition for more information. """ - self._two_tbn = bn + for inter_edges in inter_edges: + self.add_inter_edge(inter_edges[0], inter_edges[1]) - def add_transition(self, node, node_p): + def add_intra_edge(self, node, node_p): """Add a transition connecting nodes when unrolling the network. The transition is a directed edge from node X_i (`node`) to node X_i' @@ -267,7 +666,7 @@ def add_transition(self, node, node_p): ---------- node : RandomNode, String The node X_i in the next time-slice. - + node_p: RandomNode, String The corresponding node X_i' in the current time-slice. @@ -277,51 +676,175 @@ def add_transition(self, node, node_p): If a node specified in the transition model cannot be found in the corresponding networks. """ - if node_p not in self._b0.get_all_nodes(): - raise exceptions.StructureError( - 'Node "{}" is not found in B_0 network "{}".'.format( - node_p, self._b0.name)) - if node_p not in self._two_tbn.get_all_nodes(): + if node_p not in self.get_all_nodes(): raise exceptions.StructureError( - 'Node "{}" is not found in B_{->} network "{}".'.format( - node_p, self._two_tbn.name)) - if node not in self._two_tbn.get_all_nodes(): + 'Node "{}" is not found in D_{->} network "{}".') + if node not in self.get_all_nodes(): raise exceptions.StructureError( - 'Node "{}" is not found in B_{{->}} network "{}".'.format( - node, self._two_tbn.name)) - self._transitions.append((node, node_p)) + 'Node "{}" is not found in D_{{->}} network "{}".') + self._intra_edges.append((node, node_p)) - def add_transitions(self, transitions): + def add_intra_edges(self, intra_edges): """Add multiple transitions connecting nodes when unrolling the network. Parameters ---------- - transitions : [(node, node_p),] - A list of pairs, each of which represents one transition. + intra_edges : [(node, node_p),] + A list of pairs, each of which represents one transition. See add_transition for more information. - """ - for transition in transitions: - self.add_transition(transition[0], transition[1]) + """ + for intra_edges in intra_edges: + self.add_intra_edge(intra_edges[0], intra_edges[1]) - @property - def transitions(self): - """Get the transition model. + def get_intra_edges(self): + return self._intra_edges - Returns - ------- - [(node, node_p),] - A list of pairs, each of which represents one transition. - See add_transition for more information. + def get_inter_edges(self): + return self._inter_edges + + def add_transition(self, node, transition_): + """ + Adds the transition probability from the child node of a two node pair connected through a persistent edge ( + connecting two timeslices) + :param node: The child node we want to assign the CPD to + :param transition_: The transition CPD """ - return self._transitions - + + self.transitions[node] = np.array(transition_) + + def get_transition_probability(self, node): + return self.transitions[node] + + def get_future_nodes(self): + return [self.node_lookup[x] for x in self.node_lookup.keys() if "t_plus_one" in x] + + class DecisionNetwork(object): - + def __init__(self): - self.node_lookup ={} - + self.node_lookup = {} + self.partialOrdering = [] + self.random_nodes = [] + self.decision_nodes = [] + self.utility_nodes = [] + def add_node(self, node): - self.node_lookup[node.name] = node - - def add_edge(self, from_name, to_name): - self.node_lookup[to_name].add_parent(self.node_lookup[from_name]) \ No newline at end of file + if isinstance(node, RandomNode): + if node.name in self.node_lookup.keys(): + pass + if isinstance(node, DiscreteNode): + self.random_nodes.append(node) + elif isinstance(node, UtilityNode): + self.utility_nodes.append(node) + elif isinstance(node, DecisionNode): + self.decision_nodes.append(node) + else: + raise Exception("Tried to add a node which the Bayesian Decision Network can not work with") + self.node_lookup[node.name] = node + else: + raise Exception("Can only add 'Node' and its subclasses as nodes into the BayesianNetwork") + + def add_nodes(self, nodes): + for i in nodes: + if isinstance(i, RandomNode): + if i.name in self.node_lookup.keys(): + pass + if isinstance(i, DiscreteNode): + self.random_nodes.append(i) + elif isinstance(i, UtilityNode): + self.utility_nodes.append(i) + elif isinstance(i, DecisionNode): + self.decision_nodes.append(i) + else: + raise Exception("Tried to add a node which the Bayesian Decision Network can not work with") + self.node_lookup[i.name] = i + else: + raise Exception("Can only add 'Node' and its subclasses as nodes into the BayesianNetwork") + + def copy_nodes_indexed(self, nodes_to_copy, index): + for i in nodes_to_copy: + if i.get_indexed_node(index) not in self.get_all_node_names(): + if isinstance(i, RandomNode): + if isinstance(i, DiscreteNode): + new_node = DiscreteNode(i.get_indexed_node(index), i.values) + self.random_nodes.append(new_node) + self.node_lookup[new_node.name] = new_node + elif isinstance(i, UtilityNode): + new_node = UtilityNode(i.get_indexed_node(index)) + self.utility_nodes.append(new_node) + self.node_lookup[new_node.name] = new_node + + elif isinstance(i, DecisionNode): + new_node = DecisionNode(i.get_indexed_node(index), i.values) + self.decision_nodes.append(new_node) + self.node_lookup[new_node.name] = new_node + else: + raise Exception("Tried to copy a node which the Bayesian Decision Network can not work with") + + else: + raise Exception("Can only add 'Node' and its subclasses as nodes into the BayesianNetwork") + else: + raise Exception("A copy of this node already exists") + + def add_edge(self, node_from, node_to): + + self.node_lookup[node_to].add_parent(self.node_lookup[node_from]) + + def remove_edge(self, node_from, node_to): + self.node_lookup[node_to].remove_parent(self.node_lookup[node_from]) + + def get_partial_ordering(self): + return self.partialOrdering + + def set_partial_ordering(self, partial_order): + self.partialOrdering = partial_order + + def extend_partial_ordering(self, init_order, length): + """ + "Extending Partial Order when unrolling network. The exact partial order is copied for the next timeslices" + :param init_order: The initial Partial Order which is described in the d0 network + :param length: The total length of the unrolled network + """ + for i in range(0, length): + for j in init_order: + # print(type(j)) + if isinstance(j, list): + if isinstance(self.node_lookup[j[0]], DiscreteNode) or \ + isinstance(self.node_lookup[j[0]], DecisionNode): + temp = [value.split("_", 1)[0] + "_" + str(i) for value in j] + self.partialOrdering.append(temp) + else: + raise ValueError("The type of node you wanted to extend the partial " + "Ordering with is not a DiscreteNode or a DecisionNode") + # Write this if decisions are over 2 timeslices + # elif i < length - 1 or i == 0 or decisions: + # temp = [value.split("_", 1)[0] + "_" + str(i) for value in j] + # self.partialOrdering.append(temp) + elif isinstance(j, str): + # print j + # print(self.node_lookup) + if isinstance(self.node_lookup[j], DiscreteNode) or \ + isinstance(self.node_lookup[j], DecisionNode): + self.partialOrdering.append(j.split("_", 1)[0] + "_" + str(i)) + + else: + raise ValueError("The type of node you wanted to extend the partial " + "Ordering with is not a DiscreteNode or a DecisionNode") + # Write this if decisions are over 2 timeslices + # elif i < length - 1 or i == 0: + # self.partialOrdering.append(j.split("_", 1)[0] + "_" + str(i)) + + def get_random_nodes(self): + return self.random_nodes + + def get_decision_nodes(self): + return self.decision_nodes + + def get_utility_nodes(self): + return self.utility_nodes + + def get_all_nodes(self): + return self.random_nodes + self.decision_nodes + self.utility_nodes + + def get_all_node_names(self): + return self.node_lookup.keys() diff --git a/primo2/nodes.py b/primo2/nodes.py index 40d05a9..4ac725f 100644 --- a/primo2/nodes.py +++ b/primo2/nodes.py @@ -22,6 +22,8 @@ import random import numpy as np +import pandas as pd +from math import sqrt, exp, pi class RandomNode(object): @@ -30,16 +32,24 @@ class RandomNode(object): Only subclasses of this should actually be used in a network as this is underspecied for any inference algorithms. """ - + def __init__(self, nodename): self.name = nodename self.cpd = 1 self.meta = [] - + def set_cpd(self, cpd): raise NotImplementedError("Called unimplemented method.") - - + + def get_unindexed_node(self): + return self.name.split("_", 1)[0] + + def get_index(self): + return self.name.split("_", 1)[1] + # Renames a node to the scheme: "Node_index" and returns it + def get_indexed_node(self, index): + return self.name.split("_", 1)[0] + "_" + str(index) + def __eq__(self, other): """ Two random nodes are considered to be identical if they have the @@ -51,10 +61,10 @@ def __eq__(self, other): return other.name == self.name except AttributeError: return other == self.name - + def __ne__(self, other): return not self.__eq__(other) - + def __hash__(self): """ The hash of a random node is the same as the hash of its name. @@ -62,15 +72,16 @@ def __hash__(self): instantiation or their name. """ return hash(self.name) - + def __str__(self): return self.name def __repr__(self): - return self.name + return self.name + class DiscreteNode(RandomNode): - + def __init__(self, nodename, values=("True", "False")): """ Creates a new discrete node with the given name and outcomes. @@ -82,13 +93,14 @@ def __init__(self, nodename, values=("True", "False")): vales: [String,], optional List of outcome values this discrete node has. """ - super(DiscreteNode,self).__init__(nodename) - self.values = list(values) #Consider copying this if for some reason something other than strings are used in there - self.parents = dict() #Consider removing this + super(DiscreteNode, self).__init__(nodename) + self.values = list( + values) # Consider copying this if for some reason something other than strings are used in there + self.parents = dict() # Consider removing this self.parentOrder = [] self.valid = False self._update_dimensions() - + def add_parent(self, parentNode): """ Adds the given node as a parent/cause node of this node. Will invalidate the @@ -102,11 +114,11 @@ def add_parent(self, parentNode): parentNode: RandomNode The new parent/cause node """ - + self.parents[parentNode.name] = parentNode self.parentOrder.append(parentNode.name) self._update_dimensions() - + def set_values(self, new_values): """ Allows to change/set the values of this variable. This will @@ -124,9 +136,8 @@ def set_values(self, new_values): List of the new value names. """ self.values = list(new_values) - self._update_dimensions() - - + self._update_dimensions() + def _update_dimensions(self): """ Private helper function to update the dimensions of the cpd. @@ -137,7 +148,7 @@ def _update_dimensions(self): dimensions.append(len(self.parents[parentName].values)) self.cpd = np.zeros(dimensions) self.valid = False - + def remove_parent(self, parentNode): """ Removes the given parent/cause node from this node. @@ -153,7 +164,7 @@ def remove_parent(self, parentNode): del self.parents[parentNode] self.parentOrder.remove(parentNode.name) self._update_dimensions() - + def set_cpd(self, cpd): """ Allows to set the conditional probability density(table) of this @@ -172,7 +183,7 @@ def set_cpd(self, cpd): raise ValueError("The dimensions of the given cpd do not match the dependency structure of the node.") self.cpd = np.copy(cpd) self.valid = True - + def set_probability(self, valueName, prob, parentValues=None): """ Function that allows to set parts of the cpt of this variable. @@ -202,10 +213,10 @@ def set_probability(self, valueName, prob, parentValues=None): index = [self.values.index(valueName)] if self.values else [] except ValueError: raise ValueError("This node as no value {}.".format(valueName)) - + if not parentValues: parentValues = {} - + for parentName in self.parentOrder: if parentName in parentValues: try: @@ -214,9 +225,9 @@ def set_probability(self, valueName, prob, parentValues=None): raise ValueError("Parent {} does not have values {}.".format(parentName, parentValues[parentName])) else: index.append(slice(len(self.parents[parentName].values))) - + self.cpd[tuple(index)] = prob - + def get_probability(self, value, parentValues=None): """ Function to return the probability(ies) for a given value of this @@ -248,31 +259,37 @@ def get_probability(self, value, parentValues=None): index = [[self.values.index(value)]] if self.values else [] except ValueError: raise ValueError("This node as no value {}.".format(value)) - + if not parentValues: parentValues = {} - + for parentName in self.parentOrder: if parentName in parentValues: if isinstance(parentValues[parentName], list): try: index.append([self.parents[parentName].values.index(v) for v in parentValues[parentName]]) except ValueError: - raise ValueError("There is no conditional probability for parent {}, values {} in node {}.".format(parentName, parentValues[parentName], self.name)) + raise ValueError( + "There is no conditional probability for parent {}, values {} in node {}.".format( + parentName, parentValues[parentName], self.name)) else: try: index.append([self.parents[parentName].values.index(parentValues[parentName])]) except ValueError: - raise ValueError("There is no conditional probability for parent {}, value {} in node {}.".format(parentName, parentValues[parentName], self.name)) + raise ValueError( + "There is no conditional probability for parent {}, value {} in node {}.".format(parentName, + parentValues[ + parentName], + self.name)) else: index.append(range(len(self.parents[parentName].values))) - + # use np.ix_ to construct the appropriate index array! index = np.ix_(*index) tmp = self.cpd[index] res = np.squeeze(tmp) - return res #np.squeeze(np.copy(self.cpd[np.ix_(*index])) - + return res # np.squeeze(np.copy(self.cpd[np.ix_(*index])) + def _get_single_probability(self, value, parentValues=None): """ Fast-path function to return the probability for a given value of this @@ -301,19 +318,24 @@ def _get_single_probability(self, value, parentValues=None): index = [self.values.index(value)] if self.values else [] except ValueError: raise ValueError("This node as no value {}.".format(value)) - + if not parentValues: parentValues = {} - + for parentName in self.parentOrder: try: index.append(self.parents[parentName].values.index(parentValues[parentName])) except KeyError: - raise KeyError("parentValues need to specify a value for parent {} of node: {}.".format(parentName, self.name)) + raise KeyError( + "parentValues need to specify a value for parent {} of node: {}.".format(parentName, self.name)) except ValueError: - raise ValueError("There is no conditional probability for parent {}, value {} in node {}.".format(parentName, parentValues[parentName], self.name)) + raise ValueError( + "There is no conditional probability for parent {}, value {} in node {}.".format(parentName, + parentValues[ + parentName], + self.name)) return self.cpd[tuple(index)] - + def get_markov_prob(self, outcome, children, state, forward=False): """ Computes the markov probability of the given outcome of this random variable, @@ -347,7 +369,7 @@ def get_markov_prob(self, outcome, children, state, forward=False): for child in children: prob *= child._get_single_probability(state[child.name], state) return prob - + def sample_value(self, currentState, children, forward=False): """ Returns a value drawn from the probability density given by this node @@ -374,21 +396,21 @@ def sample_value(self, currentState, children, forward=False): String Name of the value this DiscreteNode has most likely adopted. """ - #Compute probabilities of possible outcomes + # Compute probabilities of possible outcomes weights = [] for outcome in self.values: adaptedState = dict(currentState) adaptedState[self] = outcome weights.append(self.get_markov_prob(outcome, children, adaptedState, forward)) - - #Perform roulette-wheel-sampling: - rndVal = random.random()* sum(weights) + + # Perform roulette-wheel-sampling: + rndVal = random.random() * sum(weights) s = 0 for i in range(len(weights)): s += weights[i] if s >= rndVal: return self.values[i] - + def sample_local(self, currentValue): """ Returns a randomly chosen possible outcome of this node. @@ -414,14 +436,14 @@ class UtilityNode(RandomNode): 1. A DiscreteNode does not have values itself 2. The CPD does not represent probabilities but rather utilities! """ - + def __init__(self, nodeName): - super(UtilityNode,self).__init__(nodeName) + super(UtilityNode, self).__init__(nodeName) self.parents = {} self.parentOrder = [] - self.utilities = self.cpd # Use different name for cpd + self.utilities = self.cpd # Use different name for cpd self.valid = False - + def add_parent(self, parentNode): """ Adds the given node as a parent/cause node of this node. @@ -439,7 +461,7 @@ def add_parent(self, parentNode): self.parents[parentNode.name] = parentNode self.parentOrder.append(parentNode.name) self._update_dimensions() - + def _update_dimensions(self): """ Private helper function to update the dimensions of the cpd. @@ -450,7 +472,7 @@ def _update_dimensions(self): dimensions.append(len(self.parents[parentName].values)) self.cpd = np.zeros(dimensions) self.valid = False - + def get_utility(self, parentValues): """ Returns the utilities for a specific parent outcome. @@ -470,11 +492,15 @@ def get_utility(self, parentValues): try: index.append(self.parents[parentName].values.index(parentValues[parentName])) except KeyError: - raise KeyError("parentValues need to specify a value for parent {} of node: {}.".format(parentName, self.name)) + raise KeyError( + "parentValues need to specify a value for parent {} of node: {}.".format(parentName, self.name)) except ValueError: - raise ValueError("There is no utility for parent {}, value {} in node {}.".format(parentName, parentValues[parentName], self.name)) + raise ValueError("There is no utility for parent {}, value {} in node {}.".format(parentName, + parentValues[ + parentName], + self.name)) return self.cpd[tuple(index)] - + def set_utilities(self, utilities): """ Allows to set the utilities directly as a utility table. @@ -494,7 +520,7 @@ def set_utilities(self, utilities): "match the dependency structure of the node.") self.cpd = np.copy(utilities) self.valid = True - + def set_utility(self, utility, parentValues): """ Allows to specify the utility for the specified parent value @@ -518,9 +544,10 @@ def set_utility(self, utility, parentValues): raise ValueError("Parent {} does not have values {}.".format(parentName, parentValues[parentName])) else: index.append(slice(len(self.parents[parentName].values))) - + self.cpd[tuple(index)] = utility + class DecisionNode(RandomNode): """ A DecisionNode is a RandomNode where the cpd represents a decision-rule @@ -528,18 +555,18 @@ class DecisionNode(RandomNode): In most single-agent cases, this will be deterministic mapping to exactly 1 outcome. """ - - + def __init__(self, nodeName, decisions=["Yes", "No"]): - super(DecisionNode,self).__init__(nodeName) + super(DecisionNode, self).__init__(nodeName) self.name = nodeName self.values = list(decisions) self.state = None self.parentOrder = [] self.parents = {} + self.decision_rule = False self.deterministic = True self._update_dimensions() - + self.action_norm = {} def add_parent(self, parentNode): """ @@ -554,7 +581,7 @@ def add_parent(self, parentNode): parentNode: RandomNode The new parent/cause node """ - + self.parents[parentNode.name] = parentNode self.parentOrder.append(parentNode.name) self._update_dimensions() @@ -568,12 +595,12 @@ def _update_dimensions(self): for parentName in self.parentOrder: dimensions.append(len(self.parents[parentName].values)) self.cpd = np.zeros(dimensions) - self.valid = False - + self.decision_rule = False + def set_decision(self, decision): """ Sets the state of this decisionNode to the given decision. - This will equivalate a deterministic decision rule where only + This will equate a deterministic decision rule where only the given decision has a probability of 1. This invalidates any prior decision assignments. @@ -588,20 +615,64 @@ def set_decision(self, decision): index = [self.values.index(decision)] if self.values else [] except ValueError: raise ValueError("This node as no value {}.".format(decision)) - + self.cpd[tuple(index)] = 1 - + def fully_mixed(self): self.cpd = np.ones(self.cpd.shape) self.cpd /= len(self.values) - + + def set_action_norm(self, action_norm): + self.action_norm = pd.DataFrame(action_norm) + + def set_cpd(self, cpd): + """ + Allows to set the conditional probability density(table) of this + node directly. Will check that the dimensions of the given cpd + is comform to the current dependency structure but will not perform + any tests on the actual values. + + Parameters + ---------- + cpd : np.array + Table containing the conditional probabilities. Each variable + is represented by a dimension in the size of the number of its + outcomes. + """ + if np.shape(self.cpd) != np.shape(cpd): + raise ValueError("The dimensions of the given cpd do not match the dependency structure of the node.") + self.cpd = np.copy(cpd) + self.decision_rule = True + + def set_best_action(self, skill_believe): + """ + Invalidates the current CPD and sets the decisionNode to the action based on the previous skill_believe + :param skill_believe: The current skill believe + :return: The reduced CPD with only the best chosen action + """ + + def get_normal(x, mu, var): + return (1 / (sqrt(2 * pi * var))) * exp(-((x - mu) ** 2) / (2 * var)) + + best_value = 0 + best_action = None + + for i in self.action_norm: + current_value = get_normal(skill_believe, self.action_norm.loc[('loc', i)], + self.action_norm.loc[('scale', i)]) + if current_value > best_value: + best_value = current_value + best_action = i + + self.set_decision(best_action) + + def check_decision_rule(self): + return self.decision_rule + if __name__ == "__main__": - d1 = DiscreteNode("a") - d1.set_cpd(np.array([0.3,0.7])) - d2 = DiscreteNode("b", []) - d2.add_parent(d1) - d2.set_cpd(np.array([100,10])) - - print(d2.cpd) - \ No newline at end of file + costs = UtilityNode("costs") + costs_ = UtilityNode("costs_") + costs_.add_parent(costs) + + print(costs_.utilities) diff --git a/primo2/tests/Factor_test.py b/primo2/tests/Factor_test.py index a41b30d..4e41fdc 100644 --- a/primo2/tests/Factor_test.py +++ b/primo2/tests/Factor_test.py @@ -23,7 +23,7 @@ import unittest import numpy as np from primo2.inference.factor import Factor -from primo2.nodes import DiscreteNode, DecisionNode, UtilityNode +from primo2.nodes import DiscreteNode class FactorTest(unittest.TestCase): @@ -41,12 +41,6 @@ def setUp(self): cpt3 = np.array([0.8,0.2]) self.n3.set_cpd(cpt3) - self.un = UtilityNode("gains") - - self.un.add_parent(self.n3) - - self.un.set_utilities(np.array([100,10])) - def test_create_from_discrete_node_error(self): with self.assertRaises(TypeError) as cm: f = Factor.from_node("Node1") @@ -68,13 +62,6 @@ def test_create_from_discrete_node_with_parents(self): for v in self.n2.parents[p].values: self.assertTrue(v in f.values[p]) - - def test_create_from_utility_node(self): - f = Factor.from_utility_node(self.un) - - self.assertEqual(f.variableOrder, ["Node3"]) - np.testing.assert_array_equal(f.potentials, np.array([100,10])) - def test_create_from_samples(self): from collections import OrderedDict samples = [{"A":"True", "B":"True"}, {"A":"True", "B":"False"}, @@ -327,26 +314,5 @@ def test_normalize(self): self.assertEqual(np.sum(f1.potentials), 1.0) - def test_joint_factor_discrete_node(self): - (pF, uF) = Factor.joint_factor(self.n2) - - np.testing.assert_array_almost_equal(pF.potentials, self.n2.cpd) - np.testing.assert_array_almost_equal(uF.potentials, np.zeros(pF.potentials.shape)) - - def test_joint_factor_decision_node(self): - decisionNode = DecisionNode("dNode", decisions=["yes", "no"]) - (pF, uF) = Factor.joint_factor(decisionNode) - - np.testing.assert_array_almost_equal(pF.potentials, np.zeros(2)) - np.testing.assert_array_almost_equal(uF.potentials, np.zeros(pF.potentials.shape)) - - def test_joint_factor_utility_node(self): - utilityNode = UtilityNode("uNode") - utilityNode.set_utilities = np.array([10,100]) - (pF, uF) = Factor.joint_factor(utilityNode) - - np.testing.assert_array_almost_equal(pF.potentials, np.ones(uF.potentials.shape)) - np.testing.assert_array_almost_equal(uF.potentials, utilityNode.cpd) - if __name__ == "__main__": unittest.main() \ No newline at end of file diff --git a/primo2/tests/Node_test.py b/primo2/tests/Node_test.py index 64eb94b..735db21 100644 --- a/primo2/tests/Node_test.py +++ b/primo2/tests/Node_test.py @@ -403,74 +403,6 @@ def test_get_markov_prob(self): class ContinousNode(unittest.TestCase): pass -class UtilityNode(unittest.TestCase): - - def test_get_utility_error(self): - n = nodes.UtilityNode("Node1") - n2 = nodes.DiscreteNode("Node2", ["Value3", "Value4", "Value5"]) - n.add_parent(n2) - with self.assertRaises(ValueError) as cm: - n.get_utility({"Node2": ["Value6"]}) - self.assertEqual(str(cm.exception), "There is no utility for parent {}, value {} in node {}.".format("Node2", "['Value6']", "Node1")) - - - def test_get_utility_multiple_parents(self): - wet = nodes.UtilityNode("wet_grass") - sprink = nodes.DiscreteNode("sprinkler") - rain = nodes.DiscreteNode("rain") - wet.add_parent(sprink) - wet.add_parent(rain) - utilities = np.array([[10,20],[100,5]]) - wet.set_utilities(utilities) - np.testing.assert_array_almost_equal(wet.get_utility({"sprinkler":"True", "rain":"False"}),utilities[0,1]) - - def test_get_utility_single_parent_value(self): - n = nodes.UtilityNode("Node1") - n2 = nodes.DiscreteNode("Node2", ["Value3", "Value4", "Value5"]) - n.add_parent(n2) - utilities = np.array([5,10,35]) - n.set_utilities(utilities) - # Get only the specified probability - self.assertEqual(n.get_utility({"Node2": "Value4"}),10) - - def test_get_utility_single_parent_value_error(self): - n = nodes.UtilityNode("Node1") - n2 = nodes.DiscreteNode("Node2", ["Value3", "Value4", "Value5"]) - n.add_parent(n2) - utilities = np.array([5,10,35]) - n.set_utilities(utilities) - with self.assertRaises(ValueError) as cm: - n.get_utility({"Node2": "Value6"}) - self.assertEqual(str(cm.exception), "There is no utility for parent {}, value {} in node {}.".format("Node2", "Value6", "Node1")) - - def test_set_utilities_wrong_dimension(self): - wet = nodes.UtilityNode("wet_grass") - sprink = nodes.DiscreteNode("sprinkler") - rain = nodes.DiscreteNode("rain") - wet.add_parent(sprink) - wet.add_parent(rain) - utilities = np.array([[10,20],[100]]) - with self.assertRaises(ValueError) as cm: - wet.set_utilities(utilities) - self.assertEqual(str(cm.exception), "The dimensions of the given " \ - "utility table do not match the dependency structure of the node.") - - def test_set_utility(self): - n = nodes.UtilityNode("Node1") - n2 = nodes.DiscreteNode("Node2", ["Value3", "Value4", "Value5"]) - n.add_parent(n2) - n.set_utility(100, {"Node2":"Value3"}) - self.assertEqual(n.get_utility({"Node2": "Value3"}),100) - - - def test_set_utility_unknown_parentValue(self): - n = nodes.UtilityNode("Node1") - n2 = nodes.DiscreteNode("Node2", ["Value3", "Value4", "Value5"]) - n.add_parent(n2) - with self.assertRaises(ValueError) as cm: - n.set_utility(100, {"Node2":"Value6"}) - self.assertEqual(str(cm.exception), "Parent {} does not have values {}.".format("Node2", "Value6")) - """TODO for Decision Networks!!!""" class DecisionNode(unittest.TestCase): diff --git a/primo2/tests/tests.py b/primo2/tests/tests.py index 3108a1e..431e6b5 100644 --- a/primo2/tests/tests.py +++ b/primo2/tests/tests.py @@ -69,7 +69,7 @@ #print (f_sp * f_wgEm).marginalize("sprinkler").potentials print "f_wgEm variableorder: ", f_wgEm.variableOrder -print "f_wgEm variables: ", f_wgEm.variables +# print "f_wgEm variables: ", f_wgEm.variables print "f_wgEm: ", f_wgEm.potentials #print "f_sp: ", f_sp.potentials print "sp*wg: ", (f_sp * f_wgEm).potentials