#!/usr/bin/env python
# coding: utf-8

# # 分類のためのモデル
# 
# 前の章で、データとモデルとアルゴリズムの仕組みを紹介したときの具体例として、線形回帰を取り上げた。ここでは、分類課題を見据えたモデルを実装し、前章の内容を補完する。用語については、「__分類__」と「__識別__」を同義的に扱う。

# __目次__
# 
# - <a href="#classifier">識別機のベースクラス</a>
# - <a href="#multiclass">多クラスのロジスティック回帰</a>
# 
# ___

# <a id="classifier"></a>
# ## 識別機のベースクラス
# 
# 実数値を予測する回帰課題と並んで、離散的なラベルを予測する分類課題が大変重要である。その予測をするシステムのことを主に「識別機」と呼ぶ。前に用意したモデルのベースクラス`Model`を踏襲して、識別機の新しい範疇を作るため、以下のように`Classifier`を実装する。

# In[1]:


import models

class Classifier(models.Model):
    '''
    Generic classifier model, an object with methods
    for both training and evaluating classifiers.
    '''

    def __init__(self, data=None, name=None):
        super(Classifier, self).__init__(name=name)

        # If given data, collect information about labels.
        if data is not None:
            self.labels = self.get_labels(data=data) # all unique labels.
            self.nc = self.labels.size # number of unique labels.


    def onehot(self, y):
        '''
        A function for encoding y into a one-hot vector.
        Inputs:
        - y is a (k,1) array, taking values in {0,1,...,nc-1}.
        '''
        nc = self.nc
        k = y.shape[0]
        C = np.zeros((k,nc), dtype=y.dtype)

        for i in range(k):
            j = y[i,0] # assumes y has only one column.
            C[i,j] = 1

        return C
        

    def get_labels(self, data):
        '''
        Get all the (unique) labels that appear in the data.
        '''
        A = (data.y_tr is None)
        B = (data.y_te is None)

        if (A and B):
            raise ValueError("No label data provided!")
        else:
            if A:
                out_labels = np.unique(data.y_te)
            elif B:
                out_labels = np.unique(data.y_tr)
            else:
                out_labels = np.unique(np.concatenate((data.y_tr,
                                                       data.y_te), axis=0))
            count = out_labels.size
            return out_labels.reshape((count,1))


    def classify(self, X):
        '''
        Must be implemented by sub-classes.
        '''
        raise NotImplementedError


    def class_perf(self, y_est, y_true):
        '''
        Given class label estimates and true values,
        compute the fraction of correct classifications
        made for each label, yielding typical binary
        classification performance metrics.

        Input:
        y_est and y_true are (k x 1) matrices of labels.

        Output:
        Returns a dictionary with two components, (1) being
        the fraction of correctly classified labels, and
        (2) being a dict of per-label precison/recall/F1
        scores. 
        '''
        
        # First, get the classification rate.
        k = y_est.size
        num_correct = (y_est == y_true).sum()
        frac_correct = num_correct / k
        frac_incorrect = 1.0 - frac_correct

        # Then, get precision/recall for each class.
        prec_rec = { i:None for i in range(self.nc) } # initialize

        for c in range(self.nc):

            idx_c = (y_true == c)
            idx_notc = (idx_c == False)

            TP = (y_est[idx_c] == c).sum()
            FN = idx_c.sum() - TP
            FP = (y_est[idx_notc] == c).sum()
            TN = idx_notc.sum() - FP

            # Precision.
            if (TP == 0 and FP == 0):
                prec = 0
            else:
                prec = TP / (TP+FP)

            # Recall.
            if (TP == 0 and FN == 0):
                rec = 0
            else:
                rec = TP / (TP+FN)

            # F1 (harmonic mean of precision and recall).
            if (prec == 0 or rec == 0):
                f1 = 0
            else:
                f1 = 2 * prec * rec / (prec + rec)

            prec_rec[c] = {"P": prec,
                           "R": rec,
                           "F1": f1}

        return {"rate": frac_incorrect,
                "PRF1": prec_rec}

    # Need to re-implement l_tr, l_te, g_tr, and
    # g_te in order to utilize C_*, the one-hot
    # representation. Doing it this way lets us
    # get around having to compute it every time
    # we evaluate a loss/grad.
    def l_tr(self, w, data, n_idx=None, lamreg=None):
        if n_idx is None:
            return self.l_imp(w=w, X=data.X_tr,
                              C=self.C_tr,
                              lamreg=lamreg)
        else:
            return self.l_imp(w=w, X=data.X_tr[n_idx,:],
                              C=self.C_tr[n_idx,:],
                              lamreg=lamreg)
    
    def l_te(self, w, data, n_idx=None, lamreg=None):
        if n_idx is None:
            return self.l_imp(w=w, X=data.X_te,
                              C=self.C_te,
                              lamreg=lamreg)
        else:
            return self.l_imp(w=w, X=data.X_te[n_idx,:],
                              C=self.C_te[n_idx,:],
                              lamreg=lamreg)

    def g_tr(self, w, data, n_idx=None, lamreg=None):
        if n_idx is None:
            return self.g_imp(w=w, X=data.X_tr,
                              C=self.C_tr,
                              lamreg=lamreg)
        else:
            return self.g_imp(w=w, X=data.X_tr[n_idx,:],
                              C=self.C_tr[n_idx,:],
                              lamreg=lamreg)
    
    def g_te(self, w, data, n_idx=None, lamreg=None):
        if n_idx is None:
            return self.g_imp(w=w, X=data.X_te,
                              C=self.C_te,
                              lamreg=lamreg)
        else:
            return self.g_imp(w=w, X=data.X_te[n_idx,:],
                              C=self.C_te[n_idx,:],
                              lamreg=lamreg)


