ComboCode.cc.ComboCode

31 32 ''' 33 The interface with which to run the ComboCode package. 34 35 ''' 36

37 - def __init__(self,inputfilename):

38 39 ''' 40 Initializing a ComboCode instance. 41 42 Once this is done, you only need to run startSession(). Then all 43 methods in this class will be called according to your inputfile. Only 44 run separate methods of the class if you know what you are doing! 45 46 Input is read and parsed, and the parameter objects (Star()) are set 47 48 The instrument and data objects, and the plotting manager are set. 49 50 The .spec file is updated here, if requested. 51 52 The inputfile can be given on the command line as: 53 python ComboCode.py inputComboCode.dat 54 55 In the python or ipython shell you can do: 56 >>> import ComboCode 57 >>> cc = ComboCode.ComboCode('/home/robinl/ComboCode/input/inputComboCode.dat') 58 59 @param inputfilename: The name of the inputfile. 60 @type inputfilename: string 61 62 ''' 63 64 self.inputfilename = inputfilename 65 self.readInput() 66 self.setGlobalPars() 67 self.setOutputFolders() 68 self.setPacs() 69 self.setSpire() 70 self.setSed() 71 self.setRadio() 72 self.setPlotManager() 73 self.setVarPars() 74 self.createStarGrid() 75 #- Only the extra transition pars will differ across the grid, so grab 76 #- the transition list from one of the Star() objects 77 if self.update_spec: 78 Transition.updateLineSpec(self.star_grid[0]['GAS_LINES']) 79 self.finished = False

80 81

82 - def startSession(self):

83 84 ''' 85 Start a ComboCode session, based on the input read upon initialisation. 86 87 The supercomputer and model managers are set and ran. 88 89 The plot manager, statistics module, fitter modules are ran if 90 requested. 91 92 The session ends by printing some info about the Star() objects. 93 94 Once started, the ComboCode object cannot be started again. You will 95 have to re-initialize. This will change in the future. 96 97 ''' 98 99 if not self.finished: 100 self.setVicManager() 101 self.setModelManager() 102 self.finished = True 103 self.runModelManager() 104 self.finalizeVic() 105 self.runChemistry() 106 self.runPlotManager() 107 self.runStatistics() 108 self.doContDiv() 109 #self.appendResults() 110 if self.write_dust_density: 111 [star.writeDensity() for star in self.star_grid] 112 self.printStarInfo() 113 else: 114 print 'This CC session is already finished. Please, create a new one.'

115 116 117

118 - def setGlobalPars(self):

119 120 ''' 121 Set the global parameters for this CC session. 122 123 ''' 124 125 default_global = [('mcmax',1),('gastronoom',1),('sphinx',1),('vic',0),\ 126 ('iterations',2),('plot_iterative',0),\ 127 ('vic_account','vsc30226'),('statistics',0),\ 128 ('vic_time_per_sphinx',30),('vic_credits',None),\ 129 ('append_results',0),('write_dust_density',0),\ 130 ('replace_db_entry',0),('update_spec',0),\ 131 ('path_gastronoom',''),('path_mcmax',''),\ 132 ('path_chemistry',''),\ 133 ('print_model_info',1),('stat_chi2','diff'),\ 134 ('contdiv_features',[]),('cfg_contdiv',''),\ 135 ('show_contdiv',0),('skip_cooling',0),\ 136 ('recover_sphinxfiles',0),('stat_print',0),\ 137 ('stat_lll_p',None),('stat_method','clipping'),\ 138 ('star_name','model'),('single_session',0),\ 139 ('stat_lll_vmin',0.0),('chemistry',0),\ 140 ('stat_lll_vmax',0.0), ('print_check_t',1),\ 141 ('chemstats',0),('chemstats_molecules',[])] 142 global_pars = dict([(k,self.processed_input.pop(k.upper(),v)) 143 for k,v in default_global]) 144 self.__dict__.update(global_pars) 145 self.__setStarName() 146 self.vic = 0 147 if not self.gastronoom or not self.mcmax: self.iterations = 1 148 if (not self.path_mcmax and self.mcmax): 149 raise IOError('Please define PATH_MCMAX in your inputfile.') 150 if (not self.path_gastronoom and self.gastronoom): 151 raise IOError('Please define PATH_GASTRONOOM in your inputfile.')

152 153 154

155 - def __setStarName(self):

156 157 ''' 158 Set star_name for the ComboCode object as a tuple. 159 160 Typically this is only one name for a standard modelling session, but 161 can be made multiple names as well for a statistical study. 162 163 The ComboCode object keeps track of all the data in dicts. 164 165 ''' 166 167 if self.multiplicative_grid.has_key('STAR_NAME') \ 168 or self.additive_grid.has_key('STAR_NAME'): 169 raise IOError('STAR_NAME incorrectly defined. Use & for grids.') 170 171 if isinstance(self.star_name,str): 172 self.star_name = (self.star_name,)

173 174 175

176 - def setPacs(self):

177 178 ''' 179 Collect the PACS relevant parameters from the inputfile and set the 180 PACS object. 181 182 ''' 183 184 pacs = self.processed_input.pop('PACS',0) 185 redo_convolution = self.processed_input.pop('PACS_REDO_CONVOLUTION',0) 186 searchstring = self.processed_input.pop('PACS_SEARCHSTRING','') 187 oversampling = self.processed_input.pop('PACS_OVERSAMPLING','') 188 intrinsic = self.processed_input.pop('PACS_INTRINSIC',1) 189 linefit = self.processed_input.pop('PACS_LINEFIT','') 190 191 #-- If PACS is not requested, put self.pacs to None. Still popping the 192 # PACS specific keywords to avoid clutter in the Star() objects. 193 # In case of a pure model, no data are available anyway. 194 self.pacs = dict() 195 for sn in self.star_name: 196 if not pacs or sn == 'model': 197 self.pacs[sn] = None 198 continue 199 self.pacs[sn] = Pacs.Pacs(star_name=sn,\ 200 path=self.path_gastronoom,\ 201 redo_convolution=redo_convolution,\ 202 oversampling=oversampling,\ 203 intrinsic=intrinsic,\ 204 path_linefit=linefit) 205 self.pacs[sn].setData(searchstring=searchstring)

