1
2
3 """
4 Tools for inter- and extrapolation.
5
6 Author: R. Lombaert
7
8 """
9
10 from scipy import array, hstack
11 from scipy import exp
12 from scipy.optimize import leastsq
13 from scipy import isnan
14
15 from cc.plotting import Plotting2
16
17
19
20 """
21 Calculate residuals between function evaluation and given y-values.
22
23 @param p: parameters for function evaluation
24 @type p: list
25 @param x: the x-grid evaluated by function
26 @type x: array
27 @param y: The given y-values
28 @type y: array
29
30 @keyword func: type of function f(x), can be ['power','exp','power10',\
31 'square','linear'] for now.
32
33 (default: 'power')
34 @type func: string
35 @return: residuals as a function of the x-grid
36 @rtype: array
37
38 """
39
40 x,y,p = array(x),array(y),array(p)
41 err = y - pEval(x,p,func)
42 return err
43
44
45
46 -def pEval(x,p,func='power',x_unit=None):
47
48 """
49 Evaluate a function prescription.
50
51 @param x: The x-grid for function evaluation
52 @type x: array
53 @param p: parameters for function evaluation
54 @type p: list
55
56 @keyword func: type of function f(x), can be ['power','exp','power10',\
57 'square','linear'] for now.
58
59 (default: 'power')
60 @type func: string
61 @keyword x_unit: the unit of the x array, only relevant if func=='planck'.
62 Units other than micron not yet available.
63
64 (default: None)
65 @type x_unit: string
66
67 @return: f(x)
68 @rtype: array
69
70 """
71
72 x,p = array(x),array(p)
73 if func.lower() == 'power':
74 return p[0] + p[2]*(x+p[1])**(p[3])
75 if func.lower() == 'exp':
76 return p[0] + p[2]*(exp(p[3]*(x+p[1])))
77 elif func.lower() == 'power10':
78 return p[0] + p[2]*(10**(p[3]*(x+p[1])))
79 elif func.lower() == 'square':
80 return p[0] + p[2]*(x+p[1]) + p[3]*(x+p[1])**2
81 elif func.lower() == 'linear':
82 return p[0] + p[2]*(x + p[1])
83
84
85
86 -def fitFunction(x_in,y_in,x_out,show=0,func='linear',initial=[1,0,0.5,-1.5],\
87 maxfev=20000):
88
89 '''
90 Fit a template function to input data.
91
92 @param x_in: The input x values
93 @type x_in: list/array
94 @param y_in: The input y values
95 @type y_in: list/array
96 @param x_out: The output x grid
97 @type x_out: list/array
98
99 @keyword show: show the fit result
100
101 (default: 0)
102 @type show: bool
103 @keyword func: The function template (see pEval)
104
105 (default: 'linear')
106 @type func: string
107 @keyword initial: initialisation parameters for your function, number of
108 parameters depends on type of function (see pEval)
109
110 (default: [1,0,0.5,-1.5])
111 @type initial: list(float)
112 @keyword maxfev: Maximum number of calls to the function for the leastsq
113 method. Zero means (1+number of elements in x_in)*100
114
115 (default: 20000)
116 @return: The output y values
117 @rtype: list/array
118
119 '''
120
121 p_lsq = leastsq(getResiduals,initial,args=(x_in,y_in,func),\
122 maxfev=20000)[0]
123 y_out = pEval(x_out,p_lsq,func=func)
124 if show:
125 print p_lsq
126 Plotting2.plotCols(x=[x_in,hstack([x_out,x_in])],\
127 y=[y_in,hstack([y_out,pEval(x_in,p_lsq,func=func)])],\
128 xlogscale=1,ylogscale=1)
129 return y_out
130
131
132
134
135 """
136 Linear interpolation between two points.
137
138 @param x_in: 2 x-points for interpolation
139 @type x_in: list
140 @param y_in: 2 y-points for interpolation corresponding to x_in
141 @type y_in: list
142 @param x: x value for interpolation
143 @type x: float/array
144
145 @return: f(x), where f is linear interpolation from x_in and y_in
146 @rtype: float/array
147
148 """
149
150 x = array(x)
151 try:
152 return y_in[0] + (y_in[0]-y_in[1])/(x_in[0]-x_in[1]) * (x-x_in[0])
153 except ZeroDivisionError:
154 print 'Identical x-coordinates were submitted: Division by zero. ' + \
155 'Aborting.'
156 return
157