From 37110de14d9a0d5c15de68f97a111a4376bbb6fb Mon Sep 17 00:00:00 2001 From: Donato Meoli Date: Mon, 29 Jul 2019 18:31:38 +0200 Subject: [PATCH 1/3] added map coloring SAT problem (#1092) * changed queue to set in AC3 Changed queue to set in AC3 (as in the pseudocode of the original algorithm) to reduce the number of consistency-check due to the redundancy of the same arcs in queue. For example, on the harder1 configuration of the Sudoku CSP the number consistency-check has been reduced from 40464 to 12562! * re-added test commented by mistake * added the mentioned AC4 algorithm for constraint propagation AC3 algorithm has non-optimal worst case time-complexity O(cd^3 ), while AC4 algorithm runs in O(cd^2) worst case time * added doctest in Sudoku for AC4 and and the possibility of choosing the constant propagation algorithm in mac inference * removed useless doctest for AC4 in Sudoku because AC4's tests are already present in test_csp.py * added map coloring SAT problems * fixed typo errors and removed unnecessary brackets * reformulated the map coloring problem * Revert "reformulated the map coloring problem" This reverts commit 20ab0e5afa238a0556e68f173b07ad32d0779d3b. * Revert "fixed typo errors and removed unnecessary brackets" This reverts commit f743146c43b28e0525b0f0b332faebc78c15946f. * Revert "added map coloring SAT problems" This reverts commit 9e0fa550e85081cf5b92fb6a3418384ab5a9fdfd. * Revert "removed useless doctest for AC4 in Sudoku because AC4's tests are already present in test_csp.py" This reverts commit b3cd24c511a82275f5b43c9f176396e6ba05f67e. * Revert "added doctest in Sudoku for AC4 and and the possibility of choosing the constant propagation algorithm in mac inference" This reverts commit 6986247481a05f1e558b93b2bf3cdae395f9c4ee. * Revert "added the mentioned AC4 algorithm for constraint propagation" This reverts commit 03551fbf2aa3980b915d4b6fefcbc70f24547b03. * added map coloring SAT problem * fixed build error * Revert "added map coloring SAT problem" This reverts commit 93af259e4811ddd775429f8a334111b9dd9e268c. * Revert "fixed build error" This reverts commit 6641c2c861728f3d43d3931ef201c6f7093cbc96. * added map coloring SAT problem * removed redundant parentheses --- csp.py | 73 ++++++++------ logic.py | 232 +++++++++++++++++++++++++++++--------------- tests/test_csp.py | 28 +++--- tests/test_logic.py | 56 ++++++----- 4 files changed, 244 insertions(+), 145 deletions(-) diff --git a/csp.py b/csp.py index ee59d4a6b..e1ee53a89 100644 --- a/csp.py +++ b/csp.py @@ -74,10 +74,12 @@ def unassign(self, var, assignment): def nconflicts(self, var, val, assignment): """Return the number of conflicts var=val has with other variables.""" + # Subclasses may implement this more efficiently def conflict(var2): return (var2 in assignment and not self.constraints(var, val, var2, assignment[var2])) + return count(conflict(v) for v in self.neighbors[var]) def display(self, assignment): @@ -153,6 +155,7 @@ def conflicted_vars(self, current): return [var for var in self.variables if self.nconflicts(var, current[var], current) > 0] + # ______________________________________________________________________________ # Constraint Propagation with AC-3 @@ -183,6 +186,7 @@ def revise(csp, Xi, Xj, removals): revised = True return revised + # ______________________________________________________________________________ # CSP Backtracking Search @@ -208,6 +212,7 @@ def num_legal_values(csp, var, assignment): return count(csp.nconflicts(var, val, assignment) == 0 for val in csp.domains[var]) + # Value ordering @@ -221,6 +226,7 @@ def lcv(var, assignment, csp): return sorted(csp.choices(var), key=lambda val: csp.nconflicts(var, val, assignment)) + # Inference @@ -245,6 +251,7 @@ def mac(csp, var, value, assignment, removals): """Maintain arc consistency.""" return AC3(csp, {(X, var) for X in csp.neighbors[var]}, removals) + # The search, proper @@ -274,6 +281,7 @@ def backtrack(assignment): assert result is None or csp.goal_test(result) return result + # ______________________________________________________________________________ # Min-conflicts hillclimbing search for CSPs @@ -302,6 +310,7 @@ def min_conflicts_value(csp, var, current): return argmin_random_tie(csp.domains[var], key=lambda val: csp.nconflicts(var, val, current)) + # ______________________________________________________________________________ @@ -356,7 +365,7 @@ def build_topological(node, parent, neighbors, visited, stack, parents): visited[node] = True for n in neighbors[node]: - if(not visited[n]): + if not visited[n]: build_topological(n, node, neighbors, visited, stack, parents) parents[node] = parent @@ -366,9 +375,9 @@ def build_topological(node, parent, neighbors, visited, stack, parents): def make_arc_consistent(Xj, Xk, csp): """Make arc between parent (Xj) and child (Xk) consistent under the csp's constraints, by removing the possible values of Xj that cause inconsistencies.""" - #csp.curr_domains[Xj] = [] + # csp.curr_domains[Xj] = [] for val1 in csp.domains[Xj]: - keep = False # Keep or remove val1 + keep = False # Keep or remove val1 for val2 in csp.domains[Xk]: if csp.constraints(Xj, val1, Xk, val2): # Found a consistent assignment for val1, keep it @@ -393,8 +402,9 @@ def assign_value(Xj, Xk, csp, assignment): # No consistent assignment available return None + # ______________________________________________________________________________ -# Map-Coloring Problems +# Map Coloring Problems class UniversalDict: @@ -446,27 +456,27 @@ def parse_neighbors(neighbors, variables=None): return dic -australia = MapColoringCSP(list('RGB'), - 'SA: WA NT Q NSW V; NT: WA Q; NSW: Q V; T: ') - -usa = MapColoringCSP(list('RGBY'), - """WA: OR ID; OR: ID NV CA; CA: NV AZ; NV: ID UT AZ; ID: MT WY UT; - UT: WY CO AZ; MT: ND SD WY; WY: SD NE CO; CO: NE KA OK NM; NM: OK TX AZ; - ND: MN SD; SD: MN IA NE; NE: IA MO KA; KA: MO OK; OK: MO AR TX; - TX: AR LA; MN: WI IA; IA: WI IL MO; MO: IL KY TN AR; AR: MS TN LA; - LA: MS; WI: MI IL; IL: IN KY; IN: OH KY; MS: TN AL; AL: TN GA FL; - MI: OH IN; OH: PA WV KY; KY: WV VA TN; TN: VA NC GA; GA: NC SC FL; - PA: NY NJ DE MD WV; WV: MD VA; VA: MD DC NC; NC: SC; NY: VT MA CT NJ; - NJ: DE; DE: MD; MD: DC; VT: NH MA; MA: NH RI CT; CT: RI; ME: NH; - HI: ; AK: """) - -france = MapColoringCSP(list('RGBY'), - """AL: LO FC; AQ: MP LI PC; AU: LI CE BO RA LR MP; BO: CE IF CA FC RA - AU; BR: NB PL; CA: IF PI LO FC BO; CE: PL NB NH IF BO AU LI PC; FC: BO - CA LO AL RA; IF: NH PI CA BO CE; LI: PC CE AU MP AQ; LO: CA AL FC; LR: - MP AU RA PA; MP: AQ LI AU LR; NB: NH CE PL BR; NH: PI IF CE NB; NO: - PI; PA: LR RA; PC: PL CE LI AQ; PI: NH NO CA IF; PL: BR NB CE PC; RA: - AU BO FC PA LR""") +australia_csp = MapColoringCSP(list('RGB'), """SA: WA NT Q NSW V; NT: WA Q; NSW: Q V; T: """) + +usa_csp = MapColoringCSP(list('RGBY'), + """WA: OR ID; OR: ID NV CA; CA: NV AZ; NV: ID UT AZ; ID: MT WY UT; + UT: WY CO AZ; MT: ND SD WY; WY: SD NE CO; CO: NE KA OK NM; NM: OK TX AZ; + ND: MN SD; SD: MN IA NE; NE: IA MO KA; KA: MO OK; OK: MO AR TX; + TX: AR LA; MN: WI IA; IA: WI IL MO; MO: IL KY TN AR; AR: MS TN LA; + LA: MS; WI: MI IL; IL: IN KY; IN: OH KY; MS: TN AL; AL: TN GA FL; + MI: OH IN; OH: PA WV KY; KY: WV VA TN; TN: VA NC GA; GA: NC SC FL; + PA: NY NJ DE MD WV; WV: MD VA; VA: MD DC NC; NC: SC; NY: VT MA CT NJ; + NJ: DE; DE: MD; MD: DC; VT: NH MA; MA: NH RI CT; CT: RI; ME: NH; + HI: ; AK: """) + +france_csp = MapColoringCSP(list('RGBY'), + """AL: LO FC; AQ: MP LI PC; AU: LI CE BO RA LR MP; BO: CE IF CA FC RA + AU; BR: NB PL; CA: IF PI LO FC BO; CE: PL NB NH IF BO AU LI PC; FC: BO + CA LO AL RA; IF: NH PI CA BO CE; LI: PC CE AU MP AQ; LO: CA AL FC; LR: + MP AU RA PA; MP: AQ LI AU LR; NB: NH CE PL BR; NH: PI IF CE NB; NO: + PI; PA: LR RA; PC: PL CE LI AQ; PI: NH NO CA IF; PL: BR NB CE PC; RA: + AU BO FC PA LR""") + # ______________________________________________________________________________ # n-Queens Problem @@ -503,16 +513,16 @@ def __init__(self, n): CSP.__init__(self, list(range(n)), UniversalDict(list(range(n))), UniversalDict(list(range(n))), queen_constraint) - self.rows = [0]*n - self.ups = [0]*(2*n - 1) - self.downs = [0]*(2*n - 1) + self.rows = [0] * n + self.ups = [0] * (2 * n - 1) + self.downs = [0] * (2 * n - 1) def nconflicts(self, var, val, assignment): """The number of conflicts, as recorded with each assignment. Count conflicts in row and in up, down diagonals. If there is a queen there, it can't conflict with itself, so subtract 3.""" n = len(self.variables) - c = self.rows[val] + self.downs[var+val] + self.ups[var-val+n-1] + c = self.rows[val] + self.downs[var + val] + self.ups[var - val + n - 1] if assignment.get(var, None) == val: c -= 3 return c @@ -560,6 +570,7 @@ def display(self, assignment): print(str(self.nconflicts(var, val, assignment)) + ch, end=' ') print() + # ______________________________________________________________________________ # Sudoku @@ -646,9 +657,12 @@ def show_cell(cell): return str(assignment.get(cell, '.')) def abut(lines1, lines2): return list( map(' | '.join, list(zip(lines1, lines2)))) + print('\n------+-------+------\n'.join( '\n'.join(reduce( abut, map(show_box, brow))) for brow in self.bgrid)) + + # ______________________________________________________________________________ # The Zebra Puzzle @@ -716,6 +730,7 @@ def zebra_constraint(A, a, B, b, recurse=0): (A in Smokes and B in Smokes)): return not same raise Exception('error') + return CSP(variables, domains, neighbors, zebra_constraint) diff --git a/logic.py b/logic.py index 6aacc4f95..24736c1a9 100644 --- a/logic.py +++ b/logic.py @@ -30,7 +30,7 @@ unify Do unification of two FOL sentences diff, simp Symbolic differentiation and simplification """ - +from csp import parse_neighbors, UniversalDict from utils import ( removeall, unique, first, argmax, probability, isnumber, issequence, Expr, expr, subexpressions @@ -42,11 +42,11 @@ import random from collections import defaultdict + # ______________________________________________________________________________ class KB: - """A knowledge base to which you can tell and ask sentences. To create a KB, first subclass this class and implement tell, ask_generator, and retract. Why ask_generator instead of ask? @@ -106,6 +106,7 @@ def retract(self, sentence): if c in self.clauses: self.clauses.remove(c) + # ______________________________________________________________________________ @@ -319,6 +320,7 @@ def pl_true(exp, model={}): else: raise ValueError("illegal operator in logic expression" + str(exp)) + # ______________________________________________________________________________ # Convert to Conjunctive Normal Form (CNF) @@ -368,6 +370,7 @@ def move_not_inwards(s): if s.op == '~': def NOT(b): return move_not_inwards(~b) + a = s.args[0] if a.op == '~': return move_not_inwards(a.args[0]) # ~~A ==> A @@ -445,6 +448,7 @@ def collect(subargs): collect(arg.args) else: result.append(arg) + collect(args) return result @@ -468,6 +472,7 @@ def disjuncts(s): """ return dissociate('|', [s]) + # ______________________________________________________________________________ @@ -481,7 +486,7 @@ def pl_resolution(KB, alpha): while True: n = len(clauses) pairs = [(clauses[i], clauses[j]) - for i in range(n) for j in range(i+1, n)] + for i in range(n) for j in range(i + 1, n)] for (ci, cj) in pairs: resolvents = pl_resolve(ci, cj) if False in resolvents: @@ -505,6 +510,7 @@ def pl_resolve(ci, cj): clauses.append(associate('|', dnew)) return clauses + # ______________________________________________________________________________ @@ -560,7 +566,6 @@ def pl_fc_entails(KB, q): """ wumpus_world_inference = expr("(B11 <=> (P12 | P21)) & ~B11") - """ [Figure 7.16] Propositional Logic Forward Chaining example """ @@ -572,9 +577,11 @@ def pl_fc_entails(KB, q): Definite clauses KB example """ definite_clauses_KB = PropDefiniteKB() -for clause in ['(B & F)==>E', '(A & E & F)==>G', '(B & C)==>F', '(A & B)==>D', '(E & F)==>H', '(H & I)==>J', 'A', 'B', 'C']: +for clause in ['(B & F)==>E', '(A & E & F)==>G', '(B & C)==>F', '(A & B)==>D', '(E & F)==>H', '(H & I)==>J', 'A', 'B', + 'C']: definite_clauses_KB.tell(expr(clause)) + # ______________________________________________________________________________ # DPLL-Satisfiable [Figure 7.17] @@ -665,7 +672,7 @@ def unit_clause_assign(clause, model): if model[sym] == positive: return None, None # clause already True elif P: - return None, None # more than 1 unbound variable + return None, None # more than 1 unbound variable else: P, value = sym, positive return P, value @@ -684,6 +691,7 @@ def inspect_literal(literal): else: return literal, True + # ______________________________________________________________________________ # Walk-SAT [Figure 7.18] @@ -714,95 +722,169 @@ def sat_count(sym): count = len([clause for clause in clauses if pl_true(clause, model)]) model[sym] = not model[sym] return count + sym = argmax(prop_symbols(clause), key=sat_count) model[sym] = not model[sym] # If no solution is found within the flip limit, we return failure return None + +# ______________________________________________________________________________ +# Map Coloring Problems + + +def MapColoringSAT(colors, neighbors): + """Make a SAT for the problem of coloring a map with different colors + for any two adjacent regions. Arguments are a list of colors, and a + dict of {region: [neighbor,...]} entries. This dict may also be + specified as a string of the form defined by parse_neighbors.""" + if isinstance(neighbors, str): + neighbors = parse_neighbors(neighbors) + colors = UniversalDict(colors) + clauses = [] + for state in neighbors.keys(): + clause = [expr(state + '_' + c) for c in colors[state]] + clauses.append(clause) + for t in itertools.combinations(clause, 2): + clauses.append([~t[0], ~t[1]]) + visited = set() + adj = set(neighbors[state]) - visited + visited.add(state) + for n_state in adj: + for col in colors[n_state]: + clauses.append([expr('~' + state + '_' + col), expr('~' + n_state + '_' + col)]) + return associate('&', map(lambda c: associate('|', c), clauses)) + + +australia_sat = MapColoringSAT(list('RGB'), """SA: WA NT Q NSW V; NT: WA Q; NSW: Q V; T: """) + +france_sat = MapColoringSAT(list('RGBY'), + """AL: LO FC; AQ: MP LI PC; AU: LI CE BO RA LR MP; BO: CE IF CA FC RA + AU; BR: NB PL; CA: IF PI LO FC BO; CE: PL NB NH IF BO AU LI PC; FC: BO + CA LO AL RA; IF: NH PI CA BO CE; LI: PC CE AU MP AQ; LO: CA AL FC; LR: + MP AU RA PA; MP: AQ LI AU LR; NB: NH CE PL BR; NH: PI IF CE NB; NO: + PI; PA: LR RA; PC: PL CE LI AQ; PI: NH NO CA IF; PL: BR NB CE PC; RA: + AU BO FC PA LR""") + +usa_sat = MapColoringSAT(list('RGBY'), + """WA: OR ID; OR: ID NV CA; CA: NV AZ; NV: ID UT AZ; ID: MT WY UT; + UT: WY CO AZ; MT: ND SD WY; WY: SD NE CO; CO: NE KA OK NM; NM: OK TX AZ; + ND: MN SD; SD: MN IA NE; NE: IA MO KA; KA: MO OK; OK: MO AR TX; + TX: AR LA; MN: WI IA; IA: WI IL MO; MO: IL KY TN AR; AR: MS TN LA; + LA: MS; WI: MI IL; IL: IN KY; IN: OH KY; MS: TN AL; AL: TN GA FL; + MI: OH IN; OH: PA WV KY; KY: WV VA TN; TN: VA NC GA; GA: NC SC FL; + PA: NY NJ DE MD WV; WV: MD VA; VA: MD DC NC; NC: SC; NY: VT MA CT NJ; + NJ: DE; DE: MD; MD: DC; VT: NH MA; MA: NH RI CT; CT: RI; ME: NH; + HI: ; AK: """) + + # ______________________________________________________________________________ # Expr functions for WumpusKB and HybridWumpusAgent -def facing_east (time): +def facing_east(time): return Expr('FacingEast', time) -def facing_west (time): + +def facing_west(time): return Expr('FacingWest', time) -def facing_north (time): + +def facing_north(time): return Expr('FacingNorth', time) -def facing_south (time): + +def facing_south(time): return Expr('FacingSouth', time) -def wumpus (x, y): + +def wumpus(x, y): return Expr('W', x, y) + def pit(x, y): return Expr('P', x, y) + def breeze(x, y): return Expr('B', x, y) + def stench(x, y): return Expr('S', x, y) + def wumpus_alive(time): return Expr('WumpusAlive', time) + def have_arrow(time): return Expr('HaveArrow', time) + def percept_stench(time): return Expr('Stench', time) + def percept_breeze(time): return Expr('Breeze', time) + def percept_glitter(time): return Expr('Glitter', time) + def percept_bump(time): return Expr('Bump', time) + def percept_scream(time): return Expr('Scream', time) + def move_forward(time): return Expr('Forward', time) + def shoot(time): return Expr('Shoot', time) + def turn_left(time): return Expr('TurnLeft', time) + def turn_right(time): return Expr('TurnRight', time) + def ok_to_move(x, y, time): return Expr('OK', x, y, time) -def location(x, y, time = None): + +def location(x, y, time=None): if time is None: return Expr('L', x, y) else: return Expr('L', x, y, time) + # Symbols def implies(lhs, rhs): return Expr('==>', lhs, rhs) + def equiv(lhs, rhs): return Expr('<=>', lhs, rhs) + # Helper Function def new_disjunction(sentences): t = sentences[0] - for i in range(1,len(sentences)): + for i in range(1, len(sentences)): t |= sentences[i] return t @@ -812,62 +894,59 @@ def new_disjunction(sentences): class WumpusKB(PropKB): """ - Create a Knowledge Base that contains the atemporal "Wumpus physics" and temporal rules with time zero. + Create a Knowledge Base that contains the a temporal "Wumpus physics" and temporal rules with time zero. """ - def __init__(self,dimrow): + def __init__(self, dimrow): super().__init__() self.dimrow = dimrow - self.tell( ~wumpus(1, 1) ) - self.tell( ~pit(1, 1) ) + self.tell(~wumpus(1, 1)) + self.tell(~pit(1, 1)) - for y in range(1, dimrow+1): - for x in range(1, dimrow+1): + for y in range(1, dimrow + 1): + for x in range(1, dimrow + 1): pits_in = list() wumpus_in = list() - if x > 1: # West room exists + if x > 1: # West room exists pits_in.append(pit(x - 1, y)) wumpus_in.append(wumpus(x - 1, y)) - if y < dimrow: # North room exists + if y < dimrow: # North room exists pits_in.append(pit(x, y + 1)) wumpus_in.append(wumpus(x, y + 1)) - if x < dimrow: # East room exists + if x < dimrow: # East room exists pits_in.append(pit(x + 1, y)) wumpus_in.append(wumpus(x + 1, y)) - if y > 1: # South room exists + if y > 1: # South room exists pits_in.append(pit(x, y - 1)) wumpus_in.append(wumpus(x, y - 1)) self.tell(equiv(breeze(x, y), new_disjunction(pits_in))) self.tell(equiv(stench(x, y), new_disjunction(wumpus_in))) - - ## Rule that describes existence of at least one Wumpus + # Rule that describes existence of at least one Wumpus wumpus_at_least = list() - for x in range(1, dimrow+1): + for x in range(1, dimrow + 1): for y in range(1, dimrow + 1): wumpus_at_least.append(wumpus(x, y)) self.tell(new_disjunction(wumpus_at_least)) - - ## Rule that describes existence of at most one Wumpus - for i in range(1, dimrow+1): - for j in range(1, dimrow+1): - for u in range(1, dimrow+1): - for v in range(1, dimrow+1): - if i!=u or j!=v: + # Rule that describes existence of at most one Wumpus + for i in range(1, dimrow + 1): + for j in range(1, dimrow + 1): + for u in range(1, dimrow + 1): + for v in range(1, dimrow + 1): + if i != u or j != v: self.tell(~wumpus(i, j) | ~wumpus(u, v)) - - ## Temporal rules at time zero + # Temporal rules at time zero self.tell(location(1, 1, 0)) - for i in range(1, dimrow+1): + for i in range(1, dimrow + 1): for j in range(1, dimrow + 1): self.tell(implies(location(i, j, 0), equiv(percept_breeze(0), breeze(i, j)))) self.tell(implies(location(i, j, 0), equiv(percept_stench(0), stench(i, j)))) @@ -881,7 +960,6 @@ def __init__(self,dimrow): self.tell(~facing_south(0)) self.tell(~facing_west(0)) - def make_action_sentence(self, action, time): actions = [move_forward(time), shoot(time), turn_left(time), turn_right(time)] @@ -895,7 +973,7 @@ def make_percept_sentence(self, percept, time): # Glitter, Bump, Stench, Breeze, Scream flags = [0, 0, 0, 0, 0] - ## Things perceived + # Things perceived if isinstance(percept, Glitter): flags[0] = 1 self.tell(percept_glitter(time)) @@ -912,7 +990,7 @@ def make_percept_sentence(self, percept, time): flags[4] = 1 self.tell(percept_scream(time)) - ## Things not perceived + # Things not perceived for i in range(len(flags)): if flags[i] == 0: if i == 0: @@ -926,15 +1004,14 @@ def make_percept_sentence(self, percept, time): elif i == 4: self.tell(~percept_scream(time)) - def add_temporal_sentences(self, time): if time == 0: return t = time - 1 - ## current location rules - for i in range(1, self.dimrow+1): - for j in range(1, self.dimrow+1): + # current location rules + for i in range(1, self.dimrow + 1): + for j in range(1, self.dimrow + 1): self.tell(implies(location(i, j, time), equiv(percept_breeze(time), breeze(i, j)))) self.tell(implies(location(i, j, time), equiv(percept_stench(time), stench(i, j)))) @@ -956,15 +1033,15 @@ def add_temporal_sentences(self, time): if j != self.dimrow: s.append(location(i, j + 1, t) & facing_south(t) & move_forward(t)) - ## add sentence about location i,j + # add sentence about location i,j self.tell(new_disjunction(s)) - ## add sentence about safety of location i,j + # add sentence about safety of location i,j self.tell( equiv(ok_to_move(i, j, time), ~pit(i, j) & ~wumpus(i, j) & wumpus_alive(time)) ) - ## Rules about current orientation + # Rules about current orientation a = facing_north(t) & turn_right(t) b = facing_south(t) & turn_left(t) @@ -990,16 +1067,15 @@ def add_temporal_sentences(self, time): s = equiv(facing_south(time), a | b | c) self.tell(s) - ## Rules about last action + # Rules about last action self.tell(equiv(move_forward(t), ~turn_right(t) & ~turn_left(t))) - ##Rule about the arrow + # Rule about the arrow self.tell(equiv(have_arrow(time), have_arrow(t) & ~shoot(t))) - ##Rule about Wumpus (dead or alive) + # Rule about Wumpus (dead or alive) self.tell(equiv(wumpus_alive(time), wumpus_alive(t) & ~percept_scream(time))) - def ask_if_true(self, query): return pl_resolution(self, query) @@ -1007,13 +1083,12 @@ def ask_if_true(self, query): # ______________________________________________________________________________ -class WumpusPosition(): +class WumpusPosition: def __init__(self, x, y, orientation): self.X = x self.Y = y self.orientation = orientation - def get_location(self): return self.X, self.Y @@ -1029,18 +1104,19 @@ def set_orientation(self, orientation): def __eq__(self, other): if other.get_location() == self.get_location() and \ - other.get_orientation()==self.get_orientation(): + other.get_orientation() == self.get_orientation(): return True else: return False + # ______________________________________________________________________________ class HybridWumpusAgent(Agent): """An agent for the wumpus world that does logical inference. [Figure 7.20]""" - def __init__(self,dimentions): + def __init__(self, dimentions): self.dimrow = dimentions self.kb = WumpusKB(self.dimrow) self.t = 0 @@ -1048,15 +1124,14 @@ def __init__(self,dimentions): self.current_position = WumpusPosition(1, 1, 'UP') super().__init__(self.execute) - def execute(self, percept): self.kb.make_percept_sentence(percept, self.t) self.kb.add_temporal_sentences(self.t) temp = list() - for i in range(1, self.dimrow+1): - for j in range(1, self.dimrow+1): + for i in range(1, self.dimrow + 1): + for j in range(1, self.dimrow + 1): if self.kb.ask_if_true(location(i, j, self.t)): temp.append(i) temp.append(j) @@ -1071,8 +1146,8 @@ def execute(self, percept): self.current_position = WumpusPosition(temp[0], temp[1], 'RIGHT') safe_points = list() - for i in range(1, self.dimrow+1): - for j in range(1, self.dimrow+1): + for i in range(1, self.dimrow + 1): + for j in range(1, self.dimrow + 1): if self.kb.ask_if_true(ok_to_move(i, j, self.t)): safe_points.append([i, j]) @@ -1080,14 +1155,14 @@ def execute(self, percept): goals = list() goals.append([1, 1]) self.plan.append('Grab') - actions = self.plan_route(self.current_position,goals,safe_points) + actions = self.plan_route(self.current_position, goals, safe_points) self.plan.extend(actions) self.plan.append('Climb') if len(self.plan) == 0: unvisited = list() - for i in range(1, self.dimrow+1): - for j in range(1, self.dimrow+1): + for i in range(1, self.dimrow + 1): + for j in range(1, self.dimrow + 1): for k in range(self.t): if self.kb.ask_if_true(location(i, j, k)): unvisited.append([i, j]) @@ -1097,13 +1172,13 @@ def execute(self, percept): if u not in unvisited_and_safe and s == u: unvisited_and_safe.append(u) - temp = self.plan_route(self.current_position,unvisited_and_safe,safe_points) + temp = self.plan_route(self.current_position, unvisited_and_safe, safe_points) self.plan.extend(temp) if len(self.plan) == 0 and self.kb.ask_if_true(have_arrow(self.t)): possible_wumpus = list() - for i in range(1, self.dimrow+1): - for j in range(1, self.dimrow+1): + for i in range(1, self.dimrow + 1): + for j in range(1, self.dimrow + 1): if not self.kb.ask_if_true(wumpus(i, j)): possible_wumpus.append([i, j]) @@ -1112,8 +1187,8 @@ def execute(self, percept): if len(self.plan) == 0: not_unsafe = list() - for i in range(1, self.dimrow+1): - for j in range(1, self.dimrow+1): + for i in range(1, self.dimrow + 1): + for j in range(1, self.dimrow + 1): if not self.kb.ask_if_true(ok_to_move(i, j, self.t)): not_unsafe.append([i, j]) temp = self.plan_route(self.current_position, not_unsafe, safe_points) @@ -1133,19 +1208,17 @@ def execute(self, percept): return action - def plan_route(self, current, goals, allowed): problem = PlanRoute(current, goals, allowed, self.dimrow) return astar_search(problem).solution() - def plan_shot(self, current, goals, allowed): shooting_positions = set() for loc in goals: x = loc[0] y = loc[1] - for i in range(1, self.dimrow+1): + for i in range(1, self.dimrow + 1): if i < x: shooting_positions.add(WumpusPosition(i, y, 'EAST')) if i > x: @@ -1157,7 +1230,7 @@ def plan_shot(self, current, goals, allowed): # Can't have a shooting position from any of the rooms the Wumpus could reside orientations = ['EAST', 'WEST', 'NORTH', 'SOUTH'] - for loc in goals: + for loc in goals: for orientation in orientations: shooting_positions.remove(WumpusPosition(loc[0], loc[1], orientation)) @@ -1186,7 +1259,7 @@ def translate_to_SAT(init, transition, goal, time): # Symbol claiming state s at time t state_counter = itertools.count() for s in states: - for t in range(time+1): + for t in range(time + 1): state_sym[s, t] = Expr("State_{}".format(next(state_counter))) # Add initial state axiom @@ -1206,11 +1279,11 @@ def translate_to_SAT(init, transition, goal, time): "Transition_{}".format(next(transition_counter))) # Change the state from s to s_ - clauses.append(action_sym[s, action, t] |'==>'| state_sym[s, t]) - clauses.append(action_sym[s, action, t] |'==>'| state_sym[s_, t + 1]) + clauses.append(action_sym[s, action, t] | '==>' | state_sym[s, t]) + clauses.append(action_sym[s, action, t] | '==>' | state_sym[s_, t + 1]) # Allow only one state at any time - for t in range(time+1): + for t in range(time + 1): # must be a state at any time clauses.append(associate('|', [state_sym[s, t] for s in states])) @@ -1363,6 +1436,7 @@ def standardize_variables(sentence, dic=None): standardize_variables.counter = itertools.count() + # ______________________________________________________________________________ @@ -1404,6 +1478,7 @@ def fol_fc_ask(KB, alpha): """A simple forward-chaining algorithm. [Figure 9.3]""" # TODO: Improve efficiency kb_consts = list({c for clause in KB.clauses for c in constant_symbols(clause)}) + def enum_subst(p): query_vars = list({v for clause in p for v in variables(clause)}) for assignment_list in itertools.product(kb_consts, repeat=len(query_vars)): @@ -1466,8 +1541,8 @@ def fol_bc_and(KB, goals, theta): P11, P12, P21, P22, P31, B11, B21 = expr('P11, P12, P21, P22, P31, B11, B21') wumpus_kb.tell(~P11) -wumpus_kb.tell(B11 | '<=>' | ((P12 | P21))) -wumpus_kb.tell(B21 | '<=>' | ((P11 | P22 | P31))) +wumpus_kb.tell(B11 | '<=>' | (P12 | P21)) +wumpus_kb.tell(B21 | '<=>' | (P11 | P22 | P31)) wumpus_kb.tell(~B11) wumpus_kb.tell(B21) @@ -1497,6 +1572,7 @@ def fol_bc_and(KB, goals, theta): 'Enemy(Nono, America)' ])) + # ______________________________________________________________________________ # Example application (not in the book). @@ -1527,7 +1603,7 @@ def diff(y, x): elif op == '/': return (v * diff(u, x) - u * diff(v, x)) / (v * v) elif op == '**' and isnumber(x.op): - return (v * u ** (v - 1) * diff(u, x)) + return v * u ** (v - 1) * diff(u, x) elif op == '**': return (v * u ** (v - 1) * diff(u, x) + u ** v * Expr('log')(u) * diff(v, x)) diff --git a/tests/test_csp.py b/tests/test_csp.py index c34d42540..a7564a395 100644 --- a/tests/test_csp.py +++ b/tests/test_csp.py @@ -10,16 +10,16 @@ def test_csp_assign(): var = 10 val = 5 assignment = {} - australia.assign(var, val, assignment) + australia_csp.assign(var, val, assignment) - assert australia.nassigns == 1 + assert australia_csp.nassigns == 1 assert assignment[var] == val def test_csp_unassign(): var = 10 assignment = {var: 5} - australia.unassign(var, assignment) + australia_csp.unassign(var, assignment) assert var not in assignment @@ -330,22 +330,22 @@ def test_forward_checking(): def test_backtracking_search(): - assert backtracking_search(australia) - assert backtracking_search(australia, select_unassigned_variable=mrv) - assert backtracking_search(australia, order_domain_values=lcv) - assert backtracking_search(australia, select_unassigned_variable=mrv, + assert backtracking_search(australia_csp) + assert backtracking_search(australia_csp, select_unassigned_variable=mrv) + assert backtracking_search(australia_csp, order_domain_values=lcv) + assert backtracking_search(australia_csp, select_unassigned_variable=mrv, order_domain_values=lcv) - assert backtracking_search(australia, inference=forward_checking) - assert backtracking_search(australia, inference=mac) - assert backtracking_search(usa, select_unassigned_variable=mrv, + assert backtracking_search(australia_csp, inference=forward_checking) + assert backtracking_search(australia_csp, inference=mac) + assert backtracking_search(usa_csp, select_unassigned_variable=mrv, order_domain_values=lcv, inference=mac) def test_min_conflicts(): - assert min_conflicts(australia) - assert min_conflicts(france) + assert min_conflicts(australia_csp) + assert min_conflicts(france_csp) - tests = [(usa, None)] * 3 + tests = [(usa_csp, None)] * 3 assert failure_test(min_conflicts, tests) >= 1 / 3 australia_impossible = MapColoringCSP(list('RG'), 'SA: WA NT Q NSW V; NT: WA Q; NSW: Q V; T: ') @@ -418,7 +418,7 @@ def test_parse_neighbours(): def test_topological_sort(): root = 'NT' - Sort, Parents = topological_sort(australia, root) + Sort, Parents = topological_sort(australia_csp, root) assert Sort == ['NT', 'SA', 'Q', 'NSW', 'V', 'WA'] assert Parents['NT'] == None diff --git a/tests/test_logic.py b/tests/test_logic.py index 378f1f0fc..fe9a9c5e3 100644 --- a/tests/test_logic.py +++ b/tests/test_logic.py @@ -1,10 +1,12 @@ import pytest + from logic import * -from utils import expr_handle_infix_ops, count, Symbol +from utils import expr_handle_infix_ops, count definite_clauses_KB = PropDefiniteKB() -for clause in ['(B & F)==>E', '(A & E & F)==>G', '(B & C)==>F', '(A & B)==>D', '(E & F)==>H', '(H & I)==>J', 'A', 'B', 'C']: - definite_clauses_KB.tell(expr(clause)) +for clause in ['(B & F)==>E', '(A & E & F)==>G', '(B & C)==>F', '(A & B)==>D', '(E & F)==>H', '(H & I)==>J', 'A', 'B', + 'C']: + definite_clauses_KB.tell(expr(clause)) def test_is_symbol(): @@ -47,7 +49,7 @@ def test_extend(): def test_subst(): - assert subst({x: 42, y:0}, F(x) + y) == (F(42) + 0) + assert subst({x: 42, y: 0}, F(x) + y) == (F(42) + 0) def test_PropKB(): @@ -55,7 +57,7 @@ def test_PropKB(): assert count(kb.ask(expr) for expr in [A, C, D, E, Q]) is 0 kb.tell(A & E) assert kb.ask(A) == kb.ask(E) == {} - kb.tell(E |'==>'| C) + kb.tell(E | '==>' | C) assert kb.ask(C) == {} kb.retract(E) assert kb.ask(E) is False @@ -94,14 +96,15 @@ def test_is_definite_clause(): def test_parse_definite_clause(): assert parse_definite_clause(expr('A & B & C & D ==> E')) == ([A, B, C, D], E) assert parse_definite_clause(expr('Farmer(Mac)')) == ([], expr('Farmer(Mac)')) - assert parse_definite_clause(expr('(Farmer(f) & Rabbit(r)) ==> Hates(f, r)')) == ([expr('Farmer(f)'), expr('Rabbit(r)')], expr('Hates(f, r)')) + assert parse_definite_clause(expr('(Farmer(f) & Rabbit(r)) ==> Hates(f, r)')) == ( + [expr('Farmer(f)'), expr('Rabbit(r)')], expr('Hates(f, r)')) def test_pl_true(): assert pl_true(P, {}) is None assert pl_true(P, {P: False}) is False - assert pl_true(P | Q, {P: True}) is True - assert pl_true((A | B) & (C | D), {A: False, B: True, D: True}) is True + assert pl_true(P | Q, {P: True}) + assert pl_true((A | B) & (C | D), {A: False, B: True, D: True}) assert pl_true((A & B) & (C | D), {A: False, B: True, D: True}) is False assert pl_true((A & B) | (A & C), {A: False, B: True, C: True}) is False assert pl_true((A | B) & (C | D), {A: True, D: False}) is None @@ -131,28 +134,28 @@ def test_dpll(): assert (dpll_satisfiable(A & ~B & C & (A | ~D) & (~E | ~D) & (C | ~D) & (~A | ~F) & (E | ~F) & (~D | ~F) & (B | ~C | D) & (A | ~E | F) & (~A | E | D)) == {B: False, C: True, A: True, F: False, D: True, E: False}) - assert dpll_satisfiable(A & B & ~C & D) == {C: False, A: True, D: True, B: True} - assert dpll_satisfiable((A | (B & C)) |'<=>'| ((A | B) & (A | C))) == {C: True, A: True} or {C: True, B: True} - assert dpll_satisfiable(A |'<=>'| B) == {A: True, B: True} + assert dpll_satisfiable(A & B & ~C & D) == {C: False, A: True, D: True, B: True} + assert dpll_satisfiable((A | (B & C)) | '<=>' | ((A | B) & (A | C))) == {C: True, A: True} or {C: True, B: True} + assert dpll_satisfiable(A | '<=>' | B) == {A: True, B: True} assert dpll_satisfiable(A & ~B) == {A: True, B: False} assert dpll_satisfiable(P & ~P) is False def test_find_pure_symbol(): - assert find_pure_symbol([A, B, C], [A|~B,~B|~C,C|A]) == (A, True) - assert find_pure_symbol([A, B, C], [~A|~B,~B|~C,C|A]) == (B, False) - assert find_pure_symbol([A, B, C], [~A|B,~B|~C,C|A]) == (None, None) + assert find_pure_symbol([A, B, C], [A | ~B, ~B | ~C, C | A]) == (A, True) + assert find_pure_symbol([A, B, C], [~A | ~B, ~B | ~C, C | A]) == (B, False) + assert find_pure_symbol([A, B, C], [~A | B, ~B | ~C, C | A]) == (None, None) def test_unit_clause_assign(): - assert unit_clause_assign(A|B|C, {A:True}) == (None, None) - assert unit_clause_assign(B|C, {A:True}) == (None, None) - assert unit_clause_assign(B|~A, {A:True}) == (B, True) + assert unit_clause_assign(A | B | C, {A: True}) == (None, None) + assert unit_clause_assign(B | C, {A: True}) == (None, None) + assert unit_clause_assign(B | ~A, {A: True}) == (B, True) def test_find_unit_clause(): - assert find_unit_clause([A|B|C, B|~C, ~A|~B], {A:True}) == (B, False) - + assert find_unit_clause([A | B | C, B | ~C, ~A | ~B], {A: True}) == (B, False) + def test_unify(): assert unify(x, x, {}) == {} @@ -175,9 +178,9 @@ def test_tt_entails(): assert tt_entails(P & Q, Q) assert not tt_entails(P | Q, Q) assert tt_entails(A & (B | C) & E & F & ~(P | Q), A & E & F & ~P & ~Q) - assert not tt_entails(P |'<=>'| Q, Q) - assert tt_entails((P |'==>'| Q) & P, Q) - assert not tt_entails((P |'<=>'| Q) & ~P, Q) + assert not tt_entails(P | '<=>' | Q, Q) + assert tt_entails((P | '==>' | Q) & P, Q) + assert not tt_entails((P | '<=>' | Q) & ~P, Q) def test_prop_symbols(): @@ -231,12 +234,13 @@ def test_move_not_inwards(): def test_distribute_and_over_or(): - def test_entailment(s, has_and = False): + def test_entailment(s, has_and=False): result = distribute_and_over_or(s) if has_and: assert result.op == '&' assert tt_entails(s, result) assert tt_entails(result, s) + test_entailment((A & B) | C, True) test_entailment((A | B) & C, True) test_entailment((A | B) | C, False) @@ -253,7 +257,8 @@ def test_to_cnf(): assert repr(to_cnf("a | (b & c) | d")) == '((b | a | d) & (c | a | d))' assert repr(to_cnf("A & (B | (D & E))")) == '(A & (D | B) & (E | B))' assert repr(to_cnf("A | (B | (C | (D & E)))")) == '((D | A | B | C) & (E | A | B | C))' - assert repr(to_cnf('(A <=> ~B) ==> (C | ~D)')) == '((B | ~A | C | ~D) & (A | ~A | C | ~D) & (B | ~B | C | ~D) & (A | ~B | C | ~D))' + assert repr(to_cnf( + '(A <=> ~B) ==> (C | ~D)')) == '((B | ~A | C | ~D) & (A | ~A | C | ~D) & (B | ~B | C | ~D) & (A | ~B | C | ~D))' def test_pl_resolution(): @@ -281,6 +286,7 @@ def test_ask(query, kb=None): return sorted( [dict((x, v) for x, v in list(a.items()) if x in test_variables) for a in answers], key=repr) + assert repr(test_ask('Farmer(x)')) == '[{x: Mac}]' assert repr(test_ask('Human(x)')) == '[{x: Mac}, {x: MrsMac}]' assert repr(test_ask('Rabbit(x)')) == '[{x: MrsRabbit}, {x: Pete}]' @@ -295,6 +301,7 @@ def test_ask(query, kb=None): return sorted( [dict((x, v) for x, v in list(a.items()) if x in test_variables) for a in answers], key=repr) + assert repr(test_ask('Criminal(x)', crime_kb)) == '[{x: West}]' assert repr(test_ask('Enemy(x, America)', crime_kb)) == '[{x: Nono}]' assert repr(test_ask('Farmer(x)')) == '[{x: Mac}]' @@ -316,6 +323,7 @@ def check_SAT(clauses, single_solution={}): if single_solution: # Cross check the solution if only one exists assert all(pl_true(x, single_solution) for x in clauses) assert soln == single_solution + # Test WalkSat for problems with solution check_SAT([A & B, A & C]) check_SAT([A | B, P & Q, P & B]) From 809988d70df63affcfb8df55ba079cf298534f5f Mon Sep 17 00:00:00 2001 From: tianqiyang Date: Mon, 5 Aug 2019 13:04:38 -0400 Subject: [PATCH 2/3] add chapter 18 and 19 for 4th edition (#1076) * chapter 18 learning * add chapter 19 * move init dataset in NN learner * add adam optimizer, add nn learner * remove cpt 19 for debug * change while loop in games4e * add chapter 19 * add sgd and adam optimizer * add chpt19 deep nn * add rnn * add auto encoder * add comments, correct tests * add more comments, change algorithms according to orders of chapter sections * add keras and numpy to requirements * add tf as requirement * add gc in test agent * fix agent bugs for running test_agent and test_agent_4e together * fix build error * add chapter 21 and 22 * add chapter 12 and part of 13 * remove chapter 12 and 13, add test of rl * modify rnn test * fix build error * Update utils4e.py --- DeepNeuralNet4e.py | 505 ++++++++++++++++++++++++ agents_4e.py | 2 +- games4e.py | 4 +- learning4e.py | 834 +++++++++++++++++++++++++++++++++++++++ nlp4e.py | 523 ++++++++++++++++++++++++ requirements.txt | 2 +- rl4e.py | 340 ++++++++++++++++ tests/test_agents.py | 20 +- tests/test_deepNN.py | 74 ++++ tests/test_learning4e.py | 103 +++++ tests/test_nlp4e.py | 135 +++++++ tests/test_rl4e.py | 66 ++++ utils4e.py | 6 + 13 files changed, 2600 insertions(+), 14 deletions(-) create mode 100644 DeepNeuralNet4e.py create mode 100644 learning4e.py create mode 100644 nlp4e.py create mode 100644 rl4e.py create mode 100644 tests/test_deepNN.py create mode 100644 tests/test_learning4e.py create mode 100644 tests/test_nlp4e.py create mode 100644 tests/test_rl4e.py diff --git a/DeepNeuralNet4e.py b/DeepNeuralNet4e.py new file mode 100644 index 000000000..a353df95c --- /dev/null +++ b/DeepNeuralNet4e.py @@ -0,0 +1,505 @@ +import math +import statistics +from utils4e import sigmoid, dotproduct, softmax1D, conv1D, GaussianKernel, element_wise_product, \ + vector_add, random_weights, scalar_vector_product, matrix_multiplication, map_vector +import random + +from keras import optimizers +from keras.models import Sequential +from keras.layers import Dense, SimpleRNN +from keras.layers.embeddings import Embedding +from keras.preprocessing import sequence + +# DEEP NEURAL NETWORKS. (Chapter 19) +# ________________________________________________ +# 19.2 Common Loss Functions + + +def cross_entropy_loss(X, Y): + """Example of cross entropy loss. X and Y are 1D iterable objects""" + n = len(X) + return (-1.0/n)*sum(x*math.log(y) + (1-x)*math.log(1-y) for x, y in zip(X, Y)) + + +def mse_loss(X, Y): + """Example of min square loss. X and Y are 1D iterable objects""" + n = len(X) + return (1.0/n)*sum((x-y)**2 for x, y in zip(X, Y)) + +# ________________________________________________ +# 19.3 Models +# 19.3.1 Computational Graphs and Layers + + +class Node: + """ + A node in computational graph, It contains the pointer to all its parents. + :param val: value of current node. + :param parents: a container of all parents of current node. + """ + + def __init__(self, val=None, parents=[]): + self.val = val + self.parents = parents + + def __repr__(self): + return "".format(self.val) + + +class NNUnit(Node): + """ + A single unit of a Layer in a Neural Network + :param weights: weights between parent nodes and current node + :param value: value of current node + """ + + def __init__(self, weights=None, value=None): + super(NNUnit, self).__init__(value) + self.weights = weights or [] + + +class Layer: + """ + A layer in a neural network based on computational graph. + :param size: number of units in the current layer + """ + + def __init__(self, size=3): + self.nodes = [NNUnit() for _ in range(size)] + + def forward(self, inputs): + """Define the operation to get the output of this layer""" + raise NotImplementedError + + +# 19.3.2 Output Layers + + +class OutputLayer(Layer): + """Example of a 1D softmax output layer in 19.3.2""" + def __init__(self, size=3): + super(OutputLayer, self).__init__(size) + + def forward(self, inputs): + assert len(self.nodes) == len(inputs) + res = softmax1D(inputs) + for node, val in zip(self.nodes, res): + node.val = val + return res + + +class InputLayer(Layer): + """Example of a 1D input layer. Layer size is the same as input vector size.""" + def __init__(self, size=3): + super(InputLayer, self).__init__(size) + + def forward(self, inputs): + """Take each value of the inputs to each unit in the layer.""" + assert len(self.nodes) == len(inputs) + for node, inp in zip(self.nodes, inputs): + node.val = inp + return inputs + +# 19.3.3 Hidden Layers + + +class DenseLayer(Layer): + """ + 1D dense layer in a neural network. + :param in_size: input vector size, int. + :param out_size: output vector size, int. + :param activation: activation function, Activation object. + """ + + def __init__(self, in_size=3, out_size=3, activation=None): + super(DenseLayer, self).__init__(out_size) + self.out_size = out_size + self.inputs = None + self.activation = sigmoid() if not activation else activation + # initialize weights + for node in self.nodes: + node.weights = random_weights(-0.5, 0.5, in_size) + + def forward(self, inputs): + self.inputs = inputs + res = [] + # get the output value of each unit + for unit in self.nodes: + val = self.activation.f(dotproduct(unit.weights, inputs)) + unit.val = val + res.append(val) + return res + +# 19.3.4 Convolutional networks + + +class ConvLayer1D(Layer): + """ + 1D convolution layer of in neural network. + :param kernel_size: convolution kernel size + """ + + def __init__(self, size=3, kernel_size=3): + super(ConvLayer1D, self).__init__(size) + # init convolution kernel as gaussian kernel + for node in self.nodes: + node.weights = GaussianKernel(kernel_size) + + def forward(self, features): + # Each node in layer takes a channel in the features. + assert len(self.nodes) == len(features) + res = [] + # compute the convolution output of each channel, store it in node.val. + for node, feature in zip(self.nodes, features): + out = conv1D(feature, node.weights) + res.append(out) + node.val = out + return res + +# 19.3.5 Pooling and Downsampling + + +class MaxPoolingLayer1D(Layer): + """1D max pooling layer in a neural network. + :param kernel_size: max pooling area size""" + + def __init__(self, size=3, kernel_size=3): + super(MaxPoolingLayer1D, self).__init__(size) + self.kernel_size = kernel_size + self.inputs = None + + def forward(self, features): + assert len(self.nodes) == len(features) + res = [] + self.inputs = features + # do max pooling for each channel in features + for i in range(len(self.nodes)): + feature = features[i] + # get the max value in a kernel_size * kernel_size area + out = [max(feature[i:i+self.kernel_size]) for i in range(len(feature)-self.kernel_size+1)] + res.append(out) + self.nodes[i].val = out + return res + +# ____________________________________________________________________ +# 19.4 optimization algorithms + + +def init_examples(examples, idx_i, idx_t, o_units): + """Init examples from dataset.examples.""" + + inputs, targets = {}, {} + # random.shuffle(examples) + for i, e in enumerate(examples): + # Input values of e + inputs[i] = [e[i] for i in idx_i] + + if o_units > 1: + # One-Hot representation of e's target + t = [0 for i in range(o_units)] + t[e[idx_t]] = 1 + targets[i] = t + else: + # Target value of e + targets[i] = [e[idx_t]] + + return inputs, targets + +# 19.4.1 Stochastic gradient descent + + +def gradient_descent(dataset, net, loss, epochs=1000, l_rate=0.01, batch_size=1): + """ + gradient descent algorithm to update the learnable parameters of a network. + :return: the updated network. + """ + # init data + examples = dataset.examples + + for e in range(epochs): + total_loss = 0 + random.shuffle(examples) + weights = [[node.weights for node in layer.nodes] for layer in net] + + for batch in get_batch(examples, batch_size): + + inputs, targets = init_examples(batch, dataset.inputs, dataset.target, len(net[-1].nodes)) + # compute gradients of weights + gs, batch_loss = BackPropagation(inputs, targets, weights, net, loss) + # update weights with gradient descent + weights = vector_add(weights, scalar_vector_product(-l_rate, gs)) + total_loss += batch_loss + # update the weights of network each batch + for i in range(len(net)): + if weights[i]: + for j in range(len(weights[i])): + net[i].nodes[j].weights = weights[i][j] + + if (e+1) % 10 == 0: + print("epoch:{}, total_loss:{}".format(e+1,total_loss)) + return net + + +# 19.4.2 Other gradient-based optimization algorithms + + +def adam_optimizer(dataset, net, loss, epochs=1000, rho=(0.9, 0.999), delta=1/10**8, l_rate=0.001, batch_size=1): + """ + Adam optimizer in Figure 19.6 to update the learnable parameters of a network. + Required parameters are similar to gradient descent. + :return the updated network + """ + examples = dataset.examples + + # init s,r and t + s = [[[0] * len(node.weights) for node in layer.nodes] for layer in net] + r = [[[0] * len(node.weights) for node in layer.nodes] for layer in net] + t = 0 + + # repeat util converge + for e in range(epochs): + # total loss of each epoch + total_loss = 0 + random.shuffle(examples) + weights = [[node.weights for node in layer.nodes] for layer in net] + + for batch in get_batch(examples, batch_size): + t += 1 + inputs, targets = init_examples(batch, dataset.inputs, dataset.target, len(net[-1].nodes)) + # compute gradients of weights + gs, batch_loss = BackPropagation(inputs, targets, weights, net, loss) + # update s,r,s_hat and r_gat + s = vector_add(scalar_vector_product(rho[0], s), + scalar_vector_product((1 - rho[0]), gs)) + r = vector_add(scalar_vector_product(rho[1], r), + scalar_vector_product((1 - rho[1]), element_wise_product(gs, gs))) + s_hat = scalar_vector_product(1 / (1 - rho[0] ** t), s) + r_hat = scalar_vector_product(1 / (1 - rho[1] ** t), r) + # rescale r_hat + r_hat = map_vector(lambda x: 1/(math.sqrt(x)+delta), r_hat) + # delta weights + delta_theta = scalar_vector_product(-l_rate, element_wise_product(s_hat, r_hat)) + weights = vector_add(weights, delta_theta) + total_loss += batch_loss + # update the weights of network each batch + for i in range(len(net)): + if weights[i]: + for j in range(len(weights[i])): + net[i].nodes[j].weights = weights[i][j] + + if (e+1) % 10 == 0: + print("epoch:{}, total_loss:{}".format(e+1,total_loss)) + return net + +# 19.4.3 Back-propagation + + +def BackPropagation(inputs, targets, theta, net, loss): + """ + The back-propagation algorithm for multilayer networks in only one epoch, to calculate gradients of theta + :param inputs: A batch of inputs in an array. Each input is an iterable object. + :param targets: A batch of targets in an array. Each target is an iterable object. + :param theta: parameters to be updated. + :param net: a list of predefined layer objects representing their linear sequence. + :param loss: a predefined loss function taking array of inputs and targets. + :return: gradients of theta, loss of the input batch. + """ + + assert len(inputs) == len(targets) + o_units = len(net[-1].nodes) + n_layers = len(net) + batch_size = len(inputs) + + gradients = [[[] for _ in layer.nodes] for layer in net] + total_gradients = [[[0]*len(node.weights) for node in layer.nodes] for layer in net] + + batch_loss = 0 + + # iterate over each example in batch + for e in range(batch_size): + i_val = inputs[e] + t_val = targets[e] + + # Forward pass and compute batch loss + for i in range(1, n_layers): + layer_out = net[i].forward(i_val) + i_val = layer_out + batch_loss += loss(t_val, layer_out) + + # Initialize delta + delta = [[] for _ in range(n_layers)] + + previous = [layer_out[i]-t_val[i] for i in range(o_units)] + h_layers = n_layers - 1 + # Backward pass + for i in range(h_layers, 0, -1): + layer = net[i] + derivative = [layer.activation.derivative(node.val) for node in layer.nodes] + delta[i] = element_wise_product(previous, derivative) + # pass to layer i-1 in the next iteration + previous = matrix_multiplication([delta[i]], theta[i])[0] + # compute gradient of layer i + gradients[i] = [scalar_vector_product(d, net[i].inputs) for d in delta[i]] + + # add gradient of current example to batch gradient + total_gradients = vector_add(total_gradients, gradients) + + return total_gradients, batch_loss + +# 19.4.5 Batch normalization + + +class BatchNormalizationLayer(Layer): + """Example of a batch normalization layer.""" + def __init__(self, size, epsilon=0.001): + super(BatchNormalizationLayer, self).__init__(size) + self.epsilon = epsilon + # self.weights = [beta, gamma] + self.weights = [0, 0] + self.inputs = None + + def forward(self, inputs): + # mean value of inputs + mu = sum(inputs) / len(inputs) + # standard error of inputs + stderr = statistics.stdev(inputs) + self.inputs = inputs + res = [] + # get normalized value of each input + for i in range(len(self.nodes)): + val = [(inputs[i] - mu)*self.weights[0]/math.sqrt(self.epsilon + stderr**2)+self.weights[1]] + res.append(val) + self.nodes[i].val = val + return res + + +def get_batch(examples, batch_size=1): + """split examples into multiple batches""" + for i in range(0, len(examples), batch_size): + yield examples[i: i+batch_size] + +# example of NNs + + +def neural_net_learner(dataset, hidden_layer_sizes=[4], learning_rate=0.01, epochs=100, optimizer=gradient_descent, batch_size=1): + """Example of a simple dense multilayer neural network. + :param hidden_layer_sizes: size of hidden layers in the form of a list""" + + input_size = len(dataset.inputs) + output_size = len(dataset.values[dataset.target]) + + # initialize the network + raw_net = [InputLayer(input_size)] + # add hidden layers + hidden_input_size = input_size + for h_size in hidden_layer_sizes: + raw_net.append(DenseLayer(hidden_input_size, h_size)) + hidden_input_size = h_size + raw_net.append(DenseLayer(hidden_input_size, output_size)) + + # update parameters of the network + learned_net = optimizer(dataset, raw_net, mse_loss, epochs, l_rate=learning_rate, batch_size=batch_size) + + def predict(example): + n_layers = len(learned_net) + + layer_input = example + layer_out = example + + # get the output of each layer by forward passing + for i in range(1, n_layers): + layer_out = learned_net[i].forward(layer_input) + layer_input = layer_out + + return layer_out.index(max(layer_out)) + + return predict + + +def perceptron_learner(dataset, learning_rate=0.01, epochs=100): + """ + Example of a simple perceptron neural network. + """ + input_size = len(dataset.inputs) + output_size = len(dataset.values[dataset.target]) + + # initialize the network, add dense layer + raw_net = [InputLayer(input_size), DenseLayer(input_size, output_size)] + # update the network + learned_net = gradient_descent(dataset, raw_net, mse_loss, epochs, l_rate=learning_rate) + + def predict(example): + + layer_out = learned_net[1].forward(example) + return layer_out.index(max(layer_out)) + + return predict + +# ____________________________________________________________________ +# 19.6 Recurrent neural networks + + +def simple_rnn_learner(train_data, val_data, epochs=2): + """ + rnn example for text sentimental analysis + :param train_data: a tuple of (training data, targets) + Training data: ndarray taking training examples, while each example is coded by embedding + Targets: ndarry taking targets of each example. Each target is mapped to an integer. + :param val_data: a tuple of (validation data, targets) + :return: a keras model + """ + + total_inputs = 5000 + input_length = 500 + + # init data + X_train, y_train = train_data + X_val, y_val = val_data + + # init a the sequential network (embedding layer, rnn layer, dense layer) + model = Sequential() + model.add(Embedding(total_inputs, 32, input_length=input_length)) + model.add(SimpleRNN(units=128)) + model.add(Dense(1, activation='sigmoid')) + model.compile(loss='binary_crossentropy', optimizer='adam', metrics=['accuracy']) + + # train the model + model.fit(X_train, y_train, validation_data=(X_val, y_val), epochs=epochs, batch_size=128, verbose=2) + + return model + + +def keras_dataset_loader(dataset, max_length=500): + """ + helper function to load keras datasets + :param dataset: keras data set type + :param max_length: max length of each input sequence + """ + # init dataset + (X_train, y_train), (X_val, y_val) = dataset + if max_length > 0: + X_train = sequence.pad_sequences(X_train, maxlen=max_length) + X_val = sequence.pad_sequences(X_val, maxlen=max_length) + return (X_train[10:], y_train[10:]), (X_val, y_val), (X_train[:10], y_train[:10]) + + +def auto_encoder_learner(inputs, encoding_size, epochs=200): + """simple example of linear auto encoder learning producing the input itself. + :param inputs: a batch of input data in np.ndarray type + :param encoding_size: int, the size of encoding layer""" + + # init data + input_size = len(inputs[0]) + + # init model + model = Sequential() + model.add(Dense(encoding_size, input_dim=input_size, activation='relu', kernel_initializer='random_uniform',bias_initializer='ones')) + model.add(Dense(input_size, activation='relu', kernel_initializer='random_uniform', bias_initializer='ones')) + # update model with sgd + sgd = optimizers.SGD(lr=0.01) + model.compile(loss='mean_squared_error', optimizer=sgd, metrics=['accuracy']) + + # train the model + model.fit(inputs, inputs, epochs=epochs, batch_size=10, verbose=2) + + return model diff --git a/agents_4e.py b/agents_4e.py index 606e3e25a..3734ee91d 100644 --- a/agents_4e.py +++ b/agents_4e.py @@ -514,7 +514,7 @@ def add_thing(self, thing, location=(1, 1), exclude_duplicate_class_items=False) def is_inbounds(self, location): """Checks to make sure that the location is inbounds (within walls if we have walls)""" x, y = location - return not (x < self.x_start or x >= self.x_end or y < self.y_start or y >= self.y_end) + return not (x < self.x_start or x > self.x_end or y < self.y_start or y > self.y_end) def random_location_inbounds(self, exclude=None): """Returns a random location that is inbounds (within walls if we have walls)""" diff --git a/games4e.py b/games4e.py index f32259175..84e082c1a 100644 --- a/games4e.py +++ b/games4e.py @@ -210,12 +210,12 @@ def backprop(n, utility): root = MCT_Node(state=state) - while N > 0: + for _ in range(N): leaf = select(root) child = expand(leaf) result = simulate(game, child.state) backprop(child, result) - N -= 1 + max_state = max(root.children, key=lambda p: p.N) return root.children.get(max_state) diff --git a/learning4e.py b/learning4e.py new file mode 100644 index 000000000..68a2d5c48 --- /dev/null +++ b/learning4e.py @@ -0,0 +1,834 @@ +from utils4e import ( + removeall, unique, mode, argmax_random_tie, isclose, dotproduct, weighted_sample_with_replacement, + num_or_str, normalize, clip, print_table, open_data, probability, random_weights +) + +import copy +import heapq +import math +import random + +from statistics import mean, stdev +from collections import defaultdict + +# Learn to estimate functions from examples. (Chapters 18) +# ______________________________________________________________________________ +# 18.2 Supervised learning. +# define supervised learning dataset and utility functions/ + + +def mean_boolean_error(X, Y): + return mean(int(x != y) for x, y in zip(X, Y)) + + +class DataSet: + """A data set for a machine learning problem. It has the following fields: + + d.examples A list of examples. Each one is a list of attribute values. + d.attrs A list of integers to index into an example, so example[attr] + gives a value. Normally the same as range(len(d.examples[0])). + d.attrnames Optional list of mnemonic names for corresponding attrs. + d.target The attribute that a learning algorithm will try to predict. + By default the final attribute. + d.inputs The list of attrs without the target. + d.values A list of lists: each sublist is the set of possible + values for the corresponding attribute. If initially None, + it is computed from the known examples by self.setproblem. + If not None, an erroneous value raises ValueError. + d.distance A function from a pair of examples to a nonnegative number. + Should be symmetric, etc. Defaults to mean_boolean_error + since that can handle any field types. + d.name Name of the data set (for output display only). + d.source URL or other source where the data came from. + d.exclude A list of attribute indexes to exclude from d.inputs. Elements + of this list can either be integers (attrs) or attrnames. + + Normally, you call the constructor and you're done; then you just + access fields like d.examples and d.target and d.inputs.""" + + def __init__(self, examples=None, attrs=None, attrnames=None, target=-1, + inputs=None, values=None, distance=mean_boolean_error, + name='', source='', exclude=()): + """Accepts any of DataSet's fields. Examples can also be a + string or file from which to parse examples using parse_csv. + Optional parameter: exclude, as documented in .setproblem(). + >>> DataSet(examples='1, 2, 3') + + """ + self.name = name + self.source = source + self.values = values + self.distance = distance + self.got_values_flag = bool(values) + + # Initialize .examples from string or list or data directory + if isinstance(examples, str): + self.examples = parse_csv(examples) + elif examples is None: + self.examples = parse_csv(open_data(name + '.csv').read()) + else: + self.examples = examples + + # Attrs are the indices of examples, unless otherwise stated. + if self.examples is not None and attrs is None: + attrs = list(range(len(self.examples[0]))) + + self.attrs = attrs + + # Initialize .attrnames from string, list, or by default + if isinstance(attrnames, str): + self.attrnames = attrnames.split() + else: + self.attrnames = attrnames or attrs + self.setproblem(target, inputs=inputs, exclude=exclude) + + def setproblem(self, target, inputs=None, exclude=()): + """Set (or change) the target and/or inputs. + This way, one DataSet can be used multiple ways. inputs, if specified, + is a list of attributes, or specify exclude as a list of attributes + to not use in inputs. Attributes can be -n .. n, or an attrname. + Also computes the list of possible values, if that wasn't done yet.""" + self.target = self.attrnum(target) + exclude = list(map(self.attrnum, exclude)) + if inputs: + self.inputs = removeall(self.target, inputs) + else: + self.inputs = [a for a in self.attrs + if a != self.target and a not in exclude] + if not self.values: + self.update_values() + self.check_me() + + def check_me(self): + """Check that my fields make sense.""" + assert len(self.attrnames) == len(self.attrs) + assert self.target in self.attrs + assert self.target not in self.inputs + assert set(self.inputs).issubset(set(self.attrs)) + if self.got_values_flag: + # only check if values are provided while initializing DataSet + list(map(self.check_example, self.examples)) + + def add_example(self, example): + """Add an example to the list of examples, checking it first.""" + self.check_example(example) + self.examples.append(example) + + def check_example(self, example): + """Raise ValueError if example has any invalid values.""" + if self.values: + for a in self.attrs: + if example[a] not in self.values[a]: + raise ValueError('Bad value {} for attribute {} in {}' + .format(example[a], self.attrnames[a], example)) + + def attrnum(self, attr): + """Returns the number used for attr, which can be a name, or -n .. n-1.""" + if isinstance(attr, str): + return self.attrnames.index(attr) + elif attr < 0: + return len(self.attrs) + attr + else: + return attr + + def update_values(self): + self.values = list(map(unique, zip(*self.examples))) + + def sanitize(self, example): + """Return a copy of example, with non-input attributes replaced by None.""" + return [attr_i if i in self.inputs else None + for i, attr_i in enumerate(example)] + + def classes_to_numbers(self, classes=None): + """Converts class names to numbers.""" + if not classes: + # If classes were not given, extract them from values + classes = sorted(self.values[self.target]) + for item in self.examples: + item[self.target] = classes.index(item[self.target]) + + def remove_examples(self, value=''): + """Remove examples that contain given value.""" + self.examples = [x for x in self.examples if value not in x] + self.update_values() + + def split_values_by_classes(self): + """Split values into buckets according to their class.""" + buckets = defaultdict(lambda: []) + target_names = self.values[self.target] + + for v in self.examples: + item = [a for a in v if a not in target_names] # Remove target from item + buckets[v[self.target]].append(item) # Add item to bucket of its class + + return buckets + + def find_means_and_deviations(self): + """Finds the means and standard deviations of self.dataset. + means : A dictionary for each class/target. Holds a list of the means + of the features for the class. + deviations: A dictionary for each class/target. Holds a list of the sample + standard deviations of the features for the class.""" + target_names = self.values[self.target] + feature_numbers = len(self.inputs) + + item_buckets = self.split_values_by_classes() + + means = defaultdict(lambda: [0] * feature_numbers) + deviations = defaultdict(lambda: [0] * feature_numbers) + + for t in target_names: + # Find all the item feature values for item in class t + features = [[] for i in range(feature_numbers)] + for item in item_buckets[t]: + for i in range(feature_numbers): + features[i].append(item[i]) + + # Calculate means and deviations fo the class + for i in range(feature_numbers): + means[t][i] = mean(features[i]) + deviations[t][i] = stdev(features[i]) + + return means, deviations + + def __repr__(self): + return ''.format( + self.name, len(self.examples), len(self.attrs)) + +# ______________________________________________________________________________ + + +def parse_csv(input, delim=','): + r"""Input is a string consisting of lines, each line has comma-delimited + fields. Convert this into a list of lists. Blank lines are skipped. + Fields that look like numbers are converted to numbers. + The delim defaults to ',' but '\t' and None are also reasonable values. + >>> parse_csv('1, 2, 3 \n 0, 2, na') + [[1, 2, 3], [0, 2, 'na']]""" + lines = [line for line in input.splitlines() if line.strip()] + return [list(map(num_or_str, line.split(delim))) for line in lines] + +# ______________________________________________________________________________ +# 18.3 Learning decision trees + + +class DecisionFork: + """A fork of a decision tree holds an attribute to test, and a dict + of branches, one for each of the attribute's values.""" + + def __init__(self, attr, attrname=None, default_child=None, branches=None): + """Initialize by saying what attribute this node tests.""" + self.attr = attr + self.attrname = attrname or attr + self.default_child = default_child + self.branches = branches or {} + + def __call__(self, example): + """Given an example, classify it using the attribute and the branches.""" + attrvalue = example[self.attr] + if attrvalue in self.branches: + return self.branches[attrvalue](example) + else: + # return default class when attribute is unknown + return self.default_child(example) + + def add(self, val, subtree): + """Add a branch. If self.attr = val, go to the given subtree.""" + self.branches[val] = subtree + + def display(self, indent=0): + name = self.attrname + print('Test', name) + for (val, subtree) in self.branches.items(): + print(' ' * 4 * indent, name, '=', val, '==>', end=' ') + subtree.display(indent + 1) + print() # newline + + def __repr__(self): + return ('DecisionFork({0!r}, {1!r}, {2!r})' + .format(self.attr, self.attrname, self.branches)) + + +class DecisionLeaf: + """A leaf of a decision tree holds just a result.""" + + def __init__(self, result): + self.result = result + + def __call__(self, example): + return self.result + + def display(self, indent=0): + print('RESULT =', self.result) + + def __repr__(self): + return repr(self.result) + +# decision tree learning in Figure 18.5 + + +def DecisionTreeLearner(dataset): + + target, values = dataset.target, dataset.values + + def decision_tree_learning(examples, attrs, parent_examples=()): + if len(examples) == 0: + return plurality_value(parent_examples) + elif all_same_class(examples): + return DecisionLeaf(examples[0][target]) + elif len(attrs) == 0: + return plurality_value(examples) + else: + A = choose_attribute(attrs, examples) + tree = DecisionFork(A, dataset.attrnames[A], plurality_value(examples)) + for (v_k, exs) in split_by(A, examples): + subtree = decision_tree_learning( + exs, removeall(A, attrs), examples) + tree.add(v_k, subtree) + return tree + + def plurality_value(examples): + """Return the most popular target value for this set of examples. + (If target is binary, this is the majority; otherwise plurality.)""" + popular = argmax_random_tie(values[target], + key=lambda v: count(target, v, examples)) + return DecisionLeaf(popular) + + def count(attr, val, examples): + """Count the number of examples that have example[attr] = val.""" + return sum(e[attr] == val for e in examples) + + def all_same_class(examples): + """Are all these examples in the same target class?""" + class0 = examples[0][target] + return all(e[target] == class0 for e in examples) + + def choose_attribute(attrs, examples): + """Choose the attribute with the highest information gain.""" + return argmax_random_tie(attrs, + key=lambda a: information_gain(a, examples)) + + def information_gain(attr, examples): + """Return the expected reduction in entropy from splitting by attr.""" + def I(examples): + return information_content([count(target, v, examples) + for v in values[target]]) + N = len(examples) + remainder = sum((len(examples_i)/N) * I(examples_i) + for (v, examples_i) in split_by(attr, examples)) + return I(examples) - remainder + + def split_by(attr, examples): + """Return a list of (val, examples) pairs for each val of attr.""" + return [(v, [e for e in examples if e[attr] == v]) + for v in values[attr]] + + return decision_tree_learning(dataset.examples, dataset.inputs) + + +def information_content(values): + """Number of bits to represent the probability distribution in values.""" + probabilities = normalize(removeall(0, values)) + return sum(-p * math.log2(p) for p in probabilities) + +# ______________________________________________________________________________ +# 18.4 Model selection and optimization + + +def model_selection(learner, dataset, k=10, trials=1): + """[Fig 18.8] + Return the optimal value of size having minimum error + on validation set. + err_train: A training error array, indexed by size + err_val: A validation error array, indexed by size + """ + errs = [] + size = 1 + + while True: + err = cross_validation(learner, size, dataset, k, trials) + # Check for convergence provided err_val is not empty + if err and not isclose(err[-1], err, rel_tol=1e-6): + best_size = 0 + min_val = math.inf + + i = 0 + while i < size: + if errs[i] < min_val: + min_val = errs[i] + best_size = i + i += 1 + return learner(dataset, best_size) + errs.append(err) + size += 1 + + +def cross_validation(learner, size, dataset, k=10, trials=1): + """Do k-fold cross_validate and return their mean. + That is, keep out 1/k of the examples for testing on each of k runs. + Shuffle the examples first; if trials>1, average over several shuffles. + Returns Training error, Validataion error""" + k = k or len(dataset.examples) + if trials > 1: + trial_errs = 0 + for t in range(trials): + errs = cross_validation(learner, size, dataset, + k=10, trials=1) + trial_errs += errs + return trial_errs/trials + else: + fold_errs = 0 + n = len(dataset.examples) + examples = dataset.examples + random.shuffle(dataset.examples) + for fold in range(k): + train_data, val_data = train_test_split(dataset, fold * (n / k), + (fold + 1) * (n / k)) + dataset.examples = train_data + h = learner(dataset, size) + fold_errs += err_ratio(h, dataset, train_data) + + # Reverting back to original once test is completed + dataset.examples = examples + return fold_errs/k + + +def err_ratio(predict, dataset, examples=None, verbose=0): + """Return the proportion of the examples that are NOT correctly predicted. + verbose - 0: No output; 1: Output wrong; 2 (or greater): Output correct""" + examples = examples or dataset.examples + if len(examples) == 0: + return 0.0 + right = 0 + for example in examples: + desired = example[dataset.target] + output = predict(dataset.sanitize(example)) + if output == desired: + right += 1 + if verbose >= 2: + print(' OK: got {} for {}'.format(desired, example)) + elif verbose: + print('WRONG: got {}, expected {} for {}'.format( + output, desired, example)) + return 1 - (right/len(examples)) + + +def train_test_split(dataset, start=None, end=None, test_split=None): + """If you are giving 'start' and 'end' as parameters, + then it will return the testing set from index 'start' to 'end' + and the rest for training. + If you give 'test_split' as a parameter then it will return + test_split * 100% as the testing set and the rest as + training set. + """ + examples = dataset.examples + if test_split == None: + train = examples[:start] + examples[end:] + val = examples[start:end] + else: + total_size = len(examples) + val_size = int(total_size * test_split) + train_size = total_size - val_size + train = examples[:train_size] + val = examples[train_size:total_size] + + return train, val + + +def grade_learner(predict, tests): + """Grades the given learner based on how many tests it passes. + tests is a list with each element in the form: (values, output).""" + return mean(int(predict(X) == y) for X, y in tests) + + +def leave_one_out(learner, dataset, size=None): + """Leave one out cross-validation over the dataset.""" + return cross_validation(learner, size, dataset, k=len(dataset.examples)) + + +# TODO learningcurve needs to fixed +def learningcurve(learner, dataset, trials=10, sizes=None): + if sizes is None: + sizes = list(range(2, len(dataset.examples) - 10, 2)) + + def score(learner, size): + random.shuffle(dataset.examples) + return train_test_split(learner, dataset, 0, size) + return [(size, mean([score(learner, size) for t in range(trials)])) + for size in sizes] + +# ______________________________________________________________________________ +# 18.5 The theory Of learning + + +def DecisionListLearner(dataset): + """A decision list is implemented as a list of (test, value) pairs.[Figure 18.11]""" + + # TODO: where are the tests from? + def decision_list_learning(examples): + if not examples: + return [(True, False)] + t, o, examples_t = find_examples(examples) + if not t: + raise Exception + return [(t, o)] + decision_list_learning(examples - examples_t) + + def find_examples(examples): + """Find a set of examples that all have the same outcome under + some test. Return a tuple of the test, outcome, and examples.""" + raise NotImplementedError + + def passes(example, test): + """Does the example pass the test?""" + return test.test(example) + raise NotImplementedError + + def predict(example): + """Predict the outcome for the first passing test.""" + for test, outcome in predict.decision_list: + if passes(example, test): + return outcome + + predict.decision_list = decision_list_learning(set(dataset.examples)) + + return predict + +# ______________________________________________________________________________ +# 18.6 Linear regression and classification + + +def LinearLearner(dataset, learning_rate=0.01, epochs=100): + """Define with learner = LinearLearner(data); infer with learner(x).""" + idx_i = dataset.inputs + idx_t = dataset.target # As of now, dataset.target gives only one index. + examples = dataset.examples + num_examples = len(examples) + + # X transpose + X_col = [dataset.values[i] for i in idx_i] # vertical columns of X + + # Add dummy + ones = [1 for _ in range(len(examples))] + X_col = [ones] + X_col + + # Initialize random weigts + num_weights = len(idx_i) + 1 + w = random_weights(min_value=-0.5, max_value=0.5, num_weights=num_weights) + + for epoch in range(epochs): + err = [] + # Pass over all examples + for example in examples: + x = [1] + example + y = dotproduct(w, x) + t = example[idx_t] + err.append(t - y) + + # update weights + for i in range(len(w)): + w[i] = w[i] + learning_rate * (dotproduct(err, X_col[i]) / num_examples) + + def predict(example): + x = [1] + example + return dotproduct(w, x) + return predict + + +def LogisticLinearLeaner(dataset, learning_rate=0.01, epochs=100): + """Define logistic regression classifier in 18.6.5""" + idx_i = dataset.inputs + idx_t = dataset.target + examples = dataset.examples + num_examples = len(examples) + + # X transpose + X_col = [dataset.values[i] for i in idx_i] # vertical columns of X + + # Add dummy + ones = [1 for _ in range(len(examples))] + X_col = [ones] + X_col + + # Initialize random weigts + num_weights = len(idx_i) + 1 + w = random_weights(min_value=-0.5, max_value=0.5, num_weights=num_weights) + + for epoch in range(epochs): + err = [] + # Pass over all examples + for example in examples: + x = [1] + example + y = 1/(1 + math.exp(-dotproduct(w, x))) + h = [y * (1-y)] + t = example[idx_t] + err.append(t - y) + + # update weights + for i in range(len(w)): + w[i] = w[i] + learning_rate * (dotproduct(dotproduct(err,h), X_col[i]) / num_examples) + + def predict(example): + x = [1] + example + return 1/(1 + math.exp(-dotproduct(w, x))) + + return predict + +# ______________________________________________________________________________ +# 18.7 Nonparametric models + + +def NearestNeighborLearner(dataset, k=1): + """k-NearestNeighbor: the k nearest neighbors vote.""" + def predict(example): + """Find the k closest items, and have them vote for the best.""" + best = heapq.nsmallest(k, ((dataset.distance(e, example), e) + for e in dataset.examples)) + return mode(e[dataset.target] for (d, e) in best) + return predict + + +# ______________________________________________________________________________ +# 18.8 Ensemble learning + + +def EnsembleLearner(learners): + """Given a list of learning algorithms, have them vote.""" + def train(dataset): + predictors = [learner(dataset) for learner in learners] + + def predict(example): + return mode(predictor(example) for predictor in predictors) + return predict + return train + + +def RandomForest(dataset, n=5): + """An ensemble of Decision Trees trained using bagging and feature bagging.""" + + def data_bagging(dataset, m=0): + """Sample m examples with replacement""" + n = len(dataset.examples) + return weighted_sample_with_replacement(m or n, dataset.examples, [1]*n) + + def feature_bagging(dataset, p=0.7): + """Feature bagging with probability p to retain an attribute""" + inputs = [i for i in dataset.inputs if probability(p)] + return inputs or dataset.inputs + + def predict(example): + print([predictor(example) for predictor in predictors]) + return mode(predictor(example) for predictor in predictors) + + predictors = [DecisionTreeLearner(DataSet(examples=data_bagging(dataset), + attrs=dataset.attrs, + attrnames=dataset.attrnames, + target=dataset.target, + inputs=feature_bagging(dataset))) for _ in range(n)] + + return predict + + +def AdaBoost(L, K): + """[Figure 18.34]""" + + def train(dataset): + examples, target = dataset.examples, dataset.target + N = len(examples) + epsilon = 1/(2*N) + w = [1/N]*N + h, z = [], [] + for k in range(K): + h_k = L(dataset, w) + h.append(h_k) + error = sum(weight for example, weight in zip(examples, w) + if example[target] != h_k(example)) + + # Avoid divide-by-0 from either 0% or 100% error rates: + error = clip(error, epsilon, 1 - epsilon) + for j, example in enumerate(examples): + if example[target] == h_k(example): + w[j] *= error/(1 - error) + w = normalize(w) + z.append(math.log((1 - error)/error)) + return WeightedMajority(h, z) + return train + + +def WeightedMajority(predictors, weights): + """Return a predictor that takes a weighted vote.""" + def predict(example): + return weighted_mode((predictor(example) for predictor in predictors), + weights) + return predict + + +def weighted_mode(values, weights): + """Return the value with the greatest total weight. + >>> weighted_mode('abbaa', [1, 2, 3, 1, 2]) + 'b' + """ + totals = defaultdict(int) + for v, w in zip(values, weights): + totals[v] += w + return max(totals, key=totals.__getitem__) + +# _____________________________________________________________________________ +# Adapting an unweighted learner for AdaBoost + + +def WeightedLearner(unweighted_learner): + """Given a learner that takes just an unweighted dataset, return + one that takes also a weight for each example. [p. 749 footnote 14]""" + def train(dataset, weights): + return unweighted_learner(replicated_dataset(dataset, weights)) + return train + + +def replicated_dataset(dataset, weights, n=None): + """Copy dataset, replicating each example in proportion to its weight.""" + n = n or len(dataset.examples) + result = copy.copy(dataset) + result.examples = weighted_replicate(dataset.examples, weights, n) + return result + + +def weighted_replicate(seq, weights, n): + """Return n selections from seq, with the count of each element of + seq proportional to the corresponding weight (filling in fractions + randomly). + >>> weighted_replicate('ABC', [1, 2, 1], 4) + ['A', 'B', 'B', 'C'] + """ + assert len(seq) == len(weights) + weights = normalize(weights) + wholes = [int(w*n) for w in weights] + fractions = [(w*n) % 1 for w in weights] + return (flatten([x]*nx for x, nx in zip(seq, wholes)) + + weighted_sample_with_replacement(n - sum(wholes), seq, fractions)) + + +def flatten(seqs): return sum(seqs, []) + +# _____________________________________________________________________________ +# Functions for testing learners on examples +# The rest of this file gives datasets for machine learning problems. + + +orings = DataSet(name='orings', target='Distressed', + attrnames="Rings Distressed Temp Pressure Flightnum") + + +zoo = DataSet(name='zoo', target='type', exclude=['name'], + attrnames="name hair feathers eggs milk airborne aquatic " + + "predator toothed backbone breathes venomous fins legs tail " + + "domestic catsize type") + + +iris = DataSet(name="iris", target="class", + attrnames="sepal-len sepal-width petal-len petal-width class") + +# ______________________________________________________________________________ +# The Restaurant example from [Figure 18.2] + + +def RestaurantDataSet(examples=None): + """Build a DataSet of Restaurant waiting examples. [Figure 18.3]""" + return DataSet(name='restaurant', target='Wait', examples=examples, + attrnames='Alternate Bar Fri/Sat Hungry Patrons Price ' + + 'Raining Reservation Type WaitEstimate Wait') + + +restaurant = RestaurantDataSet() + + +def T(attrname, branches): + branches = {value: (child if isinstance(child, DecisionFork) + else DecisionLeaf(child)) + for value, child in branches.items()} + return DecisionFork(restaurant.attrnum(attrname), attrname, print, branches) + + +""" [Figure 18.2] +A decision tree for deciding whether to wait for a table at a hotel. +""" + +waiting_decision_tree = T('Patrons', + {'None': 'No', 'Some': 'Yes', + 'Full': T('WaitEstimate', + {'>60': 'No', '0-10': 'Yes', + '30-60': T('Alternate', + {'No': T('Reservation', + {'Yes': 'Yes', + 'No': T('Bar', {'No': 'No', + 'Yes': 'Yes'})}), + 'Yes': T('Fri/Sat', {'No': 'No', 'Yes': 'Yes'})} + ), + '10-30': T('Hungry', + {'No': 'Yes', + 'Yes': T('Alternate', + {'No': 'Yes', + 'Yes': T('Raining', + {'No': 'No', + 'Yes': 'Yes'})})})})}) + + +def SyntheticRestaurant(n=20): + """Generate a DataSet with n examples.""" + def gen(): + example = list(map(random.choice, restaurant.values)) + example[restaurant.target] = waiting_decision_tree(example) + return example + return RestaurantDataSet([gen() for i in range(n)]) + +# ______________________________________________________________________________ +# Artificial, generated datasets. + + +def Majority(k, n): + """Return a DataSet with n k-bit examples of the majority problem: + k random bits followed by a 1 if more than half the bits are 1, else 0.""" + examples = [] + for i in range(n): + bits = [random.choice([0, 1]) for i in range(k)] + bits.append(int(sum(bits) > k / 2)) + examples.append(bits) + return DataSet(name="majority", examples=examples) + + +def Parity(k, n, name="parity"): + """Return a DataSet with n k-bit examples of the parity problem: + k random bits followed by a 1 if an odd number of bits are 1, else 0.""" + examples = [] + for i in range(n): + bits = [random.choice([0, 1]) for i in range(k)] + bits.append(sum(bits) % 2) + examples.append(bits) + return DataSet(name=name, examples=examples) + + +def Xor(n): + """Return a DataSet with n examples of 2-input xor.""" + return Parity(2, n, name="xor") + + +def ContinuousXor(n): + "2 inputs are chosen uniformly from (0.0 .. 2.0]; output is xor of ints." + examples = [] + for i in range(n): + x, y = [random.uniform(0.0, 2.0) for i in '12'] + examples.append([x, y, int(x) != int(y)]) + return DataSet(name="continuous xor", examples=examples) + + +def compare(algorithms=None, datasets=None, k=10, trials=1): + """Compare various learners on various datasets using cross-validation. + Print results as a table.""" + algorithms = algorithms or [ # default list + NearestNeighborLearner, DecisionTreeLearner] # of algorithms + + datasets = datasets or [iris, orings, zoo, restaurant, SyntheticRestaurant(20), # default list + Majority(7, 100), Parity(7, 100), Xor(100)] # of datasets + + print_table([[a.__name__.replace('Learner', '')] + + [cross_validation(a, d, k, trials) for d in datasets] + for a in algorithms], + header=[''] + [d.name[0:7] for d in datasets], numfmt='%.2f') diff --git a/nlp4e.py b/nlp4e.py new file mode 100644 index 000000000..98a34e778 --- /dev/null +++ b/nlp4e.py @@ -0,0 +1,523 @@ +"""Natural Language Processing (Chapter 22)""" + +from collections import defaultdict +from utils4e import weighted_choice +import copy +import operator +import heapq +from search import Problem + + +# ______________________________________________________________________________ +# 22.2 Grammars + + +def Rules(**rules): + """Create a dictionary mapping symbols to alternative sequences. + >>> Rules(A = "B C | D E") + {'A': [['B', 'C'], ['D', 'E']]} + """ + for (lhs, rhs) in rules.items(): + rules[lhs] = [alt.strip().split() for alt in rhs.split('|')] + return rules + + +def Lexicon(**rules): + """Create a dictionary mapping symbols to alternative words. + >>> Lexicon(Article = "the | a | an") + {'Article': ['the', 'a', 'an']} + """ + for (lhs, rhs) in rules.items(): + rules[lhs] = [word.strip() for word in rhs.split('|')] + return rules + + +class Grammar: + + def __init__(self, name, rules, lexicon): + """A grammar has a set of rules and a lexicon.""" + self.name = name + self.rules = rules + self.lexicon = lexicon + self.categories = defaultdict(list) + for lhs in lexicon: + for word in lexicon[lhs]: + self.categories[word].append(lhs) + + def rewrites_for(self, cat): + """Return a sequence of possible rhs's that cat can be rewritten as.""" + return self.rules.get(cat, ()) + + def isa(self, word, cat): + """Return True iff word is of category cat""" + return cat in self.categories[word] + + def cnf_rules(self): + """Returns the tuple (X, Y, Z) for rules in the form: + X -> Y Z""" + cnf = [] + for X, rules in self.rules.items(): + for (Y, Z) in rules: + cnf.append((X, Y, Z)) + + return cnf + + def generate_random(self, S='S'): + """Replace each token in S by a random entry in grammar (recursively).""" + import random + + def rewrite(tokens, into): + for token in tokens: + if token in self.rules: + rewrite(random.choice(self.rules[token]), into) + elif token in self.lexicon: + into.append(random.choice(self.lexicon[token])) + else: + into.append(token) + return into + + return ' '.join(rewrite(S.split(), [])) + + def __repr__(self): + return ''.format(self.name) + + +def ProbRules(**rules): + """Create a dictionary mapping symbols to alternative sequences, + with probabilities. + >>> ProbRules(A = "B C [0.3] | D E [0.7]") + {'A': [(['B', 'C'], 0.3), (['D', 'E'], 0.7)]} + """ + for (lhs, rhs) in rules.items(): + rules[lhs] = [] + rhs_separate = [alt.strip().split() for alt in rhs.split('|')] + for r in rhs_separate: + prob = float(r[-1][1:-1]) # remove brackets, convert to float + rhs_rule = (r[:-1], prob) + rules[lhs].append(rhs_rule) + + return rules + + +def ProbLexicon(**rules): + """Create a dictionary mapping symbols to alternative words, + with probabilities. + >>> ProbLexicon(Article = "the [0.5] | a [0.25] | an [0.25]") + {'Article': [('the', 0.5), ('a', 0.25), ('an', 0.25)]} + """ + for (lhs, rhs) in rules.items(): + rules[lhs] = [] + rhs_separate = [word.strip().split() for word in rhs.split('|')] + for r in rhs_separate: + prob = float(r[-1][1:-1]) # remove brackets, convert to float + word = r[:-1][0] + rhs_rule = (word, prob) + rules[lhs].append(rhs_rule) + + return rules + + +class ProbGrammar: + + def __init__(self, name, rules, lexicon): + """A grammar has a set of rules and a lexicon. + Each rule has a probability.""" + self.name = name + self.rules = rules + self.lexicon = lexicon + self.categories = defaultdict(list) + + for lhs in lexicon: + for word, prob in lexicon[lhs]: + self.categories[word].append((lhs, prob)) + + def rewrites_for(self, cat): + """Return a sequence of possible rhs's that cat can be rewritten as.""" + return self.rules.get(cat, ()) + + def isa(self, word, cat): + """Return True iff word is of category cat""" + return cat in [c for c, _ in self.categories[word]] + + def cnf_rules(self): + """Returns the tuple (X, Y, Z, p) for rules in the form: + X -> Y Z [p]""" + cnf = [] + for X, rules in self.rules.items(): + for (Y, Z), p in rules: + cnf.append((X, Y, Z, p)) + + return cnf + + def generate_random(self, S='S'): + """Replace each token in S by a random entry in grammar (recursively). + Returns a tuple of (sentence, probability).""" + + def rewrite(tokens, into): + for token in tokens: + if token in self.rules: + non_terminal, prob = weighted_choice(self.rules[token]) + into[1] *= prob + rewrite(non_terminal, into) + elif token in self.lexicon: + terminal, prob = weighted_choice(self.lexicon[token]) + into[0].append(terminal) + into[1] *= prob + else: + into[0].append(token) + return into + + rewritten_as, prob = rewrite(S.split(), [[], 1]) + return (' '.join(rewritten_as), prob) + + def __repr__(self): + return ''.format(self.name) + + +E0 = Grammar('E0', + Rules( # Grammar for E_0 [Figure 22.2] + S='NP VP | S Conjunction S', + NP='Pronoun | Name | Noun | Article Noun | Digit Digit | NP PP | NP RelClause', + VP='Verb | VP NP | VP Adjective | VP PP | VP Adverb', + PP='Preposition NP', + RelClause='That VP'), + + Lexicon( # Lexicon for E_0 [Figure 22.3] + Noun="stench | breeze | glitter | nothing | wumpus | pit | pits | gold | east", + Verb="is | see | smell | shoot | fell | stinks | go | grab | carry | kill | turn | feel", # noqa + Adjective="right | left | east | south | back | smelly | dead", + Adverb="here | there | nearby | ahead | right | left | east | south | back", + Pronoun="me | you | I | it", + Name="John | Mary | Boston | Aristotle", + Article="the | a | an", + Preposition="to | in | on | near", + Conjunction="and | or | but", + Digit="0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9", + That="that" + )) + +E_ = Grammar('E_', # Trivial Grammar and lexicon for testing + Rules( + S='NP VP', + NP='Art N | Pronoun', + VP='V NP'), + + Lexicon( + Art='the | a', + N='man | woman | table | shoelace | saw', + Pronoun='I | you | it', + V='saw | liked | feel' + )) + +E_NP_ = Grammar('E_NP_', # Another Trivial Grammar for testing + Rules(NP='Adj NP | N'), + Lexicon(Adj='happy | handsome | hairy', + N='man')) + +E_Prob = ProbGrammar('E_Prob', # The Probabilistic Grammar from the notebook + ProbRules( + S="NP VP [0.6] | S Conjunction S [0.4]", + NP="Pronoun [0.2] | Name [0.05] | Noun [0.2] | Article Noun [0.15] \ + | Article Adjs Noun [0.1] | Digit [0.05] | NP PP [0.15] | NP RelClause [0.1]", + VP="Verb [0.3] | VP NP [0.2] | VP Adjective [0.25] | VP PP [0.15] | VP Adverb [0.1]", + Adjs="Adjective [0.5] | Adjective Adjs [0.5]", + PP="Preposition NP [1]", + RelClause="RelPro VP [1]" + ), + ProbLexicon( + Verb="is [0.5] | say [0.3] | are [0.2]", + Noun="robot [0.4] | sheep [0.4] | fence [0.2]", + Adjective="good [0.5] | new [0.2] | sad [0.3]", + Adverb="here [0.6] | lightly [0.1] | now [0.3]", + Pronoun="me [0.3] | you [0.4] | he [0.3]", + RelPro="that [0.5] | who [0.3] | which [0.2]", + Name="john [0.4] | mary [0.4] | peter [0.2]", + Article="the [0.5] | a [0.25] | an [0.25]", + Preposition="to [0.4] | in [0.3] | at [0.3]", + Conjunction="and [0.5] | or [0.2] | but [0.3]", + Digit="0 [0.35] | 1 [0.35] | 2 [0.3]" + )) + + +E_Chomsky = Grammar('E_Prob_Chomsky', # A Grammar in Chomsky Normal Form + Rules( + S='NP VP', + NP='Article Noun | Adjective Noun', + VP='Verb NP | Verb Adjective', + ), + Lexicon( + Article='the | a | an', + Noun='robot | sheep | fence', + Adjective='good | new | sad', + Verb='is | say | are' + )) + +E_Prob_Chomsky = ProbGrammar('E_Prob_Chomsky', # A Probabilistic Grammar in CNF + ProbRules( + S='NP VP [1]', + NP='Article Noun [0.6] | Adjective Noun [0.4]', + VP='Verb NP [0.5] | Verb Adjective [0.5]', + ), + ProbLexicon( + Article='the [0.5] | a [0.25] | an [0.25]', + Noun='robot [0.4] | sheep [0.4] | fence [0.2]', + Adjective='good [0.5] | new [0.2] | sad [0.3]', + Verb='is [0.5] | say [0.3] | are [0.2]' + )) +E_Prob_Chomsky_ = ProbGrammar('E_Prob_Chomsky_', + ProbRules( + S='NP VP [1]', + NP='NP PP [0.4] | Noun Verb [0.6]', + PP='Preposition NP [1]', + VP='Verb NP [0.7] | VP PP [0.3]', + ), + ProbLexicon( + Noun='astronomers [0.18] | eyes [0.32] | stars [0.32] | telescopes [0.18]', + Verb='saw [0.5] | \'\' [0.5]', + Preposition='with [1]' + )) + +# ______________________________________________________________________________ +# 22.3 Parsing + + +class Chart: + + """Class for parsing sentences using a chart data structure. + >>> chart = Chart(E0) + >>> len(chart.parses('the stench is in 2 2')) + 1 + """ + + def __init__(self, grammar, trace=False): + """A datastructure for parsing a string; and methods to do the parse. + self.chart[i] holds the edges that end just before the i'th word. + Edges are 5-element lists of [start, end, lhs, [found], [expects]].""" + self.grammar = grammar + self.trace = trace + + def parses(self, words, S='S'): + """Return a list of parses; words can be a list or string.""" + if isinstance(words, str): + words = words.split() + self.parse(words, S) + # Return all the parses that span the whole input + # 'span the whole input' => begin at 0, end at len(words) + return [[i, j, S, found, []] + for (i, j, lhs, found, expects) in self.chart[len(words)] + # assert j == len(words) + if i == 0 and lhs == S and expects == []] + + def parse(self, words, S='S'): + """Parse a list of words; according to the grammar. + Leave results in the chart.""" + self.chart = [[] for i in range(len(words)+1)] + self.add_edge([0, 0, 'S_', [], [S]]) + for i in range(len(words)): + self.scanner(i, words[i]) + return self.chart + + def add_edge(self, edge): + """Add edge to chart, and see if it extends or predicts another edge.""" + start, end, lhs, found, expects = edge + if edge not in self.chart[end]: + self.chart[end].append(edge) + if self.trace: + print('Chart: added {}'.format(edge)) + if not expects: + self.extender(edge) + else: + self.predictor(edge) + + def scanner(self, j, word): + """For each edge expecting a word of this category here, extend the edge.""" + for (i, j, A, alpha, Bb) in self.chart[j]: + if Bb and self.grammar.isa(word, Bb[0]): + self.add_edge([i, j+1, A, alpha + [(Bb[0], word)], Bb[1:]]) + + def predictor(self, edge): + """Add to chart any rules for B that could help extend this edge.""" + (i, j, A, alpha, Bb) = edge + B = Bb[0] + if B in self.grammar.rules: + for rhs in self.grammar.rewrites_for(B): + self.add_edge([j, j, B, [], rhs]) + + def extender(self, edge): + """See what edges can be extended by this edge.""" + (j, k, B, _, _) = edge + for (i, j, A, alpha, B1b) in self.chart[j]: + if B1b and B == B1b[0]: + self.add_edge([i, k, A, alpha + [edge], B1b[1:]]) + + +# ______________________________________________________________________________ +# CYK Parsing + + +class Tree: + def __init__(self, root, *args): + self.root = root + self.leaves = [leaf for leaf in args] + + +def CYK_parse(words, grammar): + """ [Figure 22.6] """ + # We use 0-based indexing instead of the book's 1-based. + P = defaultdict(float) + T = defaultdict(Tree) + + # Insert lexical categories for each word. + for (i, word) in enumerate(words): + for (X, p) in grammar.categories[word]: + P[X, i, i] = p + T[X, i, i] = Tree(X, word) + + # Construct X(i:k) from Y(i:j) and Z(j+1:k), shortest span first + for i, j, k in subspan(len(words)): + for (X, Y, Z, p) in grammar.cnf_rules(): + PYZ = P[Y, i, j] * P[Z, j+1, k] * p + if PYZ > P[X, i, k]: + P[X, i, k] = PYZ + T[X, i, k] = Tree(X, T[Y, i, j], T[Z, j+1, k]) + + return T + + +def subspan(N): + """returns all tuple(i, j, k) covering a span (i, k) with i <= j < k""" + for length in range(2, N+1): + for i in range(1, N+2-length): + k = i + length - 1 + for j in range(i, k): + yield (i, j, k) + +# using search algorithms in the searching part + + +class TextParsingProblem(Problem): + def __init__(self, initial, grammar, goal='S'): + """ + :param initial: the initial state of words in a list. + :param grammar: a grammar object + :param goal: the goal state, usually S + """ + super(TextParsingProblem, self).__init__(initial, goal) + self.grammar = grammar + self.combinations = defaultdict(list) # article combinations + # backward lookup of rules + for rule in grammar.rules: + for comb in grammar.rules[rule]: + self.combinations[' '.join(comb)].append(rule) + + def actions(self, state): + actions = [] + categories = self.grammar.categories + # first change each word to the article of its category + for i in range(len(state)): + word = state[i] + if word in categories: + for X in categories[word]: + state[i] = X + actions.append(copy.copy(state)) + state[i] = word + # if all words are replaced by articles, replace combinations of articles by inferring rules. + if not actions: + for start in range(len(state)): + for end in range(start, len(state)+1): + # try combinations between (start, end) + articles = ' '.join(state[start:end]) + for c in self.combinations[articles]: + actions.append(state[:start] + [c] + state[end:]) + return actions + + def result(self, state, action): + return action + + def h(self, state): + # heuristic function + return len(state) + + +def astar_search_parsing(words, gramma): + """bottom-up parsing using A* search to find whether a list of words is a sentence""" + # init the problem + problem = TextParsingProblem(words, gramma, 'S') + state = problem.initial + # init the searching frontier + frontier = [(len(state)+problem.h(state), state)] + heapq.heapify(frontier) + + while frontier: + # search the frontier node with lowest cost first + cost, state = heapq.heappop(frontier) + actions = problem.actions(state) + for action in actions: + new_state = problem.result(state, action) + # update the new frontier node to the frontier + if new_state == [problem.goal]: + return problem.goal + if new_state != state: + heapq.heappush(frontier, (len(new_state)+problem.h(new_state), new_state)) + return False + + +def beam_search_parsing(words, gramma, b=3): + """bottom-up text parsing using beam search""" + # init problem + problem = TextParsingProblem(words, gramma, 'S') + # init frontier + frontier = [(len(problem.initial), problem.initial)] + heapq.heapify(frontier) + + # explore the current frontier and keep b new states with lowest cost + def explore(frontier): + new_frontier = [] + for cost, state in frontier: + # expand the possible children states of current state + if not problem.goal_test(' '.join(state)): + actions = problem.actions(state) + for action in actions: + new_state = problem.result(state, action) + if [len(new_state), new_state] not in new_frontier and new_state != state: + new_frontier.append([len(new_state), new_state]) + else: + return problem.goal + heapq.heapify(new_frontier) + # only keep b states + return heapq.nsmallest(b, new_frontier) + + while frontier: + frontier = explore(frontier) + if frontier == problem.goal: + return frontier + return False + +# ______________________________________________________________________________ +# 22.4 Augmented Grammar + + +g = Grammar("arithmetic_expression", # A Grammar of Arithmetic Expression + rules={ + 'Number_0': 'Digit_0', 'Number_1': 'Digit_1', 'Number_2': 'Digit_2', + 'Number_10': 'Number_1 Digit_0', 'Number_11': 'Number_1 Digit_1', + 'Number_100': 'Number_10 Digit_0', + 'Exp_5': ['Number_5', '( Exp_5 )', 'Exp_1, Operator_+ Exp_4', 'Exp_2, Operator_+ Exp_3', + 'Exp_0, Operator_+ Exp_5', 'Exp_3, Operator_+ Exp_2', 'Exp_4, Operator_+ Exp_1', + 'Exp_5, Operator_+ Exp_0', 'Exp_1, Operator_* Exp_5'], # more possible combinations + 'Operator_+': operator.add, 'Operator_-': operator.sub, 'Operator_*':operator.mul, 'Operator_/': operator.truediv, + 'Digit_0': 0, 'Digit_1': 1, 'Digit_2': 2, 'Digit_3': 3, 'Digit_4': 4 + }, + lexicon={}) + +g = Grammar("Ali loves Bob", # A example grammer of Ali loves Bob example + rules={ + "S_loves_ali_bob": "NP_ali, VP_x_loves_x_bob", "S_loves_bob_ali": "NP_bob, VP_x_loves_x_ali", + "VP_x_loves_x_bob": "Verb_xy_loves_xy NP_bob", "VP_x_loves_x_ali": "Verb_xy_loves_xy NP_ali", + "NP_bob": "Name_bob", "NP_ali": "Name_ali" + }, + lexicon={ + "Name_ali":"Ali", "Name_bob": "Bob", "Verb_xy_loves_xy": "loves" + }) + + diff --git a/requirements.txt b/requirements.txt index 8032818cc..3d8754e71 100644 --- a/requirements.txt +++ b/requirements.txt @@ -9,4 +9,4 @@ ipythonblocks keras numpy tensorflow -opencv-python \ No newline at end of file +opencv-python diff --git a/rl4e.py b/rl4e.py new file mode 100644 index 000000000..5575d8173 --- /dev/null +++ b/rl4e.py @@ -0,0 +1,340 @@ +"""Reinforcement Learning (Chapter 21)""" + +from collections import defaultdict +from utils import argmax +from mdp import MDP, policy_evaluation + +import random + +# _________________________________________ +# 21.2 Passive Reinforcement Learning +# 21.2.1 Direct utility estimation + + +class PassiveDUEAgent: + """Passive (non-learning) agent that uses direct utility estimation + on a given MDP and policy. + + import sys + from mdp import sequential_decision_environment + north = (0, 1) + south = (0,-1) + west = (-1, 0) + east = (1, 0) + policy = {(0, 2): east, (1, 2): east, (2, 2): east, (3, 2): None, (0, 1): north, (2, 1): north, (3, 1): None, (0, 0): north, (1, 0): west, (2, 0): west, (3, 0): west,} + agent = PassiveDUEAgent(policy, sequential_decision_environment) + for i in range(200): + run_single_trial(agent,sequential_decision_environment) + agent.estimate_U() + agent.U[(0, 0)] > 0.2 + True + + """ + + def __init__(self, pi, mdp): + self.pi = pi + self.mdp = mdp + self.U = {} + self.s = None + self.a = None + self.s_history = [] + self.r_history = [] + self.init = mdp.init + + def __call__(self, percept): + s1, r1 = percept + self.s_history.append(s1) + self.r_history.append(r1) + ## + ## + if s1 in self.mdp.terminals: + self.s = self.a = None + else: + self.s, self.a = s1, self.pi[s1] + return self.a + + def estimate_U(self): + # this function can be called only if the MDP has reached a terminal state + # it will also reset the mdp history + assert self.a is None, 'MDP is not in terminal state' + assert len(self.s_history) == len(self.r_history) + # calculating the utilities based on the current iteration + U2 = {s: [] for s in set(self.s_history)} + for i in range(len(self.s_history)): + s = self.s_history[i] + U2[s] += [sum(self.r_history[i:])] + U2 = {k: sum(v) / max(len(v), 1) for k, v in U2.items()} + # resetting history + self.s_history, self.r_history = [], [] + # setting the new utilities to the average of the previous + # iteration and this one + for k in U2.keys(): + if k in self.U.keys(): + self.U[k] = (self.U[k] + U2[k]) / 2 + else: + self.U[k] = U2[k] + return self.U + + def update_state(self, percept): + '''To be overridden in most cases. The default case + assumes the percept to be of type (state, reward)''' + return percept + +# 21.2.2 Adaptive dynamic programming + + +class PassiveADPAgent: + + """Passive (non-learning) agent that uses adaptive dynamic programming + on a given MDP and policy. [Figure 21.2] + + import sys + from mdp import sequential_decision_environment + north = (0, 1) + south = (0,-1) + west = (-1, 0) + east = (1, 0) + policy = {(0, 2): east, (1, 2): east, (2, 2): east, (3, 2): None, (0, 1): north, (2, 1): north, (3, 1): None, (0, 0): north, (1, 0): west, (2, 0): west, (3, 0): west,} + agent = PassiveADPAgent(policy, sequential_decision_environment) + for i in range(100): + run_single_trial(agent,sequential_decision_environment) + + agent.U[(0, 0)] > 0.2 + True + agent.U[(0, 1)] > 0.2 + True + """ + + class ModelMDP(MDP): + """ Class for implementing modified Version of input MDP with + an editable transition model P and a custom function T. """ + def __init__(self, init, actlist, terminals, gamma, states): + super().__init__(init, actlist, terminals, states=states, gamma=gamma) + nested_dict = lambda: defaultdict(nested_dict) + # StackOverflow:whats-the-best-way-to-initialize-a-dict-of-dicts-in-python + self.P = nested_dict() + + def T(self, s, a): + """Return a list of tuples with probabilities for states + based on the learnt model P.""" + return [(prob, res) for (res, prob) in self.P[(s, a)].items()] + + def __init__(self, pi, mdp): + self.pi = pi + self.mdp = PassiveADPAgent.ModelMDP(mdp.init, mdp.actlist, + mdp.terminals, mdp.gamma, mdp.states) + self.U = {} + self.Nsa = defaultdict(int) + self.Ns1_sa = defaultdict(int) + self.s = None + self.a = None + self.visited = set() # keeping track of visited states + + def __call__(self, percept): + s1, r1 = percept + mdp = self.mdp + R, P, terminals, pi = mdp.reward, mdp.P, mdp.terminals, self.pi + s, a, Nsa, Ns1_sa, U = self.s, self.a, self.Nsa, self.Ns1_sa, self.U + + if s1 not in self.visited: # Reward is only known for visited state. + U[s1] = R[s1] = r1 + self.visited.add(s1) + if s is not None: + Nsa[(s, a)] += 1 + Ns1_sa[(s1, s, a)] += 1 + # for each t such that Ns′|sa [t, s, a] is nonzero + for t in [res for (res, state, act), freq in Ns1_sa.items() + if (state, act) == (s, a) and freq != 0]: + P[(s, a)][t] = Ns1_sa[(t, s, a)] / Nsa[(s, a)] + + self.U = policy_evaluation(pi, U, mdp) + ## + ## + self.Nsa, self.Ns1_sa = Nsa, Ns1_sa + if s1 in terminals: + self.s = self.a = None + else: + self.s, self.a = s1, self.pi[s1] + return self.a + + def update_state(self, percept): + """To be overridden in most cases. The default case + assumes the percept to be of type (state, reward).""" + return percept + +# 21.2.3 Temporal-difference learning + + +class PassiveTDAgent: + """The abstract class for a Passive (non-learning) agent that uses + temporal differences to learn utility estimates. Override update_state + method to convert percept to state and reward. The mdp being provided + should be an instance of a subclass of the MDP Class. [Figure 21.4] + + import sys + from mdp import sequential_decision_environment + north = (0, 1) + south = (0,-1) + west = (-1, 0) + east = (1, 0) + policy = {(0, 2): east, (1, 2): east, (2, 2): east, (3, 2): None, (0, 1): north, (2, 1): north, (3, 1): None, (0, 0): north, (1, 0): west, (2, 0): west, (3, 0): west,} + agent = PassiveTDAgent(policy, sequential_decision_environment, alpha=lambda n: 60./(59+n)) + for i in range(200): + run_single_trial(agent,sequential_decision_environment) + + agent.U[(0, 0)] > 0.2 + True + agent.U[(0, 1)] > 0.2 + True + """ + + def __init__(self, pi, mdp, alpha=None): + + self.pi = pi + self.U = {s: 0. for s in mdp.states} + self.Ns = {s: 0 for s in mdp.states} + self.s = None + self.a = None + self.r = None + self.gamma = mdp.gamma + self.terminals = mdp.terminals + + if alpha: + self.alpha = alpha + else: + self.alpha = lambda n: 1 / (1 + n) # udacity video + + def __call__(self, percept): + s1, r1 = self.update_state(percept) + pi, U, Ns, s, r = self.pi, self.U, self.Ns, self.s, self.r + alpha, gamma, terminals = self.alpha, self.gamma, self.terminals + if not Ns[s1]: + U[s1] = r1 + if s is not None: + Ns[s] += 1 + U[s] += alpha(Ns[s]) * (r + gamma * U[s1] - U[s]) + if s1 in terminals: + self.s = self.a = self.r = None + else: + self.s, self.a, self.r = s1, pi[s1], r1 + return self.a + + def update_state(self, percept): + """To be overridden in most cases. The default case + assumes the percept to be of type (state, reward).""" + return percept + +# __________________________________________ +# 21.3. Active Reinforcement Learning +# 21.3.2 Learning an action-utility function + + +class QLearningAgent: + """ An exploratory Q-learning agent. It avoids having to learn the transition + model because the Q-value of a state can be related directly to those of + its neighbors. [Figure 21.8] + + import sys + from mdp import sequential_decision_environment + north = (0, 1) + south = (0,-1) + west = (-1, 0) + east = (1, 0) + policy = {(0, 2): east, (1, 2): east, (2, 2): east, (3, 2): None, (0, 1): north, (2, 1): north, (3, 1): None, (0, 0): north, (1, 0): west, (2, 0): west, (3, 0): west,} + q_agent = QLearningAgent(sequential_decision_environment, Ne=5, Rplus=2, alpha=lambda n: 60./(59+n)) + for i in range(200): + run_single_trial(q_agent,sequential_decision_environment) + + q_agent.Q[((0, 1), (0, 1))] >= -0.5 + True + q_agent.Q[((1, 0), (0, -1))] <= 0.5 + True + """ + + def __init__(self, mdp, Ne, Rplus, alpha=None): + + self.gamma = mdp.gamma + self.terminals = mdp.terminals + self.all_act = mdp.actlist + self.Ne = Ne # iteration limit in exploration function + self.Rplus = Rplus # large value to assign before iteration limit + self.Q = defaultdict(float) + self.Nsa = defaultdict(float) + self.s = None + self.a = None + self.r = None + + if alpha: + self.alpha = alpha + else: + self.alpha = lambda n: 1. / (1 + n) # udacity video + + def f(self, u, n): + """ Exploration function. Returns fixed Rplus until + agent has visited state, action a Ne number of times. + Same as ADP agent in book.""" + if n < self.Ne: + return self.Rplus + else: + return u + + def actions_in_state(self, state): + """ Return actions possible in given state. + Useful for max and argmax. """ + if state in self.terminals: + return [None] + else: + return self.all_act + + def __call__(self, percept): + s1, r1 = self.update_state(percept) + Q, Nsa, s, a, r = self.Q, self.Nsa, self.s, self.a, self.r + alpha, gamma, terminals = self.alpha, self.gamma, self.terminals, + actions_in_state = self.actions_in_state + + if s in terminals: + Q[s, None] = r1 + if s is not None: + Nsa[s, a] += 1 + Q[s, a] += alpha(Nsa[s, a]) * (r + gamma * max(Q[s1, a1] + for a1 in actions_in_state(s1)) - Q[s, a]) + if s in terminals: + self.s = self.a = self.r = None + else: + self.s, self.r = s1, r1 + self.a = argmax(actions_in_state(s1), key=lambda a1: self.f(Q[s1, a1], Nsa[s1, a1])) + return self.a + + def update_state(self, percept): + """To be overridden in most cases. The default case + assumes the percept to be of type (state, reward).""" + return percept + + +def run_single_trial(agent_program, mdp): + """Execute trial for given agent_program + and mdp. mdp should be an instance of subclass + of mdp.MDP """ + + def take_single_action(mdp, s, a): + """ + Select outcome of taking action a + in state s. Weighted Sampling. + """ + x = random.uniform(0, 1) + cumulative_probability = 0.0 + for probability_state in mdp.T(s, a): + probability, state = probability_state + cumulative_probability += probability + if x < cumulative_probability: + break + return state + + current_state = mdp.init + while True: + current_reward = mdp.R(current_state) + percept = (current_state, current_reward) + next_action = agent_program(percept) + if next_action is None: + break + current_state = take_single_action(mdp, current_state, next_action) diff --git a/tests/test_agents.py b/tests/test_agents.py index 3c133c32a..0433396ff 100644 --- a/tests/test_agents.py +++ b/tests/test_agents.py @@ -63,7 +63,7 @@ def test_RandomAgentProgram() : list = ['Right', 'Left', 'Suck', 'NoOp'] # create a program and then an object of the RandomAgentProgram program = RandomAgentProgram(list) - + agent = Agent(program) # create an object of TrivialVacuumEnvironment environment = TrivialVacuumEnvironment() @@ -139,26 +139,26 @@ def test_ReflexVacuumAgent() : def test_SimpleReflexAgentProgram(): class Rule: - + def __init__(self, state, action): self.__state = state self.action = action - + def matches(self, state): return self.__state == state - + loc_A = (0, 0) loc_B = (1, 0) - + # create rules for a two state Vacuum Environment rules = [Rule((loc_A, "Dirty"), "Suck"), Rule((loc_A, "Clean"), "Right"), Rule((loc_B, "Dirty"), "Suck"), Rule((loc_B, "Clean"), "Left")] - + def interpret_input(state): return state - + # create a program and then an object of the SimpleReflexAgentProgram - program = SimpleReflexAgentProgram(rules, interpret_input) + program = SimpleReflexAgentProgram(rules, interpret_input) agent = Agent(program) # create an object of TrivialVacuumEnvironment environment = TrivialVacuumEnvironment() @@ -306,8 +306,8 @@ def constant_prog(percept): assert not any(map(lambda x: not isinstance(x,Thing), w.things)) #Check that gold and wumpus are not present on (1,1) - assert not any(map(lambda x: isinstance(x, Gold) or isinstance(x,WumpusEnvironment), - w.list_things_at((1, 1)))) + assert not any(map(lambda x: isinstance(x, Gold) or isinstance(x,WumpusEnvironment), + w.list_things_at((1, 1)))) #Check if w.get_world() segments objects correctly assert len(w.get_world()) == 6 diff --git a/tests/test_deepNN.py b/tests/test_deepNN.py new file mode 100644 index 000000000..0a98b7e76 --- /dev/null +++ b/tests/test_deepNN.py @@ -0,0 +1,74 @@ +from DeepNeuralNet4e import * +from learning4e import DataSet, grade_learner, err_ratio +from keras.datasets import imdb +import numpy as np + + +def test_neural_net(): + iris = DataSet(name="iris") + classes = ["setosa", "versicolor", "virginica"] + iris.classes_to_numbers(classes) + nn_adam = neural_net_learner(iris, [4], learning_rate=0.001, epochs=200, optimizer=adam_optimizer) + nn_gd = neural_net_learner(iris, [4], learning_rate=0.15, epochs=100, optimizer=gradient_descent) + tests = [([5.0, 3.1, 0.9, 0.1], 0), + ([5.1, 3.5, 1.0, 0.0], 0), + ([4.9, 3.3, 1.1, 0.1], 0), + ([6.0, 3.0, 4.0, 1.1], 1), + ([6.1, 2.2, 3.5, 1.0], 1), + ([5.9, 2.5, 3.3, 1.1], 1), + ([7.5, 4.1, 6.2, 2.3], 2), + ([7.3, 4.0, 6.1, 2.4], 2), + ([7.0, 3.3, 6.1, 2.5], 2)] + assert grade_learner(nn_adam, tests) >= 1 / 3 + assert grade_learner(nn_gd, tests) >= 1 / 3 + assert err_ratio(nn_adam, iris) < 0.21 + assert err_ratio(nn_gd, iris) < 0.21 + + +def test_cross_entropy(): + loss = cross_entropy_loss([1,0], [0.9, 0.3]) + assert round(loss,2) == 0.23 + + loss = cross_entropy_loss([1,0,0,1], [0.9,0.3,0.5,0.75]) + assert round(loss,2) == 0.36 + + loss = cross_entropy_loss([1,0,0,1,1,0,1,1], [0.9,0.3,0.5,0.75,0.85,0.14,0.93,0.79]) + assert round(loss,2) == 0.26 + + +def test_perceptron(): + iris = DataSet(name="iris") + classes = ["setosa", "versicolor", "virginica"] + iris.classes_to_numbers(classes) + perceptron = perceptron_learner(iris, learning_rate=0.01, epochs=100) + tests = [([5, 3, 1, 0.1], 0), + ([5, 3.5, 1, 0], 0), + ([6, 3, 4, 1.1], 1), + ([6, 2, 3.5, 1], 1), + ([7.5, 4, 6, 2], 2), + ([7, 3, 6, 2.5], 2)] + assert grade_learner(perceptron, tests) > 1/2 + assert err_ratio(perceptron, iris) < 0.4 + + +def test_rnn(): + data = imdb.load_data(num_words=5000) + train, val, test = keras_dataset_loader(data) + train = (train[0][:1000], train[1][:1000]) + val = (val[0][:200], val[1][:200]) + model = simple_rnn_learner(train, val) + score = model.evaluate(test[0][:200], test[1][:200], verbose=0) + acc = score[1] + assert acc >= 0.3 + + +def test_auto_encoder(): + iris = DataSet(name="iris") + classes = ["setosa", "versicolor", "virginica"] + iris.classes_to_numbers(classes) + inputs = np.asarray(iris.examples) + # print(inputs[0]) + model = auto_encoder_learner(inputs, 100) + print(inputs[0]) + print(model.predict(inputs[:1])) + diff --git a/tests/test_learning4e.py b/tests/test_learning4e.py new file mode 100644 index 000000000..e80ccdd04 --- /dev/null +++ b/tests/test_learning4e.py @@ -0,0 +1,103 @@ +import pytest +import math +import random +from utils import open_data +from learning import * + + +random.seed("aima-python") + + +def test_mean_boolean_error(): + assert mean_boolean_error([1, 1], [0, 0]) == 1 + assert mean_boolean_error([0, 1], [1, 0]) == 1 + assert mean_boolean_error([1, 1], [0, 1]) == 0.5 + assert mean_boolean_error([0, 0], [0, 0]) == 0 + assert mean_boolean_error([1, 1], [1, 1]) == 0 + + +def test_exclude(): + iris = DataSet(name='iris', exclude=[3]) + assert iris.inputs == [0, 1, 2] + + +def test_parse_csv(): + Iris = open_data('iris.csv').read() + assert parse_csv(Iris)[0] == [5.1, 3.5, 1.4, 0.2, 'setosa'] + + +def test_weighted_mode(): + assert weighted_mode('abbaa', [1, 2, 3, 1, 2]) == 'b' + + +def test_weighted_replicate(): + assert weighted_replicate('ABC', [1, 2, 1], 4) == ['A', 'B', 'B', 'C'] + + +def test_means_and_deviation(): + iris = DataSet(name="iris") + + means, deviations = iris.find_means_and_deviations() + + assert round(means["setosa"][0], 3) == 5.006 + assert round(means["versicolor"][0], 3) == 5.936 + assert round(means["virginica"][0], 3) == 6.588 + + assert round(deviations["setosa"][0], 3) == 0.352 + assert round(deviations["versicolor"][0], 3) == 0.516 + assert round(deviations["virginica"][0], 3) == 0.636 + + +def test_decision_tree_learner(): + iris = DataSet(name="iris") + dTL = DecisionTreeLearner(iris) + assert dTL([5, 3, 1, 0.1]) == "setosa" + assert dTL([6, 5, 3, 1.5]) == "versicolor" + assert dTL([7.5, 4, 6, 2]) == "virginica" + + +def test_information_content(): + assert information_content([]) == 0 + assert information_content([4]) == 0 + assert information_content([5, 4, 0, 2, 5, 0]) > 1.9 + assert information_content([5, 4, 0, 2, 5, 0]) < 2 + assert information_content([1.5, 2.5]) > 0.9 + assert information_content([1.5, 2.5]) < 1.0 + + +def test_random_forest(): + iris = DataSet(name="iris") + rF = RandomForest(iris) + tests = [([5.0, 3.0, 1.0, 0.1], "setosa"), + ([5.1, 3.3, 1.1, 0.1], "setosa"), + ([6.0, 5.0, 3.0, 1.0], "versicolor"), + ([6.1, 2.2, 3.5, 1.0], "versicolor"), + ([7.5, 4.1, 6.2, 2.3], "virginica"), + ([7.3, 3.7, 6.1, 2.5], "virginica")] + assert grade_learner(rF, tests) >= 1/3 + + +def test_random_weights(): + min_value = -0.5 + max_value = 0.5 + num_weights = 10 + test_weights = random_weights(min_value, max_value, num_weights) + assert len(test_weights) == num_weights + for weight in test_weights: + assert weight >= min_value and weight <= max_value + + +def test_adaboost(): + iris = DataSet(name="iris") + iris.classes_to_numbers() + WeightedPerceptron = WeightedLearner(PerceptronLearner) + AdaboostLearner = AdaBoost(WeightedPerceptron, 5) + adaboost = AdaboostLearner(iris) + tests = [([5, 3, 1, 0.1], 0), + ([5, 3.5, 1, 0], 0), + ([6, 3, 4, 1.1], 1), + ([6, 2, 3.5, 1], 1), + ([7.5, 4, 6, 2], 2), + ([7, 3, 6, 2.5], 2)] + assert grade_learner(adaboost, tests) > 4/6 + assert err_ratio(adaboost, iris) < 0.25 diff --git a/tests/test_nlp4e.py b/tests/test_nlp4e.py new file mode 100644 index 000000000..029cbaf22 --- /dev/null +++ b/tests/test_nlp4e.py @@ -0,0 +1,135 @@ +import pytest +import nlp + +from nlp4e import Rules, Lexicon, Grammar, ProbRules, ProbLexicon, ProbGrammar, E0 +from nlp4e import Chart, CYK_parse, subspan, astar_search_parsing, beam_search_parsing +# Clumsy imports because we want to access certain nlp.py globals explicitly, because +# they are accessed by functions within nlp.py + + +def test_rules(): + check = {'A': [['B', 'C'], ['D', 'E']], 'B': [['E'], ['a'], ['b', 'c']]} + assert Rules(A="B C | D E", B="E | a | b c") == check + + +def test_lexicon(): + check = {'Article': ['the', 'a', 'an'], 'Pronoun': ['i', 'you', 'he']} + lexicon = Lexicon(Article="the | a | an", Pronoun="i | you | he") + assert lexicon == check + + +def test_grammar(): + rules = Rules(A="B C | D E", B="E | a | b c") + lexicon = Lexicon(Article="the | a | an", Pronoun="i | you | he") + grammar = Grammar("Simplegram", rules, lexicon) + + assert grammar.rewrites_for('A') == [['B', 'C'], ['D', 'E']] + assert grammar.isa('the', 'Article') + + grammar = nlp.E_Chomsky + for rule in grammar.cnf_rules(): + assert len(rule) == 3 + + +def test_generation(): + lexicon = Lexicon(Article="the | a | an", + Pronoun="i | you | he") + + rules = Rules( + S="Article | More | Pronoun", + More="Article Pronoun | Pronoun Pronoun" + ) + + grammar = Grammar("Simplegram", rules, lexicon) + + sentence = grammar.generate_random('S') + for token in sentence.split(): + found = False + for non_terminal, terminals in grammar.lexicon.items(): + if token in terminals: + found = True + assert found + + +def test_prob_rules(): + check = {'A': [(['B', 'C'], 0.3), (['D', 'E'], 0.7)], + 'B': [(['E'], 0.1), (['a'], 0.2), (['b', 'c'], 0.7)]} + rules = ProbRules(A="B C [0.3] | D E [0.7]", B="E [0.1] | a [0.2] | b c [0.7]") + assert rules == check + + +def test_prob_lexicon(): + check = {'Article': [('the', 0.5), ('a', 0.25), ('an', 0.25)], + 'Pronoun': [('i', 0.4), ('you', 0.3), ('he', 0.3)]} + lexicon = ProbLexicon(Article="the [0.5] | a [0.25] | an [0.25]", + Pronoun="i [0.4] | you [0.3] | he [0.3]") + assert lexicon == check + + +def test_prob_grammar(): + rules = ProbRules(A="B C [0.3] | D E [0.7]", B="E [0.1] | a [0.2] | b c [0.7]") + lexicon = ProbLexicon(Article="the [0.5] | a [0.25] | an [0.25]", + Pronoun="i [0.4] | you [0.3] | he [0.3]") + grammar = ProbGrammar("Simplegram", rules, lexicon) + + assert grammar.rewrites_for('A') == [(['B', 'C'], 0.3), (['D', 'E'], 0.7)] + assert grammar.isa('the', 'Article') + + grammar = nlp.E_Prob_Chomsky + for rule in grammar.cnf_rules(): + assert len(rule) == 4 + + +def test_prob_generation(): + lexicon = ProbLexicon(Verb="am [0.5] | are [0.25] | is [0.25]", + Pronoun="i [0.4] | you [0.3] | he [0.3]") + + rules = ProbRules( + S="Verb [0.5] | More [0.3] | Pronoun [0.1] | nobody is here [0.1]", + More="Pronoun Verb [0.7] | Pronoun Pronoun [0.3]" + ) + + grammar = ProbGrammar("Simplegram", rules, lexicon) + + sentence = grammar.generate_random('S') + assert len(sentence) == 2 + + +def test_chart_parsing(): + chart = Chart(nlp.E0) + parses = chart.parses('the stench is in 2 2') + assert len(parses) == 1 + + +def test_CYK_parse(): + grammar = nlp.E_Prob_Chomsky + words = ['the', 'robot', 'is', 'good'] + P = CYK_parse(words, grammar) + assert len(P) == 5 + + grammar = nlp.E_Prob_Chomsky_ + words = ['astronomers', 'saw', 'stars'] + P = CYK_parse(words, grammar) + assert len(P) == 3 + + +def test_subspan(): + spans = subspan(3) + assert spans.__next__() == (1,1,2) + assert spans.__next__() == (2,2,3) + assert spans.__next__() == (1,1,3) + assert spans.__next__() == (1,2,3) + + +def test_text_parsing(): + words = ["the", "wumpus", "is", "dead"] + grammer = E0 + assert astar_search_parsing(words, grammer) == 'S' + assert beam_search_parsing(words, grammer) == 'S' + words = ["the", "is", "wupus", "dead"] + assert astar_search_parsing(words, grammer) == False + assert beam_search_parsing(words, grammer) == False + + +if __name__ == '__main__': + pytest.main() diff --git a/tests/test_rl4e.py b/tests/test_rl4e.py new file mode 100644 index 000000000..d9c2c672d --- /dev/null +++ b/tests/test_rl4e.py @@ -0,0 +1,66 @@ +import pytest + +from rl4e import * +from mdp import sequential_decision_environment + + +north = (0, 1) +south = (0,-1) +west = (-1, 0) +east = (1, 0) + +policy = { + (0, 2): east, (1, 2): east, (2, 2): east, (3, 2): None, + (0, 1): north, (2, 1): north, (3, 1): None, + (0, 0): north, (1, 0): west, (2, 0): west, (3, 0): west, +} + +def test_PassiveDUEAgent(): + agent = PassiveDUEAgent(policy, sequential_decision_environment) + for i in range(200): + run_single_trial(agent,sequential_decision_environment) + agent.estimate_U() + # Agent does not always produce same results. + # Check if results are good enough. + #print(agent.U[(0, 0)], agent.U[(0,1)], agent.U[(1,0)]) + assert agent.U[(0, 0)] > 0.15 # In reality around 0.3 + assert agent.U[(0, 1)] > 0.15 # In reality around 0.4 + assert agent.U[(1, 0)] > 0 # In reality around 0.2 + +def test_PassiveADPAgent(): + agent = PassiveADPAgent(policy, sequential_decision_environment) + for i in range(100): + run_single_trial(agent,sequential_decision_environment) + + # Agent does not always produce same results. + # Check if results are good enough. + #print(agent.U[(0, 0)], agent.U[(0,1)], agent.U[(1,0)]) + assert agent.U[(0, 0)] > 0.15 # In reality around 0.3 + assert agent.U[(0, 1)] > 0.15 # In reality around 0.4 + assert agent.U[(1, 0)] > 0 # In reality around 0.2 + + + +def test_PassiveTDAgent(): + agent = PassiveTDAgent(policy, sequential_decision_environment, alpha=lambda n: 60./(59+n)) + for i in range(200): + run_single_trial(agent,sequential_decision_environment) + + # Agent does not always produce same results. + # Check if results are good enough. + assert agent.U[(0, 0)] > 0.15 # In reality around 0.3 + assert agent.U[(0, 1)] > 0.15 # In reality around 0.35 + assert agent.U[(1, 0)] > 0.15 # In reality around 0.25 + + +def test_QLearning(): + q_agent = QLearningAgent(sequential_decision_environment, Ne=5, Rplus=2, + alpha=lambda n: 60./(59+n)) + + for i in range(200): + run_single_trial(q_agent,sequential_decision_environment) + + # Agent does not always produce same results. + # Check if results are good enough. + assert q_agent.Q[((0, 1), (0, 1))] >= -0.5 # In reality around 0.1 + assert q_agent.Q[((1, 0), (0, -1))] <= 0.5 # In reality around -0.1 diff --git a/utils4e.py b/utils4e.py index afb60f4f0..c66020b18 100644 --- a/utils4e.py +++ b/utils4e.py @@ -420,6 +420,12 @@ def conv1D(X, K): return res + +def GaussianKernel(size=3): + mean = (size-1)/2 + stdev = 0.1 + return [gaussian(mean, stdev, x) for x in range(size)] + def gaussian_kernel_1d(size=3, sigma=0.5): mean = (size-1)/2 return [gaussian(mean, sigma, x) for x in range(size)] From fd52c720f8880bc3e406872a192c775e63ccb3b3 Mon Sep 17 00:00:00 2001 From: tianqiyang Date: Mon, 5 Aug 2019 13:58:13 -0400 Subject: [PATCH 3/3] Add chapter 12 and 13 Baysian models (#1088) * chapter 18 learning * add chapter 19 * move init dataset in NN learner * add adam optimizer, add nn learner * remove cpt 19 for debug * change while loop in games4e * add chapter 19 * add sgd and adam optimizer * add chpt19 deep nn * add rnn * add auto encoder * add comments, correct tests * add more comments, change algorithms according to orders of chapter sections * add keras and numpy to requirements * add tf as requirement * add gc in test agent * fix agent bugs for running test_agent and test_agent_4e together * fix build error * add chapter 21 and 22 * add chapter 12 and part of 13 * remove chapter 12 and 13, add test of rl * modify rnn test * add chapter 12 and 13 * change gaussian kernel util function * fix example bugs * fix build bug --- .travis.yml | 1 + DeepNeuralNet4e.py | 3 +- probability4e.py | 758 ++++++++++++++++++++++++++++++++++++ tests/test_probability4e.py | 342 ++++++++++++++++ 4 files changed, 1103 insertions(+), 1 deletion(-) create mode 100644 probability4e.py create mode 100644 tests/test_probability4e.py diff --git a/.travis.yml b/.travis.yml index b7b23e694..25750bac9 100644 --- a/.travis.yml +++ b/.travis.yml @@ -22,6 +22,7 @@ install: - pip install tensorflow - pip install opencv-python + script: - py.test --cov=./ - python -m doctest -v *.py diff --git a/DeepNeuralNet4e.py b/DeepNeuralNet4e.py index a353df95c..b68192ba8 100644 --- a/DeepNeuralNet4e.py +++ b/DeepNeuralNet4e.py @@ -1,6 +1,7 @@ import math import statistics -from utils4e import sigmoid, dotproduct, softmax1D, conv1D, GaussianKernel, element_wise_product, \ + +from utils4e import sigmoid, dotproduct, softmax1D, conv1D, gaussian_kernel_2d, GaussianKernel, element_wise_product, \ vector_add, random_weights, scalar_vector_product, matrix_multiplication, map_vector import random diff --git a/probability4e.py b/probability4e.py new file mode 100644 index 000000000..94429f2dd --- /dev/null +++ b/probability4e.py @@ -0,0 +1,758 @@ +"""Probability models. +""" + +from utils import product, argmax, isclose, probability +from logic import extend +from math import sqrt, pi, exp +import copy +import random +from collections import defaultdict +from functools import reduce + +# ______________________________________________________________________________ +# Chapter 12 Qualifying Uncertainty +# 12.1 Acting Under Uncertainty + + +def DTAgentProgram(belief_state): + """A decision-theoretic agent. [Figure 12.1]""" + def program(percept): + belief_state.observe(program.action, percept) + program.action = argmax(belief_state.actions(), + key=belief_state.expected_outcome_utility) + return program.action + program.action = None + return program + +# ______________________________________________________________________________ +# 12.2 Basic Probability Notation + + +class ProbDist: + """A discrete probability distribution. You name the random variable + in the constructor, then assign and query probability of values. + >>> P = ProbDist('Flip'); P['H'], P['T'] = 0.25, 0.75; P['H'] + 0.25 + >>> P = ProbDist('X', {'lo': 125, 'med': 375, 'hi': 500}) + >>> P['lo'], P['med'], P['hi'] + (0.125, 0.375, 0.5) + """ + + def __init__(self, varname='?', freqs=None): + """If freqs is given, it is a dictionary of values - frequency pairs, + then ProbDist is normalized.""" + self.prob = {} + self.varname = varname + self.values = [] + if freqs: + for (v, p) in freqs.items(): + self[v] = p + self.normalize() + + def __getitem__(self, val): + """Given a value, return P(value).""" + try: + return self.prob[val] + except KeyError: + return 0 + + def __setitem__(self, val, p): + """Set P(val) = p.""" + if val not in self.values: + self.values.append(val) + self.prob[val] = p + + def normalize(self): + """Make sure the probabilities of all values sum to 1. + Returns the normalized distribution. + Raises a ZeroDivisionError if the sum of the values is 0.""" + total = sum(self.prob.values()) + if not isclose(total, 1.0): + for val in self.prob: + self.prob[val] /= total + return self + + def show_approx(self, numfmt='{:.3g}'): + """Show the probabilities rounded and sorted by key, for the + sake of portable doctests.""" + return ', '.join([('{}: ' + numfmt).format(v, p) + for (v, p) in sorted(self.prob.items())]) + + def __repr__(self): + return "P({})".format(self.varname) + +# ______________________________________________________________________________ +# 12.3 Inference Using Full Joint Distributions + + +class JointProbDist(ProbDist): + """A discrete probability distribute over a set of variables. + >>> P = JointProbDist(['X', 'Y']); P[1, 1] = 0.25 + >>> P[1, 1] + 0.25 + >>> P[dict(X=0, Y=1)] = 0.5 + >>> P[dict(X=0, Y=1)] + 0.5""" + + def __init__(self, variables): + self.prob = {} + self.variables = variables + self.vals = defaultdict(list) + + def __getitem__(self, values): + """Given a tuple or dict of values, return P(values).""" + values = event_values(values, self.variables) + return ProbDist.__getitem__(self, values) + + def __setitem__(self, values, p): + """Set P(values) = p. Values can be a tuple or a dict; it must + have a value for each of the variables in the joint. Also keep track + of the values we have seen so far for each variable.""" + values = event_values(values, self.variables) + self.prob[values] = p + for var, val in zip(self.variables, values): + if val not in self.vals[var]: + self.vals[var].append(val) + + def values(self, var): + """Return the set of possible values for a variable.""" + return self.vals[var] + + def __repr__(self): + return "P({})".format(self.variables) + + +def event_values(event, variables): + """Return a tuple of the values of variables in event. + >>> event_values ({'A': 10, 'B': 9, 'C': 8}, ['C', 'A']) + (8, 10) + >>> event_values ((1, 2), ['C', 'A']) + (1, 2) + """ + if isinstance(event, tuple) and len(event) == len(variables): + return event + else: + return tuple([event[var] for var in variables]) + + +def enumerate_joint_ask(X, e, P): + """Return a probability distribution over the values of the variable X, + given the {var:val} observations e, in the JointProbDist P. [Section 12.3] + >>> P = JointProbDist(['X', 'Y']) + >>> P[0,0] = 0.25; P[0,1] = 0.5; P[1,1] = P[2,1] = 0.125 + >>> enumerate_joint_ask('X', dict(Y=1), P).show_approx() + '0: 0.667, 1: 0.167, 2: 0.167' + """ + assert X not in e, "Query variable must be distinct from evidence" + Q = ProbDist(X) # probability distribution for X, initially empty + Y = [v for v in P.variables if v != X and v not in e] # hidden variables. + for xi in P.values(X): + Q[xi] = enumerate_joint(Y, extend(e, X, xi), P) + return Q.normalize() + + +def enumerate_joint(variables, e, P): + """Return the sum of those entries in P consistent with e, + provided variables is P's remaining variables (the ones not in e).""" + if not variables: + return P[e] + Y, rest = variables[0], variables[1:] + return sum([enumerate_joint(rest, extend(e, Y, y), P) + for y in P.values(Y)]) + +# ______________________________________________________________________________ +# 12.4 Independence + + +def is_independent(variables, P): + """ + Return whether a list of variables are independent given their distribution P + P is an instance of JoinProbDist + >>> P = JointProbDist(['X', 'Y']) + >>> P[0,0] = 0.25; P[0,1] = 0.5; P[1,1] = P[1,0] = 0.125 + >>> is_independent(['X', 'Y'], P) + False + """ + for var in variables: + event_vars = variables[:] + event_vars.remove(var) + event = {} + distribution = enumerate_joint_ask(var, event, P) + events = gen_possible_events(event_vars, P) + for e in events: + conditional_distr = enumerate_joint_ask(var, e, P) + if conditional_distr.prob != distribution.prob: + return False + return True + + +def gen_possible_events(vars, P): + """Generate all possible events of a collection of vars according to distribution of P""" + events = [] + + def backtrack(vars, P, temp): + if not vars: + events.append(temp) + return + var = vars[0] + for val in P.values(var): + temp[var] = val + backtrack([v for v in vars if v != var], P, copy.copy(temp)) + backtrack(vars, P, {}) + return events + +# ______________________________________________________________________________ +# Chapter 13 Probabilistic Reasoning +# 13.1 Representing Knowledge in an Uncertain Domain + + +class BayesNet: + """Bayesian network containing only boolean-variable nodes.""" + + def __init__(self, node_specs=None): + """ + Nodes must be ordered with parents before children. + :param node_specs: an nested iterable object, each element contains (variable name, parents name, cpt) + for each node + """ + + self.nodes = [] + self.variables = [] + node_specs = node_specs or [] + for node_spec in node_specs: + self.add(node_spec) + + def add(self, node_spec): + """ + Add a node to the net. Its parents must already be in the + net, and its variable must not. + Initialize Bayes nodes by detecting the length of input node specs + """ + if len(node_spec)>=5: + node = ContinuousBayesNode(*node_spec) + else: + node = BayesNode(*node_spec) + assert node.variable not in self.variables + assert all((parent in self.variables) for parent in node.parents) + self.nodes.append(node) + self.variables.append(node.variable) + for parent in node.parents: + self.variable_node(parent).children.append(node) + + def variable_node(self, var): + """ + Return the node for the variable named var. + >>> burglary.variable_node('Burglary').variable + 'Burglary' + """ + for n in self.nodes: + if n.variable == var: + return n + raise Exception("No such variable: {}".format(var)) + + def variable_values(self, var): + """Return the domain of var.""" + return [True, False] + + def __repr__(self): + return 'BayesNet({0!r})'.format(self.nodes) + + +class BayesNode: + """ + A conditional probability distribution for a boolean variable, + P(X | parents). Part of a BayesNet. + """ + + def __init__(self, X, parents, cpt): + """ + :param X: variable name, + :param parents: a sequence of variable names or a space-separated string. Representing the names of parent nodes. + :param cpt: the conditional probability table, takes one of these forms: + + * A number, the unconditional probability P(X=true). You can + use this form when there are no parents. + + * A dict {v: p, ...}, the conditional probability distribution + P(X=true | parent=v) = p. When there's just one parent. + + * A dict {(v1, v2, ...): p, ...}, the distribution P(X=true | + parent1=v1, parent2=v2, ...) = p. Each key must have as many + values as there are parents. You can use this form always; + the first two are just conveniences. + + In all cases the probability of X being false is left implicit, + since it follows from P(X=true). + + >>> X = BayesNode('X', '', 0.2) + >>> Y = BayesNode('Y', 'P', {T: 0.2, F: 0.7}) + >>> Z = BayesNode('Z', 'P Q', + ... {(T, T): 0.2, (T, F): 0.3, (F, T): 0.5, (F, F): 0.7}) + """ + if isinstance(parents, str): + parents = parents.split() + + # We store the table always in the third form above. + if isinstance(cpt, (float, int)): # no parents, 0-tuple + cpt = {(): cpt} + elif isinstance(cpt, dict): + # one parent, 1-tuple + if cpt and isinstance(list(cpt.keys())[0], bool): + cpt = {(v,): p for v, p in cpt.items()} + + assert isinstance(cpt, dict) + for vs, p in cpt.items(): + assert isinstance(vs, tuple) and len(vs) == len(parents) + assert all(isinstance(v, bool) for v in vs) + assert 0 <= p <= 1 + + self.variable = X + self.parents = parents + self.cpt = cpt + self.children = [] + + def p(self, value, event): + """ + Return the conditional probability + P(X=value | parents=parent_values), where parent_values + are the values of parents in event. (event must assign each + parent a value.) + >>> bn = BayesNode('X', 'Burglary', {T: 0.2, F: 0.625}) + >>> bn.p(False, {'Burglary': False, 'Earthquake': True}) + 0.375 + """ + assert isinstance(value, bool) + ptrue = self.cpt[event_values(event, self.parents)] + return ptrue if value else 1 - ptrue + + def sample(self, event): + """ + Sample from the distribution for this variable conditioned + on event's values for parent_variables. That is, return True/False + at random according with the conditional probability given the + parents. + """ + return probability(self.p(True, event)) + + def __repr__(self): + return repr((self.variable, ' '.join(self.parents))) + +# Burglary example [Figure 13 .2] + + +T, F = True, False + +burglary = BayesNet([ + ('Burglary', '', 0.001), + ('Earthquake', '', 0.002), + ('Alarm', 'Burglary Earthquake', + {(T, T): 0.95, (T, F): 0.94, (F, T): 0.29, (F, F): 0.001}), + ('JohnCalls', 'Alarm', {T: 0.90, F: 0.05}), + ('MaryCalls', 'Alarm', {T: 0.70, F: 0.01}) +]) + +# ______________________________________________________________________________ +# Section 13.2. The Semantics of Bayesian Networks +# Bayesian nets with continuous variables + + +def gaussian_probability(param, event, value): + """ + Gaussian probability of a continuous Bayesian network node on condition of + certain event and the parameters determined by the event + :param param: parameters determined by discrete parent events of current node + :param event: a dict, continuous event of current node, the values are used + as parameters in calculating distribution + :param value: float, the value of current continuous node + :return: float, the calculated probability + >>> param = {'sigma':0.5, 'b':1, 'a':{'h1':0.5, 'h2': 1.5}} + >>> event = {'h1':0.6, 'h2': 0.3} + >>> gaussian_probability(param, event, 1) + 0.2590351913317835 + """ + + assert isinstance(event, dict) + assert isinstance(param, dict) + buff = 0 + for k, v in event.items(): + # buffer varianle to calculate h1*a_h1 + h2*a_h2 + buff += param['a'][k] * v + res = 1/(param['sigma']*sqrt(2*pi)) * exp(-0.5*((value-buff-param['b'])/param['sigma'])**2) + return res + + +def logistic_probability(param, event, value): + """ + Logistic probability of a discrete node in Bayesian network with continuous parents, + :param param: a dict, parameters determined by discrete parents of current node + :param event: a dict, names and values of continuous parent variables of current node + :param value: boolean, True or False + :return: int, probability + """ + + buff = 1 + for _,v in event.items(): + # buffer variable to calculate (value-mu)/sigma + + buff *= (v-param['mu'])/param['sigma'] + p = 1 - 1/(1+exp(-4/sqrt(2*pi)*buff)) + return p if value else 1-p + + +class ContinuousBayesNode: + """ A Bayesian network node with continuous distribution or with continuous distributed parents """ + + def __init__(self, name, d_parents, c_parents, parameters, type): + """ + A continuous Bayesian node has two types of parents: discrete and continuous. + :param d_parents: str, name of discrete parents, value of which determines distribution parameters + :param c_parents: str, name of continuous parents, value of which is used to calculate distribution + :param parameters: a dict, parameters for distribution of current node, keys corresponds to discrete parents + :param type: str, type of current node's value, either 'd' (discrete) or 'c'(continuous) + """ + + self.parameters = parameters + self.type = type + self.d_parents = d_parents.split() + self.c_parents = c_parents.split() + self.parents = self.d_parents + self.c_parents + self.variable = name + self.children = [] + + def continuous_p(self, value, c_event, d_event): + """ + Probability given the value of current node and its parents + :param c_event: event of continuous nodes + :param d_event: event of discrete nodes + """ + assert isinstance(c_event, dict) + assert isinstance(d_event, dict) + + d_event_vals = event_values(d_event, self.d_parents) + if len(d_event_vals) == 1: + d_event_vals = d_event_vals[0] + param = self.parameters[d_event_vals] + if self.type == "c": + p = gaussian_probability(param, c_event, value) + if self.type == "d": + p = logistic_probability(param, c_event, value) + return p + +# harvest-buy example. Figure 13.5 + + +harvest_buy = BayesNet([ + ('Subsidy', '', 0.001), + ('Harvest', '', 0.002), + ('Cost', 'Subsidy', 'Harvest', + {True: {'sigma': 0.5, 'b': 1, 'a': {'Harvest': 0.5}}, + False: {'sigma': 0.6, 'b': 1, 'a': {'Harvest': 0.5}}}, 'c'), + ('Buys', '', 'Cost', {T: {'mu':0.5, 'sigma':0.5}, F: {'mu': 0.6, 'sigma':0.6}}, 'd'), +]) + + +# ______________________________________________________________________________ +# 13.3 Exact Inference in Bayesian Networks +# 13.3.1 Inference by enumeration + + +def enumeration_ask(X, e, bn): + """ + Return the conditional probability distribution of variable X + given evidence e, from BayesNet bn. [Figure 13.10] + >>> enumeration_ask('Burglary', dict(JohnCalls=T, MaryCalls=T), burglary + ... ).show_approx() + 'False: 0.716, True: 0.284' + """ + + assert X not in e, "Query variable must be distinct from evidence" + Q = ProbDist(X) + for xi in bn.variable_values(X): + Q[xi] = enumerate_all(bn.variables, extend(e, X, xi), bn) + return Q.normalize() + + +def enumerate_all(variables, e, bn): + """ + Return the sum of those entries in P(variables | e{others}) + consistent with e, where P is the joint distribution represented + by bn, and e{others} means e restricted to bn's other variables + (the ones other than variables). Parents must precede children in variables. + """ + + if not variables: + return 1.0 + Y, rest = variables[0], variables[1:] + Ynode = bn.variable_node(Y) + if Y in e: + return Ynode.p(e[Y], e) * enumerate_all(rest, e, bn) + else: + return sum(Ynode.p(y, e) * enumerate_all(rest, extend(e, Y, y), bn) + for y in bn.variable_values(Y)) + +# ______________________________________________________________________________ +# 13.3.2 The variable elimination algorithm + + +def elimination_ask(X, e, bn): + """ + Compute bn's P(X|e) by variable elimination. [Figure 13.12] + >>> elimination_ask('Burglary', dict(JohnCalls=T, MaryCalls=T), burglary + ... ).show_approx() + 'False: 0.716, True: 0.284' + """ + assert X not in e, "Query variable must be distinct from evidence" + factors = [] + for var in reversed(bn.variables): + factors.append(make_factor(var, e, bn)) + if is_hidden(var, X, e): + factors = sum_out(var, factors, bn) + return pointwise_product(factors, bn).normalize() + + +def is_hidden(var, X, e): + """Is var a hidden variable when querying P(X|e)?""" + return var != X and var not in e + + +def make_factor(var, e, bn): + """ + Return the factor for var in bn's joint distribution given e. + That is, bn's full joint distribution, projected to accord with e, + is the pointwise product of these factors for bn's variables. + """ + node = bn.variable_node(var) + variables = [X for X in [var] + node.parents if X not in e] + cpt = {event_values(e1, variables): node.p(e1[var], e1) + for e1 in all_events(variables, bn, e)} + return Factor(variables, cpt) + + +def pointwise_product(factors, bn): + return reduce(lambda f, g: f.pointwise_product(g, bn), factors) + + +def sum_out(var, factors, bn): + """Eliminate var from all factors by summing over its values.""" + result, var_factors = [], [] + for f in factors: + (var_factors if var in f.variables else result).append(f) + result.append(pointwise_product(var_factors, bn).sum_out(var, bn)) + return result + + +class Factor: + """A factor in a joint distribution.""" + + def __init__(self, variables, cpt): + self.variables = variables + self.cpt = cpt + + def pointwise_product(self, other, bn): + """Multiply two factors, combining their variables.""" + variables = list(set(self.variables) | set(other.variables)) + cpt = {event_values(e, variables): self.p(e) * other.p(e) + for e in all_events(variables, bn, {})} + return Factor(variables, cpt) + + def sum_out(self, var, bn): + """Make a factor eliminating var by summing over its values.""" + variables = [X for X in self.variables if X != var] + cpt = {event_values(e, variables): sum(self.p(extend(e, var, val)) + for val in bn.variable_values(var)) + for e in all_events(variables, bn, {})} + return Factor(variables, cpt) + + def normalize(self): + """Return my probabilities; must be down to one variable.""" + assert len(self.variables) == 1 + return ProbDist(self.variables[0], + {k: v for ((k,), v) in self.cpt.items()}) + + def p(self, e): + """Look up my value tabulated for e.""" + return self.cpt[event_values(e, self.variables)] + + +def all_events(variables, bn, e): + """Yield every way of extending e with values for all variables.""" + if not variables: + yield e + else: + X, rest = variables[0], variables[1:] + for e1 in all_events(rest, bn, e): + for x in bn.variable_values(X): + yield extend(e1, X, x) + +# ______________________________________________________________________________ +# 13.3.4 Clustering algorithms +# [Figure 13.14a]: sprinkler network + + +sprinkler = BayesNet([ + ('Cloudy', '', 0.5), + ('Sprinkler', 'Cloudy', {T: 0.10, F: 0.50}), + ('Rain', 'Cloudy', {T: 0.80, F: 0.20}), + ('WetGrass', 'Sprinkler Rain', + {(T, T): 0.99, (T, F): 0.90, (F, T): 0.90, (F, F): 0.00})]) + +# ______________________________________________________________________________ +# 13.4 Approximate Inference for Bayesian Networks +# 13.4.1 Direct sampling methods + + +def prior_sample(bn): + """ + Randomly sample from bn's full joint distribution. The result + is a {variable: value} dict. [Figure 13.15] + """ + event = {} + for node in bn.nodes: + event[node.variable] = node.sample(event) + return event + +# _________________________________________________________________________ + + +def rejection_sampling(X, e, bn, N=10000): + """ + Estimate the probability distribution of variable X given + evidence e in BayesNet bn, using N samples. [Figure 13.16] + Raises a ZeroDivisionError if all the N samples are rejected, + i.e., inconsistent with e. + >>> random.seed(47) + >>> rejection_sampling('Burglary', dict(JohnCalls=T, MaryCalls=T), + ... burglary, 10000).show_approx() + 'False: 0.7, True: 0.3' + """ + counts = {x: 0 for x in bn.variable_values(X)} # bold N in [Figure 13.16] + for j in range(N): + sample = prior_sample(bn) # boldface x in [Figure 13.16] + if consistent_with(sample, e): + counts[sample[X]] += 1 + return ProbDist(X, counts) + + +def consistent_with(event, evidence): + """Is event consistent with the given evidence?""" + return all(evidence.get(k, v) == v + for k, v in event.items()) + +# _________________________________________________________________________ + + +def likelihood_weighting(X, e, bn, N=10000): + """ + Estimate the probability distribution of variable X given + evidence e in BayesNet bn. [Figure 13.17] + >>> random.seed(1017) + >>> likelihood_weighting('Burglary', dict(JohnCalls=T, MaryCalls=T), + ... burglary, 10000).show_approx() + 'False: 0.702, True: 0.298' + """ + + W = {x: 0 for x in bn.variable_values(X)} + for j in range(N): + sample, weight = weighted_sample(bn, e) # boldface x, w in [Figure 14.15] + W[sample[X]] += weight + return ProbDist(X, W) + + +def weighted_sample(bn, e): + """ + Sample an event from bn that's consistent with the evidence e; + return the event and its weight, the likelihood that the event + accords to the evidence. + """ + + w = 1 + event = dict(e) # boldface x in [Figure 13.17] + for node in bn.nodes: + Xi = node.variable + if Xi in e: + w *= node.p(e[Xi], event) + else: + event[Xi] = node.sample(event) + return event, w + +# _________________________________________________________________________ +# 13.4.2 Inference by Markov chain simulation + + +def gibbs_ask(X, e, bn, N=1000): + """[Figure 13.19]""" + assert X not in e, "Query variable must be distinct from evidence" + counts = {x: 0 for x in bn.variable_values(X)} # bold N in [Figure 14.16] + Z = [var for var in bn.variables if var not in e] + state = dict(e) # boldface x in [Figure 14.16] + for Zi in Z: + state[Zi] = random.choice(bn.variable_values(Zi)) + for j in range(N): + for Zi in Z: + state[Zi] = markov_blanket_sample(Zi, state, bn) + counts[state[X]] += 1 + return ProbDist(X, counts) + + +def markov_blanket_sample(X, e, bn): + """ + Return a sample from P(X | mb) where mb denotes that the + variables in the Markov blanket of X take their values from event + e (which must assign a value to each). The Markov blanket of X is + X's parents, children, and children's parents. + """ + Xnode = bn.variable_node(X) + Q = ProbDist(X) + for xi in bn.variable_values(X): + ei = extend(e, X, xi) + # [Equation 13.12:] + Q[xi] = Xnode.p(xi, e) * product(Yj.p(ei[Yj.variable], ei) + for Yj in Xnode.children) + # (assuming a Boolean variable here) + return probability(Q.normalize()[True]) + +# _________________________________________________________________________ +# 13.4.3 Compiling approximate inference + + +class complied_burglary: + """compiled version of burglary network""" + + def Burglary(self, sample): + if sample['Alarm']: + if sample['Earthquake']: + return probability(0.00327) + else: + return probability(0.485) + else: + if sample['Earthquake']: + return probability(7.05e-05) + else: + return probability(6.01e-05) + + def Earthquake(self, sample): + if sample['Alarm']: + if sample['Burglary']: + return probability(0.0020212) + else: + return probability(0.36755) + else: + if sample['Burglary']: + return probability(0.0016672) + else: + return probability(0.0014222) + + def MaryCalls(self, sample): + if sample['Alarm']: + return probability(0.7) + else: + return probability(0.01) + + def JongCalls(self, sample): + if sample['Alarm']: + return probability(0.9) + else: + return probability(0.05) + + def Alarm(self, sample): + raise NotImplementedError diff --git a/tests/test_probability4e.py b/tests/test_probability4e.py new file mode 100644 index 000000000..1ce4d7660 --- /dev/null +++ b/tests/test_probability4e.py @@ -0,0 +1,342 @@ +from probability4e import * + + +def tests(): + cpt = burglary.variable_node('Alarm') + event = {'Burglary': True, 'Earthquake': True} + assert cpt.p(True, event) == 0.95 + event = {'Burglary': False, 'Earthquake': True} + assert cpt.p(False, event) == 0.71 + # #enumeration_ask('Earthquake', {}, burglary) + + s = {'A': True, 'B': False, 'C': True, 'D': False} + assert consistent_with(s, {}) + assert consistent_with(s, s) + assert not consistent_with(s, {'A': False}) + assert not consistent_with(s, {'D': True}) + + random.seed(21) + p = rejection_sampling('Earthquake', {}, burglary, 1000) + assert p[True], p[False] == (0.001, 0.999) + + random.seed(71) + p = likelihood_weighting('Earthquake', {}, burglary, 1000) + assert p[True], p[False] == (0.002, 0.998) + +# test ProbDist + + +def test_probdist_basic(): + P = ProbDist('Flip') + P['H'], P['T'] = 0.25, 0.75 + assert P['H'] == 0.25 + assert P['T'] == 0.75 + assert P['X'] == 0.00 + + P = ProbDist('BiasedDie') + P['1'], P['2'], P['3'], P['4'], P['5'], P['6'] = 10, 15, 25, 30, 40, 80 + P.normalize() + assert P['2'] == 0.075 + assert P['4'] == 0.15 + assert P['6'] == 0.4 + + +def test_probdist_frequency(): + P = ProbDist('X', {'lo': 125, 'med': 375, 'hi': 500}) + assert (P['lo'], P['med'], P['hi']) == (0.125, 0.375, 0.5) + + P = ProbDist('Pascal-5', {'x1': 1, 'x2': 5, 'x3': 10, 'x4': 10, 'x5': 5, 'x6': 1}) + assert (P['x1'], P['x2'], P['x3'], P['x4'], P['x5'], P['x6']) == ( + 0.03125, 0.15625, 0.3125, 0.3125, 0.15625, 0.03125) + + +def test_probdist_normalize(): + P = ProbDist('Flip') + P['H'], P['T'] = 35, 65 + P = P.normalize() + assert (P.prob['H'], P.prob['T']) == (0.350, 0.650) + + P = ProbDist('BiasedDie') + P['1'], P['2'], P['3'], P['4'], P['5'], P['6'] = 10, 15, 25, 30, 40, 80 + P = P.normalize() + assert (P.prob['1'], P.prob['2'], P.prob['3'], P.prob['4'], P.prob['5'], P.prob['6']) == ( + 0.05, 0.075, 0.125, 0.15, 0.2, 0.4) + +# test JoinProbDist + + +def test_jointprob(): + P = JointProbDist(['X', 'Y']) + P[1, 1] = 0.25 + assert P[1, 1] == 0.25 + P[dict(X=0, Y=1)] = 0.5 + assert P[dict(X=0, Y=1)] == 0.5 + + +def test_event_values(): + assert event_values({'A': 10, 'B': 9, 'C': 8}, ['C', 'A']) == (8, 10) + assert event_values((1, 2), ['C', 'A']) == (1, 2) + + +def test_enumerate_joint(): + P = JointProbDist(['X', 'Y']) + P[0, 0] = 0.25 + P[0, 1] = 0.5 + P[1, 1] = P[2, 1] = 0.125 + assert enumerate_joint(['Y'], dict(X=0), P) == 0.75 + assert enumerate_joint(['X'], dict(Y=2), P) == 0 + assert enumerate_joint(['X'], dict(Y=1), P) == 0.75 + + Q = JointProbDist(['W', 'X', 'Y', 'Z']) + Q[0, 1, 1, 0] = 0.12 + Q[1, 0, 1, 1] = 0.4 + Q[0, 0, 1, 1] = 0.5 + Q[0, 0, 1, 0] = 0.05 + Q[0, 0, 0, 0] = 0.675 + Q[1, 1, 1, 0] = 0.3 + assert enumerate_joint(['W'], dict(X=0, Y=0, Z=1), Q) == 0 + assert enumerate_joint(['W'], dict(X=0, Y=0, Z=0), Q) == 0.675 + assert enumerate_joint(['W'], dict(X=0, Y=1, Z=1), Q) == 0.9 + assert enumerate_joint(['Y'], dict(W=1, X=0, Z=1), Q) == 0.4 + assert enumerate_joint(['Z'], dict(W=0, X=0, Y=0), Q) == 0.675 + assert enumerate_joint(['Z'], dict(W=1, X=1, Y=1), Q) == 0.3 + + +def test_enumerate_joint_ask(): + P = JointProbDist(['X', 'Y']) + P[0, 0] = 0.25 + P[0, 1] = 0.5 + P[1, 1] = P[2, 1] = 0.125 + assert enumerate_joint_ask( + 'X', dict(Y=1), P).show_approx() == '0: 0.667, 1: 0.167, 2: 0.167' + + +def test_is_independent(): + P = JointProbDist(['X', 'Y']) + P[0, 0] = P[0,1] = P[1, 1] = P[1, 0] = 0.25 + assert enumerate_joint_ask( + 'X', dict(Y=1), P).show_approx() == '0: 0.5, 1: 0.5' + assert is_independent(['X','Y'], P) + +# test BayesNode + + +def test_bayesnode_p(): + bn = BayesNode('X', 'Burglary', {T: 0.2, F: 0.625}) + assert bn.p(True, {'Burglary': True, 'Earthquake': False}) == 0.2 + assert bn.p(False, {'Burglary': False, 'Earthquake': True}) == 0.375 + assert BayesNode('W', '', 0.75).p(False, {'Random': True}) == 0.25 + + +def test_bayesnode_sample(): + X = BayesNode('X', 'Burglary', {T: 0.2, F: 0.625}) + assert X.sample({'Burglary': False, 'Earthquake': True}) in [True, False] + Z = BayesNode('Z', 'P Q', {(True, True): 0.2, (True, False): 0.3, + (False, True): 0.5, (False, False): 0.7}) + assert Z.sample({'P': True, 'Q': False}) in [True, False] + +# test continuous variable bayesian net + + +def test_gaussian_probability(): + param = {'sigma': 0.5, 'b': 1, 'a': {'h': 0.5}} + event = {'h': 0.6} + assert gaussian_probability(param, event, 1) == 0.6664492057835993 + + +def test_logistic_probability(): + param = {'mu': 0.5, 'sigma': 0.1} + event = {'h': 0.6} + assert logistic_probability(param, event, True) == 0.16857376940725355 + assert logistic_probability(param, event, False) == 0.8314262305927465 + + +def test_enumeration_ask(): + assert enumeration_ask( + 'Burglary', dict(JohnCalls=T, MaryCalls=T), + burglary).show_approx() == 'False: 0.716, True: 0.284' + assert enumeration_ask( + 'Burglary', dict(JohnCalls=T, MaryCalls=F), + burglary).show_approx() == 'False: 0.995, True: 0.00513' + assert enumeration_ask( + 'Burglary', dict(JohnCalls=F, MaryCalls=T), + burglary).show_approx() == 'False: 0.993, True: 0.00688' + assert enumeration_ask( + 'Burglary', dict(JohnCalls=T), + burglary).show_approx() == 'False: 0.984, True: 0.0163' + assert enumeration_ask( + 'Burglary', dict(MaryCalls=T), + burglary).show_approx() == 'False: 0.944, True: 0.0561' + + +def test_elimination_ask(): + assert elimination_ask( + 'Burglary', dict(JohnCalls=T, MaryCalls=T), + burglary).show_approx() == 'False: 0.716, True: 0.284' + assert elimination_ask( + 'Burglary', dict(JohnCalls=T, MaryCalls=F), + burglary).show_approx() == 'False: 0.995, True: 0.00513' + assert elimination_ask( + 'Burglary', dict(JohnCalls=F, MaryCalls=T), + burglary).show_approx() == 'False: 0.993, True: 0.00688' + assert elimination_ask( + 'Burglary', dict(JohnCalls=T), + burglary).show_approx() == 'False: 0.984, True: 0.0163' + assert elimination_ask( + 'Burglary', dict(MaryCalls=T), + burglary).show_approx() == 'False: 0.944, True: 0.0561' + + +# test sampling + + +def test_prior_sample(): + random.seed(42) + all_obs = [prior_sample(burglary) for x in range(1000)] + john_calls_true = [observation for observation in all_obs if observation['JohnCalls'] == True] + mary_calls_true = [observation for observation in all_obs if observation['MaryCalls'] == True] + burglary_and_john = [observation for observation in john_calls_true if observation['Burglary'] == True] + burglary_and_mary = [observation for observation in mary_calls_true if observation['Burglary'] == True] + assert len(john_calls_true) / 1000 == 46 / 1000 + assert len(mary_calls_true) / 1000 == 13 / 1000 + assert len(burglary_and_john) / len(john_calls_true) == 1 / 46 + assert len(burglary_and_mary) / len(mary_calls_true) == 1 / 13 + + +def test_prior_sample2(): + random.seed(128) + all_obs = [prior_sample(sprinkler) for x in range(1000)] + rain_true = [observation for observation in all_obs if observation['Rain'] == True] + sprinkler_true = [observation for observation in all_obs if observation['Sprinkler'] == True] + rain_and_cloudy = [observation for observation in rain_true if observation['Cloudy'] == True] + sprinkler_and_cloudy = [observation for observation in sprinkler_true if observation['Cloudy'] == True] + assert len(rain_true) / 1000 == 0.476 + assert len(sprinkler_true) / 1000 == 0.291 + assert len(rain_and_cloudy) / len(rain_true) == 376 / 476 + assert len(sprinkler_and_cloudy) / len(sprinkler_true) == 39 / 291 + + +def test_rejection_sampling(): + random.seed(47) + assert rejection_sampling( + 'Burglary', dict(JohnCalls=T, MaryCalls=T), + burglary, 10000).show_approx() == 'False: 0.7, True: 0.3' + assert rejection_sampling( + 'Burglary', dict(JohnCalls=T, MaryCalls=F), + burglary, 10000).show_approx() == 'False: 1, True: 0' + assert rejection_sampling( + 'Burglary', dict(JohnCalls=F, MaryCalls=T), + burglary, 10000).show_approx() == 'False: 0.987, True: 0.0128' + assert rejection_sampling( + 'Burglary', dict(JohnCalls=T), + burglary, 10000).show_approx() == 'False: 0.982, True: 0.0183' + assert rejection_sampling( + 'Burglary', dict(MaryCalls=T), + burglary, 10000).show_approx() == 'False: 0.965, True: 0.0348' + + +def test_rejection_sampling2(): + random.seed(42) + assert rejection_sampling( + 'Cloudy', dict(Rain=T, Sprinkler=T), + sprinkler, 10000).show_approx() == 'False: 0.56, True: 0.44' + assert rejection_sampling( + 'Cloudy', dict(Rain=T, Sprinkler=F), + sprinkler, 10000).show_approx() == 'False: 0.119, True: 0.881' + assert rejection_sampling( + 'Cloudy', dict(Rain=F, Sprinkler=T), + sprinkler, 10000).show_approx() == 'False: 0.951, True: 0.049' + assert rejection_sampling( + 'Cloudy', dict(Rain=T), + sprinkler, 10000).show_approx() == 'False: 0.205, True: 0.795' + assert rejection_sampling( + 'Cloudy', dict(Sprinkler=T), + sprinkler, 10000).show_approx() == 'False: 0.835, True: 0.165' + + +def test_likelihood_weighting(): + random.seed(1017) + assert likelihood_weighting( + 'Burglary', dict(JohnCalls=T, MaryCalls=T), + burglary, 10000).show_approx() == 'False: 0.702, True: 0.298' + assert likelihood_weighting( + 'Burglary', dict(JohnCalls=T, MaryCalls=F), + burglary, 10000).show_approx() == 'False: 0.993, True: 0.00656' + assert likelihood_weighting( + 'Burglary', dict(JohnCalls=F, MaryCalls=T), + burglary, 10000).show_approx() == 'False: 0.996, True: 0.00363' + assert likelihood_weighting( + 'Burglary', dict(JohnCalls=F, MaryCalls=F), + burglary, 10000).show_approx() == 'False: 1, True: 0.000126' + assert likelihood_weighting( + 'Burglary', dict(JohnCalls=T), + burglary, 10000).show_approx() == 'False: 0.979, True: 0.0205' + assert likelihood_weighting( + 'Burglary', dict(MaryCalls=T), + burglary, 10000).show_approx() == 'False: 0.94, True: 0.0601' + + +def test_likelihood_weighting2(): + random.seed(42) + assert likelihood_weighting( + 'Cloudy', dict(Rain=T, Sprinkler=T), + sprinkler, 10000).show_approx() == 'False: 0.559, True: 0.441' + assert likelihood_weighting( + 'Cloudy', dict(Rain=T, Sprinkler=F), + sprinkler, 10000).show_approx() == 'False: 0.12, True: 0.88' + assert likelihood_weighting( + 'Cloudy', dict(Rain=F, Sprinkler=T), + sprinkler, 10000).show_approx() == 'False: 0.951, True: 0.0486' + assert likelihood_weighting( + 'Cloudy', dict(Rain=T), + sprinkler, 10000).show_approx() == 'False: 0.198, True: 0.802' + assert likelihood_weighting( + 'Cloudy', dict(Sprinkler=T), + sprinkler, 10000).show_approx() == 'False: 0.833, True: 0.167' + + +def test_gibbs_ask(): + + g_solution = gibbs_ask('Cloudy', dict(Rain=True), sprinkler, 1000) + assert abs(g_solution.prob[False]-0.2) < 0.05 + assert abs(g_solution.prob[True]-0.8) < 0.05 + + +# The following should probably go in .ipynb: + +""" +# We can build up a probability distribution like this (p. 469): +>>> P = ProbDist() +>>> P['sunny'] = 0.7 +>>> P['rain'] = 0.2 +>>> P['cloudy'] = 0.08 +>>> P['snow'] = 0.02 + +# and query it like this: (Never mind this ELLIPSIS option +# added to make the doctest portable.) +>>> P['rain'] #doctest:+ELLIPSIS +0.2... + +# A Joint Probability Distribution is dealt with like this [Figure 13.3]: +>>> P = JointProbDist(['Toothache', 'Cavity', 'Catch']) +>>> T, F = True, False +>>> P[T, T, T] = 0.108; P[T, T, F] = 0.012; P[F, T, T] = 0.072; P[F, T, F] = 0.008 +>>> P[T, F, T] = 0.016; P[T, F, F] = 0.064; P[F, F, T] = 0.144; P[F, F, F] = 0.576 + +>>> P[T, T, T] +0.108 + +# Ask for P(Cavity|Toothache=T) +>>> PC = enumerate_joint_ask('Cavity', {'Toothache': T}, P) +>>> PC.show_approx() +'False: 0.4, True: 0.6' + +>>> 0.6-epsilon < PC[T] < 0.6+epsilon +True + +>>> 0.4-epsilon < PC[F] < 0.4+epsilon +True +""" + +if __name__ == '__main__': + pytest.main()