1
2
3 """
4 Performing statistics on SED data. Currently includes photometry only.
5
6 Author: R. Lombaert
7
8 """
9
10 import os
11 from scipy.interpolate import interp1d
12 import operator
13 import types
14 import numpy as np
15
16 import cc.path
17 from cc.tools.io import DataIO
18 from cc.data import Sed
19 from cc.modeling.codes import MCMax
20 from cc.modeling.tools import Reddening
21 from cc.statistics import BasicStats
22 from cc.statistics.Statistics import Statistics
23
24
25
27
28 """
29 Environment with several tools to perform statistics on SED photometry.
30
31 Then run several methods to set data, and models. e.g. (requires a
32 star_grid to be defined):
33
34 >>> sed = Sed.Sed('rscl')
35 >>> sedstats = SedStats.SedStats(star_name='rscl',path_code='codeSep2015')
36 >>> sedstats.setInstrument(sed=sed)
37 >>> sedstats.setModels(star_grid=star_grid)
38 >>> sedstats.setModelPhotometry()
39
40 Alternatively, you can take the SedStats object from a ComboCode object
41 given that STATISTICS=1 in the inputfile for CC. Then the above is
42 already done for you. Assuming the star_name is 'rscl':
43
44 >>> model = ComboCode.ComboCode('inputfile.dat')
45 >>> model.startSession()
46 >>> sedstats = model.sedstats['rscl']
47
48 Now you can go ahead and calculate chi^2 values and write them to a
49 file:
50
51 >>> chi2 = sedstats.calcChi2()
52 >>> sedstats.writeChi2('myfile.dat')
53
54 The writeChi2 method writes the latest chi2 calculation to a file with
55 given filename. Several options are available for selecting and sorting
56 photometry during the calculation, and sorting and adding extra columns
57 of model parameters when writing the results to a file. Check the
58 docstrings of each method for more information.
59
60 """
61
62 - def __init__(self,star_name,code='MCMax',path_code='codeSep2015'):
63
64 """
65 Initializing an instance of SedStats.
66
67 @param star_name: Star name from Star.dat
68 @type star_name: string
69
70 @keyword code: the code used for producing your output
71
72 (default: 'MCMax')
73 @type code: string
74 @keyword path_code: Output folder in the code's home folder
75
76 (default: 'codeSep2015')
77 @type path_code: string
78
79 """
80
81 super(SedStats,self).__init__(star_name=star_name,\
82 code=code,path_code=path_code)
83
84
85
86
87 self.mphot_ivs = []
88 self.mphot_other = dict()
89
90
91
92 self.dphot_ivs = dict([('wave',np.empty(0)),\
93 ('phot',np.empty(0)),\
94 ('ephot',np.empty(0))])
95 self.dphot_other = dict()
96 self.photbands = np.empty(0)
97 self.sed = None
98
99 self.chi2 = np.empty(0)
100
101
102
104
105 '''
106 Set and read the data objects for this statistics module.
107
108 The Sed object is made on the spot if sed is None. Extra arguments for
109 making the Sed can be passed through the method in kwargs.
110
111 @keyword sed: The Sed object. In case of default, it is made on the spot
112 for given star with extra arguments kwargs if needed
113
114 (default: None)
115 @type sed: Sed()
116
117 '''
118
119 super(SedStats,self).setInstrument(instrument_name='SED',\
120 instrument_instance=sed,**kwargs)
121 self.sed = self.instrument
122 self.photbands = sed.photbands
123
124 for (dt,fn),tdata in sorted([dset
125 for dset in self.sed.data.items()
126 if 'PHOT' in dset[0][0].upper()]):
127 if dt.lower() == 'photometric_ivs':
128 self.dphot_ivs['wave'] = tdata[0]
129 self.dphot_ivs['phot'] = tdata[1]
130 self.dphot_ivs['ephot'] = tdata[2]
131 else:
132 self.dphot_other[fn] = dict()
133 self.dphot_other[fn]['wave'] = tdata[0]
134 self.dphot_other[fn]['phot'] = tdata[1]
135 self.dphot_other[fn]['ephot'] = tdata[2]
136
137
139
140 '''
141 Prepare the model photometry to be compared with the data.
142
143 Two kinds: The IvS photometry with proper photometric bands, and other
144 photometry to be compared with the interpolated model spectrum.
145
146 '''
147
148
149
150 mids = [s['LAST_MCMAX_MODEL'] for s in self.star_grid]
151 if not mids:
152 print "No successfully calculated MCMax models found."
153 return
154
155 self.mwave = []
156 self.mflux = []
157 for model_id,s in zip(mids,self.star_grid):
158 dpath = os.path.join(cc.path.mout,'models',model_id)
159 fn_spec = 'spectrum{:04.1f}.dat'.format(s['RT_INCLINATION'])
160 w,f = MCMax.readModelSpectrum(dpath,s['RT_SPEC'],fn_spec)
161 if s['REDDENING']:
162 print 'Reddening models to correct for interstellar extinction.'
163 ak = self.sed.getAk(s['DISTANCE'],s['REDDENING_MAP'],\
164 s['REDDENING_LAW'])
165 f = Reddening.redden(w,f,ak,law=s['REDDENING_LAW'])
166 self.mwave.append(w)
167 self.mflux.append(f)
168
169 if self.photbands.size:
170 self.mphot_ivs = [Sed.calcPhotometry(w,f,self.photbands)
171 for w,f in zip(self.mwave,self.mflux)]
172
173 for fn in self.dphot_other.keys():
174 self.mphot_other[fn] = []
175 if self.dphot_other.keys():
176 for w,f in zip(self.mwave,self.mflux):
177 interp = interp1d(w,f)
178 for fn in self.dphot_other.keys():
179 finter = interp(self.dphot_other[fn]['wave'])
180 self.mphot_other[fn].append(finter)
181
182
183
184 - def calcChi2(self,ndf=0,fns=None,cwave=0.0,phot_ivs=1,sort=1,\
185 chi2_method='diff'):
186
187 '''
188 Calculate, save and return the chi^2 values for a given set of models.
189
190 For now, only photometry (both with and without photometric bands
191 available, the latter being associated with the Photometric_IvS file)
192 is taken into account.
193
194 You can specify the number of degrees of freedom for the reduced chi^2
195 in case you are comparing models probing a grid across ndf different
196 parameters.
197
198 The method overrides previously calculated chi2 in the object. Write the
199 chi^2 to a file before running the method if you want to save the values
200 with self.writeChi2()
201
202
203 @keyword ndf: Number of degrees of freedom. Default in case of
204 calculating for one single model. Typically the number of
205 variable grid parameters in a grid calculation.
206
207 (default: 0)
208 @type ndf: int
209 @keyword fns: The photometric filenames to be included in the chi^2 calc
210 Default if all are to be included. Empty list if none are
211 to be included.
212
213 (default: None)
214 @type fns: list(str)
215 @keyword cwave: A cut-off wavelength in micron used to calculate the
216 chi^2 stat for all photometry above or equal to the
217 wavelength in micron. Use a negative value to grab all
218 wavelengths lower than cwave. Default if no cut-off.
219
220 (default: 0.0)
221 @type cwave: float
222 @keyword phot_ivs: Include the Photometric_IvS photometry as well.
223
224 (default: 1)
225 @type phot_ivs: bool
226 @keyword sort: Sort the chi^2 values from lowest to highest.
227
228 (default: 1)
229 @type sort: bool
230 @keyword chi2_method: Method for calculating chi^2. Can be diff or
231 log (see BasicStats for more information)
232
233 (default: 'diff')
234 @type chi2_method: str
235
236 @return: The list of chi^2 values with the same length as the model grid
237 of valid MCMax models.
238 @rtype: list
239
240 '''
241
242
243
244 mphot = [[] for s in self.star_grid]
245 dphot = []
246 ephot = []
247
248 if fns is None:
249 fns = self.dphot_other.keys()
250 else:
251 fns = [fn if os.path.split(fn)[0] else os.path.join(cc.path.dsed,fn)
252 for fn in fns]
253 fns = [fn for fn in fns if fn in self.dphot_other.keys()]
254
255
256
257 if not fns and not phot_ivs:
258 fns = self.dphot_other.keys()
259
260 for fn in fns:
261
262 if cwave >= 0.0:
263 keep = self.dphot_other[fn]['wave'] >= cwave
264 else:
265 keep = self.dphot_other[fn]['wave'] < -cwave
266
267
268 if not keep.any():
269 continue
270
271
272 for i,iphot in enumerate(self.mphot_other[fn]):
273 mphot[i].extend(iphot[keep])
274 dphot.extend(self.dphot_other[fn]['phot'][keep])
275 ephot.extend(self.dphot_other[fn]['ephot'][keep])
276
277 if phot_ivs and self.photbands.size:
278 if cwave >= 0.0:
279 keep = self.dphot_ivs['wave'] >= cwave
280 else:
281 keep = self.dphot_ivs['wave'] < -cwave
282
283
284 if keep.any():
285 for i,iphot in enumerate(self.mphot_ivs):
286 mphot[i].extend(iphot[keep])
287 dphot.extend(self.dphot_ivs['phot'][keep])
288 ephot.extend(self.dphot_ivs['ephot'][keep])
289
290
291 if not dphot:
292 self.chi2 = np.empty(0)
293 return self.chi2
294
295 self.chi2 = []
296 for iphot in mphot:
297 self.chi2.append(BasicStats.calcChiSquared(dphot,iphot,ephot,ndf,\
298 chi2_method))
299 self.chi2 = np.array(self.chi2)
300
301 return np.sort(self.chi2) if sort else self.chi2
302
303
304
306
307 '''
308 Return the star grid for which the chi^2 is calculated.
309
310 This grid may differ from the original ComboCode grid because the local
311 star_grid is filtered for unavailable MCMax models (either not
312 calculated or not loaded).
313
314 Each object in the star_grid is itself a dictionary that contains all
315 the parameters of the model.
316
317 @keyword sort: Sort the star_grid according to the chi^2 values from
318 lowest to highest. Requires calcChi2 to be ran first.
319
320 (default: 1)
321 @type sort: bool
322
323 @return: The star grid
324 @rtype: list(Star())
325
326 '''
327
328 if not self.chi2.size:
329 sort = 0
330 if sort:
331 isort = np.argsort(self.chi2)
332
333 return self.star_grid[isort] if sort else self.star_grid
334
335
336
337 - def writeChi2(self,fn,sort=1,parameters=[]):
338
339 '''
340 Write the Chi^2 values to a file. Lists the model id in the first column
341 with the chi^2 value in the second.
342
343 The chi^2 values can be requested to be sorted.
344
345 Parameters from the Star() objects can be added as additional columns.
346 Given parameters must be valid.
347
348 @param fn: The output filename
349 @type fn: str
350
351 @keyword sort: Sort the star_grid according to the chi^2 values from
352 lowest to highest. Requires calcChi2 to be ran first.
353
354 (default: 1)
355 @type sort: bool
356 @keyword parameters: The additional model parameters to be added as
357 columns in the file.
358
359 (default: [])
360 @type parameters: list(str)
361
362 '''
363
364
365 if not self.chi2.size:
366 return
367
368
369 comments = ['# '] + ['ID','RedChi^2'] + parameters + ['\n']
370 DataIO.writeFile(filename=fn,input_lines=comments,delimiter='\t')
371
372
373 cols = [[s['LAST_MCMAX_MODEL'] for s in self.getStarGrid(sort=sort)]]
374 if sort:
375 isort = np.argsort(self.chi2)
376 cols.append(self.chi2[isort])
377 else:
378 cols.append(self.chi2)
379
380
381 for par in parameters:
382 cols.append([s[par] for s in self.getStarGrid(sort=sort)])
383
384
385 DataIO.writeCols(filename=fn,cols=cols,mode='a')
386