Coverage for /wheeldirectory/casa-6.7.0-12-py3.10.el8/lib/py/lib/python3.10/site-packages/casatasks/spxfit.py: 89%

27 statements  

« prev     ^ index     » next       coverage.py v7.6.4, created at 2024-11-01 07:19 +0000

1##################### generated by xml-casa (v2) from spxfit.xml #################### 

2##################### 9bd3b69a5c400df6ba6fa4d3f84602ec ############################## 

3from __future__ import absolute_import 

4import numpy 

5from casatools.typecheck import CasaValidator as _val_ctor 

6_pc = _val_ctor( ) 

7from casatools.coercetype import coerce as _coerce 

8from casatools.errors import create_error_string 

9from .private.task_spxfit import spxfit as _spxfit_t 

10from casatasks.private.task_logging import start_log as _start_log 

11from casatasks.private.task_logging import end_log as _end_log 

12from casatasks.private.task_logging import except_log as _except_log 

13 

14class _spxfit: 

15 """ 

16 spxfit ---- Fit a 1-dimensional model(s) to an image(s) or region for determination of spectral index. 

17 

18 --------- parameter descriptions --------------------------------------------- 

19 

20 imagename Name of the input image(s) 

21 box Rectangular region to select in direction plane. Default is to use the entire direction plane. 

22 region Region selection. Default is to use the full image. 

23 chans Channels to use. Default is to use all channels. 

24 stokes Stokes planes to use. Default is to use all Stokes planes. 

25 axis The profile axis. Default: use the spectral axis if one exists, axis 0 otherwise (<0). 

26 mask Mask to use. Default is none. 

27 minpts Minimum number of unmasked points necessary to attempt fit. 

28 multifit If true, fit a profile along the desired axis at each pixel in the specified region. If false, average the non-fit axis pixels and do a single fit to that average profile. Default False. 

29 spxtype Type of function to fit. "plp" = power logarithmic polynomial, "ltp" = logarithmic transformed polynomial. 

30 spxest REQUIRED. Initial estimates as array of numerical values for the spectral index function coefficients. eg [1.5, -0.8] if fitting a plp function thought to be close to 1.5*(x/div)**(-0.8) or [0.4055, -0.8] if fitting an lpt function thought to be close to ln(1.5) - 0.8*ln(x/div). 

31 spxfix Fix the corresponding spectral index function coefficients during the fit. True means hold fixed. 

32 div Divisor (numerical value or quantity) to use in the logarithmic terms of the plp or ltp function. 0 means calculate a useful value on the fly. 

33 spxsol Name of the spectral index function coefficient solution image to write. 

34 spxerr Name of the spectral index function coefficient error image to write. 

35 model Name of model image. Default: do not write the model image (""). 

36 residual Name of residual image. Default: do not write the residual image (""). 

37 wantreturn Should a record summarizing the results be returned? 

38 stretch Stretch the mask if necessary and possible? 

39 logresults Output results to logger? 

40 logfile File in which to log results. Default is not to write a logfile. 

41 append Append results to logfile? Logfile must be specified. Default is to append. False means overwrite existing file if it exists. 

42 sigma Standard deviation array or image name(s). 

43 outsigma Name of output image used for standard deviation. Ignored if sigma is empty. 

44 RETURNS record 

45 

46 --------- examples ----------------------------------------------------------- 

47 

48  

49  

50 This task fits a power logarithmic polynomial or a logarithmic tranformed polynomial to one dimensional profiles for determination of spectral indices. 

51  

52 PARAMETER SUMMARY 

53 imagename Name of the input image(s). More than one image name can be given as an 

54 array, in which case the images are concatenated along the specified axis 

55 and the resultant image is what is used by the fitter. In this case, 

56 all images must have the same dimensions along all axes other than the fit axis. 

57 box Rectangular region to select in direction plane. 

58 Default is to use the entire direction plane. 

59 region Region selection. Default is to use the full image. 

60 chans Channels to use. Default is to use all channels. 

61 stokes Stokes planes to use. Default is to use all Stokes planes. 

62 axis Axis along which to do the fit(s). <0 means use the spectral axis or the 

63 zeroth axis if a spectral axis is not present. 

64 mask Mask to use. Default is none. 

65 stretch Stretch the input mask if necessary and possible? Only used if a mask is specified. 

66  

67 minpts Minimum number of points necessary to attempt a fit. 

68 multifit Fit models at each pixel in region (true) or average profiles and fit a single model (false). 

69 spxtype Type of function to fit. "plp" => power logarithmic polynomial, "ltp" => logarithmic 

70 transformed polynomial. 

71 spxest REQUIRED. Initial estimates as array of numerical values for the spectral index 

72 function coefficients. eg [1.5, -0.8] if fitting a plp function thought to be close to 

73 1.5*(x/div)**(-0.8), or [0.4055, -0.8] if fitting an lpt function thought to be close to 

74 ln(1.5) - 0.8*ln(x/div). 

75 spxfix Fix the corresponding spx function coefficients during the fit. True=>hold fixed 

76 div Divisor (numerical value or quantity) to use in the logarithmic terms of the plp or ltp 

77 function. 0 => calculate a useful value on the fly. 

78 spxsol Name of the function coeffecients solution image to write. 

79 spxerr Name of the function coeffecients error image to write. 

80 model Name of model image to write. 

81 residual Name of residual image to write. 

82 wantreturn If true, return a record summarizing the fit results, if false, return false. 

83 stretch Stretch the mask if necessary and possible? 

84 logresults Output results to logger? 

85 logfile File in which to log results. Default is not to write a logfile. 

86 append Append results to logfile? Logfile must be specified. Default is to append. False means overwrite existing file if it exists. 

87 sigma Standard deviation numerical array, image name (string), or string array of names of images which will be 

88 concatenated to create the sigma image that is used by the fitter. 

89 outsigma Name of output image used for standard deviation. Ignored if sigma is empty. 

90  

91 This application performs a non-linear, least squares fits using the Levenberg-Marquardt algorithm of either a power logarithmic polynomial or a 

92 logarithmic tranformed polynomial to pixel values along a specified axis of an image or images. A description of the fitting algorithm may be found 

93 in AIPS++ Note 224 (http://www.astron.nl/casacore/trunk/casacore/doc/notes/224.html) and in Numerical Recipes by W.H. Press et al., Cambridge 

94 University Press. These functions are most often used for fitting the spectral index and higher order terms of a spectrum. A power logarithmic 

95 polynomial (plp) has the form 

96  

97 y = c0*(x/div)**(c1 + c2*ln(x/div) + c3*ln(x/div)**2 + ... + cn*ln(x/div)**(n - 1)) 

98  

99 and a logarithmic transformed polynomial (ltp) is simply the result of this equation after taking the natural log of both sides so that it has the form 

100  

101 ln(y) = c0 + c1*ln(x/div) + c2*ln(x/div)**2 + ... + cn*ln(x/div)**n 

102  

103 Because the logarithm of the ordinate values must be taken before fitting a logarithmic transformed polynomial, 

104 all non-positive pixel values are effectively masked for the purposes of fitting. 

105  

106 The coefficients of the two forms are equal to each other except that c0 in the second equation is equal to 

107 ln(c0) of the first. In the case of fitting a spectral index, which is traditionally represented as alpha, is 

108 equal to c1. 

109  

110 In both cases, div is a numerical value used to scale abscissa values so they are closer to unity when they are sent to the fitter. This generally 

111 improves the probability that the fit will converge. This parameter may be specified via the div parameter. A value of 0 

112 (the default) indicates that the application should determine a reasonable value for div, which is determined via 

113  

114 div = 10**int(log10(sqrt(min(x)*max(x)))) 

115  

116 where min(x) and max(x) are the minimum and maximum abscissa values, respectively. 

117  

118 So, for example, if S(nu) is proportional to nu**alpha and you expect alpha to be near -0.8 and the value of S(nu) is 1.5 at 

119 1e9 Hz and your image(s) have spectral units of Hz, you would specify spxest=[1.5, -0.8] and div=1e9 when fitting a plp function, 

120 or spxest=[0.405, -0.8] and div=1e9 if fitting an ltp function close to ln(1.5) - 0.8*ln(x/div). 

121  

122  

123 A CAUTIONARY NOTE 

124 Note that the likelihood of getting a reliable solution increases with the number of good data points as well as the goodness 

125 of the initial estimate. It is possible that the first solution found might not be the best one, and 

126 so, if a solution is found, it is recommended that the fit be repeated using the solution of the previous fit as the 

127 initial estimatE for the new fit. This process should be repeated until the solutions from one fit to the next differ only insignificantly. 

128 The convergent solution is very likely the best solution. 

129  

130 AXIS 

131 The axis parameter indicates on which axis profiles should be fit; a value <0 indicates the spectral axis should be used, 

132 or if one does not exist, that the zeroth axis should be used. 

133  

134 MINIMUM NUMBER OF PIXELS 

135 The minpts parameter indicates the minimum number of unmasked pixels that must be present in order for a fit 

136 to be attempted. When multifit=T, positions with too few good points will be masked in any output images. 

137  

138 ONE FIT OF REGION AVERAGE OR PIXEL BY PIXEL FIT 

139 The multifit parameter indicates if profiles should be fit at each pixel in the selected region (true), or if the profiles in that region should be 

140 averaged and the fit done to that average profile (false). 

141  

142 FUNCTION TYPE 

143 Which function to fit is specified in the spxtype parameter. Only two values (case insensitive) are supported. A value of 

144 "plp" indicates that a power logarithmic polynomial should be fit. A value of "ltp" indicates a logarithmic transformed 

145 polynomial should be fit. 

146  

147 INCLUDING STANDARD DEVIATIONS OF PIXEL VALUES 

148 If the standard deviations of the pixel values in the input image are known and they vary in the image (eg they are higher for pixels 

149 near the edge of the band), they can be included in the sigma parameter. This parameter takes either an array or an image name. The 

150 array or image must have one of three shapes: 1. the shape of the input image, 2. the same dimensions as the input image with the lengths 

151 of all axes being one except for the fit axis which must have length corresponding to its length in the input image, or 3. be one 

152 dimensional with lenght equal the the length of the fit axis in the input image. In cases 2 and 3, the array or pixels in sigma will 

153 be replicated such that the image that is ultimately used is the same shape as the input image. The values of sigma must be non-negative. 

154 It is only the relative values that are important. A value of 0 means that pixel should not be used in the fit. Other than that, if pixel 

155 A has a higher standard deviation than pixel B, then pixel A is noisier than pixel B and will receive a lower weight when the fit is done. 

156 The weight of a pixel is the usual 

157  

158 weight = 1/(sigma*sigma) 

159  

160 In the case of multifit=F, the sigma values at each pixel along the fit axis in the hyperplane perpendicular to the fit axis which includes 

161 that pixel are averaged and the resultant averaged standard deviation spectrum is the one used in the fit. Internally, sigma values are normalized 

162 such that the maximum value is 1. This mitigates a known overflow issue. 

163  

164 One can write the normalized standard deviation image used in the fit by specifying its name in outsigma. This image can then be 

165 used as sigma for subsequent runs. 

166  

167 RETURNED DICTIONARY STRUCTURE 

168 The returned dictionary has a (necessarily) complex structure. First, there are keys "xUnit" and "yUnit" whose values are 

169 the abscissa unit and the ordinate unit described by simple strings. Next there are arrays giving a broad overview of the 

170 fit quality. These arrays have the shape of the specified region collapsed along the fit axis with the axis corresponding to the fit 

171 axis having length of 1: 

172  

173 attempted: a boolean array indicating which fits were attempted (eg if too few unmasked points, a fit will not be attempted). 

174 converged: a boolean array indicating which fits converged. False if the fit was not attempted. 

175 valid: a boolean array indicating which solutions fall within the specified valid ranges of parameter space. Any solution for 

176 which a value or error is NaN is automatically marked as invalid. 

177 niter: an int array indicating the number of iterations for each profile, <0 if the fit did not converge 

178 direction: a string array containing the world direction coordinate for each profile 

179  

180 There is a "type" array having number of dimensions equal to the number of dimensions in the above arrays plus one. The shape of 

181 the first n-1 dimensions is the same as the shape of the above arrays. The length of the last dimension is equal to the number of 

182 components fit. The values of this array are all "POWER LOGARITHMIC POLYNOMIAL" or "LOGARITHMIC TRANSFORMED POLYNOMIAL", depending 

183 on which type function was fit. 

184  

185 There will be a subdictionary accessible via the "plp" or "ltp" key (depending on which type of function was fit) which will have 

186 subkeys "solution" and "error" which will each have an array of values. Each of these arrays will have one more dimension than the overview arrays 

187 described above. The shape of the first n-1 dimensions will be the same as the shape of the overview arrays described above, while the 

188 final dimension will have length equal to the number of parameters that were fit. Along this axis will be the 

189 corresponding fit result or associated error (depending on the array's associated key) of the fit. In cases where 

190 the fit was not attempted or did not converge, a value of NAN will be present. 

191  

192 OUTPUT IMAGES 

193 In addition to the returned dictionary, optionally one or more of any combination of output images can be written. 

194 The model and residual parameters indicate the names of the model and residual images to be written; empty values inidcate that these images 

195 should not be written. 

196  

197 The parameters spxsol and spxerr are the names of the solution and error images to write, respectively. In cases 

198 where more than one coefficient are fit, the image names will be appended with an underscore followed by the relevant 

199 coefficient number ("_0", "_1", etc). These images contain the arrays for the associated parameter solutions or errors 

200 described in previous sections. Pixels for which fits were not attempted, did not converge, or converged but have values 

201 of NaN (not a number) or INF (infinity) will be masked as bad. 

202  

203 LPT vs PLP 

204 Ultimately, the choice of which functional form to use in determining the spectral index is up to the user and should be based 

205 on the scientific goals. However, below is a summary of one user's experience and preferences as an example: 

206  

207 If the weights are known or can be determined from the images (eg. the source-free image rms and a fractional calibration error) then I 

208 favor a weighted fit using the non-linear (power-law) model. An unweighted fit using the non-linear model will, in general, give far 

209 too much leverage to large flux values. 

210  

211 If the weights are unknown or will not be considered by the fitting algorithm, then I prefer the log-transformed polynomial model. However, 

212 this does not work well in low signal-to-noise regions. A conservative mask could be created such that only high s/n areas are fit, 

213 but this could hinder many science objectives. 

214  

215 EXAMPLES 

216  

217 res = spxfit(imagename=["im0.im","im1.im"], multifit=true, spxtype="plp", spxest=[0.5,2,0.1], div="1GHz", spxsol="myplpsolutions.im") 

218 

219 

220 """ 

