Merge f3788ff into f3156e4

dolphin/analysis/output.py

@@ -204,7 +204,8 @@ def plot_model_overview(self, lens_name, model_id=None,
                             data_cmap='cubehelix', residual_cmap='RdBu_r',
                             convergence_cmap='afmhot',
                             magnification_cmap='viridis',
-                            v_min=None, v_max=None, print_results=False):
+                            v_min=None, v_max=None, print_results=False,
+                            show_source_light=False):
         """
         Plot the model, residual, reconstructed source, convergence,
         and magnification profiles. Either `model_id` or `kwargs_result`
@@ -270,8 +271,16 @@ def plot_model_overview(self, lens_name, model_id=None,
                                             v_min=-3)
         model_plot.source_plot(ax=axes[1, 0], deltaPix_source=0.02, numPix=100,
                                band_index=band_index, v_max=v_max, v_min=v_min)
-        model_plot.convergence_plot(ax=axes[1, 1], band_index=band_index,
-                                    cmap=convergence_cmap)
+        if not show_source_light:
+            model_plot.convergence_plot(ax=axes[1, 1], band_index=band_index,
+                                        cmap=convergence_cmap)
+        else:
+            model_plot.decomposition_plot(ax=axes[1, 1],
+                                          text='Source light convolved',
+                                          source_add=True,
+                                          band_index=band_index,
+                                          v_max=v_max, v_min=v_min)
+
         model_plot.magnification_plot(ax=axes[1, 2],
                                       band_index=band_index,
                                       cmap=magnification_cmap)
@@ -394,7 +403,8 @@ def get_reshaped_emcee_chain(self, lens_name, model_id, walker_ratio,
         return chain
 
     def plot_mcmc_trace(self, lens_name, model_id, walker_ratio,
-                        burn_in=-100, verbose=True, fig_width=16):
+                        burn_in=-100, verbose=True, fig_width=16,
+                        show_variables=[]):
         """
         Plot the trace of MCMC walkers.
 
@@ -421,44 +431,60 @@ def plot_mcmc_trace(self, lens_name, model_id, walker_ratio,
         num_params = self.num_params_mcmc
         num_step = chain.shape[1]
 
+        if len(show_variables) == 0:
+            variable_list = np.arange(num_params)
+        else:
+            variable_list = []
+            for i in show_variables:
+                variable_list.append(self.params_mcmc.index(i))
+
         mean_pos = np.zeros((num_params, num_step))
         median_pos = np.zeros((num_params, num_step))
         std_pos = np.zeros((num_params, num_step))
         q16_pos = np.zeros((num_params, num_step))
         q84_pos = np.zeros((num_params, num_step))
 
         # chain = np.empty((nwalker, nstep, ndim), dtype = np.double)
-        for i in np.arange(num_params):
+        for i in variable_list:
             for j in np.arange(num_step):
                 mean_pos[i][j] = np.mean(chain[:, j, i])
                 median_pos[i][j] = np.median(chain[:, j, i])
                 std_pos[i][j] = np.std(chain[:, j, i])
                 q16_pos[i][j] = np.percentile(chain[:, j, i], 16.)
                 q84_pos[i][j] = np.percentile(chain[:, j, i], 84.)
 
-        fig, ax = plt.subplots(num_params, sharex='all',
+        fig, ax = plt.subplots(len(variable_list), sharex='all',
                                figsize=(fig_width, int(fig_width/8) *
-                                        num_params))
+                                        len(variable_list)))
         last = num_step
         medians = []
 
-        for i in range(num_params):
+        for n, i in enumerate(variable_list):
             if verbose:
                 print(self.params_mcmc[i],
                       '{:.4f} ± {:.4f}'.format(median_pos[i][last - 1],
                                                (q84_pos[i][last - 1] -
                                                 q16_pos[i][last - 1]) / 2))
-            # ax[i].plot(mean_pos[i][:3000], c='b')
-            ax[i].plot(median_pos[i][:last], c='g')
-            # ax[i].axhline(np.mean(mean_pos[i][burnin:2900]), c='b')
-            ax[i].axhline(np.median(median_pos[i][burn_in:last]), c='r', lw=1)
-            ax[i].fill_between(np.arange(last), q84_pos[i][:last],
-                               q16_pos[i][:last], alpha=0.4)
-            # ax[i].fill_between(np.arange(last), mean_pos[i][:last] \
-            # +std_pos[i][:last], mean_pos[i][:last]-std_pos[i][:last],
-            # alpha=0.4)
-            ax[i].set_ylabel(self.params_mcmc[i], fontsize=10)
-            ax[i].set_xlim(0, last)
+            if len(show_variables) != 1:
+                # ax[i].plot(mean_pos[i][:3000], c='b')
+                ax[n].plot(median_pos[i][:last], c='g')
+                # ax[i].axhline(np.mean(mean_pos[i][burnin:2900]), c='b')
+                ax[n].axhline(np.median(median_pos[i][burn_in:last]),
+                              c='r', lw=1)
+                ax[n].fill_between(np.arange(last), q84_pos[i][:last],
+                                   q16_pos[i][:last], alpha=0.4)
+                # ax[i].fill_between(np.arange(last), mean_pos[i][:last] \
+                # +std_pos[i][:last], mean_pos[i][:last]-std_pos[i][:last],
+                # alpha=0.4)
+                ax[n].set_ylabel(self.params_mcmc[i], fontsize=8)
+                ax[n].set_xlim(0, last)
+            else:
+                ax.plot(median_pos[i][:last], c='g')
+                ax.axhline(np.median(median_pos[i][burn_in:last]), c='r', lw=1)
+                ax.fill_between(np.arange(last), q84_pos[i][:last],
+                                q16_pos[i][:last], alpha=0.4)
+                ax.set_ylabel(self.params_mcmc[i], fontsize=8)
+                ax.set_xlim(0, last)
 
             medians.append(np.median(median_pos[i][burn_in:last]))
 