206 207 208

209 - def setSpire(self):

210 211 ''' 212 Collect the SPIRE relevant parameters from the inputfile and set the SPIRE 213 object. 214 215 ''' 216 217 spire = self.processed_input.pop('SPIRE',0) 218 searchstring = self.processed_input.pop('SPIRE_SEARCHSTRING','') 219 resolution = self.processed_input.pop('SPIRE_RESOLUTION',0) 220 intrinsic = self.processed_input.pop('SPIRE_INTRINSIC',1) 221 oversampling = self.processed_input.pop('SPIRE_OVERSAMPLING',0) 222 linefit = self.processed_input.pop('SPIRE_LINEFIT','') 223 224 #-- If SPIRE is not requested, put self.spire to None. Still popping the 225 # SPIRE specific keywords to avoid clutter in the Star() objects. 226 # In case of a pure model, no data are available anyway. 227 self.spire = dict() 228 for sn in self.star_name: 229 if not spire or sn == 'model': 230 self.spire[sn] = None 231 continue 232 self.spire[sn] = Spire.Spire(star_name=sn,\ 233 path=self.path_gastronoom,\ 234 resolution=resolution,\ 235 intrinsic=intrinsic,\ 236 oversampling=oversampling,\ 237 path_linefit=linefit) 238 self.spire[sn].setData(searchstring=searchstring)

239 240 241

242 - def setSed(self):

243 244 ''' 245 Collect the SED data and create an Sed() object. 246 247 ''' 248 249 sed = self.processed_input.pop('SED',0) 250 remove = self.processed_input.pop('SED_PHOT_REMOVE','') 251 252 if not remove: remove = [] 253 elif isinstance(remove,str): remove = [remove] 254 else: remove = list(remove) 255 256 #-- If SED is not requested, put self.spire to None. Still popping the 257 # SED specific keywords to avoid clutter in the Star() objects. 258 # In case of a pure model, no data are available anyway. 259 self.sed = dict() 260 for sn in self.star_name: 261 if not sed or sn == 'model': 262 self.sed[sn] = None 263 continue 264 self.sed[sn] = Sed.Sed(star_name=sn,remove=remove)

265 266 267

268 - def setRadio(self):

269 270 ''' 271 Collect the relevant radio data for the requested star. Only done if 272 the pathname to the data is given. 273 274 If a database is not present, it is created. 275 276 If the radio_autosearch flag is on, transitions are automatically 277 generated based on the available data. Note that in this case, N_QUAD 278 from Star() is taken. 279 280 ''' 281 282 radio = self.processed_input.pop('RADIO',0) 283 radio_autosearch = self.processed_input.pop('RADIO_AUTOSEARCH',0) 284 285 if not radio: 286 self.radio = {sn: None for sn in self.star_name} 287 return 288 289 #-- If RADIO is not requested, put self.radio to None. Still popping the 290 # RADIO specific keywords to avoid clutter in the Star() objects. 291 # In case of a pure model, the radio db doesn't have an entry anyway. 292 radio_db = Radio.Radio() 293 294 self.radio = {sn: radio_db.get(sn, None) for sn in self.star_name} 295 self.radio_trans = {sn: None for sn in self.star_name} 296 297 if not radio_autosearch: 298 return 299 300 #-- Get the transition definitions (are in the correct format 301 # automatically, due to the methods in Radio.py). Check entry info 302 # is still ran, eg to get rid of bad 0 offset type sets. 303 #-- 11 entries, n_quad is added in the Star() object (to allow for 304 # proper n_quad gridding). N_quad specification through manual 305 # TRANSITION definition is still possible, but then overrides the 306 # default N_QUAD value. 307 radio_trans = sorted(set([tr.replace('TRANSITION=','',1) 308 for sn in self.star_name 309 if self.radio[sn] 310 for tr in self.radio[sn].keys()])) 311 radio_trans = DataIO.checkEntryInfo(radio_trans,11,'TRANSITION') 312 313 #-- Doubles of transitions not possible from db. radio_trans always 314 # a tuple. 315 #-- In case doubles after merge: they will be filtered out by Star() 316 if self.additive_grid.has_key('TRANSITION'): 317 otrl = self.additive_grid['TRANSITION'] 318 ntrl = [tl + radio_trans for tl in otrl] 319 self.additive_grid['TRANSITION'] = ntrl 320 elif self.processed_input.has_key('TRANSITION'): 321 self.processed_input['TRANSITION'] += radio_trans 322 else: 323 self.processed_input['TRANSITION'] = radio_trans

324 325 326

327 - def addRadioData(self,sn):

