Подбор важности фич для k-nearest neighbors (ну или других гиперпараметров) спуском похожим на градиентный

import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
%matplotlib inline
import warnings
warnings.filterwarnings('ignore')

df = pd.read_csv('data/telecom_churn.csv')
df['Voice mail plan'] = pd.factorize(df['Voice mail plan'])[0]
df['International plan'] = pd.factorize(df['International plan'])[0]
df['Churn'] = df['Churn'].astype('int32')
states = df['State']
y = df['Churn']
df.drop(['State','Churn'], axis = 1, inplace=True)
df.head()

%%time
from sklearn.model_selection import train_test_split, StratifiedKFold
from sklearn.tree import DecisionTreeClassifier
from sklearn.neighbors import KNeighborsClassifier
from sklearn.model_selection import GridSearchCV, cross_val_score
from sklearn.metrics import accuracy_score

X_train, X_holdout, y_train, y_holdout = train_test_split(df.values, y, test_size=0.3,
random_state=17)

tree = DecisionTreeClassifier(random_state=17, max_depth=5)
knn = KNeighborsClassifier(n_neighbors=10)

tree_params = {'max_depth': range(1,11), 'max_features': range(4,19)}
tree_grid = GridSearchCV(tree, tree_params, cv=10, n_jobs=-1, verbose=False)
tree_grid.fit(X_train, y_train)
tree_grid.best_params_, tree_grid.best_score_, accuracy_score(y_holdout, tree_grid.predict(X_holdout))

%%time
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler

knn_pipe = Pipeline([('scaler', StandardScaler()), ('knn', KNeighborsClassifier(n_jobs=-1))])
knn_params = {'knn__n_neighbors': range(1, 10)}
knn_grid = GridSearchCV(knn_pipe, knn_params, cv=10, n_jobs=-1, verbose=False)

knn_grid.fit(X_train, y_train)
knn_grid.best_params_, knn_grid.best_score_, accuracy_score(y_holdout, knn_grid.predict(X_holdout))

%%time
feature_importances = pd.DataFrame({'features': df.columns, 'importance':tree_grid.best_estimator_.feature_importances_})
print(feature_importances.sort_values(by=['importance'], inplace=False, ascending=False))

scaler = StandardScaler().fit(X_train)
X_train_transformed = scaler.transform(X_train)
X_train_transformed = X_train_transformed * np.array(feature_importances['importance'])

X_holdout_transformed = scaler.transform(X_holdout)
X_holdout_transformed = X_holdout_transformed * np.array(feature_importances['importance'])

knn_grid = GridSearchCV(KNeighborsClassifier(n_jobs=-1), {'n_neighbors': range(1, 11, 2)}, cv=5, n_jobs=-1, verbose=False)
knn_grid.fit(X_train_transformed, y_train)
print (knn_grid.best_params_, knn_grid.best_score_, accuracy_score(y_holdout, knn_grid.predict(X_holdout_transformed)))

from sklearn.base import TransformerMixin
class StandardAndPoorScaler(StandardScaler, TransformerMixin):
    #normalization = None
    def __init__(self, copy=True, with_mean=True, with_std=True, normalization = None):
        #print("new StandardAndPoorScaler(normalization=", normalization.shape if normalization is not None else normalization, ") // ", type(self))
        self.normalization = normalization
        super().__init__(copy, with_mean, with_std)
    def fit(self, X, y=None):
        #print(type(self),".fit(",X.shape, ",", y.shape if y is not None else "<null>",")")
        super().fit(X, y)
        return self
    def partial_fit(self, X, y=None):
        #print(type(self),".partial_fit(",X.shape, ",", y.shape if y is not None else "<null>)")
        super().partial_fit(X, y)
        if self.normalization is None:
            self.normalization = np.ones((X.shape[1]))
        elif type(self.normalization) != np.ndarray:
            self.normalization = np.array(self.normalization)
        if X.shape[1] != self.normalization.shape[0]:
            raise "X.shape[1]="+X.shape[1]+" in equal self.scale.shape[0]="+self.normalization.shape[0]
    def transform(self, X, copy=None):
        #print(type(self),".transform(",X.shape,",",copy,").self.normalization", self.normalization)
        Xresult = super().transform(X, copy)
        Xresult *= self.normalization
        return Xresult
    def _reset(self):
        #print(type(self),"._reset()")
        super()._reset()
    
