ComboCode.cc.ivs.sigproc.evaluate

1 """ 2 Evaluate model functions. 3 4 Most functions accept a domain of interest and a set of parameters, and return 5 a model function (e.g., a sine, gaussian etc...). Parameters are usually a 6 record array containing named fields for the parameters. 7 8 Some nonlinear optimizers expect and return parameters as a 1D array of the 9 independent parameters. Therefore, all functions also accept these 1D numpy 10 arrays, and they are converted by L{check_input} to record arrays. The latter 11 function actually accepts both record arrays and 1D arrays, and acts as a 12 translator between the two parameter representations. 13 """ 14 import logging 15 import functools 16 17 import numpy as np 18 from numpy import pi,cos,sin,sqrt,tan,arctan 19 from scipy.interpolate import splev 20 from scipy.special import erf,jn 21 from scipy.stats import distributions 22 # from ivs.timeseries import keplerorbit 23 24 logger = logging.getLogger('SIGPROC.EVAL')

25 26 -def check_input(fctn):

27 """ 28 Decorator to check input of evaluating model, and transforming input from 29 record array to 1D numpy array for optimization, and vice versa. 30 """ 31 @functools.wraps(fctn) 32 def check(*args,**kwargs): 33 args = list(args) 34 parameters = np.asarray(args[1]) 35 if not parameters.dtype.names: 36 name = fctn.__name__ 37 for rem in ['_diffcorr']: 38 name = name.replace(rem,'') 39 parameters = globals()[name+'_preppars'](parameters) 40 args[1] = parameters 41 return fctn(*args,**kwargs)

42 return check 43

44 #{ Evaluating model fits 45 46 -def RSS(data,model):

47 """ 48 Residual sum of squares. 49 50 @param data: data 51 @type data: numpy array 52 @param model: model 53 @type model: numpy array 54 @rtype: float 55 @return: residual sum of squares 56 """ 57 return np.sum( (model-data)**2 )

58

59 60 -def varred(data,model):

61 """ 62 Calculate variance reduction. 63 64 @param data: data 65 @type data: numpy array 66 @param model: model 67 @type model: numpy array 68 @rtype: percentage 69 @return: variance reduction (percentage) 70 """ 71 residus = data-model 72 variance_reduction = 1.-np.sum(residus**2.)/sum((data-data.mean())**2.) 73 return variance_reduction*100.

74

75 -def bic(data,model,k=0,n=None):

76 """ 77 BIC for time regression 78 79 This one is based on the MLE of the variance, 80 81 Sigma_k^2 = RSS/n 82 83 Where our maximum likelihood estimator is then (RSS/n)^n 84 85 @param data: data 86 @type data: numpy array 87 @param model: model 88 @type model: numpy array 89 @param k: number of free, independent parameters 90 @type k: integer 91 @param n: number of effective observations 92 @type n: integer or None 93 @rtype: float 94 @return: BIC 95 """ 96 if n is None: 97 n = len(data) 98 #-- calculate RSS 99 RSS_ = RSS(data,model) 100 #-- calculate BIC 101 BIC = n * np.log(RSS_/n) + k * np.log(n) 102 return BIC

103

104 -def aic(data,model,k=0,n=None):

105 """ 106 AIC for time regression 107 108 This one is based on the MLE of the variance, 109 110 Sigma_k^2 = RSS/n 111 112 Where our maximum likelihood estimator is then (RSS/n)^n 113 114 @param data: data 115 @type data: numpy array 116 @param model: model 117 @type model: numpy array 118 @param k: number of free, independent parameters 119 @type k: integer 120 @param n: number of effective observations 121 @type n: integer or None 122 @rtype: float 123 @return: AIC 124 """ 125 if n is None: 126 n = len(data) 127 #-- calculate RSS 128 RSS_ = RSS(data,model) 129 #-- calculate AIC 130 AIC = n * np.log(RSS_/n) + 2*k + n 131 return AIC

132

133 -def fishers_method(pvalues):

134 """ 135 Use Fisher's test to combine probabilities. 136 137 http://en.wikipedia.org/wiki/Fisher's_method 138 139 >>> rvs1 = distributions.norm().rvs(10000) 140 >>> rvs2 = distributions.norm().rvs(10000) 141 >>> pvals1 = scipy.stats.distributions.norm.sf(rvs1) 142 >>> pvals2 = scipy.stats.distributions.norm.sf(rvs2) 143 >>> fpval = fishers_method([pvals1,pvals2]) 144 145 And make a plot (compare with the wiki page). 146 147 >>> p = pl.figure() 148 >>> p = pl.scatter(pvals1,pvals2,c=fpval,edgecolors='none',cmap=pl.cm.spectral) 149 >>> p = pl.colorbar() 150 151 152 153 154 @param pvalues: array of pvalues, the first axis will be combined 155 @type pvalues: ndarray (shape Nprob times n) 156 @return: combined pvalues 157 @rtype: ndarray (shape n) 158 """ 159 pvalues = np.asarray(pvalues) 160 X2 = -2*np.sum(np.log(pvalues),axis=0) 161 df = 2*pvalues.shape[0] 162 return distributions.chi2.sf(X2,df)

163

164 165 #} 166 167 #{ Timeseries 168 169 -def phasediagram(time,signal,nu0,D=0,t0=None,forb=None,asini=None, 170 return_indices=False,chronological=False):

