Convert_NP = lambda g: g.to_numpy().reshape(-1,1).astype('float32') #convert to one column ndarray
X = Convert_NP(x) #shape (560,1)
Y = Convert_NP(y) #shape (560,1)

X_NoBias = X #USED FOR PLOTTING PURPOSES


bias = np.ones((X.shape[0],1)).astype('float32') #column of 1's represent intercept coefficient
X_bias = np.append(X,bias,axis=1) #shape (560,2)
X = X_bias
print(f"X-dim: {X.shape}\nY-dim: {Y.shape}")

x_manual = np.linspace(Y.min()-Y.std(),Y.max()+Y.std(), Y.shape[0])
plt.scatter(X_NoBias,Y)

Resid_normed = Residual/Residual.std()

def Filliben(val):
    n = val.shape[0]
    fst = 1-(0.5**(1/n))
    midd = [(i-0.3175)/(n+0.365) for i in range(1,n-1)]
    lst = 0.5**(1/n)
    return np.append(np.append(fst,midd),lst)

Resid_normed.sort(axis=0)

print(np.append(Resid_normed[0], Resid_normed[1]))
InvertedPercentile = [stats.norm.ppf(i) for i in Filliben(Resid_normed)]
plt.scatter(InvertedPercentile,Resid_normed)

Residual_Standardized = Residual/Residual.std()
# notice the standardized residual is the z-score

normdist = np.random.normal( size=506)
print(Residual_Standardized.mean())
print(Residual_Standardized.std())

fig,ax = plt.subplots(2,2)
ax[0,0].scatter(Y_hat,Residual,facecolors="none",edgecolors="black")
ax[0,0].set_title("yhat vs e")
#yy = [stats.percentileofscore(Residual_Standardized, a, 'rank') for a in Residual_Standardized]
#yy2=np.sort(yy)
SortedResidStand = np.sort(Residual_Standardized,axis=0)
normdist = np.random.normal( size=506)
Sortednormdist = np.sort(normdist)
ax[0,1].scatter(Sortednormdist,SortedResidStand)
ax[1,0].scatter(Y_hat,np.sqrt(np.abs(Residual_Standardized)))

plt.show()

Var_Yhat = ((Y_hat - Y_hat.mean())**2).sum()/Y_hat.shape[0]
Std_Yhat = Var_Yhat**0.5
TSS = ((Y_hat - Y_hat.mean())**2).sum()

assert np.isclose(Var_Yhat, Y_hat.var())
assert np.isclose(Std_Yhat, Y_hat.std())

SE_Yhat = (Var_Yhat)/Y_hat.shape[0]
print(SE_Yhat)

import torch
Xtensor = torch.from_numpy(X)
Ytensor = torch.from_numpy(Y)

W = torch.linalg.lstsq(Xtensor, Ytensor).solution
print(W)
# print(X[:,[1]])
Xlayer = torch.from_numpy(X[:,[0]])
hidlayer = torch.from_numpy(X[:,[1]])
Ytensor = torch.from_numpy(Y)

import random
import torch
from torch import nn, optim
import math
from IPython import display

device = torch.device("cuda:0" if torch.cuda.is_available() else "cpu")
seed = 1
random.seed(seed)
torch.manual_seed(seed)

inputDim = 1
outputDim = 1
hiddenDim = 100

learning_rate = 1e-4
lambda_l2 = 1e-6

model = nn.Sequential(
    nn.Linear(inputDim, hiddenDim),
    nn.Linear(hiddenDim, outputDim)
)
model.to(device) # Convert to CUDA

# nn package also has different loss functions.
# we use MSE loss for our regression task
criterion = torch.nn.MSELoss()

# we use the optim package to apply
# stochastic gradient descent for our parameter updates
optimizer = torch.optim.SGD(model.parameters(), lr=learning_rate, weight_decay=lambda_l2) # built-in L2

# Training
for t in range(1000):
    
    # Feed forward to get the logits
    y_pred = model(Xlayer)
    
    # Compute the loss (MSE)
    loss = criterion(y_pred, Ytensor)
    print("[EPOCH]: %i, [LOSS or MSE]: %.6f" % (t, loss.item()))
    display.clear_output(wait=True)
    
    # zero the gradients before running
    # the backward pass.
    optimizer.zero_grad()
    
    # Backward pass to compute the gradient
    # of loss w.r.t our learnable params. 
    loss.backward()
    
    # Update params
    optimizer.step()

plt.scatter(Xlayer.data.cpu().numpy(), Ytensor.data.cpu().numpy())
plt.plot(Xlayer.data.cpu().numpy(), y_pred.data.cpu().numpy(), 'r-', lw=5)
plt.axis('equal');

import torch
from torch.autograd import Variable
import torch.nn.functional as F

Xm = np.delete(X,1,1)
Xtensor = torch.from_numpy(Xm)
Ytensor = torch.from_numpy(Y)
numHiddenlayer = 1
#In pytorch we dont need to artifically add a hidden layer column to our X input
class LinearRegressionModel(torch.nn.Module):
 
    def __init__(self):
        super(LinearRegressionModel, self).__init__()
        self.linearA = torch.nn.Linear(Xtensor.shape[1], numHiddenlayer,bias=True)  # One in and one out
        self.linearB = torch.nn.Linear(numHiddenlayer,Ytensor.shape[1])
    def forward(self, x):
        #y_pred = F.relu(self.linear(x))
        y_pred = self.linearB(self.linearA(x))
        #y_pred = self.linearA(x)
        return y_pred
 
# our model
our_model = LinearRegressionModel()
 
criterion = torch.nn.MSELoss()

#EXTRREMELY SENSITIVE TO LEARNING RATE 1e-4 is sweetspot
learning_rate = (1e-4)

optimizer = torch.optim.SGD(our_model.parameters(), lr=learning_rate)

#w,b = our_model.parameters()
#print(f"initial: {w,b}")
bestloss = 99999
bestepoch = -1
for epoch in range(4000): #EXTREMELY SENSITIVE TO EPOCH, 50000 is sweetspot
    epoch += 1
    # Convert numpy array to torch Variable
    #inputs = torch.from_numpy(p).requires_grad_()
    #labels = torch.from_numpy(q)

    # Clear gradients w.r.t. parameters
    optimizer.zero_grad() 

    # Forward to get output
    outputs = our_model(Xtensor)

    # Calculate Loss
    loss = criterion(outputs, Ytensor)
    currloss = loss.item()
    if currloss < bestloss:
        bestloss = currloss
        bestepoch = epoch
    # if currloss - bestloss > 200 and epoch > 50:
        # print("boom")
    # Getting gradients w.r.t. parameters
    loss.backward()

    # Updating parameters
    optimizer.step()

print('epoch {}, loss {}'.format(epoch, loss.item()))
print('epoch {}, loss {}'.format(bestepoch, bestloss))
 
#w,b = our_model.parameters()
#print(w,b)

w,b = our_model.parameters()

print(w,b)