221 

222 _info_group_ = """analysis""" 

223 _info_desc_ = """Fit a 1-dimensional model(s) to an image(s) or region for determination of spectral index.""" 

224 

225 def __call__( self, imagename='', box='', region='', chans='', stokes='', axis=int(-1), mask='', minpts=int(1), multifit=False, spxtype='plp', spxest=[ ], spxfix=[ ], div=[ ], spxsol='', spxerr='', model='', residual='', wantreturn=True, stretch=False, logresults=True, logfile='', append=True, sigma='', outsigma='' ): 

226 schema = {'imagename': {'type': 'cVariant', 'coerce': [_coerce.to_variant]}, 'box': {'type': 'cStr', 'coerce': _coerce.to_str}, 'region': {'type': 'cStr', 'coerce': _coerce.to_str}, 'chans': {'type': 'cStr', 'coerce': _coerce.to_str}, 'stokes': {'type': 'cStr', 'coerce': _coerce.to_str}, 'axis': {'type': 'cInt'}, 'mask': {'type': 'cStr', 'coerce': _coerce.to_str}, 'minpts': {'type': 'cInt'}, 'multifit': {'type': 'cBool'}, 'spxtype': {'type': 'cStr', 'coerce': _coerce.to_str}, 'spxest': {'type': 'cFloatVec', 'coerce': [_coerce.to_list,_coerce.to_floatvec]}, 'spxfix': {'type': 'cBoolVec', 'coerce': [_coerce.to_list,_coerce.to_boolvec]}, 'div': {'type': 'cVariant', 'coerce': [_coerce.to_variant]}, 'spxsol': {'type': 'cStr', 'coerce': _coerce.to_str}, 'spxerr': {'type': 'cStr', 'coerce': _coerce.to_str}, 'model': {'type': 'cStr', 'coerce': _coerce.to_str}, 'residual': {'type': 'cStr', 'coerce': _coerce.to_str}, 'wantreturn': {'type': 'cBool'}, 'stretch': {'type': 'cBool'}, 'logresults': {'type': 'cBool'}, 'logfile': {'type': 'cStr', 'coerce': _coerce.to_str}, 'append': {'type': 'cBool'}, 'sigma': {'anyof': [{'type': 'cStr', 'coerce': _coerce.to_str}, {'type': 'cStrVec', 'coerce': [_coerce.to_list,_coerce.to_strvec]}, {'type': 'cFloatVec', 'coerce': [_coerce.to_list,_coerce.to_floatvec]}, {'type': 'cIntVec', 'coerce': [_coerce.to_list,_coerce.to_intvec]}]}, 'outsigma': {'type': 'cStr', 'coerce': _coerce.to_str}} 

