ivs.sed.fit

1 # -*- coding: utf-8 -*- 2 """ 3 Fit SED models to observed data using various approaches. 4 """ 5 import logging 6 import sys 7 import astropy.io.fits as pf 8 import itertools 9 import re 10 import copy 11 12 import numpy as np 13 from numpy import inf 14 import scipy 15 from scipy.interpolate import Rbf 16 from scipy.optimize import fmin,fmin_powell 17 from ivs.statistics import pca 18 from ivs.sed import model 19 from ivs.sed import filters 20 from ivs.sed.decorators import iterate_gridsearch,parallel_gridsearch 21 from ivs.sigproc import fit as sfit 22 from ivs.aux import numpy_ext 23 from ivs.aux import progressMeter 24 from ivs.aux.decorators import make_parallel 25 from ivs.units import constants 26 27 logger = logging.getLogger("SED.FIT")

28 29 #{ PCA functions 30 31 -def get_PCA_grid(colors,res=3,teffrange=(-np.inf,np.inf),loggrange=(-np.inf,np.inf), 32 ebvrange=(-np.inf,np.inf),**kwargs):

33 """ 34 Transform a flux grid to a colour grid. 35 36 The resulting grid is actually log10(flux ratio) instead of flux ratio! This 37 works better in the PCA. 38 39 Extra keyword arguments can be used to set the atmosphere model. 40 41 @parameter colors: list of desired colours (['GENEVA.U-B','STROMGREN.M1'...]) 42 @type colors: list of strings 43 @keyword res: resolution in E(B-V) 44 @type res: integer 45 @keyword teffrange: range of Teffs to use 46 @type teffrange: tuple 47 @keyword loggrange: range of loggs to use 48 @type loggrange: tuple 49 @keyword ebvrange: range of E(B-V)s to use 50 @type ebvrange: tuple 51 @return: log10(colorgrid), (teffs,loggs,ebvs) 52 @rtype: N-D numpy array,(1Darray,1Darray,1Darray) 53 """ 54 #-- read in the color parameters from the FITS file 55 gridfile = model.get_file(integrated=True,**kwargs) 56 ff = pf.open(gridfile) 57 ext = ff[1] 58 59 teff = ext.data.field('teff') 60 logg = ext.data.field('logg') 61 ebv = ext.data.field('ebv') 62 keep = (ebvrange[0]<=ebv) & (ebv<=ebvrange[1]) 63 keep = keep & (teffrange[0]<=teff) & (teff<=teffrange[1]) 64 keep = keep & (loggrange[0]<=logg) & (logg<=loggrange[1]) 65 teff,logg,ebv = teff[keep],logg[keep],ebv[keep] 66 67 A = model._get_flux_from_table(ext,colors) 68 A = A[keep] 69 70 #-- order and set the resolution of the grid 71 B = np.vstack([teff.T,logg.T,ebv.T,A.T]).T 72 B = numpy_ext.sort_order(B,[0,1,2]) 73 B = B[::res] 74 logger.info('Calculated color grid for PCA (DIM=%s, using %s)'%(B[:,3:].shape,' '.join(colors))) 75 return np.log10(B[:,3:]),(B[:,0],B[:,1],B[:,2])

76

77 78 79 80 81 82 83 84 -def get_PCA(A):

85 """ 86 Find the principal components of a color grid. 87 88 @param A: color grid obtained via C{get_PCA_grid}. 89 @type A. numpy N-D array 90 @return: PCA loadings, PCA scores, (column means, column standard devs) 91 @rtype: N-D array, N-D array, (1D array, 1D array) 92 """ 93 #-- relocate/standardize the color grid 94 means = A.mean(axis=0) 95 stds = A.std(axis=0) 96 A_st = (A-means)/stds 97 #-- find the principal components 98 # T = PCA scores 99 # P = PCA loadings 100 T,P,explained_var = pca.PCA_svd(A_st,standardize=False) 101 #-- report the explained variance per axis (and a maximum of four to show) 102 ev = explained_var*100 103 logger.info("PCA: Explained variance: %s"%(' '.join(['Axis%d=%.2f%%'%(i,j) for i,j in enumerate(ev[:4])]))) 104 return P,T,(means,stds)

105

106 107 108 109 110 111 112 113 114 -def calibrate_PCA(T,pars,function='linear'):

115 """ 116 Define the interpretations of the principal components. 117 118 T are the PCA scores, pars a list of axes (e.g. teff, logg and ebv). 119 120 @param T: PCA scores 121 @type T: N-D array 122 @param pars: real-life axes 123 @type pars: list of 1D arrays 124 @return: Rbf interpolating function 125 @rtype: n-tuple 126 """ 127 D = len(pars) # Dimension of parameter space 128 calib = [] 129 for i in range(D): 130 args = [T[:,j] for j in range(D)] + [pars[i],function] 131 calib.append(Rbf(function='linear',*args[:-1])) 132 logger.info('Calibrated first %d axes of PCA (of total %d)'%(len(pars),T.shape[1])) 133 return tuple(calib)

134

135 136 137 138 139 140 141 -def get_PCA_parameters(obsT,calib,P,means,stds,e_obsT=None,mc=None):

142 """ 143 Derive fundamental parameters of a sample of observations given a PCA. 144 145 Monte Carlo simulations are only available when you give 1 target. 146 147 @param obsT: observed colours of the sample 148 @type obsT: list of lists 149 @param calib: PCA calibration obtained via C{calibrate_PCA} 150 @type calib: list of Rbf functions 151 @param P: PCA loadings 152 @type P: N-D array 153 @param means: means of PCA grid 154 @type means: numpy array 155 @param stds: standard deviations of PCA grid 156 @type stds: numpy array 157 @param mc: number of MonteCarlo simulations 158 @type mc: integer 159 @return: list of fundamental parameters 160 @rtype: list of lists 161 """ 162 #-- put obsT in same format as is used for PCA 163 obsT = np.asarray(obsT) 164 obsT = np.log10(obsT) 165 obsT = np.dot((obsT-means)/stds,P.T) 166 167 #-- this function is made for a sample of observations: if only one 168 # observations is given, make it appear to be part of a sample 169 if not len(obsT.shape)==2: 170 obsT = np.array([obsT]) 171 172 #-- if we want MC simulations, expand the array 173 if mc is not None: 174 if e_obsT is None: 175 e_obsT = 0.01*obsT[0] 176 obsT_ = np.array([obsT[0]+np.random.normal(size=len(obsT[0]),scale=e_obsT) for i in xrange(mc)]) 177 obsT_[0] = obsT[0] 178 obsT = obsT_ 179 180 #-- prepare output of parameters 181 pars = np.zeros((len(obsT),len(calib))) 182 for i in range(len(calib)): 183 args = [obsT[:,j] for j in range(len(calib))] 184 pars[:,i] = calib[i](*args) 185 return pars

186

187 188 #} 189 190 #{ Grid search 191 192 -def stat_chi2(meas,e_meas,colors,syn,full_output=False, **kwargs):

