a lot of additions

tardisnuclear/ejecta.py

@@ -1,3 +1,5 @@
+from collections import OrderedDict
+
 from astropy import units as u
 import pandas as pd
 from pyne.material import Material
@@ -60,6 +62,10 @@ def from_masses(cls, **kwargs):
     def __init__(self, mass_msol, composition):
         self.mass_g = mass_msol * msun_to_cgs
         self.material = Material(self._normalize_composition(composition))
+        self._pad_material()
+        atomic_masses = self.get_masses()
+        self.n_per_g = [1 / atomic_masses[item]
+                            for item in self.get_all_children_nuc_name()]
 
 
     @property
@@ -82,9 +88,9 @@ def isotopes(self):
         return [nucname.name(id) for id in self.keys()]
 
     def get_decay_const(self):
-        return {nuc_name:data.decay_const(nuc_id)
+        return OrderedDict((nuc_name, data.decay_const(nuc_id))
                 for nuc_id, nuc_name in zip(self.get_all_children(),
-                                            self.get_all_children_nuc_name())}
+                                            self.get_all_children_nuc_name()))
 
     def get_half_life(self):
         return [data.half_life(nuc_id) for nuc_id in self.keys()]
@@ -106,7 +112,7 @@ def get_child(nuc_id):
             children_set.add(nuc_id)
             get_child(nuc_id)
 
-        return children_set
+        return sorted(children_set)
 
     def get_all_children_nuc_name(self):
         return [nucname.name(nuc_id) for nuc_id in self.get_all_children()]
@@ -119,11 +125,41 @@ def _normalize_composition(composition):
                               for key, value in composition.items()}
         return normed_composition
 
+    def _pad_material(self):
+        for isotope in self.get_all_children_nuc_name():
+            try:
+                self.material[isotope]
+            except KeyError:
+                self.material[isotope] = 0.0
+
 
     def decay(self, time):
         new_material= self.material.decay(time.to(u.s).value)
         return self.__class__(self.mass_g / msun_to_cgs, new_material)
 
+    def decay_epochs(self, epochs):
+        epochs = u.Quantity(epochs, u.day)
+        isotope_children = self.get_all_children()
+        columns = self.get_all_children_nuc_name()
+        decayed_fractions = pd.DataFrame(index=epochs.value, columns=columns,
+                                       dtype=np.float64)
+        day_to_s = 24 * 3600
+        for epoch, row in decayed_fractions.iterrows():
+            new_material = self.material.decay(epoch * day_to_s)
+            fractions = [0.0 if key not in new_material else new_material[key]
+                         for key in isotope_children]
+            row.values[:] = fractions
+        return decayed_fractions
+
+    def get_decayed_numbers(self, epochs):
+        epochs = u.Quantity(epochs, u.day)
+
+        fractions = self.decay_epochs(epochs)
+
+        return fractions * self.mass_g * self.n_per_g
+
+
+
     def __repr__(self):
         return self.material.__str__()
 

tardisnuclear/io/read_nndc.py

@@ -136,6 +136,10 @@ def _sanititze_table(df):
         if 'energy' in df.columns:
             df.energy = df.energy.apply(lambda x: u.Quantity(
                 float(x.split()[0]), u.keV).to(u.erg).value)
+        if 'end_point_energy' in df.columns:
+            df.end_point_energy = df.end_point_energy.apply(lambda x: u.Quantity(
+                float(x.split()[0]), u.keV).to(u.erg).value)
+
         if 'intensity' in df.columns:
             df.intensity = df.intensity.apply(lambda x:
                                               float(x.split('%')[0])/100.)
@@ -164,13 +168,13 @@ class ElectronTableParser(BaseParser):
 class BetaPlusTableParser(BaseParser):
     html_name = 'Beta+'
     name = 'beta_plus'
