DF = pd.read_table('../regression/classical/models_rrab.dat', sep='\t')
DF.drop(columns=['Mode'], inplace=True)
DF = DF[DF['M'] < 9]
pd.set_option('display.max_columns', 500)

DF['Period'] = np.log10(DF['Period'])
DF.rename(columns={'Period':'logP'}, inplace=True)

DF['W']  = DF['I_M0'] - 1.55 * (DF['V_M0'] - DF['I_M0'])
DF['VI'] = DF['V_M0'] - DF['I_M0']

DF['logR'] = np.log10(np.sqrt(10**DF['logL']) * (DF['Teff']/5777)**-2)

# five models apparently failed to produce correct periods; let's drop them 
DF = DF.drop(np.argmin(DF['logP']))
DF = DF.drop(np.argmin(DF['logP']))
DF = DF.drop(np.argmin(DF['logP']))
DF = DF.reset_index(drop=True)

DF = DF[np.logical_not(np.logical_and(DF['logP'] > 0.1, DF['logR'] < 0.8))]
DF = DF.reset_index(drop=True)

print(DF.columns)
DF

plt.plot(DF['Z'], DF['Y'], '.')
plt.xlabel('Z')
plt.ylabel('Y')
plt.show()

plt.suptitle('RRab H–R Diagram', y=0.99)
plt.scatter(DF['Teff'], DF['logL'], c=DF['M'], cmap='viridis')
plt.xlabel('$T_{\mathrm{eff}}/K$')
plt.ylabel('$\log\;L/$L$_\odot$')
plt.colorbar(label='$M/$M$_\odot$')
plt.xlim([7100, 4900])

plt.savefig('RRL-HR.pdf', bbox_inches='tight')
plt.show()

plt.suptitle('RRab P–M Relation', y=0.99)
plt.scatter(DF['logP'], DF['M'], c=DF['logR'], cmap='viridis')
plt.xlabel('$\log\;P/$days')
plt.ylabel('$M/$M$_\odot$')
plt.colorbar(label='$\log\;R/$R$_\odot$')

plt.savefig('RRL-PMR.pdf', bbox_inches='tight')
plt.show()

plt.suptitle('RRab P–R Relation')
plt.scatter(DF['logP'], DF['logR'], c=DF['M'], cmap='viridis')
plt.ylabel('$\log\;R/$R$_\odot$')
plt.xlabel('$\log\;P/$days')
plt.colorbar(label='$M/$M$_\odot$')

plt.savefig('RRL-PRM.pdf', bbox_inches='tight')
plt.show()

plt.suptitle('RRab P–L Relation')
plt.scatter(DF['logP'], DF['logL'], c=DF['logR'], cmap='viridis')
plt.ylabel('$\log\;L/$L$_\odot$')
plt.xlabel('$\log\;P/$days')
plt.colorbar(label='$\log\;R/$R$_\odot$')

plt.savefig('PL.pdf', bbox_inches='tight')
plt.show()

plt.scatter(DF['I_A2'], DF['Teff'], c=DF['logP'])
plt.xlabel('$I$-band $A_2$')
plt.ylabel('$T_{\\rm{eff}}/$K')
plt.colorbar(label='$\log\;P/$days')

plt.savefig('A2Teff.pdf', bbox_inches='tight')
plt.show()

plt.scatter(DF['V_A2'], DF['Teff'], c=DF['logP'])
plt.xlabel('$V$-band $A_2$')
plt.ylabel('$T_{\\rm{eff}}/$K')
plt.colorbar(label='$\log\;P/$days')

plt.show()

for penalty in [0, 0.00001, 0.0001, 0.001, 0.01]:
    ml = TransformedTargetRegressor(regressor=Pipeline(
        [('scaler', RobustScaler()), 
         ('mlp', MLPRegressor(hidden_layer_sizes=(100, 100), 
                              max_iter=2000, activation='relu', solver='lbfgs', tol=1e-7,
                              alpha=penalty))]),
         transformer=RobustScaler())
    print(penalty, cross_val_score(ml, DF[y_columns], DF[X_columns], cv=KFold(n_splits=2, shuffle=True, random_state=0)))

