Coverage for /wheeldirectory/casa-6.7.0-12-py3.10.el8/lib/py/lib/python3.10/site-packages/casatasks/spxfit.py: 56%
27 statements
« prev ^ index » next coverage.py v7.6.4, created at 2024-10-31 17:39 +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
14class _spxfit:
15 """
16 spxfit ---- Fit a 1-dimensional model(s) to an image(s) or region for determination of spectral index.
18 --------- parameter descriptions ---------------------------------------------
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 [1;42mRETURNS[1;m record
46 --------- examples -----------------------------------------------------------
50 This task fits a power logarithmic polynomial or a logarithmic tranformed polynomial to one dimensional profiles for determination of spectral indices.
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.
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.
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
97 y = c0*(x/div)**(c1 + c2*ln(x/div) + c3*ln(x/div)**2 + ... + cn*ln(x/div)**(n - 1))
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
101 ln(y) = c0 + c1*ln(x/div) + c2*ln(x/div)**2 + ... + cn*ln(x/div)**n
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.
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.
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
114 div = 10**int(log10(sqrt(min(x)*max(x))))
116 where min(x) and max(x) are the minimum and maximum abscissa values, respectively.
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).
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.
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.
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.
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).
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.
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
158 weight = 1/(sigma*sigma)
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.
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.
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:
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
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.
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.
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.
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.
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:
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.
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.
215 EXAMPLES
217 res = spxfit(imagename=["im0.im","im1.im"], multifit=true, spxtype="plp", spxest=[0.5,2,0.1], div="1GHz", spxsol="myplpsolutions.im")
220 """
222 _info_group_ = """analysis"""
223 _info_desc_ = """Fit a 1-dimensional model(s) to an image(s) or region for determination of spectral index."""
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
240spxfit = _spxfit( )