<jupyter_start><jupyter_text>Catboost regression optmized by Model Goodness  Package import<jupyter_code>import pandas as pd
import seaborn as sns

import sys
sys.path.insert(0,r'C:\Users\eduar\OneDrive\PhD\UTuning')

from UTuning import scorer, plots
from sklearn.model_selection import train_test_split

from catboost import CatBoostRegressor
df=pd.read_csv(r'C:\Users\eduar\OneDrive\PhD\UTuning\dataset\unconv_MV.csv')<jupyter_output><empty_output><jupyter_text>Split train test<jupyter_code>y=df['Production'].values
X=df[['Por','LogPerm','Brittle','TOC']].values
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.33)
print(X_train.shape,y_train.shape)<jupyter_output><empty_output><jupyter_text>Model definition and regression<jupyter_code>model = CatBoostRegressor(iterations=500, learning_rate=0.2, loss_function='RMSEWithUncertainty',
                          verbose=False, random_seed=0)

model.fit(X_train,y_train)

estimates = model.predict(X_test)
estimates.shape
plots.error_line(estimates[:,0],y_test,estimates[:,1],Frac = 1)<jupyter_output><empty_output><jupyter_text>Model evaluation  Scikit-learn hyperparameter optmization<jupyter_code>from sklearn.model_selection import RandomizedSearchCV
from sklearn.metrics import make_scorer
#from sklearn.ensemble import GradientBoostingRegressor
max_depth = list(np.arange(10, 200, step=10))
lr = list(np.arange(0.01, 0.1, step=.01))
param_grid = {
    "learning_rate": lr,
    "n_estimators": max_depth
}
def APG_calc(Truth, Pred, Sigma,n_quantiles):
    mask = np.random.choice([False, True],
                            len(Pred),
                            p=[0, 1]) # To display randomly less points [Remove , Keep] in fraction

    Pred=Pred[mask]
    perc = np.linspace(0.0, 1.00, n_quantiles)

    F = np.zeros(Pred.shape[0])
    Indicator_func = np.zeros((Pred.shape[0], perc.shape[0]))

    # range of symmetric p-probability intervals
    plow = (1 - perc) / 2
    pupp = (1 + perc) / 2
    
    for i in range(len(Pred)):
        F[i] = stats.norm.cdf(Truth,
                          loc=Pred[i],
                          scale=Sigma)
        for proba_low, proba_upp in zip(plow, pupp):
            for k in range(len(plow)):
                if plow[k] < F[i] <= pupp[k]:
                    Indicator_func[i, k] = 1
                else:
                    Indicator_func[i, k] = 0

    avgIndFunc = np.mean(Indicator_func, axis=0)

    a = np.zeros(len(avgIndFunc))
    for i in range(len(avgIndFunc)):
        if avgIndFunc[i] > perc[i] or avgIndFunc[i] == perc[i]:
            a[i] = 1
        else:
            a[i] = 0
    #print('Overall uncertainty = {0:2.2f}'.format(Sigma.mean()))
    U = Sigma.mean()

    Accuracy = integrate.simps(a, perc)

    Prec = a*(avgIndFunc-perc)
    
    Precision = 1-2*integrate.simps(Prec, perc)

    Sum = (3*a-2)*(avgIndFunc-perc)

    Goodness = 1-integrate.simps(Sum, perc)
    
    return Accuracy, Precision, Goodness, avgIndFunc, U


# def my_custom_loss_func(y_true, y_pred):
    
    

#     L = 10
#     mean = np.empty((L, len(perc)))

#     for p_interv in range(len(perc)):
#         for l in np.arange(0, L):
#             samples = random.choices(IF_array[:, p_interv],
#                                      k=IF_array.shape[0])
#             mean[l, p_interv] = np.mean(samples)
    
#     return np.mean(Pred-Truth)

def my_custom_loss_func(y_true, y_pred):
    
    #diff = np.abs(y_true - y_pred[:,0])
    
    n_quantiles=11
    perc = np.linspace(0.0, 1.00, n_quantiles)
    Samples = 10
    
    preds = y_pred

    Pred = preds[:,0]
    Sigma=np.sqrt(preds[:,1])
    Truth = y_true
    
    Pred_array = np.zeros((Sigma.shape[0],Samples))

    A_array=np.zeros(Pred.shape[0])
    P_array=np.zeros(Pred.shape[0])
    G_array=np.zeros(Pred.shape[0])
    U_array=np.zeros(Pred.shape[0])

    IF_array=np.zeros((Pred.shape[0],n_quantiles))

    for i in range(len(Pred)):
        Pred_array[i,:] = np.random.normal(loc=Pred[i],scale=Sigma[i],size=Samples)
        A,P,G,IF,U=APG_calc(Truth[i], Pred_array[i,:], Sigma[i],n_quantiles)
        U_array[i]=U
        A_array[i]=A
        P_array[i]=P
        G_array[i]=G
        IF_array[i,:] = IF

    avgIndFunc = np.mean(IF_array, axis=0)

    print('Accuracy = {0:1.2f}'.format(np.mean(A_array)))
    print('Precision = {0:1.2f}'.format(np.mean(P_array)))
    print('Goodness = {0:1.2f}'.format(np.mean(G_array)))

    
    return np.mean(G_array)
