# Last edited on 2021-07-31 15:43:01 by jstolfi

import rn
from math import sqrt

def basis_sampled_orthize(B, w):
  # Makes the elements of the basis {B} orhtonormal with the Gram-Schmidt method.
  #
  # The parameter {B} must be a list of list of floats.
  # Element {B[i][k]} is assumed to be the basis element {i} evaluated 
  # on the sampling point {k}. 
  #
  # If {w} is not {None}, {w[k]} is the weight of sampling point{k}.
  # If {w} is {None}, the weights are assumed to be all 1
  # The result will be orthonormal relative to the inner product 
  # {<x|y>_w = \sum x[k]*y[k]*w[k]}.
  
  C = []
  for Bi in B:
    # Orthize {Bi} rel to {C[0..i-1]}:
    Ci = Bi.copy()
    for Cj in C:
      c = rn.wdot(Ci, Cj, w)
      Ci = rn.mix(1, Ci, -c, Cj);
    d = sqrt(rn.wdot(Ci, Ci, w))
    assert abs(d) > 1.0e-12, "{B} is not linearly independent"
    Ci = rn.scale(1/d, Ci)
    C.append(Ci)
  return C
  # ----------------------------------------------------------------------