193 """ 194 Calculate Chi2 and compute angular diameter. 195 196 Colors and absolute fluxes are used to compute the Chi2, only absolute 197 fluxes are used to compute angular diameter. If no absolute fluxes are 198 given, the angular diameter is set to 0. 199 200 @param meas: array of measurements 201 @type meas: 1D array 202 @param e_meas: array containing measurements errors 203 @type e_meas: 1D array 204 @param colors: boolean array separating colors (True) from absolute fluxes (False) 205 @type colors: 1D boolean array 206 @param syn: synthetic fluxes and colors 207 @type syn: 1D array 208 @param full_output: set to True if you want individual chisq 209 @type full_output: boolean 210 @return: chi-square, scale, e_scale 211 @rtype: float,float,float 212 """ 213 #-- if syn represents only one measurement 214 if len(syn.shape)==1: 215 if sum(~colors) > 0: 216 if 'distance' in kwargs: 217 scale = 1/kwargs['distance']**2 218 e_scale = scale / 100 219 else: 220 ratio = (meas/syn)[~colors] 221 weights = (meas/e_meas)[~colors] 222 #-- weighted average and standard deviation 223 scale = np.average(ratio,weights=weights) 224 #print 'bla',weights.shape,ratio.shape,scale 225 e_scale = np.sqrt(np.dot(weights, (ratio-scale)**2)/weights.sum()) 226 else: 227 scale,e_scale = 0,0 228 #-- we don't need to scale the colors, only the absolute fluxes 229 chisq = np.where(colors, (syn-meas)**2/e_meas**2, (syn*scale-meas)**2/e_meas**2) 230 if full_output: 231 return chisq,meas/syn,meas/e_meas 232 else: 233 return chisq.sum(),scale,e_scale 234 #-- if syn is many measurements, we need to vectorize this: 235 else: 236 if sum(~colors) > 0: 237 if 'distance' in kwargs: 238 scale = 1/kwargs['distance']**2 239 scale = np.ones_like(syn[~colors][0,:]) * scale 240 e_scale = scale / 100 241 else: 242 ratio = (meas/syn)[~colors] 243 weights = (meas/e_meas)[~colors] 244 #-- weighted average and standard deviation 245 scale = np.average(ratio,weights=weights.reshape(-1),axis=0) 246 e_scale = np.sqrt(np.dot(weights.T, (ratio-scale)**2)/weights.sum(axis=0))[0] 247 #scale = np.average(ratio,axis=0) 248 #e_scale = np.zeros_like(scale) 249 else: 250 scale,e_scale = np.zeros(syn.shape[1]),np.zeros(syn.shape[1]) 251 #-- we don't need to scale the colors, only the absolute fluxes 252 chisq = np.where(colors.reshape(-1,1), (syn-meas)**2/e_meas**2, (syn*scale-meas)**2/e_meas**2) 253 if full_output: 254 return chisq,meas/syn,meas/e_meas 255 else: 256 return chisq.sum(axis=0),scale,e_scale

257

258 259 -def generate_grid_single_pix(photbands, points=None, clear_memory=True, **kwargs):

260 """ 261 Generate a grid of parameters. 262 """ 263 264 #-- Find the parameters provided and store them separately. 265 ranges, parameters = {}, [] 266 for key in kwargs.keys(): 267 if re.search('range$', key): 268 ranges[key] = kwargs.pop(key) 269 parameters.append(re.sub('range$', '', key)) 270 271 #-- report on the received grid 272 if not kwargs: 273 logger.info('Received grid (%s)'%model.defaults2str()) 274 else: 275 logger.info('Received custom grid (%s)'%kwargs) 276 277 #-- get the pixelgrid 278 axis_values,gridpnts,flux,colnames = \ 279 model._get_pix_grid(photbands,teffrange=(-inf,inf), 280 loggrange=(-inf,inf),ebvrange=(-inf,inf), 281 zrange=(-inf,inf),rvrange=(-inf,inf),vradrange=(0,0), 282 include_Labs=True,clear_memory=clear_memory,**kwargs) 283 284 #-- we first generate random teff-logg coordinates, since the grid is 285 # not exactly convex in these parameters. We assume it is for all the 286 # other parameters. We need to extract the teff/logg points, but make them 287 # unique 288 colnames = list(colnames) 289 teff_index = colnames.index('teff') 290 logg_index = colnames.index('logg') 291 teffs,loggs = gridpnts[:,teff_index],gridpnts[:,logg_index] 292 293 #-- get ranges for teff and logg 294 teffrange = ranges.pop('teffrange', (-inf,inf)) 295 loggrange = ranges.pop('loggrange', (-inf,inf)) 296 correctTeff, correctLogg = False, False 297 if teffrange[0] == teffrange[1]: 298 teffrange = [teffrange[0], teffrange[0]+1] 299 correctTeff = True 300 if loggrange[0] == loggrange[1]: 301 loggrange = [loggrange[0], loggrange[0]+0.01] 302 correctLogg = True 303 304 #-- we need to cut the grid to fit the teff and logg range: we replace the 305 # values for the upper and lower limit in the grid with those from the 306 # given ranges. This is a bit elaborate, but I don't see a better way 307 # of doin' it. 308 teffl_index = max(np.searchsorted(axis_values[teff_index],teffrange[0])-1,0) 309 teffu_index = min(np.searchsorted(axis_values[teff_index],teffrange[1]),len(axis_values[teff_index])-1) 310 teff_lower = axis_values[teff_index][teffl_index] 311 teff_upper = axis_values[teff_index][teffu_index] 312 cut = (teffs<teff_lower) | (teff_upper<teffs) 313 if teff_lower<teffrange[0]: teffs[teffs==teff_lower] = teffrange[0] 314 if teff_upper>teffrange[1]: teffs[teffs==teff_upper] = teffrange[1] 315 316 loggl_index = max(np.searchsorted(axis_values[logg_index],loggrange[0])-1,0) 317 loggu_index = min(np.searchsorted(axis_values[logg_index],loggrange[1]),len(axis_values[logg_index])-1) 318 logg_lower = axis_values[logg_index][loggl_index] 319 logg_upper = axis_values[logg_index][loggu_index] 320 cut = cut | (loggs<logg_lower) | (logg_upper<loggs) 321 if logg_lower<loggrange[0]: loggs[loggs==logg_lower] = loggrange[0] 322 if logg_upper>loggrange[1]: loggs[loggs==logg_upper] = loggrange[1] 323 teffs = teffs[~cut] 324 loggs = loggs[~cut] 325 326 #-- Generate a grid in logg/teff keeping in mind that this is not a rectangular space 327 gridpnts_ = numpy_ext.unique_arr(np.column_stack([teffs,loggs])) 328 329 #-- now we can generate random points: 330 sample1 = numpy_ext.random_rectangular_grid(gridpnts_,points) 331 if correctTeff: sample1[:,0] = np.array([teffrange[0] for i in sample1[:,0]]) 332 if correctLogg: sample1[:,1] = np.array([loggrange[0] for i in sample1[:,1]]) 333 334 for name in colnames: 335 if not name+'range' in ranges: ranges[name+'range'] = (-inf, inf) 336 sample2 = np.random.uniform(low =[max(ax.min(),ranges[name+'range'][0]) for ax,name in zip(axis_values,colnames) if not name in ['teff','logg']],\ 337 high=[min(ax.max(),ranges[name+'range'][1]) for ax,name in zip(axis_values,colnames) if not name in ['teff','logg']],\ 338 size=((len(sample1),len(colnames)-2))) 339 340 colnames.remove('teff') 341 colnames.remove('logg') 342 #-- return grid and column names 343 out_dict_ = {} 344 for col,name in zip(np.column_stack([sample1,sample2]).T,['teff','logg']+colnames): 345 out_dict_[name] = col 346 347 #-- Check if all collumns that were provided are also returned 348 out_dict = {} 349 for name in parameters: 350 if name in out_dict_: 351 out_dict[name] = out_dict_[name] 352 else: 353 out_dict[name] = np.array([ranges[name+'range'][0] for i in out_dict['teff']]) 354 355 return out_dict

356

357 358 -def generate_grid_pix(photbands, points=None, clear_memory=False,**kwargs):

