Skip to content
New issue

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

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

Already on GitHub? Sign in to your account

Read joeb VP output file #418

Merged
merged 9 commits into from
Oct 23, 2017
Merged

Read joeb VP output file #418

merged 9 commits into from
Oct 23, 2017

Conversation

profxj
Copy link
Contributor

@profxj profxj commented Oct 9, 2017

As titled.

Returns a list of AbsComponents but @jnburchett or @ntejos
may want more.

Test but no docs.

Copy link
Collaborator

@jnburchett jnburchett left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Our separate codes to accomplish the same task are quite different in some important ways. I've developed my code while dealing with the huge, and at times challenging, CASBaH lists. For one, I split up the step from joebvp to AbsLines and AbsLines to AbsComponents, as I think that may be useful for other applications. Perhaps we should iterate on this a bit? Your code is definitely cleaner, so the right answer is probably in between.

@@ -41,3 +50,67 @@ def abssys_from_json(filename):

# Return
return abs_sys


def read_joebvp(filename, coord, llist=None):
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should add 'specfile' argument. Will explain below.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

added. optional for now.

nflags, nidx = np.unique(vp_data['nflag'], return_index=True)

# Subset by nflag; Build components
for nflag in nflags:
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we shouldn't pick out components by the 'nflag' values. I see the logic and while this will work in many cases, 'nflag' can also take on values 0 (let parameter vary without being tied to other lines) or 1 (fix parameter value). My approach chooses components based on 'zsys' and then 'trans' (which should probably be called 'species', but anyway...). However, this isn't perfect either, because one may already have declared multiple components with the same systemic redshift but different velocities during the fit. The best approach might involve employing the clustering analysis to group different transitions of the same species and redshift by velocity.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

egads, I worry about the clustering.

Let me try zsys + trans.

for nflag in nflags:
mt_lines = np.where(vp_data['nflag'] == nflag)[0]
if len(np.unique(vp_data['vel'][mt_lines])) != 1:
pdb.set_trace()
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is immediately tripped by the CASBaH parameter file I tried loading because we will occasionally let the velocity centroid of a line flop around a little in the fit due to the COS wavelength uncertainty. Another reason the clustering analysis is perhaps useful here.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

gack. am still resisting to the clustering.
as a human, we had every intent of grouping
these lines together as a component, right?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

anyways, have made this an optional check and
defaulted to False

zlim = [vp_data['zsys'][idx] +
vp_data[vkey][idx] * (1 + vp_data['zsys'][idx]) / ckms
for vkey in ['vlim1', 'vlim2']]
absline = AbsLine(vp_data['restwave'][idx] * u.AA, z=z_fit, zlim=zlim, linelist=llist)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You could just declare the limits in velocity space too, right?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think so. Not for a SpectralLine
But we could add that sometime.

absline.attrib['sig_b'] = vp_data['sigbval'][idx]
absline.attrib['z'] = z_fit
absline.attrib['sig_z'] = (1 + vp_data['z_comp'][idx]) * vp_data['sigvel'][idx] / ckms
absline.attrib['specfile'] = vp_data['specfile'][idx]
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What's in the 'specfile' is probably not going to be especially useful, since it records the path to file at the time of the fit. What is useful (and even necessary for some of my analysis) is setting absline.analy['spec']. For this reason, the function should allow the user to specify the spectrum, either an XSpectrum1D object or a path to the file.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

handled

absline.attrib['coord'] = coord
absline.attrib['flag_N'] = 1
absline.attrib['logN'] = vp_data['col'][idx]
absline.attrib['sig_logN'] = vp_data['sigcol'][idx]
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

absline.attrib['N'] and 'sig_N' should also be set, as these are needed for AbsComponent.synthesize_column()

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ok, done.

@profxj profxj requested a review from ntejos October 19, 2017 11:38
Copy link
Contributor

@ntejos ntejos left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please rename method and perhaps moving it to isgm.utils? otherwise move joebvp_from_components() to isgm.io.

Please add more information to docstring and, if possible, also to the .rst documentation.

setattr(abscomp, key, abscomp.attrib[key])
# Errors must be in first line!
assert abscomp.sig_logN > 0.
comps.append(abscomp)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

"reliability" and "comment" columns should also be included.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not sure how to do this. Going to leave for @ntejos



def read_joebvp(filename, coord, llist=None, specfile=None, chk_vel=False):
"""
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Wasn't a version of this already coded in package joebvp?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't know. @jnburchett ?

----------
filename : str
joeB VP filename
coord : SkyCoord
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Explain coord is supposed to be the coordinate of the QSO sightline for the AbsComponents...

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

added

filename : str
joeB VP filename
coord : SkyCoord
llist : LineList, optional
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Explain how llist will be used if given, and what happens otherwise.

joeB VP filename
coord : SkyCoord
llist : LineList, optional
specfile : str, optional
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Explain specfile is the name of the joebvp file.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

added

if chk_vel:
if len(np.unique(vp_data['vel'][mt_lines])) != 1:
pdb.set_trace()
z_fit = vp_data['zsys'][mt_lines[0]] + vp_data['vel'][mt_lines[0]] * (1 + vp_data['zsys'][mt_lines[0]]) / ckms
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

use ltu.z_from_v() instead.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

good call

absline.attrib['z'] = z_fit
absline.attrib['sig_z'] = (1 + vp_data['z_comp'][idx]) * vp_data['sigvel'][idx] / ckms
if specfile is None:
absline.attrib['specfile'] = vp_data['specfile'][idx]
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this something we want? If the joebvp file already has the specfile, why should we allow the user to change it?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@jnburchett had a good reason (see above)

@@ -0,0 +1,4 @@
specfile|restwave|zsys|col|sigcol|bval|sigbval|vel|sigvel|nflag|bflag|vflag|vlim1|vlim2|wobs1|wobs2|pix1|pix2|z_comp|trans
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Optional columns "comment" and "reliability" may exist.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same

@@ -41,3 +51,83 @@ def abssys_from_json(filename):

# Return
return abs_sys


def read_joebvp(filename, coord, llist=None, specfile=None, chk_vel=False):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Name should be something like: components_from_joebvp() (Note that we already have isgm.utils.joebvp_from_components(). Maybe worth using both in the test case.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

renamed a few things as suggested

@coveralls
Copy link

Coverage Status

Coverage increased (+0.1%) to 71.931% when pulling 07fac5c on ingest_joebvp into dd5c4a8 on master.

@ntejos
Copy link
Contributor

ntejos commented Oct 23, 2017

Have added a commit for reliability and comments. Please review. Also, it seems that "z_fit" is the same as "z_comp". If this is the case, would be good to leave only one of them in the code.

@profxj
Copy link
Contributor Author

profxj commented Oct 23, 2017

Merging

@profxj profxj merged commit 7fbc656 into master Oct 23, 2017
@profxj profxj deleted the ingest_joebvp branch October 23, 2017 22:28
@coveralls
Copy link

Coverage Status

Coverage increased (+0.1%) to 71.926% when pulling 5ca6851 on ingest_joebvp into dd5c4a8 on master.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

5 participants