diff --git a/.DS_Store b/.DS_Store new file mode 100644 index 00000000..e59eb705 Binary files /dev/null and b/.DS_Store differ diff --git a/.gitignore b/.gitignore new file mode 100644 index 00000000..223a9456 --- /dev/null +++ b/.gitignore @@ -0,0 +1,163 @@ +# Byte-compiled / optimized / DLL files +__pycache__/ +*.py[cod] +*$py.class + +# C extensions +*.so + +# Distribution / packaging +.Python +build/ +develop-eggs/ +dist/ +downloads/ +eggs/ +.eggs/ +lib/ +lib64/ +parts/ +sdist/ +var/ +wheels/ +share/python-wheels/ +*.egg-info/ +.installed.cfg +*.egg +MANIFEST + +# PyInstaller +# Usually these files are written by a python script from a template +# before PyInstaller builds the exe, so as to inject date/other infos into it. +*.manifest +*.spec + +# Installer logs +pip-log.txt +pip-delete-this-directory.txt + +# Unit test / coverage reports +htmlcov/ +.tox/ +.nox/ +.coverage +.coverage.* +.cache +nosetests.xml +coverage.xml +*.cover +*.py,cover +.hypothesis/ +.pytest_cache/ +cover/ + +# Translations +*.mo +*.pot + +# Django stuff: +*.log +local_settings.py +db.sqlite3 +db.sqlite3-journal + +# Flask stuff: +instance/ +.webassets-cache + +# Scrapy stuff: +.scrapy + +# Sphinx documentation +docs/_build/ + +# PyBuilder +.pybuilder/ +target/ + +# Jupyter Notebook +.ipynb_checkpoints + +# IPython +profile_default/ +ipython_config.py + +# pyenv +# For a library or package, you might want to ignore these files since the code is +# intended to run in multiple environments; otherwise, check them in: +# .python-version + +# pipenv +# According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control. +# However, in case of collaboration, if having platform-specific dependencies or dependencies +# having no cross-platform support, pipenv may install dependencies that don't work, or not +# install all needed dependencies. +#Pipfile.lock + +# poetry +# Similar to Pipfile.lock, it is generally recommended to include poetry.lock in version control. +# This is especially recommended for binary packages to ensure reproducibility, and is more +# commonly ignored for libraries. +# https://python-poetry.org/docs/basic-usage/#commit-your-poetrylock-file-to-version-control +#poetry.lock + +# pdm +# Similar to Pipfile.lock, it is generally recommended to include pdm.lock in version control. +#pdm.lock +# pdm stores project-wide configurations in .pdm.toml, but it is recommended to not include it +# in version control. +# https://pdm.fming.dev/#use-with-ide +.pdm.toml + +# PEP 582; used by e.g. github.com/David-OConnor/pyflow and github.com/pdm-project/pdm +__pypackages__/ + +# Celery stuff +celerybeat-schedule +celerybeat.pid + +# SageMath parsed files +*.sage.py + +# Environments +.env +.venv +env/ +venv/ +ENV/ +env.bak/ +venv.bak/ + +# Spyder project settings +.spyderproject +.spyproject + +# Rope project settings +.ropeproject + +# mkdocs documentation +/site + +# mypy +.mypy_cache/ +.dmypy.json +dmypy.json + +# Pyre type checker +.pyre/ + +# pytype static type analyzer +.pytype/ + +# Cython debug symbols +cython_debug/ + +# PyCharm +# JetBrains specific template is maintained in a separate JetBrains.gitignore that can +# be found at https://github.com/github/gitignore/blob/main/Global/JetBrains.gitignore +# and can be added to the global gitignore or merged into this file. For a more nuclear +# option (not recommended) you can uncomment the following to ignore the entire idea folder. +#.idea/ + +# DSstore +smh/.DS_Store diff --git a/README.md b/README.md index 630af3ae..6d4e3b45 100644 --- a/README.md +++ b/README.md @@ -9,7 +9,7 @@ Authors ------- - Andrew R. Casey (Monash) - Alex Ji (University of Chicago) - - Erika Holmbeck (Carnegie Observatories) + - Erika Holmbeck (Lawrence Livermore National Laboratory) Installation ------------ diff --git a/emh_smhr3_versions.txt b/emh_smhr3_versions.txt new file mode 100644 index 00000000..e0c23f06 --- /dev/null +++ b/emh_smhr3_versions.txt @@ -0,0 +1,94 @@ +# packages in environment at /Users/holmbeck/opt/anaconda3/envs/smhr-py3: +# +# Name Version Build Channel +appnope 0.1.2 py38hecd8cb5_1001 +astropy 5.0.4 py38h67323c0_0 +asttokens 2.0.5 pyhd3eb1b0_0 +backcall 0.2.0 pyhd3eb1b0_0 +blas 1.0 mkl +brotlipy 0.7.0 py38h9ed2024_1003 +ca-certificates 2022.5.18.1 h033912b_0 conda-forge +certifi 2022.5.18.1 py38h50d1736_0 conda-forge +cffi 1.15.0 py38hc55c11b_1 +charset-normalizer 2.0.4 pyhd3eb1b0_0 +cryptography 37.0.1 py38hf6deb26_0 +cycler 0.11.0 pyhd3eb1b0_0 +decorator 5.1.1 pyhd3eb1b0_0 +executing 0.8.3 pyhd3eb1b0_0 +freetype 2.11.0 hd8bbffd_0 +icu 68.2 he49afe7_0 conda-forge +idna 3.3 pyhd3eb1b0_0 +intel-openmp 2021.4.0 hecd8cb5_3538 +ipython 8.3.0 py38hecd8cb5_0 +jedi 0.18.1 py38hecd8cb5_1 +jpeg 9e h5eb16cf_1 conda-forge +julia 0.6.1 pypi_0 pypi +kiwisolver 1.4.2 py38he9d5cce_0 +krb5 1.19.2 hcfbf3a7_3 conda-forge +libclang 11.1.0 default_he082bbe_1 conda-forge +libcxx 12.0.0 h2f01273_0 +libedit 3.1.20191231 h0678c8f_2 conda-forge +libffi 3.3 hb1e8313_2 +libgfortran 3.0.1 h93005f0_2 +libiconv 1.16 haf1e3a3_0 conda-forge +libllvm11 11.1.0 hd011deb_2 conda-forge +libpng 1.6.37 ha441bb4_0 +libpq 13.3 hea3049e_0 conda-forge +libxml2 2.9.12 h93ec3fd_0 conda-forge +libxslt 1.1.33 h5739fc3_2 conda-forge +lz4-c 1.9.3 he49afe7_1 conda-forge +matplotlib 3.1.3 py38_0 +matplotlib-base 3.1.3 py38h9aa3819_0 +matplotlib-inline 0.1.2 pyhd3eb1b0_2 +mkl 2021.4.0 hecd8cb5_637 +mkl-service 2.4.0 py38h9ed2024_0 +mkl_fft 1.3.1 py38h4ab4a9b_0 +mkl_random 1.2.2 py38hb2f4e1b_0 +mysql-common 8.0.25 h694c41f_2 conda-forge +mysql-libs 8.0.25 h115446f_2 conda-forge +ncurses 6.3 hca72f7f_2 +nspr 4.32 hcd9eead_1 conda-forge +nss 3.69 h31e2bf1_1 conda-forge +numpy 1.22.3 py38h2e5f0a9_0 +numpy-base 1.22.3 py38h3b1a694_0 +openssl 1.1.1o hfe4f2af_0 conda-forge +packaging 21.3 pyhd3eb1b0_0 +parso 0.8.3 pyhd3eb1b0_0 +pexpect 4.8.0 pyhd3eb1b0_3 +pickleshare 0.7.5 pyhd3eb1b0_1003 +pip 21.2.4 py38hecd8cb5_0 +prompt-toolkit 3.0.20 pyhd3eb1b0_0 +ptyprocess 0.7.0 pyhd3eb1b0_2 +pure_eval 0.2.2 pyhd3eb1b0_0 +pycparser 2.21 pyhd3eb1b0_0 +pyerfa 2.0.0 py38h9ed2024_0 +pygments 2.11.2 pyhd3eb1b0_0 +pyopenssl 22.0.0 pyhd3eb1b0_0 +pyparsing 3.0.4 pyhd3eb1b0_0 +pyside2 5.15.2.1 pypi_0 pypi +pysocks 1.7.1 py38_1 +python 3.8.13 hdfd78df_0 +python-dateutil 2.8.2 pyhd3eb1b0_0 +python.app 3 py38hca72f7f_0 +python_abi 3.8 2_cp38 conda-forge +pyyaml 6.0 py38hca72f7f_1 +qt 5.12.9 h126340a_4 conda-forge +readline 8.1.2 hca72f7f_1 +requests 2.27.1 pyhd3eb1b0_0 +scipy 1.7.3 py38h8c7af03_0 +setuptools 61.2.0 py38hecd8cb5_0 +shiboken2 5.15.2.1 pypi_0 pypi +six 1.16.0 pyhd3eb1b0_1 +smh 0.2 dev_0 +sqlite 3.38.3 h707629a_0 +stack_data 0.2.0 pyhd3eb1b0_0 +tk 8.6.12 h5d9f67b_0 +tornado 6.1 py38h9ed2024_0 +traitlets 5.1.1 pyhd3eb1b0_0 +urllib3 1.26.9 py38hecd8cb5_0 +wcwidth 0.2.5 pyhd3eb1b0_0 +wheel 0.37.1 pyhd3eb1b0_0 +xz 5.2.5 hca72f7f_1 +yaml 0.2.5 h0d85af4_2 conda-forge +zlib 1.2.12 h4dc903c_2 +zstd 1.5.0 h582d3a0_0 conda-forge diff --git a/smh/.DS_Store b/smh/.DS_Store index b9f4aab4..92ab8b44 100644 Binary files a/smh/.DS_Store and b/smh/.DS_Store differ diff --git a/smh/gui/base.py b/smh/gui/base.py index 7b2affe8..6398d601 100644 --- a/smh/gui/base.py +++ b/smh/gui/base.py @@ -125,6 +125,10 @@ def __init__(self, parent, session=None, self.session = session self.label_ymin = label_ymin self.label_ymax = label_ymax + # E. Holmbeck changed colors: + self.acceptable_color = "#d92653" + self.unacceptable_color = "#37ae91" + self.comparison_spectrum = comparison_spectrum @@ -191,7 +195,8 @@ def __init__(self, parent, session=None, np.nan, np.nan, np.nan, color="blue", linestyle=':', lw=1), "model_masks": [], "nearby_lines": [], - "model_fit": self.ax_spectrum.plot([np.nan], [np.nan], c="r")[0], + "model_fit": self.ax_spectrum.plot([np.nan], [np.nan], c=self.acceptable_color)[0], + "model_none": self.ax_spectrum.plot([np.nan], [np.nan], c='teal')[0], "model_residual": self.ax_residual.plot( [np.nan], [np.nan], c="k", drawstyle="steps-mid")[0], "interactive_mask": [ @@ -225,7 +230,7 @@ def reset(self): self._lines["spectrum"].set_drawstyle(drawstyle) self._lines["comparison_spectrum"].set_drawstyle(drawstyle) for key in ["spectrum", "transitions_center_main", "transitions_center_residual", - "model_fit", "model_residual"]: + "model_fit", "model_none", "model_residual"]: self._lines[key].set_data([np.nan],[np.nan]) self.label_lines(None) def new_session(self, session): @@ -366,7 +371,7 @@ def key_press_zoom(self, event): self.ax_spectrum.set_ylim(ylim) self.draw() return None - + def spectrum_left_mouse_press(self, event): """ Listener for if mouse button pressed in spectrum or residual axis @@ -627,6 +632,9 @@ def _plot_model(self): pass selected_model = self.selected_model + try: none_x,none_y = selected_model.metadata["zero_abundance"] + except: none_x = none_y = np.nan + try: (named_p_opt, cov, meta) = selected_model.metadata["fitted_result"] @@ -643,6 +651,7 @@ def _plot_model(self): except KeyError: meta = {} self._lines["model_fit"].set_data([np.nan], [np.nan]) + self._lines["model_none"].set_data([np.nan], [np.nan]) self._lines["model_residual"].set_data([np.nan], [np.nan]) else: @@ -652,7 +661,8 @@ def _plot_model(self): self._lines["model_fit"].set_data(meta[plotxkey], meta[plotykey]) self._lines["model_fit"].set_linestyle("-" if self.selected_model.is_acceptable else "--") - self._lines["model_fit"].set_color("r" if self.selected_model.is_acceptable else "b") + self._lines["model_fit"].set_color(self.acceptable_color if self.selected_model.is_acceptable else self.unacceptable_color) + self._lines["model_none"].set_data(none_x, none_y) self._lines["model_residual"].set_data(meta["model_x"], meta["residual"]) # Model yerr. @@ -661,7 +671,7 @@ def _plot_model(self): meta["model_x"], meta["model_y"] + meta["model_yerr"], meta["model_y"] - meta["model_yerr"], - facecolor="r" if self.selected_model.is_acceptable else "b", + facecolor=self.acceptable_color if self.selected_model.is_acceptable else self.unacceptable_color, edgecolor="none", alpha=0.5) # Model masks due to nearby lines. @@ -673,9 +683,9 @@ def _plot_model(self): except IndexError: self._lines["nearby_lines"].append([ self.ax_spectrum.axvspan(np.nan, np.nan, - facecolor="b", edgecolor="none", alpha=0.25), + facecolor=self.unacceptable_color, edgecolor="none", alpha=0.25), self.ax_residual.axvspan(np.nan, np.nan, - facecolor="b", edgecolor="none", alpha=0.25) + facecolor=self.unacceptable_color, edgecolor="none", alpha=0.25) ]) patches = self._lines["nearby_lines"][-1] @@ -743,18 +753,29 @@ def __init__(self, parent, xattr, yattr, error_styles=None, linefit_styles=None, linemean_styles=None, + sigma=2.5, **kwargs): assert xattr in self.allattrs, xattr assert yattr in self.allattrs, yattr self.xattr = xattr - self.yattr = yattr + self.yattr = yattr + # E. Holmbeck skipped this for now. + ''' if error_styles is not None: assert (e_xattr is not None) or (e_yattr is not None), "Must specify e_xattr and/or e_yattr for error_styles" assert (exattr is None) or (exattr in self.allattrs), exattr assert (eyattr is None) or (eyattr in self.allattrs), eyattr self.exattr = exattr self.eyattr = eyattr - + ''' + if error_styles is not None: + self.exattr = exattr + self.eyattr = eyattr + else: + self.exattr = None + self.eyattr = None + self.sigma = sigma + super(SMHScatterplot, self).__init__(parent=parent, **kwargs) @@ -774,7 +795,7 @@ def __init__(self, parent, xattr, yattr, error_styles = [None for f in filters] else: assert len(filters)==len(error_styles) - raise NotImplementedError("Need to implement error bar graphics updating!") + #raise NotImplementedError("Need to implement error bar graphics updating!") if linefit_styles is None: linefit_styles = [None for f in filters] else: @@ -789,6 +810,9 @@ def __init__(self, parent, xattr, yattr, error_objs = [] linefit_objs = [] linemean_objs = [] + # E. Holmbeck added shading + fillmean_objs = [] + for filt, point_kw, error_kw, linefit_kw, linemean_kw in zip( filters, point_styles, error_styles, linefit_styles, linemean_styles): if point_kw is None: point_objs.append(None) @@ -802,7 +826,8 @@ def __init__(self, parent, xattr, yattr, error_objs.append( self.ax.errorbar(np.nan * np.ones(2), np.nan * np.ones(2), yerr=np.nan * np.ones((2, 2)), - fmt=None,zorder=-10, + fmt='o', # E. Holmbeck changed this. + zorder=-10, **error_kw)) if linefit_kw is None: linefit_objs.append(None) @@ -810,10 +835,15 @@ def __init__(self, parent, xattr, yattr, linefit_objs.append( self.ax.plot([np.nan], [np.nan], **linefit_kw)[0]) - if linemean_kw is None: linemean_objs.append(None) + if linemean_kw is None: + linemean_objs.append(None) + fillmean_objs.append(None) else: linemean_objs.append( self.ax.axhline(np.nan, **linemean_kw)) + # E. Holmbeck added shading + fillmean_objs.append( + self.ax.fill_between([0,0],[0,0],[0,0], color=linemean_kw["color"], alpha=0.10, lw=0)) ## Save graphic objects self._selected_points = self.ax.scatter([], [], @@ -824,7 +854,9 @@ def __init__(self, parent, xattr, yattr, self._errors = error_objs self._linefits = linefit_objs self._linemeans = linemean_objs - self._graphics = list(zip(self._filters, self._points, self._errors, self._linefits, self._linemeans)) + # E. Holmbeck added shading + self._fillmeans = fillmean_objs + self._graphics = list(zip(self._filters, self._points, self._errors, self._linefits, self._linemeans, self._fillmeans)) ## Connect Interactivity if enable_zoom: @@ -840,17 +872,38 @@ def __init__(self, parent, xattr, yattr, self.update_scatterplot() self.update_selected_points(True) + ''' + def zero_out_plot(self,linefit, linemean, fillmean): + #for plot_item in self.ax.collections+self.ax.lines: + # plot_item.remove() + if linefit is not None: + linefit.set_data([np.nan], [np.nan]) + if linemean is not None: + linemean.set_data([np.nan], [np.nan]) + if fillmean is not None: + path = fillmean.get_paths()[0] + v_x = np.hstack([0.0]*len(path.vertices)) + vertices = np.vstack([v_x,v_x]).T + path.vertices = vertices + self.draw() + ''' def sizeHint(self): return QtCore.QSize(125,100) def minimumSizeHint(self): return QtCore.QSize(10,10) def reset(self): self._selected_points.set_offsets(np.array([np.nan, np.nan]).T) - for filt, point, error, linefit, linemean in self._graphics: + for filt, point, error, linefit, linemean, fillmean in self._graphics: if point is not None: point.set_offsets(np.array([np.nan, np.nan]).T) if error is not None: pass # TODO!!! if linefit is not None: linefit.set_data([np.nan],[np.nan]) if linemean is not None: linemean.set_data([0,1],[np.nan,np.nan]) + if fillmean is not None: #pass # TODO!!! + path = fillmean.get_paths()[0] + v_x = np.hstack([0.0]*len(path.vertices)) + vertices = np.vstack([v_x,v_x]).T + path.vertices = vertices + def linkToTable(self, tableview): """ view for selection; model for data @@ -905,59 +958,105 @@ def update_scatterplot(self, redraw=False): xs, ys, exs, eys = [], [], [], [] Nrows = self.tablemodel.rowCount() if Nrows==0: return None + spectral_models = self.tablemodel.get_models_from_rows(np.arange(Nrows)) for i in range(Nrows): ix = self._ix(i, self.xcol) x = self._load_value_from_table(ix) ix = self._ix(i, self.ycol) y = self._load_value_from_table(ix) - if self.excol is None: ex = np.nan - else: + if self.excol is None: + # TODO: Add x-error, which is a function of exattr... + ex = np.nan + else: # This will never be called (we don't use x-err) ix = self._ix(i, self.excol) ex = self._load_value_from_table(ix) - if self.eycol is None: ey = np.nan + if self.eycol is None: + # E. Holmbeck changed; just make sure it's right + #ey = np.nan + ey = spectral_models[i].abundance_uncertainties else: ix = self._ix(i, self.eycol) ey = self._load_value_from_table(ix) + xs.append(x) ys.append(y) exs.append(ex) - eys.append(ey) + eys.append(np.nan if ey is None else ey) xs = np.array(xs); ys = np.array(ys); exs = np.array(exs); eys = np.array(eys) + valids = [] for filt in self._filters: valids.append([filt(sm) for sm in spectral_models]) valids = np.atleast_2d(np.array(valids, dtype=bool)) - for ifilt,(filt, point, error, linefit, linemean) in enumerate(self._graphics): + + # All four datatypes: acceptable, not acceptable, user flag, upper limit + for ifilt,(filt, point, error, linefit, linemean, fillmean) in enumerate(self._graphics): valid = valids[ifilt,:] nonzero = valid.sum() > 0 x, y = xs[valid], ys[valid] + if np.all(np.isnan(xs)): + self.reset() + self.draw() + + # E. Holmbeck added; only works for yerr though + ex, ey = exs[valid], eys[valid] if point is not None: if nonzero: point.set_offsets(np.array([x,y]).T) else: point.set_offsets(np.array([np.nan,np.nan]).T) if error is not None: ## TODO not doing anything with error bars right now - pass - if nonzero and ((linefit is not None) or (linemean is not None)): - ## TODO: Figure out how best to save and return info about the lines - ## For now, just refitting whenever needed - try: - m, b, medy, stdy, stdm, N = utils.fit_line(x, y, None) - except ValueError as e: - return np.nan, np.nan, np.nan, np.nan, np.nan, np.nan - #xlim = np.array(self.ax.get_xlim()) - xlim = np.array([x.min(), x.max()]) - if (linefit is not None) and nonzero: - linefit.set_data(xlim, m*xlim + b) - else: - linefit.set_data([np.nan], [np.nan]) - if (linemean is not None) and nonzero: - linemean.set_data([0,1], [medy, medy]) - else: - linemean.set_data([0,1], [np.nan, np.nan]) - style_utils.relim_axes(self.ax) + # E. Holmbeck added this; shamelessly stolen: https://stackoverflow.com/questions/25210723/matplotlib-set-data-for-errorbar-plot + # Still very buggy. + ln, (erry_top, erry_bot), (barsy,) = error.lines + x_base = x + y_base = y + yerr_top = y_base + ey + yerr_bot = y_base - ey + erry_top.set_xdata(x_base) + erry_bot.set_xdata(x_base) + erry_top.set_ydata(yerr_top) + erry_bot.set_ydata(yerr_bot) + new_segments_y = [np.array([[x, yt], [x,yb]]) for x, yt, yb in zip(x_base, yerr_top, yerr_bot)] + barsy.set_segments(new_segments_y) + + # Reset axes before adding lines. + xlim,ylim = style_utils.relim_axes(self.ax) + xlim = np.array(xlim) self.reset_zoom_limits() - if redraw: self.draw() + + for ifilt,(filt, point, error, linefit, linemean, fillmean) in enumerate(self._graphics): + valid = valids[ifilt,:] + nonzero = valid.sum() > 0 + if not nonzero: + continue + + x, y, ey = xs[valid], ys[valid], eys[valid] + if np.all(np.isnan(x)): continue + + m,b,medy,stdy,stdm,N = utils.fit_line(x, y, ey) + + if linefit is not None: + linefit.set_data(xlim, m*xlim + b) + if linemean is not None: + linemean.set_data([0,1], [medy, medy]) + # E. Holmbeck added; shamelessly stolen: https://stackoverflow.com/questions/16120801/matplotlib-animate-fill-between-shape + if fillmean is not None: + path = fillmean.get_paths()[0] + y0new = [medy - self.sigma*stdy]*2 + y1new = [medy + self.sigma*stdy]*2 + #self.edit_sigma.text() + v_x = np.hstack([xlim[0],xlim,xlim[-1],xlim[::-1],xlim[0]]) + v_y = np.hstack([y1new[0],y0new,y0new[-1],y1new[::-1],y1new[0]]) + vertices = np.vstack([v_x,v_y]).T + path.vertices = vertices + #_,ylim = style_utils.relim_axes(self.ax, values_only=True) + #self.ax.set_ylim(*ylim) + self.ax.set_ylim([min([ylim[0],medy - (self.sigma+1)*stdy]), + max([ylim[1],medy + (self.sigma+1)*stdy])]) + + if redraw: + self.draw() return None def update_selected_points(self, redraw=False): if self.tableview is None or self.tablemodel is None: return None @@ -1272,7 +1371,7 @@ def set_moog_option(self, key, value): if 'moog_opts' not in spectral_model.metadata: spectral_model.metadata['moog_opts'] = {key: value} else: - spectral_model.metadata['moog_opts'][key] = value + spectral_model.metadata['moog_opts'][key] = value if "fitted_result" in spectral_model.metadata: num_fit += 1 try: diff --git a/smh/gui/chemical_abundances.py b/smh/gui/chemical_abundances.py index a0e24680..eced76f4 100644 --- a/smh/gui/chemical_abundances.py +++ b/smh/gui/chemical_abundances.py @@ -43,7 +43,7 @@ QtGui2.QFont.insertSubstitution(*substitute) _QFONT = QtGui2.QFont("Helvetica Neue", 10) -_ROWHEIGHT = 20 +_ROWHEIGHT = 30 DOUBLE_CLICK_INTERVAL = 0.1 # MAGIC HACK PICKER_TOLERANCE = 10 # MAGIC HACK @@ -73,8 +73,9 @@ def __init__(self, parent): self.figure.mpl_connect("key_press_event", self.key_press_selectcheck) self.figure.mpl_connect("key_press_event", self.key_press_automask_sigma) ## Stuff for extra synthesis - self.extra_spec_1 = self.ax_spectrum.plot([np.nan],[np.nan], ls='-', color='#cea2fd', lw=1.5, zorder=9999)[0] - self.extra_spec_2 = self.ax_spectrum.plot([np.nan],[np.nan], ls='-', color='#ffb07c', lw=1.5, zorder=9999)[0] + self.extra_spec_1 = self.ax_spectrum.plot([np.nan],[np.nan], ls='-', color='#ff531a', lw=1.5, zorder=9998)[0] + self.extra_spec_2 = self.ax_spectrum.plot([np.nan],[np.nan], ls='-', color='#ffc34b', lw=1.5, zorder=9998)[0] + #self.extra_spec_none = self.ax_spectrum.plot([np.nan],[np.nan], ls='-', color='teal', lw=1.2, zorder=9999)[0] ################ # BOTTOM @@ -111,6 +112,8 @@ def __init__(self, parent): self._currently_plotted_element = "All" self.new_session_loaded() + # E. Holmbeck: don't mind me... + self.acceptables_only = True def init_tab(self): """ @@ -185,7 +188,7 @@ def _create_measurement_list(self): self.btn_fit_all.setText("Fit all EW") self.btn_fit_all.setSizePolicy(sp) self.btn_measure_all = QtGui.QPushButton(self) - self.btn_measure_all.setText("Measure all acceptable EW") + self.btn_measure_all.setText("Abundances from EW") self.btn_measure_all.setSizePolicy(sp) hbox.addWidget(self.btn_fit_all) hbox.addWidget(self.btn_measure_all) @@ -299,7 +302,7 @@ def _create_checkcombo_in_hbox(parent, text): vbox_lhs.addLayout(hbox) hbox, label, line = _create_line_in_hbox(self.tab_profile, "Fit window", - 0, 1000, 1) + -4, 1000, 1) self.edit_fit_window = line vbox_lhs.addLayout(hbox) @@ -394,7 +397,7 @@ def _create_checkcombo_in_hbox(parent, text): vbox_lhs.addLayout(hbox) hbox, label, line = _create_line_in_hbox(self.tab_synthesis, "Fit window", - 0, 1000, 1) + -4, 1000, 1) self.edit_fit_window_2 = line vbox_lhs.addLayout(hbox) @@ -493,15 +496,29 @@ def _create_checkcombo_in_hbox(parent, text): #vbox_rhs.addItem(QtGui.QSpacerItem(20,20,QtGui.QSizePolicy.Minimum, # QtGui.QSizePolicy.Minimum)) + # E. Holmbeck edited + hbox = QtGui.QHBoxLayout() + self.btn_update_abund_table = QtGui.QPushButton(self.tab_synthesis) + self.btn_update_abund_table.setText("Update Ab. Table") + #vbox_rhs.addWidget(self.btn_update_abund_table) + hbox.addWidget(self.btn_update_abund_table) self.btn_fit_synth = QtGui.QPushButton(self.tab_synthesis) self.btn_fit_synth.setText("Fit Model") - vbox_rhs.addWidget(self.btn_fit_synth) - - self.btn_update_abund_table = QtGui.QPushButton(self.tab_synthesis) - self.btn_update_abund_table.setText("Update Abundance Table") - vbox_rhs.addWidget(self.btn_update_abund_table) - + #vbox_rhs.addWidget(self.btn_fit_synth) + hbox.addWidget(self.btn_fit_synth) + vbox_rhs.addLayout(hbox) + # E. Holmbeck added + hbox = QtGui.QHBoxLayout() + self.btn_fit_all_synth = QtGui.QPushButton(self.tab_synthesis) + self.btn_fit_all_synth.setText("Resynth") + hbox.addWidget(self.btn_fit_all_synth) + # E. Holmbeck added + self.btn_fit_rest_synth = QtGui.QPushButton(self.tab_synthesis) + self.btn_fit_rest_synth.setText("Synth Rest") + hbox.addWidget(self.btn_fit_rest_synth) + vbox_rhs.addLayout(hbox) + hbox = QtGui.QHBoxLayout() self.btn_clear_masks_2 = QtGui.QPushButton(self.tab_synthesis) self.btn_clear_masks_2.setText("Clear Masks") @@ -674,6 +691,12 @@ def _connect_profile_signals(self): self.clicked_btn_clear_masks) self.btn_export_synth.clicked.connect( self.clicked_export_synthesis) + + self.btn_fit_all_synth.clicked.connect( + self.clicked_fit_all_synth) + + self.btn_fit_rest_synth.clicked.connect( + self.clicked_fit_rest_synth) self._synth_signals = [ (self.edit_view_window_2.textChanged,self.update_edit_view_window_2), @@ -700,6 +723,8 @@ def _connect_profile_signals(self): (self.checkbox_upper_limit_2.stateChanged,self.clicked_checkbox_upper_limit_2), (self.btn_find_upper_limit.clicked,self.clicked_btn_find_upper_limit), (self.btn_fit_synth.clicked,self.fit_one), + (self.btn_fit_all_synth.clicked,self.clicked_fit_all_synth), + (self.btn_fit_rest_synth.clicked,self.clicked_fit_rest_synth), (self.btn_update_abund_table.clicked,self.clicked_btn_update_abund_table), (self.btn_clear_masks_2.clicked,self.clicked_btn_clear_masks), (self.btn_export_synth.clicked,self.clicked_export_synthesis) @@ -710,6 +735,8 @@ def populate_filter_combo_box(self): box = self.filter_combo_box box.clear() box.addItem("All") + box.addItem("Synthesized Lines") + box.addItem("EW Lines") all_species = set([]) for spectral_model in self.full_measurement_model.spectral_models: @@ -725,11 +752,12 @@ def populate_filter_combo_box(self): assert species == utils.element_to_species(elem) box.addItem(elem) + # E. Holmbeck added filters for Synth and EW def filter_combo_box_changed(self): elem = self.filter_combo_box.currentText() # Update the filter self.measurement_model.beginResetModel() - if self._currently_plotted_element not in ["All", "", "None"]: + if self._currently_plotted_element not in [None, "All", "", "None"]: try: self.measurement_model.delete_filter_function(self._currently_plotted_element) except KeyError as e: @@ -737,8 +765,20 @@ def filter_combo_box_changed(self): logger.debug(e) logger.debug(self.measurement_model.filter_functions) raise - if elem in [None, "", "All"]: + if elem in [None, "", "All", "None"]: + self.element_summary_text.setText("") + elif elem == "Synthesized Lines": + self.element_summary_text.setText("") + def filter_function(model): + if isinstance(model, SpectralSynthesisModel): + return model.species + self.measurement_model.add_filter_function(elem, filter_function) + elif elem == "EW Lines": self.element_summary_text.setText("") + def filter_function(model): + if isinstance(model, ProfileFittingModel): + return model.species + self.measurement_model.add_filter_function(elem, filter_function) else: species = utils.element_to_species(elem) def filter_function(model): @@ -747,6 +787,7 @@ def filter_function(model): elif isinstance(model, SpectralSynthesisModel): return np.any([species in specie for specie in model.species]) self.measurement_model.add_filter_function(elem, filter_function) + self._currently_plotted_element = elem self.measurement_model.endResetModel() self.summarize_current_table() @@ -763,12 +804,26 @@ def calculate_FeH(self): return np.nan return self.FeH + # E. Holmbeck added the synth and ew lines; TODO: make less hacky def summarize_current_table(self): elem = self.filter_combo_box.currentText() if elem is None or elem == "" or elem == "All": N = self.measurement_model.rowCount() self.element_summary_text.setText("N={} lines".format(N)) return None + elif elem == "Synthesized Lines": + filtered = [] + for line in self.parent.session.spectral_models: + if isinstance(line, SpectralSynthesisModel): filtered.append(line) + self.element_summary_text.setText("N={} lines".format(len(filtered))) + return None + elif elem == "EW Lines": + filtered = [] + for line in self.parent.session.spectral_models: + if isinstance(line, ProfileFittingModel): filtered.append(line) + self.element_summary_text.setText("N={} lines".format(len(filtered))) + return None + summary_dict = self.parent.session.summarize_spectral_models(organize_by_element=False) species = utils.element_to_species(elem) if species not in summary_dict: @@ -795,8 +850,8 @@ def refresh_table(self): self.calculate_FeH() return None - def refresh_plots(self): - self.update_spectrum_figure(True) + def refresh_plots(self, row=None): + self.update_spectrum_figure(True, row=row) return None def fit_all_profiles(self): @@ -876,6 +931,21 @@ def measure_all(self): pass return None + def fit_none(self, spectral_model): + if not isinstance(spectral_model, SpectralSynthesisModel): return + # E. Holmbeck added a "none" line + extra_abundances = self.synth_abund_table_model.get_extra_abundances() + if extra_abundances is None: + for i, elem in enumerate(spectral_model.elements): + abundances_none = deepcopy(spectral_model.metadata["rt_abundances"]) + for i, elem in enumerate(spectral_model.elements): + abundances_none[elem] = -50.0 + + x, y = spectral_model.get_synth(abundances_none) + #self.extra_spec_none.set_data([x,y]) + spectral_model.metadata["zero_abundance"] = [x,y] + + def fit_one(self): spectral_model, proxy_index, index = self._get_selected_model(True) if spectral_model is None: return None @@ -888,11 +958,15 @@ def fit_one(self): self.measurement_view.update_row(proxy_index.row()) self.summarize_current_table() self.update_fitting_options() + self.fit_none(spectral_model) self.refresh_plots() + if self.parent.session.setting("bring_to_top_after_fit", False): self.parent.raise_() self.parent.activateWindow() self.parent.showNormal() + + self.figure.draw() return None def measure_one(self): @@ -911,8 +985,8 @@ def synthesize_current_model(self): # When we get selected model, it erases the extra_abundances. # So let's cache it then put it back in... extra_abundances = self.synth_abund_table_model.get_extra_abundances() - spectral_model, proxy_index, index = self._get_selected_model(True) + if spectral_model is None: return None spectral_model.update_fit_after_parameter_change() self.measurement_view.update_row(proxy_index.row()) @@ -942,7 +1016,11 @@ def synthesize_current_model(self): self.extra_spec_1.set_data([x,y]) x, y = spectral_model.get_synth(abundances2) self.extra_spec_2.set_data([x,y]) + # E. Holmbeck added a "none" line + self.fit_none(spectral_model) self.figure.draw() + + return None def _check_for_spectral_models(self): @@ -980,41 +1058,70 @@ def _get_selected_model(self, full_output=False): model = self.parent.session.metadata["spectral_models"][index] return (model, proxy_index, index) if full_output else model - def selected_model_changed(self): - self.update_fitting_options() - self.update_spectrum_figure(True) + def selected_model_changed(self, **kwargs): + self.update_fitting_options(**kwargs) + self.update_spectrum_figure(True, **kwargs) return None - def update_spectrum_figure(self, redraw=False, reset_limits=True): + def update_spectrum_figure(self, redraw=False, reset_limits=True, row=None, selected_model=None): ## If synthesis, label selected lines self.extra_spec_1.set_data([[np.nan], [np.nan]]) self.extra_spec_2.set_data([[np.nan], [np.nan]]) - try: - selected_model = self._get_selected_model() - except IndexError: - selected_transitions = None - label_rv = None - else: - if isinstance(selected_model, SpectralSynthesisModel): - selected_elem = self.synth_abund_table.get_selected_element() - if selected_elem is not None: - transitions = selected_model.transitions - ii = np.logical_or(transitions["elem1"] == selected_elem, - transitions["elem2"] == selected_elem) - if np.sum(ii) != 0: - selected_transitions = transitions[ii] - redraw=True # force redraw - reset_limits=False - label_rv = selected_model.metadata["manual_rv"] + #self.extra_spec_none.set_data([[np.nan], [np.nan]]) + + if selected_model is None: + try: + if row==None: + selected_model = self._get_selected_model() + else: + selected_model = self.parent.session.metadata["spectral_models"][row] + except IndexError: + selected_transitions = None + label_rv = None + + else: + if isinstance(selected_model, SpectralSynthesisModel): + selected_elem = self.synth_abund_table.get_selected_element() + if selected_elem is not None: + transitions = selected_model.transitions + ii = np.logical_or(transitions["elem1"] == selected_elem, + transitions["elem2"] == selected_elem) + if np.sum(ii) != 0: + selected_transitions = transitions[ii] + redraw=True # force redraw + reset_limits=False + label_rv = selected_model.metadata["manual_rv"] + else: + selected_transitions = None + label_rv = None else: selected_transitions = None label_rv = None else: selected_transitions = None label_rv = None + + else: + selected_elem = selected_model.elements[0] + if selected_elem is not None: + transitions = selected_model.transitions + ii = np.logical_or(transitions["elem1"] == selected_elem, + transitions["elem2"] == selected_elem) + if np.sum(ii) != 0: + selected_transitions = transitions[ii] + redraw=True # force redraw + reset_limits=True + try: + label_rv = selected_model.metadata["manual_rv"] + except KeyError: + label_rv = None + else: + selected_transitions = None + label_rv = None else: selected_transitions = None label_rv = None + # Update figure self.figure.update_spectrum_figure(redraw=redraw, reset_limits=reset_limits, @@ -1022,11 +1129,12 @@ def update_spectrum_figure(self, redraw=False, reset_limits=True): label_rv=label_rv) - def update_fitting_options(self): - try: - selected_model = self._get_selected_model() - except IndexError: - return None + def update_fitting_options(self, selected_model=None): + if selected_model is None: + try: + selected_model = self._get_selected_model() + except IndexError: + return None if selected_model is None: return None if isinstance(selected_model, ProfileFittingModel): @@ -1197,6 +1305,13 @@ def update_edit_fit_window(self): return None else: model.metadata["window"] = window + transitions = self._get_selected_model().transitions + xlim = (transitions["wavelength"][0] - window, + transitions["wavelength"][-1] + window) + self.ax_spectrum.set_xlim(xlim) + self.ax_residual.set_xlim(xlim) + self.figure.reset_zoom_limits() + self.figure.draw() return None def clicked_checkbox_continuum(self): """ The checkbox for modeling the continuum was clicked. """ @@ -1414,6 +1529,7 @@ def clicked_checkbox_upper_limit_2(self): = self.checkbox_upper_limit_2.isChecked() self.measurement_view.update_row(proxy_index.row()) self.summarize_current_table() + self.fit_none(spectral_model) self.refresh_plots() return None def clicked_btn_find_upper_limit(self): @@ -1425,15 +1541,22 @@ def clicked_btn_find_upper_limit(self): except: logger.debug("Invalid sigma for finding limit") return None + # E. Holmbeck changed start_at_current from True to False upper_limit = spectral_model.find_upper_limit(sigma=sigma, start_at_current=True) # Refresh GUI self.measurement_view.update_row(proxy_index.row()) self.summarize_current_table() self.update_fitting_options() + self.fit_none(spectral_model) self.refresh_plots() return None - def clicked_btn_update_abund_table(self): - selected_model = self._get_selected_model() + + def clicked_btn_update_abund_table(self,row=None): + if row is None or not row: + selected_model = self._get_selected_model() + else: + selected_model = self.parent.session.metadata["spectral_models"][row] + if selected_model is None: return None assert isinstance(selected_model, SpectralSynthesisModel), selected_model summary_dict = self.parent.session.summarize_spectral_models(organize_by_element=True) @@ -1496,7 +1619,96 @@ def clicked_export_synthesis(self): spectral_model.export_line_list(linelist_path) logger.info("Exported to {}, {}, {}, {}".format(synth_path, data_path, param_path, linelist_path)) return + + # E. Holmbeck added + def clicked_fit_all_synth(self, skip_message=False): + if not skip_message: + num_models = 0 + for sm in self.parent.session.metadata.get("spectral_models", []): + if isinstance(sm, SpectralSynthesisModel) and sm.is_acceptable==self.acceptables_only: num_models += 1 + + time_estimate = num_models * 10.0 + if time_estimate >= 60: + units = "minutes" + time_estimate /= 60 + else: units = "seconds" + + reply = QtGui.QMessageBox.question(self, "Warning!", + f"Are you sure you want to re-synthesize all {num_models:.0f} acceptable lines?"\ + + f" This will take about {np.around(time_estimate,0):.0f} {units:}, but this may be faster or slower depending on your machine.", + QtGui.QMessageBox.Yes, QtGui.QMessageBox.No) + + if not reply==QtGui.QMessageBox.Yes: + return None + start_time = time.time() + logger.info("Re-fitting all synth lines.") + # TODO: HACKY! + row_count = -1 + for sm in self.parent.session.metadata.get("spectral_models", []): + row_count+=1 + if not isinstance(sm, SpectralSynthesisModel): continue + if sm.is_acceptable!=self.acceptables_only: continue + #self.update_spectrum_figure(redraw=True) + logger.info("Re-fitting {:} at {:.1f}.".format(sm.elements[0], sm.wavelength)) + self.clicked_btn_update_abund_table(row=row_count) + try: + res = sm.fit() + except ValueError as e: + logger.info("Please have at least one abundance for an element, otherwise, I'll skip it!") + print(f"Please have at least one abundance for {sm.elements[0]:}, otherwise, I'll skip it!") + logger.info(e) + continue + except RuntimeError as e: + logger.info("Fitting error") + logger.info(e) + self.acceptables_only = True + return None + self.summarize_current_table() + self.update_fitting_options() + self.refresh_plots(row_count) + self.fit_none(sm) + #self.figure.draw() + + end_time = time.time() + duration = end_time - start_time + minutes,seconds = divmod(duration,60) + time_string = '' + if minutes > 0: + time_string += f"{minutes:.0f} minutes and " + time_string += f"{np.around(seconds,0):.0f} seconds" + + QtGui.QMessageBox.information(self, + "Done!", f"Done! That took {time_string:}.") + + self.acceptables_only = True + return None + + def clicked_fit_rest_synth(self): + num_models = 0 + for sm in self.parent.session.metadata.get("spectral_models", []): + if not sm.is_acceptable and isinstance(sm, SpectralSynthesisModel): num_models += 1 + + time_estimate = num_models * 10.0 + if time_estimate >= 60: + units = "minutes" + time_estimate /= 60 + else: units = "seconds" + + line_or_lines = 's' if num_models>0 else '' + + reply = QtGui.QMessageBox.question(self, "Warning!", + f"You have chosen to synthesize {num_models:.0f} remaining line{line_or_lines:}."\ + + f" This will take about {np.around(time_estimate,0):.0f} {units:}, but this may be faster or slower depending on your machine. Continue?", + QtGui.QMessageBox.Yes, QtGui.QMessageBox.No) + if not reply==QtGui.QMessageBox.Yes: + return None + + logger.info("Re-fitting all synth lines.") + self.acceptables_only = False + self.clicked_fit_all_synth(skip_message=True) + + def refresh_current_model(self): spectral_model, proxy_index, index = self._get_selected_model(True) if spectral_model is None: return None diff --git a/smh/gui/normalization.py b/smh/gui/normalization.py index 4a793ef9..9a94921e 100644 --- a/smh/gui/normalization.py +++ b/smh/gui/normalization.py @@ -277,6 +277,40 @@ def __init__(self, parent): settings_layout.addWidget(self.stitch_btn) + grid_layout = QtGui.QGridLayout() + #settings_grid_layout.addLayout(hbox, 7, 1, 1, 1) + settings_layout.addLayout(grid_layout) + self.line = QtGui.QFrame() + self.line.setFrameShape(QtGui.QFrame.HLine) + self.line.setFrameShadow(QtGui.QFrame.Sunken) + #spacer = QtGui.QSpacerItem(QtGui.QSizePolicy.MinimumExpanding, 10) + #grid_layout.addItem(spacer, 0, 0, 1, 2) + grid_layout.addWidget(self.line, 0, 0, 1, 2, alignment=QtCore.Qt.AlignBottom) + + buttons = { + "Ctrl+Drag": "Zoom", + "Click+Drag": "Add mask", + "Left/Right": "Change order", + "Up/Down": "Scale continuum by 1%", + "a": "Add a point", + "x": "Remove nearest point", + "u": "Undo last mask", + "w": "Remove nearest mask", + "c": "Clear all masks and points", + "r": "Reset zoom limits", + "d": "Do not normalize this order", + "f": "Refit", + } + + for row,(button,description) in enumerate(buttons.items(),1): + label = QtGui.QLabel(self) + label.setText(button) + label.setFont(font) + grid_layout.addWidget(label, row, 0, 1, 1, alignment=QtCore.Qt.AlignCenter) #int fromRow, int fromColumn, int rowSpan, int columnSpan, alignment + label = QtGui.QLabel(self) + label.setText(description) + grid_layout.addWidget(label, row, 1, 1, 1) #int fromRow, int fromColumn, int rowSpan, int columnSpan, alignment + # Add a spacer. settings_layout.addItem(QtGui.QSpacerItem( 40, 20, QtGui.QSizePolicy.Minimum, QtGui.QSizePolicy.Expanding)) @@ -331,14 +365,15 @@ def __init__(self, parent): # Additional point markers. self.ax_order.scatter([], [], facecolor="k", zorder=5, picker=5) - self.ax_order.set_xticklabels([]) + self.ax_order.tick_params(axis='x', labelbottom=False) self.ax_order.set_ylabel("Flux") - self.ax_order_norm = self.norm_plot.figure.add_subplot(gs[1]) + self.ax_order_norm = self.norm_plot.figure.add_subplot(gs[1], sharex=self.ax_order) self.ax_order_norm.axhline(1, linestyle=":", c="#666666", zorder=1) self.ax_order_norm.plot([np.nan], [np.nan], c='k', zorder=2) # TODO: Make (0, 1.2) a default view setting. + self.ax_order_norm.tick_params(axis='x', labelbottom=True) self.ax_order_norm.set_ylim(0, 1.2) self.ax_order_norm.set_yticks([0, 0.5, 1.0]) self.ax_order_norm.set_xlabel(u"Wavelength (Å)") @@ -528,7 +563,7 @@ def figure_key_press(self, event): # undo/remove the last mask - if event.key in ("u", "U"): + if event.key in "uU": if "exclude" in self._cache["input"]: exclude_regions = self._cache["input"]["exclude"] if len(exclude_regions) > 0: @@ -580,14 +615,47 @@ def figure_key_press(self, event): # 'x': clear all added points if event.key in "xX": - for key in ["additional_points"]: - if key in self._cache["input"]: - del self._cache["input"][key] - + if "additional_points" not in self._cache["input"]: + return True + xy = np.array([event.xdata,event.ydata]) + points_xy = self._cache["input"]["additional_points"][:,:2] + if len(points_xy)==0: + return True + diff = np.power(points_xy - xy,2) + diff = np.sum(diff, axis=1) + remove_point = np.where(diff==np.min(diff))[0] + if len(remove_point)==0: + return True + self._cache["input"]["additional_points"] = np.delete(self._cache["input"]["additional_points"], remove_point[0], axis=0) + self.fit_continuum(clobber=True) self.draw_continuum(refresh=False) self.update_continuum_mask(refresh=True) return True + + # remove nearest mask + if event.key in "wW": + if "exclude" in self._cache["input"]: + exclude_regions = self._cache["input"]["exclude"] + if len(exclude_regions)==0: + return True + + x = event.xdata + diff = np.abs(exclude_regions - x) + remove_row = np.where(diff == np.min(diff))[0] + if len(remove_row)>0: + exclude_regions = np.delete(exclude_regions, remove_row[0], axis=0) + ''' + for i,(xlow,xhigh) in enumerate(exclude_regions): + if xlow<=x<=xhigh: + exclude_regions = np.delete(exclude_regions,i,axis=0) + break + #else: + ''' + self._cache["input"]["exclude"] = exclude_regions + self.fit_continuum(clobber=True) + self.draw_continuum(refresh=False) + self.update_continuum_mask(refresh=True) def figure_mouse_press(self, event): """ @@ -790,7 +858,8 @@ def _populate_widgets(self): keys = ("function", "order", "low_sigma_clip", "high_sigma_clip", "knot_spacing", "blue_trim", "red_trim", "max_iterations") self._cache = { - "input": {} + #"input": self.parent.session.metadata['normalization'].get('cache', {}) + "input": {} # Holmbeck: HACKY!! } for key in keys: self._cache["input"][key] \ @@ -887,9 +956,32 @@ def update_rv_applied(self): return None + # E. Holmbeck added this so the masks aren't overwritten in RPA data + def update_rv_applied_keep_masks(self, rv_diff): + """ + Make updates to the view when the radial velocity applied has been + updated. Keep masks. + """ + # Holmbeck: HACKY; this should be shifting the terrestrial lines by the rv_diff + rv_applied = self.parent.session.metadata["rv"].get("rv_applied", 0.0) + # These come pre-loaded + for N in range(len(self.parent.session.input_spectra)): + if 'exclude' in self.parent.session.metadata["normalization"]['normalization_kwargs'][N]: + self.parent.session.metadata["normalization"]['normalization_kwargs'][N]['exclude'] *= 1.0 + rv_diff/c + + # Update the current order fit, and the view. + self.update_order_index() + # May 7 -- added these two lines back in. + #self.update_continuum_mask(refresh=False, v=rv_diff) + self.fit_continuum(clobber=True) + self.draw_order(refresh=False) + self.draw_continuum(refresh=True) + + return None + - def update_continuum_mask(self, refresh=False): + def update_continuum_mask(self, refresh=False, v=0): """ Draw the continuum mask (relevant for all orders). """ @@ -929,13 +1021,9 @@ def update_continuum_mask(self, refresh=False): # will also be in the rest frame. So we don't need to shift the # 'rest_wavelength' mask, but we do need to shift the 'obs_wavelength' # mask - - # Get the applied velocity to shift some masks. - try: - rv_applied = self.parent.session.metadata["rv"]["rv_applied"] - except (AttributeError, KeyError): - rv_applied = 0 + # Get the applied velocity to shift some masks. + rv_applied = self.parent.session.metadata["rv"].get("rv_applied", 0.0) # ----------------------------------------------------------------- # E. Holmbeck added read-in BCV from header try: @@ -948,8 +1036,8 @@ def update_continuum_mask(self, refresh=False): _ =self.parent.session.metadata["normalization"]["normalization_kwargs"] masked_regions = [ #np.array(mask.get("rest_wavelength", [])), - np.array(mask.get("rest_wavelength", [])) * (1.0 - dop_shift/c), - np.array(mask.get("obs_wavelength", [])) * (1.0 - rv_applied/c), + np.array(mask.get("rest_wavelength", [])) * (1.0 + dop_shift/c), + np.array(mask.get("obs_wavelength", [])) * (1.0 + rv_applied/c), np.array(_[self.current_order_index].get("exclude", [])) ] if "pixel" in mask: @@ -1016,13 +1104,13 @@ def update_knot_spacing(self): return None - # ----------------------------------------------------------------- - # E. Holmbeck added these update functions + # ----------------------------------------------------------------- + # E. Holmbeck added these update functions; not sure that this is correct though. def update_blue_trim(self): try: - trim_region = int(self.blue_trim.text()) + trim_region = int(self.blue_trim.text()) except ValueError: - return None + return None if trim_region == 0: return None @@ -1065,9 +1153,9 @@ def update_blue_trim(self): def update_red_trim(self): try: - trim_region = int(self.red_trim.text()) + trim_region = int(self.red_trim.text()) except ValueError: - return None + return None if trim_region == 0: return None @@ -1105,7 +1193,7 @@ def update_red_trim(self): self.update_continuum_mask(refresh=True) return None - # ----------------------------------------------------------------- + # ----------------------------------------------------------------- def update_high_sigma_clip(self): """ Update the high sigma clip value. """ @@ -1204,7 +1292,8 @@ def update_order_index(self, index=None): v = session.metadata["rv"]["rv_applied"] except (AttributeError, KeyError): v = 0 - + + # TODO: Holmbeck note: this was -v/c...which is correct? self.current_order._dispersion *= (1 - v/c) # Update the view if the input settings don't match the settings used @@ -1451,7 +1540,7 @@ def fit_continuum(self, clobber, sender=None): regions = [] for v, masked_regions in mask_kinds: for region in np.array(masked_regions): - start, end = np.array(region) * (1 - v/c) + start, end = np.array(region) * (1 - v/c) # Holmbeck: + or -? if end >= self.current_order.dispersion[0] \ and self.current_order.dispersion[-1] >= start: diff --git a/smh/gui/review.py b/smh/gui/review.py index 1498844c..c76c228e 100644 --- a/smh/gui/review.py +++ b/smh/gui/review.py @@ -107,9 +107,9 @@ def change_measurement_table_species(self, species): # Reset the model (and its views) self.measurement_model.endResetModel() def refresh_plots(self): - self.plot1.update_scatterplot(False) - self.plot2.update_scatterplot(False) - self.plot3.update_scatterplot(False) + self.plot1.update_scatterplot(True) + self.plot2.update_scatterplot(True) + self.plot3.update_scatterplot(True) self.plot1.update_selected_points(True) self.plot2.update_selected_points(True) self.plot3.update_selected_points(True) @@ -119,6 +119,18 @@ def refresh_selected_points(self): self.plot3.update_selected_points(True) def selected_measurement_changed(self): self.refresh_selected_points() + if self.measurement_view is None or self.measurement_model is None: return None + #logger.debug("update_selected_points ({}, {})".format(self, redraw)) + row = self.measurement_view.selectionModel().selectedRows() + # TODO: this doesn't work yet; ideally, update the plots in the Chemical Abundances tab. + if len(row)>0: + row = row[0] + # We're only going to plot the first one anyway + row_index = row.row() + self.parent.chemical_abundances_tab.measurement_view.update_row(row_index) + spectral_model = self.measurement_model.get_models_from_rows([row_index])[0] + self.parent.chemical_abundances_tab.selected_model_changed(selected_model=spectral_model) + return None def _init_summary_table(self): @@ -175,17 +187,23 @@ def _init_scatterplots(self): ] linefit_styles = [{"color":"k","linestyle":"--","zorder":-99},None,None,None] linemean_styles = [{"color":"k","linestyle":":","zorder":-999},None,None,None] - self.plot1 = SMHScatterplot(None, "expot", "abundances", + # E. Holmbeck added error_styles: + error_styles = [{"ms":40,"markerfacecolor":"None","markeredgecolor":"k","ecolor":"k","lw":1}, + {"ms":40,"markerfacecolor":"None","markeredgecolor":"k","ecolor":"c","lw":1}, + {"ms":70,"markerfacecolor":"none","markeredgecolor":"red","ecolor":"red","lw":2}, + {"ms":0,"markerfacecolor":"none","markeredgecolor":"none","ecolor":"none"}, + ] + self.plot1 = SMHScatterplot(None, "expot", "abundances",# eyattr="abundance_uncertainties", tableview=self.measurement_view, - filters=filters, point_styles=point_styles, + filters=filters, point_styles=point_styles, error_styles=error_styles, linefit_styles=linefit_styles,linemean_styles=linemean_styles) - self.plot2 = SMHScatterplot(None, "reduced_equivalent_width", "abundances", + self.plot2 = SMHScatterplot(None, "reduced_equivalent_width", "abundances",# eyattr="abundance_uncertainties", tableview=self.measurement_view, - filters=filters, point_styles=point_styles, + filters=filters, point_styles=point_styles, error_styles=error_styles, linefit_styles=linefit_styles,linemean_styles=linemean_styles) - self.plot3 = SMHScatterplot(None, "wavelength", "abundances", + self.plot3 = SMHScatterplot(None, "wavelength", "abundances",# eyattr="abundance_uncertainties", tableview=self.measurement_view, - filters=filters, point_styles=point_styles, + filters=filters, point_styles=point_styles, error_styles=error_styles, linefit_styles=linefit_styles,linemean_styles=linemean_styles) sp = QtGui.QSizePolicy(QtGui.QSizePolicy.MinimumExpanding, diff --git a/smh/gui/rv.py b/smh/gui/rv.py index c590264d..c8ff92de 100644 --- a/smh/gui/rv.py +++ b/smh/gui/rv.py @@ -669,7 +669,7 @@ def redraw_normalized_order(self, refresh=False): rv_applied = self.parent.session.metadata["rv"]["rv_applied"] except (AttributeError, KeyError): rv_applied = 0 - + self.ax_order_norm.lines[1].set_data([ self._cache["normalized_order"].dispersion * (1 + rv_applied/c), self._cache["normalized_order"].flux, @@ -771,16 +771,25 @@ def correct_radial_velocity(self): """ Correct the radial velocity of the observed spectra. """ - + if "rv_applied" in self.parent.session.metadata["rv"]: + rv_diff = np.float(self.rv_applied.text()) + self.parent.session.metadata["rv"]["rv_applied"] + else: + rv_diff = 0.0 + + # rv_applied is set in this function self.parent.session.rv_correct(self.rv_applied.text()) - - # Redshift the normalized order. + + # Redshift the normalized order; will look for rv_applied. self.redraw_normalized_order(True) - + # Enable and update the normalization tab. self.parent.tabs.setTabEnabled(self.parent.tabs.indexOf(self) + 1, True) - self.parent.normalization_tab.update_rv_applied() + # New function E. Holmbeck added to keep masks; this shifts the masks by the rv_diff; + # Basically, if the saved rv_applied is different from what is in the text box, shift it by that much + #self.parent.normalization_tab.update_rv_applied() + self.parent.normalization_tab.update_rv_applied_keep_masks(rv_diff) + # Enable relevant menu actions. self.parent._action_fit_balmer_lines.setEnabled(True) @@ -1216,7 +1225,7 @@ def redraw_normalized_order(self, refresh=False): rv_applied = self.rv_tab.parent.session.metadata["rv"]["rv_applied"] except (AttributeError, KeyError): rv_applied = 0 - + self.ax_order_norm.lines[1].set_data([ self.rv_tab._cache["normalized_order"].dispersion * (1 + rv_applied/c), self.rv_tab._cache["normalized_order"].flux, diff --git a/smh/gui/stellar_parameters.py b/smh/gui/stellar_parameters.py index f7297376..1e513f49 100644 --- a/smh/gui/stellar_parameters.py +++ b/smh/gui/stellar_parameters.py @@ -53,6 +53,43 @@ _QFONT = QtGui2.QFont("Helvetica Neue", 10) _ROWHEIGHT = 20 + +# E. Holmbeck added toggle. Shamelessly stolen from https://stackoverflow.com/questions/56806987/switch-button-in-pyqt +class MySwitch(QtGui.QPushButton): + def __init__(self, parent = None): + super().__init__(parent) + self.setCheckable(True) + self.setMinimumWidth(66) + self.setMinimumHeight(22) + + def paintEvent(self, event): + label = "Fe I" if not self.isChecked() else "Fe II" + bg_color = QtCore.Qt.white if not self.isChecked() else QtCore.Qt.red + + radius = 10 + width = 25 + center = self.rect().center() + + painter = QtGui2.QPainter(self) + painter.setRenderHint(QtGui2.QPainter.Antialiasing) + painter.translate(center) + #painter.setBrush(QtGui2.QColor(0,0,0)) + painter.setBrush(QtCore.Qt.gray) + + pen = QtGui2.QPen(QtCore.Qt.black) + pen.setWidth(0) + painter.setPen(pen) + + # Draw background rounded box + painter.drawRoundedRect(QtCore.QRect(-2*width, -radius, 3*width, 2*radius), radius, radius) + painter.setBrush(QtGui2.QBrush(bg_color)) + sw_rect = QtCore.QRect(-radius, -radius, width + radius, 2*radius) + if not self.isChecked(): + sw_rect.moveLeft(-2*width) + painter.drawRoundedRect(sw_rect, radius, radius) + painter.drawText(sw_rect, QtCore.Qt.AlignCenter, label) + + class StellarParametersTab(QtGui.QWidget): @@ -189,6 +226,7 @@ def update_stellar_parameter_session(self): "surface_gravity": float(self.edit_logg.text()), "metallicity": float(self.edit_metallicity.text()), "microturbulence": float(self.edit_xi.text()), + #"numax": float(self.edit_numax.text()), "alpha": float(self.edit_alpha.text()) }) return True @@ -202,6 +240,7 @@ def update_stellar_parameter_labels(self): (self.edit_logg, "{0:.2f}", "surface_gravity"), (self.edit_metallicity, "{0:+.2f}", "metallicity"), (self.edit_xi, "{0:.2f}", "microturbulence"), + #(self.edit_numax, "{0:.2f}", "numax"), (self.edit_alpha, "{0:.2f}", "alpha") ] for widget, fmt, key in widget_info: @@ -250,31 +289,31 @@ def update_stellar_parameter_state_table(self): ## Fit lines try: - mchi1, bchi1, med1, eXH1, emchi1, N1 = utils.fit_line(chi1, eps1) + mchi1, bchi1, wmean1, eXH1, emchi1, N1 = utils.fit_line(chi1, eps1) except Exception as e: logger.debug(e) - mchi1, bchi1, med1, eXH1, emchi1, N1 = np.nan, np.nan, np.nan, np.nan, np.nan, len(eps1) + mchi1, bchi1, wmean1, eXH1, emchi1, N1 = np.nan, np.nan, np.nan, np.nan, np.nan, len(eps1) try: - mREW1, bREW1, med1, eXH1, emREW1, N1 = utils.fit_line(REW1, eps1) + mREW1, bREW1, wmean1, eXH1, emREW1, N1 = utils.fit_line(REW1, eps1) except Exception as e: logger.debug(e) - mREW1, bREW1, med1, eXH1, emREW1, N1 = np.nan, np.nan, np.nan, np.nan, np.nan, len(eps1) + mREW1, bREW1, wmean1, eXH1, emREW1, N1 = np.nan, np.nan, np.nan, np.nan, np.nan, len(eps1) try: - mchi2, bchi2, med2, eXH2, emchi2, N2 = utils.fit_line(chi2, eps2) + mchi2, bchi2, wmean2, eXH2, emchi2, N2 = utils.fit_line(chi2, eps2) except Exception as e: logger.debug(e) - mchi2, bchi2, med2, eXH2, emchi2, N2 = np.nan, np.nan, np.nan, np.nan, np.nan, len(eps2) + mchi2, bchi2, wmean2, eXH2, emchi2, N2 = np.nan, np.nan, np.nan, np.nan, np.nan, len(eps2) try: - mREW2, bREW2, med2, eXH2, emREW2, N2 = utils.fit_line(REW2, eps2) + mREW2, bREW2, wmean2, eXH2, emREW2, N2 = utils.fit_line(REW2, eps2) except Exception as e: logger.debug(e) - mREW2, bREW2, med2, eXH2, emREW2, N2 = np.nan, np.nan, np.nan, np.nan, np.nan, len(eps2) + mREW2, bREW2, wmean2, eXH2, emREW2, N2 = np.nan, np.nan, np.nan, np.nan, np.nan, len(eps2) ## Update table - XH1 = med1 - solar_composition(26.0) - XH2 = med2 - solar_composition(26.1) - meanXH1 = np.mean(eps1) - solar_composition(26.0) - meanXH2 = np.mean(eps2) - solar_composition(26.1) + XH1 = wmean1 - solar_composition(26.0) + XH2 = wmean2 - solar_composition(26.1) + meanXH1 = np.median(eps1) - solar_composition(26.0) + meanXH2 = np.median(eps2) - solar_composition(26.1) self.state_fe1_N.setText(u"Fe I ({})".format(N1)) self.state_fe1_XH.setText(u"{:.2f} ± {:.2f} ({:.2f})".format(XH1,eXH1,meanXH1)) self.state_fe1_dAdchi.setText(u"{:.3f} ± {:.3f}".format(mchi1, emchi1)) @@ -285,6 +324,8 @@ def update_stellar_parameter_state_table(self): self.state_fe2_dAdREW.setText(u"{:.3f} ± {:.3f}".format(mREW2, emREW2)) return None def refresh_plots(self): + self.expotfig.sigma = float(self.edit_sigma.text()) + self.rewfig.sigma = float(self.edit_sigma.text()) self.expotfig.update_scatterplot(False) self.rewfig.update_scatterplot(False) self.expotfig.update_selected_points(True) @@ -343,8 +384,8 @@ def solve_feh(self): ## use current state as initial guess logger.info("Setting [alpha/Fe]=0.4 to solve") self.update_stellar_parameter_session() - self.parent.session.optimize_feh(self.params_to_optimize) - self.parent.session.metadata["stellar_parameters"] + self.parent.session.optimize_feh(self.params_to_optimize, use_FeII=self.toggle_feII.isChecked()) + #self.parent.session.metadata["stellar_parameters"] ## refresh everything # E. Holmbeck added 'new_session' again; trying to fix update problem self.new_session_loaded() @@ -353,12 +394,30 @@ def solve_feh(self): def _init_rt_options(self, parent): + # E. Holmbeck: toggle for Fe I vs. Fe II grid_layout = QtGui.QGridLayout() + label = QtGui.QLabel(self) + label.setText("Use lines for parameters") + grid_layout.addWidget(label, 0, 0, 1, 1) #int fromRow, int fromColumn, int rowSpan, int columnSpan, alignment + self.toggle_feII = MySwitch() + self.toggle_feII.setChecked(True) + grid_layout.addWidget(self.toggle_feII, 0, 0, 1, 3)#, alignment=QtCore.Qt.AlignCenter) + self.toggle_feII.clicked.connect(self.toggle_feII.setChecked(False)) + label = QtGui.QLabel(self) + label.setText("Hold?") + grid_layout.addWidget(label, 0, 2, 1, 1) + self.line = QtGui.QFrame() + self.line.setFrameShape(QtGui.QFrame.HLine) + self.line.setFrameShadow(QtGui.QFrame.Sunken) + spacer = QtGui.QSpacerItem(QtGui.QSizePolicy.MinimumExpanding, 30) + grid_layout.addItem(spacer, 0, 0, 1, 3) + grid_layout.addWidget(self.line, 0, 0, 1, 3, alignment=QtCore.Qt.AlignBottom) + # Effective temperature. label = QtGui.QLabel(self) label.setText("Teff") label.setSizePolicy(QtGui.QSizePolicy(QtGui.QSizePolicy.MinimumExpanding, QtGui.QSizePolicy.Minimum)) - grid_layout.addWidget(label, 0, 0, 1, 1) + grid_layout.addWidget(label, 1, 0, 1, 1) self.edit_teff = QtGui.QLineEdit(self) self.edit_teff.setMinimumSize(QtCore.QSize(40, 0)) self.edit_teff.setMaximumSize(QtCore.QSize(50, 16777215)) @@ -367,13 +426,12 @@ def _init_rt_options(self, parent): self.edit_teff.setValidator( QtGui2.QDoubleValidator(3000, 8000, 0, self.edit_teff)) self.edit_teff.textChanged.connect(self._check_lineedit_state) - grid_layout.addWidget(self.edit_teff, 0, 1) + grid_layout.addWidget(self.edit_teff, 1, 1) # E. Holmbeck added checkbox - self.teff_const = QtGui.QCheckBox("Hold constant") + self.teff_const = QtGui.QCheckBox() self.teff_const.setChecked(False) self.teff_const.stateChanged.connect(lambda:self.const_param(self.teff_const,0)) - #grid_layout.addWidget(self.teff_const, 0, 2, -1) - grid_layout.addWidget(self.teff_const, 0, 2) + grid_layout.addWidget(self.teff_const, 1, 2, alignment=QtCore.Qt.AlignCenter) # Surface gravity. @@ -381,7 +439,7 @@ def _init_rt_options(self, parent): label.setText("logg") label.setSizePolicy(QtGui.QSizePolicy(QtGui.QSizePolicy.MinimumExpanding, QtGui.QSizePolicy.Minimum)) - grid_layout.addWidget(label, 1, 0, 1, 1) + grid_layout.addWidget(label, 2, 0, 1, 1) self.edit_logg = QtGui.QLineEdit(self) self.edit_logg.setMinimumSize(QtCore.QSize(40, 0)) self.edit_logg.setMaximumSize(QtCore.QSize(50, 16777215)) @@ -390,20 +448,20 @@ def _init_rt_options(self, parent): QtGui2.QDoubleValidator(-1, 6, 3, self.edit_logg)) self.edit_logg.setSizePolicy(QtGui.QSizePolicy(QtGui.QSizePolicy.MinimumExpanding, QtGui.QSizePolicy.Minimum)) self.edit_logg.textChanged.connect(self._check_lineedit_state) - grid_layout.addWidget(self.edit_logg, 1, 1) + grid_layout.addWidget(self.edit_logg, 2, 1) # E. Holmbeck added checkbox - self.logg_const = QtGui.QCheckBox("Hold constant") + self.logg_const = QtGui.QCheckBox()#"Hold constant") self.logg_const.setChecked(False) self.logg_const.stateChanged.connect(lambda:self.const_param(self.logg_const,2)) #grid_layout.addWidget(self.logg_const, 1, 2, -1) - grid_layout.addWidget(self.logg_const, 1, 2) + grid_layout.addWidget(self.logg_const, 2, 2, alignment=QtCore.Qt.AlignCenter) # Metallicity. label = QtGui.QLabel(self) label.setText("[M/H]") label.setSizePolicy(QtGui.QSizePolicy(QtGui.QSizePolicy.MinimumExpanding, QtGui.QSizePolicy.Minimum)) - grid_layout.addWidget(label, 2, 0, 1, 1) + grid_layout.addWidget(label, 3, 0, 1, 1) self.edit_metallicity = QtGui.QLineEdit(self) self.edit_metallicity.setMinimumSize(QtCore.QSize(40, 0)) self.edit_metallicity.setMaximumSize(QtCore.QSize(50, 16777215)) @@ -412,13 +470,13 @@ def _init_rt_options(self, parent): QtGui2.QDoubleValidator(-5, 1, 3, self.edit_metallicity)) self.edit_metallicity.setSizePolicy(QtGui.QSizePolicy(QtGui.QSizePolicy.MinimumExpanding, QtGui.QSizePolicy.Minimum)) self.edit_metallicity.textChanged.connect(self._check_lineedit_state) - grid_layout.addWidget(self.edit_metallicity, 2, 1) + grid_layout.addWidget(self.edit_metallicity, 3, 1) # E. Holmbeck added checkbox - self.feh_const = QtGui.QCheckBox("Hold constant") + self.feh_const = QtGui.QCheckBox() self.feh_const.setChecked(False) self.feh_const.stateChanged.connect(lambda:self.const_param(self.feh_const,3)) #grid_layout.addWidget(self.feh_const, 2, 2, -1) - grid_layout.addWidget(self.feh_const, 2, 2) + grid_layout.addWidget(self.feh_const, 3, 2, alignment=QtCore.Qt.AlignCenter) # Microturbulence. @@ -426,7 +484,7 @@ def _init_rt_options(self, parent): label.setText("vt") label.setSizePolicy(QtGui.QSizePolicy(QtGui.QSizePolicy.MinimumExpanding, QtGui.QSizePolicy.Minimum)) - grid_layout.addWidget(label, 3, 0, 1, 1) + grid_layout.addWidget(label, 4, 0, 1, 1) self.edit_xi = QtGui.QLineEdit(self) self.edit_xi.setMinimumSize(QtCore.QSize(40, 0)) self.edit_xi.setMaximumSize(QtCore.QSize(50, 16777215)) @@ -434,20 +492,42 @@ def _init_rt_options(self, parent): self.edit_xi.setValidator(QtGui2.QDoubleValidator(0, 5, 3, self.edit_xi)) self.edit_xi.setSizePolicy(QtGui.QSizePolicy(QtGui.QSizePolicy.MinimumExpanding, QtGui.QSizePolicy.Minimum)) self.edit_xi.textChanged.connect(self._check_lineedit_state) - grid_layout.addWidget(self.edit_xi, 3, 1) + grid_layout.addWidget(self.edit_xi, 4, 1) # E. Holmbeck added checkbox - self.vt_const = QtGui.QCheckBox("Hold constant") + self.vt_const = QtGui.QCheckBox() self.vt_const.setChecked(False) self.vt_const.stateChanged.connect(lambda:self.const_param(self.vt_const,1)) #grid_layout.addWidget(self.vt_const, 3, 2, -1) - grid_layout.addWidget(self.vt_const, 3, 2) + grid_layout.addWidget(self.vt_const, 4, 2, alignment=QtCore.Qt.AlignCenter) + + # Nu-max. WIP + ''' + label = QtGui.QLabel(self) + label.setText("nu-max") + label.setSizePolicy(QtGui.QSizePolicy(QtGui.QSizePolicy.MinimumExpanding, QtGui.QSizePolicy.Minimum)) + ''' + ''' + grid_layout.addWidget(label, 4, 0, 1, 1) + self.edit_numax = QtGui.QLineEdit(self) + self.edit_numax.setMinimumSize(QtCore.QSize(40, 0)) + self.edit_numax.setMaximumSize(QtCore.QSize(50, 16777215)) + self.edit_numax.setAlignment(QtCore.Qt.AlignCenter) + self.edit_numax.setValidator(QtGui2.QDoubleValidator(0, 5, 3, self.edit_numax)) + self.edit_numax.setSizePolicy(QtGui.QSizePolicy(QtGui.QSizePolicy.MinimumExpanding, QtGui.QSizePolicy.Minimum)) + self.edit_numax.textChanged.connect(self._check_lineedit_state) + grid_layout.addWidget(self.edit_numax, 4, 1) + self.use_nu = QtGui.QCheckBox("Use numax") + self.use_nu.setChecked(False) + self.use_nu.stateChanged.connect(lambda:self.nu_max(self.use_nu,1)) + grid_layout.addWidget(self.use_nu, 4, 2) + ''' # Alpha-enhancement. label = QtGui.QLabel(self) label.setText("alpha") label.setSizePolicy(QtGui.QSizePolicy(QtGui.QSizePolicy.MinimumExpanding, QtGui.QSizePolicy.Minimum)) - grid_layout.addWidget(label, 4, 0, 1, 1) + grid_layout.addWidget(label, 5, 0, 1, 1) self.edit_alpha = QtGui.QLineEdit(self) self.edit_alpha.setMinimumSize(QtCore.QSize(40, 0)) self.edit_alpha.setMaximumSize(QtCore.QSize(50, 16777215)) @@ -456,14 +536,33 @@ def _init_rt_options(self, parent): #self.edit_alpha.setValidator(QtGui.QDoubleValidator(0, 0.4, 3, self.edit_alpha)) self.edit_alpha.setSizePolicy(QtGui.QSizePolicy(QtGui.QSizePolicy.MinimumExpanding, QtGui.QSizePolicy.Minimum)) self.edit_alpha.textChanged.connect(self._check_lineedit_state) - grid_layout.addWidget(self.edit_alpha, 4, 1) + grid_layout.addWidget(self.edit_alpha, 5, 1) + + + # Sigma plotting + label = QtGui.QLabel(self) + label.setText("Sigma to plot") + label.setSizePolicy(QtGui.QSizePolicy(QtGui.QSizePolicy.MinimumExpanding, QtGui.QSizePolicy.Minimum)) + + grid_layout.addWidget(label, 6, 0, 1, 1) + self.edit_sigma = QtGui.QLineEdit(self) + self.edit_sigma.setMinimumSize(QtCore.QSize(40, 0)) + self.edit_sigma.setMaximumSize(QtCore.QSize(50, 16777215)) + self.edit_sigma.setAlignment(QtCore.Qt.AlignCenter) + self.edit_sigma.setValidator(QtGui2.QDoubleValidator(0, 5, 2, self.edit_sigma)) + self.edit_sigma.setSizePolicy(QtGui.QSizePolicy(QtGui.QSizePolicy.MinimumExpanding, QtGui.QSizePolicy.Minimum)) + self.edit_sigma.textChanged.connect(self._check_lineedit_state) + grid_layout.addWidget(self.edit_sigma, 6, 1) + self.edit_sigma.setText("2.5") self.edit_teff.returnPressed.connect(self.measure_abundances) self.edit_logg.returnPressed.connect(self.measure_abundances) self.edit_metallicity.returnPressed.connect(self.measure_abundances) self.edit_xi.returnPressed.connect(self.measure_abundances) + #self.edit_numax.returnPressed.connect(self.measure_abundances) self.edit_alpha.returnPressed.connect(self.measure_abundances) - + self.edit_sigma.returnPressed.connect(self.refresh_plots) + return grid_layout def _init_rt_buttons(self, parent): # Buttons for solving/measuring. @@ -509,7 +608,7 @@ def create_label(text, row, col, rowspan=1, colspan=1, align=QtCore.Qt.AlignCent return label # Create Header create_label("Species", 0, 0, align=QtCore.Qt.AlignLeft) - create_label(u"〈[X/H]〉", 0, 1) + create_label(u"〈[X/H]〉 weighted average (median)", 0, 1) create_label(u"∂A/∂χ", 0, 2) create_label(u"∂A/∂REW", 0, 3) @@ -588,16 +687,22 @@ def _init_scatterplots(self, parent): linemean_styles = [{"color":"k","linestyle":":","zorder":-999}, {"color":"r","linestyle":":","zorder":-999}, None,None] + # E. Holmbeck added error_styles: + error_styles = [{"ms":40,"markerfacecolor":"None","markeredgecolor":"k","ecolor":"k","lw":1}, + {"ms":40,"markerfacecolor":"None","markeredgecolor":"k","ecolor":"k","lw":1}, + {"ms":70,"markerfacecolor":"none","markeredgecolor":"red","ecolor":"red","lw":3}, + {"ms":70,"markerfacecolor":"none","markeredgecolor":"red","ecolor":"red","lw":3}, + ] self.expotfig = SMHScatterplot(None, "expot", "abundances", tableview=self.measurement_view, - filters=filters, point_styles=point_styles, + filters=filters, point_styles=point_styles, error_styles=error_styles, linefit_styles=linefit_styles,linemean_styles=linemean_styles, - do_not_select_unacceptable=True) + do_not_select_unacceptable=True, sigma=2.5) self.rewfig = SMHScatterplot(None, "reduced_equivalent_width", "abundances", tableview=self.measurement_view, - filters=filters, point_styles=point_styles, + filters=filters, point_styles=point_styles, error_styles=error_styles, linefit_styles=linefit_styles,linemean_styles=linemean_styles, - do_not_select_unacceptable=True) + do_not_select_unacceptable=True, sigma=2.5) sp = QtGui.QSizePolicy(QtGui.QSizePolicy.MinimumExpanding, QtGui.QSizePolicy.MinimumExpanding) @@ -677,6 +782,17 @@ def const_param(self,param_selected,param_index): else: self.params_to_optimize[param_index] = True + # E. Holmbeck added this function; WIP + def nu_max(self,param_selected,param_index): + teff = float(self.edit_teff.text()) + gsol = 4.438 + Tsol = 10.**3.7617 + nusol = 10**3.484 + nustar = float(self.edit_numax.text()) + newlogg = gsol + np.log10(data['numax']*np.sqrt(teff/Tsol)/nusol) + self.edit_logg.setText(newlogg) + self.parent.session.metadata["stellar_parameters"]["surface gravity"] = newlogg + class StellarParameterUncertaintiesDialog(QtGui.QDialog): def __init__(self, session, diff --git a/smh/gui/style_utils.py b/smh/gui/style_utils.py index 6adaa835..3f2603f0 100644 --- a/smh/gui/style_utils.py +++ b/smh/gui/style_utils.py @@ -139,7 +139,11 @@ def relim_axes(axes, percent=20): else: xlim = (x[0] - 1, x[0] + 1) - + + if xlim is not None: + if xlim[0] == xlim[1]: + xlim[0] -= 1 + xlim[1] += 1 if y.size > 1: ylim = [ diff --git a/smh/gui/ui_mainwindow.py b/smh/gui/ui_mainwindow.py index 86a13d2c..19f6c740 100644 --- a/smh/gui/ui_mainwindow.py +++ b/smh/gui/ui_mainwindow.py @@ -47,7 +47,7 @@ def __init__(self, parent=None, session_path=None, spectrum_filenames=None): self.open_session(session_path) self.setObjectName("smh") - self.resize(1200, 600) + self.resize(1200, 900) desktop = QtGui.QApplication.desktop() self.move(desktop.screen().rect().center() - self.rect().center()) @@ -257,6 +257,12 @@ def new_session(self, filenames=None): defaults = yaml.load(fp, yaml.FullLoader) except AttributeError: defaults = yaml.load(fp) + + # Holmbeck: HACKY, but I think ".update" overwrites normalization. + if 'normalization' in defaults: + normalization = defaults.pop('normalization') + self.session.metadata['normalization']['cache'] = normalization + self.session.metadata.update(defaults) # TODO: WE SHOULD REMOVE THIS: THE GUI SHOULD READ FROM .SETTINGS() @@ -591,7 +597,7 @@ def __init_ui__(self): # Create radial velocity tab self.rv_tab = rv.RVTab(self) - self.tabs.addTab(self.rv_tab, "Radial velocity") + self.tabs.addTab(self.rv_tab, "Radial Velocity") # Create normalization tab. self.normalization_tab = normalization.NormalizationTab(self) @@ -600,7 +606,7 @@ def __init_ui__(self): # Create stellar parameters tab. self.stellar_parameters_tab \ = stellar_parameters.StellarParametersTab(self) - self.tabs.addTab(self.stellar_parameters_tab, "Stellar parameters") + self.tabs.addTab(self.stellar_parameters_tab, "Stellar Parameters") # Create chemical abundances tab self.chemical_abundances_tab \ diff --git a/smh/linelists.py b/smh/linelists.py index c73c1a98..bc064bd0 100644 --- a/smh/linelists.py +++ b/smh/linelists.py @@ -628,11 +628,9 @@ def read_moog(cls,filename,moog_columns=False,**kwargs): #if not np.all(np.isnan(ew)): # print("Read {} EWs out of {} lines".format(np.sum(~np.isnan(ew)),len(ew))) colnames = colnames + ['equivalent_width'] -<<<<<<< HEAD + # E. Holmbeck removed the second float thing; why was this here..? dtypes = dtypes + [np.float32] -======= - dtypes = dtypes + [float] ->>>>>>> py38-mpl313 + #dtypes = dtypes + [float] data = data + [ew] return cls(Table(data,names=colnames,dtype=dtypes),moog_columns=moog_columns,**kwargs) diff --git a/smh/optimize_stellar_params.py b/smh/optimize_stellar_params.py index f37240de..92a1223f 100644 --- a/smh/optimize_stellar_params.py +++ b/smh/optimize_stellar_params.py @@ -96,7 +96,7 @@ def feh_optimization(func): def _decorator(request, *args, **kwargs): - params_to_optimize, previously_sampled_points, total_tolerance, individual_tolerances, use_nlte_grid = args + params_to_optimize, previously_sampled_points, total_tolerance, individual_tolerances, use_nlte_grid, use_FeII = args previously_sampled_points = np.array(previously_sampled_points) @@ -547,6 +547,7 @@ def optimize_feh(initial_guess, transitions, params_to_optimize, EWs=None, max_attempts=5, total_tolerance=1e-4, individual_tolerances=None, maxfev=30, use_nlte_grid=None, + use_FeII=False, ): """ Assumes these are all transitions you want to use for stellar parameters @@ -597,7 +598,7 @@ def optimize_feh(initial_guess, transitions, params_to_optimize, EWs=None, else: parameter_ranges["vt"] = (initial_guess[1], initial_guess[1]) - if params_to_optimize[2]: + if params_to_optimize[2]: # and not params_to_optimize[4]: # Holmbeck: I think I was trying to say, "Don't optimize logg independently, are you crazy?" But I don't know what "4" is supposed to have been. parameter_ranges["logg"] = (0, 5), else: parameter_ranges["logg"] = (initial_guess[2], initial_guess[2]) @@ -606,7 +607,7 @@ def optimize_feh(initial_guess, transitions, params_to_optimize, EWs=None, parameter_ranges["[Fe/H]"] = (-5, 0.5) else: parameter_ranges["[Fe/H]"] = (initial_guess[3], initial_guess[3]) - + # Create a mutable copy for the initial guess #solver_guess = [] @@ -625,7 +626,7 @@ def minimisation_function(stellar_parameters, *args): stellar_parameters : [teff, vt, logg, feh] """ - params_to_optimize, all_sampled_points, total_tolerance, individual_tolerances, use_nlte_grid = args + params_to_optimize, all_sampled_points, total_tolerance, individual_tolerances, use_nlte_grid, use_FeII = args # Old way: #teff, vt, logg, feh = [initial_guess[0], stellar_parameters[0], initial_guess[2], stellar_parameters[1]] @@ -647,10 +648,19 @@ def minimisation_function(stellar_parameters, *args): abundances = rt.abundance_cog(photosphere,transitions,twd=twd) transitions["abundance"] = abundances - out = utils.equilibrium_state(transitions[idx_I], - ("expot", "reduced_equivalent_width")) - dAdchi = out[26.0]['expot'][0] - dAdREW = out[26.0]['reduced_equivalent_width'][0] + # E. Holmbeck changed so that dAdREW is w.r.t. FeII + if use_FeII: + out = utils.equilibrium_state(transitions[idx_II], + ("expot", "reduced_equivalent_width")) + species=26.1 + else: + out = utils.equilibrium_state(transitions[idx_I], + ("expot", "reduced_equivalent_width")) + species=26.0 + logger.info("Using species: {:.1f}".format(species)) + + dAdchi = out[species]['expot'][0] + dAdREW = out[species]['reduced_equivalent_width'][0] dFe = np.mean(abundances[idx_I]) - np.mean(abundances[idx_II]) # E. Holmbeck changed dM to be w.r.t. FeII abundances. dM = np.mean(abundances[idx_II]) - (feh + solar_composition("Fe")) @@ -672,8 +682,8 @@ def minimisation_function(stellar_parameters, *args): for i in range(1, 1 + max_attempts): sampled_points = [] args = (params_to_optimize, sampled_points, total_tolerance, individual_tolerances, - use_nlte_grid) - + use_nlte_grid, use_FeII) + try: results = fsolve(minimisation_function, solver_guess, args=args, fprime=utils.approximate_feh_jacobian, col_deriv=1, epsfcn=0, xtol=1e-10, full_output=1, maxfev=maxfev) @@ -702,7 +712,13 @@ def minimisation_function(stellar_parameters, *args): #ier, mesg = results[-len(sampled_points)/2:] ier, mesg = results[-2:] - final_parameters = (results[0]*params_to_optimize) + (initial_guess*(-1*params_to_optimize + 1)) + # E. Holmbeck updated this + final_parameters = [] #(results[0]*params_to_optimize) + (initial_guess*(-1*params_to_optimize + 1)) + iter_results = iter(results[0]) + for pi in range(len(params_to_optimize)): + if params_to_optimize[pi]: final_parameters.append(next(iter_results)) + else: final_parameters.append(initial_guess[pi]) + final_parameters[0] = int(np.round(final_parameters[0])) # Effective temperature final_parameters_result = results[1]["fvec"] num_moog_iterations = len(sampled_points) diff --git a/smh/radiative_transfer/.DS_Store b/smh/radiative_transfer/.DS_Store new file mode 100644 index 00000000..ebd503d6 Binary files /dev/null and b/smh/radiative_transfer/.DS_Store differ diff --git a/smh/radiative_transfer/moog/utils.py b/smh/radiative_transfer/moog/utils.py index 74f43c18..15eaefa1 100644 --- a/smh/radiative_transfer/moog/utils.py +++ b/smh/radiative_transfer/moog/utils.py @@ -55,7 +55,7 @@ def twd_path(twd=None,**kwargs): #kwds.setdefault("dir", "/tmp/") #kwds.setdefault("prefix", "smh-") twd = mkdtemp(**kwds) - if len(twd) > 30: + if len(twd) > 35: logger.warn( "Temporary working directory should be as short as possible to "\ "prevent MOOG(SILENT) from falling over. Current length ({0}): "\ diff --git a/smh/session.py b/smh/session.py index d81c1dc7..8be5d9f8 100644 --- a/smh/session.py +++ b/smh/session.py @@ -550,7 +550,6 @@ def setting(self, key_tree, default_return_value=None): :param key_tree: A tuple containing a tree of dictionary keys. """ - if isinstance(key_tree, string_types): key_tree = [key_tree] @@ -570,7 +569,7 @@ def setting(self, key_tree, default_return_value=None): try: for key in key_tree: default = default[key] - except KeyError: + except (KeyError, TypeError) as e: return default_return_value else: @@ -597,22 +596,21 @@ def update_default_setting(self, key_tree, value): default = yaml.load(fp, yaml.FullLoader) except AttributeError: default = yaml.load(fp) - - branch = defaults + + # Holmbeck:changed "defaults" to "default" + branch = default for key in key_tree[:-1]: branch.setdefault(key, {}) branch = branch[key] branch[key_tree[-1]] = value with open(self._default_settings_path, "w") as fp: - fp.write(yaml.dump(defaults)) + fp.write(yaml.dump(default)) return True - - def _get_overlap_order(self, wavelength_regions, template_spectrum=None): """ Find the order (and order index) that most overlaps with the template @@ -794,7 +792,9 @@ def rv_correct(self, rv): # ----------------------------------------------------------------- # E. Holmbeck: calculate the bcv if it doesn't exist - if "barycentric_correction" in self.metadata["rv"]: + bcv = self.metadata["rv"].get("barycentric_correction", np.nan) + #if "barycentric_correction" in self.metadata["rv"]: + if ~np.isnan(bcv): return names = self._input_spectra_paths @@ -847,6 +847,7 @@ def rv_correct(self, rv): def stitch_and_stack(self, **kwargs): normalized_orders = [] + for i, (spectrum, continuum) \ in enumerate(zip(self.input_spectra, self.metadata["normalization"]["continuum"])): @@ -1660,6 +1661,8 @@ def summarize_spectral_models(self, spectral_models=None, organize_by_element=Fa I cannot imagine why you'd set it to False unless debugging :param what_fe: 1 or 2 depending on Fe I or Fe II + # Holmbeck added: + If None, use I or II depending on the ionization state of the species X. """ what_key_type = "element" if organize_by_element else "species" @@ -1805,10 +1808,13 @@ def export_spectral_model_measurements(self, filepath): ## We'll eventually put in upper limits too. spectral_models = self.metadata.get("spectral_models", []) # Erika added EW sigma to output - linedata = np.zeros((len(spectral_models), 8)) + np.nan + names=["species", "wavelength", "expot", "loggf", "EW", "e_EW", "logeps", "e_logeps"] + dtypes=[(name, "f4") for name in names[:7]] + [(names[7], "U10")] + # Convoluted, but we need to be able to print out upper limits + linedata = np.array(np.zeros(len(spectral_models)) + np.nan, dtype=dtypes) for i,spectral_model in enumerate(spectral_models): # TODO include upper limits - if not spectral_model.is_acceptable or spectral_model.is_upper_limit: continue + if not spectral_model.is_acceptable: continue# or not spectral_model.is_upper_limit: continue if isinstance(spectral_model, SpectralSynthesisModel): assert len(spectral_model.elements) == 1, spectral_model.elements wavelength = spectral_model.wavelength @@ -1816,13 +1822,15 @@ def export_spectral_model_measurements(self, filepath): expot = spectral_model.expot loggf = spectral_model.loggf EW = np.nan - e_EW = 0. + e_EW = np.nan logeps = spectral_model.abundances[0] - try: - logeps_err = spectral_model.metadata["2_sigma_abundance_error"]/2.0 - except: - logeps_err = np.nan - print("exporting synth",wavelength,species) + if spectral_model.is_upper_limit: + logeps_err = "<" + elif "2_sigma_abundance_error" in spectral_model.metadata: + logeps_err = f"{spectral_model.metadata['2_sigma_abundance_error']/2.0:6.3f}" + else: + logeps_err = "nan" + #print("exporting synth",wavelength,species) elif isinstance(spectral_model, ProfileFittingModel): line = spectral_model.transitions[0] wavelength = line['wavelength'] @@ -1835,22 +1843,30 @@ def export_spectral_model_measurements(self, filepath): # Erika added EW sigma to output e_EW = max(1000.*np.abs(spectral_model.metadata["fitted_result"][2]["equivalent_width"][1:])) logeps = spectral_model.abundances[0] - logeps_err = spectral_model.abundance_uncertainties or np.nan + logeps_err = f"{spectral_model.abundance_uncertainties or np.nan:6.3f}" except Exception as e: print(e) EW = np.nan e_EW = np.nan logeps = np.nan - logeps_err = np.nan + logeps_err = "nan" if EW is None: EW = np.nan if logeps is None: logeps = np.nan else: raise NotImplementedError # Erika added EW sigma to output - linedata[i,:] = [species, wavelength, expot, loggf, EW, e_EW, logeps, logeps_err] + try: + linedata[i] = tuple([species, wavelength, expot, loggf, EW, e_EW, logeps, logeps_err]) + except: + import pdb + pdb.set_trace() + + # This was an extra check for bad lines/upper limits #ii_bad = np.logical_or(np.isnan(linedata[:,5]), np.isnan(linedata[:,4])) - ii_bad = np.isnan(linedata[:,5]) - linedata = linedata[~ii_bad,:] + #ii_bad = np.isnan(linedata[:,6]) # Used to be 5. + #linedata = linedata[~ii_bad,:] + ii_bad = np.isnan(linedata['logeps']) + linedata = np.delete(linedata, ii_bad) if len(linedata) == 0: raise RuntimeError("No lines have abundances measured!") @@ -1864,7 +1880,8 @@ def _export_latex_measurement_table(self, filepath, linedata): raise NotImplementedError def _export_ascii_measurement_table(self, filepath, linedata): # Erika added EW sigma to output - names = ["species", "wavelength", "expot", "loggf", "EW", "e_EW", "logeps", "e_logeps"] + #names = ["species", "wavelength", "expot", "loggf", "EW", "e_EW", "logeps", "e_logeps"] + names = linedata.dtype.names tab = astropy.table.Table(linedata, names=names) tab.sort(["species","wavelength","expot"]) tab["wavelength"].format = ".3f" @@ -1873,7 +1890,7 @@ def _export_ascii_measurement_table(self, filepath, linedata): tab["EW"].format = "6.2f" tab["e_EW"].format = "6.2f" tab["logeps"].format = "6.3f" - tab["e_logeps"].format = "6.3f" + #tab["e_logeps"].format = "s" tab.write(filepath, format="ascii.fixed_width_two_line") return True diff --git a/smh/smh_plotting.py b/smh/smh_plotting.py index cbd88ea2..3d2dd55f 100644 --- a/smh/smh_plotting.py +++ b/smh/smh_plotting.py @@ -47,7 +47,7 @@ def make_summary_plot(summary_figure, normalized_spectrum, (0.33999999999999997, 0.82879999999999987, 0.86), (0.37119999999999997, 0.33999999999999997, 0.86), (0.86, 0.33999999999999997, 0.82879999999999987)] - + hls_colors = ['426b81', '7dac9f', 'd18076', '903d6f'] # Load spectra spectra_objects = {} spectra_colors = {} diff --git a/smh/specutils/motions.py b/smh/specutils/motions.py index 7ca0f99f..172204dc 100644 --- a/smh/specutils/motions.py +++ b/smh/specutils/motions.py @@ -400,6 +400,8 @@ def corrections_from_headers(headers): :type headers: A dictionary-like object. """ + # WARNING: VERY HACKY! + alt_obs = headers.get("ALT_OBS", headers.get("SITEALT", None)) lat_obs = headers.get("LAT_OBS", headers.get("SITELAT", None)) long_obs = headers.get("LONG_OBS", headers.get("SITELONG", None)) @@ -428,7 +430,7 @@ def corrections_from_headers(headers): observatory = observatories_dictionary[origin] alt_obs = observatory["elevation"] lat_obs = observatory["latitude"] - + # Get the RA/DEC. ra = headers.get("RA", None) # Assuming degrees dec = headers.get("DEC", None) @@ -446,11 +448,16 @@ def corrections_from_headers(headers): # Try to get the mid-point directly. try: mjd = Time("{0}T{1}".format(headers["UTMID"], headers["UT-MID"])).mjd - - except (IndexError, KeyError): + except KeyError: + try: + mjd = Time("{0}T{1}".format(headers["DATE-OBS"], headers["UT-MID"])).mjd + except KeyError: + mjd = None + + #except (IndexError, KeyError): + if mjd is None: # Try and calculate it from UT-START/UT-DATE keys #raise - try: utdate_key = [_ for _ in ("UT-DATE", "DATE-OBS") if _ in headers][0] if 'T' in headers[utdate_key]: @@ -465,33 +472,28 @@ def corrections_from_headers(headers): headers[utstart_key]), format="isot", scale="utc") except IndexError: - raise KeyError("cannot find all time keys: UTSTART/UTDATE") + raise KeyError("cannot find time keys: UTSTART/UTDATE") try: utend_key = [_ for _ in ("UTEND", "UT-END") if _ in headers][0] + try: + ut_end = Time("{0}T{1}".format(headers[utdate_key].replace(":", "-"), + headers[utend_key]), format="isot", scale="utc") + except: + ut_end = Time("{0}T{1}".format(headers[utdate_key].replace("/", "-"), + headers[utend_key]), format="isot", scale="utc") except IndexError: try: exp_time = headers['EXPTIME'] + ut_end = ut_start + (exp_time)*u.s except: mjd = ut_start.mjd logging.warn( "Calculating celestial corrections based on the UT-START only") - try: - ut_end = Time("{0}T{1}".format(headers[utdate_key].replace(":", "-"), - headers[utend_key]), format="isot", scale="utc") + mjd = (ut_end - ut_start).jd/2 + ut_start.mjd except: - try: - ut_end = Time("{0}T{1}".format(headers[utdate_key].replace("/", "-"), - headers[utend_key]), format="isot", scale="utc") - except: - ut_end = ut_start + (exp_time)*u.s - - # Get the MJD of the mid-point of the observation. - try: mjd - except: - mjd = (ut_end - ut_start).jd/2 + ut_start.mjd # Calculate the correction. # --------------------------------------------------------------------- diff --git a/smh/specutils/spectrum.py b/smh/specutils/spectrum.py index 42c79db2..910959c3 100644 --- a/smh/specutils/spectrum.py +++ b/smh/specutils/spectrum.py @@ -119,7 +119,8 @@ def read(cls, path, **kwargs): cls.read_alex_spectrum, cls.read_ceres, cls.read_multispec, - cls.read_galah, + cls.read_neid, #E. Holmbeck added + cls.read_galah ) failure_exceptions = [] @@ -502,7 +503,53 @@ def read_fits_multispec(cls, path, flux_ext=None, ivar_ext=None, return (dispersion, flux, ivar, metadata) + # E. Holmbeck added NEID capabilities + @classmethod + def read_neid(cls, path, **kwargs): + image = fits.open(path) + metadata = OrderedDict() + if 'SCIWAVE' not in image: + image = image[0] + + for key, value in image['SCIWAVE'].header.items(): + if key in metadata: + metadata[key] += value + else: + metadata[key] = value + metadata["smh_read_path"] = path + + dispersion = -np.ones((metadata['NAXIS2'], metadata['NAXIS1'])) + for order in range(1,metadata['NAXIS2']+1): + crval = metadata[f'CRVAL{order:.0f}'] + cdelt = metadata[f'CDELT{order:.0f}'] + poly_type = metadata[f'PS{order:.0f}_0'] + poly_param_order = metadata[f'PS{order:.0f}_1'] + #pixels = np.arange(*eval(metadata[f'PS{order:.0f}_2'])) + echelle_order = metadata[f'PS{order:.0f}_3'] + leg_range = eval(metadata[f'PS{order:.0f}_4']) + #poly_params = [metadata[f'PV{order:.0f}_{param:.0f}'] for param in range(poly_param_order)] + poly_params = [] + for param in range(poly_param_order): + value = metadata[f'PV{order:.0f}_{param:.0f}'] + if value is None: poly_params.append(np.nan) + else: poly_params.append(value) + + x = np.linspace(leg_range[0], leg_range[1], metadata['NAXIS1']) + if poly_type=='Legendre' or 'legendre' in poly_type: + from scipy.special import legendre + for li in range(poly_param_order): + dispersion[order-1] += poly_params[li]*legendre(li)(x) + + else: print(poly_type) + + flux = np.array(image['SCIFLUX'].data) + var = np.array(image['SCIVAR'].data) + dispersion = np.atleast_2d(dispersion) + flux = np.atleast_2d(flux) + ivar = np.atleast_2d(var**-1) + return (dispersion, flux, ivar, metadata) + @classmethod def read_fits_spectrum1d(cls, path, **kwargs): """ @@ -658,7 +705,7 @@ def write(self, filename, clobber=True, output_verify="warn"): if not filename.endswith('fits'): a = np.array([self.dispersion, self.flux, self.ivar]).T - np.savetxt(filename, a, fmt="%.4f".encode('ascii')) + np.savetxt(filename, a, fmt="%.6f".encode('ascii')) return else: @@ -1639,7 +1686,7 @@ def stitch(spectra, new_dispersion=None, full_output=False): finite = np.isfinite(common_flux * common_ivar) common_flux[~finite] = 0 common_ivar[~finite] = 0 - + numerator = np.sum(common_flux * common_ivar, axis=0) denominator = np.sum(common_ivar, axis=0) flux, ivar = (numerator/denominator, denominator) diff --git a/smh/utils.py b/smh/utils.py index a0995d99..f037c3e8 100644 --- a/smh/utils.py +++ b/smh/utils.py @@ -180,16 +180,67 @@ def equilibrium_state(transitions, columns=("expot", "rew"), group_by="species", def fit_line(x, y, yerr=None): - if yerr is not None: raise NotImplementedError("Does not fit with error bars yet") finite = np.isfinite(x) & np.isfinite(y) if finite.sum()==0: return np.nan, np.nan, np.nan, np.nan, np.nan, 0 + x, y = x[finite], y[finite] + if len(x)==0: + return np.nan, np.nan, np.nan, np.nan, np.nan, 0 + + # E. Holmbeck added + if np.all(yerr) is None or np.all(np.isnan(yerr)): yerr=None + if yerr is not None: + # E. Holmbeck tried it by hand, but failed. + yerr = np.array(yerr[finite], dtype=float) + xbar = np.mean(x) + n = float(len(x)) + + if len(set(x))>1: + ((m,b_bar), pcov) = optimize.curve_fit(lambda x,m,b: m*x + b, x-xbar, y, sigma=yerr) + m_stderr = pcov[0,0] + b = b_bar - m*xbar + else: + m = 0.0 + m_stderr = 0.0 + b = None + + weights = np.power(yerr,-2) + ymean = np.average(y, weights=weights) + if b is None: b = ymean + yvar = np.sqrt(np.average(np.power(y-ymean,2), weights=weights)) + # Added average measurement uncertainty; TODO: CHECK THIS + weighted_uncertainty = np.sqrt(1.0/sum(weights)) + total_yuncertainty = np.sqrt(weighted_uncertainty**2 + yvar**2) + return m, b, ymean, total_yuncertainty, m_stderr, n + ''' + #raise NotImplementedError("Does not fit with error bars yet") + ybar = np.mean(y) + sxx = np.sum(np.power(x - xbar, 2))/n + syy = np.sum(np.power(y - ybar, 2))/n + sxy = np.sum((x-xbar)*(y-ybar))/n + delta = np.array(np.power(yerr,2), dtype=float) + beta1 = syy - delta*sxx + np.sqrt(np.power(syy-delta*sxx,2) + 4*delta*np.power(sxy,2)) + beta1 /= 2.*sxy + beta0 = ybar - beta1*xbar + xstar = x + beta1*(y - beta0 - beta1*x)/(np.power(beta1,2) + delta) + #ystar = beta0 + beta1*xstar + ystar = beta0 + beta1*x # HERE! + m, b_bar, r, p, m_stderr = stats.linregress(x-xbar, ystar) # HERE! + b = b_bar - m*xbar + ''' + xbar = np.mean(x) x = x - xbar - m, b_bar, r, p, m_stderr = stats.linregress(x, y) - b = b_bar - m*xbar - return m, b, np.median(y), np.std(y), m_stderr, len(x) + if len(set(x))>1: + m, b_bar, r, p, m_stderr = stats.linregress(x, y) + b = b_bar - m*xbar + else: + m = 0.0 + m_stderr = 0.0 + b = np.mean(y) + #return m, b, np.median(y), np.std(y), m_stderr, len(x) + return m, b, np.mean(y), np.std(y), m_stderr, len(x) def spectral_model_conflicts(spectral_models, line_list): """