359 """ 360 Generate a grid of parameters for 1 or more stars. Based on the generate_grid_single_pix 361 method. The radius of the components is based on the masses if given, otherwise on the 362 radrange arguments. If masses are given the radrange arguments are ignored, meaning that 363 the returned radius can be outside those limits. If only 1 component is provided no radius 364 will be returned. 365 366 This function can handle multiple components in the same way as a single component. 367 parameter ranges are provided as <parname><component>range=(...,...). fx: 368 teffrange = (5000, 10000), loggrange = (3.5, 4.5), teff2range = (10000, 20000), 369 logg2range = (5.5, 6.0) 370 371 For the first (or only) component both no component number (teffrange) or 1 as component 372 number (teff1range) can be used. 373 374 Returns a dictionary with for each parameter that has a range provided, an array of values. 375 In the case of multiple components, the radius will be returned even is no radius ranges 376 are provided. 377 """ 378 379 #-- Find all ranges and the number of components 380 radiusrange = [] 381 ranges, parameters, components = {}, set(), set() 382 for key in kwargs.keys(): 383 if re.search('range$', key) and not re.search('^rad\d?',key): 384 ranges[key] = kwargs.pop(key) 385 name, comp = re.findall('(.*?)(\d?)range$', key)[0] 386 parameters.add(name) 387 components.add(comp) 388 elif re.search('range$', key) and re.search('^rad\d?',key): 389 radiusrange.append(kwargs.pop(key)) 390 391 #-- If only one component we can directly return the grid 392 if len(components) == 1: 393 kwargs.update(ranges) 394 return generate_grid_single_pix(photbands, points=points, clear_memory=clear_memory, **kwargs) 395 396 #-- For each component get the grid from grid_single_pix 397 pars, npoints = {}, +inf 398 for i, (comp, grid) in enumerate(zip(components, model.defaults_multiple)): 399 ranges_ = {} 400 for par in parameters: 401 ranges_[par+'range'] = ranges[par+comp+'range'] if par+comp+'range' in ranges else ranges[par+'range'] 402 403 kwargs.update(ranges_) 404 kwargs.update(grid) 405 grid_ = generate_grid_single_pix(photbands, points=points, clear_memory=clear_memory, **kwargs) 406 407 #-- prepare a permutation so different blocks are not clustered together 408 permutation = np.random.permutation(len(grid_['teff'])) 409 410 for key in grid_.keys(): 411 npoints = min(npoints,len(grid_[key])) 412 pars[key+comp] = grid_[key][permutation] 413 414 #-- The generate_grid_single_pix method does not guarantee the number of points. 415 # So force that all arrays have the same length. 416 for key in pars.keys(): 417 pars[key] = pars[key][:npoints] 418 419 #-- Check that ebv, z and rv is the same for each component 420 for comp in components: 421 if 'ebv' in parameters: pars['ebv'+comp] = pars['ebv'] 422 if 'z' in parameters: pars['z'+comp] = pars['z'] 423 if 'rv' in parameters: pars['rv'+comp] = pars['rv'] 424 425 #-- Check if we are dealing with a binary or not and set the radii accordingly 426 if 'masses' in kwargs and kwargs['masses'] != None: 427 #-- The radius of the stars is calculated based on logg and the provided masses 428 masses = kwargs['masses'] 429 G, Msol, Rsol = constants.GG_cgs, constants.Msol_cgs, constants.Rsol_cgs 430 for i, comp in enumerate(components): 431 pars['rad'+comp] = np.sqrt(G*masses[i]*Msol/10**pars['logg'+comp])/Rsol 432 else: 433 #-- We have random different radii for the stars 434 if radiusrange == []: radiusrange = [(0.1,10) for i in components] 435 for i, comp in enumerate(components): 436 pars['rad'+comp] = np.random.uniform(low=radiusrange[i][0], high=radiusrange[i][1], size=npoints) 437 438 return pars

439

440 441 -def generate_grid_single(photbands,teffrange=(-inf,inf),loggrange=(-inf,inf), 442 ebvrange=(-inf,inf),zrange=(-inf,inf), 443 points=None,res=None,clear_memory=True,**kwargs):

