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 "<Node {}>".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')
+        <DataSet(): 1 examples, 3 attributes>
+        """
+        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 '<DataSet({}): {:d} examples, {:d} attributes>'.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 '<Grammar {}>'.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 '<Grammar {}>'.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)]