p07.py

from abc import ABC
import matplotlib.pyplot as plt
import numpy as np


def check_dims(slopes):
    def wrapper(self):
        result = slopes(self)
        for i, r in enumerate(result):
            assert r.shape[0] == len(self.output), \
                f'{self.__class__.__name__}(' \
                f'{",".join([str(x.output) for x in self.inputs])}' \
                f').slopes[{i}].shape != ' \
                f'({len(self.output)}, {len(self.inputs[i].output)})'
        return result
    return wrapper


class Operand:
    def __add__(self, other):
        return Add(self, other)

    def __sub__(self, other):
        return Subtract(self, other)

    def __mul__(self, other):
        return Multiply(self, other)

    def __pow__(self, power):
        return Power(self, power=power)

    def __truediv__(self, other):
        return Divide(self, other)

    def __neg__(self):
        return Negate(self)


class Node(Operand):
    """
    Main node class with multiple inputs, 1 output, and "toString" method.
    Subclasses must implement eval method.
    """
    def __init__(self, *inputs):
        self.inputs = inputs
        self.output = None
        self.name = self.__class__.__name__

    def __call__(self):
        self.output = self.eval()
        return self.output

    def __str__(self):
        child_names = [str(node) for node in self.inputs]
        return f'{self.name}({",".join(child_names)})'

    def __repr__(self):
        return self.__str__()

    def outputs(self):
        child_outputs = [node.outputs() for node in self.inputs]
        return f'{self.output}{self.name}({",".join(child_outputs)})'

    def find(self, criterion):
        """
        Recursively search the graph to find the node with the name
        :param criterion: a function that returns True for those nodes to return.
        :return: list of nodes that satisfy the criterion.
        """
        items = set()
        if criterion(self):
            items.add(self)
        for child in self.inputs:
            items |= child.find(criterion)
        return items

    def eval(self):
        """
        Computes the input values from its children and uses them to compute
        and return this node's output.
        :return: the output value for this node.
        """
        raise NotImplementedError()

    def __getitem__(self, item):
        return self.inputs[item]


class Differentiable(Node, ABC):
    def __init__(self, *inputs):
        super().__init__(*inputs)
        self.output_slope = None

    #!!!TAKEN OUT AS PART OF P7!!!#
    # def __call__(self):
    #     result = super().__call__()
    #     self.output_slope = np.zeros_like(result)
    #     return result

    def output_slopes(self):
        """
        Return a string representation of the node and its children showing
        their slopes.
        :return: the string
        """
        child_slopes = [node.output_slopes() for node in self.inputs]
        return f'{self.output_slope}{self.name}({",".join(child_slopes)})'

    def update(self, eta):
        """
        Update all descendents. Parameters should be leaf nodes that update
        themselves differently.
        :param eta: Learning rate
        :return: None
        """
        for node in self.inputs:
            node.update(eta)

    def backprop(self, output_slope=None):
        """
        Compute and propagate the input slope to the descendents.

        :param output_slope: The output slope computed from ancestors.
        """
        if output_slope is None:
            self()  # call the graph and reset output_slopes to zero
            output_slope = np.ones_like(self.output)
        if self.output_slope is None:
            self.output_slope = np.zeros_like(self.output)
        #cast to a numpy array
        output_slope = np.array(output_slope, dtype=float)
        #accumulate output slopes
        self.output_slope += output_slope
        for curr, slope in zip(self.inputs, self.slopes()):
            curr.backprop(np.matmul(output_slope, slope))


    def slopes(self):
        """
        Compute the derivative of the output w.r.t. each input. Use each child's
        stored output rather than evaluating from scratch.
        :return: A list of Jacobian matrices, one for each input.
        """
        raise NotImplementedError()


class Variable(Differentiable):
    """
    Make a Variable node that has no inputs and evaluates to its 'output'
    value. Include a 'set' method to change its value.
    """
    def __init__(self, value, name):
        super().__init__()
        self.output = np.array(value, dtype=float)
        assert self.output.ndim == 1
        self.name = name

    def __str__(self):
        return f'{self.name}={self.output}'

    def eval(self):
        return self.output

    def slopes(self):
        return ()

    def set(self, value):
        self.output = np.array(value, dtype=float)
        assert self.output.ndim == 1

    def get(self):
        return self.output