328 329 ''' 330 Add radio data to Transition() objects in all Star() objects. 331 332 Only done if RADIO_PATH is given and if a file named radio_data.db is 333 present in the given folder. 334 335 This method can be called multiple times per session for different stars 336 in which case the data in the transition objects will be replaced. 337 338 self.radio_trans maintains a list of "sample transitions" that contain 339 data. These function as data blueprints for transitions in the star_grid 340 341 @param sn: The star name for which to add data to the transitions 342 @type sn: str 343 344 ''' 345 346 #-- If not data available, reset all transitions so no data linger for a 347 # different star. 348 if not self.radio[sn]: 349 for star in self.star_grid: 350 for tr in star['GAS_LINES']: 351 tr.resetData() 352 return 353 354 #-- Extracted transitions (copies) are saved separately, so we only read 355 # data once at most. Hence we can continuously swap between star as we 356 # progress. 357 if not self.radio_trans[sn]: 358 #-- First extract all transitions and add/read data for them. 359 trans_sel = Transition.extractTransFromStars(self.star_grid,\ 360 dtype='resolved') 361 362 #-- Note that this is a list of copied transitions. 363 for trans in trans_sel: 364 trstr = trans.getInputString(include_nquad=0) 365 if trstr in self.radio[sn].keys(): 366 #-- No need to replace. Trans extraction already resets data 367 trans.addDatafile(self.radio[sn][trstr]) 368 trans.readData() 369 370 #-- self.radio_trans was defined in setRadio. 371 self.radio_trans[sn] = trans_sel 372 373 #-- We can now add these data to all Star() objects. 374 for star in self.star_grid: 375 for trans in self.radio_trans[sn]: 376 mtrans = star.getTransition(trans) 377 if mtrans: mtrans.setData(trans,replace=1)

378 379 380

381 - def setVarPars(self):

382 383 ''' 384 Define the list of variable parameters in this CC session. 385 386 ''' 387 388 self.var_pars = [k for k in self.multiplicative_grid.keys() + \ 389 self.additive_grid.keys() 390 if k[:5] != 'T_MAX']

391 392

393 - def readInput(self):

394 395 ''' 396 Read input for ComboCode and return list. 397 398 The MOLECULE, TRANSITION and R_POINTS_MASS_LOSS parameter formats are 399 checked for errors in this method. If erroneous, an IOError is raised. 400 401 ''' 402 403 multi_keys = ['MOLECULE','TRANSITION','R_POINTS_MASS_LOSS'] 404 input_dict = DataIO.readDict(self.inputfilename,convert_floats=1,\ 405 convert_ints=1,multi_keys=multi_keys) 406 #-- keywords in multi_keys require different method 407 self.processed_input = dict() 408 self.multiplicative_grid = dict() 409 self.additive_grid = dict() 410 molecules = input_dict.pop('MOLECULE',[]) 411 transitions = input_dict.pop('TRANSITION',[]) 412 r_points_mass_loss = input_dict.pop('R_POINTS_MASS_LOSS',[]) 413 for k,v in input_dict.items(): 414 #-- Fortran input is not case sensitive. Note: use the dict 415 # value v, not input_dict[k] because of this transformation. 416 k = k.upper() 417 #-- Determine delimiter 418 try: 419 if v.find('&') != -1: delimiter = '&' 420 elif v.find(';') != -1: delimiter = ';' 421 elif v.find(',') != -1: delimiter = ',' 422 elif v.find(':') != -1: delimiter = ':' 423 #-- * while no ; or , or : means multiple values for ONE model 424 # Only * star for multiplicative grid makes no sense (just 425 # give the value without delimiter) 426 elif v.find('*') != -1: delimiter = '&' 427 else: delimiter = ' ' 428 except AttributeError: 429 #-- v is already a float, so can't use .find on it => no grids 430 #-- no need to check the rest, continue on with the next k/v pair 431 self.processed_input[k] = v 432 continue 433 #-- Expanding '*' entries: Assumes the value is first, the count 434 # second. Can't be made flexible, because in some cases the value 435 # cannot be discerned from the count (because both are low-value 436 # integers) 437 newv = delimiter.join(\ 438 [len(value.split('*')) > 1 439 and delimiter.join([value.split('*')[0]]*\ 440 int(value.split('*')[1])) 441 or value 442 for value in v.split(delimiter)]) 443 #-- Add entries to processed_input, the multiplicative grid or the 444 #-- additive grid, depending on the type of delimiter. 445 if delimiter == ' ': 446 self.processed_input[k] = v 447 elif delimiter == ',': 448 newv = [float(value) 449 for value in newv.split(',')] 450 newv = Gridding.makeGrid(*newv) 451 self.multiplicative_grid[k] = newv 452 else: 453 try: 454 if delimiter == '&': 455 newv = tuple([float(value.rstrip()) 456 for value in newv.split('&')]) 457 self.processed_input[k] = newv 458 elif delimiter == ';': 459 newv = [value.rstrip() == '%' and '%' or float(value.rstrip()) 460 for value in newv.split(';')] 461 self.multiplicative_grid[k] = newv 462 elif delimiter == ':': 463 newv = [value.rstrip() == '%' and '%' or float(value.rstrip()) 464 for value in newv.split(':')] 465 self.additive_grid[k] = newv 466 except ValueError: 467 if delimiter == '&': 468 newv = tuple([value.rstrip() 469 for value in newv.split('&')]) 470 self.processed_input[k] = newv 471 elif delimiter == ';': 472 newv = [value.rstrip() 473 for value in newv.split(';')] 474 self.multiplicative_grid[k] = newv 475 elif delimiter == ':': 476 newv = [value.rstrip() 477 for value in newv.split(':')] 478 self.additive_grid[k] = newv 479 #-- Make sure R_POINTS_MASS_LOSS, Molecule and Transition input makes 480 #-- sense and is correct 481 if molecules: 482 molecules = DataIO.checkEntryInfo(molecules,20,'MOLECULE') 483 if isinstance(molecules,list): 484 self.additive_grid['MOLECULE'] = molecules 485 else: 486 self.processed_input['MOLECULE'] = molecules 487 if r_points_mass_loss: 488 r_points_mass_loss = DataIO.checkEntryInfo(r_points_mass_loss,4,\ 489 'R_POINTS_MASS_LOSS') 490 if isinstance(r_points_mass_loss,list): 491 self.additive_grid['R_POINTS_MASS_LOSS'] = r_points_mass_loss 492 else: 493 self.processed_input['R_POINTS_MASS_LOSS'] = r_points_mass_loss 494 if transitions: 495 nk = 12 496 497 #-- Check if N_QUAD is a grid parameter: remove any manual N_QUAD 498 # definitions to allow gridding to work. N_QUAD is taken from 499 # Star() by default in this case. 500 gkeys = self.multiplicative_grid.keys() + self.additive_grid.keys() 501 if 'N_QUAD' in gkeys: 502 transitions = [' '.join(tr.split()[:-1]) for tr in transitions] 503 nk = 11 504 505 #-- If N_QUAD is a gridding key, transitions will always be a tuple. 506 transitions = DataIO.checkEntryInfo(transitions,nk,'TRANSITION') 507 if isinstance(transitions,list): 508 self.additive_grid['TRANSITION'] = transitions 509 else: 510 self.processed_input['TRANSITION'] = transitions

