SloppyCell.Residuals

1 import scipy 2

3 -class Residual:

4 - def __init__(self, key):

5 self.key = key

6

7 - def GetKey(self):

8 """ 9 A unique (and hashable) identifier for this residual. 10 """ 11 return key

12

13 - def GetValue(self, predictions, internalVars, params):

14 """ 15 The value of the residual give a current set of parameters and 16 resulting predictions. 17 """ 18 raise NotImplementedError

19

20 - def GetRequiredVarsByCalc(self):

21 """ 22 The variables that need to be calculated to evaluate this residual. 23 24 Should return a nested dictionary of the form: 25 {calc name: {dependent var name: [independent var values]}} 26 """ 27 raise NotImplementedError

28

29 - def dp(self, predictions, internalVars, params):

30 """ 31 Partial derivative of the residual with respect to any parameters. 32 33 Should return a dictionary of the form: 34 {parameter name: derivative} 35 """ 36 raise NotImplementedError

37

38 - def dy(self, predictions, internalVars, params):

39 """ 40 Partial derivative of the residual with respect to any calculated 41 variables. 42 43 Should return a dictionary of the form: 44 {calculation name: {variable name: {x value: deriv}}} 45 """ 46 raise NotImplementedError

47

48 - def dintVars(self, predictions, internalVars, params):

49 """ 50 Partial derivative of the residual with respect to any internal 51 variables. 52 53 Should return a dictionary of the form: 54 {type of internal var: {expt name: {variable name: derivative}}} 55 56 XXX: This form of nesting is only appropriate for scale factors. 57 """ 58 raise NotImplementedError

59

60 - def Dp(self, predictions, senspredictions, internalVars, internalVarsDerivs, 61 params):

62 """ 63 Total derivatives with respect to all parameters of the residual. 64 65 Should return a list with the derivatives in the same order as params. 66 67 XXX: This only works with internvalVars that are indexed like scale 68 factors. 69 """ 70 derivs_wrt_p = [] 71 for pname in params.keys(): 72 deriv = 0 73 74 # This first term is dres/dy * dy/dp 75 dres_dy = self.dy(predictions, internalVars, params) 76 for calcKey in dres_dy.keys(): 77 for yKey in dres_dy[calcKey].keys(): 78 for xVal in dres_dy[calcKey][yKey].keys(): 79 dres_dy_this = dres_dy[calcKey][yKey][xVal] 80 # We default to 0 if parameter not in senspredictions. 81 # (It may not be involved in a given calculation.) 82 dy_dp_this = senspredictions[calcKey][yKey][xVal].get(pname, 0) 83 deriv += dres_dy_this * dy_dp_this 84 85 # This term is dres/dintVar * dintVar/dp 86 dres_dintVars = self.dintVars(predictions, internalVars, params) 87 for intVar_type in dres_dintVars.keys(): 88 for exptKey in dres_dintVars[intVar_type].keys(): 89 for yKey in dres_dintVars[intVar_type][exptKey].keys(): 90 dres_dintVar_this = dres_dintVars[intVar_type][exptKey][yKey] 91 dintVar_dp_this = internalVarsDerivs[intVar_type][exptKey][yKey][pname] 92 deriv += dres_dintVar_this * dintVar_dp_this 93 94 # Finally dres/dp 95 dres_dp = self.dp(predictions, internalVars, params) 96 deriv += dres_dp.get(pname, 0) 97 98 derivs_wrt_p.append(deriv) 99 100 return derivs_wrt_p

101

102 -class ScaledErrorInFit(Residual):

103 - def __init__(self, key, depVarKey, calcKey, indVarValue, depVarMeasurement, 104 depVarSigma, exptKey):

105 Residual.__init__(self, key) 106 self.yKey = depVarKey 107 self.calcKey = calcKey 108 self.xVal = indVarValue 109 self.yMeas = depVarMeasurement 110 self.ySigma = depVarSigma 111 self.exptKey = exptKey

112

113 - def GetRequiredVarsByCalc(self):

114 return {self.calcKey: {self.yKey: [self.xVal]}}