# __練習問題__
# 
# 0. 次のメソッドの役割を説明すること。`onehot()`, `get_labels()`, `class_perf()`.
# 0. `class_perf`が何を返すか説明すること。

# <a id="multiclass"></a>
# ## 多クラスのロジスティック回帰
# 
# 簡明で使い勝手の良い多クラス識別モデルとして、多クラスロジスティック回帰は代表的である。要素を整理しておこう。
# 
# - データ： 事例とラベルからなるペア$(x,y)$. このデータは$x \in \mathbb{R}^{d}$および$y \in \{0,1,\ldots,K\}$となっており、クラスの数は$K+1$である。
# 
# - モデル： 事例$x$を与えられたときの各クラスの条件付き確率を、パラメータ$w_{0},w_{1},\ldots,w_{K} \in \mathbb{R}^{d}$を用いながら下記の通りにモデル化する。
# 
# \begin{align*}
# P\{y = j | x\} = \frac{\exp(w_{j}^{T}x)}{\sum_{k=0}^{K}\exp(w_{k}^{T}x)}.
# \end{align*}
# 
# - 損失関数： 負の対数尤度$(-1) \log P\{y | x\}$をサンプル全体に対して最小にする（詳細は後ほど）。

# 我々の目的はあくまで実装することであるから、計算の手助けとなるような表記を導入する。特に、損失関数を求めやすい形で表わすと便利である。まず、独立に得られた$n$個の標本$(x_{1},y_{1}), \ldots, (x_{n},y_{n})$を基に、各$i=1,\ldots,n$と$j=0,\ldots,K$に対して、以下のように定義する。
# 
# \begin{align*}
# p_{ij} & = P\{y_{i} = j | x_{i}\}\\
# c_{ij} & = I\{y_{i} = j\}.
# \end{align*}
# 
# インデックス$j$を追って、すべての要素をベクトルに並べ、以下のように記す。
# 
# \begin{align*}
# p_{i} & = (p_{i0},\ldots,p_{iK})\\
# c_{i} & = (c_{i0},\ldots,c_{iK}).
# \end{align*}
# 
# この$c_{i}$は元のラベル$y_{i}$のone-hotベクトル表現で、使い勝手が良いものである。
# 
# この表記にしたがって、$i$番目の事例のクラス$j$である確率は$p_{ij}^{c_{ij}}$と書ける。事例ごとの損失関数を定義するならば、以下のごとく負の対数尤度を使うことが自然である。
# 
# \begin{align*}
# L(w;x_{i},y_{i}) & = (-1)\log \prod_{j=0}^{K} p_{ij}^{c_{ij}}\\
# & = (-1)\sum_{j=0}^{K} c_{ij} \log p_{ij}\\
# & = (-1)\sum_{j=0}^{K} c_{ij}\left(w_{j}^{T} x_{i} - \log\left( \sum_{k=0}^{K}\exp(w_{k}^{T}x_{i}) \right) \right)\\
# & = \log\left( \sum_{k=0}^{K}\exp(w_{k}^{T}x_{i}) \right) - \sum_{j=0}^{K} c_{ij}w_{j}^{T} x_{i}.
# \end{align*}
# 
# これまでと同様に、ERM学習則（経験期待損失最小化）に従い、このロスのサンプル平均を目的関数とする（但し、$n$をかけている）。
# 
# \begin{align*}
# \min_{w} \sum_{i=1}^{n} L(w;x_{i},y_{i}) \to \hat{w}.
# \end{align*}
# 
# $w$の関数として、この目的関数は微分可能であり、凸関数でもある。実際に最小化に手をつけるとき、勾配降下法はしばしば使われる。この$L(w;x,y)$の勾配ベクトルを$(w,x_{i},y_{i})$という点において、$w$について求めると、いささか詳細を省きながら下記のように書ける。
# 
# \begin{align*}
# \nabla L(w;x_{i},y_{i}) = (p_{i}-c_{i}) \otimes x_{i}.
# \end{align*}
# 
# この演算子$\otimes$はクロネッカー積と呼ばれ、2つの行列に対して定義されている。たとえば、$U$が$a \times b$で、$V$が$c \times d$であるとすると、これらのクロネッカー積は
# 
# \begin{align*}
# U \otimes V =
# \begin{bmatrix}
# u_{1,1}V & \cdots & u_{1,b}V\\
# \vdots & \ddots & \vdots\\
# u_{a,1}V & \cdots & u_{a,b}V
# \end{bmatrix}
# \end{align*}
# 
# という形を取る。言うまでもなく、各$u_{i,j1}V$は$c \times d$の区分（ブロック）行列である。よって、全体として、この積$U \otimes V$の寸法は$ac \times bd$となる。今回の$\nabla L(w;x_{i},y_{i})$については、$(p_{i}-c_{i})$の寸法が$1 \times (K+1)$で、$x_{i}$の寸法が$1 \times d$であるため、勾配ベクトルは$1 \times d(K+1)$となる。
# 
# クロネッカー積は簡単なループで自分で実装できるのだが、効率的なものはNumpyに標準搭載されているので、それを使うと何の苦労もなくロジスティック回帰の実装ができる。
# 
# __パラメータの数について：__ 制御の対象となるパラメータを$w = (w_{0},\ldots,w_{K}$として、$d(k+1)$を決めていくことはまったく問題ない。ただし、確率をモデル化していることから、$K$個のクラスの条件付き確率が定まれば、必然的に残りの1クラスの確率も決まる。たとえば、すでに$w_{1},\ldots,w_{K}$の候補が定まっているのであれば、当然$p_{i1},\ldots,p_{iK}$も定まるので、残っている0番目のクラスの条件付き確率は自ずと$p_{i0} = 1 - \sum_{j=1}^{K}p_{ij}$で決まってくる。そう考えると、最初から$w_{0}=(0,\ldots,0)$と固定して、その他の$dK$個のパラメータだけ決めれば良い、ということになる。`Classifier`を踏襲する`LogisticReg`を以下の通りに実装しているが、パラメータの数は$dK$である。