511 512 513

514 - def getStars(self):

515 516 ''' 517 Return the list of Star() objects for this ComboCode session. 518 519 @return: The parameter Star() objects are returned. 520 @rtype: list[Star()] 521 522 ''' 523 524 return self.star_grid

525 526 527

528 - def createStarGrid(self):

529 530 ''' 531 Create a list of Star() objects based on the inputfile that has been 532 parsed with cc.readInput(). 533 534 The list of Star() objects is saved in self.star_grid, and is accessed 535 through cc.getStarGrid(). 536 537 ''' 538 539 base_star = Star.Star(example_star=self.processed_input,\ 540 path_gastronoom=self.path_gastronoom,\ 541 path_mcmax=self.path_mcmax,\ 542 print_check_t=self.print_check_t) 543 if self.additive_grid: 544 grid_lengths = [len(v) for v in self.additive_grid.values()] 545 if len(set(grid_lengths)) != 1: 546 raise IOError('The explicit parameter declaration using <:> '+\ 547 'has a variable amount of options (including ' +\ 548 'the R_GRID_MASS_LOSS definition). Aborting...') 549 else: 550 additive_dicts = [dict([(key,grid[index]) 551 for key,grid in self.additive_grid.items()]) 552 for index in xrange(grid_lengths[0])] 553 self.star_grid = [Star.Star(example_star=base_star,\ 554 path_gastronoom=self.path_gastronoom,\ 555 path_mcmax=self.path_mcmax,\ 556 extra_input=d,\ 557 print_check_t=self.print_check_t) 558 for d in additive_dicts] 559 else: 560 self.star_grid = [base_star] 561 for key,grid in self.multiplicative_grid.items(): 562 self.star_grid = [Star.Star(path_gastronoom=self.path_gastronoom,\ 563 path_mcmax=self.path_mcmax,\ 564 example_star=star,\ 565 extra_input=dict([(key,value)]),\ 566 print_check_t=self.print_check_t) 567 for star in self.star_grid 568 for value in grid] 569 for star in self.star_grid: 570 star.normalizeDustAbundances() 571 if self.processed_input.has_key('LAST_MCMAX_MODEL'): 572 del self.processed_input['LAST_MCMAX_MODEL'] 573 if self.processed_input.has_key('LAST_GASTRONOOM_MODEL'): 574 del self.processed_input['LAST_GASTRONOOM_MODEL'] 575 #-- Dust abundance is deleted as it is not changed during the session. 576 # Hence, it does not need to be re-updated after mutable input is 577 # removed. The keys need to be deleted to avoid inconsistencies in the 578 # star.normalizeDustAbundances() method. Some A_*** values may not be 579 # variable, while others are. Yet if any of them are variable, and 580 # have to be rescaled, then in principle they are ALL variable. This 581 # can create a mess, and therefore it is safer to just remove them 582 # from the input dictionary. 583 for akey in [k for k in self.processed_input.keys() if k[0:2] == 'A_']: 584 del self.processed_input[akey]

585 586 587

588 - def setOutputFolders(self):

589 590 ''' 591 Set the output folders. 592 593 If the folders do not already exist, they are created. 594 595 The locations are saved in cc.path for later use, but this is generally 596 only done inside a ComboCode session. Each module sets these themselves 597 598 ''' 599 600 paths = ['gastronoom','mcmax','chemistry'] 601 folders = [self.path_gastronoom,self.path_mcmax,self.path_chemistry] 602 names = ['gout','mout','cout'] 603 for p,f,n in zip(paths,folders,names): 604 if not getattr(self,p): continue 605 path = os.path.join(getattr(cc.path,p),f) 606 setattr(cc.path,n,path) 607 DataIO.testFolderExistence(path)

608 609 610

611 - def setVicManager(self):

612 613 ''' 614 Set up the VIC manager. 615 616 ''' 617 618 if self.vic and self.gastronoom and self.sphinx : 619 self.vic_manager = Vic.Vic(path=self.path_gastronoom,\ 620 account=self.vic_account,\ 621 time_per_sphinx=self.vic_time_per_sphinx,\ 622 credits_acc=self.vic_credits,\ 623 recover_sphinxfiles=self.recover_sphinxfiles) 624 if self.update_spec: 625 self.vic_manager.updateLineSpec() 626 else: 627 self.vic_manager = None

628 629

630 - def setModelManager(self):

