Skip to content

Commit

Permalink
Included optional scaling parameters in FitModel
Browse files Browse the repository at this point in the history
  • Loading branch information
Tomas Stolker committed Feb 22, 2020
1 parent ba7e695 commit 8611aac
Show file tree
Hide file tree
Showing 10 changed files with 198 additions and 83 deletions.
57 changes: 48 additions & 9 deletions species/analysis/fit_model.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -70,10 +70,10 @@ def lnprior(param,

def lnlike(param,
modelpar,
modelphot,
objphot,
distance,
spectrum,
modelphot,
modelspec):
"""
Internal function for the likelihood probability.
Expand All @@ -84,15 +84,18 @@ def lnlike(param,
Parameter values.
modelpar : list(str, )
Parameter names.
modelphot : list(species.read.read_model.ReadModel, )
objphot : list(tuple(float, float), )
List with the fluxes and uncertainties of the object.
distance : float
Distance (pc).
spectrum : dict
Dictionary with the spectrum stored as wavelength (micron), flux (W m-2 micron-1),
and error (W m-2 micron-1), and optionally the covariance matrix and the inverse of
the covariance matrix.
modelphot : list(species.read.read_model.ReadModel, )
List with the interpolated synthetic photometry.
modelspec : list(species.read.read_model.ReadModel, )
List with the interpolated synthetic spectra.
Returns
-------
Expand All @@ -101,9 +104,13 @@ def lnlike(param,
"""

paramdict = {}
spec_scaling = {}

for i, item in enumerate(modelpar):
if item == 'radius':
radius = param[i]
elif item in spectrum:
spec_scaling[item] = param[i]
else:
paramdict[item] = param[i]

Expand All @@ -118,7 +125,12 @@ def lnlike(param,

if spectrum is not None:
for i, item in enumerate(spectrum.keys()):
flux = scaling * modelspec[i].spectrum_interp(list(paramdict.values()))[0, :]
if item in spec_scaling:
scaling_new = scaling * spec_scaling[item]
flux = scaling_new * modelspec[i].spectrum_interp(list(paramdict.values()))[0, :]

else:
flux = scaling * modelspec[i].spectrum_interp(list(paramdict.values()))[0, :]

if spectrum[item][2] is not None:
spec_diff = spectrum[item][0][:, 1] - flux
Expand All @@ -134,11 +146,11 @@ def lnlike(param,
def lnprob(param,
bounds,
modelpar,
modelphot,
objphot,
distance,
prior,
spectrum,
modelphot,
modelspec):
"""
Internal function for the posterior probability.
Expand All @@ -151,8 +163,8 @@ def lnprob(param,
Parameter boundaries.
modelpar : list(str, )
Parameter names.
modelphot : list('species.read.read_model.ReadModel, )
objphot : list(tuple(float, float), )
List with the fluxes and uncertainties of the object.
distance : float
Distance (pc).
prior : tuple(str, float, float)
Expand All @@ -161,7 +173,10 @@ def lnprob(param,
used if set to None.
spectrum : dict
Wavelength (micron), apparent flux (W m-2 micron-1), and flux error (W m-2 micron-1).
modelphot : list(species.read.read_model.ReadModel, )
List with the interpolated synthetic fluxes.
modelspec : list(species.read.read_model.ReadModel, )
List with the interpolated synthetic spectra.
Returns
-------
Expand All @@ -177,10 +192,10 @@ def lnprob(param,
else:
ln_prob = ln_prior + lnlike(param,
modelpar,
modelphot,
objphot,
distance,
spectrum,
modelphot,
modelspec)

if np.isnan(ln_prob):
Expand Down Expand Up @@ -224,7 +239,10 @@ def __init__(self,
case, the radius range is set to 0-5 Rjup. It is also possible to specify the bounds
for a subset of the parameters, for example, ``{'radius': (0.5, 10.)}``. Restricting
the boundaries will decrease the computation time with the interpolation prior to the
MCMC sampling.
MCMC sampling. An additional scaling parameter can be fitted for each spectrum in which
case, the boundaries have to be provided with the database tag of the spectrum.
For example, ``{'sphere_ifs': (0.5, 2.)}`` if the spectrum is stored as ``sphere_ifs``
with :func:`~species.data.database.Database.add_object`.
inc_phot : bool
Include photometric data in the fit.
inc_spec : bool
Expand Down Expand Up @@ -314,6 +332,11 @@ def __init__(self,
self.modelpar = readmodel.get_parameters()
self.modelpar.append('radius')

if self.spectrum is not None:
for item in self.spectrum.keys():
if item in self.bounds:
self.modelpar.append(item)

def run_mcmc(self,
nwalkers,
nsteps,
Expand Down Expand Up @@ -347,6 +370,11 @@ def run_mcmc(self,

sigma = {'teff': 5., 'logg': 0.01, 'feh': 0.01, 'fsed': 0.01, 'co': 0.01, 'radius': 0.01}

if self.spectrum is not None:
for item in self.spectrum.keys():
if item in self.bounds:
sigma[item] = 0.01

print('Running MCMC...')

ndim = len(self.bounds)
Expand All @@ -367,19 +395,30 @@ def run_mcmc(self,
lnprob,
args=([self.bounds,
self.modelpar,
self.modelphot,
self.objphot,
self.distance[0],
prior,
self.spectrum,
self.modelphot,
self.modelspec]))

sampler.run_mcmc(initial, nsteps, progress=True)

species_db = database.Database()

spec_labels = []

if self.spectrum is not None:
for item in self.spectrum.keys():
if item in self.bounds:
spec_labels.append(item)

else:
spec_labels = None

species_db.add_samples(sampler=sampler,
spectrum=('model', self.model),
tag=tag,
modelpar=self.modelpar,
distance=self.distance[0])
distance=self.distance[0],
spec_labels=spec_labels)
37 changes: 26 additions & 11 deletions species/data/database.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -727,7 +727,8 @@ def add_samples(self,
spectrum,
tag,
modelpar,
distance=None):
distance=None,
spec_labels=None):
"""
Parameters
----------
Expand All @@ -742,6 +743,9 @@ def add_samples(self,
List with the model parameter names.
distance : float
Distance to the object (pc). Not used if set to None.
spec_labels : list(str, )
List with the spectrum labels that are used for fitting an additional scaling
parameter.
Returns
-------
Expand Down Expand Up @@ -773,9 +777,17 @@ def add_samples(self,
if distance:
dset.attrs['distance'] = float(distance)

count_scaling = 0

for i, item in enumerate(modelpar):
dset.attrs[f'parameter{i}'] = str(item)

if spec_labels is not None and item in spec_labels:
dset.attrs[f'scaling{count_scaling}'] = str(item)
count_scaling += 1

dset.attrs['nscaling'] = int(count_scaling)

mean_accep = np.mean(sampler.acceptance_fraction)
dset.attrs['acceptance'] = float(mean_accep)
print(f'Mean acceptance fraction: {mean_accep:.3f}')
Expand Down Expand Up @@ -831,23 +843,17 @@ def get_probable_sample(self,
# max_prob = probability[index_max]
max_sample = samples[index_max]

# print(f'Maximum probability: {max_prob:.2f}')
# print(f'Most probable sample:', end='', flush=True)

prob_sample = {}

for i in range(nparam):
par_key = dset.attrs[f'parameter{i}']
par_value = max_sample[i]

prob_sample[par_key] = par_value
# print(f' {par_key}={par_value:.2f}', end=')

if dset.attrs.__contains__('distance'):
prob_sample['distance'] = dset.attrs['distance']

# print()

h5_file.close()

return prob_sample
Expand Down Expand Up @@ -875,6 +881,11 @@ def get_median_sample(self,
dset = h5_file[f'results/mcmc/{tag}/samples']

nparam = dset.attrs['nparam']
nscaling = dset.attrs['nscaling']

scaling = []
for i in range(nscaling):
scaling.append(dset.attrs[f'scaling{i}'])

samples = np.asarray(dset)
samples = samples[:, burnin:, :]
Expand All @@ -885,7 +896,6 @@ def get_median_sample(self,
for i in range(nparam):
par_key = dset.attrs[f'parameter{i}']
par_value = np.percentile(samples[:, i], 50.)

median_sample[par_key] = par_value

if dset.attrs.__contains__('distance'):
Expand Down Expand Up @@ -923,12 +933,12 @@ def get_mcmc_spectra(self,
Boxes with the randomly sampled spectra.
"""

print('The smooth_spectrum function has not been fully tested. [WARNING]')

h5_file = h5py.File(self.database, 'r')
dset = h5_file[f'results/mcmc/{tag}/samples']

nparam = dset.attrs['nparam']
nscaling = dset.attrs['nscaling']

spectrum_type = dset.attrs['type']
spectrum_name = dset.attrs['spectrum']

Expand All @@ -952,6 +962,10 @@ def get_mcmc_spectra(self,
for i in range(nparam):
param.append(dset.attrs[f'parameter{i}'])

scaling = []
for i in range(nscaling):
scaling.append(dset.attrs[f'scaling{i}'])

if spectrum_type == 'model':
if spectrum_name == 'planck':
readmodel = read_planck.ReadPlanck(wavel_range)
Expand All @@ -966,7 +980,8 @@ def get_mcmc_spectra(self,
for i in tqdm.tqdm(range(samples.shape[0]), desc='Getting MCMC spectra'):
model_param = {}
for j in range(samples.shape[1]):
model_param[param[j]] = samples[i, j]
if param[j] not in scaling:
model_param[param[j]] = samples[i, j]

if distance:
model_param['distance'] = distance
Expand Down
2 changes: 1 addition & 1 deletion species/data/petitcode.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -444,7 +444,7 @@ def add_petitcode_hot_cloudy(input_path,
logg_val = float(filename[19:23])
feh_val = float(filename[28:32])
co_ratio_val = float(filename[36:40])
fsed_vale = float(filename[46:50])
fsed_val = float(filename[46:50])

if teff_range is not None:
if teff_val < teff_range[0] or teff_val > teff_range[1]:
Expand Down
31 changes: 26 additions & 5 deletions species/plot/plot_mcmc.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,8 @@
import matplotlib as mpl
import matplotlib.pyplot as plt

from matplotlib.ticker import ScalarFormatter

from species.data import database
from species.util import plot_util

Expand Down Expand Up @@ -144,11 +146,15 @@ def plot_posterior(tag,
for key, value in box.median_sample.items():
print(f' - {key} = {value:.2f}')

print(f'Plotting the posterior: {output}...', end='', flush=True)

samples = box.samples
par_val = tuple(box.prob_sample.values())

print(f'Maximum posterior sample:')
for key, value in box.prob_sample.items():
print(f' - {key} = {value:.2f}')

print(f'Plotting the posterior: {output}...', end='', flush=True)

labels = plot_util.update_labels(box.parameters)

ndim = samples.shape[-1]
Expand All @@ -165,13 +171,28 @@ def plot_posterior(tag,
if i >= j:
ax = axes[i, j]

ax.xaxis.set_major_formatter(ScalarFormatter(useOffset=False))
ax.yaxis.set_major_formatter(ScalarFormatter(useOffset=False))

if j == 0 and i != 0:
labelleft = True
else:
labelleft = False

if i == ndim-1:
labelbottom = True
else:
labelbottom = False

ax.tick_params(axis='both', which='major', colors='black', labelcolor='black',
direction='in', width=1, length=5, labelsize=12, top=True,
bottom=True, left=True, right=True)
bottom=True, left=True, right=True, labelleft=labelleft,
labelbottom=labelbottom, labelright=False, labeltop=False)

ax.tick_params(axis='both', which='minor', colors='black', labelcolor='black',
direction='in', width=1, length=3, labelsize=12, top=True,
bottom=True, left=True, right=True)
bottom=True, left=True, right=True, labelleft=labelleft,
labelbottom=labelbottom, labelright=False, labeltop=False)

if limits is not None:
ax.set_xlim(limits[j])
Expand Down Expand Up @@ -256,7 +277,7 @@ def plot_photometry(tag,

ax.get_xaxis().set_label_coords(0.5, -0.26)

plt.savefig(f'{os.getcwd()}/{output}', bbox_inches='tight')
plt.savefig(os.getcwd()+'/'+output, bbox_inches='tight')
plt.clf()
plt.close()

Expand Down
Loading

0 comments on commit 8611aac

Please sign in to comment.