444 """ 445 Generate grid points at which to fit an interpolated grid of SEDs. 446 447 If C{points=None}, the points are chosen on the predefined grid points. 448 Otherwise, C{points} grid points will be generated, uniformly distributed 449 between the ranges defined by C{teffrange}, C{loggrange} and C{ebvrange}. If 450 you set the resolution to C{2}, one out of every two points will be selected. 451 452 Extra keyword arguments can be used to give more details on the atmosphere 453 models to use. 454 455 Colors are automatically detected. 456 457 You can fix one parameter e.g. via setting teffrange=(10000,10000). 458 459 >>> photbands = ['GENEVA.G','GENEVA.B-V'] 460 461 Start a figure: 462 463 >>> p = pl.figure() 464 >>> rows,cols = 2,4 465 466 On the grid points, but only one in every 100 points (otherwise we have over 467 a million points): 468 469 >>> teffs,loggs,ebvs,zs = generate_grid(photbands,res=100) 470 471 >>> p = pl.subplot(rows,cols,1) 472 >>> p = pl.scatter(teffs,loggs,c=ebvs,s=(zs+5)*10,edgecolors='none',cmap=pl.cm.spectral) 473 >>> p = pl.xlim(pl.xlim()[::-1]);p = pl.ylim(pl.ylim()[::-1]) 474 >>> p = pl.xlabel('Teff');p = pl.ylabel('Logg') 475 476 >>> p = pl.subplot(rows,cols,1+cols) 477 >>> p = pl.scatter(ebvs,zs,c=teffs,s=(loggs+2)*10,edgecolors='none',cmap=pl.cm.spectral) 478 >>> p = pl.xlabel('E(B-V)');p = pl.ylabel('Z') 479 480 Randomly distributed over the grid's ranges: 481 482 >>> teffs,loggs,ebvs,zs = generate_grid(photbands,points=10000) 483 484 >>> p = pl.subplot(rows,cols,2) 485 >>> p = pl.scatter(teffs,loggs,c=ebvs,s=(zs+5)*10,edgecolors='none',cmap=pl.cm.spectral) 486 >>> p = pl.xlim(pl.xlim()[::-1]);p = pl.ylim(pl.ylim()[::-1]) 487 >>> p = pl.xlabel('Teff');p = pl.ylabel('Logg') 488 489 >>> p = pl.subplot(rows,cols,2+cols) 490 >>> p = pl.scatter(ebvs,zs,c=teffs,s=(loggs+2)*10,edgecolors='none',cmap=pl.cm.spectral) 491 >>> p = pl.xlabel('E(B-V)');p = pl.ylabel('Z') 492 493 Confined to a small area in the grid's range: 494 495 >>> teffs,loggs,ebvs,zs = generate_grid(photbands,teffrange=(8000,10000),loggrange=(4.1,4.2),zrange=(0,inf),ebvrange=(1.,2),points=10000) 496 497 >>> p = pl.subplot(rows,cols,3) 498 >>> p = pl.scatter(teffs,loggs,c=ebvs,s=(zs+5)*10,edgecolors='none',cmap=pl.cm.spectral) 499 >>> p = pl.xlim(pl.xlim()[::-1]);p = pl.ylim(pl.ylim()[::-1]) 500 >>> p = pl.xlabel('Teff');p = pl.ylabel('Logg') 501 502 >>> p = pl.subplot(rows,cols,3+cols) 503 >>> p = pl.scatter(ebvs,zs,c=teffs,s=(loggs+2)*10,edgecolors='none',cmap=pl.cm.spectral) 504 >>> p = pl.xlabel('E(B-V)');p = pl.ylabel('Z') 505 506 Confined to a small area in the grid's range with some parameters fixed: 507 508 >>> teffs,loggs,ebvs,zs = generate_grid(photbands,teffrange=(8765,8765),loggrange=(4.1,4.2),zrange=(0,0),ebvrange=(1,2),points=10000) 509 510 >>> p = pl.subplot(rows,cols,4) 511 >>> p = pl.scatter(teffs,loggs,c=ebvs,s=(zs+5)*10,edgecolors='none',cmap=pl.cm.spectral) 512 >>> p = pl.xlim(pl.xlim()[::-1]);p = pl.ylim(pl.ylim()[::-1]) 513 >>> p = pl.xlabel('Teff');p = pl.ylabel('Logg') 514 515 >>> p = pl.subplot(rows,cols,4+cols) 516 >>> p = pl.scatter(ebvs,zs,c=teffs,s=(loggs+2)*10,edgecolors='none',cmap=pl.cm.spectral) 517 >>> p = pl.xlabel('E(B-V)');p = pl.ylabel('Z') 518 519 ]include figure]]ivs_sed_fit_grids.png] 520 521 @param photbands: a list of photometric passbands, corresponding each 522 measurement 523 @type photbands: list of strings 524 @param teffrange: range of temperatures to use 525 @type teffrange: 2-tuple 526 @param loggrange: range of surface gravities to use 527 @type loggrange: 2-tuple 528 @param ebvrange: range of reddenings to use 529 @type ebvrange: 2-tuple 530 @param points: points to sample (when None, predefined grid points are used) 531 @type points: int 532 @param res: resolution of the original grid (the higher, the coarser) 533 @type res: int 534 @keyword clear_memory: flag to clear memory from previously loaded SED tables. 535 If you set it to False, you can easily get an overloaded memory! 536 @type clear_memory: boolean 537 @return: record array containing the searched grid, chi-squares and scale 538 factors 539 @rtype: record array 540 """ 541 #test 542 logger.info('Grid search with parameters teffrange=%s, loggrange=%s, ebvrange=%s, zrange=%s, points=%s'%(teffrange,loggrange,ebvrange,zrange,points)) 543 544 #-- we first get/set the grid. Calling this function means it will be 545 # memoized, so that we can safely thread (and don't have to memoize for 546 # each thread). We also have an exact view of the size of the grid here... 547 markers,(unique_teffs,unique_loggs,unique_ebvs,unique_zs),gridpnts,flux = \ 548 model._get_itable_markers(photbands,ebvrange=(-np.inf,np.inf), 549 zrange=(-np.inf,np.inf),include_Labs=True, 550 clear_memory=clear_memory,**kwargs) 551 552 if not kwargs: 553 logger.info('Received grid (%s)'%model.defaults2str()) 554 else: 555 logger.info('Received custom grid (%s)'%kwargs) 556 teffs,loggs,ebvs,zs = gridpnts.T 557 558 #-- We need to avoid having only one grid point! If nessessary the grid needs to be 559 # broader to get points in the entire interval. 560 index1 = teffrange[0] in unique_teffs and unique_teffs.searchsorted(teffrange[0]) or \ 561 max(0,unique_teffs.searchsorted(teffrange[0])-1) 562 index2 = teffrange[1] in unique_teffs and unique_teffs.searchsorted(teffrange[1]) or \ 563 min(len(unique_teffs),unique_teffs.searchsorted(teffrange[1])) 564 unique_teffs = unique_teffs[index1:index2+1] 565 index1 = max(0,unique_teffs.searchsorted(loggrange[0])-1) 566 index2 = unique_teffs.searchsorted(loggrange[1])+1 567 unique_loggs = unique_loggs[index1:index2+1] 568 569 #-- if we gave a certain number of points, we need to choose our points 570 # randomly in the grid: the grid is usually not a square, so we have to 571 # subdivide it 572 if points: 573 #-- generate appropriate evaluation points uniformly within the predefined 574 # edges 575 #-- first list all effective temperatures and their minimum logg in the grid 576 teff_min_logg = np.zeros((len(unique_teffs),2)) 577 for i,iteff in enumerate(unique_teffs): 578 teff_min_logg[i] = iteff,(loggs[teffs==iteff]).min() 579 #-- we have one square per logg: calculate their sizes 580 unique_min_loggs = sorted(list(set(teff_min_logg[:,1]))) 581 limits_and_sizes = [] 582 for index,unique_min_logg in enumerate(unique_min_loggs): 583 min_teff = teff_min_logg[:,0][teff_min_logg[:,1]==unique_min_logg].min() 584 #-- we need to avoid having gaps in the grid: 585 if index>0: 586 min_teff = max_teff 587 max_teff = teff_min_logg[:,0][teff_min_logg[:,1]==unique_min_logg].max() 588 min_logg = unique_min_logg 589 max_logg = loggs.max() 590 #-- we're at too low temperatures 591 if max_teff<teffrange[0]: continue 592 else: 593 min_teff = max(teffrange[0],min_teff) 594 max_teff = min(teffrange[1],max_teff) 595 #-- we're at too low surface gravities: 596 min_logg = max(loggrange[0],min_logg) 597 max_logg = min(loggrange[1],max_logg) 598 #-- make sure there are points defined even if some range in parameters 599 # equals zero 600 if (max_teff-min_teff)>1 and (max_logg-min_logg)>0.01: 601 size = (max_teff-min_teff)*(max_logg-min_logg) 602 elif (max_teff-min_teff)>1: 603 size = (max_teff-min_teff) 604 elif (max_logg-min_logg)>1: 605 size = (max_logg-min_logg) 606 else: 607 size = int(float(points)/(len(unique_min_loggs))) 608 if size==0: size=2 609 #-- sizes of ebv and z: 610 zrange_ = max( zrange[0],min(unique_zs)), min( zrange[1],max(unique_zs)) 611 ebvrange_ = max(ebvrange[0],min(unique_ebvs)), min(ebvrange[1],max(unique_ebvs)) 612 limits_and_sizes.append([(min_teff,max_teff),(min_logg,max_logg),ebvrange_,zrange_,size]) 613 total_size = sum([row[-1] for row in limits_and_sizes]) 614 #-- in the following case, we fall in between the grid points. We correct 615 # for this 616 if len(limits_and_sizes)==0: 617 total_size = points 618 limits_and_sizes = [[teffrange,loggrange,ebvrange,zrange,points]] 619 logger.debug('Limits and sizes of boxes:'+str(limits_and_sizes)) 620 teffs,loggs,ebvs,zs = np.hstack([np.random.uniform(low=[lims[0][0],lims[1][0],lims[2][0],lims[3][0]], 621 high=[lims[0][1],lims[1][1],lims[2][1],lims[3][1]], 622 size=(int(lims[-1]/total_size*points),4)).T for lims in limits_and_sizes]) 623 keep = (teffrange[0]<=teffs) & (teffs<=teffrange[1]) &\ 624 (loggrange[0]<=loggs) & (loggs<=loggrange[1]) &\ 625 (ebvrange[0]<=ebvs) & (ebvs<=ebvrange[1]) &\ 626 (zrange[0]<=zs) & (zs<=zrange[1]) 627 teffs,loggs,ebvs,zs = teffs[keep],loggs[keep],ebvs[keep],zs[keep] 628 629 if res: 630 teffs,loggs,ebvs,zs = teffs[::res],loggs[::res],ebvs[::res],zs[::res] 631 logger.info('Evaluating %d points in parameter space'%(len(teffs))) 632 633 return teffs,loggs,ebvs,zs

634

635 -def generate_grid(photbands,teffrange=((-inf,inf),(-inf,inf)), 636 loggrange=((-inf,inf),(-inf,inf)),ebvrange=(-inf,inf), 637 zrange=((-inf,inf),(-inf,inf)), 638 radiusrange=((1,1),(0.1,10.)),grids=None, 639 points=None,res=None,clear_memory=False, 640 type='single',primary_hottest=False, **kwargs):

