SloppyCell.Model

1 """ 2 Model class that unites theory with data. 3 """ 4 5 import logging 6 logger = logging.getLogger('Model_mod') 7 8 import copy 9 import sets 10 11 import scipy 12 13 import SloppyCell 14 import SloppyCell.Residuals as Residuals 15 import SloppyCell.Collections as Collections 16 import SloppyCell.Utility as Utility 17 import KeyedList_mod as KeyedList_mod 18 KeyedList = KeyedList_mod.KeyedList 19 20 _double_epsilon_ = scipy.finfo(scipy.float_).eps 21

22 -class Model:

23 """ 24 A Model object connects a set of experimental data with the objects used to 25 model that data. 26 27 Most importantly, a Model can calculate a cost for a given set of 28 parameters, characterizing how well those parameters fit the data contained 29 within the model. 30 """ 31 32 imag_cutoff = 1e-8 33

34 - def __init__(self, expts, calcs):

35 """ 36 expts A sequence of Experiments to be fit to. 37 calcs A sequence of calculation objects referred to by the 38 Experiments. 39 """ 40 41 self.calcVals = {} 42 self.calcSensitivityVals = {} 43 self.internalVars = {} 44 self.internalVarsDerivs = {} 45 self.residuals = KeyedList() 46 47 if isinstance(expts, list): 48 expts = Collections.ExperimentCollection(expts) 49 elif isinstance(expts, dict): 50 expts = Collections.ExperimentCollection(expts.values()) 51 self.SetExperimentCollection(expts) 52 53 if isinstance(calcs, list): 54 calcs = Collections.CalculationCollection(calcs) 55 elif isinstance(calcs, dict): 56 calcs = Collections.CalculationCollection(calcs.values()) 57 self.SetCalculationCollection(calcs) 58 59 self.observers = KeyedList()

60

61 - def compile(self):

62 """ 63 Compile all the calculations contained within the Model. 64 """ 65 for calc in self.get_calcs().values(): 66 calc.compile()

67

68 - def copy(self):

69 return copy.deepcopy(self)

70

71 - def __repr__(self):

72 return 'CalculationCollection(%s)' % repr(self.items())

73

74 - def get_params(self):

75 """ 76 Return a copy of the current model parameters 77 """ 78 return self.calcColl.GetParameters()

79

80 - def get_ICs(self):

81 """ 82 Get the initial conditions currently present in a model 83 for dynamic variables that are not assigned variables. 84 85 Outputs: 86 KeyedList with keys (calcName,varName) --> initialValue 87 """ 88 ics=KeyedList() 89 for calcName, calc in self.calcColl.items(): 90 for varName in calc.dynamicVars.keys(): 91 if varName in calc.assignedVars.keys(): continue 92 ics.set( (calcName,varName), calc.get_var_ic(varName)) 93 return ics

94

95 - def set_ICs(self, ics):

96 """ 97 Sets the initial conditions into the model. Uses the input 98 format defined by 'getICs'. 99 100 Inputs: 101 ics -- Initial conditions to set in KeyedList form: 102 keys: (calcName, varName) --> intialValue 103 104 Outputs: 105 None 106 """ 107 for (calcName, varName), initialValue in ics.items(): 108 self.calcColl.get(calcName).set_var_ic(varName, initialValue)

109

110 - def _evaluate(self, params, T=1):

111 """ 112 Evaluate the cost for the model, returning the intermediate residuals, 113 and chi-squared. 114 115 (Summing up the residuals is a negligible amount of work. This 116 arrangment makes notification of observers much simpler.) 117 """ 118 self.params.update(params) 119 self.CalculateForAllDataPoints(params) 120 self.ComputeInternalVariables(T) 121 122 resvals = [res.GetValue(self.calcVals, self.internalVars, self.params) 123 for res in self.residuals.values()] 124 125 # Occasionally it's useful to use residuals with a sqrt(-1) in them, 126 # to get negative squares. Then, however, we might get small imaginary 127 # parts in our results, which this shaves off. 128 chisq = scipy.real_if_close(scipy.sum(scipy.asarray(resvals)**2), 129 tol=self.imag_cutoff) 130 if scipy.isnan(chisq): 131 logger.warn('Chi^2 is NaN, converting to Infinity.') 132 chisq = scipy.inf 133 cost = 0.5 * chisq 134 135 entropy = 0 136 for expt, sf_ents in self.internalVars['scaleFactor_entropies'].items(): 137 for group, ent in sf_ents.items(): 138 entropy += ent 139 140 self._notify(event = 'evaluation', 141 resvals = resvals, 142 chisq = chisq, 143 cost = cost, 144 free_energy = cost-T*entropy, 145 entropy = entropy, 146 params = self.params) 147 148 return resvals, chisq, cost, entropy

149

150 - def res(self, params):

