# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# ## Imports

# +
import itertools
import numpy as np
import pandas as pd
import seaborn as sns
import scipy.stats as st
import matplotlib.pyplot as plt

from sklearn import linear_model
from sklearn.model_selection import KFold, train_test_split
from sklearn.ensemble import RandomForestRegressor
# -

# ## Данные

df = pd.read_csv('Очищенные обработанные данные.csv', sep=',', index_col=0)#, parse_dates=['DT'])
df.head()

TARGETS = ['химшлак последний Al2O3', 'химшлак последний CaO',
           'химшлак последний R', 'химшлак последний SiO2']

# убираем колонки с "будущим"
last_analysis_cols = [col for col in df.columns if (col.find('последний') != -1) & (col not in TARGETS)]

last_analysis_cols

df.drop(last_analysis_cols, axis=1, inplace=True)

len(df.columns)

# ## Feature Engineering

pearson_corr = df.drop(TARGETS, axis=1).corr(method='pearson', min_periods=50)
pos_corr = np.abs(np.nan_to_num(pearson_corr.values, 0.0))
strong_corr_pairs = []
for i in range(1, len(pos_corr)):
    if np.max(pos_corr[i, :i]) >= 0.7:
        j = np.argmax(pos_corr[i, :i])
        strong_corr_pairs.append((i, j, pearson_corr.values[i, j]))

cols = df.drop(TARGETS, axis=1).columns
cols_dict = dict(list(zip(range(len(cols)), cols)))

# сильно коррелирующие фичи
strong_corr_df = pd.DataFrame(strong_corr_pairs).replace(cols_dict).sort_values(2, ascending=False)
strong_corr_df

pairs = strong_corr_df[[0,1]].values
pairs

# #### Imputter

# +
from sklearn.impute import KNNImputer

imputer = KNNImputer(n_neighbors=2)
X = pd.DataFrame(imputer.fit_transform(df.drop(TARGETS, axis=1)), columns=df.drop(TARGETS, axis=1).columns)
# -

df[X.columns] = X.values


# Попробуем отношение.

# +
def safe_division(x, y):
    if (x != x) | (y != y) | (y == 0):
        return np.nan
    return x / y

for pair in pairs:
    new_col = pair[0]+'_'+pair[1]+'_ratio'
    df[new_col] = df.apply(lambda x: safe_division(x[pair[0]],x[pair[1]]), axis=1)
# -

import itertools
t_features = ['t вып-обр', 't обработка', 't под током', 't продувка']
t_combinations = list(itertools.combinations(t_features, 2))

for pair in t_combinations:
    new_col = pair[0]+'_'+pair[1]+'_sub'
    df[new_col] = df.apply(lambda x: abs(x[pair[0]]-x[pair[1]]), axis=1)

df.drop(['произв количество обработок'], axis=1, inplace=True)

df.shape

# ## Feature Importance

for col in df.columns:
    if (df[col].nunique() <= 50) & (col not in TARGETS):
        df[col] = df[col].astype('category')

NUMERICAL = df.select_dtypes(exclude=['object', 'datetime64']).columns.tolist()
ORDINAL = df.select_dtypes(include=['category']).columns.tolist()
for tar in TARGETS:
    if tar in NUMERICAL:
        NUMERICAL.remove(tar)
    if tar in ORDINAL:
        ORDINAL.remove(tar)

# #### Correlations

# +
correlations = dict()

for tar in TARGETS:
    correlations[tar] = dict()
    for col in NUMERICAL:
        correlations[tar][col] = df[col].corr(df[tar])
# -

corrs = pd.DataFrame(correlations)
corrs[(abs(corrs) > 0.5).any(1)]

golden_features = list(corrs[(abs(corrs) > 0.3).any(1)].index)

# ## Model

df[TARGETS].count() / df.shape[0]


def make_train_test_for_target(df, target, cols_exclude_from_sample=[]):
    new_df = df[~df[target].isna()]
    return new_df.drop(set(cols_exclude_from_sample)-set([target]), axis=1)


# ### Linear Regression

from sklearn.linear_model import LinearRegression, ElasticNet

target = 'химшлак последний Al2O3'

data = make_train_test_for_target(df, target, TARGETS)

# +
import scipy.stats as ss

fig = plt.figure()
ax1 = fig.add_subplot(211)
x = data[target].values
prob = ss.probplot(x, dist=ss.norm, plot=ax1)
ax1.set_xlabel('')
ax1.set_title('Probplot against normal distribution')

ax2 = fig.add_subplot(212)
xt, _ = ss.boxcox(x)
prob = ss.probplot(xt, dist=ss.norm, plot=ax2)
ax2.set_title('Probplot after Box-Cox transformation')
plt.show()
# -

