Examples#

Fitting with HD122563#

In this notebook, we use LOTUS to analyse a well know metal poor star HD122563 to give user a complete process when handling with real metal poor stars with high resolution EW observations. Data come from Frebel et al. 2013. For metal-poor star, we use exp_cutoff=2.7, which means that FeI lines with excitation potential smaller than 2.7 ev will not be considered in the following analysis:

[1]:

import os
import lotus_nlte
fp = os.path.dirname(os.path.realpath(lotus_nlte.__file__))
obs_path = fp + "/test/data/HD122563.csv"

[2]:

from lotus_nlte.gcogs.multigcogs import PolyMultiGCOG
star = "HD122563"
stellar_type = "K/giant/very_metal_poor"
exp_cutoff = 2.7
cal = "nlte"
mgcog = PolyMultiGCOG(star=star, stellar_type=stellar_type, exp_cutoff=exp_cutoff,
                      obs_path=obs_path,cal=cal)
mgcog.pipelines()

All optimizations are based on exist interpolated models and you don't need to interpolate them!
Hypersurface of line 3846.80A with ep=3.25ev of element FeI doesn't have enough points for interpolation
Hypersurface of line 4007.27A with ep=2.76ev of element FeI doesn't have enough points for interpolation
Hypersurface of line 4021.87A with ep=2.76ev of element FeI doesn't have enough points for interpolation
Hypersurface of line 4044.61A with ep=2.83ev of element FeI doesn't have enough points for interpolation
Hypersurface of line 4058.22A with ep=3.21ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 4062.44A with ep=2.85ev of element FeI doesn't have enough points for interpolation
Hypersurface of line 4067.98A with ep=3.21ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 4070.77A with ep=3.24ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 4073.76A with ep=3.27ev of element FeI doesn't have enough points for interpolation
Hypersurface of line 4076.63A with ep=3.21ev of element FeI doesn't have enough points for interpolation
Hypersurface of line 4098.18A with ep=3.24ev of element FeI doesn't have enough points for interpolation
Hypersurface of line 4109.06A with ep=3.29ev of element FeI doesn't have enough points for interpolation
Hypersurface of line 4109.80A with ep=2.85ev of element FeI doesn't have enough points for interpolation
Hypersurface of line 4114.44A with ep=2.83ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 4120.21A with ep=2.99ev of element FeI doesn't have enough points for interpolation
Hypersurface of line 4121.80A with ep=2.83ev of element FeI doesn't have enough points for interpolation
Hypersurface of line 4132.90A with ep=2.85ev of element FeI doesn't have enough points for interpolation
Hypersurface of line 4134.68A with ep=2.83ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 4143.41A with ep=3.05ev of element FeI doesn't have enough points for interpolation
Hypersurface of line 4153.90A with ep=3.40ev of element FeI doesn't have enough points for interpolation
Hypersurface of line 4154.50A with ep=2.83ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 4154.81A with ep=3.37ev of element FeI doesn't have enough points for interpolation
Hypersurface of line 4156.80A with ep=2.83ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 4157.78A with ep=3.42ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 4158.79A with ep=3.43ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 4175.64A with ep=2.85ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 4181.76A with ep=2.83ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 4182.38A with ep=3.02ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 4184.89A with ep=2.83ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 4195.33A with ep=3.33ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 4196.21A with ep=3.40ev of element FeI doesn't have enough points for interpolation
Hypersurface of line 4199.10A with ep=3.05ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 4217.55A with ep=3.43ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 4238.81A with ep=3.40ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 4247.43A with ep=3.37ev of element FeI doesn't have enough points for interpolation
Hypersurface of line 4388.41A with ep=3.60ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 4422.57A with ep=2.85ev of element FeI doesn't have enough points for interpolation
Hypersurface of line 4443.19A with ep=2.86ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 4454.38A with ep=2.83ev of element FeI doesn't have enough points for interpolation
Hypersurface of line 4476.02A with ep=2.85ev of element FeI doesn't have enough points for interpolation
Hypersurface of line 4484.22A with ep=3.60ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 4490.08A with ep=3.02ev of element FeI doesn't have enough points for interpolation
Hypersurface of line 4647.43A with ep=2.95ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 4678.85A with ep=3.60ev of element FeI doesn't have enough points for interpolation
Hypersurface of line 4691.41A with ep=2.99ev of element FeI doesn't have enough points for interpolation
Hypersurface of line 4707.27A with ep=3.24ev of element FeI doesn't have enough points for interpolation
Hypersurface of line 4710.28A with ep=3.02ev of element FeI doesn't have enough points for interpolation
Hypersurface of line 4736.77A with ep=3.21ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 4786.81A with ep=3.00ev of element FeI doesn't have enough points for interpolation
Hypersurface of line 4789.65A with ep=3.53ev of element FeI doesn't have enough points for interpolation
Hypersurface of line 4871.32A with ep=2.87ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 4872.14A with ep=2.88ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 4890.76A with ep=2.88ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 4891.49A with ep=2.85ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 4903.31A with ep=2.88ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 4918.99A with ep=2.85ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 4920.50A with ep=2.83ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 4938.81A with ep=2.88ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 4946.39A with ep=3.37ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 4966.09A with ep=3.33ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 5001.87A with ep=3.88ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 5006.12A with ep=2.83ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 5014.94A with ep=3.94ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 5022.24A with ep=3.98ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 5028.13A with ep=3.56ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 5068.77A with ep=2.94ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 5074.75A with ep=4.22ev of element FeI doesn't have enough points for interpolation
Hypersurface of line 5090.77A with ep=4.26ev of element FeI doesn't have enough points for interpolation
Hypersurface of line 5125.12A with ep=4.22ev of element FeI doesn't have enough points for interpolation
Hypersurface of line 5133.69A with ep=4.18ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 5162.27A with ep=4.18ev of element FeI doesn't have enough points for interpolation
Hypersurface of line 5191.45A with ep=3.04ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 5192.34A with ep=3.00ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 5217.39A with ep=3.21ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 5232.94A with ep=2.94ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 5242.49A with ep=3.63ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 5263.31A with ep=3.27ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 5266.56A with ep=3.00ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 5281.79A with ep=3.04ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 5283.62A with ep=3.24ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 5302.30A with ep=3.28ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 5324.18A with ep=3.21ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 5339.93A with ep=3.27ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 5364.87A with ep=4.45ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 5365.40A with ep=3.56ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 5367.47A with ep=4.42ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 5369.96A with ep=4.37ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 5383.37A with ep=4.31ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 5389.48A with ep=4.42ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 5393.17A with ep=3.24ev of element FeI doesn't have enough points for interpolation
Hypersurface of line 5410.91A with ep=4.47ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 5415.20A with ep=4.39ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 5424.07A with ep=4.32ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 5569.62A with ep=3.42ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 5572.84A with ep=3.40ev of element FeI doesn't have enough points for interpolation
Hypersurface of line 5576.09A with ep=3.43ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 5586.76A with ep=3.37ev of element FeI doesn't have enough points for interpolation
Hypersurface of line 5615.64A with ep=3.33ev of element FeI doesn't have enough points for interpolation
Hypersurface of line 5624.54A with ep=3.42ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 5658.82A with ep=3.40ev of element FeI doesn't have enough points for interpolation
Hypersurface of line 5662.52A with ep=4.18ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 5753.12A with ep=4.26ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 6246.32A with ep=3.60ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 6301.50A with ep=3.65ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 6336.84A with ep=3.69ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 6400.00A with ep=3.60ev of element FeI doesn't have enough points for interpolation
Hypersurface of line 6411.65A with ep=3.65ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 6592.91A with ep=2.73ev of element FeI has already existed and passes test of interpolation
Hypersurface of line 4178.86A with ep=2.58ev of element FeII has already existed and passes test of interpolation
Hypersurface of line 4233.17A with ep=2.58ev of element FeII has already existed and passes test of interpolation
Hypersurface of line 4416.82A with ep=2.78ev of element FeII has already existed and passes test of interpolation
Hypersurface of line 4489.19A with ep=2.83ev of element FeII doesn't have enough points for interpolation
Hypersurface of line 4491.41A with ep=2.86ev of element FeII has already existed and passes test of interpolation
Hypersurface of line 4508.28A with ep=2.86ev of element FeII has already existed and passes test of interpolation
Hypersurface of line 4515.34A with ep=2.84ev of element FeII has already existed and passes test of interpolation
Hypersurface of line 4520.22A with ep=2.81ev of element FeII doesn't have enough points for interpolation
Hypersurface of line 4541.52A with ep=2.86ev of element FeII doesn't have enough points for interpolation
Hypersurface of line 4555.89A with ep=2.83ev of element FeII has already existed and passes test of interpolation
Hypersurface of line 4576.34A with ep=2.84ev of element FeII has already existed and passes test of interpolation
Hypersurface of line 4583.84A with ep=2.81ev of element FeII has already existed and passes test of interpolation
Hypersurface of line 4620.52A with ep=2.83ev of element FeII has already existed and passes test of interpolation
Hypersurface of line 4731.44A with ep=2.89ev of element FeII doesn't have enough points for interpolation
Hypersurface of line 4923.93A with ep=2.89ev of element FeII has already existed and passes test of interpolation
Hypersurface of line 4993.35A with ep=2.81ev of element FeII has already existed and passes test of interpolation
Hypersurface of line 5018.45A with ep=2.89ev of element FeII doesn't have enough points for interpolation
Hypersurface of line 5197.58A with ep=3.23ev of element FeII has already existed and passes test of interpolation
Hypersurface of line 5234.63A with ep=3.22ev of element FeII has already existed and passes test of interpolation
Hypersurface of line 5276.00A with ep=3.20ev of element FeII has already existed and passes test of interpolation
Hypersurface of line 5284.08A with ep=2.89ev of element FeII doesn't have enough points for interpolation
Hypersurface of line 5414.07A with ep=3.22ev of element FeII has already existed and passes test of interpolation
Hypersurface of line 5534.83A with ep=3.25ev of element FeII has already existed and passes test of interpolation
Hypersurface of line 6247.55A with ep=3.89ev of element FeII has already existed and passes test of interpolation
Hypersurface of line 6432.68A with ep=2.89ev of element FeII has already existed and passes test of interpolation
Hypersurface of line 6456.38A with ep=3.90ev of element FeII has already existed and passes test of interpolation