171 """ 172 Construct a phasediagram, using frequency nu0. 173 174 Possibility to include a linear frequency shift D. 175 176 If needed, the indices can be returned. To 'unphase', just do: 177 178 Example usage: 179 180 >>> original = np.random.normal(size=10) 181 >>> indices = np.argsort(original) 182 >>> inv_indices = np.argsort(indices) 183 >>> all(original == np.take(np.take(original,indices),inv_indices)) 184 True 185 186 Optionally, the zero time point can be given, it will be subtracted from all 187 time points 188 189 Joris De Ridder, Pieter Degroote 190 191 @param time: time points [0..Ntime-1] 192 @type time: ndarray 193 @param signal: observed data points [0..Ntime-1] 194 @type signal: ndarray 195 @param nu0: frequency to compute the phasediagram with (inverse unit of time) 196 @type nu0: float 197 @param t0: zero time point, if None: t0 = time[0] 198 @type t0: float 199 @param D: frequency shift 200 @type D: float 201 @param chronological: set to True if you want to have a list of every phase 202 separately 203 @type chronological: boolean 204 @param return_indices: flag to return indices 205 @type return_indices: boolean 206 @return: phase points (sorted), corresponding observations 207 @rtype: ndarray,ndarray(, ndarray) 208 """ 209 cc = 173.144632674 # AU per day 210 if (t0 is None): t0 = 0. 211 if forb is None: 212 phase = np.fmod(nu0 * (time-t0) + D/2. * (time-t0)**2,1.0) 213 else: 214 alpha = nu0*asini/cc 215 phase = np.fmod(nu0 * (time-t0) + alpha*(sin(2*pi*forb*time) - sin(2*pi*forb*t0)),1.0) 216 #phase = np.fmod(nu0 * (time-t0) - asini/cc/(2*np.pi*forb)*np.cos(2*np.pi*forb*(time-t0)),1.0) 217 phase = np.where(phase<0,phase+1,phase) 218 219 #-- keep track of the phase number, and prepare lists to add the points 220 if chronological: 221 chainnr = np.floor((time-t0)*nu0) 222 myphase = [[]] 223 mysignl = [[]] 224 currentphase = chainnr[0] 225 #-- divide the signal into separate phase bins 226 for i in xrange(len(signal)): 227 if chainnr[i]==currentphase: 228 myphase[-1].append(phase[i]) 229 mysignl[-1].append(signal[i]) 230 else: 231 myphase[-1] = np.array(myphase[-1]) 232 mysignl[-1] = np.array(mysignl[-1]) 233 currentphase = chainnr[i] 234 myphase.append([phase[i]]) 235 mysignl.append([signal[i]]) 236 myphase[-1] = np.array(myphase[-1]) 237 mysignl[-1] = np.array(mysignl[-1]) 238 #-- if we 239 else: 240 indices = np.argsort(phase) 241 242 #-- possibly only return phase and signal 243 if not return_indices and not chronological: 244 return phase[indices], signal[indices] 245 #-- maybe return phase,signal and indices 246 elif not chronological: 247 return phase[indices], signal[indices], indices 248 #-- or return seperate phase bins 249 else: 250 return myphase,mysignl

251

252 253 254 @check_input 255 -def sine(times, parameters):

256 """ 257 Creates a harmonic function based on given parameters. 258 259 Parameter fields: C{const, ampl, freq, phase}. 260 261 This harmonic function is of the form 262 263 f(p1,...,pm;x1,...,x_n) = S{sum} c_i + a_i sin(2 S{pi} (f_i+ S{phi}_i)) 264 265 where p_i are parameters and x_i are observation times. 266 267 Parameters can be given as a record array containing columns 'ampl', 'freq' 268 and 'phase', and optionally also 'const' (the latter array is summed up so 269 to reduce it to one number). 270 271 C{pars = np.rec.fromarrays([consts,ampls,freqs,phases],names=('const','ampl','freq','phase'))} 272 273 Parameters can also be given as normal 1D numpy array. In that case, it will 274 be transformed into a record array by the C{sine_preppars} function. The 275 array you give must be 1D, optionally contain a constant as the first 276 parameter, followed by 3 three numbers per frequency: 277 278 C{pars = [const, ampl1, freq1, phase1, ampl2, freq2, phase2...]} 279 280 Example usage: We construct a synthetic signal with two frequencies. 281 282 >>> times = np.linspace(0,100,1000) 283 >>> const = 4. 284 >>> ampls = [0.01,0.004] 285 >>> freqs = [0.1,0.234] 286 >>> phases = [0.123,0.234] 287 288 The signal can be made via a record array 289 290 >>> parameters = np.rec.fromarrays([[const,0],ampls,freqs,phases],names = ('const','ampl','freq','phase')) 291 >>> mysine = sine(times,parameters) 292 293 or a normal 1D numpy array 294 295 >>> parameters = np.hstack([const,np.ravel(np.column_stack([ampls,freqs,phases]))]) 296 >>> mysine = sine(times,parameters) 297 298 Of course, the latter is easier if you have your parameters in a list: 299 300 >>> freq1 = [0.01,0.1,0.123] 301 >>> freq2 = [0.004,0.234,0.234] 302 >>> parameters = np.hstack([freq1,freq2]) 303 >>> mysine = sine(times,parameters) 304 305 @param times: observation times 306 @type times: numpy array 307 @param parameters: record array containing amplitudes ('ampl'), frequencies ('freq') 308 phases ('phase') and optionally constants ('const'), or 1D numpy array (see above). 309 B{Warning:} record array should have shape (n,)! 310 @type parameters: record array or 1D array 311 @return: sine signal (same shape as C{times}) 312 @rtype: array 313 """ 314 if 'const' in parameters.dtype.names: 315 total_fit = parameters['const'].sum() 316 else: 317 total_fit = 0. 318 for i in xrange(len(parameters)): 319 total_fit += parameters['ampl'][i]*sin(2*pi*(parameters['freq'][i]*times+parameters['phase'][i])) 320 return total_fit

321

