Coverage for frank/debris_fitters.py: 78%

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 

9# the Free Software Foundation, either version 3 of the License, or 

10# (at your option) any later version. 

11# 

12# This program is distributed in the hope that it will be useful, 

13# but WITHOUT ANY WARRANTY; without even the implied warranty of 

14# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 

15# GNU General Public License for more details. 

16# 

17# You should have received a copy of the GNU General Public License 

18# along with this program. If not, see <https://www.gnu.org/licenses/> 

19# 

20"""This module contains methods for fitting a radial brightness profile to a set 

21 of deprojected visibities. Routines in this file assume that the emission is 

22 optically thin with a Gaussian vertical structure. 

23""" 

24from frank.radial_fitters import FourierBesselFitter, FrankFitter 

25 

26class FourierBesselDebrisFitter(FourierBesselFitter): 

27 """ 

28 Fourier-Bessel series optically-thin model for fitting visibilities. 

29 

30 The brightness model is :math:`I(R, z) = I(R) exp(-z^2/2H(R)^2)`, where 

31 :math:`H(R)` is the (known) scale-height. 

32 

33 Parameters 

34 ---------- 

35 Rmax : float, unit = arcsec 

36 Radius of support for the functions to transform, i.e., 

37 f(r) = 0 for R >= Rmax 

38 N : int 

39 Number of collocation points 

40 geometry : SourceGeometry object 

41 Geometry used to deproject the visibilities before fitting 

42 scale_height : function R --> H 

43 Specifies the thickness of disc as a function of radius. Both 

44 units should be in arcsec. 

45 nu : int, default = 0 

46 Order of the discrete Hankel transform (DHT) 

47 block_data : bool, default = True 

48 Large temporary matrices are needed to set up the data. If block_data 

49 is True, we avoid this, limiting the memory requirement to block_size 

50 elements. 

51 block_size : int, default = 10**5 

52 Size of the matrices if blocking is used 

53 verbose : bool, default = False 

54 Whether to print notification messages 

55 """ 

56 

57 def __init__(self, Rmax, N, geometry, scale_height, nu=0, block_data=True, 

58 block_size=10 ** 5, verbose=True): 

59 

60 # All functionality is provided by the base class.  

61 # FourierBesselDebrisFitter is just a sub-set of FourierBesselFitter 

62 super(FourierBesselDebrisFitter, self).__init__( 

63 Rmax, N, geometry, nu=nu, block_data=block_data, 

64 assume_optically_thick=False, scale_height=scale_height, 

65 block_size=block_size, verbose=verbose 

66 ) 

67 

68class FrankDebrisFitter(FrankFitter): 

69 """ 

70 Fit a Gaussian process model using the Discrete Hankel Transform of 

71 Baddour & Chouinard (2015). 

72 

73 The brightness model is :math:`I(R, z) = I(R) exp(-z^2/2H(R)^2)`, where 

74 :math:`H(R)` is the (known) scale-height. 

75 

76 The GP model is based upon Oppermann et al. (2013), which use a maximum 

77 a posteriori estimate for the power spectrum as the GP prior for the 

78 real-space coefficients. 

79 

80 Parameters 

81 ---------- 

82 Rmax : float, unit = arcsec 

83 Radius of support for the functions to transform, i.e., f(r) = 0 for 

84 R >= Rmax. 

85 N : int 

86 Number of collaction points 

87 geometry : SourceGeometry object 

88 Geometry used to deproject the visibilities before fitting 

89 scale_height : function R --> H 

90 Specifies the thickness of disc as a function of radius. Both 

91 units should be in arcsec. 

92 nu : int, default = 0 

93 Order of the discrete Hankel transform, given by J_nu(r) 

94 block_data : bool, default = True 

95 Large temporary matrices are needed to set up the data. If block_data 

96 is True, we avoid this, limiting the memory requirement to block_size 

97 elements 

98 block_size : int, default = 10**5 

99 Size of the matrices if blocking is used 

100 alpha : float >= 1, default = 1.05 

101 Order parameter of the inverse gamma prior for the power spectrum 

102 coefficients 

103 p_0 : float >= 0, default = None, unit=Jy^2 

104 Scale parameter of the inverse gamma prior for the power spectrum 

105 coefficients. If not provided p_0 = 1e-15 (method="Normal") or  

106 1e-35 (method="LogNormal") will be used. 

107 weights_smooth : float >= 0, default = 1e-4 

108 Spectral smoothness prior parameter. Zero is no smoothness prior 

109 tol : float > 0, default = 1e-3 

110 Tolerence for convergence of the power spectrum iteration 

111 method : string, default="Normal" 

112 Model used for the brightness reconstrution. This must be one of 

113 "Normal" of "LogNormal". 

114 I_scale : float, default = 1e5, unit= Jy/Sr 

115 Brightness scale. Only used in the LogNormal model. Notet the  

116 LogNormal model produces I(Rmax) = I_scale. 

117 max_iter: int, default = 2000 

118 Maximum number of fit iterations 

119 check_qbounds: bool, default = True 

120 Whether to check if the first (last) collocation point is smaller 

121 (larger) than the shortest (longest) deprojected baseline in the dataset 

122 store_iteration_diagnostics: bool, default = False 

123 Whether to store the power spectrum parameters and brightness profile 

124 for each fit iteration 

125 verbose: 

126 Whether to print notification messages 

127 convergence_failure: string, default = 'raise' 

128 Decide what to do when the frank model does not converge within max_iter. 

129 Should be one of: 

130 'raise' : raise an error 

131 'warn' : print a warning message and continue 

132 'ignore' : Ignore the error. 

133 

134 References 

135 ---------- 

136 Baddour & Chouinard (2015) 

137 DOI: https://doi.org/10.1364/JOSAA.32.000611 

138 Oppermann et al. (2013) 

139 DOI: https://doi.org/10.1103/PhysRevE.87.032136 

140 """ 

141 

142 def __init__(self, Rmax, N, geometry, scale_height, nu=0, block_data=True, 

143 block_size=10 ** 5, alpha=1.05, p_0=None, weights_smooth=1e-4, 

144 tol=1e-3, method='Normal', I_scale=1e5, max_iter=2000, check_qbounds=True, 

145 store_iteration_diagnostics=False, verbose=True, 

146 convergence_failure='raise'): 

147 

148 # All functionality is provided by the base class. FrankDebrisFitter is just a  

149 # sub-set of FrankFitter 

150 super(FrankDebrisFitter, self).__init__( 

151 Rmax, N, geometry, nu=nu, block_data=block_data, 

152 block_size=block_size, alpha=alpha, p_0=p_0, weights_smooth=weights_smooth, 

153 tol=tol, method=method, I_scale=I_scale, max_iter=max_iter, 

154 check_qbounds=check_qbounds, store_iteration_diagnostics=store_iteration_diagnostics, 

155 assume_optically_thick=False, scale_height=scale_height, verbose=verbose, 

156 convergence_failure=convergence_failure)