1 """Modelling project and associated classes.
2
3 For aiding model estimation using (hybrid) dynamical systems.
4
5 Robert Clewley, September 2007.
6 """
7
8 import copy
9 import sys, traceback
10 import numpy as npy
11
12
13 import Model, Generator, Events, ModelSpec, ModelConstructor, Symbolic, \
14 Trajectory
15 import utils, common, parseUtils
16 from errors import *
17
18
19
20
21
22 _classes = ['GenTransform', 'ModelTransform', 'ModelManager', 'feature',
23 'ql_feature_node', 'ql_feature_leaf', 'qt_feature_node',
24 'qt_feature_leaf', 'binary_feature', 'always_feature',
25 'condition', 'context', 'ModelLibrary', 'GeneratorInterface',
26 'ModelInterface', 'extModelInterface', 'intModelInterface']
27
28 _functions = ['extract_from_model']
29
30 __all__ = _classes + _functions
31
32
33
34
36 try:
37 r = copy.copy(model._mspecdict)
38 except AttributeError:
39 raise ValueError("Incompatible Model argument -- it contains "
40 "no ModelSpec info")
41 for key in ['algparams', 'mspecdict', 'targetlang']:
42 assert key in r
43
44 return r
45
46
47
48
49
51 """Store a set of related candidate model types, and within each, represent
52 various relevant "dimensions" along which the model can be augmented
53 structurally."""
54 - def __init__(self, name, spectype, indepdomain, depdomain,
55 pars=None, description=''):
56 self.name = name
57 self.spectype = spectype
58
59 self.instances = {}
60 self.indepdomain = indepdomain
61 self.depdomain = depdomain
62 self.pars = pars
63 self.description = ''
64
66 return self.instances[name]
67
69 if not isinstance(specs, common._seq_types):
70 specs = [specs]
71 for spec in specs:
72 if isinstance(spec, self.spectype):
73 self.instances[spec.name] = spec
74 spec.library_tag = self.name
75 else:
76 raise PyDSTool_TypeError("Spec of wrong type")
77
79 return "Model Library %s: %s"%(self.name, self.description)
80
81
305
306
307
387
388
389
390
392 """Model management and repository class."""
393
395 assert isinstance(name, str)
396 self.proj_name = name
397
398 self._mReg = MReg()
399
400 self.trans = None
401
402 self.instances = {}
403
405 if name in self._mReg:
406 return self._mReg.descs[name]
407 else:
408 raise KeyError('Model %s does not exist in registry'%name)
409
411 if name in self._mReg:
412 return self._mReg[name]
413 else:
414 raise KeyError('Model %s does not exist in registry'%name)
415
416 - def add(self, model_desc):
417 if not isinstance(model_desc, ModelConstructor.MDescriptor):
418 raise TypeError("Invalid model descriptor")
419 if not model_desc.validate():
420 raise ValueError("Model definition not successfully validated")
421 if model_desc not in self._mReg:
422 self._mReg.add(model_desc)
423 else:
424 raise KeyError('Model with this name already exists in registry')
425
427 if name in self._mReg:
428 del(self._mReg[name])
429 else:
430 raise KeyError('Model with this name does not exist in registry')
431
432 __delitem__ = remove
433
435 """Open a model transformation transaction"""
436 if self.trans is None:
437 self.trans = ModelTransform(name, self.__getitem__(name))
438 return self._mReg.descs[name]
439 else:
440 raise AssertionError("A transaction is already open")
441
443 if self.trans is None:
444 raise AssertionError("No transaction open")
445 else:
446 self.trans = None
447
449 if self.trans is None:
450 raise AssertionError("No transaction open")
451 else:
452 self.add(self.trans.commit(new_name))
453 self.trans = None
454
455 - def build(self, name, icvalues=None, parvalues=None,
456 inputs=None, tdata=None):
457 try:
458 mdesc = copy.deepcopy(self._mReg[name])
459 except KeyError:
460 raise KeyError("No such model description")
461 for gd in mdesc.generatorspecs.values():
462 gd.modelspec.flattenSpec(ignoreInputs=True, force=True)
463 filt_keys = ('userevents', 'userfns', 'unravelInfo',
464 'inputs', 'checklevel', 'activateAllBounds',
465 'generatorspecs', 'indepvar',
466 'parvalues', 'icvalues', 'reuseTerms',
467 'withJac', 'withJacP', 'tdata',
468 'abseps', 'eventtol', 'eventPars',
469 'withStdEvts', 'stdEvtArgs')
470 if icvalues is not None:
471 mdesc.icvalues.update(icvalues)
472 if parvalues is not None:
473 mdesc.parvalues.update(parvalues)
474 if inputs is not None:
475 mdesc.inputs.update(inputs)
476 if tdata is not None:
477 mdesc.tdata = tdata
478 if not mdesc.isinstantiable(True):
479 raise ValueError("Model description incomplete: not instantiable")
480
481
482
483
484 mc = ModelConstructor.ModelConstructor(mdesc.name,
485 **common.filteredDict(dict(mdesc), filt_keys))
486 assert len(mdesc.generatorspecs) > 0, "No Generator descriptions found"
487 for gdname, gd in mdesc.generatorspecs.iteritems():
488 if gd.userEvents is not None:
489 mc.addEvents(gdname, gd.userEvents)
490 if gd.userFunctions is not None:
491 mc.addFunctions(gdname, gd.userFunctions)
492 if gd.userEventMaps is not None:
493 for em in gd.userEventMaps:
494 try:
495
496 evname, target, evmap = em
497 except ValueError:
498
499 evname, target = em
500 evmap = None
501 mc.mapEvent(gdname, evname, target, evmap)
502 model = mc.getModel()
503 self._mReg[name] = model
504
505 self.instances = self._mReg.instances
506
508 if verbose == 0:
509 outputStr = 'MProject: '+self.proj_name
510 elif verbose > 0:
511 outputStr = 'MProject: '+self.proj_name
512 if len(self._mReg):
513 for m in self._mReg:
514 outputStr += "\n" + m._infostr(verbose-1)
515 else:
516 outputStr += 'No models in MProject '+self.proj_name
517 return outputStr
518
521
522 __str__ = __repr__
523
524 - def info(self, verboselevel=1):
526
527
528
529
530
531
532
534 """End users of concrete sub-classes provide (required) evaluate method and
535 (optional) prepare, finish methods."""
536 - def __init__(self, name, description='', pars=None, ref_traj=None):
537 self.name = name
538 self.description = description
539 if pars is None:
540 self.pars = common.args()
541 elif isinstance(pars, dict):
542 self.pars = common.args(**pars)
543 elif isinstance(pars, common.args):
544 self.pars = pars
545 else:
546 raise PyDSTool_TypeError("Invalid type for pars argument")
547 if 'verbose_level' not in self.pars:
548 self.pars.verbose_level = 0
549 if 'debug' not in self.pars:
550 self.pars.debug = False
551
552 if 'penalty' not in self.pars:
553 self.pars.penalty = 2000.
554 self.ref_traj = ref_traj
555 self.results = common.args()
556 self.super_pars = common.args()
557 self.super_results = common.args()
558
559
560 try:
561 self._local_init()
562 except AttributeError:
563 pass
564 self.subfeatures = []
565
568
570 try:
571 res = self.name == other.name
572 except AttributeError:
573 return False
574 if hasattr(self, 'subfeatures'):
575 if hasattr(other, 'subfeatures'):
576 res = res and self.subfeatures == other.subfeatures
577 else:
578 return False
579 elif hasattr(other, 'subfeatures'):
580 return False
581 return res
582
584 return not self == other
585
587 """Internal function for finding index in trajectory meshpoints
588 at which containment first failed. Defaults to returning None and
589 must be overridden by a class that has access to a mesh."""
590 return None
591
593 try:
594 self.prepare(target)
595 satisfied = self.evaluate(target)
596 except KeyboardInterrupt:
597 raise
598 except:
599 display_error = self.pars.verbose_level > 0 or self.pars.debug
600 if display_error:
601 exceptionType, exceptionValue, exceptionTraceback = sys.exc_info()
602 print "******************************************"
603 print "Problem evaluating feature:", self.name
604 print " ", exceptionType, exceptionValue
605 for line in traceback.format_exc().splitlines()[-12:-1]:
606 print " " + line
607 print " originally on line:", traceback.tb_lineno(exceptionTraceback)
608 if self.pars.debug:
609 raise
610 else:
611 print "(Proceeding as 'unsatisfied')\n"
612 satisfied = False
613 if hasattr(self, 'metric'):
614 self.metric.results = self.pars.penalty * \
615 npy.ones((self.metric_len,), float)
616 for sf in self.subfeatures:
617 if hasattr(sf, 'metric'):
618 sf.metric.results = self.pars.penalty * \
619 npy.ones((sf.metric_len,), float)
620 if satisfied:
621 self.finish(target)
622 self.results.satisfied = satisfied
623 return satisfied
624
626 """May or may not be used by the feature. If not used,
627 it will be ignored."""
628 raise NotImplementedError("Override in concrete sub-class")
629
631 raise NotImplementedError("Override in concrete sub-class")
632
634 """Operations to prepare for testing (optional).
635 Override in concrete sub-class if desired"""
636 pass
637
639 """Operations to complete only if evaluate was True (optional).
640 Override in concrete sub-class if desired"""
641 pass
642
645
647 """User-definable by overriding in a sub-class"""
648 pass
649
651 """NOT YET IMPLEMENTED. Would test for reachability of a feature?
652 This is not easy and may be unnecessary!"""
653 return True
654
656 try:
657 self.metric.results = None
658 except AttributeError:
659
660 pass
661
662
663
665 """Abstract super-class for feature leaf nodes.
666 """
670
672 """Update feats and sizes lists in place with metric info, if any.
673 """
674 try:
675 sizes.append(self.metric_len)
676
677 self.metric
678 except AttributeError:
679
680 return
681 else:
682 feats.append(self)
683
685 return "Feature %s"%self.name
686
687 __repr__ = __str__
688
689
691 """Abstract super-class for feature regular nodes (supporting sub-features).
692 """
693 - def __init__(self, name, description='', pars=None,
694 ref_traj=None, subfeatures=None):
695 """Sub-features is an ordered sequence of QL or QT feature instances
696 which are (by default) evaluated in this order on a trajectory segment
697 unless evaluation method is overridden.
698
699 For more sophisticated use of sub-features, they should be provided as
700 a dictionary mapping feature names to the feature instance.
701 """
702 feature.__init__(self, name, description, pars, ref_traj)
703 if subfeatures is None:
704 self.subfeatures = ()
705 self._namemap = {}
706 elif isinstance(subfeatures, (list, tuple)):
707 for sf in subfeatures:
708 assert isinstance(sf, feature), \
709 "Only define quantitative or qualitative features"
710 self.subfeatures = subfeatures
711 self._namemap = dict(zip([sf.name for sf in subfeatures],
712 subfeatures))
713 elif isinstance(subfeatures, dict):
714 for sfname, sf in subfeatures.items():
715 assert isinstance(sf, feature), \
716 "Only define quantitative or qualitative features"
717 self.subfeatures = subfeatures
718 self._namemap = subfeatures
719 else:
720 raise TypeError("List or dictionary of sub-features expected")
721
722
724 """Update feats and sizes lists in place with metric info, if any.
725 """
726 try:
727 sizes.append(self.metric_len)
728
729 self.metric
730 except AttributeError:
731
732 pass
733 else:
734 feats.append(self)
735
736 for sf in self._namemap.values():
737 sf._residual_info(feats, sizes)
738
740 s = "Feature %s "%self.name
741 if len(self._namemap.keys()) > 0:
742 s += "- " + str(self._namemap.keys())
743 return s
744
745 __repr__ = __str__
746
748 """Return named sub-feature"""
749 return self._namemap[featname]
750
752
753 if 'verbose_level' in self.pars:
754 v = max([0,self.pars.verbose_level - 1])
755 if isinstance(sf, common._seq_types):
756 for sf_i in sf:
757 sf_i.pars.verbose_level = v
758 else:
759 sf.pars.verbose_level = v
760 if 'debug' in self.pars:
761 if isinstance(sf, common._seq_types):
762 for sf_i in sf:
763 sf_i.pars.debug = self.pars.debug
764 else:
765 sf.pars.debug = self.pars.debug
766
768 """Default method: evaluate sub-features in order (assumes they
769 are stored as a list).
770
771 Can override with more sophisticated method (e.g. for use with a
772 dictionary of sub-features).
773 """
774
775 satisfied = True
776
777
778 for sf in self.subfeatures:
779 try:
780 self.propagate_verbosity(sf)
781 except KeyboardInterrupt:
782 raise
783 except:
784 if not isinstance(self.subfeatures, common._seq_types):
785 raise TypeError("You must override the evaluate method for "
786 "dictionary-based sub-features")
787 else:
788 raise
789 sf.super_pars.update(self.pars)
790 sf.super_results.update(self.results)
791 sf.reset_metric()
792 if self.pars.verbose_level > 1:
793 print "feature_node.evaluate: sf=", sf
794 error_raised = False
795 try:
796 new_result = sf(target)
797 except KeyboardInterrupt:
798 raise
799 except:
800
801
802 new_result = False
803 error_raised = True
804 if sf.pars.debug:
805 raise
806
807
808 satisfied = satisfied and new_result
809 if error_raised:
810 print " ... error raised"
811 if hasattr(self, 'metric'):
812
813 if sf.metric.results is None:
814 sf.metric.results = self.pars.penalty * \
815 npy.ones((sf.metric_len,),float)
816 else:
817 self.results.update(sf.results)
818 return satisfied
819
821 """May or may not be used by the feature. If not used, it will be
822 ignored."""
823 self.ref_traj = ref_traj
824 self.postprocess_ref_traj()
825 if isinstance(self.subfeatures, dict):
826 sfs = self.subfeatures.values()
827 else:
828 sfs = self.subfeatures
829 for sf in sfs:
830 self.propagate_verbosity(sf)
831 sf.super_pars.update(self.pars)
832 sf.super_results.update(self.results)
833 sf.set_ref_traj(ref_traj)
834 self.pars.update(sf.pars)
835
836
838 """Qualitative feature (leaf node).
839 Add description to note assumptions used for defining feature.
840
841 input: a trajectory segment
842 output: a vector of boolean-valued event detections (non-terminal events or even
843 non-linked python events) or other function tests (e.g. existence of a fixed point)
844 stored in a list.
845 """
846
847
849 """Quantitative feature (leaf node).
850 Add description to note assumptions used for defining feature.
851
852 input: a trajectory segment
853 output: a vector of boolean-valued tolerance tests on the discrepancies between ideal and
854 actual features defined by a list of function tests
855 e.g. a test returns (residual of ideal-actual) < tolerance
856 """
857
858
860 """Qualitative feature (regular node).
861 Add description to note assumptions used for defining feature.
862
863 input: a trajectory segment
864 output: a vector of boolean-valued event detections (non-terminal events or even
865 non-linked python events) or other function tests (e.g. existence of a fixed point)
866 stored in a list.
867 """
868
869
871 """Quantitative feature (regular node).
872 Add description to note assumptions used for defining feature.
873
874 input: a trajectory segment
875 output: a vector of boolean-valued tolerance tests on the discrepancies between ideal and
876 actual features defined by a list of function tests
877 e.g. a test returns (residual of ideal-actual) < tolerance
878 """
879
880
881
882
884 """Model context condition, made up of a boolean composition of wanted and
885 unwanted features.
886 This is specified by a dictionary of feature objects mapping to True
887 (wanted feature) or False (unwanted feature).
888 """
889
890 - def __init__(self, feature_composition_dict):
891
892
893 self.namemap = {}
894 try:
895 for f, c in feature_composition_dict.iteritems():
896 assert isinstance(c, bool), \
897 "Feature composition dictionary requires boolean values"
898 assert isinstance(f, (ql_feature_leaf, qt_feature_leaf,
899 ql_feature_node, qt_feature_node)), \
900 "Only define quantitative or qualitative features"
901 self.namemap[f.name] = f
902 except AttributeError:
903 raise TypeError("Dictionary of features to Booleans expected")
904 self.fcd = feature_composition_dict
905 self.results = common.args()
906
908 try:
909 return self.namemap == other.namemap
910 except AttributeError:
911 return False
912
914 try:
915 return self.namemap != other.namemap
916 except AttributeError:
917 return True
918
920 return self.namemap.keys()
921
923 return self.namemap.values()
924
926 return self.namemap.items()
927
929 return self.namemap[name]
930
932 return name in self.namemap
933
935 """Set reference trajectory for the features (if used, otherwise will
936 be ignored or overridden in feature _local_init methods).
937 """
938 for f,c in self.fcd.iteritems():
939 f.set_ref_traj(ref_traj)
940
942 """Apply conditions to trajectory segments
943 and returns True if all are satisfied."""
944 satisfied = True
945 for f,c in self.fcd.iteritems():
946
947 new = f(target) == c
948 satisfied = satisfied and new
949 self.results[f.name] = f.results
950 return satisfied
951
952 __call__ = evaluate
953
955 s = "Condition "
956 if len(self.namemap.keys()) > 0:
957 s += "- " + str(self.namemap.keys())
958 return s
959
960 __repr__ = __str__
961
963 min_ix = npy.Inf
964 for f in self.fcd.keys():
965 f_ix = f._find_idx()
966 if f_ix is not None and f_ix < min_ix:
967 min_ix = f_ix
968 if npy.isfinite(min_ix):
969 return min_ix
970 else:
971 return None
972
974 """Update metric information used for residual / objective function,
975 from all sub-features."""
976
977 feats = []
978 sizes = []
979 for f in self.fcd.keys():
980 f._residual_info(feats, sizes)
981 return {'features': dict(zip(feats,sizes)),
982 'total_size': sum(sizes)}
983
984 - def collate_results(self, result_name, merge_lists=False,
985 feature_names=None):
986 res = []
987 if feature_names is None:
988 feature_list = self.fcd.keys()
989 else:
990 feature_list = [self.namemap[f] for f in feature_names]
991 for f in feature_list:
992 try:
993 resval = getattr(f.results, result_name)
994 except AttributeError:
995
996 continue
997 else:
998 if merge_lists and isinstance(resval, list):
999 res.extend(resval)
1000 else:
1001 res.append(resval)
1002 return res
1003
1004
1005 -class context(object):
1006 """A collection of related model interfaces that apply to a model.
1007 interface_pairs are a list of ModelInterface instance (test) and class (ref) pairs,
1008 the latter to be instantiated on a model.
1009
1010 Set the debug_mode attribute at any time, or as the optional argument at initializiation,
1011 to ensure that any exceptions that arise from interacting model interfaces and their
1012 features are fully passed back to the caller of the context containing them.
1013 """
1014 - def __init__(self, interface_pairs, debug_mode=False):
1015 self.interfaces = dict(interface_pairs)
1016 self.debug_mode = debug_mode
1017
1018
1019
1020 metric_features = {}
1021 res_feature_list = []
1022 tot_size = 0
1023 for test_mi, ref_mi_class in self.interfaces.iteritems():
1024
1025 metric_features[test_mi] = test_mi.conditions._residual_info()
1026 tot_size += metric_features[test_mi]['total_size']
1027 res_feature_list.extend([(test_mi, f) for f in \
1028 metric_features[test_mi]['features'].keys()])
1029 self.metric_features = metric_features
1030 self.res_feature_list = res_feature_list
1031 self.res_len = tot_size
1032
1033 self.ref_interface_instances = []
1034
1035 self.reset_weights()
1036
1037 - def reset_weights(self, old_weights=None):
1038 """Reset weights to unity, unless old_weights array
1039 is given, in which case reset to that.
1040 """
1041 if old_weights is None:
1042 self.weights = npy.ones(self.res_len, 'float')
1043 else:
1044 self.weights = old_weights
1045 self.weight_index_mapping = {}
1046 self.feat_weights = {}
1047 ix = 0
1048 for test_mi, feat_dict in self.metric_features.iteritems():
1049 self.weight_index_mapping[test_mi] = {}
1050 for f, size in feat_dict['features'].iteritems():
1051 self.weight_index_mapping[test_mi][f] = (ix, ix+size)
1052
1053 self.feat_weights[(test_mi, f)] = self.weights[ix]
1054 ix += size
1055
1056 - def set_single_feat_weights(self, feat, weight=1):
1057 """Set weights for a single feature given as an (interface, feature)
1058 pair, setting all others to zero."""
1059 wdict = {}
1060 for test_mi, feat_dict in self.metric_features.iteritems():
1061 if test_mi != feat[0]:
1062 continue
1063 w = {}.fromkeys(feat_dict['features'].keys(), 0)
1064 if feat[1] in w:
1065 w[feat[1]] = weight
1066 wdict[test_mi] = w
1067 self.set_weights(wdict)
1068
1069 - def set_weights(self, weight_dict):
1070 """Update weights with a dictionary keyed by test_mi, whose values are
1071 either:
1072 (1) dicts of feature -> scalar weight.
1073 (2) a scalar which will apply to all features of that model interface
1074 Features and model interfaces must correspond to those declared for the
1075 context.
1076 """
1077 for test_mi, fs in weight_dict.iteritems():
1078 try:
1079 flist = self.metric_features[test_mi]['features'].keys()
1080 except KeyError:
1081 raise AssertionError("Invalid test model interface")
1082 if isinstance(fs, common._num_types):
1083 feat_dict = {}.fromkeys(flist, fs)
1084 elif isinstance(fs, dict):
1085 assert npy.alltrue([isinstance(w, common._num_types) for \
1086 w in fs.values()]), "Invalid scalar weight"
1087 assert npy.alltrue([f in flist for f in fs.keys()]), \
1088 "Invalid features given for this test model interface"
1089 feat_dict = fs
1090 for f, w in feat_dict.iteritems():
1091 self.feat_weights[(test_mi, f)] = w
1092
1093 start_ix, end_ix = self.weight_index_mapping[test_mi][f]
1094 self.weights[start_ix:end_ix] = w
1095
1096 - def show_res_info(self, resvec):
1097 """Show detail of feature -> residual mapping for a given residual
1098 vector."""
1099 i = 0
1100 for test_mi, feat_dict in self.metric_features.iteritems():
1101 print "Test model interface:", test_mi
1102 for f in feat_dict['features']:
1103 if self.feat_weights[(test_mi, f)] == 0:
1104 continue
1105 ix0, ix1 = self.weight_index_mapping[test_mi][f]
1106 len_w = ix1-ix0
1107 f_str = " "+f.name
1108
1109 extra_space_w = " "*max([0, 13-len(f_str)])
1110 extra_space_unw = " "*max([0, len(f_str)-13])
1111 print f_str + extra_space_w , resvec[i:i+len_w]
1112 try:
1113 print " unweighted:" + extra_space_unw, \
1114 resvec[i:i+len_w]/self.weights[ix0:ix1]
1115 except ZeroDivisionError:
1116 print " (unweighted values unavailable)"
1117 i += len_w
1118
1119 - def _map_to_features(self, x):
1120 """Utility to map 1D array x onto the model interface's
1121 features with non-zero weights, returning a dictionary.
1122
1123 x is assumed to have correct length.
1124 """
1125 out = {}
1126 i = 0
1127 for test_mi, feat_dict in self.metric_features.iteritems():
1128 for f in feat_dict['features']:
1129 if self.feat_weights[(test_mi, f)] == 0:
1130 continue
1131 ix0, ix1 = self.weight_index_mapping[test_mi][f]
1132 len_w = ix1-ix0
1133 try:
1134 out[test_mi][f] = x[i:i+len_w]
1135 except KeyError:
1136 out[test_mi] = {f: x[i:i+len_w]}
1137 i += len_w
1138 return out
1139
1140 - def evaluate(self, model):
1141 """Evaluate whole context on a model instance, returning a single
1142 Boolean.
1143 """
1144 result = True
1145
1146
1147 self.ref_interface_instances = []
1148 for test_mi, ref_mi_class in self.interfaces.iteritems():
1149
1150 ref_mi = ref_mi_class(model)
1151 self.ref_interface_instances.append(ref_mi)
1152 try:
1153 new_result = test_mi(ref_mi)
1154 except KeyboardInterrupt:
1155 raise
1156 except:
1157 if self.debug_mode:
1158 raise
1159 else:
1160 print "******************************************"
1161 print "Problem evaluating interface", test_mi, "on ", ref_mi
1162 print " ", sys.exc_info()[0], sys.exc_info()[1]
1163 new_result = False
1164
1165
1166 result = result and new_result
1167 return result
1168
1169 - def residual(self, model, include_raw=False):
1170 """Evaluate whole context on a model instance, returning an array
1171 of residual error between quantitative features in the model trajectory
1172 and their target values.
1173
1174 Residual array will be weighted if one was set. Any weights set to zero
1175 will cause those features to *not appear* in the residual.
1176
1177 Provide include_raw=True argument to also return the raw, unweighted residual.
1178 (Mainly for internal use.)
1179 """
1180
1181
1182 self.evaluate(model)
1183 raw_residual = npy.concatenate(tuple([mf[1].metric.results for \
1184 mf in self.res_feature_list]))
1185 residual = process_raw_residual(raw_residual, self.weights)
1186 if include_raw:
1187 return residual, raw_residual
1188 else:
1189 return residual
1190
1191
1193 ixs = npy.nonzero(weights)
1194 residual = (weights*raw_residual)[ixs]
1195 nan_ixs = npy.where(npy.asarray(npy.isnan(residual),int))
1196 for ix in nan_ixs:
1197 residual[ix] = 100.
1198 return residual
1199
1200
1201
1202
1204 """Use this for a single vector field model that uses discrete
1205 event mappings."""
1208
1210 """Use this as a binary switch feature, toggled
1211 by a given variable name 'varname' that is supplied
1212 in the pars dict at initialization."""
1214 try:
1215 pts = target.test_traj.sample(coords=[self.pars.varname])
1216 except AttributeError:
1217 raise AttributeError("No variable name given for switch")
1218 except KeyboardInterrupt:
1219 raise
1220 except:
1221 print "Failed to find trajectory values for given variable name: %s"%self.pars.varname
1222 raise
1223 self.results.output = pts
1224 return all(self.results.output==1)
1225
1227 if self.results.satisfied:
1228
1229 return None
1230 res = self.results.output
1231 if res[0] == 1:
1232 adjusted_res = list((res - 1) != 0)
1233 else:
1234 if 1 not in res:
1235
1236 raise RuntimeError
1237 adjusted_res = list(res != 0)
1238
1239
1240
1241 return adjusted_res.index(True)
1242
1243
1244
1245
1247 """Generic and abstract interface class for dynamical systems."""
1248
1249
1250
1251
1252
1253 _setkeys = ['globalt0', 'tdata', 'pars', 'algparams', 'inputs']
1254
1255 _querykeys = ['pars', 'parameters', 'events', 'submodels',
1256 'ics', 'initialconditions', 'vars', 'variables',
1257 'auxvariables', 'auxvars', 'vardomains', 'abseps']
1258
1260 raise NotImplementedError("Only call this on a concrete sub-class")
1261
1262 - def query(self, querykey=''):
1263 return self.model.query(querykey)
1264
1265
1267 """Wrapper for Generator (for non-hybrid models) that shares similar API
1268 with ModelInterface for use in HybridModel objects."""
1269
1270 - def __init__(self, model, FScompatibleNames=None,
1271 FScompatibleNamesInv=None):
1272 """model argument must be a Generator only"""
1273 self.model = model
1274 if FScompatibleNames is None:
1275 if self.model._FScompatibleNames is None:
1276 self.model._FScompatibleNames = symbolMapClass()
1277 else:
1278 self.model._FScompatibleNames = FScompatibleNames
1279 if FScompatibleNamesInv is None:
1280 if self.model._FScompatibleNamesInv is None:
1281 self.model._FScompatibleNamesInv = symbolMapClass()
1282 else:
1283 self.model._FScompatibleNamesInv = FScompatibleNamesInv
1284 self.eventstruct = Events.EventStruct()
1285
1286
1287 - def get(self, key, ics=None, t0=0):
1288
1289 return self.model.get(key)
1290
1291 - def set(self, key, value, ics=None, t0=0):
1292 if key in self._setkeys:
1293 self.model.set(**{key:value})
1294 else:
1295 raise KeyError("Invalid or unsupported 'set' key: %s"%key)
1296
1297 - def Rhs(self, t, xdict, pdict):
1298 """Direct access to a generator's Rhs function."""
1299 return self.model.Rhs(t, xdict, pdict)
1300
1301 - def Jacobian(self, t, xdict, pdict, idict=None):
1302 """Direct access to a generator's Jacobian function (if defined)."""
1303 return self.model.Jacobian(t, xdict, pdict)
1304
1306 """Direct access to a generator's JacobianP function (if defined)."""
1307 return self.model.JacobianP(t, xdict, pdict)
1308
1310 """Direct access to a generator's MassMatrix function (if defined)."""
1311 return self.model.MassMatrix(t, xdict, pdict)
1312
1313 - def AuxVars(self, t, xdict, pdict):
1314 """Direct access to a generator's auxiliary variables
1315 definition (if defined)."""
1316 return self.model.AuxVars(t, xdict, pdict)
1317
1318
1320 """Model constraints expressed as a uni-directional interface to another
1321 formal system model:
1322 - Made up of conditions imposed on the other system's test trajectory.
1323 - Defines evaluation criteria for any view (e.g. from experimental data and
1324 test conditions).
1325 This is an abstract superclass for the 'internal' and 'external'
1326 sub-classes.
1327 """
1328 _trajname = 'test_iface_traj'
1329
1334
1335
1337 """initiator is a ModelInterface or GeneratorInterface object"""
1338 if self._initiator_cache is None:
1339 if ics is None:
1340 raise ValueError("Must pass initial conditions")
1341 else:
1342 initiator = self.model._findTrajInitiator(None, 0,
1343 t0, dict(ics))[0]
1344 self._initiator_cache = (ics, t0, initiator)
1345 else:
1346 if npy.alltrue(self._initiator_cache[0] == ics) and \
1347 self._initiator_cache[1] == t0:
1348 ics, t0, initiator = self._initiator_cache
1349 elif ics is None:
1350 raise ValueError("Must pass initial conditions")
1351 else:
1352
1353 initiator = self.model._findTrajInitiator(None, 0,
1354 t0, dict(ics))[0]
1355 self._initiator_cache = (ics, t0, initiator)
1356 return (ics, t0, initiator)
1357
1358 - def set(self, key, value, ics=None, t0=0):
1359
1360 self.model.set(**{key:value})
1361 ics, t0, initiator = self._get_initiator_cache(ics, t0)
1362
1363 if key in self._setkeys:
1364 initiator.set(key, value, ics, t0)
1365 initiator.model.set(**{key:value})
1366 else:
1367 raise KeyError("Invalid or unsupported 'set' key: %s"%key)
1368
1369 - def get(self, key, ics=None, t0=0):
1370 ics, t0, initiator = self._get_initiator_cache(ics, t0)
1371
1372 try:
1373 return initiator.get(key, ics, t0)
1374 except AttributeError:
1375 raise ValueError("Invalid or unsupported 'get' key: %s"%key)
1376
1377 - def Rhs(self, t, xdict, pdict):
1378 """Direct access to a generator's Rhs function."""
1379 ics_ignore, t_ignore, ds = self._get_initiator_cache(xdict, t)
1380 try:
1381 return self.model.Rhs(ds._supermodel.name, t, xdict, pdict)
1382 except AttributeError:
1383
1384 return self.model.Rhs(t, xdict, pdict)
1385
1386 - def Jacobian(self, t, xdict, pdict, idict=None):
1387 """Direct access to a generator's Jacobian function (if defined)."""
1388 ics_ignore, t_ignore, ds = self._get_initiator_cache(xdict, t)
1389 try:
1390 return self.model.Jacobian(ds._supermodel.name, t, xdict, pdict)
1391 except AttributeError:
1392
1393 return self.model.Jacobian(t, xdict, pdict)
1394
1396 """Direct access to a generator's JacobianP function (if defined)."""
1397 ics_ignore, t_ignore, ds = self._get_initiator_cache(xdict, t)
1398 try:
1399 return self.model.JacobianP(ds._supermodel.name, t, xdict, pdict)
1400 except AttributeError:
1401
1402 return self.model.JacobianP(t, xdict, pdict)
1403
1405 """Direct access to a generator's MassMatrix function (if defined)."""
1406 ics_ignore, t_ignore, ds = self._get_initiator_cache(xdict, t)
1407 try:
1408 return self.model.MassMatrix(ds._supermodel.name, t, xdict, pdict)
1409 except AttributeError:
1410
1411 return self.model.MassMatrix(t, xdict, pdict)
1412
1413 - def AuxVars(self, t, xdict, pdict):
1414 """Direct access to a generator's auxiliary variables
1415 definition (if defined)."""
1416 ics_ignore, t_ignore, ds = self._get_initiator_cache(xdict, t)
1417 try:
1418 return ds.model.AuxVars(ds._supermodel.name, t, xdict, pdict)
1419 except AttributeError:
1420
1421 return ds.model.AuxVars(t, xdict, pdict)
1422
1424
1425
1426 if conditions is None:
1427 self.conditions = None
1428 else:
1429
1430 self.conditions = conditions
1431 try:
1432 self.conditions.set_ref_traj(traj)
1433 except AttributeError:
1434 raise
1435
1436 - def evaluate(self, target, force=False):
1437 """Evaluate interface consistency against target interface's trajectory
1438 on specified conditions.
1439
1440 Optional force argument forces model to recompute its test trajectory,
1441 e.g. because of a known change in model parameters, ics, etc.
1442 """
1443 assert isinstance(target, ModelInterface), \
1444 "Target argument must be another interface object"
1445 if len(self.compatibleInterfaces) > 0 and \
1446 target.__class__.__name__ not in self.compatibleInterfaces \
1447 and not npy.sometrue([common.compareBaseClass(target, ctype) \
1448 for ctype in self.compatibleInterfaces]):
1449 raise ValueError("Target interface not of compatible type")
1450 try:
1451 self.conditions
1452 except AttributeError:
1453 self.setup_conditions(conditions, self.get_test_traj())
1454 force = force or target.test_traj is None
1455 if force:
1456
1457 target.get_test_traj(force=force)
1458 self.prepare_conditions(target)
1459 try:
1460 result = self.conditions(target)
1461 except KeyError:
1462 raise KeyError("Condition evaluation failed")
1463 return result
1464
1465 __call__ = evaluate
1466
1467 - def postprocess_test_traj(self, test_traj):
1468 """Called by another interface via get_test_traj.
1469 Default post-processing of test trajectory is the identity
1470 function, i.e. no processing.
1471
1472 Override this method to return a processed version of the
1473 trajectory or perform other post-computation clean-up, e.g.
1474 prepare auxiliary feature/condition-related information based
1475 on end state of trajectory so that HybridModel can use it to
1476 decide on next hybrid state to switch to.
1477 """
1478 return test_traj
1479
1481 """Called automatically by evaluate. Override with user-defined access
1482 to the target interface or processing of trajectory after return of the
1483 target's test trajectory.
1484 """
1485 pass
1486
1487
1489 """Interface providing internal evaluation criteria between models.
1490 Optional conditions (object) argument used to specify these criteria.
1491 """
1492 - def __init__(self, model, conditions=None, compatibleInterfaces=None,
1493 test_traj=None):
1494 """Set model that generates test trajectories from which the dictionary
1495 of conditions can be imposed on a connected model.
1496
1497 If no conditions are specified then the model is trivially wrapped in
1498 an "empty" interface.
1499
1500 Optionally, a dummy test traj can be supplied in case of a dummy interface
1501 for a trivial condition test that does not need to evaluate a trajectory
1502 to determine the result.
1503 """
1504 ModelInterface.__init__(self)
1505 assert isinstance(model, Model.Model), "Invalid Model object passed"
1506 self.model = model
1507
1508 self.test_traj = test_traj
1509
1511 """Cause recomputation of test trajectory if not already present in
1512 model, returning boolean for whether recomputation was performed.
1513 """
1514 info = self.model.current_defining_args()
1515
1516 new_args = self.initialize_model()
1517 if new_args is not None:
1518 info.update(new_args)
1519 if self.model.has_exact_traj(self._trajname, info):
1520
1521
1522 return False
1523 else:
1524 try:
1525 self.compute_traj(need_init=False, new_args=new_args)
1526 except KeyboardInterrupt:
1527 raise
1528 except:
1529 print "Model interface compute_traj method for model " + \
1530 "'%s' failed" % self.model.name
1531 print sys.exc_info()[0], sys.exc_info()[1]
1532 return False
1533 else:
1534 return True
1535
1537 return self.test_traj is not None
1538
1540 if need_init:
1541 new_args = self.initialize_model()
1542 if new_args is not None and len(new_args) > 0:
1543 old_info = self.model.current_defining_args()
1544 self.model.set(**new_args)
1545 else:
1546 old_info = None
1547 self.model.compute(trajname=self._trajname, force=True)
1548 if old_info is not None:
1549
1550 self.model.set(**dict(old_info))
1551
1553 """Called by another interface.
1554 Return model's test trajectory, using any post-processing
1555 specified by user-defined process_test_traj method.
1556
1557 Use force option if model is known to have changed and trajectory
1558 needs refreshing.
1559 """
1560 if force and not isinstance(self.test_traj, Trajectory.Trajectory):
1561 self.compute_traj()
1562 recomputed = True
1563 else:
1564 recomputed = self.ensure_has_test_traj()
1565 if recomputed or self.test_traj is None:
1566 self.test_traj = \
1567 self.postprocess_test_traj(self.model[self._trajname])
1568 return self.test_traj
1569
1571 """Return any unique model-specific settings here, as a dictionary with
1572 keys that can include initial conditions, parameters, tdata, algorithmic
1573 parameters. Use the same keys that are suitable for a call to the
1574 Model.set method, i.e. 'pars', 'ics', 'tdata', and 'algparams'.
1575
1576 Override in a sub-class to use. This method will be called
1577 before any trajectory computation of the model.
1578 """
1579 pass
1580
1581
1583 """Interface from a trajectory of numerical data and test conditions
1584 providing external evaluation criteria for a model.
1585 Optional conditions (object) argument used to specify these criteria.
1586 """
1587 - def __init__(self, traj=None, conditions=None, compatibleInterfaces=None):
1595
1597 """Do any user-defined preprocessing to the given trajectory, including
1598 converting it to a different type of trajectory.
1599 """
1600 self.test_traj = self.postprocess_test_traj(traj)
1601
1602 self.conditions.set_ref_traj(self.test_traj)
1603
1605 """Never needs to recompute trajectory as it is fixed, so always
1606 returns False.
1607 """
1608 assert self.has_test_traj(), "Test trajectory missing"
1609 return False
1610
1613
1615 """Called by another interface.
1616 Optional force argument is ignored for this class, as the
1617 trajectory is fixed."""
1618 return self.test_traj
1619
1620
1621
1622
1623
1624
1625
1626
1627
1628
1629
1630
1631
1632
1633
1634
1635 -class MReg(object):
1636 """Registry class for Model descriptors and instances in Model projects.
1637 For internal use by PyDSTool."""
1638
1640
1641 self.descs = {}
1642
1643 self.instances = {}
1644
1645 - def add(self, descriptor):
1646 """descriptor expected to be an MDescriptor object.
1647 """
1648 if isinstance(descriptor, ModelConstructor.MDescriptor):
1649 self.descs[descriptor.name] = descriptor
1650 self.instances[descriptor.name] = {}
1651 else:
1652 raise TypeError("Only MDescriptor objects valid for MReg class")
1653
1655 try:
1656 self.instances[name] = model_instance
1657 except KeyError:
1658 raise ValueError("No such model descriptor")
1659
1661 return name in self.descs
1662
1664 return self.descs[name]
1665
1667 del(self.descs[name])
1668 del(self.instances[name])
1669
1670 remove = __delitem__
1671
1672 - def query(self, querykey, value):
1673 """Return info about stored model specifications.
1674 Valid query keys: 'orig_name', 'in_description'
1675 """
1676 assert isinstance(querykey, str), \
1677 ("Query argument must be a single string")
1678 _keylist = ['orig_name', 'in_description']
1679 if querykey not in _keylist:
1680 print 'Valid query keys are:', _keylist
1681 raise TypeError('Query key '+querykey+' is not valid')
1682 if querykey == 'orig_name':
1683 res = []
1684 for name, regentry in self.descs.iteritems():
1685 if regentry.orig_name == value:
1686 res.append(name)
1687 return res
1688 if querykey == 'in_description':
1689 res = []
1690 for name, regentry in self.descs.iteritems():
1691 if value in regentry.description:
1692 res.append(name)
1693 return res
1694