# score = make_scorer(my_custom_loss_func,greater_is_better=True)
# model=CatBoostRegressor(loss_function='RMSEWithUncertainty',
#                         verbose = False)
# clf=model.fit(X_train, y_train)

# #preds = model.predict(X_test)
#preds.shape
# score(clf,X_test,y_test)
# model = CatBoostRegressor(iterations=100,
#                           learning_rate=0.2,
#                           loss_function='RMSEWithUncertainty',
#                           verbose=False,
#                           random_seed=0)
scorer = make_scorer(my_custom_loss_func,greater_is_better=True)

model=CatBoostRegressor(loss_function='RMSEWithUncertainty',
                        verbose=False)
#model = GradientBoostingRegressor()

random_cv = RandomizedSearchCV(model,
                               param_grid,
                               cv=2,
                               n_iter=10,
                               n_jobs=-1,
                               scoring = scorer
)
_ = random_cv.fit(X_train, y_train)

print(random_cv.cv_results_)
import pandas as pd
df = pd.DataFrame(random_cv.cv_results_)
df.info()
random_cv.best_params_<jupyter_output><empty_output><jupyter_text>Functions<jupyter_code>def error_line(Mean,Truth,STD):
    '''
    Simple function to draw an error line plot. 
    It takes three arrays of the same length, the predicted value (Mean), the truth value (Truth) and the standard deviation (STD).
    '''
    xline = [0,12000]#
    yline = [0,12000]#
    plt.figure(figsize=(10, 6))
    mask = np.random.choice([False, True], len(Mean), p=[0.0, 1]) # To display randomly less points [Remove , Keep] in fraction
    plt.errorbar(Mean[mask], Truth[mask], xerr=STD[mask], 
                 fmt='k.',
                 ecolor='k')
    plt.plot(xline, yline, '-k')
    
    plt.xlabel('Predicted value, $\hat{y}$')
    plt.ylabel('True value, $y$ ')

    plt.show()
    
def get_GAP(Pred, Sigma, Truth, n_quantiles):
    '''
    This function takes the dataframe and calculates the indicator function and the average
    we then use this information to estimate the accuracy, goodness and precision of the model.
    '''
    perc = np.linspace(0.0, 1.00, n_quantiles)
    F = np.zeros(Pred.shape[0])
    Indicator_func = np.zeros((Pred.shape[0], perc.shape[0]))

    # range of symmetric p-probability intervals
    plow = (1 - perc) / 2
    pupp = (1 + perc) / 2
        
    for i in range(len(Pred)):
        F[i] = stats.norm.cdf(Truth[i],
                              loc=Pred[i],
                              scale=Sigma[i])
        for proba_low, proba_upp in zip(plow, pupp):
            for k in range(len(plow)):
                if plow[k] < F[i] <= pupp[k]:
                    Indicator_func[i, k] = 1
                else:
                    Indicator_func[i, k] = 0

    avgIndFunc = np.mean(Indicator_func, axis=0)
    
    a = np.zeros(len(avgIndFunc))
    for i in range(len(avgIndFunc)):
        if avgIndFunc[i] > perc[i] or avgIndFunc[i] == perc[i]:
            a[i] = 1
        else:
            a[i] = 0
    print(f'Overall uncertainty = {Sigma.mean():.6f}')
    U = Sigma.mean()
    
    ##% Goodness, Precision and Accuracy
    

    Accuracy = integrate.simps(a, perc)

    Prec = a*(avgIndFunc-perc)
    Precision = 1-2*integrate.simps(Prec, perc)

    Sum = (3*a-2)*(avgIndFunc-perc)
    Goodness = 1-integrate.simps(Sum, perc)

    return Goodness, Precision, Accuracy, U, Indicator_func,perc
    