115

116 - def GetValue(self, predictions, internalVars, params):

117 scale_factor = internalVars['scaleFactors'][self.exptKey][self.yKey] 118 raw_pred_val = predictions[self.calcKey][self.yKey][self.xVal] 119 return (scale_factor * raw_pred_val - self.yMeas)/self.ySigma

120

121 - def dp(self, predictions, internalVars, params):

122 return {}

123

124 - def dy(self, predictions, internalVars, params):

125 scale_factor = internalVars['scaleFactors'][self.exptKey][self.yKey] 126 deriv = scale_factor / self.ySigma 127 return {self.calcKey: {self.yKey: {self.xVal: deriv}}}

128

129 - def dintVars(self, predictions, internalVars, params):

130 raw_pred_val = predictions[self.calcKey][self.yKey][self.xVal] 131 deriv = raw_pred_val / self.ySigma 132 return {'scaleFactors': {self.exptKey: {self.yKey: deriv}}}

133

134 -class PriorInLog(Residual):

135 - def __init__(self, key, pKey, logPVal, sigmaLogPVal):

136 Residual.__init__(self, key) 137 self.pKey = pKey 138 self.logPVal = logPVal 139 self.sigmaLogPVal = sigmaLogPVal

140

141 - def GetValue(self, predictions, internalVars, params):

142 return (scipy.log(params.get(self.pKey)) - self.logPVal) / self.sigmaLogPVal

143

144 - def dp(self, predictions, internalVars, params):

145 return {self.pKey: 1./(params.get(self.pKey) * self.sigmaLogPVal)}

146

147 - def dy(self, predictions, internalVars, params):

148 return {}

149

150 - def dintVars(self, predictions, internalVars, params):

151 return {}

152

153 -class PeriodCheckResidual(Residual):

154 - def __init__(self, key, calcKey, depVarKey, indVarValue, depVarMeasurement, 155 depVarSigma):

156 Residual.__init__(self, key) 157 self.cKey = calcKey 158 self.yKey = depVarKey 159 self.xVal = indVarValue 160 self.yMeas = depVarMeasurement 161 self.ySigma = depVarSigma

162

163 - def GetRequiredVarsByCalc(self):

164 # Need the points between xVal and xVal + 2*periods 165 return {self.cKey: {self.yKey: [self.xVal, self.xVal+2.0*self.yMeas]}}

166

167 - def GetValue(self, predictions, internalVars, params):

168 # Find the period 169 traj = predictions[self.cKey][self.yKey] 170 times = traj.keys() 171 times.sort() 172 173 maximums=[] 174 for index in range(1,len(times)-1): 175 if times[index]<self.xVal\ 176 or (self.xVal+2.0*self.yMeas)<times[index]: continue 177 t1,t2,t3 = times[index-1:index+2] 178 y1,y2,y3 = traj[t1], traj[t2], traj[t3] 179 if y1<y2 and y3<y2: 180 # Use a quadratic approximation to find the maximum 181 (a,b,c)=scipy.dot(scipy.linalg.inv([[t1**2,t1,1], 182 [t2**2,t2,1], 183 [t3**2,t3,1]]),[y1,y2,y3]) 184 maximums.append(-b/2/a) 185 186 if len(maximums)<2: 187 theoryVal = 2*self.yMeas 188 else: 189 theoryVal = maximums[1] - maximums[0] 190 191 return (theoryVal - self.yMeas)/self.ySigma

192

193 -class AmplitudeCheckResidual(Residual):

194 - def __init__(self, key, calcKey, depVarKey, indVarValue0, indVarValue1, period, 195 depVarSigma, exptKey):

196 Residual.__init__(self, key) 197 self.cKey = calcKey 198 self.yKey = depVarKey 199 self.xVal = indVarValue0 200 self.xTestVal = indVarValue1 201 self.period = period 202 self.ySigma = depVarSigma 203 self.exptKey = exptKey

204

205 - def GetRequiredVarsByCalc(self):

206 # Need the points between xVal and xVal + 2*periods 207 return {self.cKey: {self.yKey: [self.xVal, self.xVal+self.period, self.xTestVal, self.xTestVal+self.period]}}

