1 """
2 Methods for generating "perfect data" hessians.
3 """
4 __docformat__ = "restructuredtext en"
5
6 import copy
7 import logging
8 logger = logging.getLogger('RxnNets.PerfectData')
9
10 import sets
11
12 import scipy
13 import scipy.integrate
14
15 import Dynamics
16
17 from SloppyCell import HAVE_PYPAR, my_rank, my_host, num_procs
18 if HAVE_PYPAR:
19 import pypar
20
22 """
23 Return a trajectory with func applied to each variable stored in the
24 trajectory
25 """
26 if only_nonderivs:
27 keys = [key for key in self.key_column.keys()
28 if not isinstance(key, tuple)]
29 else:
30 keys = None
31
32 ret_traj = traj.copy_subset(keys)
33 for var, col in traj.key_column.items():
34 vals = func(traj, var)
35 ret_traj.values[:,col] = vals
36
37 return ret_traj
38
41 """
42 Update the typical var values for a group of networks.
43
44 Find the maximum of each variable over the integrations. In each network
45 the typical value is set to fraction of that maximum. If that maximum value
46 is less than cutoff, the typical value is set to 1.
47
48 networks List of networks to work with
49 int_times List of corresponding integration endpoints
50 (ie. [(0, 100), (0, 50)])
51 fraction Relative size of typical value, compared to max value over
52 the integrations.
53 rtol Relative tolerance for the integrations.
54 cutoff Values below this are considered to be zero
55 """
56 max_vals = {}
57 for net, times in zip(networks, int_times):
58 traj = Dynamics.integrate(net, times, rtol=rtol, fill_traj=True)
59 for var_id in net.variables.keys():
60 curr_max = max(traj.get_var_traj(var_id))
61 max_vals[var_id] = max(curr_max, max_vals.get(var_id, 0))
62
63 for var_id, val in max_vals.items():
64 for net in networks:
65 if net.variables.has_key(var_id):
66 if val > cutoff:
67 net.set_var_typical_val(var_id, val*fraction)
68 else:
69 net.set_var_typical_val(var_id, 1.0)
70
71 return max_vals
72
74 """
75 This is an uncertainty that is fraction of the variable's typical value.
76 """
77 def sigmaFunc(traj, data_id):
78 sigma = traj.get_var_typical_val(data_id) * fraction
79 if sigma < cutoff:
80 logger.warn('sigma < cutoff value (%g) for variable %s! '
81 'Taking sigma = 1.' % (cutoff, data_id))
82 sigma = 1
83 return sigma
84 return sigmaFunc
85
86 -def discrete_data(net, params, pts, interval, vars=None, random=False,
87 uncert_func=typ_val_uncert(0.1, 1e-14)):
88 """
89 Return a set of data points for the given network generated at the given
90 parameters.
91
92 net Network to generate data for
93 params Parameters for this evaluation of the network
94 pts Number of data points to output
95 interval Integration interval
96 vars Variables to output data for, defaults to all species in net
97 random If False data points are distributed evenly over interval
98 If True they are spread randomly and uniformly over each
99 variable
100 uncert_func Function that takes in a trajectory and a variable id and
101 returns what uncertainty should be assumed for that variable,
102 either as a scalar or a list the same length as the trajectory.
103 """
104
105 if vars == None:
106 vars = net.species.keys()
107
108
109 var_times = {}
110 for var in vars:
111 if random:
112 var_times[var] = scipy.rand(pts) * (interval[1]-interval[0]) + interval[0]
113 else:
114 var_times[var] = scipy.linspace(interval[0], interval[1], pts)
115
116
117 int_times = sets.Set(scipy.ravel(var_times.values()))
118 int_times.add(0)
119 int_times = list(int_times)
120 int_times.sort()
121
122
123 traj = Dynamics.integrate(net, int_times, params=params, fill_traj=False)
124
125
126 data = {}
127 for var, times in var_times.items():
128 var_data = {}
129 data[var] = var_data
130
131
132 var_uncerts = uncert_func(traj, var)
133 for time in times:
134 val = traj.get_var_val(var, time)
135 if scipy.isscalar(var_uncerts):
136 uncert = var_uncerts
137 else:
138 index = traj._get_time_index(time)
139 uncert = var_uncerts[index]
140 var_data[time] = (val, uncert)
141 return data
142
143 -def hessian_log_params(sens_traj, data_ids=None, opt_ids=None,
144 fixed_sf=False, return_dict=False,
145 uncert_func=typ_val_uncert(1.0, 1e-14)):
146 """
147 Calculate the "perfect data" hessian in log parameters given a sensitivity
148 trajectory.
149
150 sens_traj Sensitivity trajectory of Network being considered.
151 data_ids A sequence of variable id's to assume we have data for. If
152 data_ids is None, all dynamic and assigned variables will be
153 used.
154 opt_ids A sequence of parameter id's to calculate derivatives with
155 respect to. The hessian is (len(opt_ids) x len(opt_ids)).
156 If opt_ids is None, all optimizable variables are considered.
157 fixed_sf If True, calculate the hessian assuming fixed scale factors.
158 return_dict If True, returned values are (hess, hess_dict). hess_dict is a
159 dictionary keyed on the elements of data_ids; each corresponding
160 value is the hessian assuming data only on a single variable.
161 hess is the sum of all these hessians
162 uncert_func Function that takes in a trajectory and a variable id and
163 returns what uncertainty should be assumed for that variable,
164 either as a scalar or a list the same length as the trajectory.
165 """
166 if data_ids is None:
167 data_ids = sens_traj.dynamicVarKeys + sens_traj.assignedVarKeys
168 if opt_ids is None:
169 opt_ids = sens_traj.optimizableVarKeys
170
171 data_sigmas = {}
172 for data_id in data_ids:
173 ds = uncert_func(sens_traj, data_id)
174 if scipy.isscalar(ds):
175 ds = scipy.zeros(len(sens_traj), scipy.float_) + ds
176 data_sigmas[data_id] = ds
177
178 vars_assigned = [data_ids[node::num_procs] for node in range(num_procs)]
179 for worker in range(1, num_procs):
180 logger.debug('Sending to worker %i.' % worker)
181
182
183
184 vars_needed = sets.Set(sens_traj.optimizableVarKeys)
185 vars_needed.union_update(vars_assigned[worker])
186 for var in vars_assigned[worker]:
187 vars_needed.union_update([(var, ov) for ov in opt_ids])
188 worker_traj = sens_traj.copy_subset(vars_needed)
189
190 worker_ds = dict([(var, data_sigmas[var])
191 for var in vars_assigned[worker]])
192 command = 'PerfectData.compute_sf_LMHessian_conts(sens_traj, data_ids,'\
193 'data_sigmas, opt_ids, fixed_sf)'
194 args = {'sens_traj': worker_traj, 'data_ids': vars_assigned[worker],
195 'data_sigmas': worker_ds, 'opt_ids': opt_ids,
196 'fixed_sf': fixed_sf}
197 pypar.send((command, args), worker)
198
199 hess_dict = compute_sf_LMHessian_conts(sens_traj, vars_assigned[0],
200 data_sigmas, opt_ids, fixed_sf)
201
202 for worker in range(1, num_procs):
203 logger.debug('Receiving from worker %i.' % worker)
204 hess_dict.update(pypar.receive(worker))
205
206 hess = scipy.sum(hess_dict.values(), axis=0)
207 if return_dict:
208 return hess, hess_dict
209 else:
210 return hess
211
214 hess_dict = {}
215 for data_id in data_ids:
216 if fixed_sf:
217 sf_deriv = dict([(id, 0) for id in opt_ids])
218 else:
219 sf_deriv = get_sf_derivs(sens_traj, data_id, data_sigmas[data_id],
220 opt_ids)
221 hess_dict[data_id] = computeLMHessianContribution(sens_traj, data_id,
222 data_sigmas[data_id],
223 opt_ids, sf_deriv)
224
225 return hess_dict
226
228
229
230 times = traj.get_times()
231 eventIndices = scipy.compress(scipy.diff(times) == 0,
232 scipy.arange(len(times)))
233 intervals = zip([0] + list(eventIndices + 1),
234 list(eventIndices + 1) + [len(times)])
235
236 return intervals
237
239 scaleFactorDerivs = {}
240 intTheorySq = 0
241 for start, end in get_intervals(traj):
242 y = traj.get_var_traj(dataId)[start:end]
243 times = traj.get_times()[start:end]
244 sigma = data_sigma[start:end]
245
246 value = scipy.integrate.simps((y/sigma)**2, times,
247 even='last')
248 intTheorySq += value
249
250 for optId in optIds:
251 optValue = abs(traj.get_var_traj(optId)[start:end])
252 sens = traj.get_var_traj((dataId, optId))[start:end]*optValue
253
254 numerator = scipy.integrate.simps(sens*y/sigma**2, times,
255 even='last')
256
257 scaleFactorDerivs.setdefault(optId, 0)
258 scaleFactorDerivs[optId] += -numerator
259
260 for optId in optIds:
261 if intTheorySq != 0:
262 scaleFactorDerivs[optId] /= intTheorySq
263 else:
264 scaleFactorDerivs[optId] = 0
265
266 return scaleFactorDerivs
267
270 LMHessian = scipy.zeros((len(optIds), len(optIds)), scipy.float_)
271
272
273 for start, end in get_intervals(traj):
274 times = traj.timepoints[start:end]
275 y = traj.getVariableTrajectory(dataId)[start:end]
276 sigma = data_sigma[start:end]
277
278 for optIndex1, optId1 in enumerate(optIds):
279
280
281 optValue = abs(traj.get_var_traj(optId1)[start:end])
282 sens1 = traj.get_var_traj((dataId, optId1))[start:end]*optValue
283 dB1 = scaleFactorDerivs[optId1]
284 for jj, optId2 in enumerate(optIds[optIndex1:]):
285 optIndex2 = jj + optIndex1
286
287 optValue = abs(traj.get_var_traj(optId2)[start:end])
288 sens2 = traj.get_var_traj((dataId, optId2))[start:end]*optValue
289 dB2 = scaleFactorDerivs[optId2]
290
291 integrand = (sens1 + dB1 * y) * (sens2 + dB2 * y)/ sigma**2
292
293
294
295
296 value = scipy.integrate.simps(integrand, times,
297 even='last')
298
299 LMHessian[optIndex1][optIndex2] += value
300 if optIndex1 != optIndex2:
301 LMHessian[optIndex2][optIndex1] += value
302
303 LMHessian /= (traj.timepoints[-1] - traj.timepoints[0])
304
305 return LMHessian
306