<jupyter_start><jupyter_text>Imports<jupyter_code>import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import xarray as xr
import scipy as sp
import os
from sklearn.decomposition import PCA
from sklearn.linear_model import LinearRegression
from sklearn.linear_model import QuantileRegressor
from sklearn.metrics import mean_absolute_error
from sklearn.metrics import mean_squared_error
import glob
hr_eigs = np.array( pd.read_csv('/ocean/arandhawa/MOAD/analysis-armaan/notebooks/eigvecs_hrdps_annual.csv', header = None) )
hr_PCs = np.array( pd.read_csv('/ocean/arandhawa/MOAD/analysis-armaan/notebooks/PCs_hrdps_annual.csv', header = None) )
print(hr_PCs.shape)
print(hr_eigs.shape)
psl_eigs = np.array( pd.read_csv('/ocean/arandhawa/MOAD/analysis-armaan/notebooks/eigvecs_psl.csv', header = None) )
psl_PCs = np.array( pd.read_csv('/ocean/arandhawa/MOAD/analysis-armaan/notebooks/PCs_psl.csv', header = None) )
print(psl_PCs.shape)
print(psl_eigs.shape)
cw_eigs = np.array( pd.read_csv('/ocean/arandhawa/MOAD/analysis-armaan/notebooks/eigvecs_winds.csv', header = None) )
cw_PCs = np.array( pd.read_csv('/ocean/arandhawa/MOAD/analysis-armaan/notebooks/PCs_winds.csv', header = None) )
print(cw_PCs.shape)
print(cw_eigs.shape)
##data must be converted into a 2D matrix for pca analysis
##transform takes a 3D data array (time, a, b) -> (a*b, time)
##(the data grid is flattened a column using numpy.flatten)

def transform(xarr):
    arr = np.array(xarr) ##converting to numpy array
    arr = arr.reshape(arr.shape[0], arr.shape[1]*arr.shape[2]) ##reshaping from size (a, b, c) to (a, b*c)
    arr = arr.transpose()
    return arr

def reverse(mat, orig_shape):
    arr = np.copy(mat)
    arr = arr.transpose()
    arr = arr.reshape(orig_shape[0], orig_shape[1], orig_shape[2]) ##reshaping back to original array shape
    return arr
comb_eigs = np.empty((50, 2904))
for i in range(25):
    comb_eigs[2*i] = cw_eigs[i]
    comb_eigs[2*i + 1] = psl_eigs[i]
def get_scores(num_scores, eigvecs, arr, row):  
    scores = np.empty(num_scores)
    for i in range(num_scores):
        vec = eigvecs[0 : i+1 , :].T
        reg = QuantileRegressor(quantile = 0.5, alpha = 0, solver = 'highs').fit(vec, arr[row,:])
        scores[i] = reg.score(vec, arr[row,:])
    return scores
n = 5
fig, ax = plt.subplots(1, n, figsize=(5*n, 5))

for i in range(n):
    ax[i].plot(get_scores(30, cw_eigs, hr_eigs, i), label = "winds")
    ax[i].plot(get_scores(30, psl_eigs, hr_eigs, i), label = "pressure")
    ax[i].plot(get_scores(30, comb_eigs, hr_eigs, i), label = "combined")
    ax[i].set_title("Fitting CANRCM to Mode {} of HRDPS".format(i))
    ax[i].set_xlabel("CANRCM modes used")
    ax[i].set_ylabel("Correlation Coefficient ($R^2$)")
    ax[i].legend()
    