322 323 324 325 # 326 # 327 # def sine_freqshift(times,parameters,t0=None): 328 # """ 329 # Creates a sine function with a linear frequency shift. 330 # 331 # Parameter fields: C{const, ampl, freq, phase, D}. 332 # 333 # Similar to C{sine}, but with extra column 'D', which is the linear frequency 334 # shift parameter. 335 # 336 # Only accepts record arrays as parameters (for the moment). 337 # 338 # @param times: observation times 339 # @type times: numpy array 340 # @param parameters: record array containing amplitudes ('ampl'), frequencies ('freq'), 341 # phases ('phase'), frequency shifts ('D') and optionally constants ('const') 342 # @type parameters: record array 343 # @return: sine signal with frequency shift (same shape as C{times}) 344 # @rtype: array 345 # """ 346 # #-- take epoch into account 347 # if t0 is None: t0 = times[0] 348 # 349 # #-- fit the signal 350 # signal = np.zeros(len(times)) 351 # if 'const' in parameters.dtype.names: 352 # signal += np.sum(parameters['const']) 353 # 354 # for par in parameters: 355 # frequency = (par['freq'] + par['D']/2.*(times-t0)) * (times-t0) 356 # signal += par['ampl'] * sin(2*pi*(frequency + par['phase'])) 357 # 358 # return signal 359 # 360 # 361 # def sine_orbit(times,parameters,t0=None,nmax=10): 362 # """ 363 # Creates a sine function with a sinusoidal frequency shift. 364 # 365 # Parameter fields: C{const, ampl, freq, phase, asini, omega}. 366 # 367 # Similar to C{sine}, but with extra column 'asini' and 'forb', which are the 368 # orbital parameters. forb in cycles/day or something similar, asini in au. 369 # 370 # For eccentric orbits, add longitude of periastron 'omega' (radians) and 371 # 'ecc' (eccentricity). 372 # 373 # Only accepts record arrays as parameters (for the moment). 374 # 375 # @param times: observation times 376 # @type times: numpy array 377 # @param parameters: record array containing amplitudes ('ampl'), frequencies ('freq'), 378 # phases ('phase'), frequency shifts ('D') and optionally constants ('const') 379 # @type parameters: record array 380 # @return: sine signal with frequency shift (same shape as C{times}) 381 # @rtype: array 382 # """ 383 # #-- take epoch into account 384 # if t0 is None: t0 = times[0] 385 # 386 # def ane(n,e): return 2.*np.sqrt(1-e**2)/e/n*jn(n,n*e) 387 # def bne(n,e): return 1./n*(jn(n-1,n*e)-jn(n+1,n*e)) 388 # 389 # #-- fit the signal 390 # signal = np.zeros(len(times)) 391 # names = parameters.dtype.names 392 # if 'const' in names: 393 # signal += np.sum(parameters['const']) 394 # 395 # cc = 173.144632674 # speed of light in AU/d 396 # for par in parameters: 397 # alpha = par['freq']*par['asini']/cc 398 # if not 'ecc' in names:# or par['ecc']==0: 399 # frequency = par['freq']*(times-t0) + \ 400 # alpha*(sin(2*pi*par['forb']*times) - sin(2*pi*par['forb']*t0)) 401 # else: 402 # e,omega = par['ecc'],par['omega'] 403 # ns = np.arange(1,nmax+1,1) 404 # ans,bns = np.array([[ane(n,e),bne(n,e)] for n in ns]).T 405 # ksins = sqrt(ans**2*cos(omega)**2+bns**2*sin(omega)**2) 406 # thns = arctan(bns/ans*tan(omega)) 407 # tau = -np.sum(bns*sin(omega)) 408 # frequency = par['freq']*(times-t0) + \ 409 # alpha*(np.sum(np.array([ksins[i]*sin(2*pi*ns[i]*par['forb']*(times-t0)+thns[i]) for i in range(nmax)]),axis=0)+tau) 410 # signal += par['ampl'] * sin(2*pi*(frequency + par['phase'])) 411 # 412 # return signal 413 414 415 416 417 -def periodic_spline(times,parameters,t0=None):

418 """ 419 Evaluate a periodic spline function 420 421 Parameter fields: C{freq, knots, coeffs, degree}, as returned from the 422 corresponding fitting function L{fit.periodic_spline}. 423 424 @param times: observation times 425 @type times: numpy array 426 @param parameters: [frequency,"cj" spline coefficients] 427 @return: sine signal with frequency shift (same shape as C{times}) 428 @rtype: array 429 """ 430 #-- get time offset 431 if t0 is None: 432 t0 = times[0] 433 434 mytrend = 0. 435 #-- reconstruct entire time series, instead of phased version 436 for i in range(len(parameters)): 437 frequency = parameters[i]['freq'] 438 #-- reconstruct coefficients array 439 coefficients = (parameters[i]['knots'],parameters[i]['coeffs'],parameters[i]['degree']) 440 phases = np.fmod((times-t0) * frequency,1.0) 441 mytrend += splev(phases,coefficients) 442 443 return mytrend

444