def accuracy_plot(Indicator_func,perc):
    '''
    This function takes the indicator function and percentiles to draw the accuracy plot.
    '''
    
    L = 100  
    mean = np.empty((L, len(perc)))
    std = np.empty_like(mean)
    avgIndFunc = np.mean(Indicator_func, axis=0)
    for p_interv in range(len(perc)):
        for l in np.arange(0, L):
            samples = random.choices(Indicator_func[:, p_interv],
                                     k=Indicator_func.shape[0])
            mean[l, p_interv] = np.mean(samples)

    plt.figure(figsize=(10, 6))
    plt.plot(perc, avgIndFunc,'-ok',markersize=5)
    plt.plot(perc,np.round(avgIndFunc+np.std(mean, axis=0), 3),'--k')
    plt.plot(perc,np.round(avgIndFunc-np.std(mean, axis=0), 3),'--k')
    plt.plot([0, 1],[0, 1],'-k')
    plt.ylabel(r"$\overline{\xi (p)}$")
    plt.xlabel('Probability interval $p$')
    plt.ylim(0,1)
    plt.xlim(0,1)
    plt.show()

def Error_Acc(Pred,Truth,Sigma,n_quantiles,seed):
    np.random.seed(seed)
    perc = np.linspace(0.0, 1.00, n_quantiles)
    
    mask = np.random.choice([False, True],
                            len(Pred),
                            p=[0.95, 0.05]) # To display randomly less points [Remove , Keep] in fraction

    Truth=Truth[mask]
    Pred=Pred[mask]
    Sigma=Sigma[mask]
    
    F = np.zeros(Pred.shape[0])
    Indicator_func = np.zeros((Pred.shape[0], perc.shape[0]))

    # range of symmetric p-probability intervals
    plow = (1 - perc) / 2
    pupp = (1 + perc) / 2
        
    for i in range(len(Pred)):
        F[i] = stats.norm.cdf(Truth[i],
                              loc=Pred[i],
                              scale=Sigma[i])
        for proba_low, proba_upp in zip(plow, pupp):
            for k in range(len(plow)):
                if plow[k] < F[i] <= pupp[k]:
                    Indicator_func[i, k] = 1
                else:
                    Indicator_func[i, k] = 0

    avgIndFunc = np.mean(Indicator_func, axis=0)
    
    a = np.zeros(len(avgIndFunc))
    for i in range(len(avgIndFunc)):
        if avgIndFunc[i] > perc[i] or avgIndFunc[i] == perc[i]:
            a[i] = 1
        else:
            a[i] = 0
    print('Overall uncertainty = {0:2.2f}'.format(Sigma.mean()))
    U = Sigma.mean()
    
    Accuracy = integrate.simps(a, perc)
    print(a)

    Prec = a*(avgIndFunc-perc)
    print(Prec)
    Precision = 1-2*integrate.simps(Prec, perc)

    Sum = (3*a-2)*(avgIndFunc-perc)
    print(Sum)
    Goodness = 1-integrate.simps(Sum, perc)
    
    print('Accuracy = {0:1.2f}'.format(Accuracy))
    print('Precision = {0:1.2f}'.format(Precision))
    print('Goodness = {0:1.2f}'.format(Goodness))
    
    L = 100 
    mean = np.empty((L, len(perc)))

    for p_interv in range(len(perc)):
        for l in np.arange(0, L):
            samples = random.choices(Indicator_func[:, p_interv],
                                     k=Indicator_func.shape[0])
            mean[l, p_interv] = np.mean(samples)

    fig,(ax1,ax2,ax3)=plt.subplots(1,3,figsize=(16,4))
    
    xline = [0,max(Pred.max(),Truth.max())+max(Pred.max(),Truth.max())*0.1]#
    yline = [0,xline[1]]#

    ax1.errorbar(Pred, Truth, xerr=Sigma, 
                 fmt='k.',
                 ecolor='k')
    ax1.plot(xline, yline, '-k')
    ax1.set_xlabel('Predicted value, $\hat{y}$')
    ax1.set_ylabel('True value, $y$ ')
    
    ax2.plot(perc, avgIndFunc,'-ok',markersize=5)
    ax2.plot(perc,np.round(avgIndFunc+np.std(mean, axis=0), 3),'--k')
    ax2.plot(perc,np.round(avgIndFunc-np.std(mean, axis=0), 3),'--k')
    ax2.plot([0, 1],[0, 1],'-k')
    ax2.set_ylabel(r"$\overline{\xi (p)}$")
    ax2.set_xlabel('Probability interval $p$')
    ax2.set_ylim(0,1)
    ax2.set_xlim(0,1)
    
    ax2.plot(perc, avgIndFunc,'-ok',markersize=5)

