Changeset 8127:3d677acb367c in orange


Ignore:
Timestamp:
07/29/11 14:29:29 (3 years ago)
Author:
ales_erjavec <ales.erjavec@…>
Branch:
default
Convert:
14f9dba1ec2203e76fbc436d9d2d67e85474a4f4
Message:

Made the multi label EarthLearner/Classifier the default (removed the ML postfix).
Commented out the old Learner/Classifier pair.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • orange/Orange/regression/earth.py

    r8121 r8127  
    11"""\ 
    2 EARTH (Multivariate adaptive regression splines - MARS) (``earth``) 
    3  
    4   
    5 .. autoclass :: EarthLearner 
    6  
    7 .. autoclass :: EarthClassifier 
     2==================================================== 
     3Multivariate Adaptive Regression Splines (``earth``) 
     4==================================================== 
     5 
     6`Multivariate adaptive regression splines (MARS)`_ is a non-parametric 
     7regression method that extends a linear model with non-linear 
     8interactions. 
     9 
     10This module borrows the implementation of the technique from the `Earth R  
     11package`_ by Stephen Milborrow.  
     12 
     13.. _`Multivariate adaptive regression splines (MARS)`: http://en.wikipedia.org/wiki/Multivariate_adaptive_regression_splines 
     14.. _`Earth R package`: http://cran.r-project.org/web/packages/earth/index.html 
     15 
     16.. autoclass:: EarthLearner 
     17 
     18.. autoclass:: EarthClassifier 
     19 
     20.. autoclass:: EarthLearnerML 
     21 
     22.. autoclass:: EarthClassifierML 
    823 
    924""" 
     
    1530import numpy 
    1631 
    17 from Orange.misc.render import contextmanager 
    18  
    19 @contextmanager  
    20 def member_set(obj, name, val): 
    21     """ A context manager that sets member ``name`` on ``obj`` to ``val`` 
    22     and restores the previous value on exit.  
    23     """ 
    24     old_val = getattr(obj, name, val) 
    25     setattr(obj, name, val) 
    26     yield 
    27     setattr(obj, name, old_val) 
    28      
    29 def base_matrix(data, best_set, dirs, cuts): 
    30     """ Return the base matrix for the earth model. 
    31      
    32     """ 
    33     data = numpy.asarray(data) 
    34     best_set = numpy.asarray(best_set) 
    35     dirs = numpy.asarray(dirs) 
    36     cuts = numpy.asarray(cuts) 
    37      
    38     bx = numpy.zeros((data.shape[0], best_set.shape[0])) 
    39     bx[:, 0] = 1.0 # The intercept 
    40     for termi in range(1, best_set.shape[0]): 
    41         term_dirs = dirs[termi] 
    42         term_cuts = cuts[termi] 
    43          
    44         dir_p1 = numpy.where(term_dirs == 1)[0] 
    45         dir_m1 = numpy.where(term_dirs == -1)[0] 
    46         dir_2 = numpy.where(term_dirs == 2)[0] 
    47          
    48         x1 = data[:, dir_p1] - term_cuts[dir_p1] 
    49         x2 = term_cuts[dir_m1] - data[:, dir_m1] 
    50         x3 = data[:, dir_2] 
    51          
    52         x1 = numpy.where(x1 > 0.0, x1, 0.0) 
    53         x2 = numpy.where(x2 > 0.0, x2, 0.0) 
    54          
    55         X = numpy.hstack((x1, x2, x3)) 
    56         X = numpy.cumprod(X, axis=1) 
    57         bx[:, termi] = X[:, -1] if X.size else 0.0 
    58          
    59     return bx 
    60  
    61  
    62 class EarthLearner(BaseEarthLearner): 
    63     """  
    64      
    65     """ 
    66     def __new__(cls, data=None, weightId=None, **kwargs): 
    67         self = BaseEarthLearner.__new__(cls, **kwargs) 
    68         if data is not None: 
    69             self.__init__(**kwargs) 
    70             return self.__call__(data, weightId) 
    71         return self 
    72      
    73     def __init__(self, store_examples=True, **kwargs): 
    74         self.store_examples = store_examples 
    75         for key, val in kwargs.items(): 
    76             setattr(self, key, val) 
    77      
    78     def __call__(self, data, weightId=None): 
    79         if not data.domain.class_var: 
    80             raise ValueError("No class var in the domain.") 
    81          
    82         with member_set(self, "prune", False): 
    83             # We overwrite the prune argument (will do the pruning in python). 
    84             base_clsf =  BaseEarthLearner.__call__(self, data, weightId) 
    85          
    86         if self.prune: 
    87             (best_set, betas, rss, subsets, rss_per_subset, 
    88              gcv_per_subset) = self.pruning_pass(base_clsf, data) 
    89              
    90             return EarthClassifier(base_clsf, data if self.store_examples else None, 
    91                                    best_set=best_set, dirs=base_clsf.dirs, 
    92                                    cuts=base_clsf.cuts, 
    93                                    betas=betas, 
    94                                    subsets=subsets, 
    95                                    rss_per_subset=rss_per_subset, 
    96                                    gcv_per_subset=gcv_per_subset) 
    97         else: 
    98             return EarthClassifier(base_clsf, data if self.store_examples else None) 
    99      
    100      
    101     def pruning_pass(self, base_clsf, examples): 
    102         """ Prune the terms constructed in the forward pass. 
    103         (Pure numpy reimplementation) 
    104         """ 
    105         n_terms = numpy.sum(base_clsf.best_set) 
    106         n_eff_params = n_terms + self.penalty * (n_terms - 1) / 2.0 
    107         data, y, _ = examples.to_numpy_MA() 
    108         data = data.filled(0.0) 
    109         best_set = numpy.asarray(base_clsf.best_set, dtype=bool) 
    110          
    111         bx = base_matrix(data, base_clsf.best_set, 
    112                          base_clsf.dirs, base_clsf.cuts, 
    113                          ) 
    114          
    115         bx_used = bx[:, best_set] 
    116         subsets, rss_per_subset = subsets_selection_xtx(bx_used, y) # TODO: Use leaps like library 
    117         gcv_per_subset = [gcv(rss, bx.shape[0], i + self.penalty * (i - 1) / 2.0) \ 
    118                               for i, rss in enumerate(rss_per_subset, 1)] 
    119         gcv_per_subset = numpy.array(gcv_per_subset) 
    120          
    121         best_i = numpy.argmin(gcv_per_subset[1:]) + 1 # Ignore the intercept 
    122         best_ind = subsets[best_i, :best_i + 1] 
    123         bs_i = 0 
    124         for i, b in enumerate(best_set): 
    125             if b: 
    126                 best_set[i] = bs_i in best_ind 
    127                 bs_i += 1 
    128                  
    129         bx_subset = bx[:, best_set] 
    130         betas, rss, rank, s = numpy.linalg.lstsq(bx_subset, y) 
    131         return best_set, betas, rss, subsets, rss_per_subset, gcv_per_subset 
    132      
    133          
    134 class EarthClassifier(Orange.core.ClassifierFD): 
    135     def __init__(self, base_classifier=None, examples=None, best_set=None, 
    136                  dirs=None, cuts=None, betas=None, subsets=None, 
    137                  rss_per_subset=None, 
    138                  gcv_per_subset=None): 
    139         self._base_classifier = base_classifier 
    140         self.examples = examples 
    141         self.domain = base_classifier.domain 
    142         self.class_var = base_classifier.class_var 
    143          
    144         best_set = base_classifier.best_set if best_set is None else best_set 
    145         dirs = base_classifier.dirs if dirs is None else dirs 
    146         cuts = base_classifier.cuts if cuts is None else cuts 
    147         betas = base_classifier.betas if betas is None else betas 
    148          
    149         self.best_set = numpy.asarray(best_set, dtype=bool) 
    150         self.dirs = numpy.array(dirs, dtype=int) 
    151         self.cuts = numpy.array(cuts) 
    152         self.betas = numpy.array(betas) 
    153          
    154         self.subsets = subsets 
    155         self.rss_per_subset = rss_per_subset 
    156         self.gcv_per_subset = gcv_per_subset 
    157          
    158     @property 
    159     def num_terms(self): 
    160         """ Number of terms in the model (including the intercept). 
    161         """ 
    162         return numpy.sum(numpy.asarray(self.best_set, dtype=int)) 
    163      
    164     @property 
    165     def max_terms(self): 
    166         """ Maximum number of terms (as specified in the learning step). 
    167         """ 
    168         return self.best_set.size 
    169      
    170     @property 
    171     def num_preds(self): 
    172         """ Number of predictors (variables) included in the model. 
    173         """ 
    174         return len(self.used_attributes(term)) 
    175      
    176     def __call__(self, example, what=Orange.core.GetValue): 
    177         value = self.predict(example) 
    178         if isinstance(self.class_var, Orange.data.variable.Continuous): 
    179             value = self.class_var(value) 
    180         else: 
    181             value = self.class_var(int(round(value))) 
    182              
    183         dist = Orange.statistics.distribution.Distribution(self.class_var) 
    184         dist[value] = 1.0 
    185         if what == Orange.core.GetValue: 
    186             return value 
    187         elif what == Orange.core.GetProbabilities: 
    188             return dist 
    189         else: 
    190             return (value, dist) 
    191      
    192     def base_matrix(self, examples=None): 
    193         """ Return the base matrix (bx) of the Earth model for the table. 
    194         If table is not supplied the base matrix of the training examples  
    195         is returned. 
    196          
    197          
    198         :param examples: Input examples for the base matrix. 
    199         :type examples: Orange.data.Table  
    200          
    201         """ 
    202         if examples is None: 
    203             examples = self.examples 
    204              
    205         if examples is None: 
    206             raise ValueError("base matrix is only available if 'store_examples=True'") 
    207          
    208         if isinstance(examples, Orange.data.Table): 
    209             data, cls, w = examples.to_numpy_MA() 
    210             data = data.filled(0.0) 
    211         else: 
    212             data = numpy.asarray(examples) 
    213              
    214         return base_matrix(data, self.best_set, self.dirs, self.cuts) 
    215      
    216     def _anova_order(self): 
    217         """ Return indices that sort the terms into the 'ANOVA' format. 
    218         """ 
    219         terms = [([], 0)] # intercept 
    220         for i, used in enumerate(self.best_set[1:], 1): 
    221             attrs = sorted(self.used_attributes(i)) 
    222             if used and attrs: 
    223                 terms.append((attrs, i)) 
    224         terms = sotred(terms, key=lambda t:(len(t[0]), t[0])) 
    225         return [i for _, i in terms] 
    226      
    227     def format_model(self, percision=3, indent=3): 
    228         return format_model(self, percision, indent) 
    229      
    230     def print_model(self, percision=3, indent=3): 
    231         print self.format_model(percision, indent) 
    232          
    233     def predict(self, ex): 
    234         """ Return the predicted value (float) for example. 
    235         """ 
    236         x = Orange.data.Table(self.domain, [ex]) 
    237         x, c, w = x.to_numpy_MA() 
    238         x = x.filled(0.0)[0] 
    239          
    240         bx = numpy.ones(self.best_set.shape) 
    241         betas = numpy.zeros_like(self.betas) 
    242         betas[0] = self.betas[0] 
    243         beta_i = 0 
    244         for termi in range(1, len(self.best_set)): 
    245             dirs = self.dirs[termi] 
    246             cuts = self.cuts[termi] 
    247             dir_p1 = numpy.where(dirs == 1)[0] 
    248             dir_m1 = numpy.where(dirs == -1)[0] 
    249             dir_2 = numpy.where(dirs == 2)[0] 
    250              
    251             x1 = x[dir_p1] - cuts[dir_p1] 
    252             x2 = cuts[dir_m1] - x[dir_m1] 
    253             x3 = x[dir_2] 
    254              
    255             x1 = numpy.maximum(x1, 0.0) 
    256             x2 = numpy.maximum(x2, 0.0) 
    257  
    258             X = numpy.hstack((x1, x2, x3)) 
    259             X = numpy.cumprod(X) 
    260              
    261             bx[termi] = X[-1] if X.size else 0.0 
    262             if self.best_set[termi] != 0: 
    263                 beta_i += 1 
    264                 betas[beta_i] = self.betas[beta_i] 
    265  
    266         return numpy.sum(bx[self.best_set] * betas) 
    267              
    268     def used_attributes(self, term=None): 
    269         """ Return a list of used attributes. If term (index) is given 
    270         return only attributes used in that single term. 
    271          
    272         """ 
    273         if term is None: 
    274             terms = numpy.where(self.best_set)[0] 
    275         else: 
    276             terms = [term] 
    277         attrs = set() 
    278         for ti in terms: 
    279             attri = numpy.where(self.dirs[ti] != 0.0)[0] 
    280             attrs.update([self.domain.attributes[i] for i in attri]) 
    281         return attrs 
    282          
    283     def evimp(self, used_only=True): 
    284         """ Return the estimated variable importance. 
    285         """ 
    286         return evimp(self, used_only) 
    287          
    288     def __reduce__(self): 
    289         return (EarthClassifier, (self._base_classifier, self.examples, 
    290                                   self.best_set, self.dirs, self.cuts, 
    291                                   self.betas, self.subsets, 
    292                                   self.rss_per_subset, self.gcv_per_subset), 
    293                 {}) 
    294                                   
    295      
    296 def gcv(rss, n, n_effective_params): 
    297     """ Return the generalized cross validation. 
    298      
    299     .. math: gcv = rss / (n * (1 - n_effective_params / n) ^ 2) 
    300      
    301     :param rss: Residual sum of squares. 
    302     :param n: Number of training examples. 
    303     :param n_effective_params: Number of effective paramaters. 
    304       
    305     """ 
    306     return  rss / (n * (1 - n_effective_params / n) ** 2) 
    307      
    308  
    309 def subsets_selection_xtx(X, Y): 
    310     """ A numpy implementation of EvalSubsetsUsingXtx in the Earth package.  
    311     """ 
    312     X = numpy.asarray(X) 
    313     Y = numpy.asarray(Y) 
    314      
    315     var_count= X.shape[1] 
    316     rss_vec = numpy.zeros(var_count) 
    317     working_set = range(var_count) 
    318     subsets = numpy.zeros((var_count, var_count), dtype=int) 
    319      
    320     for subset_size in reversed(range(var_count)): 
    321         subsets[subset_size, :subset_size + 1] = working_set 
    322         X_work = X[:, working_set] 
    323         b, res, rank, s = numpy.linalg.lstsq(X_work, Y) 
    324         if res.size > 0: 
    325             rss_vec[subset_size] = numpy.sum(res) 
    326         else: 
    327             rss_vec[subset_size] = numpy.sum((Y - numpy.dot(X_work, b)) ** 2) 
    328              
    329         XtX = numpy.dot(X_work.T, X_work) 
    330         iXtX = numpy.linalg.pinv(XtX) 
    331         diag = numpy.diag(iXtX) 
    332          
    333         if subset_size == 0: 
    334             break 
    335          
    336         delta_rss = b ** 2 / diag 
    337         delete_i = numpy.argmin(delta_rss[1:]) + 1 # Keep the intercept 
    338         del working_set[delete_i] 
    339     return subsets, rss_vec 
    340  
    341  
    342 """ 
    343 Multilabel 
    344 """ 
    345  
    346 def is_label_attr(attr): 
    347     """ Is attribute a label. 
    348     """ 
    349     return attr.attributes.has_key("label") 
    350      
    351 def data_label_indices(domain): 
    352     """ Return the indices of label attributes in data. 
    353     """ 
    354     return numpy.where(data_label_mask(domain))[0] 
    355  
    356 def data_label_mask(domain): 
    357     """ Return an array of booleans indicating whether a variable in the 
    358     domain is a label. 
    359     """ 
    360     is_label = map(is_label_attr, domain.variables) 
    361     if domain.class_var: 
    362         is_label[-1] = True 
    363     return numpy.array(is_label, dtype=bool) 
    364  
    365 class EarthLearnerML(Orange.core.LearnerFD): 
    366     """ Multilabel (multi response) earth learner. 
     32class EarthLearner(Orange.core.LearnerFD): 
     33    """ Earth learner class. 
    36734    """ 
    36835    def __new__(cls, examples=None, weight_id=None, **kwargs): 
     
    37643    def __init__(self, degree=1, terms=21, penalty= None, thresh=1e-3, 
    37744                 min_span=0, new_var_penalty=0, fast_k=20, fast_beta=1, 
    378                  pruned_terms=None, scale_resp=False, store_examples=True, **kwds): 
     45                 pruned_terms=None, scale_resp=False, store_examples=True, 
     46                 multi_label=False, **kwds): 
    37947        """  
    38048        .. todo:: min_span, prunning_method 
    381          
    38249        """ 
    38350        self.degree = degree 
     
    39461        self.scale_resp = scale_resp 
    39562        self.store_examples = store_examples 
     63        self.multi_label = multi_label 
    39664        self.__dict__.update(kwds) 
    39765         
    39866    def __call__(self, examples, weight_id=None): 
    399         label_mask = data_label_mask(examples.domain) 
     67        if self.multi_label: 
     68            label_mask = data_label_mask(examples.domain) 
     69        else: 
     70            label_mask = numpy.zeros(len(examples.domain.variables), 
     71                                     dtype=bool) 
     72            label_mask[-1] = bool(examples.domain.class_var) 
     73             
    40074        if not any(label_mask): 
    40175            raise ValueError("The domain has no response variable.") 
     76          
    40277        data = examples.to_numpy_MA("Ac")[0] 
    40378        y = data[:, label_mask] 
     
    40580         
    40681        # TODO: y scaling 
    407         n_terms, used, bx, dirs, cuts = _forward_pass(x, y, 
     82        n_terms, used, bx, dirs, cuts = forward_pass(x, y, 
    40883            degree=self.degree, terms=self.terms, penalty=self.penalty, 
    40984            thresh=self.thresh, fast_k=self.fast_k, fast_beta=self.fast_beta, 
     
    42499        betas, res, rank, s = numpy.linalg.lstsq(bx_used, y) 
    425100         
    426         return EarthClassifierML(examples.domain, used, dirs, cuts, betas.T, 
    427                                  subsets, rss_per_subset, gcv_per_subset, 
    428                                  examples=examples if self.store_examples else None, 
    429                                  label_mask=label_mask) 
    430      
    431  
    432 class EarthClassifierML(Orange.core.ClassifierFD): 
    433     """ Multi label Earth classifier. 
     101        return EarthClassifier(examples.domain, used, dirs, cuts, betas.T, 
     102                               subsets, rss_per_subset, gcv_per_subset, 
     103                               examples=examples if self.store_examples else None, 
     104                               label_mask=label_mask, multi_flag=self.multi_label) 
     105     
     106 
     107class EarthClassifier(Orange.core.ClassifierFD): 
     108    """ Earth classifier. 
    434109    """ 
    435110    def __init__(self, domain, best_set, dirs, cuts, betas, subsets=None, 
    436111                 rss_per_subset=None, gcv_per_subset=None, examples=None, 
    437                  label_mask=None, **kwargs): 
    438         self.multi_flag = 1 
     112                 label_mask=None, multi_flag=False, **kwargs): 
     113        self.multi_flag = multi_flag 
    439114        self.domain = domain 
     115        self.class_var = domain.class_var 
    440116        self.best_set = best_set 
    441117        self.dirs = dirs 
     
    461137            probs.append(dist) 
    462138             
     139        if not self.multi_flag: 
     140            vals, probs = vals[0], probs[0] 
     141             
    463142        if result_type == Orange.core.GetValue: 
    464143            return vals 
    465         elif result_type == Orange.core.GetProbabilities: 
    466             return zip(vals, probs) 
     144        elif result_type == Orange.core.GetBoth: 
     145            return zip(vals, probs) if self.multi_flag else (vals, probs) 
    467146        else: 
    468147            return probs 
     
    537216                dict(self.__dict__)) 
    538217 
    539      
     218""" 
     219Utility functions 
     220----------------- 
     221""" 
     222 
     223from Orange.misc.render import contextmanager 
     224 
     225@contextmanager  
     226def member_set(obj, name, val): 
     227    """ A context manager that sets member ``name`` on ``obj`` to ``val`` 
     228    and restores the previous value on exit.  
     229    """ 
     230    old_val = getattr(obj, name, val) 
     231    setattr(obj, name, val) 
     232    yield 
     233    setattr(obj, name, old_val) 
     234     
     235     
     236def base_matrix(data, best_set, dirs, cuts): 
     237    """ Return the base matrix for the earth model. 
     238     
     239    :param data: Input data 
     240    :type data: :class:`numpy.ndarray` 
     241     
     242    :param best_set: A array of booleans indicating used terms. 
     243    :type best_set: :class:`numpy.ndarray` 
     244     
     245    :param dirs: Earth model's dirs members 
     246    :type dirs: :class:`numpy.ndarray` 
     247     
     248    :param cuts: Earth model's cuts members 
     249    :type cuts: :class:`numpy.ndarray` 
     250     
     251    """ 
     252    data = numpy.asarray(data) 
     253    best_set = numpy.asarray(best_set) 
     254    dirs = numpy.asarray(dirs) 
     255    cuts = numpy.asarray(cuts) 
     256     
     257    bx = numpy.zeros((data.shape[0], best_set.shape[0])) 
     258    bx[:, 0] = 1.0 # The intercept 
     259    for termi in range(1, best_set.shape[0]): 
     260        term_dirs = dirs[termi] 
     261        term_cuts = cuts[termi] 
     262         
     263        dir_p1 = numpy.where(term_dirs == 1)[0] 
     264        dir_m1 = numpy.where(term_dirs == -1)[0] 
     265        dir_2 = numpy.where(term_dirs == 2)[0] 
     266         
     267        x1 = data[:, dir_p1] - term_cuts[dir_p1] 
     268        x2 = term_cuts[dir_m1] - data[:, dir_m1] 
     269        x3 = data[:, dir_2] 
     270         
     271        x1 = numpy.where(x1 > 0.0, x1, 0.0) 
     272        x2 = numpy.where(x2 > 0.0, x2, 0.0) 
     273         
     274        X = numpy.hstack((x1, x2, x3)) 
     275        X = numpy.cumprod(X, axis=1) 
     276        bx[:, termi] = X[:, -1] if X.size else 0.0 
     277         
     278    return bx 
     279 
     280     
     281def gcv(rss, n, n_effective_params): 
     282    """ Return the generalized cross validation. 
     283     
     284    .. math: gcv = rss / (n * (1 - n_effective_params / n) ^ 2) 
     285     
     286    :param rss: Residual sum of squares. 
     287    :param n: Number of training examples. 
     288    :param n_effective_params: Number of effective paramaters. 
     289      
     290    """ 
     291    return  rss / (n * (1 - n_effective_params / n) ** 2) 
     292     
     293 
     294def subsets_selection_xtx_numpy(X, Y): 
     295    """ A numpy implementation of EvalSubsetsUsingXtx in the Earth package.  
     296    """ 
     297    X = numpy.asarray(X) 
     298    Y = numpy.asarray(Y) 
     299     
     300    var_count= X.shape[1] 
     301    rss_vec = numpy.zeros(var_count) 
     302    working_set = range(var_count) 
     303    subsets = numpy.zeros((var_count, var_count), dtype=int) 
     304     
     305    for subset_size in reversed(range(var_count)): 
     306        subsets[subset_size, :subset_size + 1] = working_set 
     307        X_work = X[:, working_set] 
     308        b, res, rank, s = numpy.linalg.lstsq(X_work, Y) 
     309        if res.size > 0: 
     310            rss_vec[subset_size] = numpy.sum(res) 
     311        else: 
     312            rss_vec[subset_size] = numpy.sum((Y - numpy.dot(X_work, b)) ** 2) 
     313             
     314        XtX = numpy.dot(X_work.T, X_work) 
     315        iXtX = numpy.linalg.pinv(XtX) 
     316        diag = numpy.diag(iXtX) 
     317         
     318        if subset_size == 0: 
     319            break 
     320         
     321        delta_rss = b ** 2 / diag 
     322        delete_i = numpy.argmin(delta_rss[1:]) + 1 # Keep the intercept 
     323        del working_set[delete_i] 
     324    return subsets, rss_vec 
     325 
     326 
     327""" 
     328Multi-label utility functions 
     329""" 
     330 
     331def is_label_attr(attr): 
     332    """ Is attribute a label. 
     333    """ 
     334    return attr.attributes.has_key("label") 
     335     
     336def data_label_indices(domain): 
     337    """ Return the indices of label attributes in data. 
     338    """ 
     339    return numpy.where(data_label_mask(domain))[0] 
     340 
     341def data_label_mask(domain): 
     342    """ Return an array of booleans indicating whether a variable in the 
     343    domain is a label. 
     344    """ 
     345    is_label = map(is_label_attr, domain.variables) 
     346    if domain.class_var: 
     347        is_label[-1] = True 
     348    return numpy.array(is_label, dtype=bool) 
     349 
    540350""" 
    541351ctypes interface to ForwardPass and EvalSubsetsUsingXtx. 
    542  
    543352""" 
    544353         
     
    576385     ] 
    577386     
    578 def _forward_pass(x, y, degree=1, terms=21, penalty=None, thresh=0.001, 
     387def forward_pass(x, y, degree=1, terms=21, penalty=None, thresh=0.001, 
    579388                  fast_k=21, fast_beta=1, new_var_penalty=2): 
    580     """ Do forward pass. 
     389    """ Do earth forward pass. 
    581390    """ 
    582391    import ctypes, orange 
     
    627436     ] 
    628437 
    629 def _subset_selection_xtx(X, Y): 
     438def subset_selection_xtx(X, Y): 
    630439    """ Subsets selection using EvalSubsetsUsingXtx in the Earth package. 
    631440    """ 
     
    662471     
    663472    """ 
    664     subsets, rss_vec = _subset_selection_xtx(bx, y) 
     473    subsets, rss_vec = subset_selection_xtx(bx, y) 
    665474     
    666475    cases, terms = bx.shape 
     
    680489def evimp(model, used_only=True): 
    681490    """ Return the estimated variable importance for the model. 
     491     
     492    :param model: Earth model. 
     493    :type model: `EarthClassifier` 
     494     
    682495    """ 
    683496    if model.subsets is None: 
     
    763576    """ Return a formated string representation of the model. 
    764577    """ 
    765     if model.class_var: 
    766         r_names = [model.class_var.name] 
    767         betas = [model.betas] 
    768     elif hasattr(model, "label_mask"): 
    769         mask = model.label_mask 
    770         r_vars = [v for v, m in zip(model.domain.variables, 
    771                                     model.label_mask) 
    772                   if m] 
    773         r_names = [v.name for v in r_vars] 
    774         betas = model.betas 
     578    mask = model.label_mask 
     579    r_vars = [v for v, m in zip(model.domain.variables, 
     580                                model.label_mask) 
     581              if m] 
     582    r_names = [v.name for v in r_vars] 
     583    betas = model.betas 
    775584         
    776585    resp = [] 
     
    820629             if d != 0] 
    821630    return " * ".join(knots) 
     631 
     632 
     633 
     634#class _EarthLearner(BaseEarthLearner): 
     635#    """ An earth learner.  
     636#    """ 
     637#    def __new__(cls, data=None, weightId=None, **kwargs): 
     638#        self = BaseEarthLearner.__new__(cls, **kwargs) 
     639#        if data is not None: 
     640#            self.__init__(**kwargs) 
     641#            return self.__call__(data, weightId) 
     642#        return self 
     643#     
     644#    def __init__(self, max_degree=1, max_terms=21, new_var_penalty=0.0, 
     645#                 threshold=0.001, prune=True, penalty=None, fast_k=20, 
     646#                 fast_beta=0.0, store_examples=True, **kwargs): 
     647#        """ Initialize the learner instance. 
     648#         
     649#        :param max_degree: 
     650#        """ 
     651#        self.max_degree = max_degree 
     652#        self.max_terms = max_terms 
     653#        self.new_var_penalty = new_var_penalty 
     654#        self.threshold = threshold 
     655#        self.prune = prunes 
     656#        if penaty is None: 
     657#            penalty = 2.0 if degree > 1 else 3.0 
     658#        self.penalty = penalty 
     659#        self.fast_k = fast_k 
     660#        self.fast_beta = fast_beta 
     661#        self.store_examples = store_examples 
     662#         
     663#        for key, val in kwargs.items(): 
     664#            setattr(self, key, val) 
     665#     
     666#    def __call__(self, data, weightId=None): 
     667#        if not data.domain.class_var: 
     668#            raise ValueError("No class var in the domain.") 
     669#         
     670#        with member_set(self, "prune", False): 
     671#            # We overwrite the prune argument (will do the pruning in python). 
     672#            base_clsf =  BaseEarthLearner.__call__(self, data, weightId) 
     673#         
     674#        if self.prune: 
     675#            (best_set, betas, rss, subsets, rss_per_subset, 
     676#             gcv_per_subset) = self.pruning_pass(base_clsf, data) 
     677#             
     678#            return _EarthClassifier(base_clsf, data if self.store_examples else None, 
     679#                                   best_set=best_set, dirs=base_clsf.dirs, 
     680#                                   cuts=base_clsf.cuts, 
     681#                                   betas=betas, 
     682#                                   subsets=subsets, 
     683#                                   rss_per_subset=rss_per_subset, 
     684#                                   gcv_per_subset=gcv_per_subset) 
     685#        else: 
     686#            return _EarthClassifier(base_clsf, data if self.store_examples else None) 
     687#     
     688#     
     689#    def pruning_pass(self, base_clsf, examples): 
     690#        """ Prune the terms constructed in the forward pass. 
     691#        (Pure numpy reimplementation) 
     692#        """ 
     693#        n_terms = numpy.sum(base_clsf.best_set) 
     694#        n_eff_params = n_terms + self.penalty * (n_terms - 1) / 2.0 
     695#        data, y, _ = examples.to_numpy_MA() 
     696#        data = data.filled(0.0) 
     697#        best_set = numpy.asarray(base_clsf.best_set, dtype=bool) 
     698#         
     699#        bx = base_matrix(data, base_clsf.best_set, 
     700#                         base_clsf.dirs, base_clsf.cuts, 
     701#                         ) 
     702#         
     703#        bx_used = bx[:, best_set] 
     704#        subsets, rss_per_subset = subsets_selection_xtx(bx_used, y) # TODO: Use leaps like library 
     705#        gcv_per_subset = [gcv(rss, bx.shape[0], i + self.penalty * (i - 1) / 2.0) \ 
     706#                              for i, rss in enumerate(rss_per_subset, 1)] 
     707#        gcv_per_subset = numpy.array(gcv_per_subset) 
     708#         
     709#        best_i = numpy.argmin(gcv_per_subset[1:]) + 1 # Ignore the intercept 
     710#        best_ind = subsets[best_i, :best_i + 1] 
     711#        bs_i = 0 
     712#        for i, b in enumerate(best_set): 
     713#            if b: 
     714#                best_set[i] = bs_i in best_ind 
     715#                bs_i += 1 
     716#                 
     717#        bx_subset = bx[:, best_set] 
     718#        betas, rss, rank, s = numpy.linalg.lstsq(bx_subset, y) 
     719#        return best_set, betas, rss, subsets, rss_per_subset, gcv_per_subset 
     720#     
     721#         
     722#class _EarthClassifier(Orange.core.ClassifierFD): 
     723#    def __init__(self, base_classifier=None, examples=None, best_set=None, 
     724#                 dirs=None, cuts=None, betas=None, subsets=None, 
     725#                 rss_per_subset=None, 
     726#                 gcv_per_subset=None): 
     727#        self._base_classifier = base_classifier 
     728#        self.examples = examples 
     729#        self.domain = base_classifier.domain 
     730#        self.class_var = base_classifier.class_var 
     731#         
     732#        best_set = base_classifier.best_set if best_set is None else best_set 
     733#        dirs = base_classifier.dirs if dirs is None else dirs 
     734#        cuts = base_classifier.cuts if cuts is None else cuts 
     735#        betas = base_classifier.betas if betas is None else betas 
     736#         
     737#        self.best_set = numpy.asarray(best_set, dtype=bool) 
     738#        self.dirs = numpy.array(dirs, dtype=int) 
     739#        self.cuts = numpy.array(cuts) 
     740#        self.betas = numpy.array(betas) 
     741#         
     742#        self.subsets = subsets 
     743#        self.rss_per_subset = rss_per_subset 
     744#        self.gcv_per_subset = gcv_per_subset 
     745#         
     746#    @property 
     747#    def num_terms(self): 
     748#        """ Number of terms in the model (including the intercept). 
     749#        """ 
     750#        return numpy.sum(numpy.asarray(self.best_set, dtype=int)) 
     751#     
     752#    @property 
     753#    def max_terms(self): 
     754#        """ Maximum number of terms (as specified in the learning step). 
     755#        """ 
     756#        return self.best_set.size 
     757#     
     758#    @property 
     759#    def num_preds(self): 
     760#        """ Number of predictors (variables) included in the model. 
     761#        """ 
     762#        return len(self.used_attributes(term)) 
     763#     
     764#    def __call__(self, example, what=Orange.core.GetValue): 
     765#        value = self.predict(example) 
     766#        if isinstance(self.class_var, Orange.data.variable.Continuous): 
     767#            value = self.class_var(value) 
     768#        else: 
     769#            value = self.class_var(int(round(value))) 
     770#             
     771#        dist = Orange.statistics.distribution.Distribution(self.class_var) 
     772#        dist[value] = 1.0 
     773#        if what == Orange.core.GetValue: 
     774#            return value 
     775#        elif what == Orange.core.GetProbabilities: 
     776#            return dist 
     777#        else: 
     778#            return (value, dist) 
     779#     
     780#    def base_matrix(self, examples=None): 
     781#        """ Return the base matrix (bx) of the Earth model for the table. 
     782#        If table is not supplied the base matrix of the training examples  
     783#        is returned. 
     784#         
     785#         
     786#        :param examples: Input examples for the base matrix. 
     787#        :type examples: Orange.data.Table  
     788#         
     789#        """ 
     790#        if examples is None: 
     791#            examples = self.examples 
     792#             
     793#        if examples is None: 
     794#            raise ValueError("base matrix is only available if 'store_examples=True'") 
     795#         
     796#        if isinstance(examples, Orange.data.Table): 
     797#            data, cls, w = examples.to_numpy_MA() 
     798#            data = data.filled(0.0) 
     799#        else: 
     800#            data = numpy.asarray(examples) 
     801#             
     802#        return base_matrix(data, self.best_set, self.dirs, self.cuts) 
     803#     
     804#    def _anova_order(self): 
     805#        """ Return indices that sort the terms into the 'ANOVA' format. 
     806#        """ 
     807#        terms = [([], 0)] # intercept 
     808#        for i, used in enumerate(self.best_set[1:], 1): 
     809#            attrs = sorted(self.used_attributes(i)) 
     810#            if used and attrs: 
     811#                terms.append((attrs, i)) 
     812#        terms = sotred(terms, key=lambda t:(len(t[0]), t[0])) 
     813#        return [i for _, i in terms] 
     814#     
     815#    def format_model(self, percision=3, indent=3): 
     816#        return format_model(self, percision, indent) 
     817#     
     818#    def print_model(self, percision=3, indent=3): 
     819#        print self.format_model(percision, indent) 
     820#         
     821#    def predict(self, ex): 
     822#        """ Return the predicted value (float) for example. 
     823#        """ 
     824#        x = Orange.data.Table(self.domain, [ex]) 
     825#        x, c, w = x.to_numpy_MA() 
     826#        x = x.filled(0.0)[0] 
     827#         
     828#        bx = numpy.ones(self.best_set.shape) 
     829#        betas = numpy.zeros_like(self.betas) 
     830#        betas[0] = self.betas[0] 
     831#        beta_i = 0 
     832#        for termi in range(1, len(self.best_set)): 
     833#            dirs = self.dirs[termi] 
     834#            cuts = self.cuts[termi] 
     835#            dir_p1 = numpy.where(dirs == 1)[0] 
     836#            dir_m1 = numpy.where(dirs == -1)[0] 
     837#            dir_2 = numpy.where(dirs == 2)[0] 
     838#             
     839#            x1 = x[dir_p1] - cuts[dir_p1] 
     840#            x2 = cuts[dir_m1] - x[dir_m1] 
     841#            x3 = x[dir_2] 
     842#             
     843#            x1 = numpy.maximum(x1, 0.0) 
     844#            x2 = numpy.maximum(x2, 0.0) 
     845# 
     846#            X = numpy.hstack((x1, x2, x3)) 
     847#            X = numpy.cumprod(X) 
     848#             
     849#            bx[termi] = X[-1] if X.size else 0.0 
     850#            if self.best_set[termi] != 0: 
     851#                beta_i += 1 
     852#                betas[beta_i] = self.betas[beta_i] 
     853# 
     854#        return numpy.sum(bx[self.best_set] * betas) 
     855#             
     856#    def used_attributes(self, term=None): 
     857#        """ Return a list of used attributes. If term (index) is given 
     858#        return only attributes used in that single term. 
     859#         
     860#        """ 
     861#        if term is None: 
     862#            terms = numpy.where(self.best_set)[0] 
     863#        else: 
     864#            terms = [term] 
     865#        attrs = set() 
     866#        for ti in terms: 
     867#            attri = numpy.where(self.dirs[ti] != 0.0)[0] 
     868#            attrs.update([self.domain.attributes[i] for i in attri]) 
     869#        return attrs 
     870#         
     871#    def evimp(self, used_only=True): 
     872#        """ Return the estimated variable importance. 
     873#        """ 
     874#        return evimp(self, used_only) 
     875#         
     876#    def __reduce__(self): 
     877#        return (EarthClassifier, (self._base_classifier, self.examples, 
     878#                                  self.best_set, self.dirs, self.cuts, 
     879#                                  self.betas, self.subsets, 
     880#                                  self.rss_per_subset, self.gcv_per_subset), 
     881#                {}) 
     882                                  
Note: See TracChangeset for help on using the changeset viewer.