641 """ 642 Generate grid points at which to fit an interpolated grid of SEDs. 643 644 If C{points=None}, the points are chosen on the predefined grid points. 645 Otherwise, C{points} grid points will be generated, uniformly distributed 646 between the ranges defined by C{teffrange}, C{loggrange} and C{ebvrange}. If 647 you set the resolution to C{2}, one out of every two points will be selected. 648 649 Setting C{primary_hottest} makes sure that, when doing a binary fit, that 650 the the effective temperature of the first star is hotter than the second. 651 652 Extra keyword arguments can be used to give more details on the atmosphere 653 models to use. 654 655 Colors are automatically detected. 656 657 You can fix one parameter e.g. via setting teffrange=(10000,10000). 658 659 >>> photbands = ['GENEVA.G','GENEVA.B-V'] 660 661 Start a figure: 662 663 >>> p = pl.figure() 664 >>> rows,cols = 2,4 665 666 On the grid points, but only one in every 100 points (otherwise we have over 667 a million points): 668 669 >>> teffs,loggs,ebvs,zs = generate_grid(photbands,res=100) 670 671 >>> p = pl.subplot(rows,cols,1) 672 >>> p = pl.scatter(teffs,loggs,c=ebvs,s=(zs+5)*10,edgecolors='none',cmap=pl.cm.spectral) 673 >>> p = pl.xlim(pl.xlim()[::-1]);p = pl.ylim(pl.ylim()[::-1]) 674 >>> p = pl.xlabel('Teff');p = pl.ylabel('Logg') 675 676 >>> p = pl.subplot(rows,cols,1+cols) 677 >>> p = pl.scatter(ebvs,zs,c=teffs,s=(loggs+2)*10,edgecolors='none',cmap=pl.cm.spectral) 678 >>> p = pl.xlabel('E(B-V)');p = pl.ylabel('Z') 679 680 Randomly distributed over the grid's ranges: 681 682 >>> teffs,loggs,ebvs,zs = generate_grid(photbands,points=10000) 683 684 >>> p = pl.subplot(rows,cols,2) 685 >>> p = pl.scatter(teffs,loggs,c=ebvs,s=(zs+5)*10,edgecolors='none',cmap=pl.cm.spectral) 686 >>> p = pl.xlim(pl.xlim()[::-1]);p = pl.ylim(pl.ylim()[::-1]) 687 >>> p = pl.xlabel('Teff');p = pl.ylabel('Logg') 688 689 >>> p = pl.subplot(rows,cols,2+cols) 690 >>> p = pl.scatter(ebvs,zs,c=teffs,s=(loggs+2)*10,edgecolors='none',cmap=pl.cm.spectral) 691 >>> p = pl.xlabel('E(B-V)');p = pl.ylabel('Z') 692 693 Confined to a small area in the grid's range: 694 695 >>> teffs,loggs,ebvs,zs = generate_grid(photbands,teffrange=(8000,10000),loggrange=(4.1,4.2),zrange=(0,inf),ebvrange=(1.,2),points=10000) 696 697 >>> p = pl.subplot(rows,cols,3) 698 >>> p = pl.scatter(teffs,loggs,c=ebvs,s=(zs+5)*10,edgecolors='none',cmap=pl.cm.spectral) 699 >>> p = pl.xlim(pl.xlim()[::-1]);p = pl.ylim(pl.ylim()[::-1]) 700 >>> p = pl.xlabel('Teff');p = pl.ylabel('Logg') 701 702 >>> p = pl.subplot(rows,cols,3+cols) 703 >>> p = pl.scatter(ebvs,zs,c=teffs,s=(loggs+2)*10,edgecolors='none',cmap=pl.cm.spectral) 704 >>> p = pl.xlabel('E(B-V)');p = pl.ylabel('Z') 705 706 Confined to a small area in the grid's range with some parameters fixed: 707 708 >>> teffs,loggs,ebvs,zs = generate_grid(photbands,teffrange=(8765,8765),loggrange=(4.1,4.2),zrange=(0,0),ebvrange=(1,2),points=10000) 709 710 >>> p = pl.subplot(rows,cols,4) 711 >>> p = pl.scatter(teffs,loggs,c=ebvs,s=(zs+5)*10,edgecolors='none',cmap=pl.cm.spectral) 712 >>> p = pl.xlim(pl.xlim()[::-1]);p = pl.ylim(pl.ylim()[::-1]) 713 >>> p = pl.xlabel('Teff');p = pl.ylabel('Logg') 714 715 >>> p = pl.subplot(rows,cols,4+cols) 716 >>> p = pl.scatter(ebvs,zs,c=teffs,s=(loggs+2)*10,edgecolors='none',cmap=pl.cm.spectral) 717 >>> p = pl.xlabel('E(B-V)');p = pl.ylabel('Z') 718 719 ]include figure]]ivs_sed_fit_grids.png] 720 721 @param photbands: a list of photometric passbands, corresponding each 722 measurement 723 @type photbands: list of strings 724 @param teffrange: range of temperatures to use 725 @type teffrange: 2-tuple 726 @param loggrange: range of surface gravities to use 727 @type loggrange: 2-tuple 728 @param ebvrange: range of reddenings to use 729 @type ebvrange: 2-tuple 730 @param points: points to sample (when None, predefined grid points are used) 731 @type points: int 732 @param res: resolution of the original grid (the higher, the coarser) 733 @type res: int 734 @keyword clear_memory: flag to clear memory from previously loaded SED tables. 735 If you set it to False, you can easily get an overloaded memory! 736 @type clear_memory: boolean 737 @return: record array containing the searched grid, chi-squares and scale 738 factors 739 @rtype: record array 740 """ 741 742 #-- Select the grid 743 # but remove metallicity, as it will be fitted! 744 if type=='single': 745 #--Single grid, uses the basic function 746 teffs,loggs,ebvs,zs = generate_grid_single(photbands,teffrange=teffrange, 747 loggrange=loggrange,ebvrange=ebvrange, 748 zrange=zrange,points=points,clear_memory=clear_memory) 749 radii = np.ones(len(teffs))#[1 for i in teffs] 750 return teffs,loggs,ebvs,zs,radii 751 752 #-- first collect the effetive temperatures, loggs, ebvs, zs for the 753 # different stars in the multiple system 754 pars = [] 755 if grids is None: 756 grids = model.defaults_multiple 757 #-- but remove metallicity, as it will be fitted! 758 for grid in grids: 759 if 'z' in grid: 760 thrash = grid.pop('z') 761 for i,grid in enumerate(grids): 762 #-- it is possible that we want certain parameters to be the same for 763 # all components 764 teffrange_ = hasattr(teffrange[0],'__iter__') and teffrange[i] or teffrange 765 loggrange_ = hasattr(loggrange[0],'__iter__') and loggrange[i] or loggrange 766 ebvrange_ = hasattr(ebvrange[0],'__iter__') and ebvrange[i] or ebvrange 767 zrange_ = hasattr(zrange[0],'__iter__') and zrange[i] or zrange 768 pars += list(generate_grid_single(photbands,teffrange=teffrange_, 769 loggrange=loggrange_,ebvrange=ebvrange_, 770 zrange=zrange_,points=points,**grid)) 771 #-- the L{generate_grid} method does not guarantee the number of points. 772 # We have to strip some points if the arrays don't have the same shape 773 nmin = np.min([len(i) for i in pars]) 774 pars = [i[:nmin] for i in pars] 775 pars = np.array(pars) 776 #-- permute parameters so that the different blocks from the generate_grid 777 # are not clustered together 778 for i in range(0,len(pars),4): 779 permutation = np.random.permutation(len(pars[0])) 780 pars[i:i+4] = pars[i:i+4,permutation] 781 #-- make arrays of the output parameters 782 teffs,loggs,ebvs,zs = pars[0::4].T,pars[1::4].T,pars[2::4].T,pars[3::4].T 783 784 #-- keep in mind that we probably want all the members in the system to have 785 # the same value for the interstellar reddening and metallicity, though 786 # this is not obligatory 787 if not hasattr(teffrange[0],'__iter__'): teffs = np.column_stack([teffs[:,0]]*len(grids)) 788 if not hasattr(loggrange[0],'__iter__'): loggs = np.column_stack([loggs[:,0]]*len(grids)) 789 #if not hasattr(ebvrange[0],'__iter__'): ebvs = np.column_stack([ebvs[:,0]]*len(grids)) 790 #if not hasattr(zrange[0],'__iter__'): zs = np.column_stack([zs[:,0]]*len(grids)) 791 ebvs = np.column_stack([ebvs[:,0]]*len(grids)) 792 zs = np.column_stack([zs[:,0]]*len(grids)) 793 794 if type=='binary': 795 #-- The radius of the stars is calculated based on logg and the provided masses 796 masses = 'masses' in kwargs and kwargs['masses'] or (1,1) 797 G = constants.GG_cgs 798 Msol = constants.Msol_cgs 799 Rsol = constants.Rsol_cgs 800 radius1 = np.sqrt(G*masses[0]*Msol/10**loggs[:,0])/Rsol 801 radius2 = np.sqrt(G*masses[1]*Msol/10**loggs[:,1])/Rsol 802 #radii = radius2/radius1 803 804 #radii = np.array([np.ones(len(radii)),radii]).T 805 radii = np.array([radius1,radius2]).T 806 807 #-- maybe we need to lower the temperatures of the secondary, so that 808 # it is not hotter than the primary effective temperature 809 if primary_hottest: 810 wrong = teffs[:,1]>teffs[:,0] 811 teffs[wrong,1] = np.random.uniform(low=teffs[:,1].min()*np.ones(sum(wrong)),\ 812 high=teffs[wrong,0]) 813 logger.info('Ensured primary is the hottest (%d/%d corrections)'%(sum(wrong),len(wrong))) 814 815 816 elif type=='multiple': 817 #-- We have random different radii for the stars 818 radii = np.array([10**np.random.uniform(low=radiusrange[0][0], high=radiusrange[0][1], size=(len(teffs))), 819 10**np.random.uniform(low=radiusrange[1][0], high=radiusrange[1][1], size=(len(teffs)))]).T 820 #radii = 10**np.random.uniform(low=[np.log10(i[0]) for i in radiusrange], 821 #high=[np.log10(i[1]) for i in radiusrange],size=(len(teffs),2)) 822 823 return teffs,loggs,ebvs,zs,radii