for idx, name in enumerate(X_columns): 
    print(name, int(np.round(((lm_err[idx,:] >= ml_err[idx,:]).sum()) / lm_err.shape[1] * 100, 0))) 

labels = ['$M$', '$\log\;R$', '$\log\;L$', '$T_{\\rm{eff}}$', '$M_I$', '$M_V$', '$V-I$', '$W$']#, 'Y', 'Z']

plt.rcParams["font.size"] = 18
x = np.arange(len(labels))  # the label locations
width = 0.45  # the width of the bars

fig, ax = plt.subplots(figsize=(6.64306*2,4.10318716*0.9))
plt.suptitle('RR Lyrae', y=0.99)
rects1 = ax.bar(x - width/2, all_r2s[0,:], width, label='Machine learning')
rects2 = ax.bar(x + width/2, all_r2s[1,:], width, label='Linear model')

plt.axhline(1, c='gray', linewidth=1.25, linestyle='--')
plt.axhline(0, c='gray', linewidth=1.25, linestyle='--')

ax.set_ylabel('$R^2$')
ax.set_xticks(x)
ax.set_xticklabels(labels)
ax.set_ylim([-0.05, 1.15])
ax.legend(loc='lower right')

def autolabel(rects):
    """Attach a text label above each bar in *rects*, displaying its height."""
    for rect in rects:
        height = rect.get_height()
        ax.annotate('%.3f' % height,
                    xy=(rect.get_x() + rect.get_width() / 2, height),
                    xytext=(0, 3),  # 3 points vertical offset
                    textcoords="offset points",
                    ha='center', va='bottom')

plt.rcParams["font.size"] = 12
autolabel(rects1)
autolabel(rects2)
plt.rcParams["font.size"] = 18
#autolabel(rects3)

plt.savefig('RRL-R2-ML-LM.pdf')
plt.show()

labels = ['$M$', '$\log\;R$', '$\log\;L$', '$T_{\\rm{eff}}$', '$M_I$', '$M_V$', '$V-I$', '$W$']#, 'Y', 'Z']

x = np.arange(len(labels))
width = 0.45

fig, ax = plt.subplots(figsize=(6.64306*2,4.10318716*0.9))
plt.suptitle('RR Lyrae', y=0.99)
rects1 = ax.bar(x, (all_sds[1,:] - all_sds[0,:])/all_sds[1,:], width)

ax.set_ylabel('Accuracy improvement\n$(\sigma_{\\rm{LM}} - \sigma_{\\rm{ML}}) / \sigma_{\\rm{LM}}$')
ax.set_xticks(x)
ax.set_xticklabels(labels)

ax.set_ylim([0, 0.8])
ax.set_yticks([0, 0.2, 0.4, 0.6, 0.8])
locs, labels = plt.yticks()
ax.set_yticklabels([str(int(loc*100)) + '%' for loc in locs])

def autolabel(rects):
    for rect in rects:
        height = rect.get_height()
        ax.annotate('%.1f' % (height*100) + '%',
                    xy=(rect.get_x() + rect.get_width() / 2, height),
                    xytext=(0, 3), 
                    textcoords="offset points",
                    ha='center', va='bottom')

autolabel(rects1)

plt.savefig('RRL-sds-ML-LM.pdf')
plt.show()

ml_stds = []
ml_r2s = []

for idx, name in enumerate(X_columns):
    std = np.array(np.array(ml_pred[:,idx]) - DF[name]).std()
    print(name, 'std: %.3f' % std)
    ml_stds += [std]
    
    SS_tot = np.sum((DF[name] - DF[name].mean())**2)
    SS_res = np.sum((ml_pred[:,idx] - DF[name].values)**2)
    R2 = 1-SS_res/SS_tot
    ml_r2s += [R2]
    
    print(name, 'R2: %.3f' % R2)

