1 from PyDSTool import *
2 from PyDSTool.PyCont.misc import getFlowMaps, getFlowJac, getLeftEvecs
3
4 __all__ = ['adjointPRC', 'rotate_phase']
5
6 -def adjointPRC(model, limit_cycle, vname, freepar, numIntervals=500,
7 numCollocation=7, method='standard', spike_est=None,
8 saveData=False, doPlot=False, verbosity=0, fname=None,
9 force=False):
10 """Create adjoint PRC from a limit cycle using AUTO, given an
11 ODE_Generator or Model object.
12
13 spike_est (default None) provides an estimate of the spike time in
14 the limit cycle, if applicable. If set, this can be useful to
15 ensure the PRC is phase-shifted to start at the associated zero
16 (in the variable given by vname).
17
18 method: 'standard' is currently the only working method for calculating
19 the eigenvetors corresponding to the unit eigenvalue. 'cyclic' selects
20 the cyclic method.
21
22 Optionally save the PRC data and make plots.
23 Returns a dictionary containing the PRC, cycles with independent
24 variable values scaled to [0,1], the Jacobian, the flow maps,
25 the eigenvector corresponding to the unit eigenvalue, and
26 other algorithmic information.
27 """
28
29 assert isinstance(freepar, str), "Free parameter must be a single string"
30
31 try:
32 MDCont = ContClass(model)
33 except:
34 print "Problem with model argument"
35 raise
36
37
38 if limit_cycle is None:
39
40 PCargs = args(name='EQ1', type='EP-C')
41 PCargs.freepars = [freepar]
42 PCargs.StepSize = 1e-2
43 PCargs.MaxNumPoints = 300
44 PCargs.LocBifPoints = 'all'
45 PCargs.NumIntervals = numIntervals
46 PCargs.NumCollocation = numCollocation
47 PCargs.FuncTol = 1e-10
48 PCargs.VarTol = 1e-10
49 PCargs.SaveEigen = True
50 PCargs.SaveFlow = True
51 PCargs.SaveJacobian = True
52 PCargs.verbosity = 2
53 MDCont.newCurve(PCargs)
54
55 print 'Computing curve...'
56 start = clock()
57 MDCont['EQ1'].forward()
58 print 'done in %.3f seconds!' % (clock()-start)
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74 PCargs.name = 'LC1'
75 PCargs.type = 'LC-C'
76 PCargs.verbosity = verbosity
77
78
79 PCargs.initpoint = 'EQ1:H1'
80 PCargs.MinStepSize = 8e-3
81 PCargs.MaxStepSize = 8e-3
82 PCargs.StepSize = 8e-3
83 PCargs.MaxNumPoints = 4
84 PCargs.LocBifPoints = []
85 PCargs.StopAtPoints = 'B'
86 PCargs.NumIntervals = numIntervals
87 PCargs.NumCollocation = numCollocation
88 PCargs.FuncTol = 1e-10
89 PCargs.VarTol = 1e-10
90 PCargs.NumSPOut = 3
91 PCargs.SolutionMeasures = 'all'
92 PCargs.SaveEigen = True
93 PCargs.SaveFlow = True
94 PCargs.SaveJacobian = True
95 MDCont.newCurve(PCargs)
96
97
98 else:
99 try:
100 cycle = limit_cycle[MDCont.model.query('variables')]
101 except:
102 raise ValueError("Model variables must match cycle coordinate names")
103 assert vname in cycle.coordnames, "vname argument must be a coordinate"
104
105 assert cycle[-1] == cycle[0], "First and last point of cycle must be equal"
106
107 assert isinstance(limit_cycle, Pointset), \
108 "cycle argument must be a Pointset"
109 assert remain(MDCont.model.query('variables'),limit_cycle.coordnames)==[], \
110 "Model variables must match cycle coordinate names"
111
112
113 PCargs = args(name = 'LC1')
114 PCargs.type = 'LC-C'
115 PCargs.verbosity = verbosity
116 PCargs.freepars = [freepar]
117 PCargs.initcycle = copy(cycle)
118 PCargs.MinStepSize = 1e-6
119 PCargs.MaxStepSize = 1e-6
120 PCargs.StepSize = 1e-6
121 PCargs.MaxNumPoints = 4
122 PCargs.LocBifPoints = []
123 PCargs.StopAtPoints = 'B'
124 PCargs.NumIntervals = numIntervals
125 PCargs.NumCollocation = numCollocation
126 PCargs.FuncTol = 1e-10
127 PCargs.VarTol = 1e-10
128 PCargs.NumSPOut = 3
129 PCargs.SolutionMeasures = 'all'
130 PCargs.SaveEigen = True
131 PCargs.SaveFlow = True
132 PCargs.SaveJacobian = True
133 MDCont.newCurve(PCargs)
134
135
136
137 MDCont['LC1'].forward()
138 pt = MDCont['LC1'].getSpecialPoint('RG1')
139 if verbosity > 1:
140 print "\nRegular point chosen:", pt, "\n"
141 J = getFlowJac(pt, verbose=verbosity>1)
142
143 n = J.shape[0]
144
145 if verbosity > 1:
146 print "Computing flow maps..."
147 maps, ntst = getFlowMaps(n, pt, 'RG', method=method)
148
149 cycle_true = pt.labels['RG']['cycle']
150 nint = MDCont['LC1'].NumIntervals
151 ncol = MDCont['LC1'].NumCollocation
152
153 assert isinstance(MDCont.model, NonHybridModel), \
154 "Only use a single vector field for the model"
155 pars = MDCont.model.query('pars')
156 idxs = ncol*arange(nint)
157 flow_vecs = array([MDCont.model.Rhs(0, x, pars) for x in \
158 cycle_true[idxs]])
159
160 if verbosity > 1:
161 print "Computing left eigenvetors corresponding to the unit ", \
162 "eigenvalue along flow"
163 evec1 = getLeftEvecs(n, ntst, maps, flow_vecs, method=method,
164 verbose=verbosity>1)
165 t = cycle_true['t'][idxs]
166 PRC = Pointset(indepvararray=t/t[-1],
167 coordarray=evec1.T/t[-1],
168 coordnames=cycle.coordnames)
169
170 if spike_est is not None:
171 assert spike_est >= 0 and spike_est <= t[-1]
172 spike_est_phase = spike_est*1./cycle_true['t'][-1]
173 res = PRC.find(spike_est_phase)
174 if isinstance(res, tuple):
175 spike_est_ix = res[0]
176 else:
177 spike_est_ix = res
178 prcv = PRC[vname]
179 phase0_ix = find_zero_phase(prcv, spike_est_ix)
180 PRC = rotate_phase(PRC, phase0_ix)
181 phase0_t = PRC['t'][phase0_ix]*t[-1]
182 if verbosity > 1:
183 print "spike_est_ix =", spike_est_ix
184 print "phase0_ix =", phase0_ix
185 print "phase0_t =", phase0_t
186
187
188
189 cycle_true_ph0_ix = phase0_ix*ncol
190 cycle_true = rotate_phase(cycle_true, cycle_true_ph0_ix)
191
192 cycle.indepvararray /= cycle['t'][-1]
193 T = cycle_true['t'][-1]
194 cycle_true.indepvararray /= T
195
196 if doPlot:
197 cyc_offset = (max(cycle_true[vname])+min(cycle_true[vname]))/2.
198 cyc_scale = abs(max(cycle_true[vname])-min(cycle_true[vname]))
199 PRC_scale = abs(max(PRC[vname])-min(PRC[vname]))
200 cyc_rescale = PRC_scale/cyc_scale
201 plt.figure()
202 plot(PRC['t'], PRC[vname], 'r', linewidth=2)
203 plot(cycle_true['t'], cyc_rescale*(cycle_true[vname]-cyc_offset))
204 plt.title('PRC overlay on (cycle%+.2e)*%.2e'%(cyc_offset,
205 cyc_rescale))
206 show()
207
208 save_objs = {'starting cycle': cycle, 'PRC cycle': cycle_true,
209 'PyCont args': PCargs, 'starting point': pt,
210 'PRC': PRC, 'nint': nint, 'ncol': ncol,
211 'Jacobian': J, 'maps': maps, 'evec1': evec1,
212 'period': T}
213 if saveData:
214
215 if fname is None:
216 savefile = MDCont.model.name + '_adjointPRC.pkl'
217 else:
218 savefile = str(fname)
219 try:
220
221 saveObjects(objlist=save_objs, filename=savefile, force=force)
222 except ValueError:
223 print "File already exists -- not saving return objects"
224
225 return_objs = save_objs
226 return_objs['PyCont'] = MDCont
227
228 return return_objs
229
230
232 """Phase shift a pointset (assumed to be a cycle) about index
233 phase0_ix, i.e. 'rotate' it to put the point at phase0_ix at the
234 beginning of the pointset.
235
236 NOTE: Does not update any label info that might be attached to pts!
237 """
238 assert phase0_ix > 0 and phase0_ix < len(pts), "phase 0 index out of range"
239 try:
240 t0 = pts['t'][phase0_ix]
241 parameterized = True
242 except PyDSTool_KeyError:
243 parameterized = False
244 pts_0 = pts[phase0_ix:]
245 if parameterized:
246 pts_0.indepvararray -= t0
247 pts_1 = pts[:phase0_ix]
248 if parameterized:
249 pts_1.indepvararray += (pts['t'][-1]-pts['t'][phase0_ix-1])
250 pts_0.extend(pts_1)
251 return pts_0
252
253
255 """prc should be of a single variable of interest"""
256 len_prc = len(prc)
257 assert est_ix < len(prc)
258 assert est_ix >= 0
259 zero_ixs = { 1: len_prc,
260 -1:-len_prc}
261 for dirn in [1, -1]:
262 ix = est_ix + dirn
263 while ix >= 0 and ix < len_prc:
264 if prc[ix]*prc[est_ix] < 0:
265 zero_ixs[dirn] = ix
266 break
267 else:
268 ix += dirn
269 dist_forward = zero_ixs[1] - est_ix
270 dist_backward = est_ix - zero_ixs[-1]
271 if dist_forward < len_prc:
272 phase0_ix = zero_ixs[1]
273 if dist_backward < dist_forward:
274 phase0_ix = zero_ixs[-1]
275 else:
276 phase0_ix = est_ix
277 return phase0_ix
278