Source code for model

#!/usr/bin/env python
# coding: utf-8
# $Id: model.py
# Author: Daniel R. Reese <daniel.reese@obspm.fr>
# Copyright (C) Daniel R. Reese and contributors
# Copyright license: GNU GPL v3.0
#
#   AIMS is free software: you can redistribute it and/or modify
#   it under the terms of the GNU General Public License as published by
#   the Free Software Foundation, either version 3 of the License, or
#   (at your option) any later version.
#
#   This program is distributed in the hope that it will be useful,
#   but WITHOUT ANY WARRANTY; without even the implied warranty of
#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#   GNU General Public License for more details.
#
#   You should have received a copy of the GNU General Public License
#   along with AIMS.  If not, see <http://www.gnu.org/licenses/>.
#

"""
A module which contains various classes relevant to the grid of models:

- :py:class:`Model`: a model
- :py:class:`Track`: an evolutionary track
- :py:class:`Model_grid`: a grid of models

These different classes allow the program to store a grid of models and perform
a number of operations, such as:

- retrieving model properties
- interpolate within the grid models
- sort the models within a given evolutionary track
- ...

"""

__docformat__ = 'restructuredtext'

import os
import sys

# AIMS configuration:
if 'AIMS_configure.py' in os.listdir(os.getcwd()):
    sys.path = [os.getcwd(), *sys.path]
import AIMS_configure as config

print(f'AIMS_configure.py path in model.py: {config.__file__}')
# packages from within AIMS:
import constants
import utilities
import aims_fortran

# various packages needed for AIMS
import math
import numpy as np
import random
import matplotlib

if (config.backend is not None):
    matplotlib.use(config.backend)
import matplotlib.pyplot as plt
from scipy.spatial import Delaunay
from bisect import bisect_right

# maximum filename length
# filename_dtype = "U1000"
filename_dtype = "str"

# Strings for age parameters
age_str = "Age"
""" Name of the physical age parameter used in the interpolation """

age_adim_str = "Age_adim"
""" Name of the non-dimensional age parameter used in the interpolation """

# Global parameters:
tol = 1e-6
""" tolerance level for points slightly outside the grid """

eps = 1e-6
""" relative tolerance on parameters used for setting up evolutionary tracks """

log0 = -1e150
""" value which is returned when log(0) is calculated (rather than causing an error) """

ntype = np.int16
""" type used for the n values """

ltype = np.int8  # int8 ranges from -128 to 127
""" type used for the l values """

ftype = np.float64
"""  type used for the frequencies """

gtype = np.float64
""" type used for grid data """

modetype = [('n', ntype), ('l', ltype), ('freq', ftype), ('inertia', ftype)]
""" structure for modes """

# quantities related to user-defined parameters:
user_params_index = {}
""" dictionary which will supply the appropriate index for the user-defined parameters"""

user_params_latex = {}
""" dictionary which will supply the appropriate latex name for the user-defined parameters"""

# integer indices for various global quantities which will be stored in a np array:
nglb = 9 + len(config.user_params)
""" total number of global quantities in a model (see :py:data:`Model.glb`)."""

nlin = 6 + len(config.user_params)
"""
total number of global quantities which are interpolated in a linear way (see
:py:func:`combine_models`).  These quantities are numbered 0:nlin-1
"""

iage_adim = 0
""" index of the parameter corresponding to dimensionless age in the :py:data:`Model.glb` array """

iage = 1
""" index of the parameter corresponding to age in the :py:data:`Model.glb` array """

imass = 2
""" index of the parameter corresponding to mass in the :py:data:`Model.glb` array """

itemperature = 3
""" index of the parameter corresponding to temperature in the :py:data:`Model.glb` array """

iz0 = 4
"""
index of the parameter corresponding to the initial metallicity
the :py:data:`Model.glb` array
"""

ix0 = 5
"""
index of the parameter corresponding to the initial hydrogen content
in the :py:data:`Model.glb` array
"""

ifreq_ref = 6 + len(config.user_params)
"""
index of the parameter corresponding to the reference frequency
(used to non-dimensionalise the pulsation frequencies of the model)
in the :py:data:`Model.glb` array
"""

iradius = 7 + len(config.user_params)
""" index of the parameter corresponding to radius in the :py:data:`Model.glb` array """

iluminosity = 8 + len(config.user_params)
""" index of the parameter corresponding to luminosity in the :py:data:`Model.glb` array """