class Parameter(Variable):
    def __init__(self, value, name, l1=0., l2=0.):
        super().__init__(value, name)
        self.l1 = l1
        self.l2 = l2
    
    def penalty(self):
        """
        Return a penalty node that computes the L1 and L2 regularization penalty for this parameter.
        :return: the root of a computational graph that computes the penalty.
        """
        penalty = (Constant([self.l1]) * Sum(Abs(self))) + (Constant([self.l2]) * Sum(self ** 2))
        return penalty
        
    
    def update(self, eta):
        #check if the slope is None to ensure 1 update per node.
        if self.output_slope is not None:   
        #that derivative is stored in self.output_slope and calculated using backprop.
            self.output -= eta * self.output_slope
            self.output_slope = np.zeros_like(self.output)


class Constant(Variable):
    def __init__(self, value):
        super().__init__(value, str(np.array(value)))

    def __str__(self):
        return f'{self.name}'

    def set(self, value):
        raise NotImplementedError("Can't set the value of a Constant")


class UnaryOperator(Differentiable, ABC):
    def __init__(self, *args):
        super().__init__(*args)
        self.a, = args


class BinaryOperator(Differentiable, ABC):
    def __init__(self, *args):
        super().__init__(*args)
        self.a, self.b = args


class Add(BinaryOperator):
    def eval(self):
        a = self.a()
        b = self.b()
        return a + b

    @check_dims
    def slopes(self):
        a = self.a.output
        b = self.b.output
        #matrix dimensions are 2*len(self.a) by len(self.a)
        #all as and bs lined up horizontally where each row is y1, y2, y3
        return [np.eye(len(a)), np.eye(len(b))]

class Subtract(BinaryOperator):
    def eval(self):
        a = self.a()
        b = self.b()
        return a - b

    @check_dims
    def slopes(self):
        a = self.a.output
        b = self.b.output
        return [np.eye(len(a)), -np.eye(len(b))]

class Multiply(BinaryOperator):
    def eval(self):
        a = self.a()
        b = self.b()
        return a * b

    @check_dims
    def slopes(self):
        a = self.a.output
        b = self.b.output
        return [np.eye(len(a)) * b, np.eye(len(b)) * a]

class Divide(BinaryOperator):
    def eval(self):
        a = self.a()
        b = self.b()
        return a / b
    
    @check_dims
    def slopes(self):
        a = self.a.output
        b = self.b.output
        jacobian1 = np.eye(len(b)) * 1/ b 
        jacobian2 = np.eye(len(a)) * (-a / (b ** 2))
        jacobian1[np.isnan(jacobian1)] = 0
        jacobian2[np.isnan(jacobian2)] = 0
        return [jacobian1, jacobian2]

class Negate(UnaryOperator):
    def eval(self):
        a = self.a()
        return -a

    @check_dims
    def slopes(self):
        a = self.a.output
        return [-np.eye(len(a))]

class Power(UnaryOperator):
    def __init__(self, *inputs, power=1):
        super().__init__(*inputs)
        self.p = power
        self.name += str(power)

    def eval(self):
        a = self.a()
        return a ** self.p

    @check_dims
    def slopes(self):
        a = self.a.output
        return [np.eye(len(a)) * self.p * (a ** (self.p-1))]


class Sigmoid(UnaryOperator):          
    def eval(self):
        a = self.a()
        idx = a >= 0
        result = np.zeros_like(a)
        result[idx] = 1 / (1 + np.exp(-a[idx]))
        result[~idx] = np.exp(a[~idx]) / (1 + np.exp(a[~idx]))
        return result
    
    @check_dims
    def slopes(self):
        a = self.a.output
        b = np.where(a >= 0, np.eye(len(a)) * (np.exp(-a) / ((1 + np.exp(-a)) ** 2)), np.eye(len(a)) * (np.exp(a) / ((1 + np.exp(a)) ** 2)))

        return [np.eye(len(a)) * b]


class Log(UnaryOperator):
    def eval(self):
        a = self.a()
        idx = a <= 0
        #if idx, return 0, else, np.log(x)
        result = np.zeros_like(a)
        result[idx] = 0
        result[~idx] = np.log(a[~idx])
        return result
        
    @check_dims
    def slopes(self):
        a = self.a.output
        b = np.where(a <= 0, 0, 1/a)
        return [np.eye(len(a)) * b]
        
        
class Sum(UnaryOperator):
    def eval(self):
        a = self.a()
        return np.array([np.sum(a)])
        
    @check_dims
    def slopes(self):
        a = self.a.output
        return [np.ones((1, len(a)))]


class Abs(UnaryOperator):
    def eval(self):
        a = self.a()
        return np.array(np.absolute(a))

    @check_dims
    def slopes(self):
        a = self.a.output
        return [np.eye(len(a)) * np.sign(a)]
    
    
