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

# +
import numpy as np
from sklearn.model_selection import cross_val_score
from sklearn.multioutput import MultiOutputRegressor
import xgboost as xgb
from sklearn.linear_model import MultiTaskElasticNet, Ridge, RidgeCV, ElasticNet, ElasticNetCV, Lasso, MultiTaskLassoCV, LinearRegression
from sklearn.metrics import mean_squared_error, accuracy_score, r2_score, explained_variance_score, mean_absolute_percentage_error
from sklearn.model_selection import RandomizedSearchCV
from sklearn.model_selection import KFold
from scipy.io import loadmat, savemat
from sklearn.preprocessing import StandardScaler
from clip import load
from big_spose_sleep import create_clip_img_transform
import glob, pickle
import torch
from os.path import exists
from scipy.stats import pearsonr
from ridge import ridge, ridge_corr, bootstrap_ridge

from sklearn.ensemble import GradientBoostingClassifier
#https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.ElasticNet.html#sklearn.linear_model.ElasticNet

from torchvision import datasets, transforms

from sklearn.datasets import make_regression
from sklearn.model_selection import train_test_split

from sklearn.multioutput import MultiOutputRegressor
from sklearn.ensemble import GradientBoostingRegressor, HistGradientBoostingRegressor


# -

def calc_corrs(timesrs_test, timesrs_pred): 
    corr_by_ft = {}

    ft_n = timesrs_test.shape[1]
    corr_by_ft = np.array([ pearsonr(timesrs_test[:,ft_i], 
                                     timesrs_pred[:,ft_i])[0] for ft_i in range(ft_n) ]).astype('float32')
    return corr_by_ft


# ## Extract CLIP features

# +
thingsroot = "/Users/katja/Documents/Data/THINGS/"

thingsimgfns = thingsroot+"/images/{}/*.*"

regr_data_fn = "../../data_spose_to_clip.mat"

reextract_data = False

thingscats = []
with open(thingsroot+"THINGS_unique_IDs.txt", 'r') as handle:
    lines = handle.readlines()
    for line in lines:
        line = line.strip()
        if len(line)>0:
            thingscats.append(line)

assert(len(thingscats)==1854)

spose_cat_emb = np.loadtxt(thingsroot+"spose_embedding_49d_sorted.txt")

assert(len(thingscats)==spose_cat_emb.shape[0])

clip_perceptor, _ = load('ViT-B/32', jit = False)
clip_transform = create_clip_img_transform(224)

# +
thingsimgs = datasets.ImageFolder(thingsroot+'/images/', transform=clip_transform)
thingsloader = torch.utils.data.DataLoader(thingsimgs, batch_size=1, shuffle=False)

print("Number of THINGS images found:", len(thingsloader))

print(clip_transform)

num_imgs = len(thingsloader)

# +
# Looping through it, get a batch on each loop, check data shape

show_shapes = False
if show_shapes: 
    show_num = 50
    for i, (thingsimg, catID) in enumerate(thingsloader): 
        print(thingscats[catID])
        print(thingsimg.shape)
        if i > show_num:
            break

# +
# %%time
#https://www.datacamp.com/community/tutorials/xgboost-in-python

if not exists(regr_data_fn) and reextract_data:
    # CPU times: user 1h 7min 42s, sys: 2min 53s, total: 1h 10min 36s
    # Wall time: 1h 8min 37s

    x_spose_vecs = np.zeros([num_imgs,49])  # X: SPoSE vectors (same for each cat)
    y_clip_vecs = np.zeros([num_imgs,512])  # Y: clip vectors

    # Load data
    img_i = 0
    for i, (thingsimg, catID) in enumerate(thingsloader):

        x_spose_vecs[i,:] = spose_cat_emb[catID,:] 
        y_clip_vecs[i,:] = clip_perceptor.encode_image(thingsimg).detach().numpy().squeeze()

        img_i += 1

    savemat(regr_data_fn, {"x_spose":x_spose_vecs, "y_clip":y_clip_vecs})
else: 
    regr_data = loadmat(regr_data_fn)
    x_spose_vecs = regr_data["x_spose"]
    y_clip_vecs = regr_data["y_clip"]

# +
print("X shape:", x_spose_vecs.shape)
print("Y shape:", y_clip_vecs.shape)

xspose_train, xspose_test, yclip_train, yclip_test = train_test_split( x_spose_vecs, 
                                                                       y_clip_vecs, 
                                                                       test_size=0.10, 
                                                                       random_state=42)

xclip_train, xclip_test, yspose_train, yspose_test = train_test_split( y_clip_vecs, 
                                                                       x_spose_vecs, 
                                                                       test_size=0.10, 
                                                                       random_state=42)
# -