151 """ 152 Return the residual values of the model fit given a set of parameters 153 """ 154 155 return self._evaluate(params)[0]

156

157 - def res_log_params(self, log_params):

158 """ 159 Return the residual values given the logarithm of the parameters 160 """ 161 return self.res(scipy.exp(log_params))

162

163 - def res_dict(self, params):

164 """ 165 Return the residual values of the model fit given a set of parameters 166 in dictionary form. 167 """ 168 return dict(zip(self.residuals.keys(), self.res(params)))

169

170 - def chisq(self, params):

171 """ 172 Return the sum of the squares of the residuals for the model 173 """ 174 return self._evaluate(params)[1]

175

176 - def redchisq(self, params):

177 """ 178 Return chi-squared divided by the number of degrees of freedom 179 180 Question: Are priors to be included in the N data points? 181 How do scale factors change the number of d.o.f.? 182 """ 183 return self.chisq(params)/(len(self.residuals) - len(self.params))

184

185 - def cost(self, params):

186 """ 187 Return the cost (1/2 chisq) of the model 188 """ 189 return self._evaluate(params)[2]

190

191 - def cost_log_params(self, log_params):

192 """ 193 Return the cost given the logarithm of the input parameters 194 """ 195 return self.cost(scipy.exp(log_params))

196

197 - def free_energy(self, params, T):

198 temp, temp, c, entropy = self._evaluate(params, T=T) 199 return c - T * entropy

200

201 - def _notify(self, **args):

202 """ 203 Call all observers with the given arguments. 204 """ 205 for obs in self.observers: 206 obs(**args)

207

208 - def attach_observer(self, obs_key, observer):

209 """ 210 Add an observer to be notified by this Model. 211 """ 212 self.observers.set(obs_key, observer)

213

214 - def detach_observer(self, obs_key):

215 """ 216 Remove an observer from the Model. 217 """ 218 self.observers.remove_by_key(obs_key)

219

220 - def get_observers(self):

221 """ 222 Return the KeyedList of observers for this model. 223 """ 224 return self.observers

225

226 - def reset_observers(self):

227 """ 228 Call reset() for all attached observers. 229 """ 230 for obs in self.observers: 231 if hasattr(obs, 'reset'): 232 obs.reset()

233 234 resDict = res_dict 235 # ... 236 237

238 - def AddResidual(self, res):

239 self.residuals.setByKey(res.key, res)

240

241 - def Force(self, params, epsf, relativeScale=False, stepSizeCutoff=None):

242 """ 243 Force(parameters, epsilon factor) -> list 244 245 Returns a list containing the numerical gradient of the cost with 246 respect to each parameter (in the parameter order of the 247 CalculationCollection). Each element of the gradient is: 248 cost(param + eps) - cost(param - eps)/(2 * eps). 249 If relativeScale is False then epsf is the stepsize used (it should 250 already be multiplied by typicalValues before Jacobian is called) 251 If relativeScale is True then epsf is multiplied by params. 252 The two previous statements hold for both scalar and vector valued 253 epsf. 254 """ 255 256 force = [] 257 params = scipy.array(params) 258 259 if stepSizeCutoff==None: 260 stepSizeCutoff = scipy.sqrt(_double_epsilon_) 261 262 if relativeScale is True: 263 eps = epsf * abs(params) 264 else: 265 eps = epsf * scipy.ones(len(params),scipy.float_) 266 267 for i in range(0,len(eps)): 268 if eps[i] < stepSizeCutoff: 269 eps[i] = stepSizeCutoff 270 271 for index, param in enumerate(params): 272 paramsPlus = params.copy() 273 paramsPlus[index] = param + eps[index] 274 costPlus = self.cost(paramsPlus) 275 276 paramsMinus = params.copy() 277 paramsMinus[index] = param - eps[index] 278 costMinus = self.cost(paramsMinus) 279 280 force.append((costPlus-costMinus)/(2.0*eps[index])) 281 282 return force

283

284 - def gradient_sens(self, params):

285 """ 286 Return the gradient of the cost, d_cost/d_param as a KeyedList. 287 288 This method uses sensitivity integration, so it only applies to 289 ReactionNetworks. 290 """ 291 self.params.update(params) 292 293 # The cost is 0.5 * sum(res**2), 294 # so the gradient is sum(res * dres_dp) 295 296 jac_dict = self.jacobian_sens(params) 297 res_dict = self.res_dict(params) 298 299 force = scipy.zeros(len(params), scipy.float_) 300 for res_key, res_val in res_dict.items(): 301 res_derivs = jac_dict.get(res_key) 302 force += res_val * scipy.asarray(res_derivs) 303 304 gradient = self.params.copy() 305 gradient.update(force) 306 307 return gradient

308

309 - def gradient_log_params_sens(self, log_params):

