danielandreasen/fasma

A python interface around MOOG, the stellar spectra code

Hi there 👋

I'm a former astronomer with research interests in stellar high precision spectroscopy. Now I'm a bioinformatician with an interest in especially somatic sequencing.

fasma's Issues

Bumping slope dependent

The amount we bump the parameters should dependent on how far away we are from a solution. E.g. at a high EP slope, we can bump quite a lot, but close to a solution, just give small bumps, so we don't diverge from a solution with a bump.

Speed up interpolation

At the moment I think we are interpolating some columns which are not necessary. The two last columns seem to always be 0 and are not used. And then there is a column with vt which is always constant.
Better to add this columns manually in the end instead of wasting time calling griddata.

Remove outliers in EW method

It could be nice to remove outliers one-by-one in the EW method minimization.

At the end of the minimization, check for outliers. If any exists beyond 3 sigma, remove the one furthest away and redo the minimization from last found best parameters. Continue until all outliers are removed.

Better estimate on rejt for ARES

We should find a way to estimate the rejt parameter in the ARES configuration file. This affect the derived metallicity if it is misplaced!

Use multiple line lists with ARES

Until now we have just used the ARES mode with one line list (usually the one for Fe lines and then the EW method). But it could be nice to prepare two line lists, first one for the EW method, and then one for maybe abundances of many elements.
It can be implemented something along these lines:

normal_optical.ares spectrum.fits extra:elements_optical.ares

The name of the first line list should be as it always is, the second line list can be called something like: spectrum_secondary.moog or if we know it is for abundance: spectrum_abund.moog.
This will make the interface better as mentioned in issue #46.

Configuration file for each driver

Instead of having a StarMe.cfg file we should have something like EW.cfg for the EW driver, synth.cfg for the synthesis driver, and abund.cfg for the abundances driver. The drivers were introduced in commit bc0f0ee and a8f5c9b.

Interpolate Teff, logg and [Fe/H] differently

If we have 4 different temperatures, the models for these should be interpolated first with a cubic spline (or something similar), and then the logg and [Fe/H] with a linear interpolation.

PHOENIX atmospheres

How about some PHOENIX atmospheres too. They seem to be the future anyway.

ewDriver

Every iteration I get this warning:

sh: \rm: command not found

We do not use this command in any function for this driver. Do you have any idea of what this is?
e.g.:

sh: \rm: command not found
 i    Teff    logg    [Fe/H]    vt    EPslope    RWslope    |FeI-FeII|
----------------------------------------------------------------------
sh: \rm: command not found
   1  5769    4.47    +0.01    0.99   -0.005     -0.025       0.02
sh: \rm: command not found

Use spectra for ARES and synth not in "spectra" folder

It could be really nice for the user to use spectra for ARES and synthesis which is not in the spectra folder. If we just do some checks on the the filepath, it should be simple to implement.

ARES crashes

Sometimes ARES crashes, and we don't really do anything about it, and just let it crash with a reference to logARES.txt. The typical solution is to remove the last line it tried fit. Maybe this can be done automatically and logged/printed to terminal/GUI.

MOOG version

The output line list so far is for MOOG2013. Is there a specific format for MOOG2014? In that case do you have the function to format to MOOG2014 so I can add it to the other drivers? Probably it will be added in the ARES driver with MOOG version as an option.

Auto fix vt

Let's add an option where vt can be auto fixed in the EW driver. This should probably be after the outlier removal part, and if the slope is positive and vt=9.99 or the slope is negative with vt=0.0, then autofix it, and rerun the minimization routine.

Linelists

We have to think which linelists to provide.
For stellar parameters we have the Sousa 2007 and Tsantaki 2013 but they need to be calibrated for the specific model atmospheres. For that we need to decide which parameters are optimal in the ewfind.par. E.g. I think we should switch from damping 2 (Blackwell) to 1 (Barklem - more recent values).
We should do the same for the list of Adibekyan 2012.

Find best parameters for a given line list

At the moment our minimization routine stops when convergence is met. However, it is likely that this does not return the best possible parameters available for a given line list.
A way around this, is to make a refinement after convergence is found. E.g. look in a small parameter grid around the converged parameters (+/-25K, +/-0.5km/s, +/-0.1 dex). Then return the best set of parameters.