lm_stds = []
lm_r2s = []

for idx, name in enumerate(X_columns):
    std = np.array(np.array(lm_pred[:,idx]) - DF[name]).std()
    print(name, 'std: %.3f' % std)
    lm_stds += [std]
    
    SS_tot = np.sum((DF[name] - DF[name].mean())**2)
    SS_res = np.sum((lm_pred[:,idx] - DF[name].values)**2)
    R2 = 1-SS_res/SS_tot
    lm_r2s += [R2]
    
    print(name, 'R2: %.3f' % R2)

var = 'M'
lims = [0.4, 0.9]

plt.figure(figsize=(6.64306,4.10318716*2))
plt.suptitle('Machine Learning: RR Lyrae', y=0.995)

ax = plt.subplot(2,1,1)
plt.scatter(DF[var], ml_pred[:,X_columns.index(var)], c=DF['logL'])
plt.plot([lims[0],lims[1]], [lims[0],lims[1]], 'k--', alpha=0.75, zorder=0)
plt.ylabel('Predicted Mass [M$_\odot$]')
plt.xlim(lims)
plt.ylim(lims)
plt.colorbar(label='$\log\;L/$L$_\odot$', pad=0.02)
plt.text(0.764, 0.42, '$R^2$: %.3f' % ml_r2s[X_columns.index(var)])

ax2 = plt.subplot(2,1,2, sharex=ax)
plt.axhline(0, c='k', ls='--', linewidth=1.25)
plt.scatter(DF[var], ml_pred[:,X_columns.index(var)] - DF[var], c=DF['logL'])
plt.xlabel('Model Mass [M$_\odot$]')
plt.ylabel('Errors [M$_\odot$]')
plt.colorbar(label='$\log\;L/$L$_\odot$', pad=0.02)
plt.xlim(lims)
plt.ylim([-0.3, 0.3])
plt.text(0.41, -0.28, 'St. dev: %.3f M$_\odot$' % ml_stds[X_columns.index(var)])

plt.savefig('RRL-ML-mass.pdf')#, bbox_inches='tight')
plt.show()

var = 'logR'
lims = [0.55, 1.15]

fig = plt.figure(figsize=(6.64306,4.10318716*2))
plt.suptitle('Linear Model: RR Lyrae', y=0.995)

ax = plt.subplot(2,1,1)
plt.scatter(DF[var], lm_pred[:,X_columns.index(var)], c=DF['M'])
plt.plot([lims[0],lims[1]], [lims[0],lims[1]], 'k--', alpha=0.75, zorder=0)
plt.ylabel('Predicted Radius [$\log$ R$_\odot$]')
plt.xlim(lims)
plt.ylim(lims)
plt.colorbar(label='Mass [M$_\odot$]', pad=0.02)
plt.text(0.98, 0.572, '$R^2$: %.3f' % lm_r2s[X_columns.index(var)])

ax2 = plt.subplot(2,1,2, sharex=ax)
plt.plot([lims[0], lims[1]], [0, 0], 'k--', zorder=0, linewidth=1.25)
plt.scatter(DF[var], lm_pred[:,X_columns.index(var)] - DF[var], c=DF['M'])
plt.colorbar(label='Mass [M$_\odot$]', pad=0.02)
plt.xlabel('Model Radius [$\log$ R$_\odot$]')
plt.ylabel('Errors [$\log$ R$_\odot$]')
plt.xlim(lims)
plt.ylim([-0.08, 0.08])
plt.text(0.56, 0.065, 'St. dev: %.3f $\log$ R$_\odot$' % lm_stds[X_columns.index(var)])

plt.savefig('RRL-LM-radius.pdf')#, bbox_inches='tight')
plt.show()