310 """ 311 Return the gradient of the cost wrt log parameters, d_cost/d_log_param 312 as a KeyedList. 313 314 This method uses sensitivity integration, so it only applies to 315 ReactionNetworks. 316 """ 317 # We just need to multiply dcost_dp by p. 318 params = scipy.exp(log_params) 319 gradient = self.gradient_sens(params) 320 gradient_log = gradient.copy() 321 gradient_log.update(scipy.asarray(gradient) * scipy.asarray(params)) 322 323 return gradient_log

324

325 - def CalculateForAllDataPoints(self, params):

326 """ 327 CalculateForAllDataPoints(parameters) -> dictionary 328 329 Gets a dictionary of measured independent variables indexed by 330 calculation from the ExperimentCollection and passes it to the 331 CalculationCollection. The returned dictionary is of the form: 332 dictionary[experiment][calculation][dependent variable] 333 [independent variabled] -> calculated value. 334 """ 335 self.params.update(params) 336 337 varsByCalc = self.GetExperimentCollection().GetVarsByCalc() 338 self.calcVals = self.GetCalculationCollection().Calculate(varsByCalc, 339 params) 340 return self.calcVals

341

342 - def CalculateSensitivitiesForAllDataPoints(self, params):

343 """ 344 CalculateSensitivitiesForAllDataPoints(parameters) -> dictionary 345 346 Gets a dictionary of measured independent variables indexed by 347 calculation from the ExperimentCollection and passes it to the 348 CalculationCollection. The returned dictionary is of the form: 349 dictionary[experiment][calculation][dependent variable] 350 [independent variabled][parameter] -> sensitivity. 351 """ 352 varsByCalc = self.GetExperimentCollection().GetVarsByCalc() 353 self.calcVals, self.calcSensitivityVals =\ 354 self.GetCalculationCollection().CalculateSensitivity(varsByCalc, 355 params) 356 return self.calcSensitivityVals

357

358 - def ComputeInternalVariables(self, T=1):

359 sf, sf_ents = self.compute_scale_factors(T) 360 self.internalVars['scaleFactors'] = sf 361 self.internalVars['scaleFactor_entropies'] = sf_ents

362

363 - def compute_scale_factors(self, T):

364 """ 365 Compute the scale factors for the current parameters and return a dict. 366 367 The dictionary is of the form dict[exptId][varId] = scale_factor 368 """ 369 scale_factors = {} 370 scale_factor_entropies = {} 371 for exptId, expt in self.GetExperimentCollection().items(): 372 scale_factors[exptId], scale_factor_entropies[exptId] =\ 373 self._compute_sf_and_sfent_for_expt(expt, T) 374 375 return scale_factors, scale_factor_entropies

376

377 - def _compute_sf_and_sfent_for_expt(self, expt, T):

378 # Compute the scale factors for a given experiment 379 scale_factors = {} 380 scale_factor_entropies = {} 381 382 exptData = expt.GetData() 383 expt_integral_data = expt.GetIntegralDataSets() 384 fixed_sf = expt.get_fixed_sf() 385 386 sf_groups = expt.get_sf_groups() 387 388 for group in sf_groups: 389 # Do any of the variables in this group have fixed scale factors? 390 fixed = sets.Set(group).intersection(sets.Set(fixed_sf.keys())) 391 fixedAt = sets.Set([fixed_sf[var] for var in fixed]) 392 393 # We'll need to index the scale factor entropies on the *group* 394 # that shares a scale factor, since we only have one entropy per 395 # shared scale factor. So we need to index on the group of 396 # variables. We sort the group and make it hashable to avoid any 397 # double-counting. 398 hash_group = expt._hashable_group(group) 399 if len(fixedAt) == 1: 400 value = fixedAt.pop() 401 for var in group: 402 scale_factors[var] = value 403 scale_factor_entropies[hash_group] = 0 404 continue 405 elif len(fixedAt) > 1: 406 raise ValueError('Shared scale factors fixed at ' 407 'inconsistent values in experiment ' 408 '%s!' % expt.GetName()) 409 410 # Finally, compute the scale factor for this group 411 theoryDotData, theoryDotTheory = 0, 0 412 # For discrete data 413 for calc in exptData: 414 # Pull out the vars we have measured for this calculation 415 for var in sets.Set(group).intersection(sets.Set(exptData[calc].keys())): 416 for indVar, (data, error) in exptData[calc][var].items(): 417 theory = self.calcVals[calc][var][indVar] 418 theoryDotData += (theory * data) / error**2 419 theoryDotTheory += theory**2 / error**2 420 421 # Now for integral data 422 for dataset in expt_integral_data: 423 calc = dataset['calcKey'] 424 theory_traj = self.calcVals[calc]['full trajectory'] 425 data_traj = dataset['trajectory'] 426 uncert_traj = dataset['uncert_traj'] 427 interval = dataset['interval'] 428 T = interval[1] - interval[0] 429 for var in group.intersection(sets.Set(dataset['vars'])): 430 TheorDotT = self._integral_theorytheory(var, theory_traj, 431 uncert_traj, 432 interval) 433 theoryDotTheory += TheorDotT/T 434 TheorDotD = self._integral_theorydata(var, theory_traj, 435 data_traj, 436 uncert_traj, 437 interval) 438 theoryDotData += TheorDotD/T 439 440 # Now for the extrema data 441 for ds in expt.scaled_extrema_data: 442 calc = ds['calcKey'] 443 if ds['type'] == 'max': 444 var = ds['var'] + '_maximum' 445 elif ds['type'] == 'min': 446 var = ds['var'] + '_minimum' 447 data, error = ds['val'], ds['sigma'] 448 theory = self.calcVals[calc][var]\ 449 [ds['minTime'],ds['maxTime']][1] 450 theoryDotData += (theory * data) / error**2 451 theoryDotTheory += theory**2 / error**2 452 453 for var in group: 454 if theoryDotTheory != 0: 455 scale_factors[var] = theoryDotData/theoryDotTheory 456 else: 457 scale_factors[var] = 1 458 entropy = expt.compute_sf_entropy(hash_group, theoryDotTheory, 459 theoryDotData, T) 460 scale_factor_entropies[hash_group] = entropy 461 462 return scale_factors, scale_factor_entropies

