Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1# -*- coding: utf-8 -*- 

2""" 

3 

4Created on Mon Mar 18 15:48:23 2013 

5Author: Josef Perktold 

6 

7TODO: 

8 - test behavior if nans or infs are encountered during the evaluation. 

9 now partially robust to nans, if increasing can be determined or is given. 

10 - rewrite core loop to use for...except instead of while. 

11 

12""" 

13import numpy as np 

14from scipy import optimize 

15 

16DEBUG = False 

17 

18 

19# based on scipy.stats.distributions._ppf_single_call 

20def brentq_expanding(func, low=None, upp=None, args=(), xtol=1e-5, 

21 start_low=None, start_upp=None, increasing=None, 

22 max_it=100, maxiter_bq=100, factor=10, 

23 full_output=False): 

24 '''find the root of a function in one variable by expanding and brentq 

25 

26 Assumes function ``func`` is monotonic. 

27 

28 Parameters 

29 ---------- 

30 func : callable 

31 function for which we find the root ``x`` such that ``func(x) = 0`` 

32 low : float or None 

33 lower bound for brentq 

34 upp : float or None 

35 upper bound for brentq 

36 args : tuple 

37 optional additional arguments for ``func`` 

38 xtol : float 

39 parameter x tolerance given to brentq 

40 start_low : float (positive) or None 

41 starting bound for expansion with increasing ``x``. It needs to be 

42 positive. If None, then it is set to 1. 

43 start_upp : float (negative) or None 

44 starting bound for expansion with decreasing ``x``. It needs to be 

45 negative. If None, then it is set to -1. 

46 increasing : bool or None 

47 If None, then the function is evaluated at the initial bounds to 

48 determine wether the function is increasing or not. If increasing is 

49 True (False), then it is assumed that the function is monotonically 

50 increasing (decreasing). 

51 max_it : int 

52 maximum number of expansion steps. 

53 maxiter_bq : int 

54 maximum number of iterations of brentq. 

55 factor : float 

56 expansion factor for step of shifting the bounds interval, default is 

57 10. 

58 full_output : bool, optional 

59 If full_output is False, the root is returned. If full_output is True, 

60 the return value is (x, r), where x is the root, and r is a 

61 RootResults object. 

62 

63 

64 Returns 

65 ------- 

66 x : float 

67 root of the function, value at which ``func(x) = 0``. 

68 info : RootResult (optional) 

69 returned if ``full_output`` is True. 

70 attributes: 

71 

72 - start_bounds : starting bounds for expansion stage 

73 - brentq_bounds : bounds used with ``brentq`` 

74 - iterations_expand : number of iterations in expansion stage 

75 - converged : True if brentq converged. 

76 - flag : return status, 'converged' if brentq converged 

77 - function_calls : number of function calls by ``brentq`` 

78 - iterations : number of iterations in ``brentq`` 

79 

80 

81 Notes 

82 ----- 

83 If increasing is None, then whether the function is monotonically 

84 increasing or decreasing is inferred from evaluating the function at the 

85 initial bounds. This can fail if there is numerically no variation in the 

86 data in this range. In this case, using different starting bounds or 

87 directly specifying ``increasing`` can make it possible to move the 

88 expansion in the right direction. 

89 

90 If 

91 

92 ''' 

93 # TODO: rtol is missing, what does it do? 

94 left, right = low, upp # alias 

95 

96 # start_upp first because of possible sl = -1 > upp 

97 if upp is not None: 

98 su = upp 

99 elif start_upp is not None: 

100 if start_upp < 0: 

101 raise ValueError('start_upp needs to be positive') 

102 su = start_upp 

103 else: 

104 su = 1. 

105 

106 if low is not None: 

107 sl = low 

108 elif start_low is not None: 

109 if start_low > 0: 

110 raise ValueError('start_low needs to be negative') 

111 sl = start_low 

112 else: 

113 sl = min(-1., su - 1.) 

114 

115 # need sl < su 

116 if upp is None: 

117 su = max(su, sl + 1.) 

118 

119 # increasing or not ? 

120 if ((low is None) or (upp is None)) and increasing is None: 

121 assert sl < su # check during development 

122 f_low = func(sl, *args) 

123 f_upp = func(su, *args) 

124 

125 # special case for F-distribution (symmetric around zero for effect 

126 # size) 

127 # chisquare also takes an indefinite time (did not wait see if it 

128 # returns) 

129 if np.max(np.abs(f_upp - f_low)) < 1e-15 and sl == -1 and su == 1: 

130 sl = 1e-8 

131 f_low = func(sl, *args) 

132 increasing = (f_low < f_upp) 

133 if DEBUG: 

134 print('symm', sl, su, f_low, f_upp) 

135 

136 # possibly func returns nan 

137 delta = su - sl 

138 if np.isnan(f_low): 

139 # try just 3 points to find ``increasing`` 

140 # do not change sl because brentq can handle one nan bound 

141 for fraction in [0.25, 0.5, 0.75]: 

142 sl_ = sl + fraction * delta 

143 f_low = func(sl_, *args) 

144 if not np.isnan(f_low): 

145 break 

146 else: 

147 raise ValueError('could not determine whether function is ' + 

148 'increasing based on starting interval.' + 

149 '\nspecify increasing or change starting ' + 

150 'bounds') 

151 if np.isnan(f_upp): 

152 for fraction in [0.25, 0.5, 0.75]: 

153 su_ = su + fraction * delta 

154 f_upp = func(su_, *args) 

155 if not np.isnan(f_upp): 

156 break 

157 else: 

158 raise ValueError('could not determine whether function is' + 

159 'increasing based on starting interval.' + 

160 '\nspecify increasing or change starting ' + 

161 'bounds') 

162 

163 increasing = (f_low < f_upp) 

164 

165 if DEBUG: 

166 print('low, upp', low, upp, func(sl, *args), func(su, *args)) 

167 print('increasing', increasing) 

168 print('sl, su', sl, su) 

169 

170 if not increasing: 

171 sl, su = su, sl 

172 left, right = right, left 

173 

174 n_it = 0 

175 if left is None and sl != 0: 

176 left = sl 

177 while func(left, *args) > 0: 

178 # condition is also false if func returns nan 

179 right = left 

180 left *= factor 

181 if n_it >= max_it: 

182 break 

183 n_it += 1 

184 # left is now such that func(left) < q 

185 if right is None and su != 0: 

186 right = su 

187 while func(right, *args) < 0: 

188 left = right 

189 right *= factor 

190 if n_it >= max_it: 

191 break 

192 n_it += 1 

193 # right is now such that func(right) > q 

194 

195 if n_it >= max_it: 

196 # print('Warning: max_it reached') 

197 # TODO: use Warnings, Note: brentq might still work even with max_it 

198 f_low = func(sl, *args) 

199 f_upp = func(su, *args) 

200 if np.isnan(f_low) and np.isnan(f_upp): 

201 # can we still get here? 

202 raise ValueError('max_it reached' + 

203 '\nthe function values at boths bounds are NaN' + 

204 '\nchange the starting bounds, set bounds' + 

205 'or increase max_it') 

206 

207 res = optimize.brentq(func, left, right, args=args, 

208 xtol=xtol, maxiter=maxiter_bq, 

209 full_output=full_output) 

210 if full_output: 

211 val = res[0] 

212 info = res[1] 

213 info.iterations_expand = n_it 

214 info.start_bounds = (sl, su) 

215 info.brentq_bounds = (left, right) 

216 info.increasing = increasing 

217 return val, info 

218 else: 

219 return res