824

825 #{ Fitting: grid search 826 827 -def igrid_search_pix(meas,e_meas,photbands, constraints={},**kwargs):

828 """ 829 Run over gridpoints and evaluate model C{model_func} via C{stat_func}. 830 831 The measurements are defined via C{meas, e_meas, photbands, colors} and 832 should be 1d arrays of equal length. C{colors} should be a boolean 833 array, C{photbands} should be a string array. 834 835 The grid points are defined via C{args}. C{args} should be a tuple of 836 1x? dimensional arrays of equal length. For single stars, this is 837 typically effective temperatures, loggs, reddenings and metallicities. 838 For multiple systems, (at least some of the) previously mentioned 839 parameters are typically doubled, and radius ratios are added. Remember 840 to specify the C{model_func} to match single or multiple systems. 841 842 At each grid point, the pre-calculated photometry will be retrieved via 843 the keyword C{model_func} and compared to the measurements via the function 844 definded via C{stat_func}. This function should be of the same form as 845 L{stat_chi2}. 846 847 Extra arguments are passed to L{parallel_gridsearch} for parallelization 848 and to {model_func} for further specification of grids etc. 849 850 The index array is returned to trace the results after parallelization. 851 852 @param meas: the measurements that have to be compared with the models 853 @type meas: 1D numpy array of floats 854 @param e_meas: errors on the measurements 855 @type e_meas: 1D numpy array of floats 856 @param photbands: names of the photometric passbands 857 @type photbands: 1D numpy array of strings 858 @keyword model_func: function to translate parameters to synthetic (model) data 859 @type model_func: function 860 @keyword stat_func: function to evaluate the fit 861 @type stat_func: function 862 @return: (chi squares, scale factors, error on scale factors, absolute 863 luminosities (R=1Rsol) 864 @rtype: array 865 """ 866 model_func = kwargs.pop('model_func',model.get_itable_pix) 867 stat_func = kwargs.pop('stat_func',stat_chi2) 868 colors = np.array([filters.is_color(photband) for photband in photbands],bool) 869 #-- run over the grid, retrieve synthetic fluces and compare with 870 # observations. 871 syn_flux,lumis = model_func(photbands=photbands,**kwargs) 872 chisqs,scales,e_scales = stat_func(meas.reshape(-1,1),\ 873 e_meas.reshape(-1,1),\ 874 colors,syn_flux, **constraints) 875 #-- return results 876 return chisqs,scales,e_scales,lumis

877

878 @parallel_gridsearch 879 @make_parallel 880 -def igrid_search(meas,e_meas,photbands,*args,**kwargs):

881 """ 882 Run over gridpoints and evaluate model C{model_func} via C{stat_func}. 883 884 The measurements are defined via C{meas, e_meas, photbands, colors} and 885 should be 1d arrays of equal length. C{colors} should be a boolean 886 array, C{photbands} should be a string array. 887 888 The grid points are defined via C{args}. C{args} should be a tuple of 889 1x? dimensional arrays of equal length. For single stars, this is 890 typically effective temperatures, loggs, reddenings and metallicities. 891 For multiple systems, (at least some of the) previously mentioned 892 parameters are typically doubled, and radius ratios are added. Remember 893 to specify the C{model_func} to match single or multiple systems. 894 895 At each grid point, the pre-calculated photometry will be retrieved via 896 the keyword C{model_func} and compared to the measurements via the function 897 definded via C{stat_func}. This function should be of the same form as 898 L{stat_chi2}. 899 900 Extra arguments are passed to L{parallel_gridsearch} for parallelization 901 and to {model_func} for further specification of grids etc. 902 903 The index array is returned to trace the results after parallelization. 904 905 @param meas: the measurements that have to be compared with the models 906 @type meas: 1D numpy array of floats 907 @param e_meas: errors on the measurements 908 @type e_meas: 1D numpy array of floats 909 @param photbands: names of the photometric passbands 910 @type photbands: 1D numpy array of strings 911 @keyword model_func: function to translate parameters to synthetic (model) data 912 @type model_func: function 913 @keyword stat_func: function to evaluate the fit 914 @type stat_func: function 915 @return: (chi squares, scale factors, error on scale factors, absolute 916 luminosities (R=1Rsol), index 917 @rtype: 4/5X1d array 918 """ 919 model_func = kwargs.pop('model_func',model.get_itable) 920 stat_func = kwargs.pop('stat_func',stat_chi2) 921 index = kwargs.pop('index',None) 922 fitkws = {} 923 if 'distance' in kwargs and kwargs['distance'] != None: 924 fitkws = {'distance':kwargs['distance']} 925 N = len(args[0]) 926 #-- prepare output arrays 927 chisqs = np.zeros(N) 928 scales = np.zeros(N) 929 e_scales = np.zeros(N) 930 lumis = np.zeros(N) 931 colors = np.array([filters.is_color(photband) for photband in photbands],bool) 932 #-- show a progressMeter when not parallelized 933 if index is None: 934 p = progressMeter.ProgressMeter(total=N) 935 #-- run over the grid, retrieve synthetic fluces and compare with 936 # observations. 937 for n,pars in enumerate(itertools.izip(*args)): 938 if index is None: p.update(1) 939 syn_flux,Labs = model_func(*pars,photbands=photbands,**kwargs) 940 chisqs[n],scales[n],e_scales[n] = stat_func(meas,e_meas,colors,syn_flux, **fitkws) 941 lumis[n] = Labs 942 #-- return results 943 if index is not None: 944 return chisqs,scales,e_scales,lumis,index 945 else: 946 return chisqs,scales,e_scales,lumis