463

464 - def _integral_theorytheory(self, var, theory_traj, uncert_traj, interval):

465 def integrand(t): 466 theory = theory_traj.evaluate_interpolated_traj(var, t) 467 uncert = uncert_traj.evaluate_interpolated_traj(var, t) 468 return theory**2/uncert**2

469 val, error = scipy.integrate.quad(integrand, interval[0], interval[1], 470 limit=int(1e5)) 471 return val

472

473 - def _integral_theorydata(self, var, theory_traj, data_traj, uncert_traj, 474 interval):

475 def integrand(t): 476 theory = theory_traj.evaluate_interpolated_traj(var, t) 477 data = data_traj.evaluate_interpolated_traj(var, t) 478 uncert = uncert_traj.evaluate_interpolated_traj(var, t) 479 return theory*data/uncert**2

480 val, error = scipy.integrate.quad(integrand, interval[0], interval[1], 481 limit=int(1e5)) 482 return val 483

484 - def ComputeInternalVariableDerivs(self):

485 """ 486 compute_scale_factorsDerivs() -> dictionary 487 488 Returns the scale factor derivatives w.r.t. parameters 489 appropriate for each chemical in each 490 experiment, given the current parameters. The returned dictionary 491 is of the form: internalVarsDerivs['scaleFactors'] \ 492 = dict[experiment][chemical][parametername] -> derivative. 493 """ 494 495 self.internalVarsDerivs['scaleFactors'] = {} 496 p = self.GetCalculationCollection().GetParameters() 497 498 for exptName, expt in self.GetExperimentCollection().items(): 499 self.internalVarsDerivs['scaleFactors'][exptName] = {} 500 exptData = expt.GetData() 501 502 # Get the dependent variables measured in this experiment 503 exptDepVars = sets.Set() 504 for calc in exptData: 505 exptDepVars.union_update(sets.Set(expt.GetData()[calc].keys())) 506 507 for depVar in exptDepVars: 508 self.internalVarsDerivs['scaleFactors'][exptName][depVar] = {} 509 if depVar in expt.GetFixedScaleFactors(): 510 for pname in p.keys(): 511 self.internalVarsDerivs['scaleFactors'][exptName]\ 512 [depVar][pname] = 0.0 513 continue 514 515 theoryDotData, theoryDotTheory = 0, 0 516 for calc in exptData: 517 if depVar in exptData[calc].keys(): 518 for indVar, (data, error)\ 519 in exptData[calc][depVar].items(): 520 theory = self.calcVals[calc][depVar][indVar] 521 theoryDotData += (theory * data) / error**2 522 theoryDotTheory += theory**2 / error**2 523 524 # now get derivative of the scalefactor 525 for pname in p.keys(): 526 theorysensDotData, theorysensDotTheory = 0, 0 527 528 for calc in exptData: 529 clc = self.calcColl.get(calc) 530 if depVar in exptData[calc].keys(): 531 for indVar, (data, error)\ 532 in exptData[calc][depVar].items(): 533 theory = self.calcVals[calc][depVar][indVar] 534 535 # Default to 0 if sensitivity not calculated for 536 # that parameter (i.e. it's not in the 537 # Calculation) 538 theorysens = self.calcSensitivityVals[calc][depVar][indVar].get(pname, 0.0) 539 theorysensDotData += (theorysens * data) / error**2 540 theorysensDotTheory += (theorysens * theory) / error**2 541 542 deriv_dict = self.internalVarsDerivs['scaleFactors'][exptName][depVar] 543 try: 544 deriv_dict[pname] = theorysensDotData/theoryDotTheory\ 545 - 2*theoryDotData*theorysensDotTheory/(theoryDotTheory)**2 546 except ZeroDivisionError: 547 deriv_dict[pname] = 0 548 549 return self.internalVarsDerivs['scaleFactors']