dolphin/processor/config.py

@@ -356,10 +356,8 @@ def custom_logL_addition(self, kwargs_lens=None, kwargs_source=None,
 
                 diff = min(abs(pa_light - pa_mass),
                            180 - abs(pa_light - pa_mass))
-                if diff < np.abs(max_delta):
-                    prior += 0.0
-                else:
-                    prior += -np.inf
+                if diff > np.abs(max_delta):
+                    return -10**-15
 
         # Ensure q_mass is smaller than q_light for the lensing galaxy, where
         # q is the ratio between the minor axis to the major axis of a profile
@@ -388,10 +386,8 @@ def custom_logL_addition(self, kwargs_lens=None, kwargs_source=None,
                                            kwargs_lens[0]['e2'])[1]
                 q_light = ellipticity2phi_q(kwargs_lens_light[0]['e1'],
                                             kwargs_lens_light[0]['e2'])[1]
-                if q_mass / q_light >= max_ratio:
-                    prior += 0.0
-                else:
-                    prior += -np.inf
+                if q_mass/q_light < max_ratio:
+                    return -10**-15
 
         # Provide logarithmic_prior on the source light profile
         if 'source_light_option' in self.settings and \
@@ -712,7 +708,7 @@ def get_lens_model_params(self):
                 })
 
                 sigma.append({
-                    'theta_E': .1, 'e1': 0.1, 'e2': 0.1,
+                    'theta_E': .1, 'e1': 0.01, 'e2': 0.01,
                     'gamma': 0.02, 'center_x': 0.1,
                     'center_y': 0.1
                 })
@@ -776,7 +772,7 @@ def get_lens_light_model_params(self):
                     'center_x': np.max(self.pixel_size) / 10.,
                     'center_y': np.max(self.pixel_size) / 10.,
                     'R_sersic': 0.05, 'n_sersic': 0.5,