445 446 447 # @check_input 448 # def kepler(times,parameters,itermax=8): 449 # """ 450 # Construct a radial velocity curve due to Kepler orbit(s). 451 # 452 # Parameter fields: C{gamma, P, T0, e, omega, K} 453 # 454 # @param times: observation times 455 # @type times: numpy array 456 # @param parameters: record array containing periods ('P' in days), 457 # times of periastron ('T0' in days), eccentricities ('e'), 458 # longitudes of periastron ('omega' in radians), amplitude ('K' in km/s) and 459 # systemic velocity ('gamma' in km/s) 460 # phases ('phase'), frequency shifts ('D') and optionally constants ('const') 461 # @type parameters: record array 462 # @param itermax: maximum number of iterations for solving Kepler equation 463 # @type itermax: integer 464 # @return: radial velocity curve (same shape as C{times}) 465 # @rtype: array 466 # """ 467 # if 'gamma' in parameters.dtype.names: 468 # RVfit = parameters['gamma'].sum() 469 # else: 470 # RVfit = 0 471 # for pars in parameters: 472 # p = [pars['P'],pars['T0'],pars['e'],pars['omega'],pars['K'],0] 473 # RVfit += keplerorbit.radial_velocity(p,times=times,itermax=itermax) 474 # return RVfit 475 # 476 # @check_input 477 # def kepler_diffcorr(times,parameters,itermax=8): 478 # """ 479 # Compute coefficients for differential correction method. 480 # 481 # See page 97-98 of Hilditch's book "An introduction to close binary stars". 482 # """ 483 # N = len(times) 484 # if 'gamma' in parameters.dtype.names: 485 # RV0 = parameters['gamma'].sum() 486 # else: 487 # RV0 = 0. 488 # P,T0,e,omega,K = parameters['P'],parameters['T0'],parameters['e'],\ 489 # parameters['omega'],parameters['K'] 490 # freq = 1./P 491 # x0 = T0*2*np.pi*freq 492 # omega = omega*np.ones(N) 493 # #-- calculate true anomaly 494 # E,theta = keplerorbit.true_anomaly(times*2*np.pi*freq-x0,e,itermax=itermax) 495 # #-- calculate coefficients: 496 # th_om = theta + omega 497 # sin_th_om = np.sin(th_om) 498 # cos_th_om = np.cos(th_om) 499 # cos_om = np.cos(omega) 500 # sin_om = np.sin(omega) 501 # cos_th = np.cos(theta) 502 # sin_th = np.sin(theta) 503 # c_deltaK = cos_th_om + e*cos_om 504 # c_deltae = K*(cos_om - (sin_th_om * sin_th * (2+e*cos_th))/(1-e**2)) 505 # c_deltao = -K*(sin_th_om + e*sin_om) 506 # c_deltaT = sin_th_om*(1+e*cos_th)**2*2*np.pi*K/P/(1-e**2)**1.5 507 # c_deltaP = sin_th_om*(1+e*cos_th)**2*2*np.pi*K*(times-T0)/P**2/(1-e**2)**1.5 508 # c_deltag = np.ones(N) 509 # A = np.vstack([c_deltag,c_deltaP,c_deltaT,c_deltae,c_deltao,c_deltaK]).T 510 # return A 511 512 # 513 # @check_input 514 # def binary(times,parameters,n1=None,itermax=8): 515 # """ 516 # Simultaneous evaluation of two binary components. 517 # 518 # parameter fields must be C{'P','T0','e','omega','K1', 'K2'}. 519 # """ 520 # if n1 is None: 521 # n1 = len(times)/2 522 # if 'gamma' in parameters.dtype.names: 523 # RVfit = parameters['gamma'].sum() 524 # else: 525 # RVfit = 0 526 # RVfit = RVfit*np.ones(len(times)) 527 # p1 = [parameters['P'],parameters['T0'],parameters['e'],parameters['omega'],parameters['K1'],0] 528 # p2 = [parameters['P'],parameters['T0'],parameters['e'],parameters['omega'],parameters['K2'],0] 529 # RVfit[:n1] += keplerorbit.radial_velocity(p1,times=times[:n1],itermax=itermax) 530 # RVfit[n1:] += keplerorbit.radial_velocity(p2,times=times[n1:],itermax=itermax) 531 # return RVfit 532 # 533 534 535 @check_input 536 -def box(times,parameters,t0=None):

537 """ 538 Evaluate a box transit model. 539 540 Parameter fields: C{cont, freq, ingress, egress, depth} 541 542 Parameters [[frequency,depth, fractional start, fraction end, continuum]] 543 @rtype: array 544 """ 545 if t0 is None: t0 = 0. 546 547 #-- continuum 548 parameters = np.asarray(parameters) 549 model = np.ones(len(times))*sum(parameters['cont']) 550 551 for parameter in parameters: 552 #-- set up the model, define which point is where in a 553 # phasediagram 554 phase = np.fmod((times - t0) * parameter['freq'], 1.0) 555 phase = np.where(phase<0,phase+1,phase) 556 transit_place = (parameter['ingress']<=phase) & (phase<=parameter['egress']) 557 model = np.where(transit_place,model-parameter['depth'],model) 558 return model

559

560 561 562 -def spots(times,spots,t0=None,**kwargs):

563 """ 564 Spotted model from Lanza (2003). 565 Extract from Carrington Map. 566 567 Included: 568 - presence of faculae through Q-factor (set constant to 5 or see 569 Chapman et al. 1997 or Solank & Unruh for a scaling law, since this 570 factor decreases strongly with decreasing rotation period. 571 - differential rotation through B-value (B/period = 0.2 for the sun) 572 - limb-darkening through coefficients ap, bp and cp 573 - non-equidistant time series 574 - individual starspot evolution and lifetime following a Gaussian law, 575 (see Hall & Henry (1993) for a relation between stellar radius and 576 sunspot lifetime), through an optional 3rd (time of maximum) and 577 4th parameter (decay time) 578 579 Not included: 580 - distinction between umbra's and penumbra 581 - sunspot structure 582 583 Parameters of global star: 584 - i_sun: inclination angle, i=pi/2 is edge-on, i=0 is pole on. 585 - l0: 'epoch angle', if L0=0, spot with coordinate (0,0) 586 starts in center of disk, if L0=0, spot is in center of 587 disk after half a rotation period 588 - period: equatorial period 589 - b: differential rotation b = Omega_pole - Omega_eq where Omega is in radians per time unit (angular velocity) 590 - ap,bp,cp: limb darkening 591 592 Parameters of spots: 593 - lambda: 'greenwich' coordinate (longitude) (pi means push to back of star if l0=0) 594 - theta: latitude (latitude=0 means at equator, latitude=pi/2 means on pole (B{not 595 colatitude}) 596 - A: size 597 - maxt0 (optional): time of maximum size 598 - decay (optional): exponential decay speed 599 600 601 Note: alpha = (Omega_eq - Omega_p) / Omega_eq 602 alpha = -b /2*pi * period 603 604 Use for fitting! 605 """ 606 #-- set zeropoint 607 if t0 is None: t0 = 0. 608 609 #-- simulated light curve template 610 signal = np.zeros(len(times)) 611 612 #-- run through all timepoints 613 this_signal = 0 614 for spot in spots: 615 #-- compute limb-darkening law 616 C = 1. / (0.25 * spots[0]['ap'] + 2./3.*spots[0]['bp'] + 1./2.*spots[0]['cp']) 617 #-- position and size of spot 618 lambda_i = spot['lambda'] 619 A_i = spot['A'] 620 if 'decay' in spot.dtype.names: 621 #-- include evolving spot 622 A_i *= np.exp(-((spot['maxt0'] - times)/spot['decay'])**2) 623 #-- find the velocity at this position (possible dependent on theta) 624 Omega_theta = 2*pi/spots[0]['P_eq'] + spots[0]['b'] * sin(spot['theta'])**2 625 #-- Correct the position of lambda_i according to the position 626 lambda_i += Omega_theta*(times-t0) 627 #-- Calculate mu_i = cos(psi_i) with psi_i the angle between the normal 628 # to the ith active region and the line of sight 629 mu_i = cos(spots[0]['i_sun']) * sin(spot['theta']) + sin(spots[0]['i_sun'])*cos(spot['theta'])*cos(lambda_i - spots[0]['l0']) 630 #-- Calculate total irradiance for this spot and add it to the signal 631 irradiance_spot = C*spots[0]['s0'] * A_i * mu_i * \ 632 (spots[0]['ap'] + spots[0]['bp']*mu_i + spots[0]['cp']*mu_i**2) *\ 633 ((spots[0]['cs']-1) + spots[0]['q']*(spots[0]['cf'] + spots[0]['cf_']*mu_i - 1)) 634 #-- only take visible spots into account 635 this_signal = where(mu_i<0,this_signal,this_signal+irradiance_spot) 636 this_signal += spots[0]['s0'] + spots[0]['sr'] 637 return this_signal