# In[ ]:


class LogisticReg(Classifier):
    '''
    Multi-class logistic regression model.
    '''

    def __init__(self, data=None):
        
        # Given data info, load up the (X,y) data.
        super(LogisticReg, self).__init__(data=data)
        
        # Convert original labels to a one-hot binary representation.
        if data.y_tr is not None:
            self.C_tr = self.onehot(y=data.y_tr)
        if data.y_te is not None:
            self.C_te = self.onehot(y=data.y_te)
        
        
    def classify(self, w, X):
        '''
        Given learned weights (w) and a matrix of one or
        more observations, classify them as {0,...,nc-1}.

        Input:
        w is a (d x 1) matrix of weights.
        X is a (k x numfeat) matrix of k observations.
        NOTE: k can be anything, the training/test sample size.

        Output:
        A vector of length k, housing labels in {0,...,nc-1}.
        '''
        
        k, numfeat = X.shape
        A = np.zeros((self.nc,k), dtype=np.float32)
        
        # Get activations, with last row as zeros.
        A[:-1,:] = w.reshape((self.nc-1, numfeat)).dot(X.T)
        
        # Now convert activations to conditional probabilities.
        maxes = np.max(A, axis=0) # largest score for each obs.
        A = A - maxes
        A = np.exp(A)
        A = A / A.sum(axis=0)  # (nc x k).
        
        # Assign classes with highest probability, (k x 1) array.
        return A.argmax(axis=0).reshape((k,1))


    def l_imp(self, w, X, C, lamreg=None):
        '''
        Implementation of the multi-class logistic regression
        loss function.

        Input:
        w is a (d x 1) matrix of weights.
        X is a (k x numfeat) matrix of k observations.
        C is a (k x nc) matrix giving a binarized encoding of the
        class labels for each observation; each row a one-hot vector.
        lam is a non-negative regularization parameter.
        NOTE: k can be anything, the training/test sample size.

        Output:
        A vector of length k with losses evaluated at k points.
        '''
        
        k, numfeat = X.shape
        A = np.zeros((self.nc,k), dtype=np.float64)
        
        # Get activations, with last row as zeros.
        A[:-1,:] = w.reshape((self.nc-1, numfeat)).dot(X.T)
        
        # Raw activations of all the correct weights.
        cvec = (A*C.T).sum(axis=0)
        
        # Compute the negative log-likelihoods.
        maxes = np.max(A, axis=0)
        err = (np.log(np.exp(A-maxes).sum(axis=0))+maxes)-cvec

        # Return the losses (all data points), with penalty if needed.
        if lamreg is None:
            return err
        else:
            penalty = lamreg * np.linalg.norm(W)**2
            return err + penalty
    
    
    def l_tr(self, w, data, n_idx=None, lamreg=None):
        if n_idx is None:
            return self.l_imp(w=w, X=data.X_tr,
                              C=self.C_tr,
                              lamreg=lamreg)
        else:
            return self.l_imp(w=w, X=data.X_tr[n_idx,:],
                              C=self.C_tr[n_idx,:],
                              lamreg=lamreg)
    
    def l_te(self, w, data, n_idx=None, lamreg=None):
        if n_idx is None:
            return self.l_imp(w=w, X=data.X_te,
                              C=self.C_te,
                              lamreg=lamreg)
        else:
            return self.l_imp(w=w, X=data.X_te[n_idx,:],
                              C=self.C_te[n_idx,:],
                              lamreg=lamreg)
        
        
    def g_imp(self, w, X, C, lamreg=0):
        '''
        Implementation of the gradient of the loss function used in
        multi-class logistic regression.

        Input:
        w is a (d x 1) matrix of weights.
        X is a (k x numfeat) matrix of k observations.
        C is a (k x nc) matrix giving a binarized encoding of the
        class labels for each observation; each row a one-hot vector.
        lamreg is a non-negative regularization parameter.
        NOTE: k can be anything, the training/test sample size.

        Output:
        A (k x d) matrix of gradients eval'd at k points.
        '''
        
        k, numfeat = X.shape
        A = np.zeros((self.nc,k), dtype=np.float32)
        
        # Get activations, with last row as zeros.
        A[:-1,:] = w.reshape((self.nc-1, numfeat)).dot(X.T)
        
        # Now convert activations to conditional probabilities.
        maxes = np.max(A, axis=0) # largest score for each obs.
        A = A - maxes
        A = np.exp(A)
        A = A / A.sum(axis=0)  # (nc x k).
        
        # Initialize a large matrix (k x d) to house per-point grads.
        G = np.zeros((k,w.size), dtype=w.dtype)

        for i in range(k):
            # A very tall vector (i.e., just one "axis").
            G[i,:] = np.kron(a=(A[:-1,i]-C[i,:-1]), b=X[i,:])
            # Note we carefully remove the last elements.
        
        if lamreg is None:
            return G
        else:
            return G + lamreg*2*w.T
        

    def g_tr(self, w, data, n_idx=None, lamreg=None):
        if n_idx is None:
            return self.g_imp(w=w, X=data.X_tr,
                              C=self.C_tr,
                              lamreg=lamreg)
        else:
            return self.g_imp(w=w, X=data.X_tr[n_idx,:],
                              C=self.C_tr[n_idx,:],
                              lamreg=lamreg)
    
    def g_te(self, w, data, n_idx=None, lamreg=None):
        if n_idx is None:
            return self.g_imp(w=w, X=data.X_te,
                              C=self.C_te,
                              lamreg=lamreg)
        else:
            return self.g_imp(w=w, X=data.X_te[n_idx,:],
                              C=self.C_te[n_idx,:],
                              lamreg=lamreg)



