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

2import numpy as np 

3from scipy.optimize import minimize 

4import theano 

5from theano import gof 

6import theano.tensor as tt 

7from theano.configparser import change_flags 

8 

9try: 

10 import healpy as hp 

11except ImportError: 

12 hp = None 

13 

14 

15__all__ = ["minimizeOp", "LDPhysicalOp"] 

16 

17 

18class minimizeOp(tt.Op): 

19 """Find the global minimum of the map intensity. 

20 

21 Returns the tuple `(lat, lon, I)`. 

22 

23 .. note:: 

24 This op is not very optimized. The idea here is to 

25 do a coarse grid search, find the minimum, then run 

26 a quick gradient descent to refine the result. 

27 """ 

28 

29 def __init__(self, intensity, P, ydeg, udeg, fdeg): 

30 self.intensity = intensity 

31 self.P = P 

32 self.ydeg = ydeg 

33 self.udeg = udeg 

34 self.fdeg = fdeg 

35 self._do_setup = True 

36 self.oversample = 1 

37 self.ntries = 1 

38 self.result = None 

39 

40 def setup(self, oversample=1, ntries=1): 

41 self.ntries = ntries 

42 

43 # Don't setup unless the user actually calls this function, 

44 # since there's quite a bit of overhead 

45 if self._do_setup or (oversample != self.oversample): 

46 

47 self.oversample = oversample 

48 

49 # Create the grid using healpy if available 

50 # Require at least `oversample * l ** 2 points` 

51 s = np.random.RandomState(0) 

52 if hp is None: 

53 npts = oversample * self.ydeg ** 2 

54 self.lat_grid = ( 

55 np.arccos(2 * s.rand(npts) - 1) * 180.0 / np.pi - 90.0 

56 ) 

57 self.lon_grid = (s.rand(npts) - 0.5) * 360 

58 else: 

59 nside = 1 

60 while hp.nside2npix(nside) < oversample * self.ydeg ** 2: 

61 nside += 1 

62 theta, phi = hp.pix2ang( 

63 nside=nside, ipix=range(hp.nside2npix(nside)) 

64 ) 

65 self.lat_grid = 0.5 * np.pi - theta 

66 self.lon_grid = phi - np.pi 

67 # Add a little noise for stability 

68 self.lat_grid += 1e-4 * s.randn(len(self.lat_grid)) 

69 self.lon_grid += 1e-4 * s.randn(len(self.lon_grid)) 

70 self.P_grid = self.P( 

71 tt.as_tensor_variable(self.lat_grid), 

72 tt.as_tensor_variable(self.lon_grid), 

73 ).eval() 

74 

75 # Set up the cost & grad function for the nonlinear solver 

76 u0 = np.zeros(self.udeg + 1) 

77 u0[0] = -1.0 

78 u0 = tt.as_tensor_variable(u0) 

79 f0 = np.zeros((self.fdeg + 1) ** 2) 

80 f0[0] = np.pi 

81 f0 = tt.as_tensor_variable(f0) 

82 latlon = tt.dvector() 

83 y = tt.dvector() 

84 with change_flags(compute_test_value="off"): 

85 self.I = theano.function( 

86 [latlon, y], 

87 [ 

88 self.intensity( 

89 latlon[0], latlon[1], y, u0, f0, 0.0, 0.0 

90 )[0], 

91 *theano.grad( 

92 self.intensity( 

93 latlon[0], latlon[1], y, u0, f0, 0.0, 0.0 

94 )[0], 

95 [latlon], 

96 ), 

97 ], 

98 ) 

99 self._do_setup = False 

100 

101 def make_node(self, *inputs): 

102 inputs = [tt.as_tensor_variable(i) for i in inputs] 

103 outputs = [ 

104 tt.TensorType(inputs[0].dtype, ())(), 

105 tt.TensorType(inputs[0].dtype, ())(), 

106 tt.TensorType(inputs[0].dtype, ())(), 

107 ] 

108 return gof.Apply(self, inputs, outputs) 

109 

110 def infer_shape(self, node, shapes): 

111 return [(), (), ()] 

112 

113 def perform(self, node, inputs, outputs): 

114 

115 assert self._do_setup is False, "Must run `setup()` first." 

116 

117 y = inputs[0] 

118 

119 # Initial guess 

120 I_grid = np.dot(self.P_grid, y) 

121 outputs[2][0] = np.inf 

122 

123 for n in range(self.ntries): 

124 

125 ind = np.argmin(I_grid) 

126 x0 = np.array([self.lat_grid[ind], self.lon_grid[ind]]) 

127 

128 # Minimize 

129 result = minimize(self.I, x0, args=(y,), jac=True) 

130 

131 if result.fun < outputs[2][0]: 

132 outputs[0][0] = result.x[0] 

133 outputs[1][0] = result.x[1] 

134 outputs[2][0] = result.fun 

135 

136 # Prepare for next iteration 

137 I_grid[ind] = np.inf 

138 

139 # Save 

140 self.result = result 

141 

142 

143class LDPhysicalOp(tt.Op): 

144 """ 

145 Check whether a limb darkening profile is physical using Sturm's theorem. 

146 

147 """ 

148 

149 def __init__(self, nroots): 

150 self.nroots = nroots 

151 

152 def make_node(self, *inputs): 

153 inputs = [tt.as_tensor_variable(i) for i in inputs] 

154 outputs = [tt.bscalar()] 

155 return gof.Apply(self, inputs, outputs) 

156 

157 def infer_shape(self, node, shapes): 

158 return [()] 

159 

160 def perform(self, node, inputs, outputs): 

161 u = inputs[0] 

162 outputs[0][0] = 1 

163 

164 # Ensure the function is *decreasing* toward the limb 

165 if u.sum() < -1: 

166 outputs[0][0] = 0 

167 return 

168 

169 # Sturm's theorem on the intensity to ensure positivity 

170 p = u[::-1] 

171 if self.nroots(p, 0, 1) > 0: 

172 outputs[0][0] = 0 

173 return 

174 

175 # Sturm's theorem on the derivative to ensure monotonicity 

176 p = (u[1:] * np.arange(1, len(u)))[::-1] 

177 if self.nroots(p, 0, 1) > 0: 

178 outputs[0][0] = 0 

179 return