Optimization with proper outliers removal#

[3]:

from lotus_nlte.optimize import DiffEvoStellarOptimization
from lotus_nlte.utils import generate_ranges

[4]:

# define boundary
bounds = [tuple(r) for r in generate_ranges(stellar_type)[:2]]
bounds.append((0.5, 3.0))

[5]:

# define how many iterations for your optimization at most
start_attemp = 1
max_attemps = 3

Removing outliers condition (when should we stop the iteration of optimization): the goal of the code is to search the global mimum values and its corresponding stellar parameters. During this process, the abundance derived from certain lines will have large scattering from the rest of abundance. Therefore, people need to perform removing outliers for those scatterings which can lead to better result. In this loop, you can set up the condition that if objective function is smaller than a certain tolerence, then we don’t need to perform any removing and stop the iteration before reaching the max_attemps, this condition can be:

\[\mathcal{F} = \left(\frac{s_{\chi,1}}{\sigma_{\chi,1}}\right)^2 + \left(\frac{s_{REW,1}}{\sigma_{REW,1}}\right)^2 + \left(\frac{\bar{A_1} - \bar{A_2}}{\sigma_{1-2}}\right)^2 \leqslant tol\]

\(s_{\chi,1}\) is the slope between excitation potentials with derived abundance for FeI, \(\sigma_{\chi,1}\) is its uncertainty; \(s_{REW,1}\) is the slope between reduced EWs with derived abundance for FeI, \(\sigma_{REW,1}\) is its uncertainty. \(\bar{A_i}\) is the mean abundance for each species while \(\sigma_{1-2}\) the standard deviation of the difference of the mean for FeI and FeII.