208

209 - def GetValue(self, predictions, internalVars, params):

210 times = predictions[self.cKey][self.yKey].keys() 211 if self.exptKey in internalVars['scaleFactors'].keys() \ 212 and self.yKey in internalVars['scaleFactors'][self.exptKey].keys(): 213 scale_factor = internalVars['scaleFactors'][self.exptKey][self.yKey] 214 else: 215 scale_factor = 1.0 216 217 # Get the indices of the points to use and integrate the areas 218 times = predictions[self.cKey][self.yKey].keys() 219 times.sort() 220 startIndex,endStartIndex = times.index(self.xVal), times.index(self.xVal+self.period) 221 testIndex,endTestIndex = times.index(self.xTestVal), times.index(self.xTestVal+self.period) 222 223 x,y=[],[] 224 for t in times[startIndex:endStartIndex+1]: 225 x.append(t) 226 y.append(scale_factor*predictions[self.cKey][self.yKey][t]) 227 measVal = scipy.integrate.simps(y,x) 228 229 x,y=[],[] 230 for t in times[testIndex:endTestIndex+1]: 231 x.append(t) 232 y.append(scale_factor*predictions[self.cKey][self.yKey][t]) 233 theoryVal = scipy.integrate.simps(y,x) 234 235 return (theoryVal - measVal)/self.ySigma

236

237 -class IntegralDataResidual(Residual):

238 - def __init__(self, name, var, exptKey, calc, traj, uncert_traj, interval):

239 self.name = name 240 self.var = var 241 self.exptKey = exptKey 242 self.calc = calc 243 self.traj = traj 244 self.uncert_traj = uncert_traj 245 self.interval = interval

246

247 - def GetValue(self, predictions, internalVars, params):

248 sf = internalVars['scaleFactors'][self.exptKey][self.var] 249 data_traj = self.traj 250 uncert_traj = self.uncert_traj 251 theory_traj = predictions[self.calc]['full trajectory'] 252 var = self.var 253 def integrand(t): 254 theory = theory_traj.evaluate_interpolated_traj(var, t) 255 data = data_traj.evaluate_interpolated_traj(var, t) 256 uncert = uncert_traj.evaluate_interpolated_traj(var, t) 257 return (sf*theory - data)**2/uncert**2

258 val, error = scipy.integrate.quad(integrand, 259 self.interval[0], self.interval[1], 260 limit = int(1e5)) 261 T = self.interval[1] - self.interval[0] 262 return scipy.sqrt(val/T)

263

264 -class ScaledExtremum(Residual):

265 - def __init__(self, key, var, calcKey, val, 266 sigma, exptKey, minTime=None, maxTime=None, type=None):

267 Residual.__init__(self, key) 268 self.yKey = var 269 self.calcKey = calcKey 270 self.yMeas = val 271 self.ySigma = sigma 272 self.exptKey = exptKey 273 self.minTime,self.maxTime = minTime,maxTime 274 self.last_time_result = None 275 self.type = type

276

277 - def GetRequiredVarsByCalc(self):

278 if self.type == 'max': 279 var = self.yKey + '_maximum' 280 elif self.type == 'min': 281 var = self.yKey + '_minimum' 282 283 return {self.calcKey: {var: [self.minTime, self.maxTime]}}

284

285 - def GetValue(self, predictions, internalVars, params):

286 scale_factor = internalVars['scaleFactors'][self.exptKey][self.yKey] 287 # predictions entry is time,value 288 # We store the last time result for use in plotting. 289 if self.type == 'max': 290 var = self.yKey + '_maximum' 291 elif self.type == 'min': 292 var = self.yKey + '_minimum' 293 294 self.last_time_result, raw_pred_val = \ 295 predictions[self.calcKey][var][self.minTime,self.maxTime] 296 return (scale_factor * raw_pred_val - self.yMeas)/self.ySigma

297

298 - def dp(self, predictions, internalVars, params):

299 return {}

300

Source Code for Module SloppyCell.Residuals