947

948 #} 949 950 #{ Fitting: minimizer 951 952 -def create_parameter_dict(**pars):

953 954 #-- Find all the parameters first 955 parnames = set() 956 atributes = set() 957 for key in pars.keys(): 958 name, att = re.findall("(.*)_([a-zA-Z]+)$", key)[0] 959 parnames.add(name) 960 atributes.add(att) 961 962 #-- create dictionary with the attributes 963 parnames = np.array(list(parnames)) 964 result = dict(names=parnames) 965 for att in atributes: 966 result[att] = np.array([None for i in parnames]) 967 968 #-- read the attributes 969 for key in pars.keys(): 970 name, att = re.findall("(.*)_([a-zA-Z]+)$", key)[0] 971 result[att][parnames == name] = pars[key] 972 973 return result

974

975 -def _get_info_from_minimizer(minimizers, startpars, photbands, meas, e_meas, **fitkws):

976 scales, lumis, chisqrs, nfevs, allpars = [], [], [], [], {} 977 for n in fitkws['pnames']: 978 allpars[n] = np.array([]) 979 allpars[n+'start'] = np.array([]) 980 for mini, start in zip(minimizers, startpars): 981 chisqrs.append(mini.chisqr) 982 nfevs.append(mini.nfev) 983 984 val, err = mini.model.get_parameters(full_output=False) 985 for n, v in zip(fitkws['pnames'], val): 986 allpars[n] = np.append(allpars[n], [v]) 987 allpars[n+'start'] = np.append(allpars[n+'start'], [start[n].value]) 988 989 synth, lum = mini.model.evaluate(photbands, **fitkws) 990 distance = fitkws['distance'] if 'distance' in fitkws else None 991 if distance != None: 992 scale = 1/distance**2 993 else: 994 ratio = (meas/synth[:,0]) 995 weights = (meas/e_meas) 996 scale = np.average(ratio,weights=weights) 997 lumis.append(lum[0]) 998 scales.append(scale) 999 1000 return np.array(chisqrs), np.array(nfevs), np.array(scales), np.array(lumis), allpars

1001

1002 -def _iminimize_model(varlist, x, *args, **kws):

1003 pnames = kws.pop('pnames') 1004 pars = {} 1005 for n, v in zip(pnames, varlist): 1006 pars[n] = np.array([v]) 1007 pars.update(kws) 1008 #print pars 1009 return model.get_itable_pix(wave_units=None, photbands=x, **pars)

1010

1011 -def _iminimize_residuals(synth, meas, weights=None, **kwargs):

1012 synth = synth[0][:,0] #select the flux. 1013 e_meas = 1 / weights 1014 if 'distance' in kwargs: 1015 scale = 1/kwargs['distance']**2 1016 else: 1017 ratio = (meas/synth) 1018 weights = (meas/e_meas) 1019 scale = np.average(ratio,weights=weights) 1020 return (meas - synth*scale)/e_meas

1021

1022 -def iminimize(meas,e_meas,photbands, points=None, return_minimizer=False,**kwargs):

1023 """ 1024 minimizer based on the sigproc.fit lmfit minimizer. 1025 provide the observed data, the fitting model, residual function and parameter 1026 information about the variables, and this function will return the best fit 1027 parameters together with extra information about the fit. 1028 1029 if the fitkws keyword is supplied, this dict will be made available to the 1030 model_func (fit model) during the fitting process. The order of the parameters 1031 will also be made available as the 'pnames' keyword. 1032 """ 1033 1034 kick_list = kwargs.pop('kick_list', None) 1035 constraints = kwargs.pop('constraints', dict()) 1036 fitmodel = kwargs.pop('model_func',_iminimize_model) 1037 residuals = kwargs.pop('res_func',_iminimize_residuals) 1038 epsfcn = kwargs.pop('epsfcn', 0.0005)# using ~3% step to derive jacobian. 1039 1040 #-- get the parameters 1041 parameters = create_parameter_dict(**kwargs) 1042 1043 #-- setup the fitting model 1044 pnames = parameters.pop('names') 1045 fmodel = sfit.Function(function=fitmodel, par_names=pnames) 1046 fmodel.setup_parameters(**parameters) 1047 1048 #-- fit the model to the data 1049 fitkws = dict(pnames=pnames) 1050 fitkws.update(constraints) 1051 if points == None: 1052 startpars = [copy.deepcopy(fmodel.parameters)] 1053 minimizer = sfit.minimize(photbands,meas, fmodel, weights=1/e_meas, kws=fitkws, \ 1054 resfunc=residuals, engine='leastsq', epsfcn=epsfcn,\ 1055 ftol=0.001, xtol=0.001) 1056 minimizer = [minimizer] 1057 else: 1058 minimizer, startpars, newmodels, chisqr = sfit.grid_minimize(photbands, meas, fmodel, \ 1059 weights=1/e_meas, kws=fitkws, resfunc=residuals, engine='leastsq', \ 1060 epsfcn=epsfcn, points=points, parameters=kick_list, return_all=True) 1061 1062 if return_minimizer: 1063 #-- return the actual minimizer used by the calculate ci methods 1064 return minimizer[0] 1065 else: 1066 #-- for all other users return the actual results 1067 chisqr, nfev, scale, lumis, grid = _get_info_from_minimizer(minimizer, startpars,\ 1068 photbands, meas, e_meas, **fitkws) 1069 ##-- collect all parameter info 1070 #val, err, vary, min, max, expr = fmodel.get_parameters(full_output=True) 1071 #pars = dict(name=pnames, value=val, cilow=min, cihigh=max) 1072 1073 return grid, chisqr, nfev, scale, lumis

1074

1075 -def calculate_iminimize_CI(meas,e_meas,photbands, CI_limit=0.66, **kwargs):

1076 """ 1077 Calculate the confidence intervalls for every parameter. 1078 When ci fails, the boundaries are returned as ci. 1079 """ 1080 1081 maxiter = kwargs.pop('maxiter', 100) 1082 1083 mini = iminimize(meas,e_meas,photbands, return_minimizer=True, **kwargs) 1084 val, err, vary, low, high, expr = mini.model.get_parameters(full_output=True) 1085 pnames = mini.model.par_names 1086 1087 #-- setup stat function 1088 def prob_func(Ndata, Nparas, new_chi, best_chi, Nfix=1.): 1089 k = Ndata-Nparas 1090 factor = max(best_chi/k,1) 1091 return scipy.stats.distributions.chi2.cdf(new_chi/factor,k)

1092 1093 cilow, cihigh = low.copy(), high.copy() 1094 for i,p in enumerate(pnames): 1095 try: 1096 ci = mini.calculate_CI(parameters=[p], short_output=True,\ 1097 sigma=CI_limit, maxiter=maxiter, prob_func=prob_func) 1098 logger.debug('Calculated ci for parameter %s: %s'%(p,ci) ) 1099 if ci[0] != None: cilow[i] = ci[0] 1100 if ci[1] != None: cihigh[i] = ci[1] 1101 except Exception: 1102 logger.warning('Failed to calculate CI for parameter: %s'%(p)) 1103 1104 return dict(name=pnames, value=val, cilow=cilow, cihigh=cihigh) 1105

