better plotting

measure_extinction/extdata.py

@@ -770,7 +770,7 @@ def trans_elv_alav(self, av=None, akav=0.112):
         -------
         Updates self.(exts, uncs)
         """
-        if (self.type_rel_band != "V") and (self.type_rel_band != 0.55 * u.micron):
+        if (self.type_rel_band != "V") and (self.type_rel_band != 0.55 * u.micron) and (self.type_rel_band != 5500. * u.angstrom):
             warnings.warn(
                 "attempt to normalize a non-E(lambda-V) curve with A(V)", UserWarning
             )

measure_extinction/model.py

@@ -730,7 +730,15 @@ def nll(params, memodel, *args):
 
         return (outmod, result)
 
-    def fit_sampler(self, obsdata, modinfo, nsteps=1000, burnfrac=0.1, multiproc=False):
+    def fit_sampler(
+        self,
+        obsdata,
+        modinfo,
+        nsteps=1000,
+        burnfrac=0.1,
+        save_samples=None,
+        multiproc=False,
+    ):
         """
         Run a samplier (specifically emcee) to find the detailed
         parameters including uncertainties.
@@ -749,6 +757,9 @@ def fit_sampler(self, obsdata, modinfo, nsteps=1000, burnfrac=0.1, multiproc=Fal
         burnfrac : float
             fraction of nsteps to discard as the burn in [default=0.1]
 
+        save_samples : filename
+            name of hd5 file to save the MCMC samples
+
         multiproc : boolean
             set to run the emcee in parallel (does not speed up things much) [default=False]
 
@@ -770,21 +781,21 @@ def fit_sampler(self, obsdata, modinfo, nsteps=1000, burnfrac=0.1, multiproc=Fal
         if not np.isfinite(outmod.lnprior()):
             raise ValueError("ln(prior) is not finite")
 
-        # simple function to turn the log(likelihood) into the chisqr
-        #  required as op.minimize function searches for the minimum chisqr (not max likelihood like MCMC algorithms)
-        # def lnprob(params, memodel, *args):
-        #    memodel.fit_to_parameters(params)
-        #    return memodel.lnprob(*args)
-
         # get the non-fixed initial parameters
         p0 = outmod.parameters_to_fit()
 
         # setup the sampliers
         ndim = len(p0)
         nwalkers = 2 * ndim
+
         # setting up the walkers to start "near" the inital guess
         p = [p0 * (1 + 0.01 * np.random.normal(0, 1.0, ndim)) for k in range(nwalkers)]
 
+        if save_samples:
+            # Don't forget to clear it in case the file already exists
+            save_backend = emcee.backends.HDFBackend(save_samples)
+            save_backend.reset(nwalkers, ndim)
+
         # setup and run the sampler
         if multiproc:
             with Pool() as pool:
@@ -802,6 +813,7 @@ def fit_sampler(self, obsdata, modinfo, nsteps=1000, burnfrac=0.1, multiproc=Fal
                 ndim,
                 _lnprob,
                 args=(outmod, obsdata, modinfo),
+                backend=save_backend,
             )
             sampler.run_mcmc(p, nsteps, progress=True)
 

measure_extinction/plotting/plot_ext.py

@@ -414,7 +414,7 @@ def plot_multi_extinction(
     fitmodel=False,
     HI_lines=False,
     range=None,
-    spread=False,
+    spread=0.0,
     exclude=[],
     log=False,
     text_offsets=[],
@@ -454,9 +454,9 @@ def plot_multi_extinction(
     range : list of 2 floats [default=None]
         Wavelength range to be plotted (in micron) - [min,max]
 
-    spread : boolean [default=False]
-        Whether or not to spread the extinction curves out by adding a vertical offset to each curve
-
+    spread : float [default=0]
+        Amount to addiatively spread the curves 
+ 
     exclude : list of strings [default=[]]
         List of data type(s) to exclude from the plot (e.g., "IRS", "IRAC1")
 
@@ -515,12 +515,17 @@ def plot_multi_extinction(
         text_angles = np.full(len(starpair_list), 10)
 
     for i, starpair in enumerate(starpair_list):
+        if ".fits" not in starpair:
+            fname = "%s%s_ext.fits" % (path, starpair.lower())
+        else:
+            fname = starpair 
+
         # read in the extinction curve data
-        extdata = ExtData("%s%s_ext.fits" % (path, starpair.lower()))
+        extdata = ExtData(fname)
 
         # spread out the curves if requested
-        if spread:
-            yoffset = 0.25 * i
+        if spread != 0.0:
+            yoffset = spread * i
         else:
             yoffset = 0.0
 
@@ -795,14 +800,12 @@ def main():
     # commandline parser
     parser = argparse.ArgumentParser()
     parser.add_argument(
-        "starpair_list",
-        nargs="+",
-        help='pairs of star names for which to plot the extinction curve, in the format "reddenedstarname_comparisonstarname", without spaces',
+        "starpair_list", nargs="+", help="filenames of extinction curves"
     )
     parser.add_argument(
         "--path",
         help="path to the data files",
-        default=get_datapath(),
+        default="./",
     )
     parser.add_argument("--alax", help="plot A(lambda)/A(X)", action="store_true")
     parser.add_argument(
@@ -827,8 +830,9 @@ def main():
     )
     parser.add_argument(
         "--spread",
-        help="spread the curves out over the figure; can only be used in combination with --onefig",
-        action="store_true",
+        help="spread the curves out over the figure by this amount; can only be used in combination with --onefig",
+        default=0.0,
+        type=float,
     )
     parser.add_argument(
         "--exclude",