<commit_before>import numpy as np
from scipy.special import expit
from scipy.stats import zscore
# Define a class for a glm solver
class glm:

    # Instance object
    distr='poisson'

    # Initialize with distr as poisson by default
    def __init__(self, distr='poisson'):
        self.distr = distr

    #--------------------------
    # Define the nonlinearity
    #--------------------------
    def qu(self, z):
        if(self.distr=='poisson'):
            eps = 0.1
            q = np.log(1+eps+np.exp(z))
        elif(self.distr=='normal'):
            q = z
        elif(self.distr=='binomial'):
            q = expit(z)
        return q

    #-------------------------------------------
    # Define the conditional intensity function
    #-------------------------------------------
    def lmb(self, beta0, beta, x):
        z = beta0 + np.dot(x,beta)
        l = self.qu(z)
        return l

    #---------------------------
    # Define the log-likelihood
    #---------------------------
    def logL(self, beta0, beta, x, y):
        l = self.lmb(beta0, beta, x)
        if(self.distr=='poisson'):
            logL = np.sum(y*np.log(l) - l)
        elif(self.distr=='normal'):
            logL = -0.5*np.sum((y-l)**2)
        elif(self.distr=='binomial'):
            #analytical formula
            #logL = np.sum(y*np.log(l) + (1-y)*np.log(1-l))

            #this prevents underflow
            z = beta0 + np.dot(x,beta)
            logL = np.sum(y*z - np.log(1+np.exp(z)))


        return logL

    #--------------------
    # Define the penalty
    #--------------------
    def penalty(self, alpha, beta):
        P = 0.5*(1-alpha)*np.linalg.norm(beta,2) + alpha*np.linalg.norm(beta,1)
        return P

    #-----------------------------------------------
    # Define the objective function for elastic net
    #-----------------------------------------------
    def loss(self, beta0, beta, alpha, reg_lambda, x, y):
        L = self.logL(beta0, beta, x, y)
        P = self.penalty(alpha, beta)
        J = -L + reg_lambda*P
        return J

    #------------------------------------------------------
    # Define only the differentiable part of the objective
    #------------------------------------------------------
    def L2loss(self, beta0, beta, alpha, reg_lambda, x, y):
        L = self.logL(beta0, beta, x, y)
        P = 0.5*(1-alpha)*np.linalg.norm(beta,2)
        J = -L + reg_lambda*P
        return J

    #------------------------------
    # Define the proximal operator
    #------------------------------
    def prox(self,x,l):
        #sx = [0. if np.abs(y) <= l else np.sign(y)*np.abs(abs(y)-l) for y in x]
        #return np.array(sx).reshape(x.shape)
        return np.sign(x) * (np.abs(x) - l) * (np.abs(x) > l)

    #---------------------
    # Define the gradient
    #---------------------
    def grad_L2loss(self, beta0, beta, alpha, reg_lambda, x, y):
        z = beta0 + np.dot(x, beta)

        if(self.distr=='poisson'):
            q = self.qu(z)
            s = expit(z)
            grad_beta0 = np.sum(s) - np.sum(y*s/q)
            grad_beta = np.transpose(np.dot(np.transpose(s), x) - np.dot(np.transpose(y*s/q), x)) \
                        + reg_lambda*(1-alpha)*beta# + reg_lambda*alpha*np.sign(beta)

        elif(self.distr=='normal'):
            grad_beta0 = -np.sum(y-z)
            grad_beta = -np.transpose(np.dot(np.transpose(y-z), x)) \
                        + reg_lambda*(1-alpha)*beta# + reg_lambda*alpha*np.sign(beta)

        elif(self.distr=='binomial'):
            s = expit(z)
            grad_beta0 =  np.sum(s-y)
            grad_beta = np.transpose(np.dot(np.transpose(s-y), x)) \
                        + reg_lambda*(1-alpha)*beta# + reg_lambda*alpha*np.sign(beta)
        return grad_beta0, grad_beta

    #-------------------------
    # Define the fit function
    #-------------------------
    def fit(self, x, y, reg_params, opt_params, verbose):
    # Implements batch gradient descent (i.e. vanilla gradient descent by computing gradient over entire training set)

        # Dataset shape
        n = x.shape[0]
        p = x.shape[1]

        # Regularization parameters
        reg_lambda = reg_params['reg_lambda']
        alpha = reg_params['alpha']

        # Optimization parameters
        max_iter = opt_params['max_iter']
        e = opt_params['learning_rate']

        # Initialize parameters
        beta0_hat = np.random.normal(0.0,1.0,1)
        beta_hat = np.random.normal(0.0,1.0,[p,1])
        fit = []

        # Outer loop with descending lambda
        if(verbose==True):
            print('---------------------------------------')
            print('Looping throug the regularization path')
            print('---------------------------------------')
        for l,rl in enumerate(reg_lambda):
            fit.append({'beta0': 0., 'beta': np.zeros([p,1])})
            if(verbose==True):
                print('Lambda: %6.4f')% rl

            # Warm initialize parameters
            if(l == 0):
                fit[-1]['beta0'] = beta0_hat
                fit[-1]['beta'] = beta_hat
            else:
                fit[-1]['beta0'] = fit[-2]['beta0']
                fit[-1]['beta'] = fit[-2]['beta']

            #---------------------------
            # Iterate until convergence
            #---------------------------
            no_convergence = 1
            convergence_threshold = 1e-3
            t = 0

            # Initialize parameters
            beta = np.zeros([p+1,1])
            beta[0] = fit[-1]['beta0']
            beta[1:] = fit[-1]['beta']

            # Initialize moment parameters
            g = np.zeros([p+1,1])

            # Initialize cost
            L = []<commit_msg>Multinomial regression added with scikit-learn dependency<commit_after>import numpy as np