var = 'logR'
lims = [0.55, 1.15]

plt.figure(figsize=(6.64306,4.10318716*2))
plt.suptitle('Machine Learning: RR Lyrae', y=0.995)

ax = plt.subplot(2,1,1)
plt.scatter(DF[var], ml_pred[:,X_columns.index(var)], c=DF['M'])
plt.plot([lims[0],lims[1]], [lims[0],lims[1]], 'k--', alpha=0.75, zorder=0)
plt.ylabel('Predicted Radius [$\log$ R$_\odot$]')
plt.xlim(lims)
plt.ylim(lims)
plt.colorbar(label='Mass [M$_\odot$]', pad=0.02)
plt.text(0.98, 0.572, '$R^2$: %.3f' % ml_r2s[X_columns.index(var)])

ax2 = plt.subplot(2,1,2, sharex=ax)
plt.plot([lims[0], lims[1]], [0, 0], 'k--', zorder=0, linewidth=1.25)
plt.scatter(DF[var], ml_pred[:,X_columns.index(var)] - DF[var], c=DF['M'])
plt.colorbar(label='Mass [M$_\odot$]', pad=0.02)
plt.xlabel('Model Radius [$\log$ R$_\odot$]')
plt.ylabel('Errors [$\log$ R$_\odot$]')
plt.xlim(lims)
plt.ylim([-0.08, 0.08])
plt.text(0.56, 0.065, 'St. dev: %.3f $\log$ R$_\odot$' % ml_stds[X_columns.index(var)])

plt.savefig('RRL-ML-radius.pdf')#, bbox_inches='tight')
plt.show()

var = 'logL'
lims = [1.4, 2.2]

fig = plt.figure(figsize=(6.64306,4.10318716*2))
plt.suptitle('Linear Model: RR Lyrae', y=0.995)

ax = plt.subplot(2,1,1)
plt.scatter(DF[var], lm_pred[:,X_columns.index(var)], c=DF['Teff'])
plt.plot([lims[0],lims[1]], [lims[0],lims[1]], 'k--', alpha=0.75, zorder=0)
plt.ylabel('Predicted Luminosity [$\log$ L$_\odot$]')
plt.xlim(lims)
plt.ylim(lims)
plt.colorbar(label='$T_{\\rm{eff}}/$K', pad=0.02)
plt.text(1.978, 1.42, '$R^2$: %.3f' % lm_r2s[X_columns.index(var)])

ax2 = plt.subplot(2,1,2, sharex=ax)
plt.plot([lims[0], lims[1]], [0, 0], 'k--', zorder=0, linewidth=1.25)
plt.scatter(DF['logL'], lm_pred[:,X_columns.index(var)] - DF['logL'], c=DF['Teff'])
plt.colorbar(label='$T_{\\rm{eff}}/$K', pad=0.02)
plt.xlabel('Model Luminosity [$\log$ L$_\odot$]')
plt.ylabel('Errors [$\log$ L$_\odot$]')
plt.xlim(lims)
plt.ylim([-0.3, 0.3])
plt.text(1.42, -0.28, 'St. dev: %.3f $\log$ L$_\odot$' % lm_stds[X_columns.index(var)])

plt.savefig('RRL-LM-luminosity.pdf')#, bbox_inches='tight')
plt.show()

var = 'logL'
lims = [1.4, 2.2]

plt.figure(figsize=(6.64306,4.10318716*2))
plt.suptitle('Machine Learning: RR Lyrae', y=0.995)

ax = plt.subplot(2,1,1)
plt.scatter(DF[var], ml_pred[:,X_columns.index(var)], c=DF['Teff'])
plt.plot([lims[0],lims[1]], [lims[0],lims[1]], 'k--', alpha=0.75, zorder=0)
plt.ylabel('Predicted Luminosity [$\log$ L$_\odot$]')
plt.xlim(lims)
plt.ylim(lims)
plt.colorbar(label='$T_{\\rm{eff}}/$K', pad=0.02)
plt.text(1.978, 1.42, '$R^2$: %.3f' % ml_r2s[X_columns.index(var)])

