-
Notifications
You must be signed in to change notification settings - Fork 65
/
Model.py
executable file
·3375 lines (3145 loc) · 152 KB
/
Model.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
"""General purpose (hybrid) model class, and associated hybrid trajectory
and variable classes.
Robert Clewley, March 2005.
A Model object's hybrid trajectory can be treated as a curve, or as
a mapping. Call the model object with the trajectory name, time(s), and
set the 'asmap' argument to be True to use an integer time to select the
trajectory segment. These are numbered from zero.
A trajectory value in a Model object's 'trajectories' dictionary
attribute is a HybridTrajectory object, having the following
attributes (among others):
timeInterval is the entire time interval for the trajectory.
timePartitions is a sequence of time_intervals (for each trajectory
segment in trajseq), and
trajSeq is a list of epoch or regime trajectory segments [traj_0, traj_1,
..., traj_(R-1)],
where traj_i is a callable Trajectory or HybridTrajectory object.
eventStruct is the event structure used to determine that trajectory.
events is a dictionary of event names -> list of times at which that
event took place.
modelNames is a list of the generators used for the trajectory (one per
partition).
variables is a dictionary that mimics the variables of the trajectory.
"""
# ----------------------------------------------------------------------------
## PyDSTool imports
from . import Generator, Events, ModelContext
from .utils import *
from .common import *
from .errors import *
from .Interval import *
from .Trajectory import *
from .Variable import *
from .Points import *
from .ModelSpec import *
from .Symbolic import isMultiRef
from .parseUtils import isHierarchicalName, NAMESEP, mapNames, symbolMapClass
## Other imports
import math, sys
from numpy import Inf, NaN, isfinite, sign, abs, array, arange, \
zeros, concatenate, transpose, shape
from numpy import sometrue, alltrue, any, all
import copy
from time import perf_counter
import pprint
__all__ = ['Model', 'HybridModel', 'NonHybridModel',
'boundary_containment', 'boundary_containment_by_postproc',
'boundary_containment_by_event',
'domain_test']
# ----------------------------------------------------------------------------
class boundary_containment(ModelContext.qt_feature_leaf):
# not implemented using metrics because the metrics are trivial
# and cause a lot of overhead for this often-evaluated feature
def __init__(self, name, description='', pars=None):
ModelContext.qt_feature_leaf.__init__(self, name, description, pars)
try:
pars.thresh
except AttributeError:
raise ValueError("Missing threshold specification")
try:
# tolerance for small rounding errors
pars.abseps
except AttributeError:
self.pars.abseps = 0
try:
assert pars.interior_dirn in [-1, 0, 1]
except AttributeError:
raise ValueError("Missing interior direction specification")
except AssertionError:
raise ValueError("Invalid interior direction specification value "
" use -1 for 'below', 1 for 'above', or 0 for 'discrete domain'")
try:
self.pars.coordname
except AttributeError:
# test all coords at once
self.pars.coordname = None
else:
assert isinstance(self.pars.coordname, str), \
"Coordinate name must be a string"
def evaluate(self, target):
raise NotImplementedError("Only call this method on a concrete "
"sub-class")
class boundary_containment_by_event(boundary_containment):
def __init__(self, name, description='', pars=None):
boundary_containment.__init__(self, name, description, pars)
try:
self.pars.bd_eventname
except AttributeError:
raise ValueError("Missing boundary event name")
# assume that the supplied model will correspond to the source of
# trajectories in evaluate method
try:
self.pars.model
except AttributeError:
raise ValueError("Missing model associated with event")
# have not verified that event is present in model
def evaluate(self, traj):
# verify whether event exists and was flagged in associated model
try:
evpts = traj.getEvents(self.pars.bd_eventname)
except ValueError as errinfo:
print(errinfo)
raise RuntimeError("Could not find flagged events for this trajectory")
try:
evpt = evpts[self.pars.coordname]
except KeyError:
raise ValueError("No such coordinate %s in the defined event"%self.pars.coordname)
except TypeError:
# no events of this kind were found, so passed feature eval test
# dereferencing None (unsubscriptable object)
if self.pars.abseps > 0:
# would like to re-evaluate event at its threshold+abseps, but
# leave for now
print("Warning -- Boundary containment feature %s:"%self.name)
print(" Check for uncertain case using events not implemented")
self.results.output = None
self.results.uncertain = False
else:
self.results.output = None
self.results.uncertain = False
return True
else:
# event found
if self.pars.abseps > 0:
# would like to re-evaluate event at its threshold+abseps, but
# leave for now
print("Warning -- Boundary containment feature %s:"%self.name)
print(" Check for uncertain case using events not implemented")
self.results.output = evpt[0] # only use first event (in case not Terminal)
self.results.uncertain = False
else:
self.results.output = evpt[0] # only use first event (in case not Terminal)
self.results.uncertain = False
return False
def _find_idx(self):
"""Helper function for finding index in trajectory meshpoints
at which containment first failed."""
if self.results.satisfied:
# Trajectory satisfied constraint!
return None
return len(self.results.output)
class boundary_containment_by_postproc(boundary_containment):
def evaluate(self, traj):
diffs = [p - self.pars.thresh for p in \
traj.sample(coords=self.pars.coordname)]
if self.pars.verbose_level > 1:
print("%s diffs in coord %s ="%(self.name,self.pars.coordname) + ", %s" % diffs)
res_strict = array([sign(d) \
== self.pars.interior_dirn for d in diffs])
satisfied_strict = alltrue(res_strict)
if self.pars.abseps > 0:
if self.pars.interior_dirn == 0:
# especially for discrete domains
res_loose = array([abs(d) < self.pars.abseps for d in diffs])
else:
res_loose = array([sign(d + self.pars.interior_dirn*self.pars.abseps) \
== self.pars.interior_dirn for d in diffs])
satisfied_loose = alltrue(res_loose)
self.results.output = res_loose
# if p values are *outside* thresh by up to abseps amount
# then flag this as 'uncertain' for use by domain_test class's
# transversality testing
self.results.uncertain = satisfied_loose and not satisfied_strict
return satisfied_loose
else:
self.results.output = res_strict
self.results.uncertain = sometrue(array(diffs)==0)
return alltrue(res_strict)
def _find_idx(self):
"""Helper function for finding index in trajectory meshpoints
at which containment first failed"""
if self.results.satisfied:
# Trajectory satisfied constraint!
return None
res = self.results.output
if res[0] == -1:
adjusted_res = list((res + 1) != 0)
elif res[0] == 1:
adjusted_res = list((res - 1) != 0)
else:
# starts with 0 already
adjusted_res = list(res != 0)
# find first index at which value is non-zero
# should never raise ValueError because this method is
# only run if there was a sign change found
return adjusted_res.index(True)
class domain_test(ModelContext.qt_feature_node):
def __init__(self, name, description='', pars=None):
ModelContext.qt_feature_node.__init__(self, name, description, pars)
try:
self.pars.interval
except AttributeError:
raise ValueError("Missing domain interval specification")
try:
# tolerance for small rounding errors
self.pars.abseps
except AttributeError:
if isinstance(self.pars.interval, Interval):
# Interval type passed, inherit its abseps
self.pars.abseps = self.pars.interval._abseps
else:
self.pars.abseps = 0
if not isinstance(self.pars.interval, (tuple, list)):
# assume a singleton numeric type passed
self.pars.interval = [self.pars.interval, self.pars.interval]
elif len(self.pars.interval)==1:
# singleton passed, so copy the value for both
# "endpoints" so that xlo_bc etc. below will work
self.pars.interval = [self.pars.interval[0],
self.pars.interval[0]]
self.isdiscrete = self.pars.interval[0] == self.pars.interval[1]
try:
self.pars.coordname
except AttributeError:
raise ValueError("Missing coordinate name")
try:
self.pars.derivname
except AttributeError:
raise ValueError("Missing coordinate derivative name")
# multiply interior directions by the integer value of not self.isdiscrete
# in order to set them to zero when the "interval" is actually a singleton
# value for a discrete domain. That fixes the boundary containment evaluation
# code which compares the sign of coord differences with that interior
# direction value.
xlo_bc = boundary_containment_by_postproc('x_test_lo',
description='Test x lower bound',
pars=args(thresh=self.pars.interval[0],
interior_dirn=1*int(not self.isdiscrete),
abseps=self.pars.abseps,
coordname=self.pars.coordname))
xhi_bc = boundary_containment_by_postproc('x_test_hi',
description='Test x upper bound',
pars=args(thresh=self.pars.interval[1],
interior_dirn=-1*int(not self.isdiscrete),
abseps=self.pars.abseps,
coordname=self.pars.coordname))
dxlo_bc = boundary_containment_by_postproc('dx_test_lo',
description='Test dx at lower x bound',
pars=args(thresh=0,
interior_dirn=1*int(not self.isdiscrete),
abseps=0,
coordname=self.pars.derivname))
dxhi_bc = boundary_containment_by_postproc('dx_test_hi',
description='Test dx at upper x bound',
pars=args(thresh=0,
interior_dirn=-1*int(not self.isdiscrete),
abseps=0,
coordname=self.pars.derivname))
self.subfeatures = {'x_test_lo': xlo_bc,
'dx_test_lo': dxlo_bc,
'x_test_hi': xhi_bc,
'dx_test_hi': dxhi_bc}
def evaluate(self, traj):
# Arg is a traj!
for sf in self.subfeatures.values():
self.propagate_verbosity(sf)
sf.super_pars = self.pars
sf.super_results = self.results
xlo_bc = self.subfeatures['x_test_lo']
xlo_test = xlo_bc(traj)
if xlo_bc.results.uncertain:
if self.pars.verbose_level > 0:
print("Lo bd uncertain")
if self.isdiscrete:
# accept uncertain case for discrete domain
xlo_test = True
else:
# check transversality at critical (boundary) value of domain
xlo_test = self.subfeatures['dx_test_lo'](traj)
xhi_bc = self.subfeatures['x_test_hi']
xhi_test = xhi_bc(traj)
if xhi_bc.results.uncertain:
if self.pars.verbose_level > 0:
print("Hi bd uncertain")
if self.isdiscrete:
# accept uncertain case for discrete domain
xhi_test = True
else:
# check transversality at critical (boundary) value of domain
xhi_test = self.subfeatures['dx_test_hi'](traj)
for sf in self.subfeatures.values():
self.results[sf.name] = sf.results
return xlo_test and xhi_test
def _find_idx(self):
"""Helper function for finding lowest index in trajectory meshpoints
at which domain test first failed"""
if self.results.satisfied:
# Trajectory satified domain conditions!
return None
lowest_idx = Inf
for sfname, sf in self.subfeatures.items():
if self.pars.verbose_level > 0:
print("\n %s %r" % (sfname, sf.results))
try:
res = list(self.results[sfname].output)
except AttributeError:
# dxdt transversality test was not run so ignore
continue
if sf.results.satisfied:
continue
if self.pars.verbose_level > 0:
print(res)
# Find first index at which value is non-zero.
# Will not raise ValueError because test satisfaction
# was already checked, so must have a zero crossing
idx = res.index(False)
if idx < lowest_idx:
lowest_idx = idx
return lowest_idx
# -------------------
class Model(object):
"""
General-purpose Hybrid and Non-Hybrid Model abstract class.
"""
_needKeys = ['name', 'modelInfo']
_optionalKeys = ['ics', 'mspecdict', 'verboselevel',
'norm', 'tdata', 'eventPars', 'abseps']
# query keys for 'query' method
_querykeys = ['pars', 'parameters', 'events', 'submodels',
'ics', 'initialconditions', 'vars', 'variables',
'auxvariables', 'auxvars', 'vardomains', 'pardomains',
'abseps']
# valid keys for 'set' method
_setkeys = ['pars', 'algparams', 'checklevel', 'abseps',
'ics', 'inputs', 'tdata', 'restrictDSlist',
'globalt0', 'verboselevel', 'inputs_t0']
def __init__(self, legit, *a, **kw):
# legit is a way to ensure that instances of this abstract class
# are not created directly
if not legit==True:
# use explicit comparison to True otherwise kw argument will
# eval to True, which is not what we want
raise RuntimeError("Only use HybridModel or NonHybridModel classes")
if len(a) > 0:
if len(a) == 1 and isinstance(a[0], dict):
if intersect(a[0].keys(),kw.keys()) != []:
raise ValueError("Cannot have initialization keys "
"common to both dictionary and keyword arguments")
kw.update(a[0])
else:
raise ValueError("Non-keyword arguments must be a single "
"dictionary")
try:
self.name = kw['name']
# modelInfo is a dict mapping model names --> a dict of:
# 'dsi': GeneratorInterface or ModelInterface object
# 'swRules': dict of switching rules (transitions between
# trajectory segments)
# 'globalConRules': list of global consistency DS names
# 'domainTests': dictionary of variable name -> domain
self.modelInfo = kw['modelInfo']
except KeyError:
raise KeyError('Necessary keys missing from argument')
_foundKeys = len(self._needKeys)
# by default, 'observables' are all variables that are common
# to the Model/Generator objects in self.modelInfo, and the
# 'internals' are those remaining generate obsvars, intvars...
# aux vars are those aux vars common to *all* sub-models.
# create self.obsvars, self.intvars, self.auxvars
self.defaultVars()
# trajectories is a dict of trajectory segments (in sequence)
# ... can only add one at a time!
self.trajectories = {}
self.trajectory_defining_args = {}
# registry of generators or models (depending on sub-class)
# Using registry provides a shortcut for accessing a sub-model regardless
# of whether it's a Generator or a Model class
self.registry = {}
for name, infodict in self.modelInfo.items():
# set super model tag of ds object (which is either a
# ModelInterface or Generator)
try:
infodict['dsi']._supermodel = self.name
except KeyError:
raise TypeError("Invalid modelInfo entry found with name %s"%name)
self.registry[name] = infodict['dsi'].model
self.diagnostics = Diagnostics()
# set initial conditions if specified already (if not, they must
# be specified before or during when compute() is called)
self.icdict = {}
if 'ics' in kw:
self.icdict = dict(kw['ics'])
_foundKeys += 1
else:
self.icdict = {}.fromkeys(self.allvars, NaN)
if 'tdata' in kw:
self.tdata = kw['tdata']
_foundKeys += 1
else:
self.tdata = None
if 'norm' in kw:
self._normord = kw['norm']
_foundKeys += 1
else:
self._normord = 2
if 'abseps' in kw:
# impose this absolute epsilon (small scale) on all components
self._abseps = kw['abseps']
_foundKeys += 1
else:
# All components use their defaults
self._abseps = self.query('abseps')
if 'verboselevel' in kw:
if kw['verboselevel'] in [0,1,2]:
self.verboselevel = kw['verboselevel']
else:
raise ValueError("Verbosity level value must be 0, 1, or 2")
_foundKeys += 1
else:
self.verboselevel = 0
if 'mspecdict' in kw:
self._mspecdict = kw['mspecdict']
_foundKeys += 1
else:
self._mspecdict = None
if 'eventPars' in kw:
self._eventPars = kw['eventPars']
_foundKeys += 1
else:
self._eventPars = {}
if _foundKeys < len(kw):
raise KeyError('Invalid keys found in arguments')
# If not already created, a True result for
# self.haveJacobian() means that the Model will have a
# Jacobian method made during compute(), which
# references the appropriate _auxfn_Jac function in the Generator
# objects.
#
# Use this dict to record if any external input t0 time shift
# values are linked to a parameter value, for automatic
# updating before trajectory computation. (Internal use)
self._inputt0_par_links = {}
def __len__(self):
"""Return number of sub-models"""
return len(self.registry)
def sub_models(self):
"""Return a list of all sub-model instances (model interfaces or generators)"""
return list(self.registry.values())
def _makeDefaultVarNames(self):
"""Return default observable, internal, and auxiliary variable names
from modelInfo."""
obsvars = []
auxvars = []
all_known_varnames = []
for infodict in self.modelInfo.values():
varnames = infodict['dsi'].query('variables')
all_known_varnames.extend(varnames)
auxvarnames = infodict['dsi'].query('auxvariables')
if auxvars == []:
# first ds to have auxvars, so just add them all
auxvars.extend(auxvarnames)
else:
auxvars = intersect(auxvars, auxvarnames)
if obsvars == []:
# first ds, so add them all
obsvars = varnames
else:
obsvars = intersect(obsvars, varnames)
intvars = remain(all_known_varnames, obsvars)
return (obsvars, intvars, auxvars)
def _generateParamInfo(self):
"""Record parameter info locally, for future queries.
Internal use only.
"""
# use query method in case model in registry is a wrapped Generator
# that uses _ versions of hierarchical names that are used natively
# here.
self.pars = {}
for model in self.registry.values():
try:
self.pars.update(model.query('pars'))
except AttributeError:
# no pars present
pass
def showDef(self, target=None, type=''):
"""type = 'spec', 'auxspec', 'auxfnspec', 'events', or 'modelspec'
(leave blank for the first *four* together).
'spec', 'auxspec' and 'auxfnspec', 'events' refer to the compiled
target language code for the specifications. 'modelspec' refers to
the pre-compiled abstract specifications of the model."""
if target is None:
print("Use showInfo() to find names of defined sub-models")
return
else:
showAll = type==''
try:
if type=='spec' or showAll:
self.registry[target].showSpec()
if type=='auxspec' or showAll:
self.registry[target].showAuxSpec()
if type=='auxfnspec' or showAll:
self.registry[target].showAuxFnSpec()
if type=='events' or showAll:
self.registry[target].showEventSpec()
if type=='modelspec':
if self._mspecdict is None:
raise PyDSTool_ExistError("Cannot use this function "
"for models not defined through ModelSpec")
info(self._mspecdict[target]['modelspec'].flattenSpec(\
[self.modelInfo[target]['dsi'].get('indepvariable').name]))
except KeyError:
raise ValueError("Model named %s is not known"%target)
def showSpec(self):
for ds in self.registry.values():
ds.showSpec()
def showAuxSpec(self):
for ds in self.registry.values():
ds.showAuxSpec()
def showAuxFnSpec(self):
for ds in self.registry.values():
ds.showAuxFnSpec()
def showEventSpec(self):
for ds in self.registry.values():
ds.showEventSpec()
def current_defining_args(self):
return args(pars=self.pars, ics=self.icdict,
tdata=self.tdata)
def has_exact_traj(self, trajname, info):
"""Compare self.pars, self.icdict and self.tdata
against what's stored for a previously computed trajectory,
so that re-computation can be avoided.
"""
try:
return info == self.trajectory_defining_args[trajname]
except KeyError:
# not even a trajectory of this name
return False
def _prepareCompute(self, trajname, **kw):
foundKeys = 0
if 'verboselevel' in kw:
self.set(verboselevel=kw['verboselevel'])
foundKeys += 1
else:
self.set(verboselevel=0)
if 'ics' in kw:
self.icdict = dict(kw['ics'])
foundKeys += 1
if 'pars' in kw:
self.set(pars=kw['pars'])
foundKeys += 1
if 'tdata' in kw:
tdata = kw['tdata']
foundKeys += 1
# print "tdata in kw of %s (%s):"%(self.name, type(self)), tdata
else:
tdata = self.tdata
# print "tdata from self of %s (%s):"%(self.name, type(self)), tdata
if 'force' in kw:
force_overwrite = kw['force']
foundKeys += 1
else:
force_overwrite = False
if len(kw) != foundKeys:
raise PyDSTool_KeyError('Invalid argument keys passed to compute()')
if tdata is None:
raise ValueError("tdata must be specified")
if len(tdata) == 1:
assert isinstance(tdata, float) or isinstance(tdata, int), \
'tdata must be either a single number or a pair'
t0_global = tdata[0]
t1_global = Inf
elif len(tdata) == 2:
t0_global = tdata[0]
t1_global = tdata[1]
else:
raise ValueError('tdata argument key may be either a single '
'float or a pair of floats')
if not force_overwrite:
assert trajname not in self.trajectories, \
'Trajectory name already exists'
assert self.modelInfo != {}, \
'No Generator or Model objects defined for this model'
return tdata, t0_global, t1_global, force_overwrite
def query(self, querykey=''):
"""Return info about Model set-up.
Valid query key: 'pars', 'parameters', 'events', 'submodels',
'ics', 'initialconditions', 'vars', 'variables',
'auxvars', 'auxvariables', 'vardomains', 'pardomains', 'abseps'
"""
assert isinstance(querykey, str), \
("Query argument must be a single string")
if querykey not in self._querykeys:
print('Valid query keys are: %r' % self._querykeys)
print("('events' key only queries model-level events, not those")
print(" inside sub-models)")
if querykey != '':
raise TypeError('Query key '+querykey+' is not valid')
if querykey in ['pars', 'parameters']:
result = copy.copy(self.pars)
elif querykey in ['ics', 'initialconditions']:
result = copy.copy(self.icdict)
elif querykey == 'events':
result = {}
for dsName, model in self.registry.items():
try:
result.update(model.eventstruct.events)
except AttributeError:
# ds is a ModelInterface, not a Generator
result.update(model.query('events'))
elif querykey == 'submodels':
result = self.registry
elif querykey in ['vars', 'variables']:
result = copy.copy(self.allvars)
elif querykey in ['vardomains', 'xdomains']:
result = {}
# accumulate domains from each sub-model for regular variables
for model in self.registry.values():
vardoms = model.query('vardomains')
if len(result)==0:
result.update(vardoms)
else:
for vname, vdom in result.items():
if vdom.issingleton:
# singleton
vdom_lo = vdom.get()
vdom_hi = vdom_lo
else:
# range
vdom_lo = vdom[0]
vdom_hi = vdom[1]
if vardoms[vname][0] < vdom[0]:
vdom[0] = vardoms[vname][0]
if vardoms[vname][1] > vdom[1]:
vdom[1] = vardoms[vname][1]
if vdom._abseps < result[vname]._abseps:
# have to keep abseps the tightest
# of any of the instances for safety
result[vname]._abseps = vdom._abseps
result[vname] = vdom
# remaining vars are promoted aux vars
for vname in remain(self.allvars, result.keys()):
result[vname] = Interval(vname, float, [-Inf, Inf])
elif querykey in ['pardomains', 'pdomains']:
result = {}
# accumulate domains from each sub-model for regular variables
for model in self.registry.values():
pardoms = model.query('pardomains')
if len(result)==0:
result.update(pardoms)
else:
for pname, pdom in result.items():
if pdom.issingleton:
# singleton
pdom_lo = pdom.get()
pdom_hi = pdom_lo
else:
# range
pdom_lo = pdom[0]
pdom_hi = pdom[1]
if pardoms[pname][0] < pdom[0]:
pdom[0] = pardoms[pname][0]
if pardoms[pname][1] > pdom[1]:
pdom[1] = pardoms[pname][1]
if pdom._abseps < result[pname]._abseps:
# have to keep abseps the tightest
# of any of the instances for safety
result[pname]._abseps = pdom._abseps
result[pname] = pdom
elif querykey in ['auxvars', 'auxvariables']:
result = copy.copy(self.auxvars)
elif querykey == 'abseps':
result = min([ds.query('abseps') for ds in self.registry.values()])
return result
def getEventMappings(self, dsName):
try:
return self.modelInfo[dsName]['swRules']
except KeyError:
raise NameError("Sub-model %s not found in model"%dsName)
def setPars(self, p, val):
# process multirefs first, then hierarchical names
# calls itself recursively to resolve all names
if isMultiRef(p):
# e.g to change a group of numerically indexed parameter names
# such as p0 - p9 or comp1.p0.g - comp1.p9.g
# extract numeric range of pars
# [ and ] are guaranteed to be present, from isMultiRef()
lbrace = p.find('[')
rbrace = p.find(']')
if rbrace < lbrace:
raise ValueError("Invalid multiple reference to pars")
rootname = p[:lbrace]
rangespec = p[lbrace+1:rbrace].split(',')
try:
remainder = p[rbrace+1:]
except KeyError:
# no more of p after this multireference
remainder = ''
if len(rangespec) != 2:
raise ValueError("Invalid multiple reference to pars")
loix = int(rangespec[0])
hiix = int(rangespec[1])
if loix >= hiix:
raise ValueError("Invalid multiple reference to pars")
# call setPars for each resolved name (these may include further
# multi references or hierarchical names
for ix in range(loix, hiix+1):
self.setPars(rootname+str(ix), val)
elif isHierarchicalName(p):
if self._mspecdict is None:
raise PyDSTool_ExistError("Cannot use this functionality for"
" models not defined through ModelSpec")
# find out: is root of p a valid 'type' in model spec?
# find all occurrences of last p
allFoundNames = []
for mspecinfo in self._mspecdict.values():
foundNames = searchModelSpec(mspecinfo['modelspec'], p)
# don't add duplicates
allFoundNames.extend(remain(foundNames,allFoundNames))
# if allFoundNames == [] then either the hierarchical name was
# a specific reference, or the type matching is invalid.
# either way, we can just call set(p) and let that resolve the issue
# (note that all multi-refs will have been dealt with by this point)
if allFoundNames == []:
self.set(pars={p: val})
else:
self.set(pars={}.fromkeys(allFoundNames, val))
else:
self.set(pars={p: val})
def setICs(self, p, val):
# process multirefs first, then hierarchical names
# calls itself recursively to resolve all names
if isMultiRef(p):
# e.g to change a group of numerically indexed parameter names
# such as p0 - p9 or comp1.p0.g - comp1.p9.g
# extract numeric range of pars
# [ and ] are guaranteed to be present, from isMultiRef()
lbrace = p.find('[')
rbrace = p.find(']')
if rbrace < lbrace:
raise ValueError("Invalid multiple reference to initial conditions")
rootname = p[:lbrace]
rangespec = p[lbrace+1:rbrace].split(',')
try:
remainder = p[rbrace+1:]
except KeyError:
# no more of p after this multireference
remainder = ''
if len(rangespec) != 2:
raise ValueError("Invalid multiple reference to initial conditions")
loix = int(rangespec[0])
hiix = int(rangespec[1])
if loix >= hiix:
raise ValueError("Invalid multiple reference to initial conditions")
# call setICs for each resolved name (these may include further
# multi references or hierarchical names
for ix in range(loix, hiix+1):
self.setICs(rootname+str(ix), val)
elif isHierarchicalName(p):
if self._mspecdict is None:
raise PyDSTool_ExistError("Cannot use this functionality for"
" models not defined through ModelSpec")
# find out: is root of p a valid 'type' in model spec?
# find all occurrences of last p
allFoundNames = []
for mspecinfo in self._mspecdict.values():
foundNames = searchModelSpec(mspecinfo['modelspec'], p)
# don't add duplicates
allFoundNames.extend(remain(foundNames,allFoundNames))
# if allFoundNames == [] then either the hierarchical name was
# a specific reference, or the type matching is invalid.
# either way, we can just call set(p) and let that resolve the issue
# (note that all multi-refs will have been dealt with by this point)
if allFoundNames == []:
self.set(ics={p: val})
else:
self.set(ics={}.fromkeys(allFoundNames, val))
else:
self.set(ics={p: val})
def set(self, **kw):
"""Set specific parameters of Model. These will get passed on to
all Generators/sub-models that support these keys unless the
restrictDSList argument is set (only applies to keys algparams,
checklevel, and abseps).
Permitted keys: 'pars', 'algparams', 'checklevel', 'abseps',
'ics', 'inputs', 'tdata', 'restrictDSlist',
'globalt0', 'verboselevel', 'inputs_t0'
"""
for key in kw:
if key not in self._setkeys:
raise KeyError("Not a permitted parameter argument: %s"%key + \
". Allowed keys: "+str(self._setkeys))
if 'restrictDSlist' in kw:
restrictDSlist = kw['restrictDSlist']
else:
restrictDSlist = []
# Handle initial conditions here, because compute will pass
# the values on to the appropriate sub-models when they are called.
if 'ics' in kw:
self.icdict.update(filteredDict(dict(kw['ics']), self.icdict.keys()))
if 'tdata' in kw:
self.tdata = kw['tdata']
if 'abseps' in kw:
self._abseps = kw['abseps']
if 'inputs_t0' in kw:
for model in self.registry.values():
# Propagate values to sub-models.
# If any values are strings they must refer to a parameter, so
# here we evaluate them. We keep a record of these for later
# automated updating.
t0val_dict = kw['inputs_t0']
for inp, val in t0val_dict.items():
if isinstance(val, str):
try:
new_val = self.pars[val]
except KeyError:
raise ValueError("Input t0 to parameter link "
"invalid: no such parameter "+val)
t0val_dict[inp] = new_val
self._inputt0_par_links.update({inp: val})
elif isinstance(val, _num_types):
if inp in self._inputt0_par_links:
# no longer parameter-linked
del self._inputt0_par_links[inp]
else:
raise TypeError("Invalid type for input t0 value")
try:
model.set(inputs_t0=t0val_dict)
except AssertionError:
# generator doesn't involve inputs
pass
if 'verboselevel' in kw:
if kw['verboselevel'] in [0,1,2]:
self.verboselevel = kw['verboselevel']
else:
raise ValueError("Verbosity level value must be 0, 1, or 2")
for model in self.registry.values():
# propagate to sub-models
try:
model.set(verboselevel=self.verboselevel)
except KeyError:
# generator doesn't support verboselevel
pass
if restrictDSlist == []:
restrictDSlist = list(self.registry.keys())
# For the remaining keys, must propagate parameter changes to all
# sub-models throughout modelInfo structure.
#
# Changed by WES 10FEB06 to handle problem of 'new' pars being
# added if the parameter names do not exist in any generators
dsis = list(self.modelInfo.values())
numDSs = len(dsis)
# loop over keywords
for key, value in filteredDict(kw, ['ics', 'tdata',
'inputs_t0', 'restrictDSlist', 'globalt0'],
neg=True).items():
# keep track of the number of errors on this keyword
if isinstance(value, dict):
# keep track of entry errors for this key
entry_err_attr = {}
entry_err_val = {}
for entrykey, entryval in value.items():
entry_err_attr[entrykey] = 0
entry_err_val[entrykey] = 0
else:
entry_err_attr = 0
entry_err_val = 0
# try setting pars in each sub-model
for infodict in dsis:
callparsDict = {}
# select out the ones relevant to this sub-model
ds = infodict['dsi'].model
if hasattr(ds, '_validKeys'):
options = ds._validKeys
elif hasattr(ds, '_setkeys'):
options = ds._setkeys
if key in options:
if key in ['algparams', 'checklevel', 'abseps'] and \
ds.name not in restrictDSlist:
# only apply these keys to the restricted list
continue
if isinstance(value, dict):
for entrykey, entryval in value.items():
try:
ds.set(**{key:{entrykey:entryval}})
except PyDSTool_AttributeError:
entry_err_attr[entrykey] += 1
except PyDSTool_ValueError:
entry_err_val[entrykey] += 1
except AssertionError:
# key not valid for this type of ds
pass
else:
try:
ds.set(**{key: value})
except PyDSTool_AttributeError:
entry_err_attr += 1
except PyDSTool_ValueError:
entry_err_val += 1
except AssertionError:
# key not valid for this type of ds
pass
# Check that none of the entries in the dictionary caused errors
# in each sub-model
if isinstance(value, dict):
for entrykey, entryval in value.items():
if entry_err_attr[entrykey] == numDSs:
raise PyDSTool_AttributeError('Parameter does not' +\
' exist in any sub-model: %s = %f'%(entrykey,
entryval))
if entry_err_val[entrykey] == numDSs:
raise PyDSTool_ValueError('Parameter value error in' +\
' every sub-model: %s = %f'%(entrykey, entryval))
else:
# can't think of other ways for this error to crop up
pass
else:
if entry_err_attr == numDSs:
raise PyDSTool_AttributeError('Parameter does not exist' +\
' in any sub-model: %s'%key)
if entry_err_val == numDSs:
raise PyDSTool_ValueError('Parameter value error in' +\
' every sub-model: %s'%key)
del(entry_err_attr)
del(entry_err_val)
self._generateParamInfo()
def __getitem__(self, trajname):
try:
return self.trajectories[trajname]
except KeyError:
raise ValueError('No such trajectory.')
def __delitem__(self, trajname):
self._delTraj(trajname)
def _delTraj(self, trajname):
"""Delete a named trajectory from the database."""
try:
traj = self.trajectories[trajname]
except KeyError:
# a trajectory piece may have been created without
# the top-level trajectory ever being completed
# (e.g. after an unxpected error or ^C interruption)
##raise ValueError('No such trajectory.')
l = len(trajname)
for m in self.registry.values():
# delete all matching pieces (of form trajname + '_' + <digits>)
for n in m.trajectories.keys():