1106 -def calculate_iminimize_CI2D(meas,e_meas,photbands, xpar, ypar, df=None, **kwargs):

1107 """ 1108 Calculate a 2D confidence grid for the given parameters. 1109 You can provide limits on the parameters yourself, if none are 1110 provided, the function will take care of it, but it might crash 1111 if the limits are outside the model range. 1112 returns: x vals, y vals, ci values 1113 """ 1114 maxiter = kwargs.pop('maxiter', 25) 1115 res = kwargs.pop('res', (10,10)) 1116 if type(res) == int: res = (res,res) 1117 limits = kwargs.pop('limits', None) 1118 citype = kwargs.pop('type', 'prob') 1119 1120 mini = iminimize(meas,e_meas,photbands, return_minimizer=True, **kwargs) 1121 val, err, vary, low, high, expr = mini.model.get_parameters(full_output=True) 1122 val, low, high = np.array(val, dtype=float), np.array(low, dtype=float), np.array(high, dtype=float) 1123 1124 #-- set some keywords so the process speeds up a bit: 1125 mini.userkws.update({'epsfcn':0.001, 'xtol':0.005, 'ftol':0.005}) 1126 1127 1128 name = mini.model.par_names 1129 xval, yval = val[name==xpar], val[name==ypar] 1130 1131 logger.debug('Starting calculations of CI2D with parameters:\n%s'%\ 1132 (mini.model.param2str(full_output=True))) 1133 1134 #-- if not provided choose own limits 1135 if limits != None: 1136 xmin, xmax = limits[0] 1137 ymin, ymax = limits[1] 1138 else: 1139 xmin = low[name==xpar] if low[name==xpar] != None else xval - 5*err[name==xpar] 1140 xmax = high[name==xpar] if high[name==xpar] != None else xval + 5*err[name==xpar] 1141 ymin = low[name==ypar] if low[name==ypar] != None else yval - 5*err[name==ypar] 1142 ymax = high[name==ypar] if high[name==ypar] != None else yval + 5*err[name==ypar] 1143 1144 #-- fix limits to include best fit in CI 1145 xstep = (xmax - xmin) / float(res[0]) 1146 ystep = (ymax - ymin) / float(res[1]) 1147 xmin = xval - np.floor((xval - xmin) / xstep) * xstep 1148 xmax = xval + (res[0] - 1 - np.floor((xval - xmin) / xstep)) * xstep 1149 ymin = yval - np.floor((yval - ymin) / ystep) * ystep 1150 ymax = yval + (res[1] - 1 - np.floor((yval - ymin) / ystep)) * ystep 1151 limits = [(xmin, xmax), (ymin, ymax)] 1152 1153 logger.debug('Recalculated grid limits to: \n\t%s: %s <-> %s \n\t%s %s <-> %s'%(xpar, 1154 limits[0][0], limits[0][1], ypar, limits[1][0], limits[1][1])) 1155 1156 x,y,ci_chi2 = mini.calculate_CI_2D(xpar=xpar, ypar=ypar, res=res, limits=limits, 1157 ctype='chi2') 1158 1159 #-- do the statistics 1160 N = mini.ndata 1161 if df == None: df = mini.nvarys 1162 k = N-df 1163 factor = max(mini.chisqr/k,1) 1164 ci_raw = scipy.stats.distributions.chi2.cdf(ci_chi2,k) 1165 ci_red = scipy.stats.distributions.chi2.cdf(ci_chi2/factor,k) 1166 1167 return x, y, ci_chi2, ci_raw, ci_red

1168

1169 -def iminimize2(meas,e_meas,photbands,*args,**kwargs):

1170 model_func = kwargs.pop('model_func',model.get_itable) 1171 #res_func = kwargs.pop('res_func',residual_single) # not yet defined! 1172 method = kwargs.pop('fitmethod','fmin') 1173 stat_func = kwargs.pop('stat_func',stat_chi2) 1174 1175 colors = np.array([filters.is_color(photband) for photband in photbands],bool) 1176 1177 # the help function which returns the chisquare NOTE: metallicity is not yet included in the fitting! 1178 def residual_single(parameters): 1179 syn_flux,Labs = model_func(*parameters,photbands=photbands,**kwargs) 1180 # in case any of the parameters goes out of its bounds 1181 #if isnan(syn_flux).any(): 1182 chisq,scale,e_scale = stat_func(meas,e_meas,colors,syn_flux,full_output=False) 1183 return chisq

1184 res_func = residual_single 1185 # calling the fitting function NOTE: the initial metallicity is returned in the output! 1186 1187 if method=='fmin': # fmin 1188 optpars,fopt,niter,funcalls,warnflag = fmin(res_func,np.array(args),xtol=0.0001,disp=0,full_output=True) 1189 elif method=='fmin_powell': #fmin_powell 1190 optpars = fmin_powell(res_func,np.array(args)) 1191 else: 1192 raise NotImplementedError 1193 logger.debug("Optimization finished") 1194 syn_flux,Labs = model_func(*optpars,photbands=photbands,**kwargs) 1195 1196 # when any of the parameters goes out of the bounds of the grid, syn_flux contains NaN 1197 if np.isnan(syn_flux).any(): 1198 warnflag = 3 1199 1200 stats = stat_func(meas,e_meas,colors,syn_flux,full_output=False) 1201 optpars = np.hstack([optpars,stats]) 1202 # stats: chisq, scale, e_scale 1203 return optpars,warnflag 1204 #} 1205 1206 1207 if __name__=="__main__": 1208 from ivs.aux import loggers 1209 import time 1210 import pylab as plt 1211 logger = loggers.get_basic_logger(clevel='DEBUG') 1212 1213 photbands = ['GENEVA.G','GENEVA.B-V'] 1214 1215 c0 = time.time() 1216 teffs,loggs,ebvs,zs,radii = generate_grid(photbands,teffrange=(5000,5800),loggrange=(4.20,4.70),zrange=(0,0),ebvrange=(0.05,0.08), grid='kurucz',points=10000) 1217 print 'Time: %i'%(time.time()-c0) 1218 1219 plt.figure(2) 1220 plt.scatter(teffs,loggs,c=ebvs,s=(zs+5)*10,edgecolors='none',cmap=plt.cm.spectral) 1221 plt.xlim(plt.xlim()[::-1]) 1222 plt.ylim(plt.ylim()[::-1]) 1223 plt.xlabel('Teff') 1224 plt.ylabel('Logg') 1225 plt.show() 1226 1227 sys.exit() 1228 1229 1230 import doctest 1231 import pylab as pl 1232 doctest.testmod() 1233 pl.show() 1234 1235 1236 sys.exit() 1237 from ivs.aux import loggers 1238 from pylab import * 1239 from numpy import * 1240 logger = loggers.get_basic_logger() 1241 random.seed(1111) 1242 A,grid = get_PCA_grid(['GENEVA.U-B','GENEVA.B1-B','GENEVA.B2-B','GENEVA.V-B','GENEVA.V1-B','GENEVA.G-B','2MASS.J-H','2MASS.KS-H'],ebvrange=(0,0.5),res=10) 1243 P,T,(means,stds) = get_PCA(A) 1244 calib = calibrate_PCA(T,grid,function='linear') 1245 sample_index = int(np.random.uniform(high=len(A))) 1246 sample = 10**A[sample_index] 1247 print [bla[sample_index] for bla in grid] 1248 print get_PCA_parameters(sample,calib,P,means,stds) 1249

Source Code for Module ivs.sed.fit