ax2 = plt.subplot(2,1,2, sharex=ax)
plt.plot([lims[0], lims[1]], [0, 0], 'k--', zorder=0, linewidth=1.25)
plt.scatter(DF['logL'], ml_pred[:,X_columns.index(var)] - DF['logL'], c=DF['Teff'])
plt.colorbar(label='$T_{\\rm{eff}}/$K', pad=0.02)
plt.xlabel('Model Luminosity [$\log$ L$_\odot$]')
plt.ylabel('Errors [$\log$ L$_\odot$]')
plt.xlim(lims)
plt.ylim([-0.3, 0.3])
plt.text(1.42, -0.28, 'St. dev: %.3f $\log$ L$_\odot$' % ml_stds[X_columns.index(var)])

plt.savefig('RRL-ML-luminosity.pdf')#, bbox_inches='tight')
plt.show()

var = 'W'
lims = [-2.5, 0]

fig = plt.figure(figsize=(6.64306,4.10318716*2))
plt.suptitle('Linear Model: RR Lyrae', y=0.995)

ax = plt.subplot(2,1,1)
plt.scatter(DF[var], lm_pred[:,X_columns.index(var)], c=DF['M'])
plt.plot([lims[0],lims[1]], [lims[0],lims[1]], 'k--', alpha=0.75, zorder=0)
plt.ylabel('Predicted Wesenheit Index')
plt.xlim(lims)
plt.ylim(lims)
plt.colorbar(label='$M/$M$_\odot$', pad=0.02)
plt.text(-0.7, -2.42, '$R^2$: %.3f' % lm_r2s[X_columns.index(var)])

ax2 = plt.subplot(2,1,2, sharex=ax)
plt.plot([lims[0], lims[1]], [0, 0], 'k--', zorder=0, linewidth=1.25)
plt.scatter(DF[var], lm_pred[:,X_columns.index(var)] - DF[var], c=DF['M'])
plt.colorbar(label='$M/$M$_\odot$', pad=0.02)
plt.xlabel('Model Wesenheit Index')
plt.ylabel('Errors [mag]')
plt.xlim(lims)
plt.ylim([-0.4, 0.4])
plt.text(-2.44, -0.37, 'St. dev: %.2f mag' % lm_stds[X_columns.index(var)])

plt.savefig('RRL-LM-W.pdf')#, bbox_inches='tight')
plt.show()

var = 'W'
lims = [-2.5, 0]

plt.figure(figsize=(6.64306,4.10318716*2))
plt.suptitle('Machine Learning: RR Lyrae', y=0.995)

ax = plt.subplot(2,1,1)
plt.scatter(DF[var], ml_pred[:,X_columns.index(var)], c=DF['M'])
plt.plot([lims[0],lims[1]], [lims[0],lims[1]], 'k--', alpha=0.75, zorder=0)
plt.ylabel('Predicted Wesenheit Index')
plt.xlim(lims)
plt.ylim(lims)
plt.colorbar(label='$M/$M$_\odot$', pad=0.02)
plt.text(-0.7, -2.42, '$R^2$: %.3f' % ml_r2s[X_columns.index(var)])

ax2 = plt.subplot(2,1,2, sharex=ax)
plt.plot([lims[0], lims[1]], [0, 0], 'k--', zorder=0, linewidth=1.25)
plt.scatter(DF[var], ml_pred[:,X_columns.index(var)] - DF[var], c=DF['M'])
plt.colorbar(label='$M/$M$_\odot$', pad=0.02)
plt.xlabel('Model Wesenheit Index')
plt.ylabel('Errors [mag]')
plt.xlim(lims)
plt.ylim([-0.4, 0.4])
plt.text(-2.44, -0.37, 'St. dev: %.3f mag' % ml_stds[X_columns.index(var)])

