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')
42 return check
43
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
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
99 RSS_ = RSS(data,model)
100
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
128 RSS_ = RSS(data,model)
129
130 AIC = n * np.log(RSS_/n) + 2*k + n
131 return AIC
132
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
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
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
217 phase = np.where(phase<0,phase+1,phase)
218
219
220 if chronological:
221 chainnr = np.floor((time-t0)*nu0)
222 myphase = [[]]
223 mysignl = [[]]
224 currentphase = chainnr[0]
225
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
239 else:
240 indices = np.argsort(phase)
241
242
243 if not return_indices and not chronological:
244 return phase[indices], signal[indices]
245
246 elif not chronological:
247 return phase[indices], signal[indices], indices
248
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
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
347 if t0 is None: t0 = times[0]
348
349
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
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
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
396 for par in parameters:
397 alpha = par['freq']*par['asini']/cc
398 if not 'ecc' in names:
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
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
431 if t0 is None:
432 t0 = times[0]
433
434 mytrend = 0.
435
436 for i in range(len(parameters)):
437 frequency = parameters[i]['freq']
438
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
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
494 E,theta = keplerorbit.true_anomaly(times*2*np.pi*freq-x0,e,itermax=itermax)
495
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
548 parameters = np.asarray(parameters)
549 model = np.ones(len(times))*sum(parameters['cont'])
550
551 for parameter in parameters:
552
553
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
607 if t0 is None: t0 = 0.
608
609
610 signal = np.zeros(len(times))
611
612
613 this_signal = 0
614 for spot in spots:
615
616 C = 1. / (0.25 * spots[0]['ap'] + 2./3.*spots[0]['bp'] + 1./2.*spots[0]['cp'])
617
618 lambda_i = spot['lambda']
619 A_i = spot['A']
620 if 'decay' in spot.dtype.names:
621
622 A_i *= np.exp(-((spot['maxt0'] - times)/spot['decay'])**2)
623
624 Omega_theta = 2*pi/spots[0]['P_eq'] + spots[0]['b'] * sin(spot['theta'])**2
625
626 lambda_i += Omega_theta*(times-t0)
627
628
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
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
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
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
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
826 if pars.dtype.names:
827
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
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
841 else:
842
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
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
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
877 if pars.dtype.names:
878
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
895 else:
896
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
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
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
935 if pars.dtype.names:
936
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
955 else:
956
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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