1
2
3 """
4 Interface for statistics.
5
6 Author: R. Lombaert
7
8 """
9
10 import os
11 import scipy.stats
12 import numpy as np
13
14 import cc.path
15 from cc.tools.io import DataIO
16 from cc.data.instruments import Pacs
17 from cc.data import Data
18 from cc.modeling.objects import Star
19
20
21
23
24 """
25 An evironment for various statistical calculations for transitions.
26
27 """
28
29 - def __init__(self,star_name,path_code='codeSep2010',code='GASTRoNOoM'):
30
31 """
32 Initializing an instance of Statistics.
33
34 @param star_name: Star name from Star.dat
35 @type star_name: string
36
37 @keyword path_code: Output folder in the code's home folder
38
39 (default: 'codeSep2010')
40 @type path_code: string
41 @keyword code: the code used for producing your output
42
43 (default: 'GASTRoNOoM')
44 @type code: string
45
46 """
47
48 self.star_name = star_name
49 self.star_grid = []
50 self.instrument = None
51 self.path_code = path_code
52 self.code = code
53 self.data_stats = dict()
54 self.data_info = dict()
55
56
57
58 - def setInstrument(self,instrument_name,searchstring='',\
59 data_filenames=[],instrument_instance=None,\
60 sample_transitions=[],\
61 stat_method='clipping',**kwargs):
62
63 '''
64 Set an instrument as a source of data. Define the instrument_name on
65 calling the method: PACS, SPIRE, SED, FREQ_RESO.
66
67 FREQ_RESO requires no additional keywords nor an instrument_instance.
68
69 For the other options, the instrument instance that contains the data
70 references and other info can be passed as an object, or can be set in
71 situ. If the latter, include required keywords for creating the
72 instrument object.
73
74 For PACS: resolution
75 For SPIRE: resolution, oversampling
76
77 See the respective modules __init__ methods for more information on
78 these keywords as well as other optional keys.
79
80 Spire and Pacs objects have to set data manually. Provide the
81 data_filenames or the searchstring for this. Not needed for FREQ_RESO or
82 SED.
83
84 @param instrument_name: The instrument (such as 'PACS', 'SPIRE', 'SED',
85 or 'FREQ_RESO')
86 @type instrument_name: string
87
88 @keyword data_filenames: The data filenames. If empty, auto-search is
89 done, with searchstring. Only for Pacs and
90 Spire
91
92 (default: [])
93 @type data_filenames: list[string]
94 @keyword searchstring: the searchstring conditional for the auto-search,
95 only for Pacs and Spire.
96
97 (default: '')
98 @type searchstring: string
99 @keyword instrument_instance: If None a new instance is made, otherwise
100 the instance is given here for the
101 instrument. Remember to pass the relevant
102 keywords to the method for the creation.
103
104 (default: None)
105 @type instrument_instance: Instrument()
106 @keyword stat_method: Std/Mean/RMS determination method for unresolved
107 data (Pacs and Spire) Options:
108 - 'clipping': 1-sigma clipping based on std/...
109 of full spectrum, then takes
110 std/... of clipped spectrum.
111 - 'preset': wavelength ranges for std/...
112 determination taken from Data.dat.
113
114 (default: 'clipping')
115 @type stat_method: string
116
117 '''
118 instrument_name = instrument_name.upper()
119 if instrument_name == 'PACS':
120 if not instrument_instance is None:
121 self.instrument = instrument_instance
122 self.instrument.setData(data_filenames=data_filenames,\
123 searchstring=searchstring)
124 else:
125 self.instrument = Pacs.Pacs(star_name=self.star_name,\
126 path=self.path_code,**kwargs)
127 self.instrument.setData(data_filenames=data_filenames,\
128 searchstring=searchstring)
129
130 elif instrument_name == 'SPIRE':
131 if not instrument_instance is None:
132 self.instrument = instrument_instance
133 self.instrument.setData(data_filenames=data_filenames,\
134 searchstring=searchstring)
135 else:
136 self.instrument = Spire.Spire(star_name=self.star_name,\
137 path=self.path_code,**kwargs)
138 self.instrument.setData(data_filenames=data_filenames,\
139 searchstring=searchstring)
140
141 elif instrument_name == 'SED':
142 if not instrument_instance is None:
143 self.instrument = instrument_instance
144 else:
145 self.instrument = Sed.Sed(star_name=self.star_name,**kwargs)
146
147 elif instrument_name == 'FREQ_RESO':
148 self.instrument = None
149
150
151
152 self.doDataStats(method=stat_method)
153
154
155
157
158 '''
159 Load the models and remember them.
160
161 @keyword star_grid: The parameter sets, if not given: model ids needed
162
163 (default: [])
164 @type star_grid: list[Star()]
165 @keyword models: the PACS ids, only relevant if star_grid == [] and if
166 instrument == PACS. In all other cases a star_grid is
167 required.
168
169 (default: [])
170 @type models: list[string]
171 @keyword extra_keywords: any instrument specific keywords that you need
172
173 (default: dict())
174 @type extra_keywords: dict
175
176 '''
177
178
179 if self.instrument and self.instrument.instrument == 'SED':
180 if not star_grid:
181 raise IOError('Statistics.setModels requires a ' + \
182 'star_grid to be defined for SED data.')
183 if set([s['LAST_MCMAX_MODEL'] for s in star_grid]) == set(['']):
184 return
185
186 self.star_grid = [s for s in star_grid if s['LAST_MCMAX_MODEL']]
187 self.star_grid = np.array(self.star_grid)
188
189
190 elif self.instrument:
191 instr = self.instrument.instrument.upper()
192 print '***********************************'
193 print '* Checking Sphinx models for comparison with %s.'%instr
194 if instr == 'SPIRE': models = []
195 if not star_grid and not models:
196 raise IOError('Statistics.setModels requires either ' + \
197 'star_grid or models to be defined.')
198
199 elif not star_grid:
200 star_grid = Star.makeStars(models=models,id_type='PACS',\
201 path=path_code,code=self.code)
202 self.instrument.addStarPars(star_grid)
203 if set([s['MOLECULE'] and 1 or 0 for s in star_grid]) == set([0]):
204 return
205 self.instrument.prepareSphinx(star_grid)
206 self.star_grid = [star
207 for star in star_grid
208 if not star['LAST_%s_MODEL'%instr] is None]
209
210 else:
211 if not star_grid:
212 raise IOError('Statistics.setModels requires a ' + \
213 'star_grid to be defined for freq-resolved data.')
214 if set([s['LAST_GASTRONOOM_MODEL'] for s in star_grid]) == set(['']):
215 return
216 if set([s['MOLECULE'] and 1 or 0 for s in star_grid]) == set([0]):
217 return
218 self.star_grid = star_grid
219
220
221
223
224 '''
225 Calculate mean, std, rms for the data files.
226
227 They are saved in self.data_stats.
228
229 @keyword method: Std/Mean/RMS determination method. Can be 'clipping'
230 or 'preset'. The former 1-sigma clips based on
231 mean/std/... of full spectrum, then takes mean/std/...
232 of clipped spectrum. Latter takes preset values from
233 Data.dat.
234
235 (default: 'clipping')
236 @type method: string
237
238 '''
239
240 method = method.lower()
241 if method not in ['preset','clipping']: method = 'clipping'
242 self.data_stats = dict()
243 if self.instrument \
244 and self.instrument.instrument.upper() in ['PACS','SPIRE']:
245 if not self.data_info: self.setDataInfo()
246 for filename,data_wave,data_flux,band in \
247 zip(self.instrument.data_filenames,\
248 self.instrument.data_wave_list,\
249 self.instrument.data_flux_list,\
250 self.instrument.data_ordernames):
251 these_stats = dict()
252 iband = self.data_info['bands'].index(band)
253 these_stats['sigma'] = self.data_info['sigmas'][iband]
254 if method == 'preset':
255 w_std_min = self.data_info['w_std_min'][iband]
256 w_std_max = self.data_info['w_std_max'][iband]
257 if w_std_min < data_wave[0] or w_std_max > data_wave[-1]:
258 method = 'clipping'
259 test_mean = Data.getMean(wave=data_wave,flux=data_flux,\
260 wmin=w_std_min,wmax=w_std_max)
261 if test_mean is None:
262 method = 'clipping'
263 if method == 'clipping':
264 totmean = Data.getMean(wave=data_wave,flux=data_flux)
265 totstd = Data.getStd(wave=data_wave,flux=data_flux)
266 limits = (-totstd,totstd)
267 these_stats['mean'] = Data.getMean(flux=data_flux,\
268 limits=limits)
269 these_stats['std'] = Data.getStd(flux=data_flux,\
270 limits=limits)
271 these_stats['rms'] = Data.getRMS(flux=data_flux-\
272 these_stats['mean'],\
273 limits=limits)
274 else:
275 these_stats['mean'] = Data.getMean(wave=data_wave,\
276 flux=data_flux,\
277 wmin=w_std_min,\
278 wmax=w_std_max)
279 these_stats['std'] = Data.getStd(wave=data_wave,\
280 flux=data_flux,\
281 wmin=w_std_min,\
282 wmax=w_std_max)
283 these_stats['rms'] = Data.getRMS(wave=data_wave,\
284 flux=data_flux-\
285 these_stats['mean'],\
286 wmin=w_std_min,\
287 wmax=w_std_max)
288 print '* Data statistics for %s using "%s" method:'\
289 %(filename,method)
290 if method == 'preset':
291 print '* Taken between %.2f and %.2f micron.'\
292 %(w_std_min,w_std_max)
293 print 'mean = ',these_stats['mean'],' Jy'
294 print 'std = ',these_stats['std'],' Jy'
295 print 'RMS = ',these_stats['rms'],' Jy'
296 self.data_stats[filename] = these_stats
297 print '***********************************'
298
299
300
302
303 '''
304 Read data info from the ComboCode/usr/Data.dat file. Will return
305 information such as sigma levels and wavelength range for determining
306 the std value.
307
308 '''
309
310 if not self.instrument: return
311 self.data_info = dict()
312 instbands = DataIO.getInputData(keyword='INSTRUMENT',\
313 filename='Data.dat')
314 instbands = [band.upper() for band in instbands]
315 indices = [i
316 for i,band in enumerate(instbands)
317 if band.find(self.instrument.instrument.upper()) == 0]
318 bands = [band.replace('%s_'%self.instrument.instrument.upper(),'')
319 for i,band in enumerate(instbands)
320 if i in indices]
321 w_std_min = [wmin
322 for i,wmin in enumerate(DataIO.getInputData(\
323 keyword='W_STD_MIN',filename='Data.dat'))
324 if i in indices]
325 w_std_max = [wmax
326 for i,wmax in enumerate(DataIO.getInputData(\
327 keyword='W_STD_MAX',filename='Data.dat'))
328 if i in indices]
329 sigmas = [sigma
330 for i,sigma in enumerate(DataIO.getInputData(\
331 keyword='SIGMA',filename='Data.dat'))
332 if i in indices]
333 self.data_info['sigmas'] = sigmas
334 self.data_info['bands'] = bands
335 self.data_info['w_std_min'] = w_std_min
336 self.data_info['w_std_max'] = w_std_max
337