def Error_Acc_All(Pred,Truth,Sigma,n_quantiles):
    perc = np.linspace(0.0, 1.00, n_quantiles)
    
    mask = np.random.choice([False, True],
                            len(Pred),
                            p=[0, 1]) # To display randomly less points [Remove , Keep] in fraction

    Truth=Truth[mask]
    Pred=Pred[mask]
    Sigma=Sigma[mask]
    
    F = np.zeros(Pred.shape[0])
    Indicator_func = np.zeros((Pred.shape[0], perc.shape[0]))

    # range of symmetric p-probability intervals
    plow = (1 - perc) / 2
    pupp = (1 + perc) / 2
        
    for i in range(len(Pred)):
        F[i] = stats.norm.cdf(Truth[i],
                              loc=Pred[i],
                              scale=Sigma[i])
        for proba_low, proba_upp in zip(plow, pupp):
            for k in range(len(plow)):
                if plow[k] < F[i] <= pupp[k]:
                    Indicator_func[i, k] = 1
                else:
                    Indicator_func[i, k] = 0

    avgIndFunc = np.mean(Indicator_func, axis=0)
    
    a = np.zeros(len(avgIndFunc))
    for i in range(len(avgIndFunc)):
        if avgIndFunc[i] > perc[i] or avgIndFunc[i] == perc[i]:
            a[i] = 1
        else:
            a[i] = 0
    print('Overall uncertainty = {0:2.2f}'.format(Sigma.mean()))
    U = Sigma.mean()
    
    Accuracy = integrate.simps(a, perc)
    print(a)

    Prec = a*(avgIndFunc-perc)
    print(Prec)
    Precision = 1-2*integrate.simps(Prec, perc)

    Sum = (3*a-2)*(avgIndFunc-perc)
    print(Sum)
    Goodness = 1-integrate.simps(Sum, perc)
    
    print('Accuracy = {0:1.2f}'.format(Accuracy))
    print('Precision = {0:1.2f}'.format(Precision))
    print('Goodness = {0:1.2f}'.format(Goodness))
    
    L = 100 
    mean = np.empty((L, len(perc)))

    for p_interv in range(len(perc)):
        for l in np.arange(0, L):
            samples = random.choices(Indicator_func[:, p_interv],
                                     k=Indicator_func.shape[0])
            mean[l, p_interv] = np.mean(samples)

    fig,(ax1,ax2,ax3)=plt.subplots(1,3,figsize=(16,4))
    
    xline = [0,max(Pred.max(),Truth.max())+max(Pred.max(),Truth.max())*0.1]#
    yline = [0,xline[1]]#

    ax1.errorbar(Pred, Truth, xerr=Sigma, 
                 fmt='k.',
                 ecolor='k')
    ax1.plot(xline, yline, '-k')
    ax1.set_xlabel('Predicted value, $\hat{y}$')
    ax1.set_ylabel('True value, $y$ ')
    
    ax2.plot(perc, avgIndFunc,'-ok',markersize=5)
    ax2.plot(perc,np.round(avgIndFunc+np.std(mean, axis=0), 3),'--k')
    ax2.plot(perc,np.round(avgIndFunc-np.std(mean, axis=0), 3),'--k')
    ax2.plot([0, 1],[0, 1],'-k')
    ax2.set_ylabel(r"$\overline{\xi (p)}$")
    ax2.set_xlabel('Probability interval $p$')
    ax2.set_ylim(0,1)
    ax2.set_xlim(0,1)
    
    ax2.plot(perc, avgIndFunc,'-ok',markersize=5)
    
def histogram(mc_predictions):
    '''
    From the Monte Carlo predictions we draw a random point and construct the
    histogram of predictions from the model
    '''
    
    Avg = []
    rand=np.random.randint(0,mc_predictions.shape[1])
    for i in range(mc_predictions.shape[0]):
        Avg.append(np.average(mc_predictions[i,rand]))

    Std = np.std(Avg)

    # Histograms
    n_bins = 20
    fig, axs = plt.subplots(1, 1, figsize=(10, 6))
    N, bins, patches = axs.hist(Avg,
                                bins=n_bins,
                                label='$\sigma$ = %2.5f' % Std)
    #axs.set_title('Root Mean squared error in barrels for each cell');
    fracs = N / N.max()
    norm = colors.Normalize(fracs.min(), fracs.max())
    for thisfrac, thispatch in zip(fracs, patches):
        color = plt.cm.binary(norm(thisfrac))
        thispatch.set_facecolor(color)
    plt.legend()
    plt.ylabel('Number of cases')
    plt.xlabel('MSE in predicted value')