scaler = StandardAndPoorScaler(normalization = feature_importances['importance'])
scaler.fit(X = X_train, y = None)
print(scaler.normalization)

%%time
knn_pipe = Pipeline([('scaler', StandardAndPoorScaler()), ('knn', KNeighborsClassifier(n_jobs=-1))])

knn_params = {'knn__n_neighbors': range(1, 11, 4), 'scaler__normalization': [feature_importances['importance']]}
knn_grid = GridSearchCV(knn_pipe, knn_params, cv=5, n_jobs=-1, verbose=False)

knn_grid.fit(X_train, y_train)
knn_grid.best_params_, knn_grid.best_score_, accuracy_score(y_holdout, knn_grid.predict(X_holdout))

feature_base = feature_importances['importance']
searchArea = np.array([feature_base - .05, feature_base, feature_base + .05])
searchArea[searchArea < 0] = 0
searchArea[searchArea > 1] = 1
print(searchArea[2,:] - searchArea[0,:])

import itertools

affected_props = [2,3,4]
parametrs_ranges = np.concatenate([
    np.linspace(searchArea[0,affected_props], searchArea[1,affected_props], 2, endpoint=False),
    np.linspace(searchArea[1,affected_props], searchArea[2,affected_props], 3, endpoint=True)]).transpose()

print(parametrs_ranges) # По пять вариантов значения каждого параметра. Всего 125 комбинаций
recombinations = itertools.product(parametrs_ranges[0],parametrs_ranges[1],parametrs_ranges[1])

variances = []
for item in recombinations: # Не знаю как это вот всё сделать в одну строчку, поэтому будут мои омерзительно медленные в Python циклы.
    varince = feature_base.copy()
    varince[affected_props] = item
    variances.append(varince)
print(variances[0])
print(len(variances))
# Методу knn важно только соотношение, поэтому можно не нормировать вообще или нормировать только при переходе от одной выборки к другой.

%%time
#scale = np.ones([18])
knn_pipe = Pipeline([('scaler', StandardAndPoorScaler()), ('knn', KNeighborsClassifier(n_neighbors = 7 , n_jobs=-1))])

knn_params = {'scaler__normalization': variances} # 'knn__n_neighbors': range(3, 9, 2), 
knn_grid = GridSearchCV(knn_pipe, knn_params, cv=10, n_jobs=-1, verbose=False)

knn_grid.fit(X_train, y_train)
knn_grid.best_params_, knn_grid.best_score_, accuracy_score(y_holdout, knn_grid.predict(X_holdout))

%%time
affected_property = 4
parametrs_range = np.concatenate([
    np.linspace(searchArea[0,affected_property], searchArea[1,affected_property], 29, endpoint=False),
    np.linspace(searchArea[1,affected_property], searchArea[2,affected_property], 30, endpoint=True)]).transpose()

print(searchArea[1,affected_property])
print(parametrs_range) # Cканируем по одному параметру, зато плотно.


variances = []
for item in parametrs_range: # Не знаю как это вот всё сделать в одну строчку, поэтому будут мои омерзительно медленные в Python циклы.
    varince = feature_base.copy()
    varince[affected_property] = item
    variances.append(varince)
print(variances[0])
print(len(variances))
# Методу knn важно только соотношение, поэтому можно не нормировать вообще или нормировать только при переходе от одной выборки к другой.

knn_pipe = Pipeline([('scaler', StandardAndPoorScaler()), ('knn', KNeighborsClassifier(n_neighbors = 7 , n_jobs=-1))])