227 doc = {'imagename': imagename, 'box': box, 'region': region, 'chans': chans, 'stokes': stokes, 'axis': axis, 'mask': mask, 'minpts': minpts, 'multifit': multifit, 'spxtype': spxtype, 'spxest': spxest, 'spxfix': spxfix, 'div': div, 'spxsol': spxsol, 'spxerr': spxerr, 'model': model, 'residual': residual, 'wantreturn': wantreturn, 'stretch': stretch, 'logresults': logresults, 'logfile': logfile, 'append': append, 'sigma': sigma, 'outsigma': outsigma} 

228 assert _pc.validate(doc,schema), create_error_string(_pc.errors) 

229 _logging_state_ = _start_log( 'spxfit', [ 'imagename=' + repr(_pc.document['imagename']), 'box=' + repr(_pc.document['box']), 'region=' + repr(_pc.document['region']), 'chans=' + repr(_pc.document['chans']), 'stokes=' + repr(_pc.document['stokes']), 'axis=' + repr(_pc.document['axis']), 'mask=' + repr(_pc.document['mask']), 'minpts=' + repr(_pc.document['minpts']), 'multifit=' + repr(_pc.document['multifit']), 'spxtype=' + repr(_pc.document['spxtype']), 'spxest=' + repr(_pc.document['spxest']), 'spxfix=' + repr(_pc.document['spxfix']), 'div=' + repr(_pc.document['div']), 'spxsol=' + repr(_pc.document['spxsol']), 'spxerr=' + repr(_pc.document['spxerr']), 'model=' + repr(_pc.document['model']), 'residual=' + repr(_pc.document['residual']), 'wantreturn=' + repr(_pc.document['wantreturn']), 'stretch=' + repr(_pc.document['stretch']), 'logresults=' + repr(_pc.document['logresults']), 'logfile=' + repr(_pc.document['logfile']), 'append=' + repr(_pc.document['append']), 'sigma=' + repr(_pc.document['sigma']), 'outsigma=' + repr(_pc.document['outsigma']) ] ) 

230 task_result = None 

231 try: 

232 task_result = _spxfit_t( _pc.document['imagename'], _pc.document['box'], _pc.document['region'], _pc.document['chans'], _pc.document['stokes'], _pc.document['axis'], _pc.document['mask'], _pc.document['minpts'], _pc.document['multifit'], _pc.document['spxtype'], _pc.document['spxest'], _pc.document['spxfix'], _pc.document['div'], _pc.document['spxsol'], _pc.document['spxerr'], _pc.document['model'], _pc.document['residual'], _pc.document['wantreturn'], _pc.document['stretch'], _pc.document['logresults'], _pc.document['logfile'], _pc.document['append'], _pc.document['sigma'], _pc.document['outsigma'] ) 

233 except Exception as exc: 

234 _except_log('spxfit', exc) 

235 raise 

236 finally: 

237 task_result = _end_log( _logging_state_, 'spxfit', task_result ) 

238 return task_result 

239 

240spxfit = _spxfit( ) 

241