[6]:

if_removing_outliers_tol = 1e-5
while start_attemp <= max_attemps:
    de = DiffEvoStellarOptimization(mgcog, bounds, physicaltol=[1e-5,2e-5,5e-5])
    result = de.optimize(disp=True, popsize=100, recombination=0.3, mutation=(0.8,1.2))
    #print out results when completing eahc iteration
    print(result['ScipyOptimizeResult'].x)
    print("Optimized slope of abund vs chi is {0:.3f}\n\
        Optimized slope of abund vs ew is {1:.3f}\n\
        The difference between FeI mean abund and FeI mean abund is {2:.3f}".format(result["dAdREW"][0], \
                                                                                    result["dAdchi"][0], \
                                                                                    result["deltaFe"][0]))

    abunds = de.abunds

    if (result["ScipyOptimizeResult"]['fun'] <= if_removing_outliers_tol):
        print("{} iteration is done".format(start_attemp))
        print("Converged already!")
        start_attemp = max_attemps + 1
        continue

    #perform sigma clip when the object function doesn't smaller than some pre-defined value, here it is 1e-5
    mgcog.remove_outliers(abunds)

    print("{} iteration is done".format(start_attemp))
    attemp = attemps + 1

differential_evolution step 1: f(x)= 4.98559
differential_evolution step 2: f(x)= 4.29034
differential_evolution step 3: f(x)= 4.22847
differential_evolution step 4: f(x)= 0.786923
differential_evolution step 5: f(x)= 0.786923
differential_evolution step 6: f(x)= 0.786923
differential_evolution step 7: f(x)= 0.786923
differential_evolution step 8: f(x)= 0.786923
differential_evolution step 9: f(x)= 0.786923
differential_evolution step 10: f(x)= 0.786923
differential_evolution step 11: f(x)= 0.583998
differential_evolution step 12: f(x)= 0.583998
differential_evolution step 13: f(x)= 0.583998
differential_evolution step 14: f(x)= 0.235628
differential_evolution step 15: f(x)= 0.235628
differential_evolution step 16: f(x)= 0.235628
differential_evolution step 17: f(x)= 0.195563
differential_evolution step 18: f(x)= 0.195563
differential_evolution step 19: f(x)= 0.195563
differential_evolution step 20: f(x)= 0.195563
differential_evolution step 21: f(x)= 0.195563
differential_evolution step 22: f(x)= 0.195563
differential_evolution step 23: f(x)= 0.195563
differential_evolution step 24: f(x)= 0.195563
differential_evolution step 25: f(x)= 0.13807
differential_evolution step 26: f(x)= 0.13807
differential_evolution step 27: f(x)= 0.13807
differential_evolution step 28: f(x)= 0.13807
differential_evolution step 29: f(x)= 0.13807
differential_evolution step 30: f(x)= 0.13807
differential_evolution step 31: f(x)= 0.13807
differential_evolution step 32: f(x)= 0.13807
differential_evolution step 33: f(x)= 0.13807
differential_evolution step 34: f(x)= 0.13807
differential_evolution step 35: f(x)= 0.13807
differential_evolution step 36: f(x)= 0.13807
differential_evolution step 37: f(x)= 0.13807
differential_evolution step 38: f(x)= 0.13807
differential_evolution step 39: f(x)= 0.13807
differential_evolution step 40: f(x)= 0.13807
differential_evolution step 41: f(x)= 0.13807
differential_evolution step 42: f(x)= 0.13807
differential_evolution step 43: f(x)= 0.13807
differential_evolution step 44: f(x)= 0.13807
differential_evolution step 45: f(x)= 0.13807
differential_evolution step 46: f(x)= 0.13807
differential_evolution step 47: f(x)= 0.13807
differential_evolution step 48: f(x)= 0.136432
differential_evolution step 49: f(x)= 0.136432
differential_evolution step 50: f(x)= 0.108798
differential_evolution step 51: f(x)= 0.108798
differential_evolution step 52: f(x)= 0.108798
differential_evolution step 53: f(x)= 0.108798
differential_evolution step 54: f(x)= 0.108798
differential_evolution step 55: f(x)= 0.0719884
differential_evolution step 56: f(x)= 0.0719884
differential_evolution step 57: f(x)= 0.0719884
differential_evolution step 58: f(x)= 0.0719884
differential_evolution step 59: f(x)= 0.0719884
differential_evolution step 60: f(x)= 0.0719884
differential_evolution step 61: f(x)= 0.0719884
differential_evolution step 62: f(x)= 0.0719884
differential_evolution step 63: f(x)= 0.0719884
differential_evolution step 64: f(x)= 0.0719884
differential_evolution step 65: f(x)= 0.0719884
differential_evolution step 66: f(x)= 0.0719884
differential_evolution step 67: f(x)= 0.0719884
differential_evolution step 68: f(x)= 0.0376823
differential_evolution step 69: f(x)= 0.0376823
differential_evolution step 70: f(x)= 0.0376823
differential_evolution step 71: f(x)= 0.0208826
differential_evolution step 72: f(x)= 0.0208826
differential_evolution step 73: f(x)= 0.0208826
differential_evolution step 74: f(x)= 0.0208826
differential_evolution step 75: f(x)= 0.0208826
differential_evolution step 76: f(x)= 0.0208826
differential_evolution step 77: f(x)= 0.0208826
differential_evolution step 78: f(x)= 0.00791499
differential_evolution step 79: f(x)= 0.00791499
differential_evolution step 80: f(x)= 0.00791499
differential_evolution step 81: f(x)= 0.00791499
differential_evolution step 82: f(x)= 0.00791499
differential_evolution step 83: f(x)= 0.00791499
differential_evolution step 84: f(x)= 0.00791499
differential_evolution step 85: f(x)= 0.00791499
differential_evolution step 86: f(x)= 0.00791499
differential_evolution step 87: f(x)= 0.00249586
differential_evolution step 88: f(x)= 0.00249586
differential_evolution step 89: f(x)= 0.00249586
differential_evolution step 90: f(x)= 0.00249586
differential_evolution step 91: f(x)= 0.00249586
differential_evolution step 92: f(x)= 0.00249586
differential_evolution step 93: f(x)= 0.00249586
differential_evolution step 94: f(x)= 0.00249586
differential_evolution step 95: f(x)= 0.00249586
differential_evolution step 96: f(x)= 0.00249586
differential_evolution step 97: f(x)= 0.00249586
differential_evolution step 98: f(x)= 0.00249586
differential_evolution step 99: f(x)= 0.00249586
differential_evolution step 100: f(x)= 0.00190166
differential_evolution step 101: f(x)= 0.000637633
differential_evolution step 102: f(x)= 0.000637633
differential_evolution step 103: f(x)= 0.000637633
differential_evolution step 104: f(x)= 0.000637633
differential_evolution step 105: f(x)= 0.000637633
differential_evolution step 106: f(x)= 0.000637633
differential_evolution step 107: f(x)= 0.000637633
differential_evolution step 108: f(x)= 0.000637633
differential_evolution step 109: f(x)= 0.000637633
differential_evolution step 110: f(x)= 0.000637633
differential_evolution step 111: f(x)= 0.000637633
differential_evolution step 112: f(x)= 0.000637633
differential_evolution step 113: f(x)= 0.00011929
differential_evolution step 114: f(x)= 0.00011929
differential_evolution step 115: f(x)= 0.00011929
differential_evolution step 116: f(x)= 0.00011929
differential_evolution step 117: f(x)= 0.00011929
differential_evolution step 118: f(x)= 0.00011929
differential_evolution step 119: f(x)= 0.00011929
differential_evolution step 120: f(x)= 0.00011929
differential_evolution step 121: f(x)= 0.00011929
differential_evolution step 122: f(x)= 0.00011929
differential_evolution step 123: f(x)= 0.00011929
differential_evolution step 124: f(x)= 0.00011929
differential_evolution step 125: f(x)= 0.00011929
differential_evolution step 126: f(x)= 0.00011929
differential_evolution step 127: f(x)= 0.00011929
differential_evolution step 128: f(x)= 0.00011929
differential_evolution step 129: f(x)= 0.00011929
differential_evolution step 130: f(x)= 0.00011929
differential_evolution step 131: f(x)= 0.00011929
differential_evolution step 132: f(x)= 0.00011929
differential_evolution step 133: f(x)= 0.000111215
differential_evolution step 134: f(x)= 0.000111215
differential_evolution step 135: f(x)= 0.000111215
differential_evolution step 136: f(x)= 0.000111215
differential_evolution step 137: f(x)= 0.000111215
differential_evolution step 138: f(x)= 0.000111215
differential_evolution step 139: f(x)= 0.000111215
differential_evolution step 140: f(x)= 0.000111215
differential_evolution step 141: f(x)= 0.000111215
differential_evolution step 142: f(x)= 0.000111215
differential_evolution step 143: f(x)= 0.000111215
differential_evolution step 144: f(x)= 0.000111215
differential_evolution step 145: f(x)= 0.000111215
differential_evolution step 146: f(x)= 0.000111215
differential_evolution step 147: f(x)= 0.000103347
differential_evolution step 148: f(x)= 0.000103347
differential_evolution step 149: f(x)= 0.000103347
differential_evolution step 150: f(x)= 0.000103347
differential_evolution step 151: f(x)= 0.000103347
differential_evolution step 152: f(x)= 0.000103347
differential_evolution step 153: f(x)= 0.000103347
differential_evolution step 154: f(x)= 0.000103347
differential_evolution step 155: f(x)= 0.000103347
differential_evolution step 156: f(x)= 0.000103347
differential_evolution step 157: f(x)= 0.000103347
differential_evolution step 158: f(x)= 0.000103347
differential_evolution step 159: f(x)= 0.000103347
differential_evolution step 160: f(x)= 8.30115e-05
differential_evolution step 161: f(x)= 8.30115e-05
differential_evolution step 162: f(x)= 8.30115e-05
differential_evolution step 163: f(x)= 8.30115e-05
differential_evolution step 164: f(x)= 8.30115e-05
differential_evolution step 165: f(x)= 8.30115e-05
differential_evolution step 166: f(x)= 8.30115e-05
differential_evolution step 167: f(x)= 8.30115e-05
differential_evolution step 168: f(x)= 8.30115e-05
differential_evolution step 169: f(x)= 8.30115e-05
differential_evolution step 170: f(x)= 8.30115e-05
differential_evolution step 171: f(x)= 8.30115e-05
differential_evolution step 172: f(x)= 8.30115e-05
differential_evolution step 173: f(x)= 8.30115e-05
differential_evolution step 174: f(x)= 7.51282e-05
differential_evolution step 175: f(x)= 7.51282e-05
differential_evolution step 176: f(x)= 6.55286e-05
differential_evolution step 177: f(x)= 6.55286e-05
differential_evolution step 178: f(x)= 6.55286e-05
differential_evolution step 179: f(x)= 6.55286e-05
differential_evolution step 180: f(x)= 4.07915e-05
differential_evolution step 181: f(x)= 4.07915e-05
differential_evolution step 182: f(x)= 4.07915e-05
differential_evolution step 183: f(x)= 1.66898e-05
differential_evolution step 184: f(x)= 1.66898e-05
differential_evolution step 185: f(x)= 1.66898e-05
differential_evolution step 186: f(x)= 1.66898e-05
differential_evolution step 187: f(x)= 1.66898e-05
differential_evolution step 188: f(x)= 1.66898e-05
differential_evolution step 189: f(x)= 1.66898e-05
differential_evolution step 190: f(x)= 1.66898e-05
differential_evolution step 191: f(x)= 1.66898e-05
differential_evolution step 192: f(x)= 1.66898e-05
differential_evolution step 193: f(x)= 1.66796e-05
differential_evolution step 194: f(x)= 1.66796e-05
differential_evolution step 195: f(x)= 1.66796e-05
differential_evolution step 196: f(x)= 1.66796e-05
differential_evolution step 197: f(x)= 1.66796e-05
differential_evolution step 198: f(x)= 1.66796e-05
differential_evolution step 199: f(x)= 1.66796e-05
differential_evolution step 200: f(x)= 1.66796e-05
differential_evolution step 201: f(x)= 1.66796e-05
differential_evolution step 202: f(x)= 1.66796e-05
differential_evolution step 203: f(x)= 1.66796e-05
differential_evolution step 204: f(x)= 1.66796e-05
differential_evolution step 205: f(x)= 1.00495e-05
differential_evolution step 206: f(x)= 1.00495e-05
differential_evolution step 207: f(x)= 1.00495e-05
differential_evolution step 208: f(x)= 1.00495e-05
differential_evolution step 209: f(x)= 1.00495e-05
differential_evolution step 210: f(x)= 1.11717e-06
[4.60782388e+03 1.18430819e+00 1.87197435e+00]
Optimized slope of abund vs chi is 0.000
        Optimized slope of abund vs ew is 0.000
        The difference between FeI mean abund and FeI mean abund is -0.000