550

551 - def jacobian_log_params_sens(self, log_params):

552 """ 553 Return a KeyedList of the derivatives of the model residuals w.r.t. 554 the lograithms of the parameters parameters. 555 556 The method uses the sensitivity integration. As such, it will only 557 work with ReactionNetworks. 558 559 The KeyedList is of the form: 560 kl.get(resId) = [dres/dlogp1, dres/dlogp2...] 561 """ 562 params = scipy.exp(log_params) 563 j = self.jacobian_sens(params) 564 j_log = j.copy() 565 j_log.update(scipy.asarray(j) * scipy.asarray(params)) 566 567 return j_log

568

569 - def jacobian_sens(self, params):

570 """ 571 Return a KeyedList of the derivatives of the model residuals w.r.t. 572 parameters. 573 574 The method uses the sensitivity integration. As such, it will only 575 work with ReactionNetworks. 576 577 The KeyedList is of the form: 578 kl[resId] = [dres/dp1, dres/dp2...] 579 """ 580 self.params.update(params) 581 582 # Calculate sensitivities 583 self.CalculateSensitivitiesForAllDataPoints(params) 584 self.ComputeInternalVariables() 585 self.ComputeInternalVariableDerivs() 586 587 # Calculate residual derivatives 588 deriv = [(resId, res.Dp(self.calcVals, self.calcSensitivityVals, 589 self.internalVars, self.internalVarsDerivs, 590 self.params)) 591 for (resId, res) in self.residuals.items()] 592 593 return KeyedList(deriv)

594

595 - def jacobian_fd(self, params, eps, 596 relativeScale=False, stepSizeCutoff=None):

597 """ 598 Return a KeyedList of the derivatives of the model residuals w.r.t. 599 parameters. 600 601 The method uses finite differences. 602 603 Inputs: 604 params -- Parameters about which to calculate the jacobian 605 eps -- Step size to take, may be vector or scalar. 606 relativeScale -- If true, the eps is taken to be the fractional 607 change in parameter to use in finite differences. 608 stepSizeCutoff -- Minimum step size to take. 609 """ 610 res = self.resDict(params) 611 612 orig_vals = scipy.array(params) 613 614 if stepSizeCutoff is None: 615 stepSizeCutoff = scipy.sqrt(_double_epsilon_) 616 617 if relativeScale: 618 eps_l = scipy.maximum(eps * abs(params), stepSizeCutoff) 619 else: 620 eps_l = scipy.maximum(eps * scipy.ones(len(params),scipy.float_), 621 stepSizeCutoff) 622 623 J = KeyedList() # will hold the result 624 for resId in res.keys(): 625 J.set(resId, []) 626 # Two-sided finite difference 627 for ii in range(len(params)): 628 params[ii] = orig_vals[ii] + eps_l[ii] 629 resPlus = self.resDict(params) 630 631 params[ii] = orig_vals[ii] - eps_l[ii] 632 resMinus = self.resDict(params) 633 634 params[ii] = orig_vals[ii] 635 636 for resId in res.keys(): 637 res_deriv = (resPlus[resId]-resMinus[resId])/(2.*eps_l[ii]) 638 J.get(resId).append(res_deriv) 639 640 # NOTE: after call to ComputeResidualsWithScaleFactors the Model's 641 # parameters get updated, must reset this: 642 self.params.update(params) 643 return J

644

645 - def GetJacobian(self,params):

646 """ 647 GetJacobian(parameters) -> dictionary 648 649 Gets a dictionary of the sensitivities at the time points of 650 the independent variables for the measured dependent variables 651 for each calculation and experiment. 652 Form: 653 dictionary[(experiment,calculation,dependent variable, 654 independent variable)] -> result 655 656 result is a vector of length number of parameters containing 657 the sensitivity at that time point, in the order of the ordered 658 parameters 659 660 """ 661 return self.jacobian_sens(params)

662

663 - def Jacobian(self, params, epsf, relativeScale=False, stepSizeCutoff=None):

664 """ 665 Finite difference the residual dictionary to get a dictionary 666 for the Jacobian. It will be indexed the same as the residuals. 667 Note: epsf is either a scalar or an array. 668 If relativeScale is False then epsf is the stepsize used (it should 669 already be multiplied by typicalValues before Jacobian is called) 670 If relativeScale is True then epsf is multiplied by params. 671 The two previous statements hold for both scalar and vector valued 672 epsf. 673 """ 674 return self.jacobian_fd(params, epsf, 675 relativeScale, stepSizeCutoff)