def APG_calc(Truth, Pred, Sigma,n_quantiles):

    mask = np.random.choice([False, True],
                            len(Pred),
                            p=[0, 1]) # To display randomly less points [Remove , Keep] in fraction

    #Truth=Truth[mask]
    Pred=Pred[mask]
    #Sigma=Sigma[mask]


    #n_quantiles = 11

    perc = np.linspace(0.0, 1.00, n_quantiles)

    F = np.zeros(Pred.shape[0])
    Indicator_func = np.zeros((Pred.shape[0], perc.shape[0]))

    # range of symmetric p-probability intervals
    plow = (1 - perc) / 2
    pupp = (1 + perc) / 2
    
    for i in range(len(Pred)):
        F[i] = stats.norm.cdf(Truth,
                          loc=Pred[i],
                          scale=Sigma)
        for proba_low, proba_upp in zip(plow, pupp):
            for k in range(len(plow)):
                if plow[k] < F[i] <= pupp[k]:
                    Indicator_func[i, k] = 1
                else:
                    Indicator_func[i, k] = 0

    avgIndFunc = np.mean(Indicator_func, axis=0)

    a = np.zeros(len(avgIndFunc))
    for i in range(len(avgIndFunc)):
        if avgIndFunc[i] > perc[i] or avgIndFunc[i] == perc[i]:
            a[i] = 1
        else:
            a[i] = 0
    #print('Overall uncertainty = {0:2.2f}'.format(Sigma.mean()))
    U = Sigma.mean()

    Accuracy = integrate.simps(a, perc)

    Prec = a*(avgIndFunc-perc)
    
    Precision = 1-2*integrate.simps(Prec, perc)

    Sum = (3*a-2)*(avgIndFunc-perc)

    Goodness = 1-integrate.simps(Sum, perc)

    # print('Accuracy = {0:1.2f}'.format(Accuracy))
    # print('Precision = {0:1.2f}'.format(Precision))
    # print('Goodness = {0:1.2f}'.format(Goodness))
    
    return Accuracy, Precision, Goodness, avgIndFunc, U<jupyter_output><empty_output><jupyter_text>file import<jupyter_code>df.info()
from sklearn.feature_selection import mutual_info_regression

x = df.iloc[:,[1,2,3,4,5,6]]              # separate DataFrames for predictor and response features
y = df.iloc[:,[7]]

mi = mutual_info_regression(x,np.ravel(y)) # calculate mutual information
mi /= np.max(mi)                          # calculate relative mutual information

indices = np.argsort(mi)[::-1]            # find indicies for descending order

print("Feature ranking:")                 # write out the feature importances
for f in range(x.shape[1]):
    print("%d. feature %s = %f" % (f + 1, x.columns[indices][f], mi[indices[f]]))

plt.subplot(111)                          # plot the relative mutual information 
plt.title("Mutual Information")
plt.bar(range(x.shape[1]), mi[indices],
       color="g", align="center")
plt.xticks(range(x.shape[1]), x.columns[indices],rotation=90)
plt.xlim([-1, x.shape[1]])
plt.subplots_adjust(left=0.0, bottom=0.0, right=1., top=1., wspace=0.2, hspace=0.2)
plt.show()<jupyter_output><empty_output><jupyter_text>Model definition and regression<jupyter_code># # predict mean value and data uncertainty

# model = CatBoostRegressor(iterations=100, learning_rate=0.2, loss_function='RMSEWithUncertainty',
#                           verbose=False, random_seed=0)
#                          #task_type = '0:1')

# #train_pool=Pool(X_train,y_train)

# model.fit(X_train,y_train)

# preds = model.predict(X_test)
# print(preds.shape)<jupyter_output><empty_output><jupyter_text>Model evaluation<jupyter_code># np.random.seed(0)

# n_quantiles=11

# perc = np.linspace(0.0, 1.00, n_quantiles)

# Samples = 100

# Sigma=np.sqrt(preds[:,1])

# Pred = preds[:,0]
# Truth = y_test

# Pred_array = np.zeros((Sigma.shape[0],Samples))

# A_array=np.zeros(Pred.shape[0])
# P_array=np.zeros(Pred.shape[0])
# G_array=np.zeros(Pred.shape[0])
# U_array=np.zeros(Pred.shape[0])

# IF_array=np.zeros((Pred.shape[0],n_quantiles))

# for i in range(len(Pred)):
#     Pred_array[i,:] = np.random.normal(loc=Pred[i],scale=Sigma[i],size=Samples)
#     A,P,G,IF,U=APG_calc(Truth[i], Pred_array[i,:], Sigma[i],n_quantiles)
#     U_array[i]=U
#     A_array[i]=A
#     P_array[i]=P
#     G_array[i]=G
#     IF_array[i,:] = IF

# avgIndFunc = np.mean(IF_array, axis=0)