-                    'e1': 0.1, 'e2': 0.1
+                    'e1': 0.01, 'e2': 0.01
                 })
 
                 lower.append({

dolphin/processor/files.py

@@ -281,7 +281,7 @@ def save_output_h5(self, lens_name, model_id, output):
                                             )
                     subgroup.create_dataset('param_list',
                                             data=np.array(single_output[2],
-                                                          dtype='S10')
+                                                          dtype='S25')
                                             )
                 elif single_output[0] == 'EMCEE':
                     subgroup.create_dataset('samples',
@@ -290,7 +290,7 @@ def save_output_h5(self, lens_name, model_id, output):
                                             )
                     subgroup.create_dataset('param_list',
                                             data=np.array(single_output[2],
-                                                          dtype='S10')
+                                                          dtype='S25')
                                             )
                     subgroup.create_dataset('log_likelihood',
                                             data=np.array(single_output[3],

io_directory_example/settings/lens_system3_config.yml

@@ -35,6 +35,8 @@ lens_option:
     - - theta_E
       - 1.11
       - 0.13
+  constrain_position_angle_from_lens_light: 15
+  limit_mass_eccentricity_from_light: true
 mask:
   centroid_offset:
   - - 0.0

test/test_analysis/test_output.py

@@ -111,7 +111,8 @@ def test_plot_model_overview(self):
 
         fig2 = self.output.plot_model_overview('lens_system2', 'example',
                                                v_min=-3.0, v_max=1.0,
-                                               print_results=True)
+                                               print_results=True,
+                                               show_source_light=True)
 
         plt.close(fig2)
 
@@ -146,6 +147,13 @@ def test_plot_mcmc_trace(self):
 
         plt.close(fig)
 
+        fig2 = self.output.plot_mcmc_trace('lens_system2', 'example', 2,
+                                          verbose=True, burn_in=0,
+                                          show_variables=["gamma_lens"])
+
+        plt.close(fig2)
+
+
     def test_get_reshaped_emcee_chain(self):
         """
         Test `get_reshaped_emcee_chain` method.

test/test_processor/test_config.py

@@ -234,24 +234,24 @@ def test_custom_logL_addition(self):
         :return:
         :rtype:
         """
-        # phi_m = 0 deg, q_m = 0.8
+        # Mass paramters : (phi_m = 0 deg, q_m = 0.8)
         # Satisfy both priors (phi_L = 10 deg, q_L = 0.8)
-        prior1 = self.config.custom_logL_addition(
+        prior = self.config.custom_logL_addition(
             kwargs_lens=[{'e1': 0.111, 'e2': 0.0}],
             kwargs_lens_light=[{'e1': 0.166, 'e2': 0.060}])
-        assert prior1 == 0
+        assert prior == 0
 
         # qm < qL (phi_L = 0 deg, q_L = 0.9)
-        prior2 = self.config.custom_logL_addition(
+        prior = self.config.custom_logL_addition(
             kwargs_lens=[{'e1': 0.111, 'e2': 0.0}],
             kwargs_lens_light=[{'e1': 0.0526, 'e2': 0.0}])
-        assert prior2 == -np.inf
+        assert prior == -1e-15
 
         # Angle out of sync (phi_L = 20 deg, q_L = 0.8)
-        prior3 = self.config.custom_logL_addition(
+        prior = self.config.custom_logL_addition(
             kwargs_lens=[{'e1': 0.111, 'e2': 0.0}],
             kwargs_lens_light=[{'e1': 0.0851, 'e2': 0.0714}])
-        assert prior3 == -np.inf
+        assert prior == -1e-15
 
         # Settings set to False  (phi_L = 20 deg, q_L = 0.9)
         config2 = deepcopy(self.config)
@@ -261,43 +261,51 @@ def test_custom_logL_addition(self):
             'limit_mass_eccentricity_from_light'] = False
         config2.settings['source_light_option'][
             'shapelet_scale_logarithmic_prior'] = False
-        prior4 = config2.custom_logL_addition(
+        prior = config2.custom_logL_addition(
             kwargs_lens=[{'e1': 0.111, 'e2': 0.0}],
             kwargs_lens_light=[{'e1': 0.0403, 'e2': 0.0338}])
-        assert prior4 == 0
+        assert prior == 0
 
         # Change setting data type
-        config3 = deepcopy(self.config)
-        config3.settings['lens_option'][
+        config3a = deepcopy(self.config)
+        config3a.settings['lens_option'][
             'constrain_position_angle_from_lens_light'] = True
-        config3.settings['lens_option'][
+        prior = config3a.custom_logL_addition(
+            kwargs_lens=[{'e1': 0.111, 'e2': 0.0}],
+            kwargs_lens_light=[{'e1': 0.0851, 'e2': 0.0714}])
+        assert prior == -1e-15
+
+        config3b = deepcopy(self.config)
+        config3b.settings['lens_option'][
             'limit_mass_eccentricity_from_light'] = 1.0
-        prior5 = config3.custom_logL_addition(
+        prior = config3b.custom_logL_addition(
             kwargs_lens=[{'e1': 0.111, 'e2': 0.0}],
-            kwargs_lens_light=[{'e1': 0.0403, 'e2': 0.0338}])
-        assert prior5 == -np.inf
+            kwargs_lens_light=[{'e1': 0.0526, 'e2': 0.0}])
+        assert prior == -1e-15
 
         # Raise error when settings are not bool, int or float
-        config4 = deepcopy(self.config)
-        config4.settings['lens_option'][
+        config4a = deepcopy(self.config)
+        config4a.settings['lens_option'][
             'constrain_position_angle_from_lens_light'] = "Test"
         with pytest.raises(TypeError):
-            config4.custom_logL_addition(
+            config4a.custom_logL_addition(
                 kwargs_lens=[{'e1': 0.111, 'e2': 0.0}],
-                kwargs_lens_light=[{'e1': 0.0403, 'e2': 0.0338}])
+                kwargs_lens_light=[{'e1': 0.166, 'e2': 0.060}])
 
-        config5 = deepcopy(self.config)
-        config5.settings['lens_option'][
+        config4b = deepcopy(self.config)
+        config4b.settings['lens_option'][
             'limit_mass_eccentricity_from_light'] = "Test"
         with pytest.raises(TypeError):
-            config5.custom_logL_addition(
+            config4b.custom_logL_addition(
                 kwargs_lens=[{'e1': 0.111, 'e2': 0.0}],
-                kwargs_lens_light=[{'e1': 0.0403, 'e2': 0.0338}])
+                kwargs_lens_light=[{'e1': 0.166, 'e2': 0.060}])
 
-        # Test Jeffrey's Prior
-        prior6 = self.config3.custom_logL_addition(
+        # Test logarithmic shapelets prior
+        prior = self.config3.custom_logL_addition(
+            kwargs_lens=[{'e1': 0.111, 'e2': 0.0}],
+            kwargs_lens_light=[{'e1': 0.166, 'e2': 0.060}],
             kwargs_source=[{'beta': 0.1}, {'beta': 0.1}])
-        assert round(prior6, 2) == 4.61
+        assert round(prior, 2) == 4.61
 
     def test_get_masks(self):
         """