631 632 ''' 633 Set up the model manager. 634 635 ''' 636 637 self.model_manager = MM(iterations=self.iterations,\ 638 processed_input=self.processed_input,\ 639 var_pars=self.var_pars,\ 640 path_gastronoom=self.path_gastronoom,\ 641 mcmax=self.mcmax,\ 642 gastronoom=self.gastronoom,\ 643 sphinx=self.sphinx,\ 644 iterative=self.plot_iterative,\ 645 num_model_sessions=len(self.star_grid),\ 646 vic_manager=self.vic_manager,\ 647 replace_db_entry=self.replace_db_entry,\ 648 path_mcmax=self.path_mcmax,\ 649 skip_cooling=self.skip_cooling,\ 650 recover_sphinxfiles=self.recover_sphinxfiles,\ 651 single_session=self.single_session,\ 652 )

653 654

655 - def setPlotManager(self):

656 657 ''' 658 Set up the plot manager(s) for each star name. 659 660 ''' 661 662 plot_pars = dict([(k,self.processed_input.pop(k)) 663 for k,v in self.processed_input.items() 664 if k[0:5] == 'PLOT_' or k[0:4] == 'CFG_']) 665 fn_add_star = plot_pars.pop('PLOT_FN_ADD_STAR',1) 666 self.plot_manager = {sn: PM(star_name=sn,\ 667 gastronoom=self.gastronoom,\ 668 mcmax=self.mcmax,\ 669 chemistry=self.chemistry,\ 670 path_gastronoom=self.path_gastronoom,\ 671 path_mcmax=self.path_mcmax,\ 672 path_chemistry=self.path_chemistry,\ 673 inputfilename=self.inputfilename,\ 674 pacs=self.pacs[sn],\ 675 spire=self.spire[sn],\ 676 sed=self.sed[sn],\ 677 fn_add_star=fn_add_star,\ 678 plot_pars=plot_pars) 679 for sn in self.star_name}

680 681 682

683 - def runModelManager(self):

684 685 ''' 686 Start up the modeling. 687 688 ''' 689 690 if self.gastronoom or self.mcmax: 691 print '***********************************' 692 print '** Starting grid calculation.' 693 print '***********************************' 694 for star_index, star in enumerate(self.star_grid): 695 print '***********************************' 696 print '** Model #%i out of %i requested models.'\ 697 %(star_index+1,len(self.star_grid)) 698 print '***********************************' 699 self.model_manager.startModeling(star,star_index) 700 #-- mline_done is True if in previous model an mline calculation 701 #-- was done: Only then do a progress check, because a lot of time 702 #-- has passed, but then a wait time is used to make sure the newly 703 #-- queued sphinx models after the mline model are properly queued. 704 if self.vic_manager \ 705 and self.vic_manager.getQueue() \ 706 and self.model_manager.mline_done_list[-1]: 707 print '***********************************' 708 print '** Current VIC queue:' 709 print self.vic_manager.getQueue() 710 self.vic_manager.checkProgress(wait_qstat=1) 711 712 if self.single_session: 713 if self.gastronoom: 714 self.model_manager.cool_db.sync() 715 self.model_manager.ml_db.sync() 716 self.model_manager.sph_db.sync() 717 if self.mcmax: 718 self.model_manager.mcmax_db.sync()

719 720 721

722 - def finalizeVic(self):

723 724 ''' 725 At the end of a modeling session, allow Vic to be finalized and clean up. 726 727 ''' 728 729 if not self.vic_manager is None: 730 #vic_running = vic_manager.checkProgress() 731 if self.vic_manager.getQueue(): 732 vic_running = True 733 else: 734 vic_running = False 735 while vic_running: 736 print 'VIC is not yet finished. Waiting 5 minutes before checking again.' 737 print self.vic_manager.getQueue() 738 try: 739 time.sleep(300) 740 except KeyboardInterrupt: 741 print 'Ending wait time, continuing with progress check immediately.' 742 vic_running = self.vic_manager.checkProgress() 743 self.vic_manager.finalizeVic() 744 745 if self.single_session: 746 self.model_manager.sph_db.sync()

747 748 749

750 - def runPlotManager(self):

751 752 ''' 753 Run the plotting manager. 754 755 ''' 756 757 #-- First check if any plots are requested at all 758 dds = [self.plot_manager[sn].gas_pars.keys() 759 for sn in self.star_name 760 if self.plot_manager[sn].gas_pars.keys()] + \ 761 [self.plot_manager[sn].dust_pars.keys() 762 for sn in self.star_name 763 if self.plot_manager[sn].dust_pars.keys()] + \ 764 [self.plot_manager[sn].chem_pars.keys() 765 for sn in self.star_name 766 if self.plot_manager[sn].chem_pars.keys()] 767 if not dds: return 768 769 #-- Continue with the plots 770 print '************************************************' 771 print '****** Plotting final results.' 772 print '************************************************' 773 for sn in self.star_name: 774 #-- Make sure resolved data files are correctly included. 775 self.addRadioData(sn) 776 print '** Plots for %s:'%sn 777 if self.plot_iterative: 778 print '** Plotting results for each iterative step.' 779 print '************' 780 #- star_grid_old remembers all old models if iterative is on, 781 #- meaning that every list in star_grid_old consists of Star models 782 #- associated with one level of iteration. Following line plots all 783 #- iterative steps for one star immutable parameter set 784 pm = self.plot_manager[sn] 785 for i in range(len(self.star_grid)): 786 pm.startPlotting(self.model_manager.star_grid_old[i],i+1) 787 self.plot_manager[sn].startPlotting(self.star_grid)

788 789 790

791 - def appendResults(self):

