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

# In[1]:


import numpy as np
from scipy.linalg import svd
from sklearn.datasets import load_boston
from sklearn.linear_model import RidgeCV as skRidgeCV


# ### Implementation 1
# - based on svd
# - only support fit_intercept=False
# - similar to scikit-learn gcv_mode='svd'

# In[2]:


class RidgeCV():
    def __init__(self, alphas=[0.1, 1.0, 10]):
        self.alphas = alphas

    def fit(self, X, y):
        U, S, VT = svd(X, full_matrices=False)
        cv_values = np.zeros((X.shape[0], len(self.alphas)))
        best_score, best_coef, best_alpha = None, None, None
        for i, alpha in enumerate(self.alphas):
            w = (S ** 2 + alpha) ** -1 - alpha ** -1
            G_inv = np.dot(U * w[np.newaxis, :], U.T) + np.eye(X.shape[0]) * alpha ** -1
            c = np.dot(G_inv, y)
            looe = c / np.diag(G_inv)
            errors = looe ** 2
            cv_values[:, i] = errors
            score = -np.mean(errors)
            if best_score is None or best_score < score:
                best_score = score
                best_coef = c
                best_alpha = alpha
        self.cv_values_ = cv_values
        self.best_score = best_score
        self.coef_ = np.dot(X.T, best_coef)
        self.alpha_ = best_alpha
        return self

    def predict(self, X):
        return np.dot(X, self.coef_)


# In[3]:


X, y = load_boston(return_X_y=True)
X -= X.mean(axis=0)
y -= y.mean()
alpha = [0.1, 1, 10]
reg1 = RidgeCV().fit(X, y)
reg2 = skRidgeCV(fit_intercept=False, store_cv_values=True).fit(X, y)
assert np.allclose(reg1.cv_values_, reg2.cv_values_)
assert np.allclose(reg1.coef_, reg2.coef_)
assert np.allclose(reg1.alpha_, reg2.alpha_)
pred1 = reg1.predict(X)
pred2 = reg2.predict(X)
assert np.allclose(pred1, pred2)