# print('Accuracy = {0:1.2f}'.format(np.mean(A_array)))
# print('Precision = {0:1.2f}'.format(np.mean(P_array)))
# print('Goodness = {0:1.2f}'.format(np.mean(G_array)))

# L = 10
# mean = np.empty((L, len(perc)))

# for p_interv in range(len(perc)):
#     for l in np.arange(0, L):
#         samples = random.choices(IF_array[:, p_interv],
#                                  k=IF_array.shape[0])
#         mean[l, p_interv] = np.mean(samples)
# fig,(ax1,ax2)=plt.subplots(1,2,figsize=(12,4))

# xline = [0,max(Pred.max(),Truth.max())+max(Pred.max(),Truth.max())*0.1]#
# yline = [0,xline[1]]#

# ax1.errorbar(Pred, Truth, xerr=Sigma, 
#              fmt='k.',
#              ecolor='k')
# ax1.plot(xline, yline, '-k')
# ax1.set_xlabel('Predicted value, $\hat{y}$')
# ax1.set_ylabel('True value, $y$ ')

# ax2.plot(perc, avgIndFunc,'-ok',markersize=5)
# ax2.plot(perc,np.round(avgIndFunc+np.std(mean, axis=0), 3),'--k')
# ax2.plot(perc,np.round(avgIndFunc-np.std(mean, axis=0), 3),'--k')
# ax2.plot([0, 1],[0, 1],'-k')
# ax2.set_ylabel(r"$\overline{\xi (p)}$")
# ax2.set_xlabel('Probability interval $p$')
# ax2.set_ylim(0,1)
# ax2.set_xlim(0,1)

# ax2.plot(perc, avgIndFunc,'-ok',markersize=5)<jupyter_output><empty_output><jupyter_text>Virtual ensemble<jupyter_code>def virt_ensemble(X_train,y_train, num_samples=100, iters=1000, lr=0.2):
    ens_preds = []
    
    model = CatBoostRegressor(iterations=iters, learning_rate=lr, loss_function='RMSEWithUncertainty',
                          verbose=False, random_seed=0)



    model.fit(X_train,y_train)

    
    ens_preds = model.virtual_ensembles_predict(X_test, prediction_type='VirtEnsembles', 
                                                virtual_ensembles_count=num_samples,
                                                thread_count=8)
    return np.asarray(ens_preds)
np.random.seed(0)

#Pred_array = np.zeros((Sigma.shape[0],Samples))

n_quantiles=11

perc = np.linspace(0.0, 1.00, n_quantiles)

Samples = 10

ens_preds=virt_ensemble(X_train,y_train, num_samples=Samples)

Pred_array = ens_preds[:,:,0]
# print(Pred_array[:5,1])
# print(Truth[:5])
# #plt.scatter(np.average(Pred_array,axis=1),Truth)
# plt.scatter(Pred_array[:,0],Truth)
# plt.scatter(Pred_array[:,1],Truth)
# plt.scatter(Pred_array[:,2],Truth)
# plt.scatter(Pred_array[:,3],Truth)
# plt.scatter(Pred_array[:,9],Truth)
Knowledge_u=np.sqrt(np.var(Pred_array,axis=1)) #Knowledge uncertainty
Data_u=np.sqrt(np.mean(ens_preds[:,:,1],axis=1)) #Data uncertainty
Sigma=Knowledge_u+Data_u

Truth = y_test

A_array=np.zeros(Pred_array.shape[0])
P_array=np.zeros(Pred_array.shape[0])
G_array=np.zeros(Pred_array.shape[0])
U_array=np.zeros(Pred_array.shape[0])

IF_array=np.zeros((Pred_array.shape[0],n_quantiles))

for i in range(Pred_array.shape[0]):
    #Pred_array[i,:] = np.random.normal(loc=Pred[i],scale=Sigma[i],size=Samples)
    A,P,G,IF,U=APG_calc(Truth[i], Pred_array[i,:], Sigma[i],n_quantiles)
    U_array[i]=U
    A_array[i]=A
    P_array[i]=P
    G_array[i]=G
    IF_array[i,:] = IF

avgIndFunc = np.mean(IF_array, axis=0)

print('Accuracy = {0:1.2f}'.format(np.mean(A_array)))
print('Precision = {0:1.2f}'.format(np.mean(P_array)))
print('Goodness = {0:1.2f}'.format(np.mean(G_array)))

L = 10
mean = np.empty((L, len(perc)))

for p_interv in range(len(perc)):
    for l in np.arange(0, L):
        samples = random.choices(IF_array[:, p_interv],
                                 k=IF_array.shape[0])
        mean[l, p_interv] = np.mean(samples)