knn_params = {'scaler__normalization': variances} # 'knn__n_neighbors': range(3, 9, 2), 
knn_grid = GridSearchCV(knn_pipe, knn_params, cv=10, n_jobs=-1, verbose=False)

knn_grid.fit(X_train, y_train)
knn_grid.best_params_, knn_grid.best_score_, accuracy_score(y_holdout, knn_grid.predict(X_holdout))

searchArea = np.array([np.zeros((18,)), np.ones((18,)) /18, np.ones((18,))])
print(searchArea[:,0])

history_parametrs = [searchArea[1,:].copy()]
scaler = StandardAndPoorScaler(normalization=searchArea[1,:])
scaler.fit(X_train)
knn = KNeighborsClassifier(n_neighbors = 7 , n_jobs=-1)
knn.fit(scaler.transform(X_train), y_train)
history_holdout_score = [accuracy_score(y_holdout, knn.predict(scaler.transform(X_holdout)))]

%%time
def changePropertyNormalization(affected_property, points_count = 15):
    test_range = np.concatenate([
        np.linspace(searchArea[0,affected_property], searchArea[1,affected_property], points_count//2, endpoint=False),
        np.linspace(searchArea[1,affected_property], searchArea[2,affected_property], points_count//2 + 1, endpoint=True)]).transpose()
    variances = [searchArea[1,:].copy() for i in range(test_range.shape[0])]
    for row in range(len(variances)):
        variances[row][affected_property] = test_range[row]
    
    knn_pipe = Pipeline([('scaler', StandardAndPoorScaler()), ('knn', KNeighborsClassifier(n_neighbors = 7 , n_jobs=-1))])
    knn_params = {'scaler__normalization': variances} # 'knn__n_neighbors': range(3, 9, 2), 
    knn_grid = GridSearchCV(knn_pipe, knn_params, cv=10, n_jobs=-1, verbose=False)

    knn_grid.fit(X_train, y_train)
    holdout_score = accuracy_score(y_holdout, knn_grid.predict(X_holdout))
    best_param = knn_grid.best_params_['scaler__normalization'][affected_property]
    print(affected_property,
          'property:', searchArea[1, affected_property], "=>", best_param,
          'holdout:', history_holdout_score[-1], "=>", holdout_score, '(', knn_grid.best_score_, ')')
    # А вот теперь надо решить как изменить в связи с этим диапазон поиска.
    before = searchArea[:, affected_property]
    propertySearchArea = searchArea[:, affected_property].copy()
    if best_param == propertySearchArea[0]:
        print('|<<')
        searchArea[0, affected_property] = best_param/2 if best_param > 0.01 else 0
        searchArea[2, affected_property] = (best_param + searchArea[2, affected_property])/2
        searchArea[1, affected_property] = best_param
    elif best_param == propertySearchArea[2]:
        print('>>|')
        searchArea[2, affected_property] = (best_param + 1)/2 if best_param < 0.99 else 1
        searchArea[0, affected_property] = (best_param + searchArea[0, affected_property])/2
        searchArea[1, affected_property] = best_param
    elif best_param < (propertySearchArea[0] + propertySearchArea[1])/2:
        print('<<')
        searchArea[0, affected_property] = max(propertySearchArea[0]*1.1 - .1*propertySearchArea[1], 0)
        searchArea[2, affected_property] = (best_param + propertySearchArea[2])/2
        searchArea[1, affected_property] = best_param
    elif best_param > (propertySearchArea[1] + propertySearchArea[2])/2:
        print('>>')
        searchArea[0, affected_property] = (best_param + propertySearchArea[0])/2
        searchArea[2, affected_property] = min(propertySearchArea[2]*1.1 - .1*propertySearchArea[1], 1)
        searchArea[1, affected_property] = best_param
    elif best_param < propertySearchArea[1]:
        print('<')
        searchArea[2, affected_property] = searchArea[1, affected_property]*.25 + .75*searchArea[2, affected_property]
        searchArea[1, affected_property] = best_param
    elif best_param > propertySearchArea[1]:
        print('>')
        searchArea[0, affected_property] = searchArea[1, affected_property]*.25 + .75*searchArea[0, affected_property]
        searchArea[1, affected_property] = best_param
    else:
        print('=')
        searchArea[0, affected_property] = searchArea[1, affected_property]*.25 + .75*searchArea[0, affected_property]
        searchArea[2, affected_property] = searchArea[1, affected_property]*.25 + .75*searchArea[2, affected_property]
    normalization = searchArea[1,:].sum() #Нормализовать, чтобы параметры не убегали в бесконечность.
    searchArea[:,:] /= normalization
    print(before, "=>",searchArea[:, affected_property])
    history_parametrs.append(searchArea[1,:].copy())
    history_holdout_score.append(holdout_score)
    
changePropertyNormalization(1, 9)
changePropertyNormalization(1, 9)

%%time
# Первоначальное наполнение 
searchArea = np.array([np.zeros((18,)), np.ones((18,)) /18, np.ones((18,))])
print(searchArea[:,0])

history_parametrs = [searchArea[1,:].copy()]
scaler = StandardAndPoorScaler(normalization=searchArea[1,:])
scaler.fit(X_train)
knn = KNeighborsClassifier(n_neighbors = 7 , n_jobs=-1)
knn.fit(scaler.transform(X_train), y_train)
history_holdout_score = [accuracy_score(y_holdout, knn.predict(scaler.transform(X_holdout)))]

for tick in range(40):
    for p in range(searchArea.shape[1]):
        changePropertyNormalization(p, 7)
    
print(searchArea[1,:])
print(history_holdout_score)

feature_importances['knn_importance'] = history_parametrs[-1]
diagramma = feature_importances.copy()
indexes = diagramma.index
diagramma.index = diagramma['features']
diagramma.drop('features', 1, inplace = True)
diagramma.plot(kind='bar');
plt.savefig("images/pic1.png", format = 'png')
plt.show()
feature_importances

parametrs_df = pd.DataFrame(history_parametrs)
parametrs_df['scores'] = history_holdout_score
parametrs_df.index.name = 'index'
parametrs_df.to_csv('parametrs_and_scores.csv')

parametrs_df = pd.read_csv('parametrs_and_scores.csv', index_col = 'index')
history_holdout_score = np.array(parametrs_df['scores'])
parametrs_df.drop('scores',axis=1)
history_parametrs = np.array(parametrs_df.drop('scores',axis=1))

last = history_holdout_score[-1]
steps = np.arange(0, history_holdout_score.shape[0])[history_holdout_score != last].max()
print(steps/18)

%%time
# Задаём произвольные единичные проекции:
import matplotlib.pyplot as plt
%matplotlib inline
import random
import math
random.seed(17)
property_projection = np.array([[math.sin(a), math.cos(a)] for a in [random.uniform(-math.pi, math.pi) for i in range(history_parametrs[0].shape[0])]]).transpose()
history = np.array(history_parametrs[::18]) # Полный цикл - 18 параметров.
# Мне мозгов не зватает чтобы одно на другое помножить без циклов. :(
points = np.array([(history[i] * property_projection).sum(axis=1) for i in range(history.shape[0])])
plt.plot(points[:36,0],points[0:36,1]);
plt.savefig("images/pic2.png", format = 'png')
plt.show()

plt.plot(points[4:36,0],points[4:36,1]);
plt.savefig("images/pic3.png", format = 'png')

plt.plot(points[14:36,0],points[14:36,1]);
plt.savefig("images/pic4.png", format = 'png')
plt.show()
plt.plot(points[24:36,0],points[24:36,1]);
plt.plot(points[35:,0],points[35:,1], color = 'red');
plt.savefig("images/pic5.png", format = 'png')
plt.show()

%%time
import sklearn.manifold as manifold
tsne = manifold.TSNE(random_state=19)
tsne_representation = tsne.fit_transform(history)
plt.plot(tsne_representation[:, 0], tsne_representation[:, 1])
plt.savefig("images/pic6.png", format = 'png')
plt.show();

from time import time
from scipy import linalg
# This list will contain the positions of the map points at every iteration.
positions = []
def _gradient_descent(objective, p0, it, n_iter,
                      n_iter_check=1, n_iter_without_progress=300,
                      momentum=0.8, learning_rate=200.0, min_gain=0.01,
                      min_grad_norm=1e-7, verbose=0, args=None, kwargs=None):
    # The documentation of this function can be found in scikit-learn's code.
    if args is None:
        args = []
    if kwargs is None:
        kwargs = {}

    p = p0.copy().ravel()
    update = np.zeros_like(p)
    gains = np.ones_like(p)
    error = np.finfo(np.float).max
    best_error = np.finfo(np.float).max
    best_iter = i = it

    tic = time()
    for i in range(it, n_iter):
        positions.append(p.copy()) # We save the current position.
        
        check_convergence = (i + 1) % n_iter_check == 0
        # only compute the error when needed
        kwargs['compute_error'] = check_convergence or i == n_iter - 1

        error, grad = objective(p, *args, **kwargs)
        grad_norm = linalg.norm(grad)

        inc = update * grad < 0.0
        dec = np.invert(inc)
        gains[inc] += 0.2
        gains[dec] *= 0.8
        np.clip(gains, min_gain, np.inf, out=gains)
        grad *= gains
        update = momentum * update - learning_rate * grad
        p += update

        if check_convergence:
            toc = time()
            duration = toc - tic
            tic = toc

            if verbose >= 2:
                print("[t-SNE] Iteration %d: error = %.7f,"
                      " gradient norm = %.7f"
                      " (%s iterations in %0.3fs)"
                      % (i + 1, error, grad_norm, n_iter_check, duration))

            if error < best_error:
                best_error = error
                best_iter = i
            elif i - best_iter > n_iter_without_progress:
                if verbose >= 2:
                    print("[t-SNE] Iteration %d: did not make any progress "
                          "during the last %d episodes. Finished."
                          % (i + 1, n_iter_without_progress))
                break
            if grad_norm <= min_grad_norm:
                if verbose >= 2:
                    print("[t-SNE] Iteration %d: gradient norm %f. Finished."
                          % (i + 1, grad_norm))
                break

    return p, error, i

manifold.t_sne._gradient_descent = _gradient_descent

tsne_representation = manifold.TSNE(random_state=17).fit_transform(history)
X_iter = np.dstack(position.reshape(-1, 2) for position in positions)
position_reshape = [position.reshape(-1, 2) for position in positions]
print(position_reshape[0].shape)
print('[0] min', position_reshape[0][:,0].min(),'max', position_reshape[0][:,0].max())
print('[1] min', position_reshape[1][:,0].min(),'max', position_reshape[1][:,0].max())
print('[2] min', position_reshape[2][:,0].min(),'max', position_reshape[2][:,0].max())

%%time
from sklearn.preprocessing import MinMaxScaler
minMaxScaler = MinMaxScaler()
minMaxScaler.fit_transform(position_reshape[0])
position_reshape = [minMaxScaler.fit_transform(frame) for frame in position_reshape]
position_reshape[0].min(), position_reshape[0].max()

%%time

from matplotlib.animation import FuncAnimation, PillowWriter
#plt.style.use('seaborn-pastel')

fig = plt.figure()

ax = plt.axes(xlim=(0, 1), ylim=(0, 1))
line, = ax.plot([], [], lw=3)

def init():
    line.set_data([], [])
    return line,
def animate(i):
    x = position_reshape[i][:,0]
    y = position_reshape[i][:,1]
    line.set_data(x, y)
    return line,

anim = FuncAnimation(fig, animate, init_func=init, frames=36, interval=20, blit=True, repeat_delay = 1000)
anim.save('images/animate_tsne_learning.gif', writer=PillowWriter(fps=5))