638

639 640 641 #} 642 643 644 #{ General purpose 645 646 @check_input 647 -def gauss(x, parameters):

648 """ 649 Evaluate a gaussian fit. 650 651 Parameter fields: C{'ampl','mu','sigma'} and optionally C{'const'}. 652 653 Evaluating one Gaussian: 654 655 >>> x = np.linspace(-20,20,10000) 656 >>> pars = [0.5,0.,3.] 657 >>> y = gauss(x,pars) 658 >>> p = pl.plot(x,y,'k-') 659 660 Evaluating multiple Gaussian, with and without constants: 661 662 >>> pars = [0.5,0.,3.,0.1,-10,1.,.1] 663 >>> y = gauss(x,pars) 664 >>> p = pl.plot(x,y,'r-') 665 666 @parameter x: domain 667 @type x: ndarray 668 @parameter parameters: record array. 669 @type x: numpy record array 670 @return: fitted signal 671 @rtype: array 672 """ 673 if 'const' in parameters.dtype.names: 674 C = parameters['const'].sum() 675 else: 676 C = 0. 677 678 y = C 679 for pars in parameters: 680 y += pars['ampl'] * np.exp( -(x-pars['mu'])**2 / (2.*pars['sigma']**2)) 681 682 return y

683

684 @check_input 685 -def lorentz(x,parameters):

686 """ 687 Evaluate Lorentz profile 688 689 P(f) = A / ( (x-mu)**2 + gamma**2) 690 691 Parameters fields: C{ampl,mu,gamma} and optionally C{const}. 692 693 Evaluating one Lorentzian: 694 695 >>> x = np.linspace(-20,20,10000) 696 >>> pars = [0.5,0.,3.] 697 >>> y = lorentz(x,pars) 698 >>> p = pl.figure() 699 >>> p = pl.plot(x,y,'k-') 700 701 Evaluating multiple Lorentzians, with and without constants: 702 703 >>> pars = [0.5,0.,3.,0.1,-10,1.,.1] 704 >>> y = lorentz(x,pars) 705 >>> p = pl.plot(x,y,'r-') 706 """ 707 if 'const' in parameters.dtype.names: 708 y = parameters['const'].sum() 709 else: 710 y = 0. 711 for par in parameters: 712 y += par['ampl'] / ((x-par['mu'])**2 + par['gamma']**2) 713 return y

714

715 716 717 @check_input 718 -def voigt(x,parameters):

719 """ 720 Evaluate a Voigt profile. 721 722 Field parameters: C{ampl, mu, gamma, sigma} and optionally C{const} 723 724 Function:: 725 726 z = (x + gamma*i) / (sigma*sqrt(2)) 727 V = A * Real[cerf(z)] / (sigma*sqrt(2*pi)) 728 729 >>> x = np.linspace(-20,20,10000) 730 >>> pars1 = [0.5,0.,2.,2.] 731 >>> pars2 = [0.5,0.,1.,2.] 732 >>> pars3 = [0.5,0.,2.,1.] 733 >>> y1 = voigt(x,pars1) 734 >>> y2 = voigt(x,pars2) 735 >>> y3 = voigt(x,pars3) 736 >>> p = pl.figure() 737 >>> p = pl.plot(x,y1,'b-') 738 >>> p = pl.plot(x,y2,'g-') 739 >>> p = pl.plot(x,y3,'r-') 740 741 Multiple voigt profiles: 742 743 >>> pars4 = [0.5,0.,2.,1.,0.1,-3,0.5,0.5,0.01] 744 >>> y4 = voigt(x,pars4) 745 >>> p = pl.plot(x,y4,'c-') 746 747 @rtype: array 748 @return: Voigt profile 749 """ 750 if 'const' in parameters.dtype.names: 751 y = parameters['const'].sum() 752 else: 753 y = 0. 754 for par in parameters: 755 x = x-par['mu'] 756 z = (x+1j*par['gamma'])/(par['sigma']*sqrt(2)) 757 y += par['ampl']*_complex_error_function(z).real/(par['sigma']*sqrt(2*pi)) 758 return y

759

760 @check_input 761 -def power_law(x,parameters):

