Skip to content

Commit

Permalink
make sur npix_array is returned
Browse files Browse the repository at this point in the history
also return covariance
  • Loading branch information
segasai committed Nov 15, 2024
1 parent 8dbb050 commit 63a6480
Show file tree
Hide file tree
Showing 2 changed files with 13 additions and 4 deletions.
5 changes: 3 additions & 2 deletions py/rvspecfit/desi/desi_fit.py
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,8 @@ def __str__(self):
'RV_WARN': 2, # rv is to close to the edge
'RVERR_WARN': 4, # RV error is too large
'PARAM_WARN': 8, # parameters are too close to the edge
'VSINI_WARN': 16 # vsini is too large
'VSINI_WARN': 16, # vsini is too large
'BAD_SPECTRUM': 32 # some issue with the spectrum
}


Expand Down Expand Up @@ -774,7 +775,6 @@ def get_specdata(waves,
medspec = np.nanmedian(np.abs(spec))
if not np.isfinite(medspec) or medspec == 0:
# bail out the spectrum is insane
# TODO make the logic clearer
continue
baddat = ~np.isfinite(spec + curivars)
if mask_dicroic:
Expand Down Expand Up @@ -1042,6 +1042,7 @@ def proc_desi(fname,
else:
rr_z = np.zeros(len(seqid_to_fit)) + np.nan
rr_spectype = np.ma.zeros(len(seqid_to_fit), dtype=str)

subset_ret = subset.copy()
# returned subset to deal with the fact that
# I'll skip some spectra
Expand Down
12 changes: 10 additions & 2 deletions py/rvspecfit/vel_fit.py
Original file line number Diff line number Diff line change
Expand Up @@ -395,6 +395,12 @@ def _minimum_sampler(func,


def _uncertainties_from_hessian(hessian):
"""
Take the hessian and return the uncertainties vector and
covariance matrix
Here I also protect all sorts of failures when encountering
bad hessian
"""
diag_hessian = np.diag(hessian)
inv_diag_hessian = 1. / (diag_hessian + (diag_hessian == 0))
inv_diag_hessian[diag_hessian == 0] = np.inf
Expand All @@ -419,7 +425,7 @@ def _uncertainties_from_hessian(hessian):
diag_err0[sub2] = 0
diag_err = np.sqrt(diag_err0)
diag_err[sub2] = np.nan
return diag_err
return diag_err, hessian_inv


def process(specdata,
Expand Down Expand Up @@ -630,15 +636,17 @@ def hess_func_wrap(p):
hess_step = ndf.MaxStepGenerator(base_step=hess_step)
hessian = ndf.Hessian(hess_func_wrap, step=hess_step)(
[ret['param'][_] for _ in specParamNames])
diag_err = _uncertainties_from_hessian(hessian)
diag_err, covar_mat = _uncertainties_from_hessian(hessian)
ret['param_err'] = dict(zip(specParamNames, diag_err))
ret['param_covar'] = covar_mat
ret['minimize_success'] = minimize_success

ret['yfit'] = outp['models']
ret['raw_models'] = outp['raw_models']
ret['chisq'] = outp['chisq']
ret['logl'] = outp['logl']
ret['chisq_array'] = outp['chisq_array']
ret['npix_array'] = outp['npix_array']
t6 = time.time()
logging.debug('Timings process: %.4f %.4f %.4f %.4f, %.4f %.4f' %
(t1 - t0, t2 - t1, t3 - t2, t4 - t3, t5 - t4, t6 - t5))
Expand Down

0 comments on commit 63a6480

Please sign in to comment.