792 793 ''' 794 Append results at the end of the inputfile. 795 796 ''' 797 print '** Appending results to inputfile and copying to output folders.' 798 print '***********************************' 799 #-- Check if the transition was intended to be calculated, and if it was 800 #-- successful (ie don't add if it had already been done) 801 timestring = '%.4i-%.2i-%.2ih%.2i-%.2i-%.2i'\ 802 %(time.gmtime()[0],time.gmtime()[1],time.gmtime()[2],\ 803 time.gmtime()[3],time.gmtime()[4],time.gmtime()[5]) 804 appendage = [] 805 if self.model_manager.trans_bool_list: 806 model_ids_list = [list(set([(trans.molecule.molecule,\ 807 trans.getModelId()) 808 for boolean,trans in zip(trans_bool,\ 809 star['GAS_LINES']) 810 if trans.getModelId() \ 811 and (not trans_bool \ 812 or self.append_results)])) 813 for star,trans_bool in zip(self.star_grid,\ 814 self.model_manager.trans_bool_list)] 815 #-- all unique molecules over all stars 816 molec_list = list(set([molec 817 for model_ids in model_ids_list 818 for molec,model_id in model_ids 819 if model_ids])) 820 #-- all unique modelids for every star separately 821 model_id_unique = [list(set([model_id 822 for molec,model_id in model_ids])) 823 for model_ids in model_ids_list] 824 if [modelids for modelids in model_ids_list if modelids] != []: 825 appendage += \ 826 ['#########################################',\ 827 '## Successfully calculated transition model_ids on %s:'\ 828 %timestring] 829 appendage.extend(['## molecule %s'%molec 830 for molec in molec_list]) 831 for i,(star,model_ids) in enumerate(zip(self.star_grid,\ 832 model_ids_list)): 833 if model_ids: 834 appendage += ['## For Model %i : cooling id %s'\ 835 %(i+1,star['LAST_GASTRONOOM_MODEL'])] + \ 836 ['#molecule %s #%s' %(molecule,model_id) 837 for molecule,model_id in model_ids] + \ 838 ['########'] 839 for star,model_ids in zip(self.star_grid,model_id_unique): 840 for model_id in model_ids: 841 try: 842 i = 0 843 while True: 844 dummy = DataIO.readFile(\ 845 os.path.join(cc.path.gout,'models',model_id,\ 846 os.path.split(self.inputfilename)[1]+\ 847 '_%s_%i'%(model_id,i))) 848 i += 1 849 except IOError: 850 subprocess.call(['cp %s %s'%(self.inputfilename,\ 851 os.path.join(cc.path.gout,\ 852 'models',model_id,\ 853 os.path.split(self.inputfilename)[1]+\ 854 '_%s_%i'%(model_id,i)))],shell=True) 855 if self.model_manager.mcmax_done_list: 856 model_ids = [star['LAST_MCMAX_MODEL'] 857 for star,boolean in zip(self.star_grid,\ 858 self.model_manager.mcmax_done_list) 859 if boolean or self.append_results] 860 if model_ids: 861 appendage += ['#########################################',\ 862 '## MCMax model_ids associated with this grid on %s:'\ 863 %timestring] 864 appendage += ['#%s'%model_id for model_id in model_ids] 865 for model_id in model_ids: 866 try: 867 i = 0 868 while True: 869 dummy = DataIO.readFile(os.path.join(\ 870 cc.path.mout,'models',model_id,\ 871 os.path.split(self.inputfilename)[1]+\ 872 '_%s_%i'%(model_id,i))) 873 i += 1 874 except IOError: 875 subprocess.call(['cp %s %s'%(self.inputfilename,os.path.join(\ 876 cc.path.mout,'models',model_id,\ 877 os.path.split(self.inputfilename)[1]+\ 878 '_%s_%i'%(model_id,i)))],shell=True) 879 if appendage: DataIO.writeFile(filename=self.inputfilename,\ 880 input_lines=appendage+['\n'],mode='a')

881 882 883

884 - def runStatistics(self):

885 886 ''' 887 Run the statistics module. 888 889 ''' 890 891 self.pacsstats = dict() 892 self.spirestats = dict() 893 self.unresostats = [] 894 self.resostats = dict() 895 self.sedstats = dict() 896 self.chemstats = dict() 897 898 #-- Data are handled by these objects themselves for unresolved lines. 899 for sn in self.star_name: 900 if self.statistics and not self.pacs[sn] is None: 901 ss = UnresoStats.UnresoStats(star_name=sn,\ 902 path_code=self.path_gastronoom) 903 ss.setInstrument(instrument_name='PACS',\ 904 instrument_instance=self.pacs[sn],\ 905 stat_method=self.stat_method) 906 self.pacsstats[sn] = ss 907 self.unresostats.append(ss) 908 if self.statistics and not self.spire[sn] is None: 909 ss = UnresoStats.UnresoStats(star_name=sn,\ 910 path_code=self.path_gastronoom) 911 ss.setInstrument(instrument_name='SPIRE',\ 912 instrument_instance=self.spire[sn],\ 913 stat_method=self.stat_method) 914 self.spirestats[sn] = ss 915 self.unresostats.append(ss) 916 917 for ss in self.unresostats: 918 instr = ss.instrument.instrument.upper() 919 print '************************************************' 920 print '** %s statistics for %s.'%(instr,ss.star_name) 921 print '************************************************' 922 ss.setModels(star_grid=self.star_grid) 923 ss.setLineStrengths() 924 ss.calcChiSquared(chi2_method=self.stat_chi2) 925 ss.plotRatioWav(inputfilename=self.inputfilename) 926 927 for sn in self.star_name: 928 if self.statistics and not self.sed[sn] is None: 929 print '************************************************' 930 print '** SED statistics for %s.'%sn 931 print '************************************************' 932 ss = SedStats.SedStats(star_name=sn,path_code=self.path_mcmax) 933 ss.setInstrument(sed=self.sed[sn]) 934 ss.setModels(star_grid=self.star_grid) 935 ss.setModelPhotometry() 936 ss.calcChi2(chi2_method=self.stat_chi2) 937 self.sedstats[sn] = ss 938 939 if self.statistics and self.radio[sn]: 940 #-- Make sure all resolved data properties are set. 941 self.addRadioData(sn) 942 if not self.radio_trans[sn]: 943 return 944 print '************************************************' 945 print '** Statistics for spectrally resolved lines in %s.'%sn 946 print '** Use cc_session.resostats[sn] for interactive tools.' 947 print '************************************************' 948 ss = ResoStats.ResoStats(star_name=sn,\ 949 path_code=self.path_gastronoom,\ 950 lll_p=self.stat_lll_p,\ 951 vmin=self.stat_lll_vmin,\ 952 vmax=self.stat_lll_vmax) 953 ss.setInstrument(self.radio_trans[sn]) 954 ss.setModels(star_grid=self.star_grid) 955 ss.setIntensities() 956 if self.stat_print: ss.printStats() 957 self.resostats[sn] = ss 958 #bfms = self.resostats.selectBestFitModels(mode='int') 959 #self.plot_manager.plotTransitions(star_grid=bfms,fn_suffix='BFM',force=1) 960 961 for sn in self.star_name: 962 if self.statistics and self.chemistry: 963 self.chemstats_molecules = list(self.chemstats_molecules) 964 ss = ChemStats.ChemStats(star_name=sn,\ 965 path_code=self.path_chemistry,\ 966 star_grid=self.star_grid,\ 967 molecules = self.chemstats_molecules) 968 self.chemstats[sn] = ss