676

677 - def GetJandJtJ(self,params):

678 j = self.GetJacobian(params) 679 mn = scipy.zeros((len(params),len(params)),scipy.float_) 680 681 for paramind in range(0,len(params)): 682 for paramind1 in range(0,len(params)): 683 sum = 0.0 684 for kys in j.keys(): 685 sum = sum + j.get(kys)[paramind]*j.get(kys)[paramind1] 686 687 mn[paramind][paramind1] = sum 688 return j,mn

689

690 - def GetJandJtJInLogParameters(self,params):

691 # Formula below is exact if you have perfect data. If you don't 692 # have perfect data (residuals != 0) you get an extra term when you 693 # compute d^2(cost)/(dlogp[i]dlogp[j]) which is 694 # sum_resname (residual[resname] * jac[resname][j] * delta_jk * p[k]) 695 # but can be ignored when residuals are zeros, and maybe should be 696 # ignored altogether because it can make the Hessian approximation 697 # non-positive definite 698 pnolog = scipy.exp(params) 699 jac, jtj = self.GetJandJtJ(pnolog) 700 for i in range(len(params)): 701 for j in range(len(params)): 702 jtj[i][j] = jtj[i][j]*pnolog[i]*pnolog[j] 703 704 res = self.resDict(pnolog) 705 for resname in self.residuals.keys(): 706 for j in range(len(params)): 707 # extra term --- not including it 708 # jtj[j][j] += res[resname]*jac[resname][j]*pnolog[j] 709 jac.get(resname)[j] = jac.get(resname)[j]*pnolog[j] 710 711 return jac,jtj

712

713 - def hessian_elem(self, func, f0, params, i, j, epsi, epsj, 714 relativeScale, stepSizeCutoff, verbose):

715 """ 716 Return the second partial derivative for func w.r.t. parameters i and j 717 718 f0: The value of the function at params 719 eps: Sets the stepsize to try 720 relativeScale: If True, step i is of size p[i] * eps, otherwise it is 721 eps 722 stepSizeCutoff: The minimum stepsize to take 723 """ 724 origPi, origPj = params[i], params[j] 725 726 if relativeScale: 727 # Steps sizes are given by eps*the value of the parameter, 728 # but the minimum step size is stepSizeCutoff 729 hi, hj = scipy.maximum((epsi*abs(origPi), epsj*abs(origPj)), 730 (stepSizeCutoff, stepSizeCutoff)) 731 else: 732 hi, hj = epsi, epsj 733 734 if i == j: 735 params[i] = origPi + hi 736 fp = func(params) 737 738 params[i] = origPi - hi 739 fm = func(params) 740 741 element = (fp - 2*f0 + fm)/hi**2 742 else: 743 ## f(xi + hi, xj + h) 744 params[i] = origPi + hi 745 params[j] = origPj + hj 746 fpp = func(params) 747 748 ## f(xi + hi, xj - hj) 749 params[i] = origPi + hi 750 params[j] = origPj - hj 751 fpm = func(params) 752 753 ## f(xi - hi, xj + hj) 754 params[i] = origPi - hi 755 params[j] = origPj + hj 756 fmp = func(params) 757 758 ## f(xi - hi, xj - hj) 759 params[i] = origPi - hi 760 params[j] = origPj - hj 761 fmm = func(params) 762 763 element = (fpp - fpm - fmp + fmm)/(4 * hi * hj) 764 765 params[i], params[j] = origPi, origPj 766 767 self._notify(event = 'hessian element', i = i, j = j, 768 element = element) 769 if verbose: 770 print 'hessian[%i, %i] = %g' % (i, j, element) 771 772 return element

773

774 - def hessian(self, params, epsf, relativeScale = True, 775 stepSizeCutoff = None, jacobian = None, 776 verbose = False):

