RF_ITEP_2L.py

print('----------- RF for the inverse Transmission Eigenvalue problem for picewise constant refractive index -------------------')
print('----------- N. Pallikarakis and A. Ntargraras - https://arxiv.org/abs/2212.04279 ----------------')
print('----------- Copyright (C) 2023 N. Pallikarakis ----------------------------------')

import numpy as np
import matplotlib.pyplot as plt

###########function for regression chart
def chart_regression(pred, y, sort=False):
    import pandas as pd
    import matplotlib.pyplot as plt
    t = pd.DataFrame({'pred': pred, 'y': y.flatten()})
    if sort:
        t.sort_values(by=['y'], inplace=True)
    plt.plot(t['y'].tolist(), marker='o',label='expected')
    plt.plot(t['pred'].tolist(), marker='o',label='prediction')
    plt.ylabel('output RF')
    plt.legend()
    plt.show()


#############function for perturbation rank
def perturbation_rank(model, x, y, names, regression):
    errors = []
    from sklearn import metrics
    import pandas as pd

    for i in range(x.shape[1]):
        hold = np.array(x[:, i])
        np.random.shuffle(x[:, i])

        if regression:
            pred = model.predict(x)
            error = metrics.mean_squared_error(y, pred)
        else:
            pred = model.predict_proba(x)
            error = metrics.log_loss(y, pred)

        errors.append(error)
        x[:, i] = hold

    max_error = np.max(errors)
    importance = [e / max_error for e in errors]

    data = {'name': names, 'error': errors, 'importance': importance}
    result = pd.DataFrame(data, columns=['name', 'error', 'import'
                                                          'ance'])
    result.sort_values(by=['importance'], ascending=[0], inplace=True)
    result.reset_index(inplace=True, drop=True)
    return result

######### load data
#train data of the direct transmission eigenvalue problem, from the spectral-galerkin method
train_data = np.loadtxt('classic_train.csv', delimiter=',',skiprows=1)
#original data coming from separation of variables in discs
real_data = np.loadtxt('classic_original_10.csv', delimiter=',',skiprows=1) ##ok

print(train_data.shape)
print(real_data.shape)

########## random shuffle of the train data
np.random.seed(100)
#np.random.seed(5)
#np.random.seed(105)
#np.random.seed(205)
#np.random.seed(305)
#np.random.seed(405)
#np.random.seed(505)
#np.random.seed(605)
#np.random.seed(705)
#np.random.seed(805)
#np.random.seed(10) #in place of seed 705

#np.random.seed(10)
np.random.shuffle(train_data)


######## Split them into dependent and independent variables
### training data#####
# the lowest 6 real eigenvalues
X_tr = train_data[:,1:7] #classic & mod
# the refractive index n1,n2 and discontinuiity d1
Y_tr = train_data[:,10:13] #classic


X_test = real_data[:,0:6]
Y_test = real_data[:,6:9]

print(Y_tr.shape)
print(X_tr.shape)
print(X_test.shape)
print(X_test.shape)


##### Prerocessing Step

from sklearn import preprocessing, metrics
from sklearn.model_selection import train_test_split

# Splitting train data into train and valid(test) set to scale them avoiding leaking of data
# The valid(test) set is held back in order to provide unbiased evaluation during a models hyperparameter tuning
X_train, X_val, y_train, y_val = train_test_split(X_tr, Y_tr, test_size=0.3, random_state=10)


#### Standardization, or mean removal and variance scaling
sscaler_X = preprocessing.StandardScaler()
X_train_ss = sscaler_X.fit_transform(X_train)
X_val_ss = sscaler_X.transform(X_val)
X_test_ss = sscaler_X.transform(X_test)


######################### RF ###################


X = X_train_ss
y= y_train



from sklearn.ensemble import RandomForestRegressor
rdregressor=RandomForestRegressor(n_estimators = 200,random_state = 42)

model = rdregressor.fit(X,y)
Y_predict = rdregressor.predict(X_test_ss)
Y_predict_tr = rdregressor.predict(X)
Y_predict_val = rdregressor.predict(X_val_ss)

#print('y_pred',Y_predict)
#print('y_original',Y_test)

#np.savetxt('results_2L_RF_ss.csv', Y_predict,delimiter=',')

########print predicted y vs real y
chart_regression(Y_predict[:,0].flatten(), Y_test[:,0])
chart_regression(Y_predict[:,1].flatten(), Y_test[:,1])
chart_regression(Y_predict[:,2].flatten(), Y_test[:,2])

#### R squared #################
print('----------- model scores -------------------')
print('score training set', round(rdregressor.score(X, y)*100, 2))
print('score validating set', round(rdregressor.score(X_val_ss, y_val)*100, 2))
print('score testing set', round(rdregressor.score(X_test_ss, Y_test)*100, 2))
print('----------- r2 scores -------------------')
from sklearn.metrics import r2_score
print('R squared training set r2', round(metrics.r2_score(y_train, Y_predict_tr)*100, 2))
print('R squared validating set r2', round(metrics.r2_score(y_val, Y_predict_val)*100, 2))
print('R squared testing set r2', round(metrics.r2_score(Y_test, Y_predict)*100, 2))
from sklearn.metrics import r2_score

print('----------- r2 score individually -------------------')
#print('R squared real y set', round(r2_score(Y_re,Y_predict, )*100,2))
print('R squared real n1 set', round(r2_score(Y_test[:,0],Y_predict[:,0])*100,2))
print('R squared real n2 set', round(r2_score(Y_test[:,1],Y_predict[:,1])*100,2))
print('R squared real r1 set', round(r2_score( Y_test[:,2],Y_predict[:,2])*100,2))


###########ranking of eigenvalues
#"""
print('----------- rankings -------------------')
columns = X_tr.shape[1]
names = [i for i in range(1,columns+1)]

#names = list(['k1', 'k2', 'k3', 'k4', 'k5', 'k6'])
#train set
rank = perturbation_rank(model, X, y, names, True)
print(rank)
#np.savetxt('rank_RF_2L_ss.csv', rank,delimiter=',')
#"""