969 970 971

972 - def runChemistry(self):

973 if self.chemistry: 974 print '************************************************' 975 print '** Running Chemistry ' 976 print '************************************************' 977 978 chemistry_db_path = os.path.join(cc.path.cout,'Chemistry_models.db') 979 self.chem_db = Database.Database(db_path=chemistry_db_path) 980 981 ch = Chemistry.Chemistry(path_chemistry=self.path_chemistry,\ 982 replace_db_entry = self.replace_db_entry,\ 983 db = self.chem_db,\ 984 single_session=self.single_session) 985 for i,star in enumerate(self.star_grid): 986 print '***********************************' 987 print '** Model #%i out of %i requested models.'\ 988 %(i+1,len(self.star_grid)) 989 print '***********************************' 990 991 ch.doChemistry(star) 992 993 if ch.model_id: 994 star['LAST_CHEMISTRY_MODEL'] = ch.model_id

995 996 997

998 - def doContDiv(self):

999 1000 ''' 1001 Run the Continuum Division class for this CC session. 1002 1003 ''' 1004 1005 if isinstance(self.contdiv_features,str): 1006 self.contdiv_features = [self.contdiv_features] 1007 for k in self.contdiv_features: 1008 features = ['MGS','H2O3.1'] 1009 if k in features: 1010 print '** Plotting continuum division for the %s feature.'%k 1011 all_franges = [[20.0,23.,46.,48.5],\ 1012 [2.6,2.85,3.3,3.7]] 1013 all_funcs = ['linear',\ 1014 'power'] 1015 franges = all_franges[features.index(k)] 1016 func = all_funcs[features.index(k)] 1017 self.contdiv = ContinuumDivision.ContinuumDivision(\ 1018 star_grid=self.star_grid,\ 1019 spec=[self.sed],franges=franges,\ 1020 plot=self.show_contdiv,func=func,\ 1021 cfg=self.cfg_contdiv) 1022 self.contdiv.prepareModels() 1023 self.contdiv.prepareData() 1024 self.contdiv.show() 1025 for key,val in self.contdiv.eq_width.items(): 1026 if key[0:3] == 'sws': 1027 key = self.sed.star_name_plots 1028 print 'Equivalent width for %s in %s: %f'\ 1029 %(k,self.sed.star_name_plots,val) 1030 for star in self.star_grid: 1031 model_id = star['LAST_MCMAX_MODEL'] 1032 print 'Equivalent width for %s in %s: %f'\ 1033 %(k,model_id,self.contdiv.eq_width[model_id])

1034 1035 1036

1037 - def printStarInfo(self):