plt.tight_layout()
def get_scores2(num_scores, eigvecs, arr, row):  
    scores = np.empty(num_scores)
    for i in range(num_scores):
        vec = eigvecs[0 : i+1 , :].T
        reg = QuantileRegressor(quantile = 0.5, alpha = 0, solver = 'highs').fit(vec, arr[row,:])
        scores[i] = mean_absolute_error(
    return scores
n = 5
fig, ax = plt.subplots(1, n, figsize=(5*n, 5))

for i in range(n):
    ax[i].plot(get_scores2(30, cw_eigs, hr_eigs, i), label = "winds")
    ax[i].plot(get_scores2(30, psl_eigs, hr_eigs, i), label = "pressure")
    ax[i].plot(get_scores2(30, comb_eigs, hr_eigs, i), label = "combined")
    ax[i].set_title("Fitting CANRCM to Mode {} of HRDPS".format(i))
    ax[i].set_xlabel("CANRCM modes used")
    ax[i].set_ylabel("Correlation Coefficient ($R^2$)")
    ax[i].legend()
    
plt.tight_layout()
##testing quantile regressor - regresses by minimizing absolute error rather than squared error
X = psl_eigs[0:30,:].T
y = hr_eigs[1]
print(X.shape, y.shape)

reg = QuantileRegressor(quantile = 0.5, alpha = 0, solver = 'highs').fit(X, y)
result = reg.predict(X)

r_squared = reg.score(X, y)
MAE = mean_absolute_error(y, result)
print(MAE)
plt.plot(hr_eigs[0])
plt.plot(result)
def fit_modes(vectors, num_vec, data, result_size):  
    result = np.empty((result_size, vectors.shape[1]))
    scores = np.empty(result_size)
    X = vectors[0:num_vec,:].T
        
    for i in range(result_size):
        y = data[i,:]
        reg = QuantileRegressor(quantile = 0.5, alpha = 0, solver = 'highs').fit(X, y)
        result[i] = reg.predict(X)
        scores[i] = reg.score(X, y)

    return (result, scores)
n = 10
new_eigs, scores = fit_modes(psl_eigs, n, hr_eigs, 5)

fig, ax = plt.subplots(1, 5, figsize=(20, 5))

for i in range(5):
    
    MAE = mean_absolute_error(hr_eigs[i], new_eigs[i])
    rmse = np.sqrt( mean_squared_error(hr_eigs[i], new_eigs[i]) )
    r_squared = scores[i]
    
    ax[i].plot(hr_eigs[i], label = "HRDPS")
    ax[i].plot(new_eigs[i], label = "{} pressure CANRCM modes".format(n))
    ax[i].set_ylim(-0.1, 0.1)
    ax[i].annotate("r-squared = {:.3f}".format(r_squared), (0, -0.07))
    ax[i].annotate("MAE = {:.3f}".format(MAE), (0, -0.08))
    ax[i].annotate("RMSE = {:.3f}".format(rmse), (0, -0.09))
    ax[i].set_title("Mode {}".format(i))
    ax[i].legend()
    
    print(np.sqrt( np.mean(hr_eigs[i]**2) / np.mean(new_eigs[i]**2) ))
plt.tight_layout()
n = 10
new_eigs2, scores = fit_modes(cw_eigs, n, hr_eigs, 5)

fig, ax = plt.subplots(1, 5, figsize=(20, 5))

for i in range(5):
    
    MAE = mean_absolute_error(hr_eigs[i], new_eigs[i])
    rmse = np.sqrt( mean_squared_error(hr_eigs[i], new_eigs[i]) )
    r_squared = scores[i]
    
    ax[i].plot(hr_eigs[i], label = "HRDPS")
    ax[i].plot(new_eigs[i], label = "{} wind CANRCM modes".format(n))
    ax[i].set_ylim(-0.1, 0.1)
    ax[i].annotate("r-squared = {:.3f}".format(r_squared), (0, -0.07))
    ax[i].annotate("MAE = {:.3f}".format(MAE), (0, -0.08))
    ax[i].annotate("RMSE = {:.3f}".format(rmse), (0, -0.09))
    ax[i].set_title("Mode {}".format(i))
    ax[i].legend()
    
    print(np.sqrt( np.mean(hr_eigs[i]**2) / np.mean(new_eigs[i]**2) ))
plt.tight_layout()
n = 10
new_eigs, scores = fit_modes(comb_eigs, 3, hr_eigs, 30)

fig, ax = plt.subplots(1, 5, figsize=(20, 5))

for i in range(5):
    
    MAE = mean_absolute_error(hr_eigs[i], new_eigs[i])
    rmse = np.sqrt( mean_squared_error(hr_eigs[i], new_eigs[i]) )
    r_squared = scores[i]
    
    ax[i].plot(hr_eigs[i], label = "HRDPS")
    ax[i].plot(new_eigs[i], label = "{} combined CANRCM modes".format(n))
    ax[i].set_ylim(-0.1, 0.1)
    ax[i].annotate("r-squared = {:.3f}".format(r_squared), (0, -0.07))
    ax[i].annotate("MAE = {:.3f}".format(MAE), (0, -0.08))
    ax[i].annotate("RMSE = {:.3f}".format(rmse), (0, -0.09))
    ax[i].set_title("Mode {}".format(i))
    ax[i].legend()
    
    print(np.sqrt( np.mean(hr_eigs[i]**2) / np.mean(new_eigs[i]**2) ))
plt.tight_layout()<jupyter_output>/home/arandhawa/conda_envs/analysis-armaan/lib/python3.10/site-packages/sklearn/linear_model/_quantile.py:270: ConvergenceWarning: Linear programming for QuantileRegressor did not succeed.
Status is 4: Numerical difficulties encountered.
Result message of linprog:
An optimal solution to the scaled model was found but was HiGHS Status Code 0: HighsModelStatusNOTSET in the unscaled model. For more information run with the option `disp: True`.
  warnings.warn(<jupyter_text>Conserving Energy<jupyter_code>new_eigs, scores = fit_modes(cw_eigs, 10, hr_eigs, 5)

fig, ax = plt.subplots(1, 5, figsize=(20, 5))

for i in range(5):
    
    f = np.sqrt( np.mean(hr_eigs[i]**2) / np.mean(new_eigs[i]**2) )
    new_eigs[i] = new_eigs[i]*f
    
    MAE = mean_absolute_error(hr_eigs[i], new_eigs[i])
    rmse = np.sqrt( mean_squared_error(hr_eigs[i], new_eigs[i]) )
    ##correlation_matrix = np.corrcoef(x, y)
    ##r_squared = correlation_matrix[0,1]**2
    ##r_squared = scores[i]
    ax[i].plot(hr_eigs[i], label = "HRDPS")
    ax[i].plot(new_eigs[i], label = "10 wind CANRCM modes")
    ax[i].set_ylim(-0.1, 0.1)
    ##ax[i].annotate("r-squared = {:.3f}".format(r_squared), (0, -0.07))
    ax[i].annotate("MAE = {:.3f}".format(MAE), (0, -0.08))
    ax[i].annotate("RMSE = {:.3f}".format(rmse), (0, -0.09))
    ax[i].set_title("Mode {}".format(i))
    ax[i].legend()
    
plt.tight_layout()
def equalize_energy(eigs, old_eigs):
    new_eigs = np.sqrt( (old_eigs[0:eigs.shape[0]]**2).mean(axis = 1)/(eigs**2).mean(axis = 1))* eigs.T
    return new_eigs.T<jupyter_output><empty_output><jupyter_text>Exports<jupyter_code>new_eigs, _ = fit_modes(cw_eigs, 30, hr_eigs, 30)
new_eigs = equalize_energy(new_eigs, hr_eigs)
new_eigs2.shape
np.savetxt("eigvecs_recon.csv", new_eigs2, delimiter=",")
##MAE (duplicate notebook)
##test different #'s of modes

##compare total energy, for each new mode when reconstructed from 10, 20, and 30 modes
##do the dot plots from last notebook<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17546867855509832, 'fraction_numerical': 0.025377229080932786, 'mean_word_length': 4.192284866468842, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '404321', 'license': 'Apache-2.0', 'n_tokens_mistral': 3889, 'n_tokens_neox': 3656, 'n_words': 711, 'path': 'notebooks/reconstructingdataP3.ipynb', 'repo_name': 'UBC-MOAD/analysis-armaan', 'size': '8748'}