[docs] def string_to_latex(string, prefix="", postfix=""): """ Return a fancy latex name for an input string. :param string: string that indicates for which parameter we're seeking a latex name :param prefix: optional prefix to add to the string :param postfix: optional postfix to add to the string :type string: string :type prefix: string :type postfix: string :return: a fancy latex name :rtype: string .. note:: This also works for the names of the amplitude parameters for surface corrections. """ if (string.startswith("log_")): return string_to_latex(string[4:], prefix + r"\log_{10}\left(", r"\right)" + postfix) if (string.startswith("ln_")): return string_to_latex(string[3:], prefix + r"\ln\left(", r"\right)" + postfix) if (string.startswith("exp_")): return string_to_latex(string[4:], prefix + r"\exp\left(", r"\right)" + postfix) if (string == "Mass"): return r'Mass, $%sM/M_{\odot}%s$' % (prefix, postfix) if (string == "Radius"): return r'Radius, $%sR/R_{\odot}%s$' % (prefix, postfix) if (string == "Luminosity"): return r'Luminosity, $%sL/L_{\odot}%s$' % (prefix, postfix) if (string == "Z"): return r'Metallicity, $%sZ_0%s$' % (prefix, postfix) if (string == "Y"): return r'Helium content, $%sY_0%s$' % (prefix, postfix) if (string == "X"): return r'Hydrogen content, $%sX_0%s$' % (prefix, postfix) if (string == "Ys"): return r'Helium content, $%sY_s%s$' % (prefix, postfix) if (string == "Yc"): return r'Central helium, $%sY_s%s$' % (prefix, postfix) if (string == "zsx_s"): return r'Ratio, $%s(Z/X)_s%s$' % (prefix, postfix) if (string == "zsx_0"): return r'Ratio, $%s(Z/X)_0%s$' % (prefix, postfix) if (string == "Fe_H"): return r'Iron content, $%s\mathrm{[Fe/H]}%s$' % (prefix, postfix) if (string == "Fe_H0"): return r'Initial iron content, $%s\mathrm{[Fe/H]}%s$' % (prefix, postfix) if (string == "M_H"): return r'Metal content, $%s\mathrm{[M/H]}%s$' % (prefix, postfix) if (string == "M_H0"): return r'Initial metal content, $%s\mathrm{[M/H]}_0%s$' % (prefix, postfix) if (string == "dY_dZ"): return r'Enrichment law, $%s\mathrm{dY/dZ}%s$' % (prefix, postfix) if (string == age_str): return r'Age (in Myrs), $%st%s$' % (prefix, postfix) if (string == age_adim_str): return r'Age parameter, $%s\tau%s$' % (prefix, postfix) if (string == "Teff"): return r'$%sT_{\mathrm{eff}}%s$ (in K)' % (prefix, postfix) if (string == "Dnu"): return r'Large separation (in $\mu$Hz), $%s\Delta \nu%s$' % (prefix, postfix) if (string == "numax"): return r'$%s\nu_{\mathrm{max}}%s$ (in $\mu$Hz)' % (prefix, postfix) if (string == "Rho"): return r'Density (in g.cm$^{-3}$), $%s\rho%s$' % (prefix, postfix) if (string == "g"): return r'Surface gravity (in cm.s$^{-2}$), $%sg%s$' % (prefix, postfix) if (string == "Theta_BCZ"): return r'BCZ normalised acoustic depth, $%s\tau_{\mathrm{BCZ}}/2\tau%s$' % ( prefix, postfix) if (string == "Theta_He"): return r'HeII normalised acoustic depth, $%s\tau_{\mathrm{He}}/2\tau%s$' % ( prefix, postfix) if (string == "A_surf"): return r'$%sA_{\mathrm{surf}}%s$' % (prefix, postfix) if (string == "A3_surf"): return r'$%sA_3^{\mathrm{surf}}%s$' % (prefix, postfix) if (string == "Am1_surf"): return r'$%sA_{-1}^{\mathrm{surf}}%s$' % (prefix, postfix) if (string == "alpha_surf"): return r'$%s\alpha_{\mathrm{surf}}%s$' % (prefix, postfix) if (string == "b_Kjeldsen2008"): return r'$%sb_{\mathrm{Kjeldsen\,et\,al.\,(2008)}}%s$' % (prefix, postfix) if (string == "beta_Sonoi2015"): return r'$%s\beta_{\mathrm{Sonoi\,et\,al.\,(2015)}}%s$' % (prefix, postfix) # raises a KeyError if `string` isn't a valid key return user_params_latex[string] % (prefix, postfix)
[docs] def get_surface_parameter_names(surface_option): """ Return the relevant parameter names for a given surface correction option. :param surface_option: specifies the type of surface correction. :type surface_option: string :return: names for the surface parameters :rtype: tuple of strings """ if (surface_option is None): return () if (surface_option == "Kjeldsen2008"): return ("A_surf",) if (surface_option == "Kjeldsen2008_scaling"): return ("A_surf",) if (surface_option == "Kjeldsen2008_2"): return ("A_surf", "b_Kjeldsen2008") if (surface_option == "Ball2014"): return ("A3_surf",) if (surface_option == "Ball2014_2"): return ("A3_surf", "Am1_surf") if (surface_option == "Sonoi2015"): return ("alpha_surf",) if (surface_option == "Sonoi2015_scaling"): return ("alpha_surf",) if (surface_option == "Sonoi2015_2"): return ("alpha_surf", "beta_Sonoi2015") raise ValueError("Unknown surface correction: " + surface_option)
[docs] class Model: """A class which contains a stellar model, including classical and seismic information."""
[docs] def string_to_param(self, string): """ Return a parameter for an input string. :param string: string that indicates which parameter we're seeking :type string: string :return: the value of the parameter :rtype: float """ if (string.startswith("log_")): return math.log10(self.string_to_param(string[4:])) if (string.startswith("ln_")): return math.log(self.string_to_param(string[3:])) if (string.startswith("exp_")): return math.exp(self.string_to_param(string[4:])) if (string == "Mass"): return self.glb[imass] / constants.solar_mass if (string == "Radius"): return self.glb[iradius] / constants.solar_radius if (string == "Luminosity"): return self.glb[iluminosity] / constants.solar_luminosity if (string == "Z"): return self.glb[iz0] if (string == "Y"): return 1.0 - self.glb[iz0] - self.glb[ix0] if (string == "X"): return self.glb[ix0] if (string == "Ys"): return 1.0 - self.glb[user_params_index["Zs"]] - self.glb[user_params_index["Xs"]] if (string == "Yc"): return 1.0 - self.glb[user_params_index["Zc"]] - self.glb[user_params_index["Xc"]] if (string == "zsx_s"): return self.zsx_s if (string == "zsx_0"): return self.zsx_0 if (string == "Fe_H"): return self.FeH if (string == "Fe_H0"): return self.FeH0 if (string == "M_H"): return self.MH if (string == "M_H0"): return self.MH0 if (string == "dY_dZ"): return (1.0 - self.glb[iz0] - self.glb[ix0] - constants.Yp) / self.glb[iz0] if (string == age_str): return self.glb[iage] if (string == age_adim_str): return self.glb[iage_adim] if (string == "Teff"): return self.glb[itemperature] if (string == "Dnu"): return self.find_large_separation() * self.glb[ifreq_ref] if (string == "numax"): return self.numax if (string == "Rho"): return 3.0 * self.glb[imass] / (4.0 * math.pi * self.glb[iradius] ** 3) if (string == "g"): return constants.G * self.glb[imass] / self.glb[iradius] ** 2 if (string == "Theta_BCZ"): return self.glb[user_params_index["tau_BCZ"]] / ( 2.0 * self.glb[user_params_index["tau"]]) if (string == "Theta_He"): return self.glb[user_params_index["tau_He"]] / ( 2.0 * self.glb[user_params_index["tau"]]) if (string == "beta_Sonoi2015"): return self.beta_Sonoi2015 if (string == "b_Kjeldsen2008"): return self.b_Kjeldsen2008 # raises a KeyError if `string` isn't a valid key return self.glb[user_params_index[string]]
def __init__(self, _glb, _name=None, _modes=None): """ :param _glb: 1D array of global parameters for this model. Its dimension should be greater or equal to :py:data:`nglb` :param _name: name of the model (typically the second part of its path) :param _modes: list of modes in the form of tuples (n,l,freq,inertia) which will be appended to the set of modes in the model. :type _glb: np.array :type _name: string :type _modes: list of (int, int, float, float) """ # check inputs assert (_glb[imass] >= 0.0), "A star cannot have a negative mass!" assert (_glb[iradius] >= 0.0), "A star cannot have a negative radius!" assert (_glb[iluminosity] >= 0.0), "A star cannot have a negative luminosity!" assert (_glb[iz0] >= 0.0), "A star cannot have a negative metallicity!" assert (_glb[ix0] >= 0.0), "A star cannot have a negative hydrogen abundance!" assert (_glb[iage] >= 0.0), "A star cannot have a negative age!" assert (_glb[itemperature] >= 0.0), "A star cannot have a negative temperature!" self.name = _name """Name of the model, typically the second part of its path""" self.glb = _glb """Array which will contain various global quantities""" self.glb[ifreq_ref] = 5e5 * math.sqrt(constants.G * self.glb[imass] / self.glb[iradius] ** 3) / math.pi """Characteristic frequency of the model in :math:`\\mathrm{cyclic \\, \\mu Hz}`""" """array containing the modes (n, l, freq, inertia)""" if _modes is not None: self.modes = np.array(_modes, dtype=modetype) else: self.modes = np.empty([0], dtype=modetype) def __del__(self): """ Clears a model instance. """ del self.modes del self.glb del self.name
[docs] def read_file(self, filename): """ Read in a set of modes from a file. This method will either call :py:meth:`read_file_CLES`, :py:meth:`read_file_MESA`, or :py:meth:`read_file_agsm` according to the value of the ``mode_format`` variable in ``AIMS_configure.py``. :param filename: name of the file with the modes. The format of this file is decided by the ``mode_format`` variable in ``AIMS_configure.py``. :type filename: string :return: ``True`` if at least one frequency has been discarded (see note below). :rtype: boolean .. note:: At this stage the frequencies should be expressed in :math:`\\mu\\mathrm{Hz}`. They will be non-dimensionalised in :py:func:`read_model_list_standard`. """ if (config.mode_format == "simple"): # for retrocompatibility return self.read_file_CLES(filename) if (config.mode_format == "CLES"): return self.read_file_CLES(filename) if (config.mode_format == "CLES_Mod"): return self.read_file_CLES_Mod(filename) if (config.mode_format == "MESA"): return self.read_file_MESA(filename) if (config.mode_format == "agsm"): return self.read_file_agsm(filename) if (config.mode_format == "PLATO"): return self.read_file_PLATO(filename) raise ValueError("Unrecognised format \"" + config.mode_format + "\".\n" + "Please choose another value for mode_format in AIMS_configure.py")
[docs] def read_file_CLES(self, filename): """ Read in a set of modes from a file. This uses the "simple" or "CLES" format as specified in the ``mode_format`` variable in ``AIMS_configure.py``. :param filename: name of the file with the modes. The file should contain a one-line header followed by five columns which correspond to l, n, frequency, unused, inertia. :type filename: string :return: ``True`` if at least one frequency has been discarded (see note below). :rtype: boolean .. note:: - The fourth column is discarded. - Frequencies above `config.cutoff` * :math:`\\nu_{\\mathrm{cut-off}}` are discarded. """ freqlim = config.cutoff * self.cutoff exceed_freqlim = False with open(filename) as freqfile: freqfile.readline() # skip head mode_temp = [] for line in freqfile: line = line.strip() columns = line.split() n = int(columns[1]) freq = utilities.to_float(columns[2]) # remove frequencies above AIMS_configure.cutoff*nu_{cut-off} if (freq > freqlim): exceed_freqlim = True continue if (config.npositive and (n < 0)): # remove g-modes if need be continue mode_temp.append((n, int(columns[0]), freq, utilities.to_float(columns[4]))) self.modes = np.array(mode_temp, dtype=modetype) return exceed_freqlim
[docs] def read_file_CLES_Mod(self, filename): """ Read in a set of modes from a file. This uses the "CLES_Mod" format as specified in the ``mode_format`` variable in ``AIMS_configure.py``. :param filename: name of the file with the modes. The file may contain a header (the header lines start with "#") followed by four columns which correspond to l, n, unused, frequency. :type filename: string :return: ``True`` if at least one frequency has been discarded (see note below). :rtype: boolean .. note:: - The third column is discarded. - Frequencies above `config.cutoff` * :math:`\\nu_{\\mathrm{cut-off}}` are discarded. """ freqlim = config.cutoff * self.cutoff exceed_freqlim = False with open(filename) as freqfile: mode_temp = [] for line in freqfile: if line[0] == '#': # Skip comments continue line = line.strip() columns = line.split() n = int(columns[1]) freq = utilities.to_float(columns[3]) * 10 ** 6 # freq tu muHz # remove frequencies above AIMS_configure.cutoff*nu_{cut-off} if (freq > freqlim): exceed_freqlim = True continue if (config.npositive and (n < 0)): # remove g-modes if need be continue mode_temp.append((n, int(columns[0]), freq, 1.0)) # No inertia in file self.modes = np.array(mode_temp, dtype=modetype) return exceed_freqlim
[docs] def read_file_MESA(self, filename): """ Read in a set of modes from a file. This uses the GYRE (i.e. MESA) format as specified in the ``mode_format`` variable in ``AIMS_configure.py``. :param filename: name of the file with the modes. The file should contain a seven-line header followed by various columns which contain l, n, frequency, and inertia for each pulsation mode. :type filename: string :return: ``True`` if at least one frequency has been discarded (see note below). :rtype: boolean .. note:: - Frequencies above `config.cutoff` * :math:`\\nu_{\\mathrm{cut-off}}` are discarded. """ freqlim = config.cutoff * self.cutoff exceed_freqlim = False with open(filename) as freqfile: for i in range(6): freqfile.readline() mode_temp = [] for line in freqfile: line = line.strip() columns = line.split() n = int(columns[1]) # n_pg # n = int(columns[2]) # n_p freq = utilities.to_float(columns[4]) # remove frequencies above AIMS_configure.cutoff*nu_{cut-off} if (freq > freqlim): exceed_freqlim = True continue if (config.npositive and (n < 0)): # remove g-modes if need be continue mode_temp.append((n, int(columns[0]), freq, utilities.to_float(columns[7]))) self.modes = np.array(mode_temp, dtype=modetype) return exceed_freqlim
[docs] def read_file_agsm(self, filename): """ Read in a set of modes from a file. This uses the "agsm" format as specified in the ``mode_format`` variable in ``AIMS_configure.py``. :param filename: name of the file with the modes. This file is a binary fortran "agsm" file produced by the ADIPLS code. See instructions to the ADIPLS code for a description of this format. :type filename: string :return: ``True`` if at least one frequency has been discarded (see note below). :rtype: boolean .. note:: - Frequencies above `config.cutoff` * :math:`\\nu_{\\mathrm{cut-off}}` are discarded. """ narr, larr, farr, iarr, nn, exceed_freqlim = \ aims_fortran.read_file_agsm(filename, config.npositive, config.agsm_cutoff, \ config.cutoff * self.cutoff) self.modes = np.array(list(zip(narr[0:nn], larr[0:nn], farr[0:nn], iarr[0:nn])), dtype=modetype) return exceed_freqlim
[docs] def read_file_PLATO(self, filename): """ Read in a set of modes from a file. This uses the "PLATO" format as specified in the ``mode_format`` variable in ``AIMS_configure.py``. :param filename: name of the file with the modes. The file may contain a header (the header lines start with "#") followed by four columns which correspond to l, n, frequency, inertia. :type filename: string :return: ``False`` :rtype: boolean .. note:: - No mode selection is carried out, hence the return value - np.loadtxt is used to gain speed (the grid is quite large) """ self.modes = np.loadtxt(filename, dtype=modetype, usecols=(1, 0, 2, 3)) return False
[docs] def read_modes_Aldo(self, norder, farray): """ Read in a set of modes from a table. This uses the Aldo format as specified in the ``mode_format`` variable in ``AIMS_configure.py``. :param norder: table which gives the starting n values and number of modes per l value :type norder: 2D numpy int array :param farray: table with mode frequencies and inertias as follows: frequency1, inertia1, frequency2, inertia2 ... :type farray: 1D numpy float array :return: ``True`` if at least one frequency has been discarded (see note below). :rtype: boolean .. note:: - Frequencies above `config.cutoff` * :math:`\\nu_{\\mathrm{cut-off}}` are discarded. """ freqlim = config.cutoff * self.cutoff exceed_freqlim = False mode_temp = [] ii = 0 for l in range(4): for i in range(norder[l, 1]): n = i + norder[l, 0] freq = farray[ii] # remove zero freuqency modes (these modes are not calculated): if (abs(freq) < 1e-10): ii += 2 continue # remove frequencies above AIMS_configure.cutoff*nu_{cut-off} if (freq > freqlim): exceed_freqlim = True ii += 2 continue # remove gravity modes if need be if (config.npositive and (n < 0)): ii += 2 continue # we can add this mode: mode_temp.append((n, l, freq, farray[ii + 1])) self.modes = np.array(mode_temp, dtype=modetype) ii += 2 return exceed_freqlim
[docs] def write_file_simple(self, filename): """ Write a set of modes into a file using the "simple" format as described in :py:meth:`read_file_simple`. :param filename: name of the file where the modes should be written. :type filename: string .. note:: - The output frequencies are expressed in :math:`\\mathrm{\\mu Hz}` """ with open(filename, "w") as output: # write header output.write("# %1s %3s %22s %6s %22s\n" % ("l", "n", "nu_theo (muHz)", "unused", "Inertia")) for i in range(self.modes.shape[0]): output.write(" %1d %3d %22.15e 0.0 %22.15e\n" % ( self.modes["l"][i], self.modes["n"][i], self.modes["freq"][i] * self.glb[ifreq_ref], self.modes["inertia"][i]))
[docs] def append_modes(self, modes): """ Append a list of modes to the model. :param modes: list of modes which are in the form of tuples: (n,l,freq,inertia). :type modes: list of (int, int, float, float) """ self.modes = np.append(self.modes, np.array(modes, dtype=modetype))
[docs] def sort_modes(self): """ Sort the modes by l, then n, then freq. """ # sorts by l, then n, then freq ind = np.lexsort((self.modes['freq'], self.modes['n'], self.modes['l'])) self.modes = np.array([self.modes[i] for i in ind], dtype=modetype)
[docs] def remove_duplicate_modes(self): """ Remove duplicate modes. Modes are considered to be duplicate if they have the same l and n values (regardless of frequency). :return: ``True`` if at least one mode has been removed. :rtype: boolean .. warning:: This method assumes the modes are sorted. """ ind = [] for i in range(len(self.modes) - 1, 1, -1): if ((self.modes['n'][i] == self.modes['n'][i - 1]) and \ (self.modes['l'][i] == self.modes['l'][i - 1])): ind.append(i) self.modes = np.delete(self.modes, ind) return (len(ind) == 0)
[docs] def get_age(self): """ Return age of stellar model. This is useful for sorting purposes. :return: the age of the model :rtype: float """ return self.glb[iage]
[docs] def get_freq(self, surface_option=None, a=[]): """ Obtain model frequencies, with optional frequency corrections. :param surface_option: specifies the type of surface correction. Options include: - ``None``: no corrections are applied - ``"Kjeldsen2008"``: apply a correction based on Kjeldsen et al. (2008) - ``"Kjeldsen2008_scaling"``: apply a correction based on Kjeldsen et al. (2008). The exponent is based on a scaling relation from Sonoi et al. (2015). - ``"Kjeldsen2008_2"``: apply a correction based on Kjeldsen et al. (2008). The exponent is a free parameter. - ``"Ball2014"``: apply a one-term correction based on Ball and Gizon (2014) - ``"Ball2014_2"``: apply a two-term correction based on Ball and Gizon (2014) - ``"Sonoi2015"``: apply a correction based on Sonoi et al. (2015) - ``"Sonoi2015_scaling"``: apply a correction based on Sonoi et al. (2015) The exponent is based on a scaling relation from Sonoi et al. (2015). - ``"Sonoi2015_2"``: apply a correction based on Sonoi et al. (2015) The exponent is a free parameter. :param a: amplitude parameters which intervene in the surface correction :type surface_option: string :type a: array-like :return: models frequencies (including surface corrections) :rtype: np.array .. note:: If surface_option==None or a==[], the original frequencies are returned (hence modifying them modifies the :py:class:`Model` object). """ if (surface_option is None) or (len(a) == 0): return self.modes['freq'] return self.modes['freq'] + self.get_surface_correction(surface_option, a)
[docs] def get_surface_correction(self, surface_option, a): """ Obtain corrections on model frequencies (these corrections should be *added* to the *theorectical* frequencies). :param surface_option: specifies the type of surface correction. Options include: - ``None``: no corrections are applied - ``"Kjeldsen2008"``: apply a correction based on Kjeldsen et al. (2008) - ``"Kjeldsen2008_scaling"``: apply a correction based on Kjeldsen et al. (2008). The exponent is based on a scaling relation from Sonoi et al. (2015). - ``"Kjeldsen2008_2"``: apply a correction based on Kjeldsen et al. (2008). The exponent is a free parameter. - ``"Ball2014"``: apply a one-term correction based on Ball and Gizon (2014) - ``"Ball2014_2"``: apply a two-term correction based on Ball and Gizon (2014) - ``"Sonoi2015"``: apply a correction based on Sonoi et al. (2015) - ``"Sonoi2015_scaling"``: apply a correction based on Sonoi et al. (2015) The exponent is based on a scaling relation from Sonoi et al. (2015). - ``"Sonoi2015_2"``: apply a correction based on Sonoi et al. (2015) The exponent is a free parameter. :param a: parameters which intervene in the surface correction. According to the correction they take on the following meanings: - ``"Kjeldsen2008"``: a[0]*freq**b_Kjeldsen2008 - ``"Kjeldsen2008_scaling"``: a[0]*freq**b_scaling - ``"Kjeldsen2008_2"``: a[0]*freq**a[1] - ``"Ball2014"``: a[0]*freq**3/I - ``"Ball2014_2"``: a[0]*freq**3/I + a[1]/(freq*I) - ``"Sonoi2015"``: a[0]*[1 - 1/(1 + (nu/numax)**beta_Sonoi2015)] - ``"Sonoi2015_scaling"``: a[0]*[1 - 1/(1 + (nu/numax)**beta_scaling)] - ``"Sonoi2015_2"``: a[0]*[1 - 1/(1 + (nu/numax)**a[1])] :type surface_option: string :type a: array-like :return: surface corrections on the model frequencies :rtype: np.array .. note:: The array operations lead to the creation of a new array with the result, which avoids modifications of the original frequencies and inertias. """ # easy exit: if (surface_option is None): return np.zeros(len(self.modes), dtype=ftype) freq = self.glb[ifreq_ref] * self.modes['freq'] / self.numax if (surface_option == "Kjeldsen2008"): return a[0] * freq ** config.b_Kjeldsen2008 if (surface_option == "Kjeldsen2008_scaling"): return a[0] * freq ** self.b_Kjeldsen2008 if (surface_option == "Kjeldsen2008_2"): return a[0] * freq ** a[1] if (surface_option == "Ball2014"): return a[0] * freq ** 3 / self.modes['inertia'] if (surface_option == "Ball2014_2"): return a[0] * freq ** 3 / self.modes['inertia'] \ + a[1] / (freq * self.modes['inertia']) if (surface_option == "Sonoi2015"): return a[0] * (1.0 - 1.0 / (1.0 + freq ** config.beta_Sonoi2015)) if (surface_option == "Sonoi2015_scaling"): return a[0] * (1.0 - 1.0 / (1.0 + freq ** self.beta_Sonoi2015)) if (surface_option == "Sonoi2015_2"): return a[0] * (1.0 - 1.0 / (1.0 + freq ** a[1])) raise ValueError("Unknown surface correction: " + surface_option)
@property def b_Kjeldsen2008(self): """ Return the exponent for the Kjeldsen et al. (2008) surface correction recipe, as calculated based on the Sonoi et al. (2015) scaling relation. :return: the Kjeldsen et al. exponent :rtype: float """ return 10.0 ** (-3.16 * self.string_to_param("log_Teff") + 0.184 * self.string_to_param("log_g") + 11.7) @property def beta_Sonoi2015(self): """ Return the exponent for the Sonoi et al. (2015) surface correction recipe, as calculated based on the Sonoi et al. (2015) scaling relation. :return: the Kjeldsen et al. exponent :rtype: float """ return 10.0 ** (-3.86 * self.string_to_param("log_Teff") + 0.235 * self.string_to_param("log_g") + 14.2)
[docs] def multiply_modes(self, constant): """ Multiply the frequencies by constant. :param constant: constant by which the mode frequencies are multiplied :type constant: float """ # NOTE: inertias are non-dimensionless, so they don't change. # (this has been tested numerically) self.modes['freq'] *= constant
[docs] def find_mode(self, ntarget, ltarget): """ Find a mode with specific n and l values. :param ntarget: target n value :param ltarget: target l value :type ntarget: int :type ltarget: int :return: the frequency of the mode :rtype: float """ size = len(self.modes) for i in range(size): if ((self.modes['n'][i] == ntarget) and (self.modes['l'][i] == ltarget)): return self.modes['freq'][i] return np.nan
[docs] def find_mode_range(self): """ Find n and l ranges of the modes in the model. :return: the n and l ranges of the modes :rtype: int, int, int, int """ if (len(self.modes) < 1): return -1, -1, -1, -1 nmin = np.nanmin(self.modes['n']) nmax = np.nanmax(self.modes['n']) lmin = np.nanmin(self.modes['l']) lmax = np.nanmax(self.modes['l']) return nmin, nmax, lmin, lmax
[docs] def find_large_separation(self): """ Find large frequency separation using only radial modes. :return: the large frequency separation :rtype: float """ ind = (self.modes['l'] == 0).flat n = self.modes['n'].compress(ind) if (len(np.unique(n)) < 2): return np.nan freq = self.modes['freq'].compress(ind) coeff = np.polyfit(n, freq, 1) return coeff[0]
[docs] def find_epsilon(self, ltarget): """ Find epsilon, the constant offset in a simplified version of Tassoul's asymptotic formula: :math:`\\nu_n = \\Delta \\nu (n + \\varepsilon)` :param ltarget: target l value. Only modes with this l value will be used in obtaining epsilon. :type ltarget: int :return: the constant offset :rtype: float """ dnu = self.find_large_separation() one = n = nu = 0.0 for i in range(len(self.modes)): if (self.modes['l'][i] != ltarget): continue one += 1.0 n += self.modes['n'][i] nu += self.modes['freq'][i] if (one == 0.0): return 0.0 else: return (nu / dnu - n) / one
@property def FeH(self): """ Find [Fe/H] value for model. The conversion from (Xs,Zs) to [Fe/H] is performed using the following formula: :math:`\\mathrm{[Fe/H] = \\frac{[M/H]}{A_{FeH}} \ = \\frac{1}{A_{FeH}} \\log_{10} \ \\left(\\frac{z/x}{z_{\\odot}/x_{\\odot}} \\right)}` :return: the :math:`\\mathrm{[Fe/H]}` value :rtype: float .. note:: The relevant values are given in :py:mod:`constants` """ try: return self.MH / constants.A_FeH except ValueError: return log0 # a rather low value @property def FeH0(self): """ Find initial [Fe/H] value for model. The conversion from (X,Z) to [Fe/H] is performed using the following formula: :math:`\\mathrm{[Fe/H] = \\frac{[M/H]}{A_{FeH}} \ = \\frac{1}{A_{FeH}} \\log_{10} \ \\left(\\frac{z/x}{z_{\\odot}/x_{\\odot}} \\right)}` :return: the initial :math:`\\mathrm{[Fe/H]}` value :rtype: float .. note:: The relevant values are given in :py:mod:`constants` """ try: return self.MH0 / constants.A_FeH except ValueError: return log0 # a rather low value @property def MH(self): """ Find [M/H] value for model. The conversion from (Xs,Zs) to [M/H] is performed using the following formula: :math:`\\mathrm{[M/H] = \\log_{10} \\left(\\frac{z/x}{z_{\\odot}/x_{\\odot}} \\right)}` :return: the :math:`\\mathrm{[M/H]}` value :rtype: float .. note:: The relevant values are given in :py:mod:`constants` """ available_keys = user_params_index.keys() have_Zs = False have_Ys = False have_Xs = False have_ZXs = False if 'Zs' in available_keys: Zs = self.glb[user_params_index["Zs"]] have_Zs = True if 'Ys' in available_keys: Ys = self.glb[user_params_index["Ys"]] have_Ys = True if 'Xs' in available_keys: Xs = self.glb[user_params_index["Xs"]] have_Xs = True if 'ZXs' in available_keys: ZXs = self.glb[user_params_index["ZXs"]] have_ZXs = True if 'ZHsurf' in available_keys: ZXs = self.glb[user_params_index["ZHsurf"]] have_ZXs = True if not have_Zs: if have_Xs and have_Ys: Zs = 1 - self.glb[user_params_index["Xs"]] - self.glb[user_params_index["Ys"]] have_Zs = True if have_ZXs and have_Xs: Zs = ZXs * Xs have_Zs = True try: return math.log10(Zs * constants.solar_x / (Xs * constants.solar_z)) except ValueError: return log0 # a rather low value @property def MH0(self): """ Find initial [M/H] value for model. The conversion from (X,Z) to [M/H] is performed using the following formula: :math:`\\mathrm{[M/H] = \\log_{10} \\left(\\frac{z/x}{z_{\\odot}/x_{\\odot}} \\right)}` :return: the initial :math:`\\mathrm{[M/H]}` value :rtype: float .. note:: The relevant values are given in :py:mod:`constants` """ try: return math.log10(self.glb[iz0] * constants.solar_x / (self.glb[ix0] * constants.solar_z)) except ValueError: return log0 # a rather low value @property def zsx_s(self): """ Find the Zs/Xs value :return: the Zs/Xs value :rtype: float """ return self.glb[user_params_index["Zs"]] / self.glb[user_params_index["Xs"]] @property def zsx_0(self): """ Find the Z0/X0 value :return: the Z0/X0 value :rtype: float """ return self.glb[iz0] / self.glb[ix0] @property def numax(self): """ Find :math:`\\nu_{\\mathrm{max}}` for model. The :math:`\\nu_{\\mathrm{max}}` value is obtained from the following scaling relation: :math:`\\frac{\\nu_{\\mathrm{max}}}{\\nu_{\\mathrm{max},\\odot}} \ = \\left(\\frac{M}{M_{\\odot}}\\right) \ \\left(\\frac{R}{R_{\\odot}}\\right)^{-2} \ \\left(\\frac{T_{\\mathrm{eff}}}{T_{\\mathrm{eff},\\odot}}\\right)^{-1/2}` :return: the :math:`\\nu_{\\mathrm{max}}` value :rtype: float .. note:: The relevant values are given in :py:mod:`constants` """ return constants.solar_numax * (self.glb[imass] / constants.solar_mass) \ / ((self.glb[iradius] / constants.solar_radius) ** 2 \ * math.sqrt(self.glb[itemperature] / constants.solar_temperature)) @property def cutoff(self): """ Find :math:`\\nu_{\\mathrm{cut-off}}` for model. The :math:`\\nu_{\\mathrm{cut-off}}` value is obtained from the following scaling relation: :math:`\\frac{\\nu_{\\mathrm{cut-off}}}{\\nu_{\\mathrm{cut-off},\\odot}} \ = \\left(\\frac{M}{M_{\\odot}}\\right) \ \\left(\\frac{R}{R_{\\odot}}\\right)^2 \ \\left(\\frac{T_{\\mathrm{eff}}}{T_{\\mathrm{eff},\\odot}}\\right)^{-1/2}` :return: the :math:`\\nu_{\\mathrm{cut-off}}` value :rtype: float .. note:: The relevant values are given in :py:mod:`constants` """ return constants.solar_cutoff * (self.glb[imass] / constants.solar_mass) \ / ((self.glb[iradius] / constants.solar_radius) ** 2 \ * math.sqrt(self.glb[itemperature] / constants.solar_temperature))
[docs] def freq_sorted(self): """ Check to see if the frequencies are in ascending order for each l value. :return: ``True`` if the frequencies are in ascending order. :rtype: boolean """ for i in range(len(self.modes) - 1): # if (self.modes['l'][i] > 0): continue if (self.modes['l'][i] != self.modes['l'][i + 1]): continue if (self.modes['freq'][i] > self.modes['freq'][i + 1]): return False return True
[docs] def print_me(self): """ Print classical and seismic characteristics of the model to standard output.""" print("----- Model: " + str(self.name) + " -----") print("Mass (in M_sun): %.5f" % (self.glb[imass] / constants.solar_mass)) print("Radius (in R_sun): %.5f" % (self.glb[iradius] / constants.solar_radius)) print("Reference frequency (in uHz): %.3f" % self.glb[ifreq_ref]) print("Temperature (in K): %.1f" % self.glb[itemperature]) print("Luminosity (in L_sun): %.3g" % (self.glb[iluminosity] / constants.solar_luminosity)) print("Age (in Myrs): %.2f" % self.glb[iage]) print("Age parameter: %.2f" % self.glb[iage_adim]) print("Z: %.4f" % self.glb[iz0]) print("X: %.4f" % self.glb[ix0]) for (name, latex_name) in config.user_params: print("{0:29} {1:.5e}".format(name, self.glb[user_params_index[name]])) print("Modes (in muHz):") size = self.modes.shape[0] for i in range(size): print(" (n,l,freq,IK) = (%d, %d, %.15f, %.5e)" % \ (self.modes['n'][i], self.modes['l'][i], \ self.modes['freq'][i] * self.glb[ifreq_ref], \ self.modes['inertia'][i]))
[docs] class Track: """ An evolutionary track. """ def __init__(self, aModel, grid_params): """ :param aModel: first model to be added to evolutionary track (it does not need to be the youngest model in an evolutionary sequence). This Model is used to obtain the relevant parameters for the evolutionary track (as given by the :py:data:`grid_params` variable). :param grid_params: list of strings which are the names of the parameters which describe the evolutionary track. :type aModel: :py:class:`Model` :type grid_params: list of strings """ self.grid_params = grid_params """Names of the parameters used to construct the grid""" self.params = np.array(utilities.my_map(aModel.string_to_param, self.grid_params)) """Parameters which characterise this evolutionary track""" self.nmodes = len(aModel.modes) """Total number pulsation modes from all of the models in this evolutionary track""" self.names = [aModel.name, ] """ List of model names """ self.glb = aModel.glb.reshape((1, nglb)) """ Global properties of the models """ self.modes = aModel.modes """ array containing the modes (n, l, freq, inertia) of all of the models """ self.mode_indices = [0, len(aModel.modes)] """ starting indices in :py:data:`Track.modes` array corresponding to a given model """ del aModel # clean up to avoid using too much memory def __del__(self): """ Remove a track and clear the arrays it contained. """ del self.grid_params del self.params del self.nmodes del self.names del self.glb del self.modes del self.mode_indices
[docs] def append(self, aModel): """ Append a model to the evolutionary track. :param aModel: model which is being appended to the track :type aModel: :py:class:`Model` """ self.nmodes += len(aModel.modes) self.names.append(aModel.name) self.glb = np.append(self.glb, aModel.glb.reshape((1, nglb)), axis=0) self.modes = np.append(self.modes, aModel.modes) self.mode_indices.append(self.nmodes) del aModel # clean up to avoid using too much memory
[docs] def append_track(self, aTrack): """ Append a track to the current evolutionary track (i.e. combine the two tracks), and remove the track which has been appended. The resultant track is then sorted according to age. :param aModel: model which is being appended to the track :type aModel: :py:class:`Model` """ self.names += aTrack.names self.glb = np.concatenate((self.glb, aTrack.glb), axis=0) self.modes = np.concatenate((self.modes, aTrack.modes), axis=0) self.mode_indices += [self.nmodes + i for i in aTrack.mode_indices[1:]] self.nmodes += aTrack.nmodes self.sort()
[docs] def matches(self, aModel, params_tol=None): """ Check to see if a model matches the evolutionary track and can therefore be included in the track. :param aModel: input model being tested :type aModel: :py:class:`Model` :param params_tol: the tolerance on each parameter. ``None`` means no tolerance (i.e. the parameters have to be exactly the same). :type params_tol: float np.array :return: ``True`` only if all of the parameters of the input model match those of the evolutionary track within the provided tolerances :rtype: boolean """ params_bis = np.array(utilities.my_map(aModel.string_to_param, self.grid_params)) if (params_tol is None): return np.all(self.params == params_bis) else: return np.all(np.abs(self.params - params_bis) <= params_tol)
[docs] def sort(self): """Sort models within evolutionary track according to age.""" index = np.argsort(self.glb[:, iage]) # sort model names temp = utilities.my_map(lambda i: self.names[i], index) del self.names self.names = temp # sort global parameters temp = self.glb[index] del self.glb self.glb = temp # sort modes and mode indices: temp = np.empty(self.modes.shape, dtype=modetype) ntot = 0 temp2 = [ntot, ] for i in index: nmodes = self.mode_indices[i + 1] - self.mode_indices[i] temp[ntot:ntot + nmodes] = self.modes[self.mode_indices[i]:self.mode_indices[i + 1]] ntot += nmodes temp2.append(ntot) del self.modes del self.mode_indices self.modes = temp self.mode_indices = temp2
[docs] def is_sorted(self): """ Check to see of models are in ascending order according to age. :return: ``True`` if the models ar in order of increasing age :rtype: boolean """ return all(self.glb[i, iage] < self.glb[i + 1, iage] for i in range(len(self.names) - 1))
[docs] def is_sorted_adim(self): """ Check to see of models are in ascending order according to age. :return: ``True`` if the models ar in order of increasing age :rtype: boolean """ return all(self.glb[i, iage_adim] < self.glb[i + 1, iage_adim] for i in range(len(self.names) - 1))
[docs] def freq_sorted(self, imodel): """ Check to see if the frequencies of model i are in ascending order for each l value. :param imodel: the number of the model to be checked :type imodel: int :return: ``True`` if the frequencies are in ascending order. :rtype: boolean """ for i in range(self.mode_indices[imodel], self.mode_indices[imodel + 1] - 1): # if (self.modes['l'][i] > 0): continue if (self.modes['l'][i] != self.modes['l'][i + 1]): continue if (self.modes['freq'][i] > self.modes['freq'][i + 1]): return False return True
[docs] def find_epsilon(self, imodel, ltarget): """ Find epsilon, the constant offset in a simplified version of Tassoul's asymptotic formula: :math:`\\nu_n = \\Delta \\nu (n + \\varepsilon)` :param imodel: model number :type imodel: ltarget :param ltarget: target l value. Only modes with this l value will be used in obtaining epsilon. :type ltarget: int :return: the constant offset :rtype: float """ dnu = self.find_large_separation(imodel) one = n = nu = 0.0 for i in range(self.mode_indices[imodel], self.mode_indices[imodel + 1]): if (self.modes['l'][i] != ltarget): continue one += 1.0 n += self.modes['n'][i] nu += self.modes['freq'][i] if (one == 0.0): return 0.0 else: return (nu / dnu - n) / one
[docs] def find_large_separation(self, imodel): """ Find large frequency separation using only radial modes. :param imodel: model number :type imodel: ltarget :return: the large frequency separation :rtype: float """ istart = self.mode_indices[imodel] istop = self.mode_indices[imodel + 1] ind = (self.modes['l'][istart:istop] == 0).flat n = self.modes['n'][istart:istop].compress(ind) if (len(n) < 2): return np.nan freq = self.modes['freq'][istart:istop].compress(ind) coeff = np.polyfit(n, freq, 1) return coeff[0]
[docs] def duplicate_ages(self): """ Check to see if the evolutionary track contains models with duplicate ages. :return: ``True`` if there are duplicate age(s) :rtype: boolean .. warning:: This method should only be applied after the track has been sorted. """ return any(self.glb[i, iage] == self.glb[i + 1, iage] for i in range(len(self.names) - 1))
[docs] def remove_duplicate_ages(self): """ Removes models with duplicate ages from the evolutionary track. :return: ``True`` if models with duplicate age(s) have been removed :rtype: boolean .. warning:: This method should only be applied after the track has been sorted. """ # produce list of models to be removed: remove_models = [i for i in range(1, len(self.names)) \ if self.glb[i - 1, iage] == self.glb[i, iage]] if (len(remove_models) == 0): return False else: # remove duplicate models: mode_indices_rm = [] for i in remove_models: mode_indices_rm += list(range(self.mode_indices[i], self.mode_indices[i + 1])) for i in remove_models[::-1]: del self.names[i] self.glb = np.delete(self.glb, remove_models, 0) self.modes = np.delete(self.modes, mode_indices_rm, 0) ntot = self.mode_indices[1] new_indices = [0, ntot] # NOTE: the first model is kept (by construction) for i in range(1, len(self.mode_indices) - 1): if (i not in remove_models): ntot += self.mode_indices[i + 1] - self.mode_indices[i] new_indices.append(ntot) self.mode_indices = new_indices self.nmodes = ntot return True
[docs] def age_adim_to_age(self): """ Replace dimensionless age parameter by the physical age. .. warning:: This method should only be applied after the track has been sorted (according to age). """ self.glb[:, iage_adim] = self.glb[:, iage]
[docs] def age_adim_to_scale_age(self): """ Replace dimensionless age parameter by the physical age scaled from 0 to 1. .. warning:: This method should only be applied after the track has been sorted (according to age). """ age_start = self.glb[0, iage] age_stop = self.glb[-1, iage] self.glb[:, iage_adim] = (self.glb[:, iage] - age_start) / (age_stop - age_start)
[docs] def age_adim_to_scale_Xc(self): """ Replace dimensionless age parameter by the central hydrogen abundance, Xc, scaled from 0 to 1. .. warning:: This method should only be applied after the track has been sorted (according to age). """ # Raises a key error if Xc isn't in the grid, # which should be informative enough ixc = user_params_index["Xc"] Xc_start = self.glb[0, ixc] Xc_stop = self.glb[-1, ixc] self.glb[:, iage_adim] = (self.glb[:, ixc] - Xc_start) / (Xc_stop - Xc_start)
[docs] def interpolate_model(self, age): """ Return a model at a given age which is obtained using linear interpolation. :param age: age of desired model in :math:`\\mathrm{Myrs}` :type age: float :return: the interpolated model :rtype: :py:class:`Model` .. warning:: This method assumes the track is sorted, since it applies a binary search algorithm for increased efficiency. """ # easy exit: if (age < self.glb[0, iage]): return None if (age > self.glb[-1, iage]): return None istart = 0 istop = len(self.names) - 1 while (istop > istart + 1): itemp = (istop + istart) // 2 if (age < self.glb[itemp, iage]): istop = itemp else: istart = itemp mu = (age - self.glb[istart, iage]) \ / (self.glb[istop, iage] - self.glb[istart, iage]) return combine_models( Model(self.glb[istart], _modes=self.modes[self.mode_indices[istart]:self.mode_indices[istart + 1]]), 1.0 - mu, \ Model(self.glb[istop], _modes=self.modes[self.mode_indices[istop]:self.mode_indices[istop + 1]]), mu)
[docs] def find_combination(self, age, coef): """ Return a model combination at a given age which is obtained using linear interpolation. :param age: age of desired model in :math:`\\mathrm{Myrs}` :param coef: coefficient which multiplies this combination :type age: float :type coef: float :return: pairs composed of an interpolation coefficient and a model name :rtype: tuple of (float, string) .. warning:: This method assumes the track is sorted, since it applies a binary search algorithm for increased efficiency. """ # easy exit: if (age < self.glb[0, iage]): return None if (age > self.glb[-1, iage]): return None istart = 0 istop = len(self.names) - 1 while (istop > istart + 1): itemp = (istop + istart) // 2 if (age < self.glb[itemp, iage]): istop = itemp else: istart = itemp mu = (age - self.glb[istart, iage]) \ / (self.glb[istop, iage] - self.glb[istart, iage]) return ((coef * (1.0 - mu), self.names[istart]), (coef * mu, self.names[istop]))
[docs] def find_modes(self, ntarget, ltarget): """ Return two lists, one with the ages of the models and the other with the mode non-dimensional frequencies corresponding to target n and l values. This function is useful for seeing how the frequency of a particular mode changes with stellar age. :param ntarget: target n value :param ltarget: target l value :type ntarget: int :type ltarget: int :return: lists of ages and frequencies :rtype: list, list """ ages = list(self.glb[:, iage]) freqs = [] for i in range(len(self.names)): for j in range(self.mode_indices[i], self.mode_indices[i + 1]): if ((self.modes['n'][j] == ntarget) and (self.modes['l'][j] == ltarget)): freqs.append(self.modes['freq'][j]) else: freqs.append(np.nan) return ages, freqs
[docs] def find_modes_dim(self, ntarget, ltarget): """ Return two lists, one with the ages of the models and the other with the mode dimensional frequencies corresponding to target n and l values. This function is useful for seeing how the frequency of a particular mode changes with stellar age. :param ntarget: target n value :param ltarget: target l value :type ntarget: int :type ltarget: int :return: lists of ages and frequencies :rtype: list, list """ ages = list(self.glb[:, iage]) freqs = [] for i in range(len(self.names)): for j in range(self.mode_indices[i], self.mode_indices[i + 1]): if ((self.modes['n'][j] == ntarget) and (self.modes['l'][j] == ltarget)): freqs.append(self.modes['freq'][j] * self.glb[i, ifreq_ref]) else: freqs.append(np.nan) return ages, freqs
[docs] def find_mode_range(self): """ Find n and l ranges of modes in models :return: the n and l ranges :rtype: int, int, int, int """ if (len(self.modes) < 1): return -1, -1, -1, -1 else: return np.min(self.modes['n']), np.max(self.modes['n']), \ np.min(self.modes['l']), np.max(self.modes['l'])
@property def age_range(self): """ Provides the age range for an evolutionary track. """ return abs(self.glb[-1, iage] - self.glb[0, iage]) @property def age_lower(self): """ Provides the lowest age an evolutionary track. """ return min(self.glb[0, iage], self.glb[-1, iage]) @property def age_upper(self): """ Provides the highest age for an evolutionary track. """ return max(self.glb[0, iage], self.glb[-1, iage])
[docs] def test_interpolation(self, nincr): """ Test accuracy of interpolation along evolutionary track. This method removes every other model and retrieves its frequencies by interpolation from neighbouring models. The accuracy of the interpolated frequencies and global parameters are tested by carrying out comparisons with the original models. :param nincr: increment with which to carry out the interpolation. By comparing results for different values of ``nincr``, one can evaluate how the interpolation error depends on the size of the interval over which the interpolation is carried out. :type nincr: int :return: the interpolation errors :rtype: np.array """ # initialisation nmodels = len(self.names) ndim = len(self.params) + 1 result = np.zeros((nmodels - 2 * nincr, ndim + nglb + 6), dtype=gtype) # loop through all models: for i in range(nincr, nmodels - nincr): # carry out interpolation mu = (self.glb[i, iage] - self.glb[i - nincr, iage]) \ / (self.glb[i + nincr, iage] - self.glb[i - nincr, iage]) aModel = combine_models(Model(self.glb[i - nincr], _modes=self.modes[ self.mode_indices[i - nincr]:self.mode_indices[ i - nincr + 1]]), 1.0 - mu, \ Model(self.glb[i + nincr], _modes=self.modes[ self.mode_indices[i + nincr]:self.mode_indices[ i + nincr + 1]]), mu) result[i - nincr, 0:ndim - 1] = self.params result[i - nincr, ndim - 1] = self.glb[i, iage] result[i - nincr, ndim:ndim + nglb + 6] = compare_models(aModel, Model(self.glb[i], _modes=self.modes[ self.mode_indices[ i]: self.mode_indices[ i + 1]])) return result
[docs] class Model_grid: """ A grid of models. """ def __init__(self): self.ndim = 0 """ Number of dimensions for the grid (excluding age), as based on the :py:data:`Model_grid.grid_params` variable """ self.tracks = [] """List of evolutionary tracks contained in the grid.""" self.ndx = [] """List containing track indices""" self.grid = None """Array containing the grid parameters for each evolutionary track (excluding age).""" self.tessellation = None """Object containing the tessellation of the grid used for interpolation.""" self.grid_params = None """ Set of parameters (excluding age) used to construct the grid and do interpolations. .. note:: For best interpolation results, these parameters should be comparable. """ self.prefix = None """Root folder with grid of models (including final slash).""" self.postfix = ".freq" """Last part of the filenames which contain the model frequencies (default = ".freq").""" self.user_params = config.user_params """ The set of user parameters involved in the grid. This is to avoid having a different set of user parameters in `AIMS_configure.py` """ self.distort_mat = None """ Transformation matrix used to break the Cartesian character of the grid and reduce computation time """
[docs] def read_model_list(self, filename): """ Read list of models from a file and construct a grid. If the mode format is Aldo, an entirely different strategy must be used. :param filename: name of the file with the list. :type filename: string """ if (config.mode_format == "Aldo"): self.read_model_list_Aldo(filename) else: self.read_model_list_standard(filename)
[docs] def read_model_list_standard(self, filename): """ Read list of models from a file and construct a grid. :param filename: name of the file with the list. The first line of this file should contain a prefix which is typically the root folder of the grid of models. This followed by a file with multiple columns. The first 9 contain the following information for each model: 1. the second part of the path. When concatenated with the prefix on the first line, this should give the full path to the model. 2. The stellar mass in :math:`\\mathrm{g}` 3. The stellar radius in :math:`\\mathrm{cm}` 4. The stellar luminosity in :math:`\\mathrm{g.cm^2.s^{-3}}` 5. The metallicity 6. The hydrogen content 7. The stellar age in :math:`\\mathrm{Myrs}` 8. The effective temperature in :math:`\\mathrm{K}` 9. A dimensionless age parameter The following columns contain the parameters specified in the :py:data:`AIMS_configure.user_params` variable. :type filename: string """ self.grid_params = config.grid_params # set the correct dimension: self.ndim = len(self.grid_params) # set prefix and postfix: with open(filename, "r") as listfile: line = listfile.readline().strip() columns = line.split() if len(columns) < 1: raise IOError("Erroneous first line in %s." % filename) self.prefix = columns[0] if (len(columns) > 1): self.postfix = columns[1] if (self.postfix == "None"): self.postfix = "" # read list with numpy: print("Reading list file") perm = [imass, iradius, iluminosity, iz0, ix0, iage, itemperature, iage_adim] \ + [user_params_index[name_pair[0]] for name_pair in config.user_params] invperm = [1, ] * (len(perm) + 1) for i in range(len(perm)): invperm[perm[i]] = i + 1 glb = np.loadtxt(filename, skiprows=1, usecols=invperm, dtype=gtype) names = np.loadtxt(filename, skiprows=1, usecols=(0,), dtype=filename_dtype) # sort list: this separates the evolutionary tracks (which still need # to be partitionned) print("Sorting models") grid_list = [age_str, ] + list(self.grid_params) grid_temp = np.array([utilities.my_map(Model(glb[i]).string_to_param, \ grid_list) for i in range(glb.shape[0])]) ind = np.lexsort([grid_temp[:, i] for i in range(grid_temp.shape[1])]) params_span = np.array([np.max(grid_temp[:, i]) - np.min(grid_temp[:, i]) \ for i in range(1, grid_temp.shape[1])]) del grid_temp # sanity check: for i in range(len(self.grid_params)): span = params_span[i] if np.isnan(span): raise ValueError("The values of some models parameters are NaN.") if np.isinf(span): raise ValueError("The values of some models parameters are infinite.") if span == 0.0: raise ValueError("ERROR: parameter %s is constant in your grid" % self.grid_params[i] + "and cannot be used as a grid parameter.") # create tracks: print("Creating evolutionary tracks") nmodes = 0 models_small_spectra = [] for nmodels in range(glb.shape[0]): i = ind[nmodels] aModel = Model(glb[i], _name=names[i]) exceed_freqlim = aModel.read_file(self.prefix + names[i] + self.postfix) aModel.multiply_modes(1.0 / aModel.glb[ifreq_ref]) # make frequencies non-dimensional aModel.sort_modes() aModel.remove_duplicate_modes() if ((len(self.tracks) > 0) and (self.tracks[-1].matches(aModel))): self.tracks[-1].append(aModel) else: aTrack = Track(aModel, self.grid_params) self.tracks.append(aTrack) nmodes += len(aModel.modes) if (not exceed_freqlim): models_small_spectra.append(aModel.name) if (not config.batch): print("%d %d %d" % (len(self.tracks), nmodels + 1, nmodes)) print('\033[2A') # backup two line - might not work in all terminals print("%d %d %d" % (len(self.tracks), glb.shape[0], nmodes)) del glb del names del ind # merge tracks which are too close: imerge = 0 i = 1 params_span *= eps while (i < len(self.tracks)): aModel = Model(self.tracks[i].glb[0]) for j in range(i): if (self.tracks[j].matches(aModel, params_tol=params_span)): self.tracks[j].append_track(self.tracks[i]) del self.tracks[i] imerge += 1 break else: i += 1 print("I merged %d track(s)" % (imerge)) # write list of models with spectra which are too small in a file: with open("models_small_spectra", "w") as output: for name in models_small_spectra: output.write(name + "\n") # sort tracks: for track in self.tracks: track.sort() # remove duplicate models: for track in self.tracks: if track.remove_duplicate_ages(): print("WARNING: the track %s = %s" % (str(track.grid_params), str(track.params))) print(" has models with the same age. Removing duplicate models.") # update list of indices: self.ndx = range(len(self.tracks)) # need to create grid from scratch since tracks have been sorted. self.grid = np.asarray([track.params for track in self.tracks])
[docs] def read_model_list_Aldo(self, filename): """ Read list of models in Aldo format from a file and construct a grid. A track in Aldo format contains 3 parts: 1. A three line header which is discarded. 2. A set of four lines, for l=0 to 3, which describe the mode structure. The first columns gives the start value of n, whereas the second column gives the number of modes for that particular l value. 3. A table where each line corresponds to a model and the different columns correspond to the model ID, global quantities, and pulsation frequencies and inertias. :param filename: name of the file with the list of track files. :type filename: string """ # sanity check: if (len(config.user_params) != 5): raise ValueError('user_params should have length 5 for Aldo format but got %i.' % len(config.user_params)) for key in ['Xc', 'Zc', 'Xs', 'Zs', 'Mass_true']: if not key in user_params_index: raise KeyError("%s missing from user_params variable but needed for Aldo format." % key) self.grid_params = config.grid_params # set the correct dimension: self.ndim = len(self.grid_params) # set prefix and postfix: with open(filename, "r") as listfile: line = listfile.readline().strip() columns = line.split() if (len(columns) < 1): raise IOError("Erroneous first line in %s." % filename) self.prefix = columns[0] self.postfix = "" # read models and put them into evolutionary tracks: nmodels = 0 nmodes = 0 models_small_spectra = [] norder = np.zeros((4, 2), dtype=ntype) for line in listfile: if line[0] == "#": continue line = line.strip() if len(line) == 0: continue # read track file: x0 = None z0 = None Mass0 = None with open(self.prefix + line, "r") as trackfile: # skip header: for i in range(3): trackfile.readline() # read mode structure ntot = 0 for l in range(4): columns = trackfile.readline().strip().split() norder[l, 0] = int(columns[0]) norder[l, 1] = int(columns[1]) ntot += norder[l, 1] # read models with numpy: glbs = np.loadtxt(self.prefix + line, skiprows=7, usecols=range(1, 11 + 2 * ntot), dtype=gtype) names = np.loadtxt(self.prefix + line, skiprows=7, usecols=(0,), dtype=filename_dtype) # read models for i in range(glbs.shape[0]): if (x0 is None): x0 = glbs[i, 7] if (z0 is None): z0 = glbs[i, 6] if (Mass0 is None): Mass0 = glbs[i, 1] * constants.solar_mass glb = np.empty((nglb,), dtype=gtype) glb[imass] = Mass0 glb[iradius] = glbs[i, 4] * constants.solar_radius glb[iluminosity] = constants.solar_luminosity * 10.0 ** glbs[i, 3] glb[iz0] = z0 glb[ix0] = x0 glb[iage] = glbs[i, 0] glb[itemperature] = glbs[i, 2] glb[iage_adim] = glbs[i, 0] # may be replace later on glb[user_params_index["Mass_true"]] = glbs[i, 1] glb[user_params_index["Xs"]] = glbs[i, 7] glb[user_params_index["Zs"]] = glbs[i, 6] glb[user_params_index["Xc"]] = glbs[i, 8] glb[user_params_index["Zc"]] = 1.0 - glbs[i, 9] - glbs[i, 8] aModel = Model(glb, _name=names[i]) exceed_freqlim = aModel.read_modes_Aldo(norder, glbs[i, 10:]) aModel.multiply_modes(1.0 / aModel.glb[ifreq_ref]) # make frequencies non-dimensional aModel.sort_modes() # We're assuming the models are grouped together in tracks, so # we only need to check the previous track: if len(self.tracks) > 0 and self.tracks[-1].matches(aModel): self.tracks[-1].append(aModel) else: aTrack = Track(aModel, self.grid_params) self.tracks.append(aTrack) nmodels += 1 nmodes += len(aModel.modes) if not exceed_freqlim: models_small_spectra.append(aModel.name) if not config.batch: print("%d %d %d" % (len(self.tracks), nmodels, nmodes)) print('\033[2A') # backup two line - might not work in all terminals del glbs del names print("%d %d %d" % (len(self.tracks), nmodels, nmodes)) # find span for each grid parameter prior to merging tracks: grid_temp = np.asarray([track.params for track in self.tracks]) params_span = np.array([np.max(grid_temp[:, i]) - np.min(grid_temp[:, i]) \ for i in range(1, grid_temp.shape[1])]) del grid_temp # merge tracks which are too close: imerge = 0 i = 1 params_span *= eps while (i < len(self.tracks)): aModel = Model(self.tracks[i].glb[0]) for j in range(i): if (self.tracks[j].matches(aModel, params_tol=params_span)): self.tracks[j].append_track(self.tracks[i]) del self.tracks[i] imerge += 1 break else: i += 1 print("I merged %d track(s)" % (imerge)) # write list of models with spectra which are too small in a file: with open("models_small_spectra", "w") as output: for name in models_small_spectra: output.write(name + "\n") # sort tracks: for track in self.tracks: track.sort() # sanity check: for track in self.tracks: if track.remove_duplicate_ages(): print("WARNING: the track %s = %s" % (str(track.grid_params), str(track.params))) print(" has models with the same age. Removing duplicate models.") # update list of indices: self.ndx = range(len(self.tracks)) # need to create grid from scratch since tracks have been sorted. self.grid = np.asarray([track.params for track in self.tracks])
[docs] def replace_age_adim(self): """ This replaces the dimensionless ages in the tracks according to the :py:data:`replace_age_adim` option chosen in ``AIMS_configure.py``. """ if (config.replace_age_adim is None): return elif (config.replace_age_adim == "age"): for track in self.tracks: track.age_adim_to_age() elif (config.replace_age_adim == "scale_age"): for track in self.tracks: track.age_adim_to_scale_age() elif (config.replace_age_adim == "scale_Xc"): for track in self.tracks: track.age_adim_to_scale_Xc() else: raise ValueError("Unknown option %s for replace_scale_age in AIMS_configure.py." % config.replace_age_adim)
[docs] def check_age_adim(self): """ This checks that all of the tracks are sorted according to dimensionless age parameter. """ for track in self.tracks: if not track.is_sorted_adim(): raise ValueError("ERROR: track(s) not strictly sorted according to dimensionless age")
[docs] def range(self, aParam): """ Find range of values for the input parameter. :param aParam: name of the parameter for which to find the range :type aParam: str .. warning:: The input parameter can only be one of the grid parameters or an age/mHe parameter. """ if (aParam == age_str): param_min = self.tracks[0].age_lower param_max = self.tracks[0].age_upper for track in self.tracks: if (param_min > track.age_lower): param_min = track.age_lower if (param_max < track.age_upper): param_max = track.age_upper else: i = self.grid_params.index(aParam) param_min = self.tracks[0].params[i] param_max = self.tracks[0].params[i] for track in self.tracks: if (param_min > track.params[i]): param_min = track.params[i] if (param_max < track.params[i]): param_max = track.params[i] return [param_min, param_max]
[docs] def tessellate(self): """Apply Delauny triangulation to obtain the grid tessellation.""" if (self.distort_mat is None): self.tessellation = Delaunay(self.grid) else: self.tessellation = Delaunay(np.dot(self.grid, self.distort_mat))
[docs] def plot_tessellation(self): """ Plot the grid tessellation. .. warning:: This only works for two-dimensional tessellations. """ if (self.ndim != 2): print("Only able to plot the tessellation in two dimensions.") return # find bounds: xmin = np.nanmin(self.grid[:, 0]) xmax = np.nanmax(self.grid[:, 0]) ymin = np.nanmin(self.grid[:, 1]) ymax = np.nanmax(self.grid[:, 1]) dx = xmax - xmin dy = ymax - ymin xmin -= dx * 0.03 xmax += dx * 0.03 ymin -= dy * 0.05 ymax += dy * 0.05 # plt.semilogy(self.grid[:,0],self.grid[:,1],'o') plt.plot(self.grid[:, 0], self.grid[:, 1], 'o') plt.triplot(self.grid[:, 0], self.grid[:, 1], self.tessellation.simplices.copy()) plt.xlim((xmin, xmax)) plt.ylim((ymin, ymax)) plt.xlabel(string_to_latex(self.grid_params[0])) plt.ylabel(string_to_latex(self.grid_params[1])) plt.savefig("tessellation.eps")
[docs] def remove_tracks(self, nthreshold): """ Removes stellar evolution tracks with fewer than nthreshold models. :param nthreshold: lower limit on number of models in a stellar evolutionary track :type nthreshold: int :return: True if tracks have been removed and the grid needs to be retessellated :rtype: boolean """ # easy exit: if (nthreshold <= 0): return removedTracks = False for i in range(len(self.tracks) - 1, -1, -1): if (len(self.tracks[i].names) < nthreshold): del self.tracks[i] removedTracks = True # redo tessellation if need be if (removedTracks): print("I removed evolutionary tracks from the grid!") # update list of indices: self.ndx = range(len(self.tracks)) # need to create grid from scratch since tracks have been sorted. self.grid = np.asarray([track.params for track in self.tracks]) return removedTracks
[docs] def distort_grid(self): """ Define distortion matrix with which to distort grid to break its Cartesian character prior to tessellation. This can cause find_simplex to run much faster. """ self.distort_mat = np.dot(make_scale_matrix(self.grid), make_distort_matrix(self.ndim)) print("Distortion matrix condition number: %e" % (np.linalg.cond(self.distort_mat)))
[docs] def test_interpolation(self): """ Test interpolation between different evolutionary tracks in a given grid. :return: The following four items are returned: - the interpolation errors - the first half of the partition (where the interpolation is tested) - the second half of the partition (used to carry out the interpolation) - the tessellation associated with the second half of the partition :rtype: np.array, list, list, tessellation object """ ndx1, ndx2 = self.find_partition() if (self.distort_mat is None): tessellation = Delaunay(self.grid[ndx2, :]) else: tessellation = Delaunay(np.dot(self.grid[ndx2, :], self.distort_mat)) # initialisation results = [] ndim = self.ndim + 1 for j in ndx1: nmodels = len(self.tracks[j].names) aResult = np.empty((nmodels, ndim + nglb + 6), dtype=gtype) pt = list(self.tracks[j].params) + [0.0, ] for i in range(nmodels): aModel1 = Model(self.tracks[j].glb[i], \ _modes=self.tracks[j].modes[self.tracks[j].mode_indices[i]: \ self.tracks[j].mode_indices[i + 1]]) pt[-1] = aModel1.glb[iage] aModel2 = interpolate_model(self, pt, tessellation, ndx2) aResult[i, 0:ndim] = pt if (aModel2 is None): aResult[i, ndim:ndim + nglb + 6] = np.nan else: aResult[i, ndim:ndim + nglb + 6] = compare_models(aModel1, aModel2) results.append(aResult) return results, ndx1, ndx2, tessellation
[docs] def find_partition(self): """ Find a partition of the grid for use with :py:meth:`Model_grid.test_interpolation` :return: a random partition of [0 ... n-1] into two equal halves, where n is the number of tracks in the grid :rtype: two lists of int """ ndx = list(range(len(self.tracks))) random.shuffle(ndx) nn = len(self.tracks) // 2 return ndx[:nn], ndx[nn:]
[docs] def test_freq(self): """ Test to see if frequencies in all of the models of the grid are in ascending order for each l value. :return: The following items are returned - the effective temperatures of the models with frequencies out of order - the luminosities of the models with frequencies out of order - the effective temperatures of the models with sorted frequencies - the luminosities of the models with sorted frequencies :rtype: four lists of floats """ Teffs, Lums, Teffs_out, Lums_out = [], [], [], [] for track in self.tracks: for i in range(len(track.names)): if (not track.freq_sorted(i)): print(track.names[i]) Teffs_out.append(track.glb[i, itemperature]) Lums_out.append(math.log10(track.glb[i, iluminosity] / constants.solar_luminosity)) else: Teffs.append(track.glb[i, itemperature]) Lums.append(math.log10(track.glb[i, iluminosity] / constants.solar_luminosity)) return Teffs_out, Lums_out, Teffs, Lums
[docs] def find_epsilons(self, ltarget): """ Find epsilon values in models from the grid :param ltarget: target l value for which epsilons are being obtained :type ltarget: int :return: the epsilon values :rtype: list of floats """ epsilons = [] for track in self.tracks: for i in range(len(track.names)): epsilon = track.find_epsilon(i, ltarget) if (epsilon != 0.0): epsilons.append(epsilon) return epsilons
[docs] def init_user_param_dict(): """ Initialise the dictionaries which are related to user-defined parameters. For a given parameter, these dictionaries provide the appropriate index for for the :py:data:`Model.glb` array as well as the appropriate latex name. """ i = nlin - len(config.user_params) for (name, latex_name) in config.user_params: user_params_index[name] = i user_params_latex[name] = latex_name i += 1
[docs] def make_distort_matrix(d, theta=0.157): """ Create a distortion matrix which can be used to make the grid more "tessellation-friendly", i.e. which leads to much shorter computation times for finding simplices. :param d: number of dimensions :type d: int :param theta: a small angle :type theta: float :return: a distortion matrix :rtype: 2D float array """ cost = math.cos(theta) sint = math.sin(theta) mat = np.eye(d) aux = np.empty((d, d), dtype=float) for i in range(d): for j in range(i + 1, d): aux[:, :] = 0.0 for k in range(d): if ((k == i) or (k == j)): continue aux[k, k] = 1.0 aux[i, i] = 1.0 aux[j, j] = cost aux[j, i] = 0.0 aux[i, j] = sint mat = np.dot(mat, aux) return mat
[docs] def make_scale_matrix(grid): """ Create a distortion matrix which can be used to make the grid more "tessellation-friendly", i.e. which leads to much shorter computation times for finding simplices. :param grid: set of points used in the construction of the tessellation :type d: 2D float array :return: a distortion matrix :rtype: 2D float array """ eps = 1e-3 d = grid.shape[1] mat = np.eye(d) for i in range(d): x = np.asarray(list(set(grid[:, i]))) x = np.sort(x) x = x[1:] - x[:-1] x = np.sort(x) ind = bisect_right(x, eps) mat[i, i] = 1.0 / x[ind] print("Scale factor %d: %e" % (i, mat[i, i])) return mat
[docs] def combine_models(model1, coef1, model2, coef2): """ Do linear combination of this model with another. This method returns a new model which is the weighted sum of two models for the purposes of model interpolation. The classical parameters are combined in a self-consistent way as are the frequencies. :param model1: first model :param coef1: weighting coefficient applied to first model :param model2: second model :param coef2: weighting coefficient applied to second model :type model1: :py:class:`Model` :type coef1: float :type model2: :py:class:`Model` :type coef2: float :return: the combined model :rtype: :py:class:`Model` .. warning:: One should avoid negative or zero coefficients as these could lead to undefined results. """ # find global parameters (try to be self-consistent): # this first part is simply a linear combination: glb = np.empty((nglb,), dtype=gtype) glb[0:nlin] = coef1 * model1.glb[0:nlin] + coef2 * model2.glb[0:nlin] # this next part depends on previous results: glb[iradius] = (glb[imass] / (coef1 * model1.glb[imass] / model1.glb[iradius] ** 3 + coef2 * model2.glb[imass] / model2.glb[iradius] ** 3)) ** (1.0 / 3.0) cnst1 = model1.glb[iluminosity] / (model1.glb[iradius] ** 2 * model1.glb[itemperature] ** 4) cnst2 = model2.glb[iluminosity] / (model2.glb[iradius] ** 2 * model2.glb[itemperature] ** 4) glb[iluminosity] = (coef1 * cnst1 + coef2 * cnst2) * glb[iradius] ** 2 * glb[itemperature] ** 4 # glb[ifreq_ref] will be correctly defined when the Model() constructor is invoked # interpolate spectra: size3 = min(model1.modes.shape[0], model2.modes.shape[0]) nvalues = np.empty((size3,), dtype=ntype) lvalues = np.empty((size3,), dtype=ltype) fvalues = np.empty((size3,), dtype=ftype) ivalues = np.empty((size3,), dtype=ftype) nvalues, lvalues, fvalues, ivalues, n3 = aims_fortran.combine_modes( \ coef1, model1.modes['n'], model1.modes['l'], model1.modes['freq'], model1.modes['inertia'], \ coef2, model2.modes['n'], model2.modes['l'], model2.modes['freq'], model2.modes['inertia'], \ nvalues, lvalues, fvalues, ivalues) return Model(glb, _modes=list(zip(nvalues[0:n3], lvalues[0:n3], fvalues[0:n3], ivalues[0:n3])))
[docs] def compare_models(model1, model2): """ Compare two models and find the largest frequency different for radial and non-radial modes. :param model1: first model :param model2: second model :type model1: :py:class:`Model` :type model2: :py:class:`Model` :return: a 1D array to be used in ``plot_test_interpolation.py`` with the following measurements of the differences between the two models: - ``result[0]`` = maximum error on the radial modes - ``result[1]`` = RMS error on the radial modes - ``result[2]`` = RMS error on the radial modes near :math:`\\nu_{\\mathrm{max}}` - ``result[3]`` = maximum error on the non radial modes - ``result[4]`` = RMS error on the non radial modes - ``result[5]`` = RMS error on the non radial modes near :math:`\\nu_{\\mathrm{max}}` - ``result[6+[0:nglb]]`` = errors on the global parameters :rtype: np.array """ if (config.interpolation_test_units == "microHz"): C1, C2 = model1.glb[ifreq_ref], model1.glb[ifreq_ref] elif (config.interpolation_test_units is None): C1 = C2 = 1.0 else: raise ValueError("ERROR: Unrecognised units for interpolation_test_units in AIMS_configure.py") # initialisation: n_radial = 0 n_radial_numax = 0 n_non_radial = 0 n_non_radial_numax = 0 result = np.zeros((6 + nglb,), dtype=gtype) # define frequency interval around numax: numax = 0.5 * (C1 * model1.numax / model1.glb[ifreq_ref] \ + C2 * model2.numax / model2.glb[ifreq_ref]) a = 0.8 * numax b = 1.2 * numax # compare frequency spectra: size1 = len(model1.modes) size2 = len(model2.modes) i1 = i2 = 0 while ((i1 < size1) and (i2 < size2)): if (model1.modes['l'][i1] < model2.modes['l'][i2]): i1 += 1; continue if (model1.modes['l'][i1] > model2.modes['l'][i2]): i2 += 1; continue if (model1.modes['n'][i1] < model2.modes['n'][i2]): i1 += 1; continue if (model1.modes['n'][i1] > model2.modes['n'][i2]): i2 += 1; continue # now the two modes have the same n and l values: diff = abs(C1 * model1.modes['freq'][i1] - C2 * model2.modes['freq'][i2]) avg_freq = (C1 * model1.modes['freq'][i1] + C2 * model2.modes['freq'][i2]) / 2.0 if (model1.modes['l'][i1] == 0): if (result[0] < diff): result[0] = diff diff *= diff # square diff result[1] += diff n_radial += 1 # in python, this is called an interval comparison: if (a <= avg_freq <= b): result[2] += diff n_radial_numax += 1 else: if (result[3] < diff): result[3] = diff diff *= diff # square diff result[4] += diff n_non_radial += 1 if (a <= avg_freq <= b): result[5] += diff n_non_radial_numax += 1 i1 += 1 i2 += 1 # avoid divisions by zero: if (n_radial > 0): result[1] = math.sqrt(result[1] / float(n_radial)) else: result[1] = np.nan if (n_radial_numax > 0): result[2] = math.sqrt(result[2] / float(n_radial_numax)) else: result[2] = np.nan if (n_non_radial > 0): result[4] = math.sqrt(result[4] / float(n_non_radial)) else: result[4] = np.nan if (n_non_radial_numax > 0): result[5] = math.sqrt(result[5] / float(n_non_radial_numax)) else: result[5] = np.nan # absolute differences on global parameters: result[6:6 + nglb] = np.absolute(model1.glb - model2.glb) return result
[docs] def find_interpolation_coefficients(grid, pt, tessellation, ndx): """ Find interpolation weights from the corresponding simplex. Linear interpolation weights are obtained with the simplex by finding the barycentric coordinates of the point given by ``pt``. :param grid: grid of models in which we're carrying out the interpolation :param pt: set of parameters used for finding the interpolation weights. The first part contains the grid parameters (relevant to this interpolation), whereas the last element is the age (not used here). If the provided set of parameters lies outside the grid, then ``None`` is returned instead of an interpolated model. :param tessellation: tessellation with which to carry out the interpolation. :param ndx: indices of the grid points associated with the tessellation :type grid: :py:class:`Model_grid` :type pt: array-like :type tessellation: Delaunay tessellation object :type ndx: list of int :return: lists of interpolation coefficients and tracks :rtype: list of floats, list of :py:class:`Track` """ if (pt is None): return None, None if (tessellation is None): return None, None if (grid.distort_mat is not None): pt1 = np.dot(np.asarray(pt[0:-1], dtype=gtype), grid.distort_mat) else: pt1 = np.asarray(pt[0:-1], dtype=gtype) val = tessellation.find_simplex(pt1.reshape((1, grid.ndim)))[0] # see if point is outside tessellation if (val == -1): return None, None mat = tessellation.transform[val] # make sure the transformation matrix is defined: if (math.isnan(np.sum(mat))): return None, None b = mat[:grid.ndim].dot(pt1 - mat[grid.ndim]) coefs = np.r_[b, 1.0 - b.sum()] ind = tessellation.simplices[val] # check to make sure you're not outside the grid: for coef in coefs: if (coef < -tol): return None, None # produce results, filtering out zero elements: coefs_out = [] tracks = [] for coef, i in zip(coefs, ind): # remove negative coefficients to avoid problems. if (coef > 0.0): coefs_out.append(coef) tracks.append(grid.tracks[ndx[i]]) return coefs_out, tracks
[docs] def find_ages(coefs, tracks, age): """ Find ages to which each track needs to be interpolated for a specified age. The variable :py:data:`age_interpolation` in AIMS_configure.py decides between the following options: 1. ``age_interpolation`` = ``age``: each track is simply interpolated to ``age``. 2. ``age_interpolation`` = ``scale_age``: the age of each model along each evolutionary track, including the interpolated track, is linearly mapped onto the interval [0,1]. A dimensionless parameter ``eta`` is obtained by interpolating ``age`` onto the interval [0,1], using the linear transformation associated with the interpolated track. Using the parameter eta, a corresponding age is obtained along each track. 3. ``age_interpolation`` = ``age_adim``: each track is interpolated to the appropriate dimensionless age so that the interpolated physical age reproduces the input age. This dimensionless age parameter is calculated via a fortran subroutine. .. figure:: ./figures/age_interpolation.* :figclass: align-center This diagram illustrates age interpolation for the two first options, namely ``age`` and ``scale_age``, and shows the advantages of selecting the latter. :param coefs: interpolation coefficients used to weight each track. :param tracks: evolutionary tracks involved in the interpolation. :param age: target age for the output interpolated model. :type coefs: list of floats :type tracks: list of :py:class:`Track` :type age: float :return: the relevant age for each track :rtype: list of floats .. note:: - the interpolation coefficients should add up to 1.0 - there should be as many tracks as interpolation coefficients. """ assert (len(coefs) == len(tracks)), "Mismatch between len(coefs) and len(tracks)" if (config.age_interpolation == "age"): return [age] * len(coefs) elif (config.age_interpolation == "scale_age"): age_s = 0.0 age_f = 0.0 for coef, track in zip(coefs, tracks): age_s += coef * track.glb[0, iage] age_f += coef * track.glb[-1, iage] eta = (age - age_s) / (age_f - age_s) # check to see if the age lies within the interpolated track: if (eta < 0.0): return None if (eta > 1.0): return None ages = [] for coef, track in zip(coefs, tracks): ages.append((1.0 - eta) * track.glb[0, iage] + eta * track.glb[-1, iage]) return ages elif (config.age_interpolation == "age_adim"): coefs = np.array(coefs) ntracks = len(tracks) sze = np.array([track.glb.shape[0] for track in tracks], dtype=int) max_sze = np.max(sze) age_adim_array = np.zeros((max_sze, ntracks), dtype=gtype) age_array = np.zeros((max_sze, ntracks), dtype=gtype) weights = np.zeros((ntracks,), dtype=gtype) indices = np.zeros((ntracks, 3), dtype=int) age_adim = np.nan # this is important to catch failures for i in range(ntracks): nmodels = sze[i] age_array[0:nmodels, i] = tracks[i].glb[:, iage] age_adim_array[0:nmodels, i] = tracks[i].glb[:, iage_adim] age_adim, indices, weights = aims_fortran.find_tau(age_adim_array, \ age_array, coefs, sze, age, age_adim, indices, weights) if (math.isnan(age_adim)): return None ages = np.zeros((ntracks,), dtype=gtype) for i in range(ntracks): ages[i] = weights[i] * age_array[indices[i, 0] - 1, i] \ + (1.0 - weights[i]) * age_array[indices[i, 0], i] return ages else: raise ValueError("Unrecognised age_interpolation option %s." % config.age_interpolation)
[docs] def interpolate_model(grid, pt, tessellation, ndx): """ Interpolate model in grid using provided parameters. The interpolation is carried out in two steps. First, linear interpolation according to age is carried out on each node of the simplex containing the set of parameters. This interpolation is done using the :py:class:`Track.interpolate_model` method. Then, linear interpolation is carried out within the simplex. This achieved by finding the barycentric coordinates of the model (i.e. the weights), before combining the age-interpolated models form the nodes using the :py:class:`combine_models` method. In this manner, the weights are only calculated once, thereby increasing computational efficiency. :param grid: grid of models in which we're carrying out the interpolation :param pt: set of parameters used for the interpolation. The first part contains the grid parameters, whereas the last element is the age. If the provided set of parameters lies outside the grid, then ``None`` is returned instead of an interpolated model. :param tessellation: tessellation with which to carry out the interpolation. :param ndx: indices of the grid points associated with the tessellation :type grid: :py:class:`Model_grid` :type pt: array-like :type tessellation: Delaunay tessellation object :type ndx: list of int :return: the interpolated model :rtype: :py:class:`Model` """ # find simplex interpolation coefficients coefs, tracks = find_interpolation_coefficients(grid, pt, tessellation, ndx) if (coefs is None): return None # find ages: ages = find_ages(coefs, tracks, pt[-1]) if (ages is None): return None n = len(tracks) # treat the case where there is only 1 model: if (n == 1): if (abs(coefs[0] - 1.0) > eps): print("WARNING: erroneous interpolation coefficient: " + str(coefs[0])) return tracks[0].interpolate_model(ages[0]) # treat the case where there are at least 2 models: aModel1 = tracks[0].interpolate_model(ages[0]) if (aModel1 is None): return None aModel2 = tracks[1].interpolate_model(ages[1]) if (aModel2 is None): return None aModel1 = combine_models(aModel1, coefs[0], aModel2, coefs[1]) for i in range(2, n): aModel2 = tracks[i].interpolate_model(ages[i]) if (aModel2 is None): return None aModel1 = combine_models(aModel1, 1.0, aModel2, coefs[i]) return aModel1
[docs] def find_combination(grid, pt): """ Find linear combination of models which corresponds to interpolating the model based on the provided parameters. The interpolation is carried out using the same procedure as in :py:func:`interpolate_model`. :param grid: grid of models in which we're carrying out the interpolation :param pt: set of parameters used for the interpolation. The first part contains the grid parameters, whereas the last element is the age. If the provided set of parameters lies outside the grid, then ``None`` is returned instead of an interpolated model. :type grid: :py:class:`Model_grid` :type pt: array-like :return: pairs of coefficients and model names :rtype: tuple of (float,string) """ # find simplex interpolation coefficients coefs, tracks = find_interpolation_coefficients(grid, pt, grid.tessellation, grid.ndx) if (coefs is None): return None # find ages: ages = find_ages(coefs, tracks, pt[-1]) if (ages is None): return None n = len(tracks) # combine multiple models: results = () for coef, track, age in zip(coefs, tracks, ages): if (coef < 0.0): # make sure we're not outside the grid return None result = track.find_combination(age, coef) if (result is None): return None results += result return results