1038 1039 ''' 1040 Print extra Star() info to the shell. 1041 1042 ''' 1043 1044 if len(self.star_grid)<20 and self.print_model_info: 1045 print '************************************************' 1046 for star in self.star_grid: 1047 if star['LAST_MCMAX_MODEL']: 1048 print 'Requested MCMax parameters for %s:'\ 1049 %star['LAST_MCMAX_MODEL'] 1050 if star.has_key('R_INNER_DUST'): del star['R_INNER_DUST'] 1051 print '%s = %f R_STAR (effective inner radius)'%('R_INNER_DUST',star['R_INNER_DUST']) 1052 print '%s = %s'%('TEMDUST_FILENAME',star['TEMDUST_FILENAME']) 1053 print '%s = %s'%('DUST_TEMPERATURE_FILENAME',star['DUST_TEMPERATURE_FILENAME']) 1054 print '%s = %s km/s'%('V_EXP_DUST',star['V_EXP_DUST']) 1055 1056 if not int(star['MRN_DUST']): 1057 cdcalc = ColumnDensity.ColumnDensity(star) 1058 dlist = [sp 1059 for sp in star.getDustList() 1060 if 'H2O' not in sp] 1061 print 'Dust FULL column densities (if available):' 1062 for species in dlist: 1063 print 'C_%s = %.3e g/cm^-2'\ 1064 %(species,cdcalc.dustFullColDens(species)) 1065 print 'Dust FULL number column densities (if available):' 1066 for species in dlist: 1067 print 'N_%s = %.3e cm^-2'\ 1068 %(species,cdcalc.dustFullNumberColDens(species)) 1069 print 'Dust associated molecular abundances (if available):' 1070 print 'Note: Does not use above column densities. See ColumnDensity.py.' 1071 for species in dlist: 1072 print 'A_%s/A_H2 = %.3e'\ 1073 %(species,cdcalc.dustMolecAbun(species)) 1074 if star.has_key('T_DES_H2O') or star.has_key('T_DES_CH2O')\ 1075 or star.has_key('T_DES_AH2O') \ 1076 and not int(star['MRN_DUST']): 1077 print ', '.join(['%s = %s K'%(ts,star[ts]) 1078 for ts in ['T_DES_H2O','T_DES_CH2O',\ 1079 'T_DES_AH2O'] 1080 if star.has_key(ts)]) 1081 print ', '.join(['%s = %.2f R_STAR = %.2e cm'\ 1082 %(rs,star[rs],\ 1083 star[rs]*star['R_STAR']*star.Rsun) 1084 for rs in ['R_DES_H2O','R_DES_AH2O',\ 1085 'R_DES_CH2O'] 1086 if star.has_key(rs)]) 1087 cdh2o1 = cdcalc.dustFullColDens('CH2O') 1088 cdh2o2 = cdcalc.dustFullColDens('AH2O') 1089 print 'The FULL water ice column density:' 1090 print 'C_H2O_ice = %.3e g cm-2'%(cdh2o1+cdh2o2) 1091 ncdh2o1 = cdcalc.dustFullNumberColDens('CH2O') 1092 ncdh2o2 = cdcalc.dustFullNumberColDens('AH2O') 1093 print 'The FULL water ice column number density:' 1094 print 'N_H2O_ice = %.3e cm-2'%(ncdh2o1+ncdh2o2) 1095 nh2o1 = cdcalc.dustMolecAbun('CH2O') 1096 nh2o2 = cdcalc.dustMolecAbun('AH2O') 1097 print 'The associated molecular abundance of water ice:' 1098 print 'A_H2O_ice/A_H2 = %.3e'%(nh2o1+nh2o2) 1099 print '' 1100 if star['LAST_GASTRONOOM_MODEL']: 1101 print 'Requested GASTRoNOoM parameters for %s:'%star['LAST_GASTRONOOM_MODEL'] 1102 print '%s = %s'%('DENSFILE',star['DENSFILE']) 1103 print '%s = %f Rsun = %f AU'%('R_STAR',star['R_STAR'],star['R_STAR']*star.Rsun/star.au) 1104 print '%s = %f R_STAR'%('R_OUTER_GAS',star['R_OUTER_EFFECTIVE']) 1105 print '%s = %f R_STAR'%('R_INNER_GAS',star['R_INNER_GAS']) 1106 print '%s = %f'%('DUST_TO_GAS_INITIAL',star['DUST_TO_GAS_INITIAL']) 1107 print '%s = %f'%('DUST_TO_GAS_ITERATED',star['DUST_TO_GAS_ITERATED']) 1108 print '%s (semi-empirical) = %f'%('DUST_TO_GAS',star['DUST_TO_GAS']) 1109 if star['DUST_TO_GAS_CHANGE_ML_SP']: 1110 print '%s = %s'%('DUST_TO_GAS_CHANGE_ML_SP',star['DUST_TO_GAS_CHANGE_ML_SP']) 1111 print '%s = %f'%('[N_H2O/N_H2]',star['F_H2O']) 1112 print '%s = %.2f R_STAR = %.2e cm'%('R_OH1612_NETZER',star['R_OH1612_NETZER'],star['R_OH1612_NETZER']*star.Rsun*star['R_STAR']) 1113 if star['R_OH1612_AS']: 1114 print '%s = %.2f R_STAR = %.2e cm'%('R_OH1612_OBS',star['R_OH1612'],star['R_OH1612']*star.Rsun*star['R_STAR']) 1115 print '-----------------------------------' 1116 #if star.has_key('R_DES_H2O') or star.has_key('R_DES_CH2O') or star.has_key('R_DES_AH2O'): 1117 #(nh2o,nh2o_ice,nh2,nh2o_full,nh2_full) = wa.getWaterInfo(star) 1118 #print 'Ice shell water VAPOUR column density [cm-2]:' 1119 #print '%.3e'%nh2o 1120 #print 'Ice shell H2 column density [cm-2]:' 1121 #print '%.3e'%nh2 1122 #print 'Total water vapour abundance (ortho + para) wrt H2:' 1123 #print '%.3e'%(nh2o/nh2) 1124 #print 'Ice shell water ICE MOLECULAR column density [cm-2]:' 1125 #print '%.3e'%nh2o_ice 1126 #print 'Minimum required water vapour abundance N(H2O)/N(H2) for this ice column density:' 1127 #print '%.3e'%(nh2o_ice/nh2) 1128 #print 'Ice/vapour fraction in ice shell:' 1129 #print '%.2f'%(nh2o_ice/nh2o) 1130 #print 1131 #print 'FULL shell water VAPOUR column density [cm-2]:' 1132 #print '%.3e'%nh2o_full 1133 #print 'Total water vapour abundance (ortho + para) wrt H2:' 1134 #print '%.3e'%(nh2o_full/nh2_full) 1135 #print 1136 #if not int(star['MRN_DUST'): 1137 1138 #else: 1139 #print 'No water ice present in dust model.' 1140 print '************************************************'

Source Code for Module ComboCode.cc.ComboCode