762 """ 763 Evaluate a power law. 764 765 Parameter fields: C{A, B, C}, optionally C{const} 766 767 Function: 768 769 P(f) = A / (1+ Bf)**C + const 770 771 >>> x = np.linspace(0,10,1000) 772 >>> pars1 = [0.5,1.,1.] 773 >>> pars2 = [0.5,1.,2.] 774 >>> pars3 = [0.5,2.,1.] 775 >>> y1 = power_law(x,pars1) 776 >>> y2 = power_law(x,pars2) 777 >>> y3 = power_law(x,pars3) 778 >>> p = pl.figure() 779 >>> p = pl.plot(x,y1,'b-') 780 >>> p = pl.plot(x,y2,'g-') 781 >>> p = pl.plot(x,y3,'r-') 782 783 Multiple power laws: 784 >>> pars4 = pars1+pars2+pars3 785 >>> y4 = power_law(x,pars4) 786 >>> p = pl.plot(x,y4,'c-') 787 788 @parameter x: domain 789 @type x: ndarray 790 @parameter parameters: record array. 791 @type x: numpy record array 792 @return: fitted signal 793 @rtype: array 794 """ 795 if 'const' in parameters.dtype.names: 796 y = par['const'] 797 else: 798 y = 0. 799 for par in parameters: 800 y += par['A'] / (1+ (par['B']*x)**par['C']) 801 return y

802

803 #} 804 805 806 807 #{ Convert record arrays to flat arrays and back 808 809 810 -def sine_preppars(pars):

811 """ 812 Prepare sine function parameters in correct form for evaluating/fitting. 813 814 If you input a record array, this function will output a 1D numpy array 815 containing only the independent parameters for use in nonlinear fitting 816 algorithms. 817 818 If you input a 1D numpy array, it will output a record array. 819 820 @param pars: input parameters in record array or normal array form 821 @type pars: record array or normal numpy array 822 @return: input parameters in normal array form or record array 823 @rtype: normal numpy array or record array 824 """ 825 #-- from record array to flat 826 if pars.dtype.names: 827 # with constant 828 if 'const' in pars.dtype.names: 829 converted_pars = np.zeros(3*len(pars)+1) 830 converted_pars[0] = pars['const'].sum() 831 converted_pars[1::3] = pars['ampl'] 832 converted_pars[2::3] = pars['freq'] 833 converted_pars[3::3] = pars['phase'] 834 # without constant 835 else: 836 converted_pars = np.zeros(3*len(pars)) 837 converted_pars[0::3] = pars['ampl'] 838 converted_pars[1::3] = pars['freq'] 839 converted_pars[2::3] = pars['phase'] 840 #-- from flat to record array 841 else: 842 # with constant 843 if len(pars)%3==1: 844 converted_pars = np.zeros((4,(len(pars)-1)/3)) 845 converted_pars[0,0] = pars[0] 846 converted_pars[1] = pars[1::3] 847 converted_pars[2] = pars[2::3] 848 converted_pars[3] = pars[3::3] 849 names = ['const','ampl','freq','phase'] 850 # without constant 851 else: 852 converted_pars = np.zeros((3,len(pars)/3)) 853 converted_pars[0] = pars[0::3] 854 converted_pars[1] = pars[1::3] 855 converted_pars[2] = pars[2::3] 856 names = ['ampl','freq','phase'] 857 converted_pars = np.rec.fromarrays(converted_pars,names=names) 858 return converted_pars

859

860 861 # def kepler_preppars(pars): 862 # """ 863 # Prepare Kepler orbit parameters in correct form for evaluating/fitting. 864 # 865 # If you input a record array, this function will output a 1D numpy array 866 # containing only the independent parameters for use in nonlinear fitting 867 # algorithms. 868 # 869 # If you input a 1D numpy array, it will output a record array. 870 # 871 # @param pars: input parameters in record array or normal array form 872 # @type pars: record array or normal numpy array 873 # @return: input parameters in normal array form or record array 874 # @rtype: normal numpy array or record array 875 # """ 876 # #-- from record array to flat 877 # if pars.dtype.names: 878 # # with systemic velocity 879 # if 'gamma' in pars.dtype.names: 880 # converted_pars = np.zeros(5*len(pars)+1) 881 # converted_pars[0] = pars['gamma'].sum() 882 # converted_pars[1::5] = pars['P'] 883 # converted_pars[2::5] = pars['T0'] 884 # converted_pars[3::5] = pars['e'] 885 # converted_pars[4::5] = pars['omega'] 886 # converted_pars[5::5] = pars['K'] 887 # else: 888 # converted_pars = np.zeros(5*len(pars)) 889 # converted_pars[0::5] = pars['P'] 890 # converted_pars[1::5] = pars['T0'] 891 # converted_pars[2::5] = pars['e'] 892 # converted_pars[3::5] = pars['omega'] 893 # converted_pars[4::5] = pars['K'] 894 # #-- from flat to record array 895 # else: 896 # # with constant 897 # if len(pars)%5==1: 898 # converted_pars = np.zeros((6,(len(pars)-1)/5)) 899 # converted_pars[0,0] = pars[0] 900 # converted_pars[1] = pars[1::5] 901 # converted_pars[2] = pars[2::5] 902 # converted_pars[3] = pars[3::5] 903 # converted_pars[4] = pars[4::5] 904 # converted_pars[5] = pars[5::5] 905 # names = ['gamma','P','T0','e','omega','K'] 906 # # without constant 907 # else: 908 # converted_pars = np.zeros((5,len(pars)/5)) 909 # converted_pars[0] = pars[0::5] 910 # converted_pars[1] = pars[1::5] 911 # converted_pars[2] = pars[2::5] 912 # converted_pars[3] = pars[3::5] 913 # converted_pars[4] = pars[4::5] 914 # names = ['P','T0','e','omega','K'] 915 # converted_pars = np.rec.fromarrays(converted_pars,names=names) 916 # return converted_pars 917 # 918 # 919 # def binary_preppars(pars): 920 # """ 921 # Prepare Kepler orbit parameters in correct form for evaluating/fitting. 922 # 923 # If you input a record array, this function will output a 1D numpy array 924 # containing only the independent parameters for use in nonlinear fitting 925 # algorithms. 926 # 927 # If you input a 1D numpy array, it will output a record array. 928 # 929 # @param pars: input parameters in record array or normal array form 930 # @type pars: record array or normal numpy array 931 # @return: input parameters in normal array form or record array 932 # @rtype: normal numpy array or record array 933 # """ 934 # #-- from record array to flat 935 # if pars.dtype.names: 936 # # with systemic velocity 937 # if 'gamma' in pars.dtype.names: 938 # converted_pars = np.zeros(6*len(pars)+1) 939 # converted_pars[0] = pars['gamma'].sum() 940 # converted_pars[1::6] = pars['P'] 941 # converted_pars[2::6] = pars['T0'] 942 # converted_pars[3::6] = pars['e'] 943 # converted_pars[4::6] = pars['omega'] 944 # converted_pars[5::6] = pars['K1'] 945 # converted_pars[6::6] = pars['K2'] 946 # else: 947 # converted_pars = np.zeros(6*len(pars)) 948 # converted_pars[0::6] = pars['P'] 949 # converted_pars[1::6] = pars['T0'] 950 # converted_pars[2::6] = pars['e'] 951 # converted_pars[3::6] = pars['omega'] 952 # converted_pars[4::6] = pars['K1'] 953 # converted_pars[5::6] = pars['K2'] 954 # #-- from flat to record array 955 # else: 956 # # with constant 957 # if len(pars)%6==1: 958 # converted_pars = np.zeros((7,(len(pars)-1)/6)) 959 # converted_pars[0,0] = pars[0] 960 # converted_pars[1] = pars[1::6] 961 # converted_pars[2] = pars[2::6] 962 # converted_pars[3] = pars[3::6] 963 # converted_pars[4] = pars[4::6] 964 # converted_pars[5] = pars[5::6] 965 # converted_pars[6] = pars[6::6] 966 # names = ['gamma','P','T0','e','omega','K1','K2'] 967 # # without constant 968 # else: 969 # converted_pars = np.zeros((6,len(pars)/6)) 970 # converted_pars[0] = pars[0::6] 971 # converted_pars[1] = pars[1::6] 972 # converted_pars[2] = pars[2::6] 973 # converted_pars[3] = pars[3::6] 974 # converted_pars[4] = pars[4::6] 975 # converted_pars[5] = pars[5::6] 976 # names = ['P','T0','e','omega','K1','K2'] 977 # converted_pars = np.rec.fromarrays(converted_pars,names=names) 978 # return converted_pars 979 # 980 981 982 -def box_preppars(pars):