Model atmospheres for ODFNEW

We can not fill in the gaps in the grid here with the previous grid. We need to interpolate our own grid.

ewDriver

Every iteration I get this warning:
sh: \rm: command not found

We do not use this command in any function for this driver. Do you have any idea of what this is?
e.g.:
sh: \rm: command not found

i Teff logg [Fe/H] vt EPslope RWslope |FeI-FeII|

sh: \rm: command not found
1 5769 4.47 +0.01 0.99 -0.005 -0.025 0.02
sh: \rm: command not found
...

Parallel analysis

With large samples it could be really nice to do the analysis parallel. This is true for all the drivers. Just by running on two cores, it should be twice as fast.
Some warning: The output files genereted from ARES, MOOG and MOOGme has to be labeled in order not to be confused.

Clean output on ARES run

At the moment the output from ARES is a bit messy, and not newlines between different runs.
Let's make it easier to read, and also print which spectra and linelist is used for all runs. Is this something for you @MariaTsantaki 😄

Abundances

Actually this could be easy:

0. Find a line list of element
   Let's start with iron!
1. Set teff, logg, metal, vmic and interpolate the model.
2. Run MOOG (abfind.par)
3. Read output
   Return the average abundance.
   Save output file
   Apply Vardan's outlier removal suggestion

Detailed documentation

More detailed documentation is needed.

I am developing it at the moment.

If there is any particular suggestion let me know.

More information in results.csv

A column with which options were set for a given run, would be really helpful. E.g. if some parameters were fixed, and if the logg is corrected from light curve loggs.

Broadening

When vmac, vsini are too high, the points fall outside the interval. New points need to be added

refine option

When is it optimal to use the refine option?

MOOG version greater than 2013

In MOOG 2014 versions and greater there is an extra column in the summary file. This means we read a wrong column, hence the minimization does not work.
We should fix this fast!

Warning when using normal optical linelist for cool stars

For the EW driver, we could give the user a warning, if the effective temperature is low and the normal line list (optical) is used, instead of the one from Tsantaki et al. Since the latter is a subset of the former, we could just check for a few lines we knoe should only be in the "normal" line list.

Normalization

What is the best way to match the observed (already normalized) with the synthetic spectrum?

1. find the median of the max points
2. linear fit

Atomic data calibration

I have a line e.g.:

Wavelength     ele    EP     loggf        EW
5345.588        6.0      8.85    -3.511                                             6.6

I use the recalibration.py and the driver=ewfind to calculate the new loggf values. The final value is -0.846 which corresponds to EWtheoretical = 0.0.
The difference with the measured value is 6.6. Is this good enough to stop iterations?
Also the value depends on the limits of loggf. If I set the limits +- 1.0 I get -2.513 but with +-10.0 I get -0.846.
Both values give EWtheoretical=0, but it may be important for other lines.

PS We should only accept negative loggfs.

One function for read summary.out

One function to read them all, one function to in the darkness bind them.

Bumping parameters

When bumping parameters, we should only bump parameters which are not fixed, or return the previous value.
We could pass an alpha dependent on the fix, if this is 0 (i.e. the parameter is fixed), then don't change this parameter.

Present results from abundance of many stars

Any good ideas to how we can present the results from abundance analysis of more than 1 star.
Currently, we create a csv file (abundances.csv) which is nice, but it is cramped.
A tool like preprocess.py would be nice, but I don't know how we should present the data.

Error on the flux

We should add an array for the error on the flux for synthesis.
This can be related to the SNR.

Check if files exists

This is in general, before any calculations are done, check if the file requested exists. This should be done in the drivers. If the file(s) does not exists, then continue to the next one.

I have had too many times trying to do ionization balance on a fits file now!

abundances.csv is overwritten in GUI mode

In fact, this file is overwritten every time we run abundanceDriver.py. There should be a overwrite parameter, and it should by default be off.
I'm currently working on this... It's difficult.

Reading atmospheric parameters from MOOG

A few times there are problems reading the atmospheric parameters, but we know Teff, logg, [Fe/H], and vt already, so just pass those on where needed, and problem solved!