X, y = data.drop(target, axis=1), data[target]
features = X.columns

train_X, val_X, train_y, val_y = train_test_split(X, np.log(y+0.0001), random_state=0, test_size=.4, shuffle=True)

# +
# # !pip install sklearn_pandas --user

# +
from sklearn.preprocessing import StandardScaler, MinMaxScaler
from sklearn_pandas import DataFrameMapper, gen_features

def map_features(features=[]):
    numerical_def = gen_features(
        columns=[[c] for c in NUMERICAL if c in features],
        classes=[
            {'class': StandardScaler}
        ]
    )

    ordinal_def = gen_features(
        columns=[[c] for c in ORDINAL if c in features],
        classes=[
            {'class': MinMaxScaler}
        ]
    )
    return numerical_def #+ ordinal_def


# -

mapper1 = DataFrameMapper(map_features(features), df_out=True)

# +
from sklearn.metrics import r2_score, mean_absolute_error, mean_squared_error
from sklearn.model_selection import cross_validate

def mean_absolute_percentage_error(y_true, y_pred): 
    y_true, y_pred = np.array(y_true), np.array(y_pred)
    return np.mean(np.abs((y_true - y_pred) / y_true)) * 100

def evaluate_model(train, val, tr_y, val_y, features, est):
    mapper = DataFrameMapper(features, df_out=True)
#     трансформим отдельно трейн (фиттим) и тест
    train = mapper.fit_transform(train)
    val = mapper.transform(val)

    est.fit(train, tr_y)
    pred_val = est.predict(val)
    pred_train = est.predict(train)

    return pd.DataFrame({
        'train_R2': [r2_score(tr_y, pred_train)],
        'train_MAE': [mean_absolute_error(tr_y, pred_train)], 
        'train_MSE': [mean_squared_error(tr_y, pred_train)], 
        'train_MAPe': [mean_absolute_percentage_error(tr_y, pred_train)],
        'test_R2': [r2_score(val_y, pred_val)],
        'test_MAE': [mean_absolute_error(val_y, pred_val)],
        'test_MSE': [mean_squared_error(val_y, pred_val)],
        'test_MAPe': [mean_absolute_percentage_error(val_y, pred_val)]
    })


# -

from sklearn.linear_model import LassoCV, RidgeCV
reg_Al2O3 = RidgeCV()
scores = evaluate_model(train_X, val_X, train_y, val_y, map_features(features), reg_Al2O3)
scores

_ = mapper1.fit(train_X)
reg_Al2O3.fit(mapper1.transform(train_X), train_y)

import pickle
filename = 'Al2O3.sav'
# pickle.dump(reg_Al2O3, open(filename, 'wb'))

preds = reg_Al2O3.predict(mapper1.transform(data[features]))
fig = plt.figure(figsize=(15,10))
plt.plot(np.exp(preds), label='preds')
plt.plot(data['химшлак последний Al2O3'].values,alpha=0.5, label='true')
plt.legend(fontsize='14')

# #### Следующий таргет

from imblearn.under_sampling import RandomUnderSampler
target = 'химшлак последний R'
data = make_train_test_for_target(df, target, set(TARGETS)) # добавим колонку Al2O3
data['химшлак последний Al2O3'] = np.exp(reg_Al2O3.predict(mapper1.transform(data[features]))) # добавим предыдущее предсказание
data = data[~(data['химшлак последний R'].isin([3.1, 1.9, 2.12, 2.9, 3.0]))] # эти значения встречаются по 1-2 раза
features_2 = set(data.columns) - set([target])
NUMERICAL += ['химшлак последний Al2O3']
X, y = data[features_2], data[target]

from collections import Counter
Counter(data['химшлак последний R'])

from imblearn.over_sampling import SMOTENC
# дискретизируем значения и сделаем оверсэмплинг по "классам"
ord_cols = [X.columns.get_loc(c) for c in X.columns if c in ORDINAL]
smote_nc = SMOTENC(categorical_features=ord_cols,sampling_strategy='not majority',random_state=0, k_neighbors=10)
x_rus, y_rus = smote_nc.fit_resample(X, y.astype(str))

Counter(y_rus)

train_X, val_X, train_y, val_y = train_test_split(x_rus, y_rus.astype(float), random_state=0, test_size=.4, shuffle=True, stratify=y_rus)

from sklearn.linear_model import LassoCV, RidgeCV
from sklearn.ensemble import RandomForestRegressor
reg_R = RandomForestRegressor(n_estimators=10, max_depth=5, random_state=0)
scores = evaluate_model(train_X, val_X, train_y, val_y, map_features(features_2), reg_R)
scores