plt.savefig('RRL-ML-W.pdf')#, bbox_inches='tight')
plt.show()

# (neural network with log P and lightcurve) vs (linear model with log P) with the 5x2 cv test 

np.random.seed(0) 

ml = TransformedTargetRegressor(regressor=Pipeline(
            [('scaler', RobustScaler()), 
             ('mlp', MLPRegressor(hidden_layer_sizes=(100, 100), 
                                  max_iter=2000, activation='relu', solver='lbfgs', tol=1e-7))]),
         transformer=RobustScaler())
lm = LinearRegression()

first_p1 = 0 
s2s = []

for ii in range(5):
    ml_scores = []
    lm_scores = []
    
    cv = ShuffleSplit(n_splits=2, random_state=ii, train_size=0.5, test_size=0.5)
    
    for train_index, test_index in cv.split(DF):
        ml.fit(DF.iloc[train_index][y_columns], DF.iloc[train_index][X_columns])
        lm.fit(DF.iloc[train_index][['logP']],  DF.iloc[train_index][X_columns])
        
        ml_predict = ml.predict(DF.iloc[test_index][y_columns])
        lm_predict = lm.predict(DF.iloc[test_index][['logP']])
        
        ml_scores_ = []
        lm_scores_ = []
        
        for idx, name in enumerate(X_columns):
            SS_tot = np.sum((DF.iloc[test_index][name] - DF.iloc[test_index][name].mean())**2)
            
            ml_SS_res = np.sum((ml_predict[:,idx] - DF.iloc[test_index][name].values)**2)
            ml_R2 = 1-ml_SS_res/SS_tot
            #ml_scores_ += [ml_R2]
            ml_scores_ += [ml_SS_res]
            
            lm_SS_res = np.sum((lm_predict[:,idx] - DF.iloc[test_index][name].values)**2)
            lm_R2 = 1-lm_SS_res/SS_tot
            #lm_scores_ += [lm_R2]
            lm_scores_ += [lm_SS_res]
        
        ml_scores += [ml_scores_]
        lm_scores += [lm_scores_]
    
    ml_scores = np.array(ml_scores)
    lm_scores = np.array(lm_scores)
    
    p1 = ml_scores[0,:] - lm_scores[0,:]
    p2 = ml_scores[1,:] - lm_scores[1,:]
    p_bar = (p1 + p2)/2
    s2 = (p1 - p_bar)**2 + (p2 - p_bar)**2
    
    if ii == 0:
        first_p1 = p1
    s2s += [s2]

first_p1 = np.array(first_p1)
s2s = np.array(s2s)

for idx, name in enumerate(X_columns):
    p = stats.t.sf(np.abs(first_p1[idx]) / np.sqrt(1/5. * sum(s2s[:,idx])), 5)/2. # one sided test 
    print(name, '%.2E' % p) 

# (neural network with log P and lightcurve) vs (linear model with log P) with Wilcoxon signed-rank test 

np.random.seed(0) 

ml = TransformedTargetRegressor(regressor=Pipeline(
            [('scaler', RobustScaler()), 
             ('mlp', MLPRegressor(hidden_layer_sizes=(100, 100), 
                                  max_iter=2000, activation='relu', solver='lbfgs', tol=1e-7))]),
         transformer=RobustScaler())
lm = LinearRegression()

cv = ShuffleSplit(n_splits=2, random_state=0, train_size=0.5, test_size=0.5)
train_set, test_set = list(cv.split(DF))[0]

ml.fit(DF.iloc[train_set][y_columns], DF.iloc[train_set][X_columns])
lm.fit(DF.iloc[train_set][['logP']],  DF.iloc[train_set][X_columns])

ml_predict = ml.predict(DF.iloc[test_set][y_columns])
lm_predict = lm.predict(DF.iloc[test_set][['logP']])