Interface between ewMode/synthMode and abundanceMode

The abundance mode will most likely be used after the parameters are found. But if we have to write the entire StarMe_abund.cfg from scratch for multiple stars it will be a pain, since the parameters has to be specified (otherwise solar are used).
I don't have a good idea for the interface, but we should find a way to easily write the parameters to the configuration file.

Minimization problems

Here is a link to the SO question asked for the minimization: http://stackoverflow.com/questions/29855849/minimizing-a-multivariable-function-with-scipy-derivative-not-known

fix feh

The fix_feh option does not work. Either cause the value feh is not read properly or there is a problem in the minimization.

Class based minimization

Make the minimization class based as the Readmoog is. It will be easier to maintain, and easier to read. We will be free of all the parameters been passed around, as they will be part of self.

output results extension

We have been adding the extension '.csv' to the results.
Considering that this file extension has, by convention, a meaning, i.e., comma-separated values and that our output file is actually tab separated values we should consider:
a) change the extension of the output file(e.g. to '.dat')
b) change the output file to actually have commas separating the values

Why is this important? Because of the convention, there are some programs which have options to explicitly import '.csv' file(e.g. TOPCAT).

I know that this is not a priority but, given that we have made an effort, in other places, to keep conventions, e.g. coding conventions, we should be coherent and also try to use conventions in our output files.

What do the rest of you think?

Missing grid points in the model atmosphere

There are some missing models in the grid of Kuricz '93 (and probably also others). So instead of returning the path, check it the file exists first. If not then lower/raise the temperature until the file is found.
Be sure to also return the right temperature range.

Documenting functions in a standard way

We should really document our function using some sort of standard (can be our own), where input and output are clear.
Any suggestions on a standard?

def myfunc(x, y=2.0):
    """My function that does this and that.

    Input
    -----
    x : int
      The first variable.
    y : float
      The second variable (default: 2)

    Output
    ------
    r : list
      The result from x and y
    """

Something like that maybe? It will take some time to do, but it is extremely important for future users (and us).

Make GetModels and interpolation work better together

It would be a lot easier if we could do the interpolation by

interpolator(GetModels(teff, logg, feh, vt, atmtype), save=True)
# or
interpolator(GetModels(teff, logg, feh, atmtype), vt=vt, save=True)  # vt outside GetModels
# or
interpolator(teff, logg, feh, vt, atmtype, save=True)  # Let interpolator handle GetModels

Instead of

grid = GetModels(teff, logg, feh, atmtype)
models, nt, nl, nf = grid.getmodels()
model = interpolator(models, teff=(teff, nt), logg=(logg, nl), feh=(feh, nf))
save_model(model, x)

logg from LC not working

The correction for logg we apply doesn't work right. When I switch in on in batch mode, it doesn't work, but sometimes it work in the GUI version, but I can't turn it off.

Gooey free version of MOOGme

Since there are quite a few problems with the installation of Gooey (when people don't want to use something like anaconda), maybe we should make a Gooey free version of MOOGme.py, where it is entirely called from the commandline.

Implementing synthesis

We need to implement synthesis! I think a todo list as in issue #30 is a good way to go.

• Create/get a spectrum for testing.
• Make i xi-square minimization to re-obtain parameters (try only few parameters free at the time).
• Match continuum points between observed and synthetic, i.e. normalization
• Make a MCMC sampling to obtain parameters (same as above).

Many things more. Please edit here, and add more points.

Implement ARES in the GUI

We can implement ARES as well. It is easy to write a wrapper as we do for MOOG. This is connected to issue #7 but this is the main issue thread.

Add to the GUI

Add both a synthetic mode, and an abundance mode (for analyzing different elements).
With the EW mode we will have in total 3 modes.

Grid at edges and on grid points

We can find a model grid to be interpolated, if the values are in the middle of everything. However, if we choose a parameter which is in the models, the grid should be smaller. Also need to handle grid near edges.

Better guess on initial conditions

Since the execution time on our minimization code is heavily dependent on the initial parameters, maybe we can use something like TMcalc to give a better estimate, if one is not provided.