Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

SVD for orth. projection vectors #17

Open
wants to merge 3 commits into
base: main
Choose a base branch
from
Open

SVD for orth. projection vectors #17

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
cee67b6
svd implementation of distrotion matrix and projection vectors
p-slash Oct 2, 2023
cbf2eef
Merge branch 'main' into distortion-vectors-wsvd
p-slash Oct 2, 2023
b4fc6fd
fix zero weight in get_orthonormal_vectors_svd
p-slash Oct 2, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
17 changes: 12 additions & 5 deletions lya_2pt/tracer.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,8 @@
import numpy as np

from lya_2pt.constants import ABSORBER_IGM
from lya_2pt.tracer_utils import rebin, project_deltas, get_angle_list, gram_schmidt
from lya_2pt.tracer_utils import (
rebin, project_deltas, get_angle_list, gram_schmidt, get_orthonormal_vectors_svd)


class Tracer:
Expand Down Expand Up @@ -269,18 +270,24 @@ def rebin(self, rebin_factor, dwave, absorption_line):
self.logwave_term = log_lambda - np.sum(log_lambda * weights) / self.sum_weights
self.term3_norm = (weights * self.logwave_term**2).sum()

def project(self, old_projection=True):
def project(self, old_projection=True, use_svd=True):
"""Apply projection matrix to deltas"""
assert not self.is_projected, "Tracer already projected"

if old_projection:
self.deltas = project_deltas(self.log_lambda, self.deltas, self.weights, self.order)
self.is_projected = True
return

if use_svd:
basis = get_orthonormal_vectors_svd(self.log_lambda, self.weights, self.order)
else:
basis = gram_schmidt(self.log_lambda, self.weights, self.order)

for b in basis:
self.deltas -= b * np.dot(b * self.weights, self.deltas)
proj_vec_mat = basis * self.weights
self.deltas -= basis.T.dot(proj_vec_mat.dot(self.deltas))

self.proj_vec_mat = (basis * self.weights).T
self.proj_vec_mat = proj_vec_mat.T

self.is_projected = True

Expand Down
10 changes: 10 additions & 0 deletions lya_2pt/tracer_utils.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -189,6 +189,16 @@ def get_projection_matrix(log_lambda, weights, order):
return Vh.T, np.eye(weights.size) - Vh.T @ Vh


def get_orthonormal_vectors_svd(log_lambda, weights, order):
wsqrt = np.sqrt(weights)
input_vectors_matrix = np.vander(log_lambda, order + 1).T * wsqrt
Vh = np.linalg.svd(input_vectors_matrix, full_matrices=False)[2]
s = weights != 0
Vh[~s] = 0
Vh[s] /= wsqrt[s]
return Vh


def gram_schmidt(log_lambda, weights, order):
basis = []
for n in range(order + 1):
Expand Down