fig,(ax1,ax2)=plt.subplots(1,2,figsize=(12,4))

xline = [0,max(np.mean(Pred_array,axis=1).max(),Truth.max())+max(np.mean(Pred_array,axis=1).max(),Truth.max())*0.1]#
yline = [0,xline[1]]#

ax1.errorbar(np.mean(Pred_array,axis=1), Truth, xerr=Sigma, 
             fmt='k.',
             ecolor='k')
ax1.plot(xline, yline, '-k')
ax1.set_xlabel('Predicted value, $\hat{y}$')
ax1.set_ylabel('True value, $y$ ')

ax2.plot(perc, avgIndFunc,'-ok',markersize=5)
ax2.plot(perc,np.round(avgIndFunc+np.std(mean, axis=0), 3),'--k')
ax2.plot(perc,np.round(avgIndFunc-np.std(mean, axis=0), 3),'--k')
ax2.plot([0, 1],[0, 1],'-k')
ax2.set_ylabel(r"$\overline{\xi (p)}$")
ax2.set_xlabel('Probability interval $p$')
ax2.set_ylim(0,1)
ax2.set_xlim(0,1)

ax2.plot(perc, avgIndFunc,'-ok',markersize=5)<jupyter_output><empty_output><jupyter_text>Optimization<jupyter_code>value='value'
xname='Tnumber'
yname='Lrate'
zname='Tdepth'
def virt_ensemble(X_train,y_train, num_samples=100, iters=1000, lr=0.2, depth=100):
    ens_preds = []
    
    model = CatBoostRegressor(iterations=iters,
                              learning_rate=lr,
                              depth=depth,
                              loss_function='RMSEWithUncertainty',
                              verbose=False)



    model.fit(X_train,y_train)

    
    ens_preds = model.virtual_ensembles_predict(X_test, prediction_type='VirtEnsembles', 
                                                virtual_ensembles_count=num_samples,
                                                thread_count=8)
    return np.asarray(ens_preds)

def objective(trial):
    iterations = trial.suggest_float("{0}".format(xname), 50, 1000)
    lrate = trial.suggest_float("{0}".format(yname), 0.001, 0.2)
    depth = trial.suggest_float("{0}".format(zname), 4, 16)
    
    #np.random.seed(0)
    

    n_quantiles=11

    perc = np.linspace(0.0, 1.00, n_quantiles)

    Samples = 100

    ens_preds=virt_ensemble(X_train,
                            y_train,
                            num_samples=Samples,
                            iters=int(iterations),
                            lr=lrate,
                            depth = int(depth))

    Pred_array = ens_preds[:,:,0]
    
    Knowledge_u=np.sqrt(np.var(Pred_array,axis=1)) #Knowledge uncertainty
    Data_u=np.sqrt(np.mean(ens_preds[:,:,1],axis=1)) #Data uncertainty
    Sigma=Knowledge_u+Data_u
    
    #Pred_array = np.zeros((Sigma.shape[0],Samples))
    Truth = y_test
    
    
    A_array=np.zeros(Pred_array.shape[0])
    P_array=np.zeros(Pred_array.shape[0])
    G_array=np.zeros(Pred_array.shape[0])
    U_array=np.zeros(Pred_array.shape[0])

    IF_array=np.zeros((Pred_array.shape[0],n_quantiles))

    for i in range(Pred_array.shape[0]):
        A,P,G,IF,U=APG_calc(Truth[i], Pred_array[i,:], Sigma[i],n_quantiles)
        U_array[i]=U
        A_array[i]=A
        P_array[i]=P
        G_array[i]=G
        IF_array[i,:] = IF

    avgIndFunc = np.mean(IF_array, axis=0)

    print('Accuracy = {0:1.2f}'.format(np.mean(A_array)))
    print('Precision = {0:1.2f}'.format(np.mean(P_array)))
    print('Goodness = {0:1.2f}'.format(np.mean(G_array)))

    L = 100
    mean = np.empty((L, len(perc)))

    for p_interv in range(len(perc)):
        for l in np.arange(0, L):
            samples = random.choices(IF_array[:, p_interv],
                                     k=IF_array.shape[0])
            mean[l, p_interv] = np.mean(samples)
            
    fig,(ax1,ax2)=plt.subplots(1,2,figsize=(12,4))

    xline = [0,max(np.mean(Pred_array,axis=1).max(),Truth.max())+max(np.mean(Pred_array,axis=1).max(),Truth.max())*0.1]#
    yline = [0,xline[1]]#

    ax1.errorbar(np.mean(Pred_array,axis=1), Truth, xerr=Sigma, 
                 fmt='k.',
                 ecolor='k')
    ax1.plot(xline, yline, '-k')
    ax1.set_xlabel('Predicted value, $\hat{y}$')
    ax1.set_ylabel('True value, $y$ ')

    ax2.plot(perc, avgIndFunc,'-ok',markersize=5)
    ax2.plot(perc,np.round(avgIndFunc+np.std(mean, axis=0), 3),'--k')
    ax2.plot(perc,np.round(avgIndFunc-np.std(mean, axis=0), 3),'--k')
    ax2.plot([0, 1],[0, 1],'-k')
    ax2.set_ylabel(r"$\overline{\xi (p)}$")
    ax2.set_xlabel('Probability interval $p$')
    ax2.set_ylim(0,1)
    ax2.set_xlim(0,1)

    ax2.plot(perc, avgIndFunc,'-ok',markersize=5)
    plt.show()
    
    return np.mean(G_array)