for idx, name in enumerate(X_columns):
    ml_SE = (ml_predict[:,idx] - DF.iloc[test_set][name].values)**2
    lm_SE = (lm_predict[:,idx] - DF.iloc[test_set][name].values)**2
    
    wil = sp.stats.wilcoxon(ml_SE, lm_SE)
    print(name, ml_SE.mean() < lm_SE.mean(), wil.pvalue/2.) # one sided test 

# (neural network with log P and lightcurve) vs (linear model with log P and lightcurve)

np.random.seed(0) 

ml = TransformedTargetRegressor(regressor=Pipeline(
            [('scaler', RobustScaler()), 
             ('mlp', MLPRegressor(hidden_layer_sizes=(100, 100), 
                                  max_iter=2000, activation='relu', solver='lbfgs', tol=1e-7))]),
         transformer=RobustScaler())
lm = LinearRegression()

first_p1 = 0 
s2s = []

for ii in range(5):
    ml_scores = []
    lm_scores = []
    
    cv = ShuffleSplit(n_splits=2, random_state=ii, train_size=0.5, test_size=0.5)
    
    for train_index, test_index in cv.split(DF):
        ml.fit(DF.iloc[train_index][y_columns], DF.iloc[train_index][X_columns])
        lm.fit(DF.iloc[train_index][y_columns], DF.iloc[train_index][X_columns])
        
        ml_predict = ml.predict(DF.iloc[test_index][y_columns])
        lm_predict = lm.predict(DF.iloc[test_index][y_columns])
        
        ml_scores_ = []
        lm_scores_ = []
        
        for idx, name in enumerate(X_columns):
            SS_tot = np.sum((DF.iloc[test_index][name] - DF.iloc[test_index][name].mean())**2)
            
            ml_SS_res = np.sum((ml_predict[:,idx] - DF.iloc[test_index][name].values)**2)
            ml_R2 = 1-ml_SS_res/SS_tot
            #ml_scores_ += [ml_R2]
            ml_scores_ += [ml_SS_res]
            
            lm_SS_res = np.sum((lm_predict[:,idx] - DF.iloc[test_index][name].values)**2)
            lm_R2 = 1-lm_SS_res/SS_tot
            #lm_scores_ += [lm_R2]
            lm_scores_ += [lm_SS_res]
        
        ml_scores += [ml_scores_]
        lm_scores += [lm_scores_]
    
    ml_scores = np.array(ml_scores)
    lm_scores = np.array(lm_scores)
    
    p1 = ml_scores[0,:] - lm_scores[0,:]
    p2 = ml_scores[1,:] - lm_scores[1,:]
    p_bar = (p1 + p2)/2
    s2 = (p1 - p_bar)**2 + (p2 - p_bar)**2
    
    if ii == 0:
        first_p1 = p1
    s2s += [s2]

first_p1 = np.array(first_p1)
s2s = np.array(s2s)

for idx, name in enumerate(X_columns):
    p = stats.t.sf(np.abs(first_p1[idx]) / np.sqrt(1/5. * sum(s2s[:,idx])), 5)/2. # one sided test 
    print(name, '%.2E' % p) 

np.random.seed(0) 

ml = TransformedTargetRegressor(regressor=Pipeline(
            [('scaler', RobustScaler()), 
             ('mlp', MLPRegressor(hidden_layer_sizes=(100, 100), 
                                  max_iter=2000, activation='relu', solver='lbfgs', tol=1e-7))]),
         transformer=RobustScaler())
lm = LinearRegression()

ml.fit(DF[y_columns], DF[X_columns])
lm.fit(DF[['logP']], DF[X_columns])

import joblib
joblib.dump(ml, 'RR_Lyrae.ann')

maxs = DF[y_columns].quantile(0.95)
mins = DF[y_columns].quantile(0.05)
print(maxs, '\n\n', mins)