calibration spectra

species/__init__.py

@@ -1,8 +1,12 @@
-from species.analysis.fit import FitSpectrum
+from species.analysis.fit_model import FitModel
+
+from species.analysis.fit_spectrum import FitSpectrum
 
 from species.analysis.photometry import SyntheticPhotometry, \
                                         apparent_to_absolute
 
+from species.read.read_calibration import ReadCalibration
+
 from species.read.read_filter import ReadFilter
 
 from species.read.read_model import ReadModel, \

species/analysis/fit.py → species/analysis/fit_model.py

@@ -1,13 +1,10 @@
-"""
+'''
 Text
-"""
+'''
 
-import os
 import sys
 import math
-import configparser
 
-import h5py
 import emcee
 import progress.bar
 import numpy as np
@@ -25,7 +22,7 @@ def lnprior(param,
             bounds,
             modelpar,
             prior):
-    """
+    '''
     :param param: Parameter values.
     :type param: numpy.ndarray
     :param bounds: Parameter boundaries.
@@ -39,7 +36,7 @@ def lnprior(param,
 
     :return: Log prior probability.
     :rtype: float
-    """
+    '''
 
     if prior:
 
@@ -72,7 +69,7 @@ def lnlike(param,
            synphot,
            sampling,
            distance):
-    """
+    '''
     :param param:
     :type param:
     :param modelpar:
@@ -90,7 +87,7 @@ def lnlike(param,
 
     :return: Log likelihood probability.
     :rtype: float
-    """
+    '''
 
     global MIN_CHISQ
     global MIN_PARAM
@@ -123,7 +120,7 @@ def lnprob(param,
            sampling,
            distance,
            prior):
-    """
+    '''
     :param param:
     :type param:
     :param bounds:
@@ -145,7 +142,7 @@ def lnprob(param,
 
     :return:
     :rtype:
-    """
+    '''
 
     ln_prior = lnprior(param, bounds, modelpar, prior)
 
@@ -164,18 +161,18 @@ def lnprob(param,
     return ln_prob
 
 
-class FitSpectrum:
-    """
+class FitModel:
+    '''
     Text
-    """
+    '''
 
     def __init__(self,
                  objname,
                  filters,
                  model,
                  sampling,
                  bounds):
-        """
+        '''
         :param objname: Object name in the database.
         :type objname: str
         :param filters: Filter IDs for which the photometry is selected. All available
@@ -184,17 +181,14 @@ def __init__(self,
         :name model: Atmospheric model.
         :type model: str
         :name sampling: Wavelength sampling for the computation of synthetic photometry
-                        ("specres" or "gaussian").
+                        ('specres' or 'gaussian').
         :type sampling: tuple
         :name bounds: Parameter boundaries. Full parameter range is used if None or not specified.
                       The radius parameter range is set to 0-5 Rjup if not specified.
         :type bounds: dict
 
         :return: None
-        """
-
-        self.parsec = 3.08567758147e16 # [m]
-        self.r_jup = 71492000. # [m]
+        '''
 
         self.object = read_object.ReadObject(objname)
         self.distance = self.object.get_distance()
@@ -252,10 +246,10 @@ def run_mcmc(self,
                  guess,
                  tag,
                  prior=None,
-                 ncpu=1,):
-        """
+                 ncpu=1):
+        '''
         :return: None
-        """
+        '''
 
         global MIN_CHISQ
         global MIN_PARAM
@@ -300,75 +294,11 @@ def run_mcmc(self,
 
         progbar.finish()
 
-        self.store_samples(sampler, self.model, tag, (MIN_CHISQ, MIN_PARAM))
-
-    def store_samples(self,
-                      sampler,
-                      model,
-                      tag,
-                      chisquare):
-        """
-        :param sampler: Ensemble sampler.
-        :type sampler: emcee.ensemble.EnsembleSampler
-        :param model: Atmospheric model.
-        :type model: str
-        :param tag: Database tag.
-        :type tag: str
-        :param chisquare: Maximum likelihood solution. Tuple with the chi-square value and related
-                          parameter values.
-        :type chisquare: tuple(float, float)
-
-        :return: None
-        """
-
-        config_file = os.path.join(os.getcwd(), 'species_config.ini')
-
-        config = configparser.ConfigParser()
-        config.read_file(open(config_file))
-
-        species_db = config['species']['database']
-
-        h5_file = h5py.File(species_db, 'a')
-
-        if 'results' not in h5_file:
-            h5_file.create_group('results')
-
-        if 'results/mcmc' not in h5_file:
-            h5_file.create_group('results/mcmc')
-
-        if 'results/mcmc/'+tag in h5_file:
-            del h5_file['results/mcmc/'+tag]
-
-        samples = sampler.chain
-
-        dset = h5_file.create_dataset('results/mcmc/'+tag,
-                                      data=samples,
-                                      dtype='f')
-
-        dset.attrs['model'] = str(model)
-        dset.attrs['distance'] = float(self.distance)
-        dset.attrs['nparam'] = int(len(self.modelpar))
-
-        for i, item in enumerate(self.modelpar):
-            dset.attrs['parameter'+str(i)] = str(item)
-
-        dset.attrs['min_chi'] = float(chisquare[0])
-        for i, item in enumerate(self.modelpar):
-            dset.attrs['chisquare'+str(i)] = float(chisquare[1][item])
-
-        mean_accep = np.mean(sampler.acceptance_fraction)
-        dset.attrs['acceptance'] = float(mean_accep)
-        print('Mean acceptance fraction: {0:.3f}'.format(mean_accep))
-
-        try:
-            int_auto = emcee.autocorr.integrated_time(sampler.flatchain)
-            print('Integrated autocorrelation time =', int_auto)
-
-        except emcee.autocorr.AutocorrError:
-            int_auto = None
-
-        if int_auto is not None:
-            for i, item in enumerate(int_auto):
-                dset.attrs['autocorrelation'+str(i)] = float(item)
+        species_db = database.Database()
 
-        h5_file.close()
+        species_db.add_samples(sampler=sampler,
+                               spectrum=('model', self.model),
+                               tag=tag,
+                               chisquare=(MIN_CHISQ, MIN_PARAM),
+                               modelpar=self.modelpar,
+                               distance=self.distance)

species/analysis/fit_spectrum.py

@@ -0,0 +1,171 @@
+'''
+Text
+'''
+
+import os
+import sys
+import math
+import configparser
+
+import h5py
+import emcee
+import progress.bar
+import numpy as np
+
+from species.analysis import photometry
+from species.data import database
+from species.read import read_model, read_object, read_calibration
+
+
+MIN_CHISQ = np.inf
+MIN_PARAM = None
+
+
+def lnprob(param,
+           bounds,
+           objphot,
+           specphot):
+    '''
+    :param param: Parameter values.
+    :type param: numpy.ndarray
+    :param bounds: Parameter boundaries.
+    :type bounds: tuple(float, float)
+    :param objphot:
+    :type objphot:
+    :param specphot:
+    :type specphot:
+
+    :return:
+    :rtype:
+    '''
+
+    global MIN_CHISQ
+    global MIN_PARAM
+
+    if bounds[0] <= param <= bounds[1]:
+        ln_prior = 0.
+
+    else:
+        ln_prior = -np.inf
+
+    if math.isinf(ln_prior):
+        ln_prob = -np.inf
+
+    else:
+        chisq = 0.
+        for i, item in enumerate(objphot):
+            chisq += (objphot[i][0]-param*specphot[0])**2 / objphot[i][1]**2
+
+        if chisq < MIN_CHISQ:
+            MIN_CHISQ = chisq
+            MIN_PARAM = {'scaling':param}
+
+        ln_prob = ln_prior - 0.5*chisq
+
+    return ln_prob
+
+
+class FitSpectrum:
+    '''
+    Text
+    '''
+
+    def __init__(self,
+                 objname,
+                 filters,
+                 spectrum,
+                 bounds):
+        '''
+        :param objname: Object name in the database.
+        :type objname: str
+        :param filters: Filter IDs for which the photometry is selected. All available
+                        photometry of the object is selected if set to None.
+        :type filters: tuple(str, )
+        :param spectrum: Calibration spectrum.
+        :type spectrum: str
+        :param bounds: Range of the scaling parameter (min, max).
+        :type bounds: tuple(float, float)
+
+        :return: None
+        '''
+
+        self.object = read_object.ReadObject(objname)
+
+        self.spectrum = spectrum
+        self.bounds = bounds
+
+        self.objphot = []
+        self.specphot = []
+
+        if filters is None:
+            species_db = database.Database()
+            objectbox = species_db.get_object(objname, None)
+            filters = objectbox.filter
+
+        for item in filters:
+            readcalib = read_calibration.ReadCalibration(self.spectrum, item)
+            calibspec = readcalib.get_spectrum()
+
+            synphot = photometry.SyntheticPhotometry(item)
+            spec_phot = synphot.spectrum_to_photometry(calibspec.wavelength, calibspec.flux)
+            self.specphot.append(spec_phot)
+
+            obj_phot = self.object.get_photometry(item)
+            self.objphot.append((obj_phot[2], obj_phot[3]))
+
+        self.modelpar = ['scaling']
+
+    def run_mcmc(self,
+                 nwalkers,
+                 nsteps,
+                 guess,
+                 tag):
+        '''
+        :param nwalkers: Number of walkers.
+        :type nwalkers: int
+        :param nsteps: Number of steps for each walker.
+        :type nsteps: int
+        :param guess: Guess of the scaling factor.
+        :type guess: float
+        :param tag: Database tag for the results.
+        :type tag: int
+
+        :return: None
+        '''
+
+        global MIN_CHISQ
+        global MIN_PARAM
+
+        sys.stdout.write('Running MCMC...')
+        sys.stdout.flush()
+
+        ndim = 1
+
+        initial = np.zeros((nwalkers, ndim))
+        initial[:, 0] = guess + np.random.normal(0, 1e-3*guess, nwalkers)
+
+        sampler = emcee.EnsembleSampler(nwalkers=nwalkers,
+                                        dim=ndim,
+                                        lnpostfn=lnprob,
+                                        a=2.,
+                                        args=([self.bounds,
+                                               self.objphot,
+                                               self.specphot]))
+
+        progbar = progress.bar.Bar('\rRunning MCMC...',
+                                   max=nsteps,
+                                   suffix='%(percent)d%%')
+
+        for i, _ in enumerate(sampler.sample(initial, iterations=nsteps)):
+            progbar.next()
+
+        progbar.finish()
+
+        species_db = database.Database()
+
+        species_db.add_samples(sampler=sampler,
+                               spectrum=('calibration', self.spectrum),
+                               tag=tag,
+                               chisquare=(MIN_CHISQ, MIN_PARAM),
+                               modelpar=self.modelpar,
+                               distance=None)