# __練習問題__
# 
# 0. 上記`LogisticReg`の`classify`, `l_imp`および`g_imp`というメソッドの中身を説明すること。
# 
# 0. もしクラスが2つしかなく、$y \in \{0,1\}$という状況であれば、従来のロジスティック回帰では、$d$次元ベクトル$w$を一つだけ決めるようになっている。その$w$は$x \mapsto f(w^{T}x) = P\{y = 1 | x\}$の写像に使う。この$f(u) = 1/(1+\exp(-u))$はロジスティック関数と呼ばれる。このときの勾配はかなり単純な形になる。そこで、2クラスロジスティック回帰のモデルを、上記の`np.kron`を__使わずに__、実装してみること。自作の`g_imp`と`LogisticReg`の`g_imp`の出力が数値的に一致することを自分で確かめること。
# 
# 0. なぜ`l_tr`, `l_te`, `g_tr`, `g_te`をわざわざ実装しなおしたか。
# 
# 0. `maxes`を使って、何をしているか。何のためにこの演算を行なっているか。ちなみに、任意のベクトル$\mathbf{u}$とスカラー$a$に対して、以下の等式が成り立つことを利用している。
# 
# \begin{align*}
# \text{softmax}(\mathbf{u}+a) & = \text{softmax}(\mathbf{u})\\
# \log \left( \sum_{j} \exp(u_{j}) \right) & = a + \log \left( \sum_{j} \exp(u_{j}-a) \right)
# \end{align*}
# 

# ___