mapper2 = DataFrameMapper(map_features(features_2), df_out=True)
_ = mapper2.fit(train_X[features_2])
reg_R.fit(mapper2.transform(train_X), train_y)

preds = reg_R.predict(mapper2.transform(data[features_2]))
fig = plt.figure(figsize=(15,10))
plt.plot(np.around(preds, 2), label='preds')
plt.plot(data['химшлак последний R'].values, alpha=0.5, label='true')
plt.legend(fontsize='14')

filename = 'R.sav'
# pickle.dump(reg_R, open(filename, 'wb'))

# #### Следующий таргет

df['химшлак последний SiO2'].hist() # с бимодальным распределением сложно работать

target = 'химшлак последний SiO2'
data = make_train_test_for_target(df, target, TARGETS)
data['химшлак последний Al2O3'] = np.exp(reg_Al2O3.predict(mapper1.transform(data[features]))) # добавим предыдущее предсказание
data['химшлак последний R'] = np.around(reg_R.predict(mapper2.transform(data[features_2])), 2) # добавим предыдущее предсказание
data['химшлак последний R_химшлак последний Al2O3_mul'] = data['химшлак последний Al2O3']*data['химшлак последний R'] # из логики и корреляции фич
features_4 = set(data.columns) - set([target])
# ORDINAL += ['химшлак последний R']
NUMERICAL += ['химшлак последний R_химшлак последний Al2O3_mul', 'химшлак последний R']
X, y = data[features_4], data[target]

# попытка стратификации по значениям
y_strat = y.copy()
y_strat.loc[y_strat <=24] = 22
y_strat.loc[y_strat > 24] = 26

# train_X, val_X, train_y, val_y = train_test_split(X, y, random_state=0, test_size=.4, shuffle=True, stratify=y_strat.astype(str))
train_X, val_X, train_y, val_y = train_test_split(X, y, random_state=0, test_size=.4, shuffle=True)

from sklearn.linear_model import LogisticRegression
# reg_SiO2 = RandomForestRegressor(n_estimators=50, max_depth=3)
reg_SiO2 = RidgeCV(cv=3)
scores = evaluate_model(train_X, val_X, train_y, val_y, map_features(features_4), reg_SiO2)
scores

reg_SiO2.coef_

mapper4 = DataFrameMapper(map_features(features_4), df_out=True)
_ = mapper4.fit(train_X[features_4])

reg_SiO2.fit(mapper4.transform(train_X), train_y)

preds = reg_SiO2.predict(mapper4.transform(X))
fig = plt.figure(figsize=(15,10))
# plt.plot((preds-np.mean(preds))*np.std(train_y)/np.std(preds)+np.mean(preds))
plt.plot(preds)
plt.plot(data['химшлак последний SiO2'].values, alpha=0.5)

filename = 'SiO2.sav'
# pickle.dump(reg_SiO2, open(filename, 'wb'))

# #### Последний таргет

df['химшлак последний CaO'].hist()

target = 'химшлак последний CaO'
data = make_train_test_for_target(df, target, [target]) # добавим колонку Al2O3
data['химшлак последний Al2O3'] = reg_Al2O3.predict(mapper1.transform(data[features])) # добавим предыдущее предсказание
data['химшлак последний R'] = reg_R.predict(mapper2.transform(data[features_2])) # добавим предыдущее предсказание
data['химшлак последний R_химшлак последний Al2O3_mul'] = data['химшлак последний Al2O3']*data['химшлак последний R']
data['химшлак последний SiO2'] = reg_SiO2.predict(mapper4.transform(data[features_4]))
features_3 = set(data.columns) - set([target])
NUMERICAL += ['химшлак последний SiO2']
X, y = data[features_3], data[target]

# +
fig = plt.figure()
ax1 = fig.add_subplot(211)
x = data[target].values
prob = ss.probplot(x, dist=ss.norm, plot=ax1)
ax1.set_xlabel('')
ax1.set_title('Probplot against normal distribution')

ax2 = fig.add_subplot(212)
xt, _ = ss.boxcox(x)
prob = ss.probplot(xt, dist=ss.norm, plot=ax2)
ax2.set_title('Probplot after Box-Cox transformation')
plt.show()
# -

train_X, val_X, train_y, val_y = train_test_split(X, y, random_state=0, test_size=.2, shuffle=True)

# reg_CaO = RandomForestRegressor(n_estimators=50, max_depth=3)
reg_CaO = RidgeCV(cv=3)
scores = evaluate_model(train_X, val_X, train_y, val_y, map_features(features_3), reg_CaO)
scores

reg_CaO.coef_

filename = 'CaO.sav'
# pickle.dump(reg_CaO, open(filename, 'wb'))

len(df.columns)



StarcoderdataJupyterScriptsDedupFiltered