if False: # scaling
    scaler = StandardScaler()
    xspose_train = scaler.fit_transform(xspose_train)
    xspose_test = scaler.fit_transform(xspose_test)
    
    xclip_train = scaler.fit_transform(xclip_train)
    xclip_test = scaler.fit_transform(xclip_test)

# ## xgboost model

# +
# %%time
# on hyperparameter tuning: https://github.com/KSpiliop/Fraud_Detection
# visual guide parameter tuning: 
# https://kevinvecmanis.io/machine%20learning/hyperparameter%20tuning/dataviz/python/2019/05/11/XGBoost-Tuning-Visual-Guide.html
# https://stackoverflow.com/questions/69786993/tuning-xgboost-hyperparameters-with-randomizedsearchcv

model = MultiOutputRegressor( xgb.XGBRegressor(), n_jobs=-1)

grid = {  'n_estimators': [100, 500, 900, 1100, 1500],
          'max_depth': [2, 3, 5, 10, 15],
          'learning_rate': [0.05, 0.1, 0.15, 0.20],
          'min_child_weight': [1, 2, 3, 4]    }

kfold = KFold(n_splits=5, shuffle=True, random_state=10)

grid_search = RandomizedSearchCV(model, param_grid, scoring="r2", n_jobs=-1, verbose=1, n_iter = 500, cv=kfold)
grid_result = grid_search.fit(xspose_train, yclip_train, **param_grid)

# +
# %%time
# on hyperparameter tuning: https://github.com/KSpiliop/Fraud_Detection
# visual guide parameter tuning: 
# https://kevinvecmanis.io/machine%20learning/hyperparameter%20tuning/dataviz/python/2019/05/11/XGBoost-Tuning-Visual-Guide.html
# https://stackoverflow.com/questions/69786993/tuning-xgboost-hyperparameters-with-randomizedsearchcv
# example with multioutput and randomsearchcv: 
# https://stackoverflow.com/questions/43532811/gridsearch-over-multioutputregressor

"""
regressor = xgb.XGBRegressor()
multiregressor = MultiOutputRegressor( estimator=regressor ) 

hyparam_grid = {  'estimator__n_estimators': [100, 500, 900, 1100, 1500],
                  'estimator__max_depth': [2, 3, 5, 10, 15],
                  'estimator__learning_rate': [0.05, 0.1, 0.15, 0.20],
                  'estimator__min_child_weight': [1, 2, 3, 4]   }

# Set up the random search with 4-fold cross validation
random_cv = RandomizedSearchCV( estimator=multiregressor,
                                param_distributions=hyparam_grid,
                                cv=5, n_iter=50, n_jobs=-1, 
                                scoring = 'neg_mean_absolute_error', 
                                verbose = 5, 
                                return_train_score = True,
                                random_state=42 )

random_cv.fit(xspose_train, yclip_train)

print( "Best estimater:\n", random_cv.best_estimator_ )
"""
# -

# %%time
model = MultiOutputRegressor( xgb.XGBRegressor() )
model.fit(xspose_train, yclip_train)



# +
yclip_test_pred = model.predict(xspose_test)

print("Pearson score:", np.mean(calc_corrs(yclip_test, yclip_test_pred)))

print("R2 score:", r2_score(yclip_test, yclip_test_pred, multioutput='variance_weighted'))

print("MSE:", mean_squared_error(yclip_test, yclip_test_pred))
# Pearson score: 0.70394534
# R2 score: 0.5203283862199036
# MSE: 0.052056707074646616
# -

# train on all data
modelfull = MultiOutputRegressor( xgb.XGBRegressor() )
modelfull.fit(x_spose_vecs, y_clip_vecs)

# save all data model
with open('../../xgb_spose_to_clip_defaultmodel.pkl','wb') as f:
    pickle.dump(modelfull,f)

# load for trying
with open('../../xgb_spose_to_clip_defaultmodel.pkl', 'rb') as f:
    modelfull = pickle.load(f)

yclip_test_pred = modelfull.predict(xspose_test)

print("Pearson score:", np.mean(calc_corrs(yclip_test, yclip_test_pred)))

yclip_test_pred[233]



# ## Ridge model

# +
# checking data for issues by doing CLIP-to-SPoSE (which is supposed to work)
model = Ridge()
model.fit(xclip_train, yspose_train)
yspose_test_pred = model.predict(xclip_test)

print("Default sklearn Ridge CLIP-to-SPoSE R2:", r2_score(yspose_test, yspose_test_pred))
# -

# checking data for issues by doing CLIP-to-SPoSE (which is supposed to work)
model = RidgeCV()
clf = model.fit(xclip_train, yspose_train)
print("sklearn CV-Ridge CLIP-to-SPoSE R2:", clf.score(xclip_train, yspose_train) )