-    columns = ['type', 'end_point_energy', 'energy', 'intensity', 'dose']
+    columns = ['energy', 'end_point_energy', 'intensity', 'dose']
 
 
 class BetaMinusTableParser(BaseParser):
     html_name = 'Beta-'
     name = 'beta_minus'
-    columns = ['type', 'end_point_energy', 'energy', 'intensity', 'dose']
+    columns = ['energy', 'end_point_energy', 'intensity', 'dose']
 
 
 

tardisnuclear/multinest/__init__.py

@@ -0,0 +1 @@
+

tardisnuclear/multinest/fitting.py

@@ -0,0 +1,209 @@
+import numpy as np
+from scipy import optimize
+import sys
+
+from astropy import modeling
+from itertools import chain
+from tardisnuclear.ejecta import Ejecta
+from tardisnuclear.rad_trans import SimpleLateTime
+from tardisnuclear.nuclear_data import NuclearData
+
+from scipy import stats
+from collections import OrderedDict
+import pandas as pd
+
+from astropy import units as u
+import pymultinest
+
+msun_to_cgs = u.Msun.to(u.g)
+mpc_to_cm = u.Mpc.to(u.cm)
+
+
+class BaseModel(modeling.Model):
+
+    def __call__(self, *inputs, **kwargs):
+        parameters = self._param_sets(raw=True)
+        return self.evaluate(*chain(inputs, parameters))
+
+
+class BolometricLightCurveModel(BaseModel):
+    pass
+
+class BolometricLightCurveModelIa(object):
+
+
+
+    def __init__(self, epochs, lum_dens, lum_dens_err, ni56, ni57, co55, ti44):
+        self.epochs = epochs
+        self.lum_dens = lum_dens
+        self.lum_dens_err = lum_dens_err
+        self.ejecta = Ejecta.from_masses(Ni56=ni56 * u.Msun, Ni57=ni57 * u.Msun,
+                                         Co55=co55 * u.Msun, Ti44=ti44 * u.Msun)
+        self.nuclear_data = NuclearData(self.ejecta.get_all_children_nuc_name())
+        self.rad_trans = SimpleLateTime(self.ejecta, self.nuclear_data)
+
+    def calculate_light_curve(self, ni56, ni57, co55, ti44, fraction=1.0,
+                              distance=6.4, epochs=None):
+
+        if epochs is None:
+            epochs = self.epochs
+        total_mass = ni56 + ni57 + co55 + ti44
+
+        self.ejecta.mass_g = total_mass * msun_to_cgs
+        self.ejecta['Ni56'] = ni56 / total_mass
+        self.ejecta['Ni57'] = ni57 / total_mass
+        self.ejecta['Co55'] = co55 / total_mass
+        self.ejecta['Ti44'] = ti44 / total_mass
+        luminosity_density =  self.rad_trans.total_bolometric_light_curve(epochs)
+        return (luminosity_density * fraction /
+                (4 * np.pi * (distance * mpc_to_cm)**2))
+
+
+    def calculate_individual_light_curve(self, ni56, ni57, co55, ti44, fraction=1.0,
+                              distance=6.4, epochs=None):
+
+        if epochs is None:
+            epochs = self.epochs
+        total_mass = ni56 + ni57 + co55 + ti44
+
+        self.ejecta.mass_g = total_mass * msun_to_cgs
+        self.ejecta['Ni56'] = ni56 / total_mass
+        self.ejecta['Ni57'] = ni57 / total_mass
+        self.ejecta['Co55'] = co55 / total_mass
+        self.ejecta['Ti44'] = ti44 / total_mass
+        luminosity_density =  self.rad_trans.bolometric_light_curve(epochs)
+        return (luminosity_density * fraction /
+                (4 * np.pi * (distance * mpc_to_cm)**2))
+
+
+    def fitness_function(self, ni56, ni57, co55, ti44, fraction, distance):
+
+        model_light_curve = self.calculate_light_curve(ni56, ni57, co55, ti44,
+                                                 fraction, distance)
+        return (model_light_curve.value - self.lum_dens)/self.lum_dens_err
+
+
+    def log_likelihood(self, model_param, ndim, nparam):
+        #return -5
+
+        model_param = [model_param[i] for i in xrange(6)]
+        return (-0.5 * self.fitness_function(*model_param)**2).sum()
+
+    def simple_fit(self, ni56, ni57, co55, ti44, method='Nelder-Mead'):
+        def fit_func(isotopes):
+            ni57, co55, ti44 = np.abs(isotopes)
+            mdl = self.evaluate(ni56, ni57, co55, ti44)
+            mdl *= np.mean(self.luminosity / mdl.value)
+            return ((mdl.value - self.luminosity)**2).sum()
+
+        fit = optimize.minimize(fit_func, (ni57, co55, ti44),
+                                 method=method)
+        mdl = self.evaluate(ni56, *fit.x)
+        norm_factor = np.mean(self.luminosity / mdl.value)
+        mdl *= norm_factor
+
+        return fit, norm_factor, mdl
+
+
+    def multinest_fit(self, priors, **kwargs):
+
+
+        mn_fit = pymultinest.run(self.log_likelihood, priors.prior_transform, 6,
+                                 outputfiles_basename='sn11fe/fit', **kwargs)
+
+        return mn_fit
+
+
+
+
+class MultiNestResult():
+
+
+    @classmethod
+    def from_multinest_basename(cls, basename, parameter_names):
+        """
+        Reading a MultiNest result from a basename
+
+        Parameters
+        ----------
+
+        basename: str
+            basename (path + prefix) for a multinest run
+
+        Returns
+            : ~MultinestResult
+        """
+
+        posterior_data = cls.read_posterior_data(basename, parameter_names)
+
+        return cls(posterior_data)
+
+    @classmethod
+    def from_hdf5(cls, h5_fname, key):
+        """
+        Reading a Multinest result from its generated HDF5 file
+
+        Parameters
+        ----------
+
+        h5_fname: ~str
+            HDF5 filename
+
+        key: ~str
+            group identifier in the store
+        """
+
+        posterior_data = pd.read_hdf(h5_fname, key)
+
+        return cls(posterior_data)
+
+
+    @staticmethod
+    def read_posterior_data(basename, parameter_names):
+        """
+        Reading the posterior data into a pandas dataframe
+
+        """
+        posterior_data = pd.read_csv('{0}/fit.txt'.format(basename),
+                           delim_whitespace=True,
+                           names=['posterior', 'x'] + parameter_names)
+        posterior_data.index = np.arange(len(posterior_data))
+        return posterior_data
+
+    def __init__(self, posterior_data):
+        self.posterior_data = posterior_data
+        self.parameter_names = [col_name for col_name in posterior_data.columns
+                                if col_name not in ['x', 'posterior']]
+
+    def calculate_sigmas(self, sigma):
+        sigmas = OrderedDict()
+        for parameter_name in self.parameter_names:
+            posterior_data = self.posterior_data.sort(parameter_name)
+            parameter_values, posterior_values = (posterior_data[parameter_name],
+                                                  posterior_data['posterior'])
+            posterior_cumsum = posterior_values.cumsum()
+
+            norm_distr = stats.norm(loc=0.0, scale=1.)
+
+            sigma_low = np.interp(norm_distr.cdf(-sigma), posterior_cumsum,
+                                  parameter_values)
+
+            sigma_high = np.interp(norm_distr.cdf(sigma), posterior_cumsum,
+                                  parameter_values)
+
+
+            sigmas[parameter_name] = (sigma_low, sigma_high)
+
+        return sigmas
+
+    @property
+    def mean(self):
+        if not hasattr(self, '_mean'):
+            _mean = OrderedDict([(param_name,
+                                  np.average(self.posterior_data[param_name],
+                                             weights=
+                                             self.posterior_data['posterior']))
+                                 for param_name in self.parameter_names])
+            self._mean = _mean
+
+        return self._mean