983 """ 984 Prepare sine function parameters in correct form for evaluating/fitting. 985 986 If you input a record array, this function will output a 1D numpy array 987 containing only the independent parameters for use in nonlinear fitting 988 algorithms. 989 990 If you input a 1D numpy array, it will output a record array. 991 992 @param pars: input parameters in record array or normal array form 993 @type pars: record array or normal numpy array 994 @return: input parameters in normal array form or record array 995 @rtype: normal numpy array or record array 996 """ 997 #-- from record array to flat 998 if pars.dtype.names: 999 converted_pars = np.ones(4*len(pars)+1) 1000 converted_pars[0] += sum(pars['cont']) 1001 converted_pars[1::4] = pars['freq'] 1002 converted_pars[2::4] = pars['depth'] 1003 converted_pars[3::4] = pars['ingress'] 1004 converted_pars[4::4] = pars['egress'] 1005 #-- from flat to record array 1006 else: 1007 converted_pars = np.ones((5,(len(pars)-1)/4)) 1008 converted_pars[0,0] = pars[0] 1009 converted_pars[1] = pars[1::4] 1010 converted_pars[2] = pars[2::4] 1011 converted_pars[3] = pars[3::4] 1012 converted_pars[4] = pars[4::4] 1013 names = ['cont','freq','depth','ingress','egress'] 1014 converted_pars = np.rec.fromarrays(converted_pars,names=names) 1015 return converted_pars

1016

1017 1018 -def gauss_preppars(pars):

1019 """ 1020 Prepare gauss function parameters in correct form for evaluating/fitting. 1021 1022 If you input a record array, this function will output a 1D numpy array 1023 containing only the independent parameters for use in nonlinear fitting 1024 algorithms. 1025 1026 If you input a 1D numpy array, it will output a record array. 1027 1028 @param pars: input parameters in record array or normal array form 1029 @type pars: record array or normal numpy array 1030 @return: input parameters in normal array form or record array 1031 @rtype: normal numpy array or record array 1032 """ 1033 #-- from record array to flat 1034 if pars.dtype.names: 1035 if 'const' in pars.dtype.names: 1036 converted_pars = np.zeros(3*len(pars)+1) 1037 converted_pars[-1] = pars['const'].sum() 1038 else: 1039 converted_pars = np.zeros(3*len(pars)) 1040 converted_pars[0::3] = pars['ampl'] 1041 converted_pars[1::3] = pars['mu'] 1042 converted_pars[2::3] = pars['sigma'] 1043 #-- from flat to record array 1044 else: 1045 if len(pars)%3==0: 1046 converted_pars = np.zeros((3,len(pars)/3)) 1047 names = ['ampl','mu','sigma'] 1048 else: 1049 converted_pars = np.zeros((4,(len(pars)-1)/3)) 1050 converted_pars[3,0] = pars[-1] 1051 names = ['ampl','mu','sigma','const'] 1052 converted_pars[0] = pars[0:-1:3] 1053 converted_pars[1] = pars[1::3] 1054 converted_pars[2] = pars[2::3] 1055 converted_pars = np.rec.fromarrays(converted_pars,names=names) 1056 return converted_pars

1057

1058 -def lorentz_preppars(pars):