1 iteration is done
Converged already!

Plot the equilibira for the optimal stellar parameters#

[7]:

%matplotlib notebook
from lotus_nlte.plot import plot_optimized_equilibrium
import numpy as np
idx_fei = np.where(np.array(mgcog.obs_ele) =="FeI")
idx_feii = np.where(np.array(mgcog.obs_ele) =="FeII")
REWs = np.log10(1e-3*np.array(mgcog.obs_ew)/np.array(mgcog.obs_wavelength))
chis = np.array(mgcog.obs_ep)
fit_pars = [result['dAdREW'], result['dAdchi']]
abunds = de.abunds

fig1 = plot_optimized_equilibrium("HD122563", result['stellarpars'], fit_pars, REWs1=REWs[idx_fei], REWs2=REWs[idx_feii],
                                chis1=chis[idx_fei], chis2=chis[idx_feii], abunds1=abunds[idx_fei],
                                abunds2=abunds[idx_feii])

Perform slicesampling MCMC to constrain derived stellar parameters#

Priors and their estimated uncertainty comes from last step given by differential evolution

[8]:

import numpy as np
import pymc3 as pm

from lotus_nlte.sampling import slicesampling
trace = slicesampling(de.log_likelihood, mgcog, result['ScipyOptimizeResult'].x, result['ScipyOptimizeResult'].stderrs, de.bounds)
summary = pm.summary(trace, var_names=["Teff", "logg", "vt", 'feh'])
samples = np.vstack([trace[k] for k in ["Teff", "logg", "vt", 'feh']]).T
priors = np.append(result['ScipyOptimizeResult'].x, np.mean(np.array(abunds)[idx_fei]-7.5))
summary