class Dot(BinaryOperator):
    def eval(self):
        a = self.a()
        b = self.b()
        return np.array([np.dot(a, b)])

    @check_dims
    def slopes(self):
        a = self.a.output
        b = self.b.output
        return [np.array([b]), np.array([a])]


class PolynomialFeatures(UnaryOperator):
    def __init__(self, *inputs, degree=2):
        """
        Only works for scalar inputs.
        :param inputs: list containing one input
        :param degree: the degree of the polynomial
        """
        super().__init__(*inputs)
        self.degree = degree

    def eval(self):
        a = self.a()
        res = np.array([a[0]**p for p in range(1, self.degree+1)])
        return res
    
    def slopes(self):
        a = self.a.output
        list = [n*a[0]**(n-1)for n in range (1,self.degree+1)]
        return [np.array(list).reshape((self.degree, 1))]

def stochastic_gradient_descent(cost, x, y, data, num_iter, eta):
    """
    Adapt the parameters to reduce the cost using SGD
    :param cost: the root node of the cost function
    :param x: the x-node
    :param y: the y-node
    :param data: a list of (x, y) vector (list) pairs
    :param num_iter: the number of epochs to process
    :param eta: the learning rate
    :return: list of costs (floats)
    """
    np.random.seed(1)
    costs = []
    for i in range(num_iter):
        np.random.shuffle(data)
        c = np.zeros_like(data[0][1])
        for xi, yi in data:
            x.set(np.array(xi))
            y.set(np.array(yi))
            c = c + cost.backprop()
            cost.update(eta)

        c /= len(data)
        costs.append(c)
        print(f'{i}: {c}')
    return costs


#p07
def batch_update(cost, x, y, data, eta):
    """
    A batch update takes one step of gradient descent using all data passed in as an input argument.
    The output_slope should be set to zero and then accumulated for all samples, culminating with
    one update step with the learning rate eta.

    Args:
        cost (Node): The root node of the 'cost' graph
        x (Node):  The input 'x' node
        y (Node): The output 'y' node
        data ((vector, vector)): : A batch of samples (a list of (x, y) vector (list) pairs)
        eta (float): Learning rate
    
    Returns:
        float: the average cost for the samples in the batch
    """
    c=0.0
    #1. For each sample (x,y) in data:
    for (sx, sy) in data:
        #1.1 set the value of x and y in the graph
        x.set(np.array(sx))
        y.set(np.array(sy))
        #1.2 backprop and accumulate the gradient
        c += cost()
        cost.backprop()
    #2. update each parameter by following the gradient with the learning rate eta.
        #theta = theta - eta * gradient(theta)
    cost.update(eta)
    #3. return the average cost per sample
    return c / len(data)


def minibatch_gradient_descent(cost, x, y, data, batch_size=1, num_iter=25, 
                               eta=0.01):
    """
    Adapt the parameters to reduce the cost using minibatch gradient descent.

    Args:
        cost (Node): the root node of the cost function
        x (Node):  the x-node
        y (Node): the y-node
        data ((vector, vector)): a list of (x, y) vector (list) pairs
        batch_size (int, optional): number of samples per batch. Defaults to 1.
        num_iter (int, optional): the number of epochs to process. Defaults to 25.
        eta (float, optional): the learning rate. Defaults to 0.01.

    Returns:
        [float]: list of costs (floats) that are the average cost per minibatch.
    """
    np.random.seed(1)
    costs = []
    for epoch in range(num_iter):
        np.random.shuffle(data)
        i = 0
        sublist = data[i:i+batch_size]
        while len(sublist) > 0:
            costs.append(batch_update(cost, x, y, sublist, eta))
            i += batch_size
            sublist = data[i:i+batch_size] 
        
        #print the total average cost per minibatch for this epoch
        print(f'{epoch}: {costs}')
    return costs

class Linear(UnaryOperator):
    #linear activation function is just a pass through
    def eval(self):
        a = self.a()
        return a
    def slopes(self):
        a = self.a.output
        return [np.eye(len(a))]       
        
class Tanh(UnaryOperator):
    def eval(self):
        a = self.a()
        return np.tanh(a)
    def slopes(self):
        a = self.a.output
        return [np.diagflat(1 - (np.tanh(a) ** 2))]
        
        