1059 """ 1060 Prepare Lorentz function parameters in correct form for evaluating/fitting. 1061 1062 If you input a record array, this function will output a 1D numpy array 1063 containing only the independent parameters for use in nonlinear fitting 1064 algorithms. 1065 1066 If you input a 1D numpy array, it will output a record array. 1067 1068 @param pars: input parameters in record array or normal array form 1069 @type pars: record array or normal numpy array 1070 @return: input parameters in normal array form or record array 1071 @rtype: normal numpy array or record array 1072 """ 1073 #-- from record array to flat 1074 if pars.dtype.names: 1075 if 'const' in pars.dtype.names: 1076 converted_pars = np.zeros(3*len(pars)+1) 1077 converted_pars[-1] = pars['const'].sum() 1078 else: 1079 converted_pars = np.zeros(3*len(pars)) 1080 converted_pars[0::3] = pars['ampl'] 1081 converted_pars[1::3] = pars['mu'] 1082 converted_pars[2::3] = pars['gamma'] 1083 #-- from flat to record array 1084 else: 1085 if len(pars)%3==0: 1086 converted_pars = np.zeros((3,len(pars)/3)) 1087 names = ['ampl','mu','gamma'] 1088 else: 1089 converted_pars = np.zeros((4,(len(pars)-1)/3)) 1090 converted_pars[3,0] = pars[-1] 1091 names = ['ampl','mu','gamma','const'] 1092 converted_pars[0] = pars[0:-1:3] 1093 converted_pars[1] = pars[1::3] 1094 converted_pars[2] = pars[2::3] 1095 converted_pars = np.rec.fromarrays(converted_pars,names=names) 1096 return converted_pars

1097

1098 -def voigt_preppars(pars):

1099 """ 1100 Prepare voigt function parameters in correct form for evaluating/fitting. 1101 1102 If you input a record array, this function will output a 1D numpy array 1103 containing only the independent parameters for use in nonlinear fitting 1104 algorithms. 1105 1106 If you input a 1D numpy array, it will output a record array. 1107 1108 @param pars: input parameters in record array or normal array form 1109 @type pars: record array or normal numpy array 1110 @return: input parameters in normal array form or record array 1111 @rtype: normal numpy array or record array 1112 """ 1113 #-- from record array to flat 1114 if pars.dtype.names: 1115 if 'const' in pars.dtype.names: 1116 converted_pars = np.zeros(4*len(pars)+1) 1117 converted_pars[-1] = pars['const'].sum() 1118 else: 1119 converted_pars = np.zeros(4*len(pars)) 1120 converted_pars[0::4] = pars['ampl'] 1121 converted_pars[1::4] = pars['mu'] 1122 converted_pars[2::4] = pars['sigma'] 1123 converted_pars[3::4] = pars['gamma'] 1124 #-- from flat to record array 1125 else: 1126 if len(pars)%4==0: 1127 converted_pars = np.zeros((4,len(pars)/4)) 1128 names = ['ampl','mu','sigma','gamma'] 1129 else: 1130 converted_pars = np.zeros((5,(len(pars)-1)/4)) 1131 converted_pars[4,0] = pars[-1] 1132 names = ['ampl','mu','sigma','gamma','const'] 1133 converted_pars[0] = pars[0:-1:4] 1134 converted_pars[1] = pars[1::4] 1135 converted_pars[2] = pars[2::4] 1136 converted_pars[3] = pars[3::4] 1137 converted_pars = np.rec.fromarrays(converted_pars,names=names) 1138 return converted_pars

1139

1140 1141 -def power_law_preppars(pars):

1142 """ 1143 Prepare gauss function parameters in correct form for evaluating/fitting. 1144 1145 If you input a record array, this function will output a 1D numpy array 1146 containing only the independent parameters for use in nonlinear fitting 1147 algorithms. 1148 1149 If you input a 1D numpy array, it will output a record array. 1150 1151 @param pars: input parameters in record array or normal array form 1152 @type pars: record array or normal numpy array 1153 @return: input parameters in normal array form or record array 1154 @rtype: normal numpy array or record array 1155 """ 1156 #-- from record array to flat 1157 if pars.dtype.names: 1158 if 'const' in pars.dtype.names: 1159 converted_pars = np.zeros(3*len(pars)+1) 1160 converted_pars[-1] = pars['const'].sum() 1161 else: 1162 converted_pars = np.zeros(3*len(pars)) 1163 converted_pars[0::3] = pars['A'] 1164 converted_pars[1::3] = pars['B'] 1165 converted_pars[2::3] = pars['C'] 1166 #-- from flat to record array 1167 else: 1168 if len(pars)%3==0: 1169 converted_pars = np.zeros((3,len(pars)/3)) 1170 names = ['A','B','C'] 1171 else: 1172 converted_pars = np.zeros((4,(len(pars)-1)/3)) 1173 converted_pars[3,0] = pars[-1] 1174 names = ['A','B','C','const'] 1175 converted_pars[0] = pars[0:-1:3] 1176 converted_pars[1] = pars[1::3] 1177 converted_pars[2] = pars[2::3] 1178 converted_pars = np.rec.fromarrays(converted_pars,names=names) 1179 return converted_pars

1180 #}

1181 1182 #{ Helper functions 1183 1184 -def _complex_error_function(x):

1185 """ 1186 Complex error function 1187 """ 1188 cef_value = np.exp(-x**2)*(1-erf(-1j*x)) 1189 if sum(np.isnan(cef_value))>0: 1190 logger.warning("Complex Error function: NAN encountered, results are biased") 1191 noisnans = np.compress(1-np.isnan(cef_value),cef_value) 1192 try: 1193 last_value = noisnans[-1] 1194 except: 1195 last_value = 0 1196 logger.warning("Complex Error function: all values are NAN, results are wrong") 1197 cef_value = np.where(np.isnan(cef_value),last_value,cef_value) 1198 1199 return cef_value

1200 1201 #} 1202 1203 1204 if __name__=="__main__": 1205 import doctest 1206 import pylab as pl 1207 import sys 1208 doctest.testmod() 1209 pl.show() 1210

Source Code for Module ComboCode.cc.ivs.sigproc.evaluate