# %%time
wt, corr, valphas, bscorrs, valinds = bootstrap_ridge(xspose_train, yclip_train,
                                                      xspose_test, yclip_test, 
                                                      alphas=np.logspace(-100, 100, 60),
                                                      nboots=5,
                                                      chunklen=10, nchunks=15, return_wt=True)

# +
yclip_test_pred = np.dot(xspose_test, wt)

print("Pearson score:", np.mean(calc_corrs(yclip_test, yclip_test_pred)))

print("R2 score:", r2_score(yclip_test, yclip_test_pred, multioutput='variance_weighted'))

print("MSE:", mean_squared_error(yclip_test, yclip_test_pred))
# Pearson score: 0.42788085
# R2 score: 0.23764935063586742
# MSE: 0.08273465283753259
# -

print(np.mean(corr))  # 0.68

# train on all data
wt, corr, valphas, bscorrs, valinds = bootstrap_ridge(x_spose_vecs, y_clip_vecs, 
                                                      xspose_test, yclip_test,
                                                      alphas=np.logspace(-2, 100, 100),
                                                      nboots=5,
                                                      chunklen=10, nchunks=15, return_wt=True)

# ## Model test ground

# +
# checking data for issues by doing CLIP-to-SPoSE (which is supposed to work)
model = LinearRegression()
model.fit(xspose_train, yclip_train)
yclip_test_pred = model.predict(xspose_test)

print("Default LinearRegression SPoSE-to-CLIP R2:", r2_score(yclip_test, yclip_test_pred))

# +
# %%time

# fitting
model = MultiOutputRegressor(xgb.XGBRegressor(objective='reg:squarederror'))
model.fit(xspose_train, yclip_train)
yclip_test_pred = model.predict(xspose_test)

print("Default XGBRegressor SPoSE-to-CLIP R2:", r2_score(yclip_test, yclip_test_pred))
# -

"""
params = {"objective":"reg:linear",'colsample_bytree': 0.3,'learning_rate': 0.1,
                'max_depth': 5, 'al                'max_depth': 5, 'alpha': 10}
pha': 10}

cv_results = xgb.cv(dtrain=data_dmatrix, params=params, nfold=3,
                    num_boost_round=50,early_stopping_rounds=10,metrics="rmse", as_pandas=True, seed=123)
# https://www.datacamp.com/community/tutorials/xgboost-in-python
"""

# +
# checking data for issues by doing CLIP-to-SPoSE (which is supposed to work)
model = Ridge()
model.fit(xspose_train, yclip_train)
yclip_test_pred = model.predict(xspose_test)

print("Default RidgeRegression SPoSE-to-CLIP R2:", r2_score(yclip_test, yclip_test_pred))
# -

wt.shape



np.max(corr)





print(np.mean(corr))  # 0.68

savemat('big_sleep/data/W_aridge_spose_to_clip.mat', {'W':wt} )



# checking data for issues by doing CLIP-to-SPoSE (which is supposed to work)
model = MultiTaskLassoCV()
clf = model.fit(xclip_train, yspose_train)
clf.score(xclip_train, yspose_train)

# checking data for issues by doing CLIP-to-SPoSE (which is supposed to work)
model = MultiTaskElasticNetCV()
clf = model.fit(xclip_train, yspose_train)
clf.score(xclip_train, yspose_train)

# +
# TODO: set up modeling with CLIP-to-spose

# TODO: standard scaler
# TODO: compare regression methods: https://towardsdatascience.com/quickly-test-multiple-models-a98477476f0
# TODO: determine best parameters via CV
# TODO: try alexridge

# +
# TODO: set up modeling with CLIP-to-spose



# +
model = MultiTaskElasticNet()
model.fit(xspose_train, yclip_train)
yclip_test_pred = model.predict(xspose_test)

# TODO: determine best parameters
# -

model = MultiTaskLasso()
model.fit(xspose_train, yclip_train)
yclip_test_pred = model.predict(xspose_test)

# %%time
#yclip_test_pred = MultiOutputRegressor(GradientBoostingRegressor(random_state=0)).fit(xspose_train, yclip_train).predict(xspose_test)

r2_score(yclip_test, yclip_test_pred)

# +
model = LinearRegression()  # TODO: elasticNet
# fit model
model.fit(X, y)
# make a prediction

# define the evaluation procedure
cv = RepeatedKFold(n_splits=10, n_repeats=3, random_state=1)
# evaluate the model and collect the scores
n_scores = cross_val_score(model, X, y, scoring='neg_mean_absolute_error', cv=cv, n_jobs=-1)
# force the scores to be positive
n_scores = absolute(n_scores)
# summarize performance
print('MAE: %.3f (%.3f)' % (mean(n_scores), std(n_scores)))

# -

































StarcoderdataJupyterScriptsDedupFiltered