777 """ 778 Returns the hessian of the model. 779 780 epsf: Sets the stepsize to try 781 relativeScale: If True, step i is of size p[i] * eps, otherwise it is 782 eps 783 stepSizeCutoff: The minimum stepsize to take 784 jacobian: If the jacobian is passed, it will be used to estimate 785 the step size to take. 786 vebose: If True, a message will be printed with each hessian element 787 calculated 788 """ 789 790 nOv = len(params) 791 if stepSizeCutoff is None: 792 stepSizeCutoff = scipy.sqrt(_double_epsilon_) 793 794 params = scipy.asarray(params) 795 if relativeScale: 796 eps = epsf * abs(params) 797 else: 798 eps = epsf * scipy.ones(len(params),scipy.float_) 799 800 # Make sure we don't take steps smaller than stepSizeCutoff 801 eps = scipy.maximum(eps, stepSizeCutoff) 802 803 if jacobian is not None: 804 # Turn off the relative scaling since that would overwrite all this 805 relativeScale = False 806 807 jacobian = scipy.asarray(jacobian) 808 if len(jacobian.shape) == 0: 809 resDict = self.resDict(params) 810 new_jacobian = scipy.zeros(len(params),scipy.float_) 811 for key, value in resDict.items(): 812 new_jacobian += 2.0*value*scipy.array(jacobian[0][key]) 813 jacobian = new_jacobian 814 elif len(jacobian.shape) == 2: # Need to sum up the total jacobian 815 residuals = scipy.asarray(self.res(params)) 816 # Changed by rng7. I'm not sure what is meant by "sum up the 817 # total jacobian". The following line failed due to shape 818 # mismatch. From the context below, it seems that the dot 819 # product is appropriate. 820 #jacobian = 2.0*residuals*jacobian 821 jacobian = 2.0 * scipy.dot(residuals, jacobian) 822 823 # If parameters are independent, then 824 # epsilon should be (sqrt(2)*J[i])^-1 825 factor = 1.0/scipy.sqrt(2) 826 for i in range(nOv): 827 if jacobian[i] == 0.0: 828 eps[i] = 0.5*abs(params[i]) 829 else: 830 # larger than stepSizeCutoff, but not more than 831 # half of the original parameter value 832 eps[i] = min(max(factor/abs(jacobian[i]), stepSizeCutoff), 833 0.5*abs(params[i])) 834 835 ## compute cost at f(x) 836 f0 = self.cost(params) 837 838 hess = scipy.zeros((nOv, nOv), scipy.float_) 839 840 ## compute all (numParams*(numParams + 1))/2 unique hessian elements 841 for i in range(nOv): 842 for j in range(i, nOv): 843 hess[i][j] = self.hessian_elem(self.cost, f0, 844 params, i, j, 845 eps[i], eps[j], 846 relativeScale, stepSizeCutoff, 847 verbose) 848 hess[j][i] = hess[i][j] 849 850 return hess

851

852 - def hessian_log_params(self, params, eps, 853 relativeScale=False, stepSizeCutoff=1e-6, 854 verbose=False):

855 """ 856 Returns the hessian of the model in log parameters. 857 858 eps: Sets the stepsize to try 859 relativeScale: If True, step i is of size p[i] * eps, otherwise it is 860 eps 861 stepSizeCutoff: The minimum stepsize to take 862 vebose: If True, a message will be printed with each hessian element 863 calculated 864 """ 865 nOv = len(params) 866 if scipy.isscalar(eps): 867 eps = scipy.ones(len(params), scipy.float_) * eps 868 869 ## compute cost at f(x) 870 f0 = self.cost_log_params(scipy.log(params)) 871 872 hess = scipy.zeros((nOv, nOv), scipy.float_) 873 874 ## compute all (numParams*(numParams + 1))/2 unique hessian elements 875 for i in range(nOv): 876 for j in range(i, nOv): 877 hess[i][j] = self.hessian_elem(self.cost_log_params, f0, 878 scipy.log(params), 879 i, j, eps[i], eps[j], 880 relativeScale, stepSizeCutoff, 881 verbose) 882 hess[j][i] = hess[i][j] 883 884 return hess

885

886 - def CalcHessianInLogParameters(self, params, eps, relativeScale = False, 887 stepSizeCutoff = 1e-6, verbose = False):

888 return self.hessian_log_params(params, eps, relativeScale, 889 stepSizeCutoff, verbose)

890

891 - def CalcHessian(self, params, epsf, relativeScale = True, 892 stepSizeCutoff = None, jacobian = None, verbose = False):

893 """ 894 Finite difference the residual dictionary to get a dictionary 895 for the Hessian. It will be indexed the same as the residuals. 896 Note: epsf is either a scalar or an array. 897 If relativeScale is False then epsf is the stepsize used (it should 898 already be multiplied by typicalValues before Jacobian is called) 899 If relativeScale is True then epsf is multiplied by params. 900 The two previous statements hold for both scalar and vector valued 901 epsf. 902 """ 903 return self.hessian(params, epsf, relativeScale, 904 stepSizeCutoff, jacobian, verbose)

905

906 - def CalcResidualResponseArray(self, j, h):

907 """ 908 Calculate the Residual Response array. This array represents the change 909 in a residual obtained by a finite change in a data value. 910 911 Inputs: 912 (self, j, h) 913 j -- jacobian matrix to use 914 h -- hessian matrix to use 915 916 Outputs: 917 response -- The response array 918 """ 919 j,h = scipy.asarray(j), scipy.asarray(h) 920 [m,n] = j.shape 921 response = scipy.zeros((m,m),scipy.float_) 922 ident = scipy.eye(m,typecode=scipy.float_) 923 hinv = scipy.linalg.pinv2(h,1e-40) 924 tmp = scipy.dot(hinv,scipy.transpose(j)) 925 tmp2 = scipy.dot(j,tmp) 926 response = ident - tmp2 927 928 return response