study = optuna.create_study(direction= 'maximize')

study.optimize(objective,
               n_trials=100,
               show_progress_bar=True)
StudyName = 'Maximize_Random_G'

joblib.dump(study, StudyName+'.pkl')

fig = optuna.visualization.plot_contour(study, params = ['Tnumber','Lrate'])
fig.show()
fig = optuna.visualization.plot_contour(study, params = ['Tnumber','Tdepth'])
fig.show()
fig = optuna.visualization.plot_contour(study, params = ['Lrate','Tdepth'])
fig.show()
fig = optuna.visualization.plot_contour(study, params = ['Tnumber','Lrate','Tdepth'])
fig.show()
x = np.arange(50,1000,10,dtype=float)
y = np.linspace(.001,.2,20,dtype=float)
z = np.arange(4,16,1,dtype=float)

search_space = {"{0}".format(xname): x,
                "{0}".format(yname): y,
                "{0}".format(zname): z}

study = optuna.create_study(sampler=optuna.samplers.GridSampler(search_space),
                            direction= 'maximize')

study.optimize(objective,
               #n_trials=K_space.shape[0] * F_space.shape[0],
               n_trials=len(x)*len(y)*len(z),
               show_progress_bar=True)
StudyName = 'example'

joblib.dump(study, StudyName+'.pkl')

study = joblib.load(StudyName+'.pkl')

study.trials_dataframe()

df = study.trials_dataframe().drop(['state',
                                    'datetime_start',
                                    'datetime_complete',
                                    'duration',
                                    'system_attrs_grid_id',
                                    'system_attrs_search_space',
                                    'state'], axis=1)
df.info()
import scipy.ndimage
import matplotlib.tri as tri
from scipy.ndimage.filters import gaussian_filter

z=df['{0}'.format(value)].values
x=df['params_{0}'.format(xname)].values
y=df['params_{0}'.format(yname)].values

fig, (ax1) = plt.subplots(nrows=1,figsize=(12,6))

npoints=77
smooth=4

# Create grid values first.
xi = np.linspace(x.min(), x.max(), npoints)
yi = np.linspace(y.min(), y.max(), npoints)

# Linearly interpolate the data (x, y) on a grid defined by (xi, yi).
triang = tri.Triangulation(x, y)
interpolator = tri.LinearTriInterpolator(triang, z)
Xi, Yi = np.meshgrid(xi, yi)
zi = interpolator(Xi, Yi)

zi = gaussian_filter(zi, smooth)

levels=10

ax1.contour(xi, yi, zi, levels=levels, linewidths=0.1, colors='k')
cntr1 = ax1.contourf(xi, yi, zi, levels=levels, cmap="inferno",alpha=0.99)

cbar = plt.colorbar(cntr1, ax=ax1)

cbar.set_label('ErrorValue', rotation=270,labelpad=30)

ax1.set(xlim=(df['params_{0}'.format(xname)].min(),
              df['params_{0}'.format(xname)].max()),
        ylim=(df['params_{0}'.format(yname)].min(),
              df['params_{0}'.format(yname)].max()))
ax1.scatter(x,y,s=3,color='darkgray')
ax1.set_ylabel('{0}'.format(xname))
ax1.set_xlabel('{0}'.format(yname))
plt.savefig("{0}.png".format(StudyName), dpi=600,bbox_inches='tight')

plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16279759941963992, 'fraction_numerical': 0.025093978764096813, 'mean_word_length': 2.7366929521931986, 'pattern_counts': {'":': 2, '<': 51, '<?xml version=': 0, '>': 46, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 21, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '100984', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 12066, 'n_tokens_neox': 11489, 'n_words': 1797, 'path': 'examples/Catboost_regression-ScikitTunning.ipynb', 'repo_name': 'emaldonadocruz/UTuning', 'size': '30326'}