def get_pca(
    X : np.ndarray,
    alpha = None,
    r = None
) -> tuple:
  '''
  Performs principal component analysis for a data matrix X

  Returns

  P : np.ndarray
    projection matrix

  mean : np.ndarray
    mean vector
  '''

  # solve for eigenvectors/values
  mean = np.mean(X, axis=0,keepdims=True)
  z = X - mean
  w, V = eigh(np.cov(z.T))

  # Rearrange eigenvectors according to eigenvalues in descending order
  indices = np.argsort(w)[::-1]
  w = w[indices]
  V = V[:,indices]

  # if r value(s) are not provided, decide based on alpha
  if r is None:
    # fraction of variance for every component
    variance_explained = [w[i] / w.sum() for i in range(len(w))]
    cumulative_variance_explained = np.cumsum(variance_explained)

    # if one alpha is provided, decide the dimensions accordingly
    if type(alpha) == int:
      info = cumulative_variance_explained >= alpha
      for i,v in enumerate(info):
        if v:
          r = i + 1
          break
    # if a list of alpha values are provied, calculate r for each value
    elif alpha is not None:
      r = []
      for a in alpha:
        info = cumulative_variance_explained >= a
        for i,v in enumerate(info):
          if v:
            r_a = i + 1
            break
        r.append(r_a)

  print(f"at alpha = {alpha}, r = {r}")

  # if only one r value is decided, calculate the corresponding projection matrix
  if type(r) == int:
    P = V[:,:r]
  # if multiple r values, Calculate projection matrix for each r value
  else:
    P = []
    for i in r:
      P.append(V[:,:i])

  return (P, mean)

def LDA():
  numberOfClasses = 40
  numberInEachClass = 5
  index = 0
  Zindex = 0
  sumOfEach = np.zeros((40,10304))
  totalSum = np.zeros((1,10304))
  meanOfEach = np.zeros((40,10304))
  totalMean = np.zeros((1,10304))
  Z = np.zeros((200,10304))
  Sb = np.zeros((10304,10304))
  S = np.zeros((10304,10304))
  U = np.zeros((39,10304))
  w = np.zeros((1,39))

  for i in range(1,numberOfClasses+1):
    # gets the indexes of the label with id=i
    place = np.where(label_train[:] == i)
    #print(place[0])
    # sums the 5 pics values in the same id
    for j in range(numberInEachClass):
      #print(data_train[place[0][j],:])
      sumOfEach[index][:] = sumOfEach[index][:] + data_train[place[0][j],:]
      #print("sum here " + str(sumOfEach))
    # getting the mean of each class
    meanOfEach[index][:] = sumOfEach[index][:] / numberInEachClass
    #print("mean here "+ str(meanOfEach[index][:]))

    # getting the centre class matrices Z
    for k in range(numberInEachClass):
      Z[place[0][k],:] = data_train[place[0][k],:] - meanOfEach[index][:]
    #print("printing Z " + str(Z))
    index += 1

  totalSum = np.sum(sumOfEach, axis=0)
  #print(totalSum.shape)
  # getting the overall sample mean
  totalMean = totalSum / numberOfClasses
  #print(totalMean)

  # getting the between class scatter matrix
  for i in range(numberOfClasses):
    difference = meanOfEach[i] - totalMean
    Sb += (numberInEachClass * np.dot(difference, np.transpose(difference)))

  # getting the within class scatter matrices S
  #for i in range(numberOfClasses):
  S = np.dot(np.transpose(Z),Z)

  w, U = np.linalg.eigh(np.dot(np.linalg.inv(S),Sb))
  sorted_vectors=U[w.argsort()]
  sorted_vectors = np.real(sorted_vectors)
  U=sorted_vectors[-39:]
  print(U.shape)
  return U