Coverage for frank/make_figs.py: 87%
Shortcuts on this page
r m x toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1# Frankenstein: 1D disc brightness profile reconstruction from Fourier data
2# using non-parametric Gaussian Processes
3#
4# Copyright (C) 2019-2020 R. Booth, J. Jennings, M. Tazzari
5#
6# This program is free software: you can redistribute it and/or modify
7# it under the terms of the GNU General Public License as published by
8# the Free Software Foundation, either version 3 of the License, or
9# (at your option) any later version.
10#
11# This program is distributed in the hope that it will be useful,
12# but WITHOUT ANY WARRANTY; without even the implied warranty of
13# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14# GNU General Public License for more details.
15#
16# You should have received a copy of the GNU General Public License
17# along with this program. If not, see <https://www.gnu.org/licenses/>
18#
19"""This module generates figures for a Frankenstein fit and its diagnostics.
20"""
21import os
22import numpy as np
23import matplotlib.pyplot as plt
24from matplotlib.gridspec import GridSpec
25from matplotlib.colors import PowerNorm
26import logging
28from frank.utilities import UVDataBinner
30from frank.plot import (
31 plot_deprojection_effect,
32 plot_brightness_profile,
33 plot_vis_quantity,
34 plot_vis_hist,
35 plot_iterations,
36 plot_2dsweep,
37 plot_convergence_criterion
38)
40# Suppress some benign warnings
41import warnings
42warnings.filterwarnings('ignore', '.*handles with labels found.*')
45# Global settings for plots
46cs = ['#a4a4a4', 'k', '#f781bf', '#dede00']
47cs2 = ['#3498DB', '#984ea3', '#4daf4a', '#ff7f00']
48hist_cs = ['#377eb8', '#ff7f00', '#e41a1c', '#999999', '#4daf4a', '#f781bf',
49 '#984ea3', '#dede00']
50multifit_cs = ['#e41a1c', '#999999', '#377eb8', '#ff7f00', '#4daf4a', '#f781bf',
51 '#984ea3', '#dede00']
52ms = ['x', '+', '.', '1']
55class _DoNothingContextManager(object):
56 def __init__(self, dummy_resource=None):
57 self.dummy_resource = dummy_resource
58 def __enter__(self):
59 return self.dummy_resource
60 def __exit__(self, *args):
61 pass
63def frank_plotting_style_context_manager(use_frank_style=True):
64 """Get a context manager for temporary use of frank's own plotting style"""
65 if use_frank_style:
66 frank_path = os.path.dirname(__file__)
67 style_path = os.path.join(frank_path, 'frank.mplstyle')
68 return plt.style.context(style_path)
69 else:
70 return _DoNothingContextManager()
72def use_frank_plotting_style():
73 """Set matplotlib to use frank's own plotting style"""
74 frank_path = os.path.dirname(__file__)
75 style_path = os.path.join(frank_path, 'frank.mplstyle')
76 plt.style.use(style_path)
79def make_deprojection_fig(u, v, vis, weights, geom, bin_widths, logx=False,
80 logy=False, force_style=True, save_prefix=None,
81 figsize=(8, 6)):
82 r"""
83 Produce a simple figure showing the effect of deprojection on the (u, v)
84 coordinates and visibilities
86 Parameters
87 ----------
88 u, v : array, unit = :math:`\lambda`
89 Projected (u, v) coordinates
90 vis : array, unit = Jy
91 Projected visibilities (complex: real + imag * 1j)
92 weights : array, unit = Jy^-2
93 Weights assigned to observed visibilities, of the form
94 :math:`1 / \sigma^2`
95 geom : SourceGeometry object
96 Fitted geometry (see frank.geometry.SourceGeometry)
97 bin_widths : list, unit = \lambda
98 Bin widths in which to bin the observed visibilities
99 logy : bool, default = False
100 Whether to plot the visibility distributions in log(flux)
101 force_style: bool, default = True
102 Whether to use preconfigured matplotlib rcParams in generated figure
103 save_prefix : string, default = None
104 Prefix for saved figure name. If None, the figure won't be saved
105 figsize : tuple = (width, height) of figure, unit = inch
107 Returns
108 -------
109 fig : Matplotlib `.Figure` instance
110 The produced figure, including the GridSpec
111 axes : Matplotlib `~.axes.Axes` class
112 The axes of the produced figure
113 """
115 logging.info(' Making deprojection figure')
117 # Apply the deprojection to the provided (u, v) coordinates
118 # and visibility amplitudes
119 up, vp, visp = geom.apply_correction(u, v, vis)
121 re_vis = np.real(vis)
122 re_visp = np.real(visp)
124 bsp = np.hypot(up, vp)
126 with frank_plotting_style_context_manager(force_style):
127 gs = GridSpec(2, 1, hspace=0.2, top=0.97, left=0.07, right=0.97)
128 fig = plt.figure(figsize=figsize)
130 ax0 = fig.add_subplot(gs[0])
131 ax1 = fig.add_subplot(gs[1])
132 axes = [ax0, ax1]
134 plot_deprojection_effect(u / 1e6, v / 1e6, up / 1e6, vp / 1e6,
135 re_vis * 1e3, re_visp * 1e3, axes)
137 for i in range(len(bin_widths)):
138 binned_vis = UVDataBinner(bsp, visp, weights, bin_widths[i])
139 vis_re = binned_vis.V.real
140 vis_err_re = binned_vis.error.real
142 if logy:
143 lab=r'Deprojected, >0, {:.0f} k$\lambda$ bins'.format(bin_widths[i]/1e3)
144 else:
145 lab=r'Deprojected, {:.0f} k$\lambda$ bins'.format(bin_widths[i]/1e3)
146 plot_vis_quantity(binned_vis.uv / 1e6, vis_re * 1e3, ax1, c=cs[i],
147 marker=ms[i], ls='None',
148 label=lab)
150 if logy:
151 plot_vis_quantity(binned_vis.uv / 1e6, -vis_re * 1e3, ax1, c=cs2[i],
152 marker=ms[i], ls='None',
153 label=r'Deprojected, <0, {:.0f} k$\lambda$ bins'.format(bin_widths[i]/1e3))
155 ax0.set_xlabel(r'u [M$\lambda$]')
156 ax0.set_ylabel(r'v [M$\lambda$]')
158 ax1.set_xlabel(r'Baseline [M$\lambda$]')
159 ax1.set_ylabel('Re(V) [mJy]')
160 if logx:
161 ax1.set_xscale('log')
162 if logy:
163 ax1.set_yscale('log')
164 ax1.set_ylim(bottom=1e-3)
166 ax0.legend(loc=0)
167 ax1.legend(loc=0)
169 if save_prefix:
170 plt.savefig(save_prefix + '_frank_deprojection.png', dpi=300,
171 bbox_inches='tight')
172 plt.close()
174 return fig, axes
177def make_quick_fig(u, v, vis, weights, sol, bin_widths, dist=None,
178 logx=False, force_style=True, save_prefix=None,
179 stretch='power', gamma=1.0, asinh_a=0.02, figsize=(8,6)):
180 r"""
181 Produce a simple figure showing just a Frankenstein fit, not any diagnostics
183 Parameters
184 ----------
185 u, v : array, unit = :math:`\lambda`
186 u and v coordinates of observations
187 vis : array, unit = Jy
188 Observed visibilities (complex: real + imag * 1j)
189 weights : array, unit = Jy^-2
190 Weights assigned to observed visibilities, of the form
191 :math:`1 / \sigma^2`
192 sol : _HankelRegressor object
193 Reconstructed profile using Maximum a posteriori power spectrum
194 (see frank.radial_fitters.FrankFitter)
195 bin_widths : list, unit = \lambda
196 Bin widths in which to bin the observed visibilities
197 dist : float, optional, unit = AU, default = None
198 Distance to source, used to show second x-axis in [AU]
199 logx : bool, default = False
200 Whether to plot the visibility distributions in log(baseline)
201 gamma : float, default = 1.0
202 Index of power law normalization to apply to swept profile image's
203 colormap (see matplotlib.colors.PowerNorm).
204 gamma=1.0 yields a linear colormap
205 force_style: bool, default = True
206 Whether to use preconfigured matplotlib rcParams in generated figure
207 save_prefix : string, default = None
208 Prefix for saved figure name. If None, the figure won't be saved
209 stretch : string, default = 'power'
210 Transformation to apply to the colorscale. The default 'power' is a
211 power law stretch. The other option is 'asinh', an arcsinh stretch,
212 which requires astropy.visualization.mpl_normalize.simple_norm
213 asinh_a : float, default = 0.02
214 Scale parameter for an asinh stretch
215 figsize : tuple = (width, height) of figure, unit = inch
217 Returns
218 -------
219 fig : Matplotlib `.Figure` instance
220 The produced figure, including the GridSpec
221 axes : Matplotlib `~.axes.Axes` class
222 The axes of the produced figure
223 """
225 logging.info(' Making quick figure')
227 Rmax, N, alpha, wsmooth, p0 = sol._info["Rmax"], sol._info["N"], \
228 sol._info["alpha"], sol._info["wsmooth"], sol._info["p0"]
230 with frank_plotting_style_context_manager(force_style):
231 gs = GridSpec(3, 2, hspace=0, bottom=.12)
232 gs2 = GridSpec(3, 2, hspace=.2)
233 fig = plt.figure(figsize=figsize)
234 fig.suptitle(r'Fit hyperparameters: $\alpha$={:.2f}, $w_{{smooth}}$={:.1e}, R$_{{max}}$={:.2f}, '
235 'N={:0d}, p$_0$={:.0e}'.format(alpha, wsmooth, Rmax, N, p0))
237 ax0 = fig.add_subplot(gs[0])
238 ax1 = fig.add_subplot(gs[2])
240 ax2 = fig.add_subplot(gs[1])
241 ax3 = fig.add_subplot(gs[3])
243 ax4 = fig.add_subplot(gs2[4])
244 ax5 = fig.add_subplot(gs2[5])
246 ax0.text(.5, .9, 'a)', transform=ax0.transAxes)
247 ax1.text(.5, .9, 'b)', transform=ax1.transAxes)
249 ax2.text(.5, .9, 'c)', transform=ax2.transAxes)
250 ax3.text(.5, .9, 'd)', transform=ax3.transAxes)
252 ax4.text(.5, .9, 'e)', c='w', transform=ax4.transAxes)
253 ax5.text(.5, .9, 'f)', c='w', transform=ax5.transAxes)
255 axes = [ax0, ax1, ax2, ax3, ax4, ax5]
257 plot_brightness_profile(sol.r, sol.I / 1e10, ax0, c='r', label='frank')
259 if dist:
260 ax0_5 = plot_brightness_profile(sol.r, sol.I / 1e10, ax0, dist=dist, c='r')
261 xlims = ax0.get_xlim()
262 ax0_5.set_xlim(np.multiply(xlims, dist))
264 plot_brightness_profile(sol.r, sol.I / 1e10, ax1, c='r', label='frank')
266 if hasattr(sol, 'nonneg'):
267 plot_brightness_profile(sol.r, sol.nonneg / 1e10, ax0, ls='--', c='#009933', label='non-neg.')
268 plot_brightness_profile(sol.r, sol.nonneg / 1e10, ax1, ls='--', c='#009933', label='non-neg.')
271 u_deproj, v_deproj, vis_deproj = sol.geometry.apply_correction(u, v, vis)
272 baselines = (u_deproj**2 + v_deproj**2)**.5
273 grid = np.logspace(np.log10(min(baselines.min(), sol.q[0])),
274 np.log10(max(baselines.max(), sol.q[-1])),
275 10**4)
277 for i in range(len(bin_widths)):
278 binned_vis = UVDataBinner(
279 baselines, vis_deproj, weights, bin_widths[i])
280 vis_re = binned_vis.V.real * 1e3
281 vis_err_re = binned_vis.error.real * 1e3
282 vis_fit = sol.predict_deprojected(binned_vis.uv).real * 1e3
284 for ax in [ax2, ax3]:
285 plot_vis_quantity(binned_vis.uv / 1e6, vis_re, ax,
286 vis_err_re, c=cs[i],
287 marker=ms[i], ls='None',
288 label=r'Obs., {:.0f} k$\lambda$ bins'.format(bin_widths[i]/1e3))
290 vis_fit = sol.predict_deprojected(grid).real * 1e3
291 plot_vis_quantity(grid / 1e6, vis_fit, ax2, c='r', label='frank', zorder=10)
293 plot_vis_quantity(grid / 1e6, vis_fit, ax3, c='r', label='frank', zorder=10)
295 if hasattr(sol, 'nonneg'):
296 vis_fit_nonneg = sol.predict_deprojected(grid, I=sol.nonneg).real * 1e3
297 plot_vis_quantity(grid / 1e6, vis_fit_nonneg, ax2, ls='--', c='#009933', label='non-neg.', zorder=10)
298 plot_vis_quantity(grid / 1e6, vis_fit_nonneg, ax3, ls='--', c='#009933', label='non-neg.', zorder=10)
300 # Make a guess of good y-bounds for zooming in on the visibility fit
301 zoom_ylim_guess = vis_fit.mean()
302 zoom_bounds = [-1.5 * zoom_ylim_guess, 1.5 * zoom_ylim_guess]
303 ax3.set_ylim(zoom_bounds)
305 vmax = sol.I.max()
306 if stretch == 'asinh':
307 vmin = max(0, min(sol.I))
308 from astropy.visualization.mpl_normalize import simple_norm
309 norm = simple_norm(sol.I, stretch='asinh', asinh_a=asinh_a,
310 min_cut=vmin, max_cut=vmax)
311 elif stretch == 'power':
312 vmin = 0
313 norm = PowerNorm(gamma, vmin, vmax)
314 else:
315 err = ValueError("Unknown 'stretch'. Should be one of 'power' or 'asinh'")
316 raise err
318 plot_2dsweep(sol.r, sol.I, ax=ax4, cmap='inferno', norm=norm,
319 project=False)
320 plot_2dsweep(sol.r, sol.I, ax=ax5, cmap='inferno', norm=norm,
321 project=True, geom=sol.geometry)
323 ax1.set_xlabel('r ["]')
324 ax0.set_ylabel(r'Brightness [$10^{10}$ Jy sr$^{-1}$]')
325 ax1.set_ylabel(r'Brightness [$10^{10}$ Jy sr$^{-1}$]')
326 ax1.set_yscale('log')
327 ax1.set_ylim(sol.I[sol.I > 0].min() * 0.9 / 1e10, sol.I.max() * 1.1 / 1e10)
329 ax3.set_xlabel(r'Baseline [M$\lambda$]')
330 ax2.set_ylabel('Re(V) [mJy]')
331 ax3.set_ylabel('Re(V) [mJy]')
333 if logx:
334 ax2.set_xscale('log')
335 ax3.set_xscale('log')
337 ax2.set_xlim(right=max(baselines) / 1e6 * 1.2)
338 else:
339 ax2.set_xlim(0, max(baselines) / 1e6 * 1.2)
340 xlims = ax2.get_xlim()
341 ax3.set_xlim(xlims)
343 ax4.set_xlabel('RA offset ["]')
344 ax4.set_ylabel('Dec offset ["]')
345 ax5.set_xlabel('RA offset ["]')
346 ax5.set_ylabel('Dec offset ["]')
348 plt.setp(ax0.get_xticklabels(), visible=False)
349 plt.setp(ax2.get_xticklabels(), visible=False)
351 if save_prefix:
352 plt.savefig(save_prefix + '_frank_fit_quick.png', dpi=300,
353 bbox_inches='tight')
354 plt.close()
356 return fig, axes
359def make_full_fig(u, v, vis, weights, sol, bin_widths,
360 dist=None, logx=False, force_style=True,
361 save_prefix=None, norm_residuals=False, stretch='power',
362 gamma=1.0, asinh_a=0.02, figsize=(8, 6)):
363 r"""
364 Produce a figure showing a Frankenstein fit and some useful diagnostics
366 Parameters
367 ----------
368 u, v : array, unit = :math:`\lambda`
369 u and v coordinates of observations
370 vis : array, unit = Jy
371 Observed visibilities (complex: real + imag * 1j)
372 weights : array, unit = Jy^-2
373 Weights assigned to observed visibilities, of the form
374 :math:`1 / \sigma^2`
375 sol : _HankelRegressor object
376 Reconstructed profile using maximum a posteriori power spectrum
377 (see frank.radial_fitters.FrankFitter)
378 bin_widths : list, unit = \lambda
379 Bin widths in which to bin the observed visibilities
380 dist : float, optional, unit = AU, default = None
381 Distance to source, used to show second x-axis in [AU]
382 logx : bool, default = False
383 Whether to plot the visibility distributions in log(baseline)
384 force_style: bool, default = True
385 Whether to use preconfigured matplotlib rcParams in generated figure
386 save_prefix : string, default = None
387 Prefix for saved figure name. If None, the figure won't be saved
388 norm_residuals : bool, default = False
389 Whether to normalize the residual visibilities by the data's
390 visibility amplitudes
391 stretch : string, default = 'power'
392 Transformation to apply to the colorscale. The default 'power' is a
393 power law stretch. The other option is 'asinh', an arcsinh stretch,
394 which requires astropy.visualization.mpl_normalize.simple_norm
395 gamma : float, default = 1.0
396 Index of power law normalization to apply to swept profile image's
397 colormap (see matplotlib.colors.PowerNorm).
398 gamma=1.0 yields a linear colormap
399 asinh_a : float, default = 0.02
400 Scale parameter for an asinh stretch
401 figsize : tuple = (width, height) of figure, unit = inch
403 Returns
404 -------
405 fig : Matplotlib `.Figure` instance
406 The produced figure, including the GridSpec
407 axes : Matplotlib `~.axes.Axes` class
408 The axes of the produced figure
409 """
411 logging.info(' Making full figure')
413 Rmax, N, alpha, wsmooth, p0 = sol._info["Rmax"], sol._info["N"], \
414 sol._info["alpha"], sol._info["wsmooth"], sol._info["p0"]
416 with frank_plotting_style_context_manager(force_style):
417 gs = GridSpec(3, 3, hspace=0, wspace=0.25, left=0.06, right=0.98)
418 gs1 = GridSpec(4, 3, hspace=0, wspace=0.25, left=0.06, right=0.98)
419 gs2 = GridSpec(3, 3, hspace=.35, wspace=0.25, left=0.06, right=0.98)
420 fig = plt.figure(figsize=figsize)
421 fig.suptitle(r'Fit hyperparameters: $\alpha$={:.2f}, $w_{{smooth}}$={:.1e}, R$_{{max}}$={:.2f}, '
422 'N={:0d}, p$_0$={:.0e}'.format(alpha, wsmooth, Rmax, N, p0))
424 ax0 = fig.add_subplot(gs[0])
425 ax1 = fig.add_subplot(gs[3])
426 ax2 = fig.add_subplot(gs2[6])
428 ax3 = fig.add_subplot(gs[1])
429 ax4 = fig.add_subplot(gs[4])
430 ax5 = fig.add_subplot(gs[7])
432 ax6 = fig.add_subplot(gs1[2])
433 ax7 = fig.add_subplot(gs1[5])
434 ax8 = fig.add_subplot(gs1[8])
435 ax9 = fig.add_subplot(gs1[11])
437 ax0.text(.9, .6, 'a)', transform=ax0.transAxes)
438 ax1.text(.9, .6, 'b)', transform=ax1.transAxes)
439 ax2.text(.1, .9, 'c)', c='w', transform=ax2.transAxes)
441 ax3.text(.9, .5, 'd)', transform=ax3.transAxes)
442 ax4.text(.9, .1, 'e)', transform=ax4.transAxes)
443 ax5.text(.9, .1, 'f)', transform=ax5.transAxes)
444 ax6.text(.9, .1, 'g)', transform=ax6.transAxes)
445 ax7.text(.9, .9, 'h)', transform=ax7.transAxes)
446 ax8.text(.9, .1, 'i)', transform=ax8.transAxes)
447 ax9.text(.9, .1, 'j)', transform=ax9.transAxes)
449 axes = [ax0, ax1, ax2, ax3, ax4, ax5, ax6, ax7, ax8, ax9]
451 # Plot the fitted brightness profile in linear- and log-y
452 plot_brightness_profile(sol.r, sol.I / 1e10, ax0, c='r', label='frank')
454 if dist:
455 ax0_5 = plot_brightness_profile(sol.r, sol.I / 1e10, ax0, dist=dist, c='r')
456 xlims = ax0.get_xlim()
457 ax0_5.set_xlim(np.multiply(xlims, dist))
459 plot_brightness_profile(sol.r, sol.I / 1e10, ax1, c='r', label='frank')
461 if hasattr(sol, 'nonneg'):
462 plot_brightness_profile(sol.r, sol.nonneg / 1e10, ax0, ls='--', c='#009933', label='non-neg.')
463 plot_brightness_profile(sol.r, sol.nonneg / 1e10, ax1, ls='--', c='#009933', label='non-neg.')
465 # Apply deprojection to the provided (u, v) coordinates
466 # and visibility amplitudes
467 u_deproj, v_deproj, vis_deproj = sol.geometry.apply_correction(u, v, vis)
468 baselines = (u_deproj**2 + v_deproj**2)**.5
469 # Set a grid of baselines on which to plot the visibility domain frank fit
470 grid = np.logspace(np.log10(min(baselines.min(), sol.q[0])),
471 np.log10(max(baselines.max(), sol.q[-1])), 10**4
472 )
474 # Bin the observed (real and imaginary components of the) visibilities
475 # for plotting
476 for i in range(len(bin_widths)):
477 binned_vis = UVDataBinner(baselines, vis_deproj, weights, bin_widths[i])
478 vis_re = binned_vis.V.real * 1e3
479 vis_im = binned_vis.V.imag * 1e3
480 vis_err_re = binned_vis.error.real * 1e3
481 vis_err_im = binned_vis.error.imag * 1e3
482 vis_fit = sol.predict_deprojected(binned_vis.uv).real * 1e3
483 if hasattr(sol, 'nonneg'):
484 vis_fit_nn = sol.predict_deprojected(binned_vis.uv, I=sol.nonneg).real * 1e3
486 # Determine the visiblity domain frank fit residuals (and RMS error)
487 # for Real(V)
488 resid = vis_re - vis_fit
489 if norm_residuals:
490 resid /= vis_re
491 rmse = (np.mean(resid**2))**.5
492 if hasattr(sol, 'nonneg'):
493 resid_nn = vis_re - vis_fit_nn
494 if norm_residuals:
495 resid_nn /= vis_re
496 rmse_nn = (np.mean(resid_nn**2))**.5
498 # Plot the observed, binned visibilities (with errorbars) and the residuals
499 plot_vis_quantity(binned_vis.uv / 1e6, vis_re, ax3, c=cs[i],
500 marker=ms[i], ls='None',
501 label=r'Obs., {:.0f} k$\lambda$ bins'.format(bin_widths[i]/1e3))
503 plot_vis_quantity(binned_vis.uv / 1e6, vis_re, ax4, c=cs[i],
504 marker=ms[i], ls='None',
505 label=r'Obs., {:.0f} k$\lambda$ bins'.format(bin_widths[i]/1e3))
507 plot_vis_quantity(binned_vis.uv / 1e6, vis_im, ax6, c=cs[i],
508 marker=ms[i], ls='None',
509 label=r'Obs., {:.0f} k$\lambda$ bins'.format(bin_widths[i]/1e3))
511 plot_vis_quantity(binned_vis.uv / 1e6, resid, ax5, c='r', marker=ms[i], ls='None',
512 label=r'{:.0f} k$\lambda$ bins, RMSE {:.3f} mJy'.format(bin_widths[i]/1e3, rmse))
514 if hasattr(sol, 'nonneg'):
515 plot_vis_quantity(binned_vis.uv / 1e6, resid_nn, ax5, c='#009933', marker=ms[i], ls='None',
516 label=r'{:.0f} k$\lambda$ bins, RMSE {:.3f} mJy'.format(bin_widths[i]/1e3, rmse_nn))
518 # Plot a histogram of the observed visibilties to examine how the
519 # visibility count varies with baseline
520 plot_vis_hist(binned_vis, ax8, rescale=1e6, color=hist_cs[i],
521 label=r'Obs., {:.0f} k$\lambda$ bins'.format(bin_widths[i]/1e3))
522 # Plot the binned data signal-to-noise as a function of baseline
523 plot_vis_quantity(binned_vis.uv / 1e6,
524 binned_vis.V.real**2 / binned_vis.error.real**2,
525 ax9, color=hist_cs[i], marker=ms[i], ls='None',
526 label=r'{:.0f} k$\lambda$ bins'.format(bin_widths[i]/1e3))
528 # Map the frank visibility fit to `grid`, considering only the real component
529 # that frank fits
530 vis_fit = sol.predict_deprojected(grid).real * 1e3
532 # Plot the visibility domain frank fit
533 plot_vis_quantity(grid / 1e6, vis_fit, ax3, c='r', label='frank')
534 plot_vis_quantity(grid / 1e6, vis_fit, ax4, c='r', label='frank')
536 if hasattr(sol, 'nonneg'):
537 vis_fit_nonneg = sol.predict_deprojected(grid, I=sol.nonneg).real * 1e3
538 plot_vis_quantity(grid / 1e6, vis_fit_nonneg, ax3, ls='--', c='#009933', label='non-neg.')
539 plot_vis_quantity(grid / 1e6, vis_fit_nonneg, ax4, ls='--', c='#009933', label='non-neg.')
541 # Make a guess of good y-bounds for zooming in on the visibility fit
542 zoom_ylim_guess = vis_fit.mean()
543 zoom_bounds = [-1.5 * zoom_ylim_guess, 1.5 * zoom_ylim_guess]
544 ax4.set_ylim(zoom_bounds)
546 # Plot the frank inferred power spectrum
547 plot_vis_quantity(sol.q / 1e6, sol.power_spectrum, ax7, c='k')
549 # Plot a sweep over 2\pi of the frank 1D fit
550 # (analogous to a model image of the source)
551 vmax = sol.I.max()
552 if stretch == 'asinh':
553 vmin = max(0, min(sol.I))
554 from astropy.visualization.mpl_normalize import simple_norm
555 norm = simple_norm(sol.I, stretch='asinh', asinh_a=asinh_a,
556 min_cut=vmin, max_cut=vmax)
557 elif stretch == 'power':
558 vmin = 0
559 norm = PowerNorm(gamma, vmin, vmax)
560 else:
561 err = ValueError("Unknown 'stretch'. Should be one of 'power' or 'asinh'")
562 raise err
564 plot_2dsweep(sol.r, sol.I, ax=ax2, cmap='inferno', norm=norm,
565 project=True, geom=sol.geometry)
567 ax1.set_xlabel('r ["]')
568 ax0.set_ylabel(r'Brightness [$10^{10}$ Jy sr$^{-1}$]')
569 ax1.set_ylabel(r'Brightness [$10^{10}$ Jy sr$^{-1}$]')
570 ax1.set_yscale('log')
571 ax1.set_ylim(sol.I[sol.I > 0].min() * 0.9 / 1e10, sol.I.max() * 1.1 / 1e10)
572 ax2.set_xlabel('RA offset ["]')
573 ax2.set_ylabel('Dec offset ["]')
575 axs = [ax3, ax4, ax5, ax6, ax7, ax8, ax9]
576 if logx:
577 for aa in axs:
578 aa.set_xscale('log')
579 ax3.set_xlim(right=max(baselines) / 1e6 * 1.2)
580 else:
581 ax3.set_xlim(0, max(baselines) / 1e6 * 1.2)
582 xlims = ax3.get_xlim()
583 for aa in axs:
584 aa.set_xlim(xlims)
586 ax7.set_yscale('log')
587 ax8.set_yscale('log')
589 ax8.set_yticks([1, 10, 1e2, 1e3])
591 ax9.set_yscale('log')
592 ax9.set_ylim(1e-2,200)
593 ax9.set_yticks([1e-2, 1e-1, 1, 10, 100])
594 ax9.axhline(y=1, ls='--', c='c')
596 ax3.set_ylabel('Re(V) [mJy]')
597 ax4.set_ylabel('Re(V) [mJy]')
598 if norm_residuals:
599 ax5.set_ylabel('Norm. residual')
600 else:
601 ax5.set_ylabel('Residual [mJy]')
602 ax5.set_xlabel(r'Baseline [M$\lambda$]')
604 ax6.set_ylabel('Im(V) [mJy]')
605 ax7.set_ylabel(r'Power [Jy$^2$]')
606 ax8.set_ylabel('Count')
607 ax9.set_ylabel(r'SNR=$\mu_{\rm bin}^2 / \sigma_{\rm bin}^2$')
608 ax9.set_xlabel(r'Baseline [M$\lambda$]')
610 axs = [ax0, ax3, ax4, ax6, ax7, ax8]
611 for aa in axs:
612 plt.setp(aa.get_xticklabels(), visible=False)
614 if save_prefix:
615 plt.savefig(save_prefix + '_frank_fit_full.png', dpi=300)
616 plt.close()
618 return fig, axes
621def make_diag_fig(r, q, iteration_diagnostics, iter_plot_range=None, logx=False,
622 force_style=True, save_prefix=None, figsize=(8, 6)):
623 r"""
624 Produce a diagnostic figure showing fit convergence metrics
626 Parameters
627 ----------
628 r : array
629 Radial data coordinates at which the brightness profile is defined.
630 The assumed unit (for the x-label) is arcsec
631 q : array
632 Baselines at which the power spectrum is defined.
633 The assumed unit (for the x-label) is :math:`\lambda`
634 iteration_diagnostics : dict
635 The iteration_diagnostics from FrankFitter.
636 iter_plot_range : list
637 Range of iterations in the fit over which to
638 plot brightness profile and power spectrum reconstructions
639 logx : bool, default = False
640 Whether to plot the visibility distributions in log(baseline)
641 force_style: bool, default = True
642 Whether to use preconfigured matplotlib rcParams in generated figure
643 save_prefix : string, default = None
644 Prefix for saved figure name. If None, the figure won't be saved
645 figsize : tuple = (width, height) of figure, unit = inch
647 Returns
648 -------
649 fig : Matplotlib `.Figure` instance
650 The produced figure, including the GridSpec
651 axes : Matplotlib `~.axes.Axes` class
652 The axes of the produced figure
653 """
655 logging.info(' Making diagnostic figure')
657 if iter_plot_range is None:
658 logging.info(" diag_plot is 'true' in your parameter file but"
659 " iter_plot_range is 'null' --> Defaulting to"
660 " plotting all iterations")
662 iter_plot_range = [0, iteration_diagnostics['num_iterations']]
664 else:
665 if iter_plot_range[0] > iteration_diagnostics['num_iterations']:
666 logging.info(' iter_plot_range[0] in your parameter file'
667 ' exceeds the number of fit iterations -->'
668 ' Defaulting to plotting all iterations')
670 iter_plot_range = [0, iteration_diagnostics['num_iterations']]
672 with frank_plotting_style_context_manager(force_style):
673 gs = GridSpec(2, 2, hspace=0, bottom=.35)
674 gs2 = GridSpec(3, 2, hspace=0, top=.7)
675 fig = plt.figure(figsize=figsize)
677 ax0 = fig.add_subplot(gs[0])
678 ax1 = fig.add_subplot(gs[2])
680 ax2 = fig.add_subplot(gs[1])
681 ax3 = fig.add_subplot(gs[3])
683 ax4 = fig.add_subplot(gs2[4])
685 ax0.text(.92, .5, 'a)', transform=ax0.transAxes)
686 ax1.text(.92, .1, 'b)', transform=ax1.transAxes)
688 ax2.text(.05, .5, 'c)', transform=ax2.transAxes)
689 ax3.text(.92, .5, 'd)', transform=ax3.transAxes)
691 ax4.text(.92, .9, 'e)', transform=ax4.transAxes)
693 axes = [ax0, ax1, ax2, ax3, ax4]
695 profile_iter = iteration_diagnostics['MAP']
696 profile_iter_toplot = [x / 1e10 for x in profile_iter]
697 pwr_spec_iter = iteration_diagnostics['power_spectrum']
698 num_iter = iteration_diagnostics['num_iterations']
700 plot_iterations(r, profile_iter_toplot, iter_plot_range, ax0)
702 # Plot the difference in the profile between the last 100 iterations
703 iter_plot_range_end = [max(iter_plot_range[1] - 100, 0),
704 iter_plot_range[1] - 1]
706 plot_iterations(r, np.diff(profile_iter_toplot, axis=0),
707 iter_plot_range_end, ax1,
708 cmap=plt.cm.cividis) # pylint: disable=no-member
710 plot_iterations(q, pwr_spec_iter, iter_plot_range, ax2)
712 # Plot the difference in the power spectrum between the last 100 iterations
713 plot_iterations(q, abs(np.diff(pwr_spec_iter, axis=0)),
714 iter_plot_range_end, ax3,
715 cmap=plt.cm.cividis) # pylint: disable=no-member
717 plot_convergence_criterion(profile_iter_toplot, num_iter, ax4, c='k')
719 ax0.set_ylabel(r'I [$10^{10}$ Jy sr$^{-1}$]')
720 ax1.set_ylabel(r'$I_i - I_{i-1}$ [$10^{10}$ Jy sr$^{-1}$]')
721 ax1.set_xlabel('r ["]')
723 ax2.set_ylabel(r'Power [Jy$^2$]')
724 ax3.set_ylabel(r'|PS$_i$ - PS$_{i-1}$| [Jy$^2$]')
725 ax3.set_xlabel(r'Baseline [$\lambda$]')
726 if logx:
727 ax2.set_xscale('log')
728 ax3.set_xscale('log')
729 ax2.set_yscale('log')
730 ax3.set_yscale('log')
732 ax4.set_xlabel('Fit iteration')
733 ax4.set_ylabel('Convergence criterion,\n' +
734 r'max(|$I_i - I_{i-1}$|) / max($I_i$)')
735 ax4.set_yscale('log')
737 xlims = ax2.get_xlim()
738 ax3.set_xlim(xlims)
740 plt.setp(ax0.get_xticklabels(), visible=False)
741 plt.setp(ax2.get_xticklabels(), visible=False)
743 if save_prefix:
744 plt.savefig(save_prefix + '_frank_fit_diag.png', dpi=300,
745 bbox_inches='tight')
746 plt.close()
748 return fig, axes, iter_plot_range
751def make_clean_comparison_fig(u, v, vis, weights, sol, clean_profile,
752 bin_widths, stretch='power',
753 gamma=1.0, asinh_a=0.02, MAP_convolved=None,
754 dist=None, logx=False, logy=False, ylims=None,
755 force_style=True, save_prefix=None,
756 figsize=(8, 10)):
757 r"""
758 Produce a figure comparing a frank fit to a CLEAN fit, in real space by
759 convolving the frank fit with the CLEAN beam, and in visibility space by
760 taking the discrete Hankel transform of the CLEAN profile
762 Parameters
763 ----------
764 u, v : array, unit = :math:`\lambda`
765 u and v coordinates of observations
766 vis : array, unit = Jy
767 Observed visibilities (complex: real + imag * 1j)
768 weights : array, unit = Jy^-2
769 Weights assigned to observed visibilities, of the form
770 :math:`1 / \sigma^2`
771 sol : _HankelRegressor object
772 Reconstructed profile using Maximum a posteriori power spectrum
773 (see frank.radial_fitters.FrankFitter)
774 clean_profile : dict
775 Dictionary with entries 'r' for the radial points [arcsec],
776 'I' for the brightness [Jy / sr], and optionally the negative and positive
777 brightness uncertainties 'lo_err' and 'hi_err' [Jy / sr]. If only the
778 negative uncertainty is provided, the positive uncertainty is assumed
779 equal to it
780 bin_widths : list, unit = \lambda
781 Bin widths in which to bin the observed visibilities
782 stretch : string, default = 'power'
783 Transformation to apply to the colorscale. The default 'power' is a
784 power law stretch. The other option is 'asinh', an arcsinh stretch,
785 which requires astropy.visualization.mpl_normalize.simple_norm
786 gamma : float, default = 1.0
787 Index of power law normalization to apply to swept profile image's
788 colormap (see matplotlib.colors.PowerNorm).
789 gamma=1.0 yields a linear colormap
790 asinh_a : float, default = 0.02
791 Scale parameter for an asinh stretch
792 MAP_convolved : None (default) or array, unit = Jy / sr
793 frank brightness profile convolved with a CLEAN beam
794 (see utilities.convolve_profile).
795 The assumed unit is for the x-label
796 dist : float, optional, unit = AU, default = None
797 Distance to source, used to show second x-axis in [AU]
798 logx : bool, default = False
799 Whether to plot the visibility distributions in log(baseline)
800 logy : bool, default = False
801 Whether to plot the visibility distributions in log(flux)
802 ylims : list of len(2), default = None
803 Lower and upper y-bounds for the visibility domain plot
804 force_style: bool, default = True
805 Whether to use preconfigured matplotlib rcParams in generated figure
806 save_prefix : string, default = None
807 Prefix for saved figure name. If None, the figure won't be saved
808 figsize : tuple = (width, height) of figure, unit = inch
810 Returns
811 -------
812 fig : Matplotlib `.Figure` instance
813 The produced figure, including the GridSpec
814 axes : Matplotlib `~.axes.Axes` class
815 The axes of the produced figure
816 """
818 logging.info(' Making CLEAN comparison figure')
820 with frank_plotting_style_context_manager(force_style):
821 gs = GridSpec(3, 1)
822 gs2 = GridSpec(3, 3)
824 fig = plt.figure(figsize=figsize)
826 ax0 = fig.add_subplot(gs[0])
827 ax1 = fig.add_subplot(gs[1])
829 ax2 = fig.add_subplot(gs2[6])
830 ax3 = fig.add_subplot(gs2[7])
831 ax4 = fig.add_subplot(gs2[8])
833 axes = [ax0, ax1, ax2, ax3, ax4]
835 if 'lo_err' in clean_profile.keys():
836 low_uncer = clean_profile['lo_err']
837 else:
838 low_uncer = None
839 if 'hi_err' in clean_profile.keys():
840 high_uncer = clean_profile['hi_err']
841 else:
842 high_uncer = None
844 plot_brightness_profile(clean_profile['r'], clean_profile['I'] / 1e10, ax0,
845 low_uncer=low_uncer, high_uncer=high_uncer,
846 c='b', ls='--', label='CLEAN')
848 plot_brightness_profile(sol.r, sol.I / 1e10, ax0, c='r', ls=':', label='frank')
850 if MAP_convolved is not None:
851 plot_brightness_profile(sol.r, MAP_convolved / 1e10, ax0, c='k', ls='-',
852 label='frank, convolved')
854 if dist:
855 ax0_5 = plot_brightness_profile(sol.r, sol.I / 1e10, ax0, dist=dist, c='r', ls=':')
857 u_deproj, v_deproj, vis_deproj = sol.geometry.apply_correction(u, v, vis)
858 baselines = (u_deproj**2 + v_deproj**2)**.5
859 grid = np.logspace(np.log10(min(baselines.min(), sol.q[0])),
860 np.log10(max(baselines.max(), sol.q[-1])), 10**4
861 )
863 for i in range(len(bin_widths)):
864 binned_vis = UVDataBinner(baselines, vis_deproj, weights, bin_widths[i])
865 vis_re = binned_vis.V.real * 1e3
866 vis_err_re = binned_vis.error.real * 1e3
868 if logy:
869 lab=r'Obs.>0, {:.0f} k$\lambda$ bins'.format(bin_widths[i]/1e3)
870 else:
871 lab=r'Obs., {:.0f} k$\lambda$ bins'.format(bin_widths[i]/1e3)
872 plot_vis_quantity(binned_vis.uv / 1e6, vis_re, ax1, c=cs[i],
873 marker=ms[i], ls='None',
874 label=lab)
875 if logy:
876 plot_vis_quantity(binned_vis.uv / 1e6, -vis_re, ax1, c=cs2[i],
877 marker=ms[i], ls='None',
878 label=r'Obs.<0, {:.0f} k$\lambda$ bins'.format(bin_widths[i]/1e3))
880 vis_fit = sol.predict_deprojected(grid).real * 1e3
882 # Take the discrete Hankel transform of the CLEAN profile, using the same
883 # collocation points for the DHT as those in the frank fit
884 Inu_interp = np.interp(
885 sol.r, clean_profile['r'], clean_profile['I'].real) * 1e3
886 clean_DHT = sol.predict_deprojected(grid, I=Inu_interp)
888 if logy:
889 lab = 'frank, >0'
890 lab2 = 'DHT of CLEAN, >0'
891 else:
892 lab = 'frank'
893 lab2 = 'DHT of CLEAN'
894 plot_vis_quantity(grid / 1e6, vis_fit, ax1, c='r', label=lab)
895 plot_vis_quantity(grid / 1e6, clean_DHT, ax1, c='b', label=lab2)
896 if logy:
897 plot_vis_quantity(grid / 1e6, -vis_fit, ax1, c='r', ls='--', label='frank, <0')
898 plot_vis_quantity(grid / 1e6, -clean_DHT, ax1, c='b', ls='--', label='DHT of CLEAN, <0')
900 if logx:
901 ax1.set_xscale('log')
902 ax1.set_xlim(right=1.2 * max(baselines) / 1e6)
904 if logy:
905 ax1.set_yscale('log')
906 ax1.set_ylim(bottom=1e-3)
907 if ylims:
908 ax1.set_ylim(ylims)
910 if MAP_convolved is not None:
911 vmax = max(sol.I.max(), MAP_convolved.max(), clean_profile['I'].max())
912 else:
913 vmax = max(sol.I.max(), clean_profile['I'].max())
914 if stretch == 'asinh':
915 vmin = max(0, min(sol.I))
916 from astropy.visualization.mpl_normalize import simple_norm
917 norm = simple_norm(sol.I, stretch='asinh', asinh_a=asinh_a,
918 min_cut=vmin, max_cut=vmax)
919 elif stretch == 'power':
920 vmin = 0
921 norm = PowerNorm(gamma, vmin, vmax)
922 else:
923 err = ValueError("Unknown 'stretch'. Should be one of 'power' or 'asinh'")
924 raise err
926 plot_2dsweep(sol.r, sol.I, ax=ax2, cmap='inferno', norm=norm,
927 xmax=sol.Rmax, plot_colorbar=True)
928 if MAP_convolved is not None:
929 plot_2dsweep(sol.r, MAP_convolved, ax=ax3, cmap='inferno', norm=norm,
930 xmax=sol.Rmax, plot_colorbar=True)
932 # Interpolate the CLEAN profile onto the frank grid to ensure the CLEAN
933 # swept 'image' has the same pixel resolution as the frank swept 'images'
934 from scipy.interpolate import interp1d
935 interp = interp1d(clean_profile['r'], clean_profile['I'])
936 regrid_I_clean = interp(sol.r)
937 plot_2dsweep(sol.r, regrid_I_clean, ax=ax4, cmap='inferno', norm=norm,
938 xmax=sol.Rmax, plot_colorbar=True)
940 ax0.legend(loc='best')
941 ax1.legend(loc='best', ncol=2)
943 ax0.set_xlabel('r ["]')
944 ax0.set_ylabel(r'Brightness [$10^{10}$ Jy sr$^{-1}$]')
945 ax1.set_xlabel(r'Baseline [M$\lambda$]')
946 ax1.set_ylabel(r'Re(V) [mJy]')
947 ax2.set_xlabel('RA offset ["]')
948 ax3.set_xlabel('RA offset ["]')
949 ax4.set_xlabel('RA offset ["]')
950 ax2.set_ylabel('Dec offset ["]')
952 ax0.set_xlim(right=sol.Rmax)
953 if dist:
954 xlims = ax0.get_xlim()
955 ax0_5.set_xlim(np.multiply(xlims, dist))
957 if logx:
958 ax1.set_xscale('log')
959 if logy:
960 ax1.set_yscale('log')
961 ax1.set_ylim(bottom=1e-3)
963 ax2.set_title('Unconvolved frank profile swept')
964 ax3.set_title('Convolved frank profile swept')
965 ax4.set_title('CLEAN profile swept')
967 ax0.text(.5, .9, 'a)', transform=ax0.transAxes)
968 ax1.text(.5, .9, 'b)', transform=ax1.transAxes)
969 ax2.text(.1, .9, 'c)', c='w', transform=ax2.transAxes)
970 ax3.text(.1, .9, 'd)', c='w', transform=ax3.transAxes)
971 ax4.text(.1, .9, 'e)', c='w', transform=ax4.transAxes)
973 if save_prefix:
974 plt.savefig(save_prefix + '_frank_clean_comparison.png', dpi=300)
975 plt.close()
977 return fig, axes
980def make_multifit_fig(u, v, vis, weights, sols, bin_widths, dist=None,
981 logx=False, force_style=True,
982 save_prefix=None, figsize=(8, 8)
983 ):
984 r"""
985 Produce a figure overplotting multiple fits
987 Parameters
988 ----------
989 u, v : array, unit = :math:`\lambda`
990 u and v coordinates of observations
991 vis : array, unit = Jy
992 Observed visibilities (complex: real + imag * 1j)
993 weights : array, unit = Jy^-2
994 Weights assigned to observed visibilities, of the form
995 :math:`1 / \sigma^2`
996 sols : list of _HankelRegressor objects
997 Reconstructed profile using Maximum a posteriori power spectrum
998 (see frank.radial_fitters.FrankFitter), for each of multiple fits
999 bin_widths : list, unit = \lambda
1000 Bin widths in which to bin the observed visibilities
1001 dist : float, optional, unit = AU, default = None
1002 Distance to source, used to show second x-axis in [AU]
1003 logx : bool, default = False
1004 Whether to plot the visibility distributions in log(baseline)
1005 force_style: bool, default = True
1006 Whether to use preconfigured matplotlib rcParams in generated figure
1007 save_prefix : string, default = None
1008 Prefix for saved figure name. If None, the figure won't be saved
1009 figsize : tuple = (width, height) of figure, unit = inch
1011 Returns
1012 -------
1013 fig : Matplotlib `.Figure` instance
1014 The produced figure, including the GridSpec
1015 axes : Matplotlib `~.axes.Axes` class
1016 The axes of the produced figure
1017 """
1019 logging.info(' Making multifit figure')
1021 with frank_plotting_style_context_manager(force_style):
1022 gs = GridSpec(3, 2, hspace=0)
1023 fig = plt.figure(figsize=figsize)
1025 ax0 = fig.add_subplot(gs[0])
1026 ax1 = fig.add_subplot(gs[2])
1028 ax2 = fig.add_subplot(gs[1])
1029 ax3 = fig.add_subplot(gs[3])
1030 ax4 = fig.add_subplot(gs[5])
1032 ax0.text(.9, .4, 'a)', transform=ax0.transAxes)
1033 ax1.text(.5, .9, 'b)', transform=ax1.transAxes)
1035 ax2.text(.5, .9, 'c)', transform=ax2.transAxes)
1036 ax3.text(.5, .9, 'd)', transform=ax3.transAxes)
1037 ax4.text(.5, .9, 'e)', transform=ax4.transAxes)
1039 axes = [ax0, ax1, ax2, ax3, ax4]
1041 # Assume the fitted geometry and thus deprojected baseline distribution
1042 # is the same for all fits, plotting the common dataset
1043 u_deproj, v_deproj, vis_deproj = sols[0].geometry.apply_correction(u, v, vis)
1044 baselines = (u_deproj**2 + v_deproj**2)**.5
1045 grid = np.logspace(np.log10(min(baselines.min(), sols[0].q[0])),
1046 np.log10(max(baselines.max(), sols[0].q[-1])),
1047 10**4)
1049 for i in range(len(bin_widths)):
1050 binned_vis = UVDataBinner(
1051 baselines, vis_deproj, weights, bin_widths[i])
1052 vis_re = binned_vis.V.real * 1e3
1053 vis_err_re = binned_vis.error.real * 1e3
1055 plot_vis_quantity(binned_vis.uv, vis_re, ax2, c=cs[i],
1056 marker=ms[i], ls='None', label=r'Obs., {:.0f} k$\lambda$ bins'.format(bin_widths[i]/1e3))
1058 plot_vis_quantity(binned_vis.uv, vis_re, ax3, c=cs[i],
1059 marker=ms[i], ls='None',
1060 label=r'Obs.>0, {:.0f} k$\lambda$ bins'.format(bin_widths[i]/1e3))
1061 plot_vis_quantity(binned_vis.uv, -vis_re, ax3, c=cs2[i],
1062 marker=ms[i], ls='None',
1063 label=r'Obs.<0, {:.0f} k$\lambda$ bins'.format(bin_widths[i]/1e3))
1065 # Overplot the multiple fits
1066 for ii in range(len(sols)):
1067 plot_brightness_profile(sols[ii].r, sols[ii].I / 1e10, ax0, c=multifit_cs[ii])#,
1068 # label='alpha = {}, wsmooth = {}'.format(sols[ii].info['alpha'], sols[ii].info['wsmooth'])
1069 if dist and ii == len(sols) - 1:
1070 ax0_5 = plot_brightness_profile(sols[ii].r, sols[ii].I / 1e10, ax0, dist=dist, c=multifit_cs[ii])
1071 xlims = ax0.get_xlim()
1072 ax0_5.set_xlim(np.multiply(xlims, dist))
1074 plot_brightness_profile(sols[ii].r, sols[ii].I / 1e10, ax1, c=multifit_cs[ii])
1076 vis_fit = sols[ii].predict_deprojected(grid).real * 1e3
1077 plot_vis_quantity(grid, vis_fit, ax2, c=multifit_cs[ii])
1079 plot_vis_quantity(grid, vis_fit, ax3, c=multifit_cs[ii])
1080 plot_vis_quantity(grid, -vis_fit, ax3, c=multifit_cs[ii], ls='--')
1082 plot_vis_quantity(sols[ii].q, sols[ii].power_spectrum, ax4, c=multifit_cs[ii])
1084 ax1.set_xlabel('r ["]')
1085 ax0.set_ylabel(r'Brightness [$10^{10}$ Jy sr$^{-1}$]')
1086 ax1.set_ylabel(r'Brightness [$10^{10}$ Jy sr$^{-1}$]')
1087 ax1.set_yscale('log')
1088 ax1.set_ylim(bottom=1e-3)
1090 ax4.set_xlabel(r'Baseline [$\lambda$]')
1091 ax2.set_ylabel('Re(V) [mJy]')
1092 ax3.set_ylabel('Re(V) [mJy]')
1093 ax4.set_ylabel(r'Power [Jy$^2$]')
1094 if logx:
1095 ax2.set_xscale('log')
1096 ax3.set_xscale('log')
1097 ax4.set_xscale('log')
1098 ax3.set_yscale('log')
1099 ax4.set_yscale('log')
1100 ax3.set_ylim(bottom=1e-4)
1102 xlims = ax2.get_xlim()
1103 ax3.set_xlim(xlims)
1104 ax4.set_xlim(xlims)
1106 plt.setp(ax0.get_xticklabels(), visible=False)
1107 plt.setp(ax2.get_xticklabels(), visible=False)
1109 ax0.legend()
1110 ax2.legend()
1111 ax3.legend()
1113 if save_prefix:
1114 plt.savefig(save_prefix + '_frank_multifit.png', dpi=300)
1115 plt.close()
1117 return fig, axes
1120def make_bootstrap_fig(r, profiles, force_style=True,
1121 save_prefix=None, figsize=(8, 6)):
1122 r"""
1123 Produce a figure showing a bootstrap analysis for a Frankenstein fit
1125 Parameters
1126 ----------
1127 r : array, unit = arcsec
1128 Single set of radial collocation points used in all bootstrap fits
1129 profiles : array, unit = Jy / sr
1130 Brightness profiles of all bootstrap fits
1131 force_style: bool, default = True
1132 Whether to use preconfigured matplotlib rcParams in generated figure
1133 save_prefix : string, default = None
1134 Prefix for saved figure name. If None, the figure won't be saved
1135 figsize : tuple = (width, height) of figure, unit = inch
1137 Returns
1138 -------
1139 fig : Matplotlib `.Figure` instance
1140 The produced figure, including the GridSpec
1141 axes : Matplotlib `~.axes.Axes` class
1142 The axes of the produced figure
1143 """
1145 logging.info(' Making bootstrap summary figure')
1147 with frank_plotting_style_context_manager(force_style):
1148 gs = GridSpec(2, 2, hspace=0)
1149 fig = plt.figure(figsize=figsize)
1151 ax0 = fig.add_subplot(gs[0])
1152 ax1 = fig.add_subplot(gs[2])
1154 ax2 = fig.add_subplot(gs[1])
1155 ax3 = fig.add_subplot(gs[3])
1157 axes = [ax0, ax1, ax2, ax3]
1159 ax0.text(.6, .9, 'a) Bootstrap: {} trials'.format(len(profiles)),
1160 transform=ax0.transAxes)
1161 ax1.text(.9, .7, 'b)', transform=ax1.transAxes)
1162 ax2.text(.9, .9, 'c)', transform=ax2.transAxes)
1163 ax3.text(.9, .7, 'd)', transform=ax3.transAxes)
1165 mean_profile = np.mean(profiles, axis=0)
1166 std = np.std(profiles, axis=0)
1168 plot_brightness_profile(r, mean_profile / 1e10, ax1, low_uncer=std / 1e10,
1169 color='r', alpha=.7, label='Stan. dev. of trials')
1170 plot_brightness_profile(r, mean_profile / 1e10, ax3, low_uncer=std / 1e10,
1171 color='r', alpha=.7, label='Stan. dev. of trials')
1173 for i in range(len(profiles)):
1174 plot_brightness_profile(r, profiles[i] / 1e10, ax0, c='k', alpha=.2)
1176 plot_brightness_profile(r, profiles[i] / 1e10, ax2, c='k', alpha=.2)
1178 plot_brightness_profile(r, mean_profile / 1e10, ax1,
1179 c='#35F9E1', label='Mean of trials')
1180 plot_brightness_profile(r, mean_profile / 1e10, ax3,
1181 c='#35F9E1', label='Mean of trials')
1183 ax0.set_ylabel(r'Brightness [$10^{10}$ Jy sr$^{-1}$]')
1184 ax1.set_ylabel(r'Brightness [$10^{10}$ Jy sr$^{-1}$]')
1185 ax2.set_ylabel(r'Brightness [$10^{10}$ Jy sr$^{-1}$]')
1186 ax3.set_ylabel(r'Brightness [$10^{10}$ Jy sr$^{-1}$]')
1187 ax1.set_xlabel('r ["]')
1188 ax3.set_xlabel('r ["]')
1190 ax2.set_yscale('log')
1191 ax3.set_yscale('log')
1192 ax2.set_ylim(bottom=1e-4)
1193 ax3.set_ylim(bottom=1e-4)
1195 if save_prefix:
1196 plt.savefig(save_prefix + '_frank_bootstrap.png', dpi=300,
1197 bbox_inches='tight')
1198 plt.close()
1200 return fig, axes