929

930 - def CalcParameterResponseToResidualArray(self,j,h):

931 """ 932 Calculate the parameter response to residual array. This array 933 represents the change in parameter resulting from a change in data 934 (residual). 935 936 Inputs: 937 (self, j, h) 938 j -- jacobian matrix to use 939 h -- hessian matrix to use 940 941 Outputs: 942 response -- The response array 943 """ 944 j,h = scipy.asarray(j), scipy.asarray(h) 945 [m,n] = j.shape 946 response = scipy.zeros((n,m),scipy.float_) 947 hinv = scipy.linalg.pinv2(h,1e-40) 948 response = -scipy.dot(hinv,scipy.transpose(j)) 949 950 return response

951 952 ############################################################################ 953 # Getting/Setting variables below 954

955 - def SetExperimentCollection(self, exptColl):

956 self.exptColl = exptColl 957 958 for exptKey, expt in exptColl.items(): 959 exptData = expt.GetData() 960 for calcKey, calcData in exptData.items(): 961 for depVarKey, depVarData in calcData.items(): 962 sortedData = depVarData.items() 963 sortedData.sort() 964 for indVar, (value, uncert) in sortedData: 965 resName = (exptKey, calcKey, depVarKey, indVar) 966 res = Residuals.ScaledErrorInFit(resName, depVarKey, 967 calcKey, indVar, value, 968 uncert, exptKey) 969 self.residuals.setByKey(resName, res) 970 971 # Add in the PeriodChecks 972 for period in expt.GetPeriodChecks(): 973 calcKey, depVarKey, indVarValue = period['calcKey'], \ 974 period['depVarKey'], period['startTime'] 975 resName = (exptKey, calcKey, depVarKey, indVarValue, 976 'PeriodCheck') 977 res = Residuals.PeriodCheckResidual(resName, calcKey, depVarKey, 978 indVarValue, 979 period['period'], 980 period['sigma']) 981 self.residuals.setByKey(resName, res) 982 983 # Add in the AmplitudeChecks 984 for amplitude in expt.GetAmplitudeChecks(): 985 calcKey, depVarKey = amplitude['calcKey'], \ 986 amplitude['depVarKey'] 987 indVarValue0, indVarValue1 = amplitude['startTime'],\ 988 amplitude['testTime'] 989 resName = (exptKey, calcKey, depVarKey, indVarValue0, 990 indVarValue1, 'AmplitudeCheck') 991 res = Residuals.AmplitudeCheckResidual(resName, calcKey, 992 depVarKey, indVarValue0, 993 indVarValue1, 994 amplitude['period'], 995 amplitude['sigma'], 996 exptKey) 997 self.residuals.setByKey(resName, res) 998 999 # Add in the integral data 1000 for ds in expt.GetIntegralDataSets(): 1001 for var in ds['vars']: 1002 resName = (exptKey, ds['calcKey'], var, 'integral data') 1003 res = Residuals.IntegralDataResidual(resName, var, 1004 exptKey, 1005 ds['calcKey'], 1006 ds['trajectory'], 1007 ds['uncert_traj'], 1008 ds['interval']) 1009 self.residuals.setByKey(resName, res) 1010 1011 for ds in expt.scaled_extrema_data: 1012 ds['exptKey'] = expt.name 1013 ds['key'] = '%s_%simum_%s_%s' % (ds['var'], ds['type'], 1014 str(ds['minTime']), 1015 str(ds['maxTime'])) 1016 res = Residuals.ScaledExtremum(**ds) 1017 self.AddResidual(res)

1018 1019

1020 - def get_expts(self):

1021 return self.exptColl

1022

1023 - def set_var_optimizable(self, var, is_optimizable):

1024 for calc in self.get_calcs().values(): 1025 try: 1026 calc.set_var_optimizable(var, is_optimizable) 1027 except KeyError: 1028 pass 1029 self.params = self.calcColl.GetParameters()

1030 1031 GetExperimentCollection = get_expts 1032

1033 - def SetCalculationCollection(self, calcColl):

1034 self.calcColl = calcColl 1035 self.params = calcColl.GetParameters()

1036

1037 - def get_calcs(self):

1038 return self.calcColl

1039 1040 GetCalculationCollection = get_calcs 1041

1042 - def GetScaleFactors(self):

1043 return self.internalVars['scaleFactors']

1044

1045 - def GetResiduals(self):

1046 return self.residuals

1047

1048 - def GetCalculatedValues(self):

1049 return self.calcVals

1050

1051 - def GetInternalVariables(self):

1052 return self.internalVars

1053

Source Code for Module SloppyCell.Model_mod