Multiprocess sampling (4 chains in 4 jobs)
CompoundStep
>Slice: [log_f]
>Slice: [vt]
>Slice: [logg]
>Slice: [Teff]
Sampling 4 chains: 100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████| 4800/4800 [15:36<00:00,  5.13draws/s]

[8]:
mean sd mc_error hpd_2.5 hpd_97.5 n_eff Rhat
Teff 4607.901299 57.527494 1.178904 4496.620479 4716.454628 2179.654198 1.001032
logg 1.174598 0.174566 0.002933 0.806337 1.501483 3268.420109 0.999964
vt 1.869406 0.051727 0.001060 1.756691 1.962808 2218.177535 1.000786
feh -2.768153 0.054252 0.001045 -2.869687 -2.662219 2537.737441 1.001092

In corner plot , poterior probablity distributions of each parameter are displayed and blue points show the priors:

[9]:

import corner
fig2 = corner.corner(samples, labels=["$T_{eff}$", "log$\mathit{g}$", "$\\xi_t$", '[Fe/H]'], truths=priors)

Collect abundance at the mean of posterios probability of each stellar parameter of certain line:

[10]:

wavelength = mgcog.obs_wavelength[0]
exp = mgcog.obs_ep[0]
ele = mgcog.obs_ele[0]
ew = mgcog.obs_ew[0]
model = mgcog.models[0]

Teff = summary["mean"]["Teff"]
e_Teff = summary["sd"]["Teff"]
logg = summary["mean"]["logg"]
e_logg = summary["sd"]["logg"]
vt = summary["mean"]["vt"]
e_vt = summary["sd"]["vt"]
abund = de.generate_met(model, Teff, logg, vt, ew)[0] + 7.46#solor iron abundance
print("The abundance of line {0:.2f}A with ep={1:.2f}ev of element {2:s} at Teff={3:.0f}K, logg={4:.2f} and vt={5:.2f}km/s is {6:.2f}".format(wavelength,exp, ele,
                                                           Teff, logg, vt, abund))

The abundance of line 4058.22A with ep=3.21ev of element FeI at Teff=4608K, logg=1.17 and vt=1.87km/s is 4.95

Check if the optimal objective function at stellar parameters is close with element-wise minimum of objective function#

[11]:

best_vals = {
  "Teff": [Teff, e_Teff],
  "logg": [logg, e_logg],
  "Vmic": [vt, e_vt]
}

[12]:

from lotus_nlte.plot import plot_results_brute
de.mgcog = mgcog
grid = de._set_up_meshgrids()
fig3 = plot_results_brute(result=result, grid=grid, best_vals=best_vals)

Looks good!