casa
$Rev:20696$
|
00001 # 00002 # This file was generated using xslt from its XML file 00003 # 00004 # Copyright 2008, Associated Universities Inc., Washington DC 00005 # 00006 import sys 00007 import os 00008 from casac import * 00009 import string 00010 import time 00011 import inspect 00012 import gc 00013 import numpy 00014 from odict import odict 00015 from task_specfit import specfit 00016 from task_specfit import casalog 00017 00018 class specfit_pg_: 00019 __name__ = "specfit" 00020 00021 def __init__(self) : 00022 self.__bases__ = (specfit_pg_,) 00023 self.__doc__ = self.__call__.__doc__ 00024 00025 00026 def __call__(self, imagename=None, box=None, region=None, chans=None, stokes=None, axis=None, mask=None, ngauss=None, poly=None, estimates=None, minpts=None, multifit=None, model=None, residual=None, amp=None, amperr=None, center=None, centererr=None, fwhm=None, fwhmerr=None, integral=None, integralerr=None, wantreturn=None, stretch=None, logresults=None, pampest=None, pcenterest=None, pfwhmest=None, pfix=None, gmncomps=None, gmampcon=None, gmcentercon=None, gmfwhmcon=None, gmampest=None, gmcenterest=None, gmfwhmest=None, gmfix=None, logfile=None, append=None, pfunc=None, goodamprange=None, goodcenterrange=None, goodfwhmrange=None, sigma=None, outsigma=None, async=None): 00027 00028 """Fit 1-dimensional gaussians and/or polynomial models to an image or image region 00029 00030 This task simultaneously fits one or more gaussian singlets lorentzian singlets, gaussian multiplets, and/or a polynomial to one dimensional profiles. 00031 00032 ARAMETER SUMMARY 00033 imagename Name of the input (CASA, FITS, MIRIAD) image 00034 box Direction plane box specification, "blcx, blcy, trcx, trcy". Only one box 00035 may be specified. If not specified, region is used if specified. If region 00036 is also not specified, entire directional plane unioned with any chans and 00037 stokes specification determines the region. 00038 region Region of interest. See help par.region for possible specifications. 00039 chans Optional contiguous frequency channel number specification. Not used if 00040 region is specified. Default is all channels. 00041 stokes Contiguous stokes planes specification. Not used if region is specified. 00042 Default is all stokes. 00043 axis Axis along which to do the fit(s). <0 means use the spectral axis or the 00044 zeroth axis if a spectral axis is not present. 00045 mask Mask to use. See help par.mask. Default is none. 00046 stretch Stretch the input mask if necessary and possible? Only used if a mask is specified. 00047 See help par.stretch. 00048 ngauss Maximum number of gaussians to fit. 00049 poly Order of polynomial to fit. <0 means do not fit a polynomial. 00050 estimates Name of file containing initial gaussian estimates. 00051 minpts Minimum number of points necessary to attempt a fit. 00052 multifit Fit models at each pixel in region (true) or average profiles and fit a single model (false). 00053 model Name of model image to write. 00054 residual Name of residual image to write. 00055 amp Name of amplitude solution image. Default: do not write the image ("") 00056 amperr Name of amplitude solution error image. Default: do not write the image ("") 00057 center Name of center solution image. Default: do not write the image ("") 00058 centererr Name of center solution error image. Default: do not write the image ("") 00059 fwhm Name of fwhm solution image. Default: do not write the image ("") 00060 fwhmerr Name of fwhm solution error image. Default: do not write the image ("") 00061 integral Name of integral solution image. Default: do not write the image ("") 00062 integralerr Name of integral solution error image. Default: do not write the image ("") 00063 wantreturn If true, return a record summarizing the fit results, if false, return false. 00064 stretch Stretch the mask if necessary and possible? See help par.stretch 00065 logresults Output results to logger? 00066 pampest Initial estimate of PCF profile (gaussian or lorentzian) amplitudes. 00067 pcenterest Initial estimate PCF profile centers, in pixels. 00068 pfwhmest Initial estimate PCF profile FWHMs, in pixels. 00069 pfix PCF profile parameters to fix during fit. 00070 pfunc PCF singlet functions to fit. "gaussian" or "lorentzian" (minimal match supported). Unspecified means all gaussians. 00071 gmncomps Number of components in each Gaussian multiplet to fit. 00072 gmampcon The amplitude ratio constraints for non-reference components to reference component in gaussian multiplets. 00073 gmcentercon The center offset constraints (in pixels) for non-reference components to reference component in gaussian multiplets. 00074 gmfwhmcon The FWHM ratio constraints for non-reference components to reference component in gaussian multiplets. 00075 gmampest Initial estimate of individual gaussian amplitudes in gaussian multiplets. 00076 gmcenterest Initial estimate of individual gaussian centers in gaussian multiplets, in pixels. 00077 gmfwhmest Initial estimate of individual gaussian FWHMss in gaussian multiplets, in pixels. 00078 gmfix Parameters of individual gaussians in gaussian multiplets to fix during fit. 00079 logfile File in which to log results. Default is not to write a logfile. 00080 append Append results to logfile? Logfile must be specified. Default is to append. False means overwrite existing file if it exists. 00081 goodamprange Acceptable amplitude solution range. 0 => all amplitude solutions are acceptable. [0.0] => All amplitude solutions are acceptable. 00082 goodcenterrange Acceptable center solution range in pixels relative to region start. [0.0] => all center solutions are acceptable. 00083 goodfwhmrange Acceptable FWHM solution range in pixels. [0.0] => all FWHM solutions are acceptable. 00084 sigma Standard deviation array or image name. 00085 outsigma Name of output image used for standard deviation. Ignored if sigma is empty. 00086 00087 This task simultaneously fits one or more gaussian singlets, one or more lorentzian singlets, one or more gaussian multiplets, 00088 and/or a polynomial to one dimensional profiles. A gaussian/lorentzian singlet is a gaussian/lorentzian whose parameters (amplitude, 00089 center position, and width) are all independent from any other feature that may be simultaneously fit. A gaussian multiplet is a set of two or 00090 more gaussian lines in which at least one (and possibly two or three) parameter of each line is dependent on the parameter of another, 00091 single (reference) profile in the multiplet. For example, one can specify a doublet in which the amplitude of the first line is 0.6 times the 00092 amplitude of the zeroth line and/or the center of the first line is 20 pixels from the center of the zeroth line, and/or the fwhm of the first 00093 line is identical (in pixels) to that of the zeroth line. There is no limit to the number of components one can specify in a multiplet 00094 (except of course that the number of parameters to be fit should be significantly less than the number of data points), but there can be only 00095 a single reference profile in a multiplet to which to tie constraints of parameters of the other profiles in the set. 00096 00097 AXIS 00098 The axis parameter indicates on which axis profiles should be fit; a value <0 indicates the spectral axis should be used, or if one does not exist, 00099 that the zeroth axis should be used. 00100 00101 MINIMUM NUMBER OF PIXELS 00102 The minpts parameter indicates the minimum number of unmasked pixels that must be present in order for a fit 00103 to be attempted. When multifit=T, positions with too few good points will be masked in any output images. 00104 00105 ONE FIT OF REGION AVERAGE OR PIXEL BY PIXEL FIT 00106 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 00107 averaged and the fit done to that average profile (false). 00108 00109 POLYNOMIAL FITTING 00110 The order of the polynomial to fit is specified only via the poly parameter. If poly<0, no polynomial will be fit. No initial estimates of 00111 coefficients can be specified; these are determined automatically. 00112 00113 GAUSSIAN SINGLET FITTING 00114 In the absence of an estimates file and no estimates being specified by the p*est parameters, and gmncomps=0 or is empty, the ngauss parameter 00115 indicates the maximum number of gaussian singlets that should be fit. The initial estimates of the parameters for these gaussians will be attempted 00116 automatically in this case. If it deems appropriate, the fitter will fit fewer than this number. In the case where an estimates file is supplied, 00117 ngauss is ignored (see below). ngauss is also ignored if the p*est parameters are specified or if gmncomps is not an empty array or, if an integer, 00118 is greater than zero. If estimates is not specified or the p*est parameters are not specified and ngauss=0, gmncomps is empty or 0, and poly<0, 00119 an error will occur as this indicates there is nothing to fit. 00120 00121 One can specify initial estimates of gaussian singlet parameters via an estimates file or the pampest, pcenterest, pfwhmest, and optionally, the 00122 pfix parameters. The latter is the recommended way to specify these estimates as support for estimates files may be deprecated in the future. No matter 00123 which option is used, an amplitude initial estimate must always be nonzero. A negative fwhm estimate will be silently changed to positve. 00124 00125 SPECIFYING INITIAL ESTIMATES FOR GAUSSIAN AND LORENTZIAN SINGLETS (RECOMMENDED METHOD) 00126 One may specify initial estimates via the pampest, pcenterest, and pfwhmest parameters. In the case of a single gaussian or lorentzian singlet, 00127 these parameters can be numbers. pampest must be specified in image brightness units, pcenterest must be given in the number of pixels from the 00128 zeroth pixel, and pfwhmest must be given in pixels. Optionally pfix can be specified and in the case of a single gaussian or lorentzian singlet 00129 can be a string. In it is coded which parameters should be held constant during the fix. Any combination of "p" (amplitude), "c" (center), or "f" 00130 (fwhm) is allowed; eg pfix="pc" means fix both the amplitude and center during the fit. In the case of more than one gaussian and/or lorentzian 00131 singlets, these parameters must be specified as arrays of numbers. The length of the arrays indicates the number of singlets to fit and must be 00132 the same for all the p*est parameters. 00133 00134 If no parameters are to be fixed for any of the singlets, pfix can be set to the empty string. However, if at least one parameter of one singlet 00135 is to be fixed, pfix must be an array of strings and have a length equal to the p*est arrays. Singlets which are not to have any parameters fixed 00136 should be represented as an empty string in the pfix array. So, for example, if one desires to fit three singlets and fix the fwhm of the middle 00137 one, one must specify pfix=["", "f", ""], the empty strings indicating no parameters of the zeroth and second singlet should be held constant. 00138 00139 In the case of multifit=True, the initial estimates, whether from the p*est parameters or from a file (see below), will be applied to the location 00140 of the first fit. This is normally the bottom left corner of the region selected. If masked, not enough good points to perform a fit, or the 00141 attempted fit fails, the fitting proceeds to the next pixel with the pixel value of the lowest numbered axis changing the fastest. Once a 00142 successful fit has been performed, subsequent fits will use the results of a fit for a nearest pixel for which a previous fit was successful as the 00143 initial estimate for the parameters at the current location. The fixed parameter string will be honored for every fit performed when multifit=True. 00144 00145 One specifies what type of PCF profile to fit via the pfunc parameter. A PCF function is one that can be parameterized by a peak, center, and FWHM, 00146 as both gaussian and lorentzian singlets can. If all singlets to be fit are gaussians, one can set pfunc equal to the empty string and all snglets 00147 will be assumed to be gaussians. If at least one lorentzian is to be fit, pfunc must be specified as a string (in the case of a single singlet) or 00148 an array of strings (in the case of multiple singlets). The position of each string corresponds to the positions of the initial estimates in the 00149 p*est and pfix arrays. Minimal match ("g", "G", "l", or "L") is supported. So, if one wanted to simultaneously fit two gaussian and two lorentzian 00150 singlets, the zeroth and last of which were lorentzians, one would specify pfunc=["L", "G", "G", "L"]. 00151 00152 ESTIMATES FILE FOR GAUSSIAN SINGLETS (NONRECOMMENDED METHOD) 00153 Initial estimates for gaussian singlets can be specified in an estimates file. Estimates files may be deprecated in the future in favor of the 00154 p*est parameters, so it is recommended users use those parameters instead. If an estimates file is desired to be used, the p*est parameters 00155 must be 0 or empty and mgncomps must be 0 or empty. Only gaussian singlets can be specified in an estimates file. If one desires to fit one or 00156 more gaussian multiplets and/or one or more lorentzian singlets simultaneously, the p*est parameters must be used to specify the initial parameters 00157 of all gaussian singlets to fit; one cannot use an estimates file in this case. If an estimates file is specified, a polynomial 00158 can be fit simultaneously by specifying the poly parameter. The estimates file must contain initial estimates of parameters 00159 for all gaussian singlets to be fit. The number of gaussian singlets to fit is gotten from the number of estimates in the file. The file can contain 00160 comments which are indicated by a "#" at the beginning of a line. All non-comment lines will be interpreted as initial estimates. The 00161 format of such a line is 00162 00163 [peak intensity], [center], [fwhm], [optional fixed parameter string] 00164 00165 The first three values are required and must be numerical values. The peak intensity must be expressed in image brightness units, while the 00166 center must be specified in pixels offset from the zeroth pixel, and fwhm must be specified in pixels. The fourth value is optional and if present, 00167 represents the parameter(s) that should be held constant during the fit. Any combination of the characters 'p' (peak), 'c' (center), and 'f' (fwhm) are 00168 permitted, eg "fc" means hold the fwhm and the center constant during the fit. Fixed parameters will have no error associated with them. Here is an 00169 example file: 00170 00171 # estimates file indicating that two gaussians should be fit 00172 # first guassian estimate, peak=40, center at pixel number 10.5, fwhm = 5.8 pixels, all parameters allowed to vary during 00173 # fit 00174 40, 10.5, 5.8 00175 # second gaussian, peak = 4, center at pixel number 90.2, fwhm = 7.2 pixels, hold fwhm constant 00176 4, 90.2, 7.2, f 00177 # end file 00178 00179 GAUSSIAN MULTIPLET FITTING 00180 Any number of gaussian multiplets, each containing any number of two or more components, can be simultaneously fit, optionally with a 00181 polynomial and/or any number of gaussian and/or lorentzian singlets, the only caveat being that the number of parameters to be fit should be 00182 significantly less than the number of data points. The gmncomps parameter indicates the number of multiplets to fit and the number of 00183 components in each multiplet. In the case of a single multiplet, an integer (>1) can be specified. For example, mgncomps=4 means fit a 00184 single quadruplet of gaussians. In the case of 2 or more multiplets, and array of integers (all >1) must be specified. For example, 00185 gmncomps=[2, 4, 3] means 3 seperate multiples are to be fit, the zeroth being a doublet, the first being a quadruplet, and the second 00186 being a triplet. 00187 00188 Initial estimates of all gaussians in all multiplets are specified via the gm*est parameters which must be arrays of numbers. The order 00189 starts with the zeroth component of the zeroth multiplet to the last component of the zeroth multiplet, then the zeroth component of 00190 the first multiplet to the last compoenent of the first multiplet, etc to the zeroth component of the last multiplet to the last 00191 element of the last multiplet. The zeroth element of a multiplet is defined as the reference component of that multiplet and has the special 00192 significance that it is the profile to which all constraints of all other profiles in that multiplet are referenced (see below). So, 00193 in our example of gmncomps=[2, 4, 3], gmampest, gmcenterest, and gmfwhmest must each be nine (the total number of individual gaussian 00194 profiles summed over all multiplets) element arrays. The zeroth, second, and sixth elements represent parameters of the reference profiles 00195 in the zeroth, first, and second multiplet, respectively. 00196 00197 The fixed relationships between the non-reference profile(s) and the reference profile of a multiplet are specified via the gmampcon, 00198 gmcentercon, and gmfwhmcon parameters. At least one, and any combination, of constraints can be specified for any non-reference 00199 component of a multiplet. The amplitude ratio of a non-reference line to that of the reference line is set in gmampcon. The ratio of 00200 the fwhm of a non-reference line to that of the reference line is set in gmfwhmcon. The offset in pixels of the center position of 00201 a non-reference line to that of the reference line is set in gmcentercon. In the case where a parameter is not constrained for any 00202 non-reference line of any multiplet, the value of the associated parameter must be 0. In the case of 00203 a single doublet, a constraint may be specified as a number or an array of a single number. For example, mgncomps=2 and gmampcon=0.65 00204 and gmcentercon=[32.4] means there is a single doublet to fit where the amplitude ratio of the first to the zeroth line is constained 00205 to be 0.65 and the center of the first line is constrained to be offset by 32.4 pixels from the center of the zeroth line. In cases 00206 of a total of three or more gaussians, the constraints parameters must be specified as arrays with lengths equal to the total number 00207 of gaussians summed over all multiplets minus the number of reference lines (one per multiplet, or just number of multiplets, since 00208 reference lines cannot be constrained by themselves). In the cases where an array must be specified but a component in that array 00209 does not have that constraint, 0 should be specified. Here's an example 00210 00211 gmncomps=[2, 4, 3] 00212 gmampcon= [ 0 , 0.2, 0 , 0.1, 4.5, 0 ] 00213 gcentercon=[24.2, 45.6, 92.7, 0 , -22.8, -33.5] 00214 gfwhmcon="" 00215 00216 In this case we have our previous example of one doublet, one quadruplet, and one triplet. The first component of the doublet has the constraint 00217 that its center is offset by 24.2 pixels from the zeroth (reference) component. The first component of the quadruplet is constrained to have 00218 an amplitude of 0.2 times that of the quadruplet's zeroth component and its center is constrained to be offset by 45.6 pixels from the 00219 reference component. The second component of the quadruplet is constained to have its center offset by 92.7 pixels from the associated 00220 reference component and the third component is constrained to have an amplitude of 0.1 times that of the associated reference component. 00221 The first component of the triplet is constrained to have an amplitude of 4.5 times that of its associated reference component and its center 00222 is constrained to be offset by -22.8 pixels from the reference component's center. The second component of the triplet is constrained to have 00223 its center offset by -33.5 pixels from the center of the reference component. No lines have FWHM constraints, so the empty string can be given 00224 for that parameter. Note that using 0 to indicate no constraint for line center means that one cannot specify a line centered at the same 00225 position as the reference component but having a different FWHM from the reference component. If you must specify this very unusual case, 00226 try using a very small positive (or even negative) value for the center constraint. 00227 00228 Note that when a parameter for a line is constrained, the corresponding value for that component in the corresponding gm*est array is 00229 ignored and the value of the constrained parameter is automatically used instead. So let's say, for our example above, we had specified 00230 the following estimates: 00231 00232 gmampest = [ 1, .2, 2, .1, .1, .5, 3, 2, 5] 00233 gmcenterest = [20, 10 , 30, 45.2, 609 , -233, 30, -859, 1] 00234 00235 Before any fitting is done, the constraints would be taken into account and these arrays would be implicitly rewritten as: 00236 00237 gmampest = [ 1, .2, 2, .4, .1, .2, 3, 13.5, 5 ] 00238 gmcenterest = [20, 44.2, 30, 75.6, 127.7, -233, 30, 7.2, -3.5] 00239 00240 The value of gmfwhmest would be unchanged since there are no FWHM constraints in this example. 00241 00242 In addition to be constrained by values of the reference component, parameters of individual components can be fixed. Fixed parameters 00243 are specified via the gmfix parameter. If no parameters are to be fixed, gmfix can be specified as the empty string or a zero element 00244 array. In the case where any parameter is to be fixed, gmfix must be specified as an array of strings with length equal to the total number of 00245 components summed over all multiplets. These strings encode which parameters to be fixed for the corresponding components. If 00246 a component is to have no parameters fixed, an empty string is used. In other cases one or more of any combination of parameters can 00247 be fixed using "p", "c", and/or "f" described above for fixing singlet parameters. There are a couople of special cases 00248 to be aware of. In the case where a non-reference component parameter is constrained and the corresponding reference component parameter is 00249 set as fixed, that parameter in the non-reference parameter will automatically be fixed even if it was specified not to be fixed in 00250 the gmfix array. This is the only way the constraint can be honored afterall. In the converse case of when a constrained parameter of a 00251 non-reference component is specified as fixed, but the corresponding parameter in the reference component is not specified to be fixed, 00252 an error will occur. Fixing an unconstrained parameter in a non-reference component is always legal as is fixing any combination of 00253 parameters in a reference component (with the above caveat that corresponding constrained parameters in non-reference components will 00254 be silently held fixed as well). 00255 00256 The same rules that apply to singlets when multifit=True apply to multiplets. 00257 00258 LIMITING RANGES FOR SOLUTION PARAMETERS 00259 In cases of low (or no) signal to noise spectra, it is still possible for the fit to converge, but often to a 00260 nonsensical solution. The astronomer can use her knowledge of the source to filter out obviously bogus solutions. 00261 One can limit the ranges of solution parameters to known "good" values via the goodamprange, goodcenterrange, and goodfwhmrange 00262 parameters. Any combination can be specified and the limit constraints will be ANDed together. The ranges apply to all PCF components 00263 that might be fit; choosing ranges on a component by component basis is not supported. If specified, 00264 an array of exactly two numerical values must be given to indicate the range of acceptable solution values for 00265 that parameter. goodamprange is expressed in terms of image brightness units. goodcenterrange is expressed in terms of pixels 00266 from the zeroth pixel in the specified region. goodfwhmrange is expressed in terms of pixels (only non-negative values should be 00267 given for FWHM range endpoints). In the case of a multiple-PCF fit, if any of the corresponding solutions are outside the specified 00268 ranges, the entire solution is considered to be invalid. 00269 00270 In addition, solutions for which the absolute value of the ratio of the amplitude error to the amplitude exceeds 100 or the 00271 ratio of the FWHM error to the FWHM exceeds 100 are automatically marked as invalid. 00272 00273 INCLUDING STANDARD DEVIATIONS OF PIXEL VALUES 00274 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 00275 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 00276 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 00277 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 00278 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 00279 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. 00280 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 00281 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. 00282 The weight of a pixel is the usual 00283 00284 weight = 1/(sigma*sigma) 00285 00286 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 00287 that pixel are averaged and the resultant averaged standard deviation spectrum is the one used in the fit. Internally, sigma values are normalized 00288 such that the maximum value is 1. This mitigates a known overflow issue. 00289 00290 One can write the normalized standard deviation image used in the fit but specifying its name in outsigma. This image can then be 00291 used as sigma for subsequent runs. 00292 00293 RETURNED DICTIONARY STRUCTURE 00294 The dictionary returned (if wantreturn=True) has a (necessarily) complex structure. First, there are keys "xUnit" and "yUnit" whose values are 00295 the abscissa unit and the ordinate unit described by simple strings. Next there are arrays giving a broad overview of the 00296 fit quality. These arrays have the shape of the specified region collapsed along the fit axis with the axis corresponding to the fit 00297 axis having length of 1: 00298 00299 attempted: a boolean array indicating which fits were attempted (eg if too few unmasked points, a fit will not be attempted). 00300 converged: a boolean array indicating which fits converged. False if the fit was not attempted. 00301 valid: a boolean array indicating which solutions fall within the specified valid ranges of parameter space 00302 niter: an int array indicating the number of iterations for each profile, <0 if the fit did not converge 00303 ncomps: the number of components (gaussian singlets + lorentzian singlets + gaussian multiplets + polynomial) fit for the profile, 00304 <0 if the fit did not converge 00305 00306 There is a "type" array having number of dimensions equal to the number of dimensions in the above arrays plus one. The shape of 00307 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 00308 components fit. The values of this array are strings describing the components that were fit at each possition ("POLYNOMIAL", 00309 "GAUSSIAN" in the case of gaussian singlets, "LORENTZIAN" in the case of lorentzian singlets, and ""GAUSSIAN MULTPLET"). 00310 00311 If any gaussian singlets were fit, there will be a subdictionary accessible via the "gs" key which will have subkeys "amp", "ampErr", "center", 00312 "centerErr", "fwhm", "fwhmErr, "integral", and "integralErr". Each of these arrays will have one more dimension than the overview arrays described 00313 above. The shape of the first n-1 dimensions will be the same as the shape of the arrays described above, while the final dimension will 00314 have length equal to the maximum number of gaussian singlets that were fit. Along this axis will be the 00315 corresponding fit result or associated error (depending on the array's associated key) of the fit for that singlet component number. In cases where 00316 the fit did not converge, or that particular component was excluded from the fit, a value of NAN will be present. 00317 00318 If any lorentzian singlets were fit, their solutions will be accessible via the "ls" key. These arrays follow the same rules 00319 as the "gs" arrays described above. 00320 00321 If any gaussian multiplets were fit, there will be subdictionaries accessible by keys "gm0", "gm1", ..., "gm{n-1}" where n is the number of gaussian 00322 muliplets that were fit. Each of these dictionaries will have the same arrays described above for gaussian singlets. The last dimension 00323 will have length equal to the number of components in that particular multiplet. Each pixel along the last axis will be the parameter solution 00324 value or error for that component number in the multiplet, eg the zeroth pixel along that axis contains 00325 the parameter solution or error for the reference component of the multiplet. 00326 00327 The polynomial coefficient solutions and errors are not returned, although they are logged. 00328 00329 OUTPUT IMAGES 00330 In addition to the returned dictionary, optionally one or more of any combination of output images can be written. 00331 The model and residual parameters indicate the names of the model and residual images to be written; blank values inidcate that these images 00332 should not be written. 00333 00334 One can also write none, any or all of the solution and error images for gaussian singlet, lorentzian singlet, and gaussian multiplet fits 00335 via the parameters amp, amperr, center, centererr, fwhm, fwhmerr, integral, and integralerr when doing multi-pixel fits. These images simply 00336 contain the arrays described for the associated parameter solutions or errors described in previous sections. In the case of lorentzian 00337 singlets, "_ls" is appended to the image names, in the case of gaussian multiplets, "_gm0", "_gm1", etc are appended to the image names to 00338 distinguish each multiplet. Pixels for which fits were not attempted or did not converge will be masked as bad. The last axis of these images 00339 is a linear axis and repesents component number (and is named accordingly). 00340 00341 Writing analogous images for polynomial coefficients is not supported. 00342 00343 EXAMPLE 00344 res = specif(imagename="myspectrum.im", ngauss=2, box="3,3,4,5", poly=2, multifit=true, wantreturn=True) 00345 00346 """ 00347 a=inspect.stack() 00348 stacklevel=0 00349 for k in range(len(a)): 00350 if (string.find(a[k][1], 'ipython console') > 0) or (string.find(a[k][1], '<string>') >= 0): 00351 stacklevel=k 00352 break 00353 myf=sys._getframe(stacklevel).f_globals 00354 myf['__last_task'] = 'specfit' 00355 myf['taskname'] = 'specfit' 00356 ### 00357 myf['update_params'](func=myf['taskname'],printtext=False) 00358 ### 00359 ### 00360 #Handle globals or user over-ride of arguments 00361 # 00362 function_signature_defaults=dict(zip(self.__call__.func_code.co_varnames,self.__call__.func_defaults)) 00363 useLocalDefaults = False 00364 00365 for item in function_signature_defaults.iteritems(): 00366 key,val = item 00367 keyVal = eval(key) 00368 if (keyVal == None): 00369 #user hasn't set it - use global/default 00370 pass 00371 else: 00372 #user has set it - use over-ride 00373 if (key != 'self') : 00374 useLocalDefaults = True 00375 #myf[key]=keyVal 00376 00377 myparams = {} 00378 if useLocalDefaults : 00379 for item in function_signature_defaults.iteritems(): 00380 key,val = item 00381 keyVal = eval(key) 00382 exec('myparams[key] = keyVal') 00383 if (keyVal == None): 00384 exec('myparams[key] = '+ key + ' = self.itsdefault(key)') 00385 keyVal = eval(key) 00386 if(type(keyVal) == dict) : 00387 exec('myparams[key] = ' + key + ' = keyVal[len(keyVal)-1][\'value\']') 00388 00389 else : 00390 uselessvariable = None 00391 myparams['imagename'] = imagename = myf['imagename'] 00392 myparams['box'] = box = myf['box'] 00393 myparams['region'] = region = myf['region'] 00394 myparams['chans'] = chans = myf['chans'] 00395 myparams['stokes'] = stokes = myf['stokes'] 00396 myparams['axis'] = axis = myf['axis'] 00397 myparams['mask'] = mask = myf['mask'] 00398 myparams['ngauss'] = ngauss = myf['ngauss'] 00399 myparams['poly'] = poly = myf['poly'] 00400 myparams['estimates'] = estimates = myf['estimates'] 00401 myparams['minpts'] = minpts = myf['minpts'] 00402 myparams['multifit'] = multifit = myf['multifit'] 00403 myparams['model'] = model = myf['model'] 00404 myparams['residual'] = residual = myf['residual'] 00405 myparams['amp'] = amp = myf['amp'] 00406 myparams['amperr'] = amperr = myf['amperr'] 00407 myparams['center'] = center = myf['center'] 00408 myparams['centererr'] = centererr = myf['centererr'] 00409 myparams['fwhm'] = fwhm = myf['fwhm'] 00410 myparams['fwhmerr'] = fwhmerr = myf['fwhmerr'] 00411 myparams['integral'] = integral = myf['integral'] 00412 myparams['integralerr'] = integralerr = myf['integralerr'] 00413 myparams['wantreturn'] = wantreturn = myf['wantreturn'] 00414 myparams['stretch'] = stretch = myf['stretch'] 00415 myparams['logresults'] = logresults = myf['logresults'] 00416 myparams['pampest'] = pampest = myf['pampest'] 00417 myparams['pcenterest'] = pcenterest = myf['pcenterest'] 00418 myparams['pfwhmest'] = pfwhmest = myf['pfwhmest'] 00419 myparams['pfix'] = pfix = myf['pfix'] 00420 myparams['gmncomps'] = gmncomps = myf['gmncomps'] 00421 myparams['gmampcon'] = gmampcon = myf['gmampcon'] 00422 myparams['gmcentercon'] = gmcentercon = myf['gmcentercon'] 00423 myparams['gmfwhmcon'] = gmfwhmcon = myf['gmfwhmcon'] 00424 myparams['gmampest'] = gmampest = myf['gmampest'] 00425 myparams['gmcenterest'] = gmcenterest = myf['gmcenterest'] 00426 myparams['gmfwhmest'] = gmfwhmest = myf['gmfwhmest'] 00427 myparams['gmfix'] = gmfix = myf['gmfix'] 00428 myparams['logfile'] = logfile = myf['logfile'] 00429 myparams['append'] = append = myf['append'] 00430 myparams['pfunc'] = pfunc = myf['pfunc'] 00431 myparams['goodamprange'] = goodamprange = myf['goodamprange'] 00432 myparams['goodcenterrange'] = goodcenterrange = myf['goodcenterrange'] 00433 myparams['goodfwhmrange'] = goodfwhmrange = myf['goodfwhmrange'] 00434 myparams['sigma'] = sigma = myf['sigma'] 00435 myparams['outsigma'] = outsigma = myf['outsigma'] 00436 00437 if type(gmampest)==float: gmampest=[gmampest] 00438 if type(gmcenterest)==float: gmcenterest=[gmcenterest] 00439 if type(gmfwhmest)==float: gmfwhmest=[gmfwhmest] 00440 if type(goodamprange)==float: goodamprange=[goodamprange] 00441 if type(goodcenterrange)==float: goodcenterrange=[goodcenterrange] 00442 if type(goodfwhmrange)==float: goodfwhmrange=[goodfwhmrange] 00443 00444 result = None 00445 00446 # 00447 # The following is work around to avoid a bug with current python translation 00448 # 00449 mytmp = {} 00450 00451 mytmp['imagename'] = imagename 00452 mytmp['box'] = box 00453 mytmp['region'] = region 00454 mytmp['chans'] = chans 00455 mytmp['stokes'] = stokes 00456 mytmp['axis'] = axis 00457 mytmp['mask'] = mask 00458 mytmp['ngauss'] = ngauss 00459 mytmp['poly'] = poly 00460 mytmp['estimates'] = estimates 00461 mytmp['minpts'] = minpts 00462 mytmp['multifit'] = multifit 00463 mytmp['model'] = model 00464 mytmp['residual'] = residual 00465 mytmp['amp'] = amp 00466 mytmp['amperr'] = amperr 00467 mytmp['center'] = center 00468 mytmp['centererr'] = centererr 00469 mytmp['fwhm'] = fwhm 00470 mytmp['fwhmerr'] = fwhmerr 00471 mytmp['integral'] = integral 00472 mytmp['integralerr'] = integralerr 00473 mytmp['wantreturn'] = wantreturn 00474 mytmp['stretch'] = stretch 00475 mytmp['logresults'] = logresults 00476 mytmp['pampest'] = pampest 00477 mytmp['pcenterest'] = pcenterest 00478 mytmp['pfwhmest'] = pfwhmest 00479 mytmp['pfix'] = pfix 00480 mytmp['gmncomps'] = gmncomps 00481 mytmp['gmampcon'] = gmampcon 00482 mytmp['gmcentercon'] = gmcentercon 00483 mytmp['gmfwhmcon'] = gmfwhmcon 00484 mytmp['gmampest'] = gmampest 00485 mytmp['gmcenterest'] = gmcenterest 00486 mytmp['gmfwhmest'] = gmfwhmest 00487 mytmp['gmfix'] = gmfix 00488 mytmp['logfile'] = logfile 00489 mytmp['append'] = append 00490 mytmp['pfunc'] = pfunc 00491 mytmp['goodamprange'] = goodamprange 00492 mytmp['goodcenterrange'] = goodcenterrange 00493 mytmp['goodfwhmrange'] = goodfwhmrange 00494 mytmp['sigma'] = sigma 00495 mytmp['outsigma'] = outsigma 00496 pathname='file:///'+os.environ.get('CASAPATH').split()[0]+'/share/xml/' 00497 trec = casac.utils().torecord(pathname+'specfit.xml') 00498 00499 casalog.origin('specfit') 00500 if not trec.has_key('specfit') or not casac.utils().verify(mytmp, trec['specfit']) : 00501 return False 00502 00503 00504 try : 00505 casalog.post('') 00506 casalog.post('##########################################') 00507 casalog.post('##### Begin Task: specfit #####') 00508 casalog.post('') 00509 result = specfit(imagename, box, region, chans, stokes, axis, mask, ngauss, poly, estimates, minpts, multifit, model, residual, amp, amperr, center, centererr, fwhm, fwhmerr, integral, integralerr, wantreturn, stretch, logresults, pampest, pcenterest, pfwhmest, pfix, gmncomps, gmampcon, gmcentercon, gmfwhmcon, gmampest, gmcenterest, gmfwhmest, gmfix, logfile, append, pfunc, goodamprange, goodcenterrange, goodfwhmrange, sigma, outsigma) 00510 casalog.post('') 00511 casalog.post('##### End Task: specfit #####') 00512 casalog.post('##########################################') 00513 00514 00515 # saveinputs for individule engine has no use 00516 # saveinputs should alos be removed from casa_in_py.py 00517 # 00518 # 00519 # saveinputs = myf['saveinputs'] 00520 # saveinputs('specfit', 'specfit.last', myparams) 00521 # 00522 # 00523 except Exception, instance: 00524 #print '**** Error **** ',instance 00525 pass 00526 00527 gc.collect() 00528 return result 00529 # 00530 # 00531 ## 00532 # def paramgui(self, useGlobals=True): 00533 # """ 00534 # Opens a parameter GUI for this task. If useGlobals is true, then any relevant global parameter settings are used. 00535 # """ 00536 # import paramgui 00537 # 00538 # a=inspect.stack() 00539 # stacklevel=0 00540 # for k in range(len(a)): 00541 # if (string.find(a[k][1], 'ipython console') > 0) or (string.find(a[k][1], '<string>') >= 0): 00542 # stacklevel=k 00543 # break 00544 # myf = sys._getframe(stacklevel).f_globals 00545 # 00546 # if useGlobals: 00547 # paramgui.setGlobals(myf) 00548 # else: 00549 # paramgui.setGlobals({}) 00550 # 00551 # paramgui.runTask('specfit', myf['_ip']) 00552 # paramgui.setGlobals({}) 00553 # 00554 # 00555 # 00556 # 00557 def defaults(self, param=None): 00558 a=inspect.stack() 00559 stacklevel=0 00560 for k in range(len(a)): 00561 if (string.find(a[k][1], 'ipython console') > 0) or (string.find(a[k][1], '<string>') >= 0): 00562 stacklevel=k 00563 break 00564 myf=sys._getframe(stacklevel).f_globals 00565 a = odict() 00566 a['imagename'] = '' 00567 a['box'] = '' 00568 a['region'] = '' 00569 a['chans'] = '' 00570 a['stokes'] = '' 00571 a['axis'] = -1 00572 a['mask'] = '' 00573 a['poly'] = -1 00574 a['estimates'] = '' 00575 a['minpts'] = 0 00576 a['multifit'] = False 00577 a['model'] = '' 00578 a['residual'] = '' 00579 a['wantreturn'] = True 00580 a['logresults'] = True 00581 a['gmncomps'] = 0 00582 a['gmampcon'] = '' 00583 a['gmcentercon'] = '' 00584 a['gmfwhmcon'] = '' 00585 a['gmampest'] = [0.0] 00586 a['gmcenterest'] = [0.0] 00587 a['gmfwhmest'] = [0.0] 00588 a['gmfix'] = '' 00589 a['logfile'] = '' 00590 a['goodamprange'] = [0.0] 00591 a['goodcenterrange'] = [0.0] 00592 a['goodfwhmrange'] = [0.0] 00593 a['sigma'] = '' 00594 00595 a['async']=False 00596 a['mask'] = { 00597 0:odict([{'notvalue':''}, {'stretch':False}])} 00598 a['estimates'] = { 00599 0:odict([{'value':''}, {'ngauss':1}, {'pampest':''}, {'pcenterest':''}, {'pfwhmest':''}, {'pfix':''}, {'pfunc':''}])} 00600 a['multifit'] = { 00601 0:odict([{'value':True}, {'amp':""}, {'amperr':""}, {'center':""}, {'centererr':""}, {'fwhm':""}, {'fwhmerr':""}, {'integral':""}, {'integralerr':""}])} 00602 a['gmncomps'] = { 00603 0:odict([{'notvalue':0}, {'gmampcon':''}, {'gmcentercon':''}, {'gmfwhmcon':''}, {'gmampest':0.0}, {'gmcenterest':0.0}, {'gmfwhmest':0.0}, {'gmfix':''}])} 00604 a['logfile'] = { 00605 0:odict([{'notvalue':''}, {'append':True}])} 00606 a['sigma'] = { 00607 0:odict([{'notvalue':''}, {'outsigma':''}])} 00608 00609 ### This function sets the default values but also will return the list of 00610 ### parameters or the default value of a given parameter 00611 if(param == None): 00612 myf['__set_default_parameters'](a) 00613 elif(param == 'paramkeys'): 00614 return a.keys() 00615 else: 00616 if(a.has_key(param)): 00617 #if(type(a[param]) == dict) : 00618 # return a[param][len(a[param])-1]['value'] 00619 #else : 00620 return a[param] 00621 00622 00623 # 00624 # 00625 def check_params(self, param=None, value=None): 00626 a=inspect.stack() 00627 stacklevel=0 00628 for k in range(len(a)): 00629 if (string.find(a[k][1], 'ipython console') > 0) or (string.find(a[k][1], '<string>') >= 0): 00630 stacklevel=k 00631 break 00632 myf=sys._getframe(stacklevel).f_globals 00633 00634 # print 'param:', param, 'value:', value 00635 try : 00636 if str(type(value)) != "<type 'instance'>" : 00637 value0 = value 00638 value = myf['cu'].expandparam(param, value) 00639 matchtype = False 00640 if(type(value) == numpy.ndarray): 00641 if(type(value) == type(value0)): 00642 myf[param] = value.tolist() 00643 else: 00644 #print 'value:', value, 'value0:', value0 00645 #print 'type(value):', type(value), 'type(value0):', type(value0) 00646 myf[param] = value0 00647 if type(value0) != list : 00648 matchtype = True 00649 else : 00650 myf[param] = value 00651 value = myf['cu'].verifyparam({param:value}) 00652 if matchtype: 00653 value = False 00654 except Exception, instance: 00655 #ignore the exception and just return it unchecked 00656 myf[param] = value 00657 return value 00658 00659 # 00660 # 00661 def description(self, key='specfit', subkey=None): 00662 desc={'specfit': 'Fit 1-dimensional gaussians and/or polynomial models to an image or image region', 00663 'imagename': 'Name of the input image', 00664 'box': 'Rectangular box in direction coordinate blc, trc. Default: entire image ("").', 00665 'region': 'Region of interest. See help par.region for possible specifications. Default: Do not use a region.', 00666 'chans': 'Channels to use. Channels must be contiguous. Default: all channels ("").', 00667 'stokes': 'Stokes planes to use. Planes must be contiguous. Default: all stokes ("").', 00668 'axis': 'The profile axis. Default: use the spectral axis if one exists, axis 0 otherwise (<0).', 00669 'mask': 'Mask to use. See help par.mask. Default is none..', 00670 'ngauss': 'Number of Gaussian elements. Default: 1.', 00671 'poly': 'Order of polynomial element. Default: do not fit a polynomial (<0).', 00672 'estimates': 'Name of file containing initial estimates. Default: No initial estimates ("").', 00673 'minpts': 'Minimum number of unmasked points necessary to attempt fit.', 00674 '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.', 00675 'model': 'Name of model image. Default: do not write the model image ("").', 00676 'residual': 'Name of residual image. Default: do not write the residual image ("").', 00677 'amp': 'Name of amplitude solution image. Default: do not write the image ("").', 00678 'amperr': 'Name of amplitude solution error image. Default: do not write the image ("").', 00679 'center': 'Name of center solution image. Default: do not write the image ("").', 00680 'centererr': 'Name of center solution error image. Default: do not write the image ("").', 00681 'fwhm': 'Name of fwhm solution image. Default: do not write the image ("").', 00682 'fwhmerr': 'Name of fwhm solution error image. Default: do not write the image ("").', 00683 'integral': 'Prefix of ame of integral solution image. Name of image will have gaussian component number appended. Default: do not write the image ("").', 00684 'integralerr': 'Prefix of name of integral error solution image. Name of image will have gaussian component number appended. Default: do not write the image ("").', 00685 'wantreturn': 'Should a record summarizing the results be returned?', 00686 'stretch': 'Stretch the mask if necessary and possible? See help par.stretch ', 00687 'logresults': 'Output results to logger?', 00688 'pampest': 'Initial estimate of PCF profile (gaussian or lorentzian) amplitudes.', 00689 'pcenterest': 'Initial estimate PCF profile centers, in pixels.', 00690 'pfwhmest': 'Initial estimate PCF profile FWHMs, in pixels.', 00691 'pfix': 'PCF profile parameters to fix during fit.', 00692 'gmncomps': 'Number of components in each gaussian multiplet to fit', 00693 'gmampcon': 'The amplitude ratio constraints for non-reference components to reference component in gaussian multiplets.', 00694 'gmcentercon': 'The center offset constraints (in pixels) for non-reference components to reference component in gaussian multiplets.', 00695 'gmfwhmcon': 'The FWHM ratio constraints for non-reference components to reference component in gaussian multiplets.', 00696 'gmampest': 'Initial estimate of individual gaussian amplitudes in gaussian multiplets.', 00697 'gmcenterest': 'Initial estimate of individual gaussian centers in gaussian multiplets, in pixels.', 00698 'gmfwhmest': 'Initial estimate of individual gaussian FWHMss in gaussian multiplets, in pixels.', 00699 'gmfix': 'Parameters of individual gaussians in gaussian multiplets to fix during fit.', 00700 'logfile': 'File in which to log results. Default is not to write a logfile.', 00701 'append': 'Append results to logfile? Logfile must be specified. Default is to append. False means overwrite existing file if it exists.', 00702 'pfunc': 'PCF singlet functions to fit. "gaussian" or "lorentzian" (minimal match supported). Unspecified means all gaussians.', 00703 'goodamprange': 'Acceptable amplitude solution range. [0.0] => all amplitude solutions are acceptable.', 00704 'goodcenterrange': 'Acceptable center solution range in pixels relative to region start. [0.0] => all center solutions are acceptable.', 00705 'goodfwhmrange': 'Acceptable FWHM solution range in pixels. [0.0] => all FWHM solutions are acceptable.', 00706 'sigma': 'Standard deviation array or image name.', 00707 'outsigma': 'Name of output image used for standard deviation. Ignored if sigma is empty.', 00708 00709 'async': 'If true the taskname must be started using specfit(...)' 00710 } 00711 00712 # 00713 # Set subfields defaults if needed 00714 # 00715 00716 if(desc.has_key(key)) : 00717 return desc[key] 00718 00719 def itsdefault(self, paramname) : 00720 a = {} 00721 a['imagename'] = '' 00722 a['box'] = '' 00723 a['region'] = '' 00724 a['chans'] = '' 00725 a['stokes'] = '' 00726 a['axis'] = -1 00727 a['mask'] = '' 00728 a['ngauss'] = 1 00729 a['poly'] = -1 00730 a['estimates'] = '' 00731 a['minpts'] = 0 00732 a['multifit'] = False 00733 a['model'] = '' 00734 a['residual'] = '' 00735 a['amp'] = '' 00736 a['amperr'] = '' 00737 a['center'] = '' 00738 a['centererr'] = '' 00739 a['fwhm'] = '' 00740 a['fwhmerr'] = '' 00741 a['integral'] = '' 00742 a['integralerr'] = '' 00743 a['wantreturn'] = True 00744 a['stretch'] = False 00745 a['logresults'] = True 00746 a['pampest'] = '' 00747 a['pcenterest'] = '' 00748 a['pfwhmest'] = '' 00749 a['pfix'] = '' 00750 a['gmncomps'] = 0 00751 a['gmampcon'] = '' 00752 a['gmcentercon'] = '' 00753 a['gmfwhmcon'] = '' 00754 a['gmampest'] = [0.0] 00755 a['gmcenterest'] = [0.0] 00756 a['gmfwhmest'] = [0.0] 00757 a['gmfix'] = '' 00758 a['logfile'] = '' 00759 a['append'] = True 00760 a['pfunc'] = '' 00761 a['goodamprange'] = [0.0] 00762 a['goodcenterrange'] = [0.0] 00763 a['goodfwhmrange'] = [0.0] 00764 a['sigma'] = '' 00765 a['outsigma'] = '' 00766 00767 if a.has_key(paramname) : 00768 return a[paramname] 00769 specfit_pg = specfit_pg_()