from scipy.special import expit
from scipy.stats import zscore
from sklearn.metrics import classification

def softmax(w):
    w = np.array(w)

    maxes = np.amax(w, axis=1)
    maxes = maxes.reshape(maxes.shape[0], 1)
    e = np.exp(w - maxes)
    dist = e / np.sum(e, axis=1, keepdims=True)
    return dist

# Define a class for a glm solver
class glm:

    # Instance object
    distr='poisson'

    # Initialize with distr as poisson by default
    def __init__(self, distr='poisson'):
        self.distr = distr

    #--------------------------
    # Define the nonlinearity
    #--------------------------
    def qu(self, z):
        if(self.distr=='poisson'):
            eps = 0.1
            q = np.log(1+eps+np.exp(z))
        elif(self.distr=='normal'):
            q = z
        elif(self.distr=='binomial'):
            q = expit(z)
        elif(self.distr=='multinomial'):
            q = softmax(z)
        return q

    #-------------------------------------------
    # Define the conditional intensity function
    #-------------------------------------------
    def lmb(self, beta0, beta, x):
        z = beta0 + np.dot(x,beta)
        l = self.qu(z)
        return l

    #---------------------------
    # Define the log-likelihood
    #---------------------------
    def logL(self, beta0, beta, x, y):
        l = self.lmb(beta0, beta, x)
        if(self.distr=='poisson'):
            logL = np.sum(y*np.log(l) - l)
        elif(self.distr=='normal'):
            logL = -0.5*np.sum((y-l)**2)
        elif(self.distr=='binomial'):
            #analytical formula
            #logL = np.sum(y*np.log(l) + (1-y)*np.log(1-l))

            #this prevents underflow
            z = beta0 + np.dot(x,beta)
            logL = np.sum(y*z - np.log(1+np.exp(z)))
        elif(self.distr=='multinomial'):
            logL = classification.log_loss(y, l, normalize=False)
        return logL

    #--------------------
    # Define the penalty
    #--------------------
    def penalty(self, alpha, beta):
        P = 0.5*(1-alpha)*np.linalg.norm(beta,2) + alpha*np.linalg.norm(beta,1)
        return P

    #-----------------------------------------------
    # Define the objective function for elastic net
    #-----------------------------------------------
    def loss(self, beta0, beta, alpha, reg_lambda, x, y):
        L = self.logL(beta0, beta, x, y)
        P = self.penalty(alpha, beta)
        J = -L + reg_lambda*P
        return J

    #------------------------------------------------------
    # Define only the differentiable part of the objective
    #------------------------------------------------------
    def L2loss(self, beta0, beta, alpha, reg_lambda, x, y):
        L = self.logL(beta0, beta, x, y)
        P = 0.5*(1-alpha)*np.linalg.norm(beta,2)
        J = -L + reg_lambda*P
        return J

    #------------------------------
    # Define the proximal operator
    #------------------------------
    def prox(self,x,l):
        #sx = [0. if np.abs(y) <= l else np.sign(y)*np.abs(abs(y)-l) for y in x]
        #return np.array(sx).reshape(x.shape)
        return np.sign(x) * (np.abs(x) - l) * (np.abs(x) > l)

    #---------------------
    # Define the gradient
    #---------------------
    def grad_L2loss(self, beta0, beta, alpha, reg_lambda, x, y):
        z = beta0 + np.dot(x, beta)

        if(self.distr=='poisson'):
            q = self.qu(z)
            s = expit(z)
            grad_beta0 = np.sum(s) - np.sum(y*s/q)
            grad_beta = np.transpose(np.dot(np.transpose(s), x) - np.dot(np.transpose(y*s/q), x)) \
                        + reg_lambda*(1-alpha)*beta# + reg_lambda*alpha*np.sign(beta)

        elif(self.distr=='normal'):
            grad_beta0 = -np.sum(y-z)
            grad_beta = -np.transpose(np.dot(np.transpose(y-z), x)) \
                        + reg_lambda*(1-alpha)*beta# + reg_lambda*alpha*np.sign(beta)

        elif(self.distr=='binomial'):
            s = expit(z)
            grad_beta0 =  np.sum(s-y)
            grad_beta = np.transpose(np.dot(np.transpose(s-y), x)) \
                        + reg_lambda*(1-alpha)*beta# + reg_lambda*alpha*np.sign(beta)
        elif(self.distr=='multinomial'):
            # this assumes that y is already as a one-hot encoding
            pred = self.qu(z)
            grad_beta0 = -np.sum(y - pred)
            grad_beta = -np.transpose(np.dot(np.transpose(y - pred), x)) \
                        + reg_lambda*(1-alpha)*beta


        return grad_beta0, grad_beta

    #-------------------------
    # Define the fit function
    #-------------------------
    def fit(self, x, y, reg_params, opt_params, verbose):
    # Implements batch gradient descent (i.e. vanilla gradient descent by computing gradient over entire training set)

        # Dataset shape
        n = x.shape[0]
        p = x.shape[1]

        # Regularization parameters
        reg_lambda = reg_params['reg_lambda']
        alpha = reg_params['alpha']

        # Optimization parameters
        max_iter = opt_params['max_iter']
        e = opt_params['learning_rate']

        # Initialize parameters
        if self.distr=='multinomial':
            beta0_hat = np.random.normal(0.0,1.0, y.shape[1])
            beta_hat = np.random.normal(0.0, 1.0, [p, y.shape[1]])
        else:
            beta0_hat = np.random.normal(0.0,1.0,1)
            beta_hat = np.random.normal(0.0,1.0,[p,1])
        fit = []

        # Outer loop with descending lambda
        if(verbose==True):
            print('---------------------------------------')
            print('Looping throug the regularization path')
            print('---------------------------------------')
        for l,rl in enumerate(reg_lambda):
            fit.append({'beta0': beta0_hat, 'beta':beta_hat})
            if(verbose==True):
                print('Lambda: %6.4f')% rl

            # Warm initialize parameters
            if(l == 0):
                fit[-1]['beta0'] = beta0_hat
                fit[-1]['beta'] = beta_hat
            else:
                fit[-1]['beta0'] = fit[-2]['beta0']
                fit[-1]['beta'] = fit[-2]['beta']

            #---------------------------
            # Iterate until convergence
            #---------------------------
            no_convergence = 1
            convergence_threshold = 1e-3
            t = 0

            # Initialize parameters
            if self.distr=='multinomial':
                beta = np.zeros([p + 1, y.shape[1]])
            else:
                beta = np.zeros([p+1,1])
            beta[0] = fit[-1]['beta0']
            beta[1:] = fit[-1]['beta']

            # Initialize moment parameters
            if self.distr=='multinomial':
                g = np.zeros([p + 1, y.shape[1]])
            else:
                g = np.zeros([p+1,1])

            # Initialize cost
            L = []
{'fraction_non_alphanumeric': 0.24010096229689226, 'fraction_numerical': 0.014355576589367407, 'mean_word_length': 1.2861521817526145, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4278790', 'n_tokens_mistral': 3982, 'n_tokens_neox': 3809, 'n_words': 818}