Skip to content

Commit

Permalink
indentation corrected
Browse files Browse the repository at this point in the history
  • Loading branch information
abhi0395 committed Aug 29, 2023
1 parent 5ff5d17 commit 55b04ef
Showing 1 changed file with 67 additions and 67 deletions.
134 changes: 67 additions & 67 deletions py/redrock/zscan.py
Original file line number Diff line number Diff line change
Expand Up @@ -273,88 +273,88 @@ def Tb_for_archetype(spectra, tdata, nbasis, n_nbh, nleg):

def per_camera_coeff_with_least_square(spectra, tdata, nleg, method=None, n_nbh=None):

"""
This function calculates coefficients for archetype mode in each camera using normal linear algebra matrix solver or BVLS (bounded value least square) method
BVLS described in : https://www.stat.berkeley.edu/~stark/Preprints/bvls.pdf
"""
This function calculates coefficients for archetype mode in each camera using normal linear algebra matrix solver or BVLS (bounded value least square) method
Scipy: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.lsq_linear.html
BVLS described in : https://www.stat.berkeley.edu/~stark/Preprints/bvls.pdf
Parameters
---------------------
Scipy: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.lsq_linear.html
spectra (object): target spectra object
tdata (dict): template data for model fit
nleg (int): number of Legendre polynomials
method (string): 'PCA' or 'bvls'
n_nbh (int): number of nearest best archetypes
Parameters
---------------------
spectra (object): target spectra object
tdata (dict): template data for model fit
nleg (int): number of Legendre polynomials
method (string): 'PCA' or 'bvls'
n_nbh (int): number of nearest best archetypes
Returns
--------------------
coefficients and chi2
Returns
--------------------
coefficients and chi2
"""
"""

ncam = 3 # number of cameras in DESI: b, r, z
ncam = 3 # number of cameras in DESI: b, r, z

flux = np.concatenate([s.flux for s in spectra])
weights = np.concatenate([s.ivar for s in spectra])
wflux = flux*weights
flux = np.concatenate([s.flux for s in spectra])
weights = np.concatenate([s.ivar for s in spectra])
wflux = flux*weights

nbasis = n_nbh+nleg*ncam # n_nbh : for actual physical archetype(s), nleg: number of legendre polynomials, ncamera: number of cameras
nbasis = n_nbh+nleg*ncam # n_nbh : for actual physical archetype(s), nleg: number of legendre polynomials, ncamera: number of cameras

#linear templates in each camera (only the Legendre terms will vary per camera, the lead archetype(s) will remain same for entire spectra)
#linear templates in each camera (only the Legendre terms will vary per camera, the lead archetype(s) will remain same for entire spectra)

Tb = Tb_for_archetype(spectra, tdata, nbasis, n_nbh, nleg)
M = Tb.T.dot(np.multiply(weights[:,None], Tb))
y = Tb.T.dot(wflux)
ret_zcoeff= {'alpha':[], 'b':[], 'r':[], 'z':[]}
Tb = Tb_for_archetype(spectra, tdata, nbasis, n_nbh, nleg)
M = Tb.T.dot(np.multiply(weights[:,None], Tb))
y = Tb.T.dot(wflux)
ret_zcoeff= {'alpha':[], 'b':[], 'r':[], 'z':[]}

# PCA method will use numpy Linear Algebra method to solve the best fit linear equation
if method=='pca':
try:
zcoeff = solve_matrices(M, y, solve_algorithm='PCA', use_gpu=False)
except np.linalg.LinAlgError:
return 9e+99, np.zeros(nbasis)

# BVLS implementation with scipy
if method=='bvls':
bounds = []
for i in range(nbasis):
if i in [j for j in range(n_nbh)]:
bounds.append([0.0, np.inf]) # archetype term(s), these coefficients must be positive
else:
bounds.append([-np.inf, np.inf]) # constant and slope terms in archetype method (can be positive or negative)

bounds = np.array(bounds).T
try:
res = lsq_linear(M, y, bounds=bounds, method='bvls')
zcoeff = res.x
except np.linalg.LinAlgError:
return 9e+99, np.zeros(nbasis)
# PCA method will use numpy Linear Algebra method to solve the best fit linear equation
if method=='pca':
try:
zcoeff = solve_matrices(M, y, solve_algorithm='PCA', use_gpu=False)
except np.linalg.LinAlgError:
return 9e+99, np.zeros(nbasis)

# BVLS implementation with scipy
if method=='bvls':
bounds = []
for i in range(nbasis):
if i in [j for j in range(n_nbh)]:
bounds.append([0.0, np.inf]) # archetype term(s), these coefficients must be positive
else:
bounds.append([-np.inf, np.inf]) # constant and slope terms in archetype method (can be positive or negative)

model = Tb.dot(zcoeff)
zchi2 = np.dot((flux - model)**2, weights)
bounds = np.array(bounds).T
try:
res = lsq_linear(M, y, bounds=bounds, method='bvls')
zcoeff = res.x
except np.linalg.LinAlgError:
return 9e+99, np.zeros(nbasis)

# saving leading archetype coefficients in correct order
ret_zcoeff['alpha'] = [zcoeff[k] for k in range(n_nbh)] # archetype coefficient(s)

if nleg>=1:
split_coeff = np.split(zcoeff[n_nbh:], ncam) # n_camera = 3
old_coeff = {}

# In target spectra redrock saves values as 'b', 'z', 'r'.
# So just re-ordering them here to 'b', 'r', 'z' for easier reading
model = Tb.dot(zcoeff)
zchi2 = np.dot((flux - model)**2, weights)

old_coeff['b']=split_coeff[0]
old_coeff['r']=split_coeff[2]
old_coeff['z']=split_coeff[1]
# saving leading archetype coefficients in correct order
ret_zcoeff['alpha'] = [zcoeff[k] for k in range(n_nbh)] # archetype coefficient(s)

for band in ['b', 'r', 'z']:# 3 cameras
ret_zcoeff[band] = old_coeff[band]

coeff = np.concatenate([c for c in ret_zcoeff.values()])
return zchi2, coeff
if nleg>=1:
split_coeff = np.split(zcoeff[n_nbh:], ncam) # n_camera = 3
old_coeff = {}

# In target spectra redrock saves values as 'b', 'z', 'r'.
# So just re-ordering them here to 'b', 'r', 'z' for easier reading

old_coeff['b']=split_coeff[0]
old_coeff['r']=split_coeff[2]
old_coeff['z']=split_coeff[1]

for band in ['b', 'r', 'z']:# 3 cameras
ret_zcoeff[band] = old_coeff[band]

coeff = np.concatenate([c for c in ret_zcoeff.values()])
return zchi2, coeff

def batch_dot_product_sparse(spectra, tdata, nz, use_gpu):
"""Calculate a batch dot product of the 3 sparse matrices in spectra
Expand Down

0 comments on commit 55b04ef

Please sign in to comment.