class Concatenate(Differentiable):
    #combines outputs of multiple neurons.
    def eval(self):
        inputs = [x() for x in self.inputs]
        return np.concatenate(inputs)
    
    def slopes(self):
        inputs = [x.output for x in self.inputs]
        num_outputs = len(self.output)
        slopes = []
        k = 0
        
        for x in inputs:
            num_inputs = len(x)
            jacobian = np.zeros((num_outputs, num_inputs))
            jacobian[k:k+num_inputs, :] = np.eye(num_inputs)
            slopes.append(jacobian)
            k += num_inputs
        return slopes   
    
         
class Neuron(Linear):
    def __init__(self, x, l1=0., l2=0., activation=UnaryOperator):
        """
        Evaluates activation(w.x + b)

        Args:
            x (_type_): the input to the neuron
            l1 (_type_, optional): l1 regularization coefficient for the weights. Defaults to 0..
            l2 (_type_, optional): l2 regularization coefficient for the bias. Defaults to 0..
            activation (_type_, optional): the activation function. Defaults to UnaryOperator.
        """
        
        # initialize w parameter to normal(0, 0.1) elements
        # initialize b to vector of zero
        # construct graph for a = phi(w.x + b)
        # inputs should be 'a'
        # because Neuron "is a" Linear, the output of 'a' is passed through as well as the slopes.
        # this creates an extra linear node at the output of the activation function for this neuron
        num_inputs = len(x.eval())
        np.random.seed(1)
        w = Parameter(np.array(np.random.normal(0,0.1,num_inputs)), 'w', l1, l2) #need to regularize these weights.
        b = Parameter([0], 'b') #not regularized
        a = activation(Dot(w,x) + b)
        super().__init__(a)
        #output of a neuron is y_hat
        
class Model:
    def __init__(self, num_inputs, num_outputs, **kwargs):
        """
        Construct a model with input 'x' and output 'y_hat', construct its
        loss graph and penalty graph.

        Args:
            num_inputs (int): number of inputs (length of 'x')
            num_outputs (int):  number of outputs (length of 'y')
        """
        #set fields: x, y, y_hat, and cost = loss + penalty
        self.x, self.y_hat = self.build_model(num_inputs, num_outputs, **kwargs)
        self.y = Variable(np.zeros(num_outputs), 'y')
        self.cost = self.loss(self.y_hat, self.y) + self.penalty()
    
    #abstract functions
    def build_model(self, num_inputs, num_outputs, **kwargs): raise NotImplementedError()
    def loss(self, y_hat, y): raise NotImplementedError()
    
    def params(self):
        """
        Returns:
            Node:[Parameters] -> list of parameters
        """
        return self.y_hat.find(lambda node: isinstance(node, Parameter))
    
    def penalty(self):
        """
        Returns:
            Node: returns the root node of a graph that computes the sum of penalties for 
            the parameters in this model. 
        """
        p = self.params()
        penalty = Constant([0])
        for param in p:
            penalty += param.penalty()
        return penalty
        
    
    def fit (self, data, batch_size=1, num_iter=100, eta=0.01):
        """
        Uses minibatch_gradient_descent to optimize the model and return
        a list of costs (floats).
        """
        costs = minibatch_gradient_descent(self.cost, self.x, self.y, data, batch_size, num_iter, eta)
        return costs
    
    
class LinearRegression(Model):
    def build_model(self, num_inputs, num_outputs, l1=0., l2=0.,):
        #construct the model with input 'x' and output 'y_hat'
        #create a neuron with linear with activation
        x = Variable(np.zeros(num_inputs), 'x')
        
        y_hat = Neuron(x, l1, l2, Linear)
        return x, y_hat
    
    def loss (self, y_hat, y):
        #return the root node of the graph that computes the loss function.add()
        #return loss, activation is sigmoid. only difference is activation function.
        loss = (y_hat - y) ** 2
        return loss


class LogisticRegression(Model):
    def build_model(self, num_inputs, num_outputs, l1=0., l2=0.,):
        #construct the model with input 'x' and output 'y_hat'
        x=Variable(np.zeros(num_inputs), 'x')
        y_hat = Neuron(x, l1, l2, Sigmoid)
        return x, y_hat
    
    def loss (self, y_hat, y):
        #return the root node of the graph that computes the loss function.add()
        loss = -y * Log(y_hat) - (Constant([1]) - y) * Log(Constant([1]) - y_hat)
        return loss
    
    
class NonlinearClassifier(Model):
    def build_model(self, num_inputs, num_outputs, l1=0., l2=0., num_neurons=1,  activation=Differentiable):
        x = Variable(np.zeros(num_inputs), 'x')
        
        neurons = []
        for i in range(num_neurons):
           neurons.append(Neuron(x, l1, l2, activation))
        
        y_hat = Neuron(Concatenate(*neurons), l1, l2, Sigmoid)
        return x, y_hat
    
    def loss(self, y_hat, y):
        #return the root node of the graph that computes the loss function.
        loss = -y * Log(y_hat) - (Constant([1]) - y) * Log(Constant([1]) - y_hat)
        return loss

def plot_output_1d(x, output, data, steps=101):
    old_x = x.get()

    x_data, y_data = zip(*data)
    x_data = np.array(x_data)
    y_data = np.array(y_data)

    x_min = x_data.min(axis=0, initial=None)
    x_max = x_data.max(axis=0, initial=None)

    xx = np.linspace(x_min, x_max, steps)

    yy = np.zeros(steps)
    for i in range(len(xx)):
        x.set(xx[i])
        yy[i] = output()

    x.set(old_x)

    fig, ax = plt.subplots()
    ax.plot(x_data, y_data, linestyle='None', marker='o', color='blue')
    ax.plot(xx, yy, color='green')
    ax.set(title=f'{output.name}')
    ax.set(ylabel=f'{output.name}', xlabel='x')


def plot_output_2d(x, output, data, steps=101, quiver=False):
    old_x = x.get()

    x_data, y_data = zip(*data)
    x_data = np.array(x_data)
    y_data = np.array(y_data)

    x1_min, x2_min = x_data.min(axis=0, initial=None)
    x1_max, x2_max = x_data.max(axis=0, initial=None)

    x1 = np.linspace(x1_min, x1_max, steps)
    x2 = np.linspace(x2_min, x2_max, steps)

    num_features = x_data.shape[1]

    yy = np.zeros((steps, steps))
    yy_slope = None
    if quiver:
        yy_slope = np.zeros((steps, steps, num_features))

    for i in range(len(x2)):
        for j in range(len(x1)):
            x.set([x1[j], x2[i]])
            y_hat = output()
            yy[i, j] = y_hat
            if quiver:
                output.backprop()
                yy_slope[i, j] = x.output_slope

    x.set(old_x)

    fig, ax = plt.subplots()
    im = ax.imshow(yy,
                   extent=(x1_min, x1_max, x2_min, x2_max),
                   origin='lower')
    ax.scatter(x_data[:, 0], x_data[:, 1], c=y_data,
               edgecolor='black')
    ax.set(title=f'{output.name}')
    ax.set(ylabel='x2', xlabel='x1')

    if quiver:
        ax.quiver(x1, x2, yy_slope[:, :, 0], yy_slope[:, :, 1])

    fig.subplots_adjust(right=0.8)
    cax = fig.add_axes([0.85, 0.15, 0.05, 0.7])
    fig.colorbar(im, cax=cax)

#TESTING HELPERS FOR P07
def example_one(w, b, name):
    np.random.seed(1)
    n = 1000
    m = 2
    x = np.random.uniform(-1, 1, (n, m))
    z = x.dot(w) + b
    y = (np.sign(z) + 1) // 2
    data = list(zip(x.tolist(), y.tolist()))
    model = LogisticRegression(m, 1)
    c = model.fit(data, batch_size=1, num_iter=25, eta=0.01)
    print(model.params())
    
    #plot output
    plot_output_2d(model.x, model.y_hat, data, quiver=False, steps=101)
    
    fig, ax = plt.subplots()
    ax.plot(c)
    ax.set(xlabel='epoch', ylabel='cost', yscale='linear',
           title=f'{name.title()}')
    plt.show()


def example_two(w, b, u, c, name):
    np.random.seed(1)
    n = 1000
    m = 2
    x = np.random.uniform(-1, 1, (n, m))
    z1 = x.dot(w) + b
    a1 = (np.sign(z1) + 1) // 2
    z2 = a1.dot(u) + c
    y = (np.sign(z2) + 1) // 2
    data = list(zip(x.tolist(), y.tolist()))
    model = NonlinearClassifier(m, 1, num_neurons=2, activation=Sigmoid)
    c = model.fit(data, batch_size=32, num_iter=25, eta=0.1)
    print(model.params())
    
    #plot output
    plot_output_2d(model.x, model.y_hat, data, quiver=False, steps=101)
    
    fig, ax = plt.subplots()
    ax.plot(c)
    ax.set(xlabel='epoch', ylabel='cost', yscale='linear',
           title=f'{name.title()}')
    plt.show()