Ignore:
Timestamp:
06/01/12 14:55:41 (23 months ago)
Author:
Ales Erjavec <ales.erjavec@…>
Branch:
default
Children:
10898:6b7e7d55a81b, 10903:04accc6a6f21
Message:

Code style fixup.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • Orange/regression/earth.py

    r10896 r10897  
    1010interactions. 
    1111 
    12 This module borrows the implementation of the technique from the `Earth R  
    13 package`_ by Stephen Milborrow.  
    14  
    15 .. _`Multivariate adaptive regression splines (MARS)`: http://en.wikipedia.org/wiki/Multivariate_adaptive_regression_splines 
     12This module borrows the implementation of the technique from the `Earth R 
     13package`_ by Stephen Milborrow. 
     14 
     15.. _`Multivariate adaptive regression splines (MARS)`: 
     16        http://en.wikipedia.org/wiki/Multivariate_adaptive_regression_splines 
     17 
    1618.. _`Earth R package`: http://cran.r-project.org/web/packages/earth/index.html 
    1719 
     
    5557from Orange.feature import Discrete, Continuous 
    5658from Orange.data import Table, Domain 
    57 from Orange.data.preprocess import Continuize as Preprocessor_continuize, \ 
    58                               Impute as Preprocessor_impute, \ 
    59                               PreprocessorList as Preprocessor_preprocessorList, \ 
    60                               DomainContinuizer 
     59from Orange.data.preprocess import DomainContinuizer 
    6160 
    6261import numpy 
     62 
    6363 
    6464def is_discrete(var): 
    6565    return isinstance(var, Discrete) 
    6666 
     67 
    6768def is_continuous(var): 
    6869    return isinstance(var, Continuous) 
     70 
    6971 
    7072def expand_discrete(var): 
     
    7274    variable for each value of ``var`` (if the number of values is grater 
    7375    then 2 else return only one indicator variable). 
    74      
     76 
    7577    """ 
    7678    if len(var.values) > 2: 
     
    9092    return new_vars 
    9193 
     94 
    9295def select_attrs(table, features, class_var=None, 
    9396                 class_vars=None, metas=None): 
     
    101104        domain.add_metas(metas) 
    102105    return Table(domain, table) 
    103      
     106 
     107 
    104108class EarthLearner(Orange.regression.base.BaseRegressionLearner): 
    105109    """Earth learner class. Supports both regression and classification 
    106     problems. For classification, class values are expanded into  
    107     continuous indicator columns (one for each value if the number of  
     110    problems. For classification, class values are expanded into 
     111    continuous indicator columns (one for each value if the number of 
    108112    values is grater then 2), and a multi response model is fit to these 
    109113    new columns. The resulting classifier the computes response 
    110114    values on new instances to select the final predicted class. 
    111       
     115 
    112116    """ 
    113117    def __new__(cls, instances=None, weight_id=None, **kwargs): 
     
    118122        else: 
    119123            return self 
    120          
    121     def __init__(self, degree=1, terms=21, penalty= None, thresh=1e-3, 
     124 
     125    def __init__(self, degree=1, terms=21, penalty=None, thresh=1e-3, 
    122126                 min_span=0, new_var_penalty=0, fast_k=20, fast_beta=1, 
    123127                 pruned_terms=None, scale_resp=True, store_instances=True, 
    124128                **kwds): 
    125129        """Initialize the learner instance. 
    126          
     130 
    127131        :param degree: Maximum degree (num. of hinge functions per term) 
    128132            of the terms in the model (default: 1). 
     
    133137                setting in earth R package. 
    134138        :type terms: int 
    135         :param penalty: Penalty for hinges in the GCV computation (used  
     139        :param penalty: Penalty for hinges in the GCV computation (used 
    136140            in the pruning pass). Default is 3.0 if ``degree`` is above 1, 
    137             and 2.0 otherwise.  
     141            and 2.0 otherwise. 
    138142        :type penalty: float 
    139143        :param thresh: Threshold for RSS decrease in the forward pass 
     
    155159            (default: ``True``). 
    156160        :type store_instances: bool 
    157           
     161 
    158162        .. todo:: min_span, prunning_method (need Leaps like functionality, 
    159             currently only eval_subsets_using_xtx is implemented).  
    160          
     163            currently only eval_subsets_using_xtx is implemented). 
     164 
    161165        """ 
    162          
     166 
    163167        super(EarthLearner, self).__init__() 
    164          
     168 
    165169        self.degree = degree 
    166170        self.terms = terms 
    167171        if penalty is None: 
    168172            penalty = 3 if degree > 1 else 2 
    169         self.penalty = penalty  
     173        self.penalty = penalty 
    170174        self.thresh = thresh 
    171175        self.min_span = min_span 
     
    177181        self.store_instances = store_instances 
    178182        self.__dict__.update(kwds) 
    179          
     183 
    180184        self.continuizer.class_treatment = DomainContinuizer.Ignore 
    181          
     185 
    182186    def __call__(self, instances, weight_id=None): 
    183          
    184187        expanded_class = None 
    185188        multitarget = False 
    186          
     189 
    187190        if instances.domain.class_var: 
    188191            instances = self.impute_table(instances) 
    189192            instances = self.continuize_table(instances) 
    190              
     193 
    191194            if is_discrete(instances.domain.class_var): 
    192195                # Expand a discrete class with indicator columns 
     
    206209            x_table = select_attrs(instances, instances.domain.attributes) 
    207210            y_table = select_attrs(instances, instances.domain.class_vars) 
    208              
     211 
    209212            # Impute and continuize only the x_table 
    210213            x_table = self.impute_table(x_table) 
     
    212215            domain = Domain(x_table.domain.attributes, 
    213216                            class_vars=instances.domain.class_vars) 
    214              
     217 
    215218            (x, ) = x_table.to_numpy_MA("A") 
    216219            (y, ) = y_table.to_numpy_MA("A") 
    217              
     220 
    218221            multitarget = True 
    219222        else: 
    220223            raise ValueError("Class variable expected.") 
    221          
     224 
    222225        if self.scale_resp and y.shape[1] == 1: 
    223226            sy = y - numpy.mean(y, axis=0) 
     
    225228        else: 
    226229            sy = y 
    227              
     230 
    228231        terms = self.terms 
    229232        if terms is None: 
    230233            # Automatic maximum number of terms 
    231234            terms = min(200, max(20, 2 * x.shape[1])) + 1 
    232              
     235 
    233236        n_terms, used, bx, dirs, cuts = forward_pass(x, sy, 
    234237            degree=self.degree, terms=terms, penalty=self.penalty, 
    235238            thresh=self.thresh, fast_k=self.fast_k, fast_beta=self.fast_beta, 
    236239            new_var_penalty=self.new_var_penalty) 
    237          
     240 
    238241        # discard unused terms from bx, dirs, cuts 
    239242        bx = bx[:, used] 
    240243        dirs = dirs[used, :] 
    241244        cuts = cuts[used, :] 
    242          
     245 
    243246        # pruning 
    244247        used, subsets, rss_per_subset, gcv_per_subset = \ 
    245248            pruning_pass(bx, y, self.penalty, 
    246249                         pruned_terms=self.pruned_terms) 
    247          
     250 
    248251        # Fit betas 
    249252        bx_used = bx[:, used] 
    250253        betas, res, rank, s = numpy.linalg.lstsq(bx_used, y) 
    251          
     254 
    252255        return EarthClassifier(instances.domain, used, dirs, cuts, betas.T, 
    253256                               subsets, rss_per_subset, gcv_per_subset, 
     
    275278        if self.multitarget: 
    276279            self.class_vars = domain.class_vars 
    277              
     280 
    278281        self.best_set = best_set 
    279282        self.dirs = dirs 
     
    287290        self.original_domain = original_domain 
    288291        self.__dict__.update(kwargs) 
    289          
     292 
    290293    def __call__(self, instance, result_type=Orange.core.GetValue): 
    291294        if self.multitarget and self.domain.class_vars: 
     
    295298        else: 
    296299            resp_vars = [self.class_var] 
    297              
     300 
    298301        vals = self.predict(instance) 
    299302        vals = [var(val) for var, val in zip(resp_vars, vals)] 
    300          
     303 
    301304        from Orange.statistics.distribution import Distribution 
    302          
     305 
    303306        if not self.multitarget and is_discrete(self.class_var): 
    304307            dist = Distribution(self.class_var) 
     
    307310            else: 
    308311                probs = soft_max(map(float, vals)) 
    309                  
     312 
    310313            for val, p in zip(self.class_var.values, probs): 
    311314                dist[val] = p 
     
    318321                dist[val] = 1.0 
    319322                probs.append(dist) 
    320              
     323 
    321324        if not self.multitarget: 
    322325            vals, probs = vals[0], probs[0] 
    323              
     326 
    324327        if result_type == Orange.core.GetValue: 
    325328            return vals 
     
    328331        else: 
    329332            return probs 
    330          
     333 
    331334    def base_matrix(self, instances=None): 
    332335        """Return the base matrix (bx) of the Earth model for the table. 
    333         If table is not supplied, the base matrix of the training instances  
     336        If table is not supplied, the base matrix of the training instances 
    334337        is returned. 
    335338        Base matrix is a len(instances) x num_terms matrix of computed values 
    336339        of terms in the model (not multiplied by beta) for each instance. 
    337          
     340 
    338341        :param instances: Input instances for the base matrix. 
    339         :type instances: :class:`Orange.data.Table`  
    340          
     342        :type instances: :class:`Orange.data.Table` 
     343 
    341344        """ 
    342345        if instances is None: 
     
    351354        The attributes can be used in Orange's domain translation 
    352355        (i.e. they define the proper ``get_value_from`` functions). 
    353          
     356 
    354357        """ 
    355358        terms = [] 
     
    357360        cuts = self.cuts[self.best_set] 
    358361 
    359         for dir, cut in zip(dirs[1:], cuts[1:]):  # Drop the intercept (first column)  
     362        for dir, cut in zip(dirs[1:], cuts[1:]):  # Drop the intercept (first column) 
    360363            hinge = [_format_knot(self, attr.name, dir1, cut1) \ 
    361364                     for (attr, dir1, cut1) in \ 
     
    370373    def predict(self, instance): 
    371374        """ Predict the response value(s) 
    372          
     375 
    373376        :param instance: Data instance 
    374377        :type instance: :class:`Orange.data.Instance` 
    375          
     378 
    376379        """ 
    377380        data = Orange.data.Table(self.domain, [instance]) 
     
    380383        vals = numpy.dot(bx_used, self.betas.T).ravel() 
    381384        return vals 
    382      
     385 
    383386    def used_attributes(self, term=None): 
    384387        """Return the used terms for term (index). If no term is given, 
    385388        return all attributes in the model. 
    386          
     389 
    387390        :param term: term index 
    388391        :type term: int 
    389          
     392 
    390393        """ 
    391394        if term is None: 
     
    393396                                      for i in range(self.best_set.size)], 
    394397                          set()) 
    395              
     398 
    396399        attrs = self.domain.attributes 
    397          
     400 
    398401        used_mask = self.dirs[term, :] != 0.0 
    399402        return [a for a, u in zip(attrs, used_mask) if u] 
    400      
     403 
    401404    def evimp(self, used_only=True): 
    402405        """ Return the estimated variable importances. 
    403          
     406 
    404407        :param used_only: if True return only used attributes 
    405          
    406         """   
     408 
     409        """ 
    407410        return evimp(self, used_only) 
    408      
     411 
    409412    def __reduce__(self): 
    410413        return (type(self), (self.domain, self.best_set, self.dirs, 
    411414                            self.cuts, self.betas), 
    412415                dict(self.__dict__)) 
    413          
     416 
    414417    def to_string(self, percision=3, indent=3): 
    415418        """ Return a string representation of the model. 
    416419        """ 
    417420        return format_model(self, percision, indent) 
    418          
     421 
    419422    def __str__(self): 
    420423        return self.to_string() 
     
    424427----------------- 
    425428""" 
    426      
     429 
     430 
    427431def base_matrix(data, best_set, dirs, cuts): 
    428432    """ Return the base matrix for the earth model. 
    429      
     433 
    430434    :param data: Input data 
    431435    :type data: :class:`numpy.ndarray` 
    432      
     436 
    433437    :param best_set: A array of booleans indicating used terms. 
    434438    :type best_set: :class:`numpy.ndarray` 
    435      
     439 
    436440    :param dirs: Earth model's dirs members 
    437441    :type dirs: :class:`numpy.ndarray` 
    438      
     442 
    439443    :param cuts: Earth model's cuts members 
    440444    :type cuts: :class:`numpy.ndarray` 
    441      
     445 
    442446    """ 
    443447    data = numpy.asarray(data) 
     
    445449    dirs = numpy.asarray(dirs) 
    446450    cuts = numpy.asarray(cuts) 
    447      
     451 
    448452    bx = numpy.zeros((data.shape[0], best_set.shape[0])) 
    449     bx[:, 0] = 1.0 # The intercept 
     453    bx[:, 0] = 1.0  # The intercept 
    450454    for termi in range(1, best_set.shape[0]): 
    451455        term_dirs = dirs[termi] 
    452456        term_cuts = cuts[termi] 
    453          
     457 
    454458        dir_p1 = numpy.where(term_dirs == 1)[0] 
    455459        dir_m1 = numpy.where(term_dirs == -1)[0] 
    456460        dir_2 = numpy.where(term_dirs == 2)[0] 
    457          
     461 
    458462        x1 = data[:, dir_p1] - term_cuts[dir_p1] 
    459463        x2 = term_cuts[dir_m1] - data[:, dir_m1] 
    460464        x3 = data[:, dir_2] 
    461          
     465 
    462466        x1 = numpy.where(x1 > 0.0, x1, 0.0) 
    463467        x2 = numpy.where(x2 > 0.0, x2, 0.0) 
    464          
     468 
    465469        X = numpy.hstack((x1, x2, x3)) 
    466470        X = numpy.cumprod(X, axis=1) 
    467471        bx[:, termi] = X[:, -1] if X.size else 0.0 
    468          
     472 
    469473    return bx 
    470474 
    471      
     475 
    472476def gcv(rss, n, n_effective_params): 
    473477    """ Return the generalized cross validation. 
    474      
     478 
    475479    .. math:: gcv = rss / (n * (1 - NumEffectiveParams / n) ^ 2) 
    476      
     480 
    477481    :param rss: Residual sum of squares. 
    478482    :param n: Number of training instances. 
    479483    :param n_effective_params: Number of effective paramaters. 
    480       
     484 
    481485    """ 
    482486    return  rss / (n * (1. - n_effective_params / n) ** 2) 
     
    529533 
    530534_c_forward_pass_.argtypes = \ 
    531     [ctypes.POINTER(ctypes.c_int),  #pnTerms: 
    532      ctypeslib.ndpointer(dtype=ctypes.c_bool, ndim=1),  #FullSet 
    533      ctypeslib.ndpointer(dtype=ctypes.c_double, ndim=2, flags="F_CONTIGUOUS"), #bx 
    534      ctypeslib.ndpointer(dtype=ctypes.c_int, ndim=2, flags="F_CONTIGUOUS"),    #Dirs 
    535      ctypeslib.ndpointer(dtype=ctypes.c_double, ndim=2, flags="F_CONTIGUOUS"), #Cuts 
    536      ctypeslib.ndpointer(dtype=ctypes.c_int, ndim=1),  #nFactorsInTerms 
    537      ctypeslib.ndpointer(dtype=ctypes.c_int, ndim=1),  #nUses 
    538      ctypeslib.ndpointer(dtype=ctypes.c_double, ndim=2, flags="F_CONTIGUOUS"), #x 
    539      ctypeslib.ndpointer(dtype=ctypes.c_double, ndim=2, flags="F_CONTIGUOUS"), #y 
     535    [ctypes.POINTER(ctypes.c_int),  # pnTerms: 
     536     ctypeslib.ndpointer(dtype=ctypes.c_bool, ndim=1),  # FullSet 
     537     ctypeslib.ndpointer(dtype=ctypes.c_double, ndim=2, flags="F_CONTIGUOUS"), # bx 
     538     ctypeslib.ndpointer(dtype=ctypes.c_int, ndim=2, flags="F_CONTIGUOUS"),    # Dirs 
     539     ctypeslib.ndpointer(dtype=ctypes.c_double, ndim=2, flags="F_CONTIGUOUS"), # Cuts 
     540     ctypeslib.ndpointer(dtype=ctypes.c_int, ndim=1),  # nFactorsInTerms 
     541     ctypeslib.ndpointer(dtype=ctypes.c_int, ndim=1),  # nUses 
     542     ctypeslib.ndpointer(dtype=ctypes.c_double, ndim=2, flags="F_CONTIGUOUS"), # x 
     543     ctypeslib.ndpointer(dtype=ctypes.c_double, ndim=2, flags="F_CONTIGUOUS"), # y 
    540544     ctypeslib.ndpointer(dtype=ctypes.c_double, ndim=1), # Weights 
    541      ctypes.c_int,  #nCases 
    542      ctypes.c_int,  #nResp 
    543      ctypes.c_int,  #nPred 
    544      ctypes.c_int,  #nMaxDegree 
    545      ctypes.c_int,  #nMaxTerms 
    546      ctypes.c_double,   #Penalty 
    547      ctypes.c_double,   #Thresh 
    548      ctypes.c_int,  #nFastK 
    549      ctypes.c_double,   #FastBeta 
    550      ctypes.c_double,   #NewVarPenalty 
    551      ctypeslib.ndpointer(dtype=ctypes.c_int, ndim=1),  #LinPreds 
    552      ctypes.c_bool, #UseBetaCache 
    553      ctypes.c_char_p    #sPredNames 
     545     ctypes.c_int,  # nCases 
     546     ctypes.c_int,  # nResp 
     547     ctypes.c_int,  # nPred 
     548     ctypes.c_int,  # nMaxDegree 
     549     ctypes.c_int,  # nMaxTerms 
     550     ctypes.c_double,   # Penalty 
     551     ctypes.c_double,   # Thresh 
     552     ctypes.c_int,  # nFastK 
     553     ctypes.c_double,   # FastBeta 
     554     ctypes.c_double,   # NewVarPenalty 
     555     ctypeslib.ndpointer(dtype=ctypes.c_int, ndim=1),  # LinPreds 
     556     ctypes.c_bool, # UseBetaCache 
     557     ctypes.c_char_p    # sPredNames 
    554558     ] 
    555      
     559 
     560 
    556561def forward_pass(x, y, degree=1, terms=21, penalty=None, thresh=0.001, 
    557562                  fast_k=21, fast_beta=1, new_var_penalty=2): 
     
    568573    n_cases = x.shape[0] 
    569574    n_preds = x.shape[1] 
    570      
     575 
    571576    n_resp = y.shape[1] if y.ndim == 2 else y.shape[0] 
    572      
     577 
    573578    # Output variables 
    574579    n_term = ctypes.c_int() 
     
    582587    lin_preds = numpy.zeros((n_preds,), dtype=ctypes.c_int, order="F") 
    583588    use_beta_cache = True 
    584      
     589 
    585590    # These tests are performed in ForwardPass, and if they fail the function 
    586591    # calls exit. So we must check it here and raise a exception to avoid a 
     
    618623                     n_factors_in_terms, n_uses, x, y, weights, n_cases, 
    619624                     n_resp, n_preds, degree, terms, penalty, thresh, 
    620                      fast_k, fast_beta, new_var_penalty, lin_preds,  
     625                     fast_k, fast_beta, new_var_penalty, lin_preds, 
    621626                     use_beta_cache, None) 
    622627    return n_term.value, full_set, bx, dirs, cuts 
     
    626631 
    627632_c_eval_subsets_xtx.argtypes = \ 
    628     [ctypeslib.ndpointer(dtype=ctypes.c_bool, ndim=2, flags="F_CONTIGUOUS"),   #PruneTerms 
    629      ctypeslib.ndpointer(dtype=ctypes.c_double, ndim=1),   #RssVec 
    630      ctypes.c_int,  #nCases 
    631      ctypes.c_int,  #nResp 
    632      ctypes.c_int,  #nMaxTerms 
    633      ctypeslib.ndpointer(dtype=ctypes.c_double, ndim=2, flags="F_CONTIGUOUS"),   #bx 
    634      ctypeslib.ndpointer(dtype=ctypes.c_double, ndim=2, flags="F_CONTIGUOUS"),   #y 
    635      ctypeslib.ndpointer(dtype=ctypes.c_double, ndim=1)  #WeightsArg 
     633    [ctypeslib.ndpointer(dtype=ctypes.c_bool, ndim=2, flags="F_CONTIGUOUS"),   # PruneTerms 
     634     ctypeslib.ndpointer(dtype=ctypes.c_double, ndim=1),   # RssVec 
     635     ctypes.c_int,  # nCases 
     636     ctypes.c_int,  # nResp 
     637     ctypes.c_int,  # nMaxTerms 
     638     ctypeslib.ndpointer(dtype=ctypes.c_double, ndim=2, flags="F_CONTIGUOUS"),  # bx 
     639     ctypeslib.ndpointer(dtype=ctypes.c_double, ndim=2, flags="F_CONTIGUOUS"),  # y 
     640     ctypeslib.ndpointer(dtype=ctypes.c_double, ndim=1)  # WeightsArg 
    636641     ] 
    637      
     642 
    638643_c_eval_subsets_xtx.restype = ctypes.c_int 
     644 
    639645 
    640646def subset_selection_xtx(X, Y): 
     
    645651    if Y.ndim == 1: 
    646652        Y = Y.reshape((-1, 1), order="F") 
    647          
     653 
    648654    if X.shape[0] != Y.shape[0]: 
    649655        raise ValueError("First dimensions of bx and y must be the same") 
    650          
     656 
    651657    var_count = X.shape[1] 
    652658    resp_count = Y.shape[1] 
     
    656662    rss_vec = numpy.zeros((var_count,), dtype=ctypes.c_double, order="F") 
    657663    weights = numpy.ones((cases,), dtype=ctypes.c_double, order="F") 
    658      
     664 
    659665    rval = _c_eval_subsets_xtx(subsets, rss_vec, cases, resp_count, var_count, 
    660666                        X, Y, weights) 
    661667    if rval != 0: 
    662668        raise numpy.linalg.LinAlgError("Lin. dep. terms in X") 
    663      
     669 
    664670    subsets_ind = numpy.zeros((var_count, var_count), dtype=int) 
    665671    for i, used in enumerate(subsets.T): 
    666672        subsets_ind[i, :i + 1] = numpy.where(used)[0] 
    667          
     673 
    668674    return subsets_ind, rss_vec 
    669      
     675 
     676 
    670677def subset_selection_xtx_numpy(X, Y): 
    671678    """ A numpy implementation of EvalSubsetsUsingXtx in the Earth package. 
    672679    Using numpy.linalg.lstsq 
    673      
     680 
    674681    """ 
    675682    X = numpy.asarray(X) 
    676683    Y = numpy.asarray(Y) 
    677      
     684 
    678685    var_count = X.shape[1] 
    679686    rss_vec = numpy.zeros(var_count) 
    680687    working_set = range(var_count) 
    681688    subsets = numpy.zeros((var_count, var_count), dtype=int) 
    682      
     689 
    683690    for subset_size in reversed(range(var_count)): 
    684691        subsets[subset_size, :subset_size + 1] = working_set 
     
    689696        else: 
    690697            rss_vec[subset_size] = numpy.sum((Y - numpy.dot(X_work, b)) ** 2) 
    691              
     698 
    692699        XtX = numpy.dot(X_work.T, X_work) 
    693700        iXtX = numpy.linalg.pinv(XtX) 
    694701        diag = numpy.diag(iXtX).reshape((-1, 1)) 
    695          
     702 
    696703        if subset_size == 0: 
    697704            break 
    698          
     705 
    699706        delta_rss = b ** 2 / diag 
    700707        delta_rss = numpy.sum(delta_rss, axis=1) 
    701         delete_i = numpy.argmin(delta_rss[1:]) + 1 # Keep the intercept 
     708        delete_i = numpy.argmin(delta_rss[1:]) + 1  # Keep the intercept 
    702709        del working_set[delete_i] 
    703710    return subsets, rss_vec 
    704      
     711 
     712 
    705713def subset_selection_xtx2(X, Y): 
    706714    """ Another implementation (this uses qr decomp). 
     
    714722    rss_vec = numpy.zeros((col_count,)) 
    715723    QR, k, _, jpvt = linalg.qr_decomp(X) 
    716      
     724 
    717725    if k < col_count: 
    718726        # remove jpvt[k:] from the work set. Will have zero  
     
    722730            rss_vec[len(working_set)] = float("inf") 
    723731        col_count = len(working_set) 
    724          
     732 
    725733    for subset_size in reversed(range(col_count)): 
    726734        subsets[subset_size, :subset_size + 1] = working_set 
     
    731739        iXtX = numpy.linalg.pinv(XtX) 
    732740        diag = numpy.diag(iXtX) 
    733          
     741 
    734742        if subset_size == 0: 
    735743            break 
    736          
     744 
    737745        delta_rss = b ** 2 / diag 
    738         delete_i = numpy.argmin(delta_rss[1:]) + 1 # Keep the intercept 
     746        delete_i = numpy.argmin(delta_rss[1:]) + 1  # Keep the intercept 
    739747        del working_set[delete_i] 
    740748    return subsets, rss_vec 
    741      
     749 
     750 
    742751def pruning_pass(bx, y, penalty, pruned_terms=-1): 
    743752    """ Do pruning pass 
    744      
     753 
    745754    .. todo:: pruned_terms, Leaps 
    746      
     755 
    747756    """ 
    748757    try: 
     
    750759    except numpy.linalg.LinAlgError: 
    751760        subsets, rss_vec = subset_selection_xtx_numpy(bx, y) 
    752      
     761 
    753762    cases, terms = bx.shape 
    754763    n_effective_params = numpy.arange(terms) + 1.0 
    755764    n_effective_params += penalty * (n_effective_params - 1.0) / 2.0 
    756      
     765 
    757766    gcv_vec = gcv(rss_vec, cases, n_effective_params) 
    758      
     767 
    759768    min_i = numpy.argmin(gcv_vec) 
    760769    used = numpy.zeros((terms), dtype=bool) 
    761      
     770 
    762771    used[subsets[min_i, :min_i + 1]] = True 
    763      
     772 
    764773    return used, subsets, rss_vec, gcv_vec 
    765774 
     
    767776Printing functions. 
    768777""" 
    769      
     778 
     779 
    770780def format_model(model, percision=3, indent=3): 
    771781    """ Return a formated string representation of the earth model. 
     
    777787    else: 
    778788        r_vars = [model.class_var] 
    779          
     789 
    780790    r_names = [v.name for v in r_vars] 
    781791    betas = model.betas 
    782          
     792 
    783793    resp = [] 
    784794    for name, betas in zip(r_names, betas): 
     
    786796                                     percision, indent)) 
    787797    return "\n\n".join(resp) 
     798 
    788799 
    789800def _format_response(model, resp_name, betas, percision=3, indent=3): 
     
    813824    terms = sorted(terms, key=lambda t: (len(t[0]), t[0])) 
    814825    return "\n".join([header] + [indent + t for _, t in terms]) 
    815          
     826 
     827 
    816828def _format_knot(model, name, dir, cut, percision=3): 
    817829    fmt = "%%.%if" % percision 
     
    824836    return txt 
    825837 
    826      
     838 
    827839"""\ 
    828840Variable importance estimation 
     
    831843 
    832844from collections import defaultdict 
     845 
    833846 
    834847def collect_source(vars): 
     
    837850    back to the variables in ``vars`` (assumes the default preprocessor in 
    838851    EarthLearner). 
    839      
     852 
    840853    """ 
    841854    source = defaultdict(list) 
     
    852865    return dict(source) 
    853866 
     867 
    854868def map_to_source_var(var, sources): 
    855     """  
     869    """ 
    856870    """ 
    857871    if var in sources: 
     
    866880    else: 
    867881        return None 
    868      
     882 
     883 
    869884def evimp(model, used_only=True): 
    870885    """ Return the estimated variable importance for the model. 
    871      
     886 
    872887    :param model: Earth model. 
    873888    :type model: `EarthClassifier` 
    874      
     889 
    875890    """ 
    876891    if model.subsets is None: 
    877892        raise ValueError("No subsets. Use the learner with 'prune=True'.") 
    878      
     893 
    879894    subsets = model.subsets 
    880895    n_subsets = numpy.sum(model.best_set) 
    881      
     896 
    882897    rss = -numpy.diff(model.rss_per_subset) 
    883898    gcv = -numpy.diff(model.gcv_per_subset) 
    884899    attributes = list(model.domain.variables) 
    885      
     900 
    886901    attr2ind = dict(zip(attributes, range(len(attributes)))) 
    887902    importances = numpy.zeros((len(attributes), 4)) 
    888903    importances[:, 0] = range(len(attributes)) 
    889      
     904 
    890905    for i in range(1, n_subsets): 
    891906        term_subset = subsets[i, :i + 1] 
     
    898913    imp_min = numpy.min(importances[:, [2, 3]], axis=0) 
    899914    imp_max = numpy.max(importances[:, [2, 3]], axis=0) 
    900      
     915 
    901916    #Normalize importances. 
    902917    importances[:, [2, 3]] = 100.0 * (importances[:, [2, 3]] \ 
    903918                            - [imp_min]) / ([imp_max - imp_min]) 
    904      
     919 
    905920    importances = list(importances) 
    906921    # Sort by n_subsets and gcv. 
     
    908923                         reverse=True) 
    909924    importances = numpy.array(importances) 
    910      
     925 
    911926    if used_only: 
    912         importances = importances[importances[:,1] > 0.0] 
    913      
     927        importances = importances[importances[:, 1] > 0.0] 
     928 
    914929    res = [(attributes[int(row[0])], tuple(row[1:])) for row in importances] 
    915930    return res 
     
    919934    """ Plot the variable importances as returned from 
    920935    :obj:`EarthClassifier.evimp` call. 
    921      
     936 
    922937    :: 
    923938 
     
    928943 
    929944    .. image:: files/earth-evimp.png 
    930       
     945 
    931946    The left axis is the nsubsets measure and on the right are the normalized 
    932947    RSS and GCV. 
    933      
     948 
    934949    """ 
    935950    from Orange.ensemble.bagging import BaggedClassifier 
     
    938953    elif isinstance(evimp, BaggedClassifier): 
    939954        evimp = bagged_evimp(evimp) 
    940          
     955 
    941956    import pylab 
    942957    fig = pylab.figure() 
     
    948963    l1 = axes1.plot(X, imp[:, 0], "b-",) 
    949964    axes2 = axes1.twinx() 
    950      
     965 
    951966    l2 = axes2.plot(X, imp[:, 1], "g-",) 
    952967    l3 = axes2.plot(X, imp[:, 2], "r-",) 
    953      
     968 
    954969    x_axis = axes1.xaxis 
    955970    x_axis.set_ticks(X) 
    956971    x_axis.set_ticklabels([a.name for a in attrs], rotation=90) 
    957      
     972 
    958973    axes1.yaxis.set_label_text("nsubsets") 
    959974    axes2.yaxis.set_label_text("normalized gcv or rss") 
     
    963978    fig.show() 
    964979 
    965      
     980 
    966981def bagged_evimp(classifier, used_only=True): 
    967982    """ Extract combined (average) evimp from an instance of BaggedClassifier 
    968      
     983 
    969984    Example:: 
    970985 
     
    977992        if not isinstance(object, class_): 
    978993            raise TypeError("Instance of %r expected." % (class_)) 
    979      
     994 
    980995    from Orange.ensemble.bagging import BaggedClassifier 
    981      
     996 
    982997    assert_type(classifier, BaggedClassifier) 
    983998    bagged_imp = defaultdict(list) 
     
    9871002        imp = evimp(c, used_only=used_only) 
    9881003        for attr, score in imp: 
    989             bagged_imp[attr.name].append(score) # map by name 
     1004            bagged_imp[attr.name].append(score)  # map by name 
    9901005            attrs_by_name[attr.name].append(attr) 
    991              
     1006 
    9921007    for attr, scores in bagged_imp.items(): 
    9931008        scores = numpy.average(scores, axis=0) 
    9941009        bagged_imp[attr] = tuple(scores) 
    995      
    996      
     1010 
    9971011    bagged_imp = sorted(bagged_imp.items(), 
    9981012                        key=lambda t: (t[1][0], t[1][1]), 
    9991013                        reverse=True) 
    1000      
     1014 
    10011015    bagged_imp = [(attrs_by_name[name][0], scores) for name, scores in bagged_imp] 
    1002      
     1016 
    10031017    if used_only: 
    10041018        bagged_imp = [(a, r) for a, r in bagged_imp if r[0] > 0] 
     
    10111025""" 
    10121026from Orange.feature import scoring 
    1013              
     1027 
     1028 
    10141029class ScoreEarthImportance(scoring.Score): 
    10151030    """ A subclass of :class:`Orange.feature.scoring.Score` that. 
    10161031    scores features based on their importance in the Earth 
    10171032    model using ``bagged_evimp``. 
    1018      
    1019     """ 
    1020     # Return types   
     1033 
     1034    """ 
     1035    # Return types 
    10211036    NSUBSETS = 0 
    10221037    RSS = 1 
    10231038    GCV = 2 
    1024      
     1039 
    10251040    handles_discrete = True 
    10261041    handles_continuous = True 
    10271042    computes_thresholds = False 
    10281043    needs = scoring.Score.Generator 
    1029      
     1044 
    10301045    def __new__(cls, attr=None, data=None, weight_id=None, **kwargs): 
    10311046        self = scoring.Score.__new__(cls) 
     
    10381053        else: 
    10391054            raise ValueError("Both 'attr' and 'data' arguments expected.") 
    1040          
    1041     def __init__(self, t=10, degree=2, terms=10, score_what="nsubsets", cached=True): 
     1055 
     1056    def __init__(self, t=10, degree=2, terms=10, score_what="nsubsets", 
     1057                 cached=True): 
    10421058        """ 
    10431059        :param t: Number of earth models to train on the data 
    10441060            (using BaggedLearner). 
    1045              
     1061 
    10461062        :param score_what: What to return as a score. 
    10471063            Can be one of: "nsubsets", "rss", "gcv" or class constants 
    10481064            NSUBSETS, RSS, GCV. 
    1049              
     1065 
    10501066        """ 
    10511067        self.t = t 
     
    10531069        self.terms = terms 
    10541070        if isinstance(score_what, basestring): 
    1055             score_what = {"nsubsets":self.NSUBSETS, "rss":self.RSS, 
    1056                           "gcv":self.GCV}.get(score_what, None) 
    1057                            
     1071            score_what = {"nsubsets": self.NSUBSETS, "rss": self.RSS, 
     1072                          "gcv": self.GCV}.get(score_what, None) 
     1073 
    10581074        if score_what not in range(3): 
    10591075            raise ValueError("Invalid  'score_what' parameter.") 
     
    10631079        self._cache_ref = None 
    10641080        self._cached_evimp = None 
    1065          
     1081 
    10661082    def __call__(self, attr, data, weight_id=None): 
    10671083        ref = self._cache_ref 
     
    10751091            self._cache_ref = data 
    10761092            self._cached_evimp = evimp 
    1077          
     1093 
    10781094        evimp = dict(evimp) 
    10791095        score = evimp.get(attr, None) 
    1080          
     1096 
    10811097        if score is None: 
    10821098            source = collect_source(evimp.keys()) 
     
    10891105        else: 
    10901106            return score[self.score_what] 
    1091      
     1107 
     1108 
    10921109class ScoreRSS(scoring.Score): 
    1093      
     1110 
    10941111    handles_discrete = False 
    10951112    handles_continuous = True 
    10961113    computes_thresholds = False 
    1097      
     1114 
    10981115    def __new__(cls, attr=None, data=None, weight_id=None, **kwargs): 
    10991116        self = scoring.Score.__new__(cls) 
     
    11061123        else: 
    11071124            raise ValueError("Both 'attr' and 'data' arguments expected.") 
    1108          
     1125 
    11091126    def __init__(self): 
    11101127        self._cache_data = None 
    11111128        self._cache_rss = None 
    1112          
     1129 
    11131130    def __call__(self, attr, data, weight_id=None): 
    11141131        ref = self._cache_data 
     
    11261143                subset = subsets[s_size, :s_size + 1] 
    11271144                rss[subset] += rss_diff[s_size - 1] 
    1128             rss = rss[1:] #Drop the intercept 
     1145            rss = rss[1:]  # Drop the intercept 
    11291146            self._cache_data = data 
    11301147            self._cache_rss = rss 
    1131              
     1148 
    11321149        index = list(data.domain.attributes).index(attr) 
    11331150        return rss[index] 
    1134      
    1135  
    1136 #from Orange.core import EarthLearner as BaseEarthLearner, \ 
    1137 #                        EarthClassifier as BaseEarthClassifier 
    1138 #from Orange.misc import member_set 
    1139 #  
    1140 #class _EarthLearner(BaseEarthLearner): 
    1141 #    """ An earth learner.  
    1142 #    """ 
    1143 #    def __new__(cls, data=None, weightId=None, **kwargs): 
    1144 #        self = BaseEarthLearner.__new__(cls, **kwargs) 
    1145 #        if data is not None: 
    1146 #            self.__init__(**kwargs) 
    1147 #            return self.__call__(data, weightId) 
    1148 #        return self 
    1149 #     
    1150 #    def __init__(self, max_degree=1, max_terms=21, new_var_penalty=0.0, 
    1151 #                 threshold=0.001, prune=True, penalty=None, fast_k=20, 
    1152 #                 fast_beta=0.0, store_examples=True, **kwargs): 
    1153 #        """ Initialize the learner instance. 
    1154 #         
    1155 #        :param max_degree: 
    1156 #        """ 
    1157 #        self.max_degree = max_degree 
    1158 #        self.max_terms = max_terms 
    1159 #        self.new_var_penalty = new_var_penalty 
    1160 #        self.threshold = threshold 
    1161 #        self.prune = prunes 
    1162 #        if penaty is None: 
    1163 #            penalty = 2.0 if degree > 1 else 3.0 
    1164 #        self.penalty = penalty 
    1165 #        self.fast_k = fast_k 
    1166 #        self.fast_beta = fast_beta 
    1167 #        self.store_examples = store_examples 
    1168 #         
    1169 #        for key, val in kwargs.items(): 
    1170 #            setattr(self, key, val) 
    1171 #     
    1172 #    def __call__(self, data, weightId=None): 
    1173 #        if not data.domain.class_var: 
    1174 #            raise ValueError("No class var in the domain.") 
    1175 #         
    1176 #        with member_set(self, "prune", False): 
    1177 #            # We overwrite the prune argument (will do the pruning in python). 
    1178 #            base_clsf =  BaseEarthLearner.__call__(self, data, weightId) 
    1179 #         
    1180 #        if self.prune: 
    1181 #            (best_set, betas, rss, subsets, rss_per_subset, 
    1182 #             gcv_per_subset) = self.pruning_pass(base_clsf, data) 
    1183 #             
    1184 #            return _EarthClassifier(base_clsf, data if self.store_examples else None, 
    1185 #                                   best_set=best_set, dirs=base_clsf.dirs, 
    1186 #                                   cuts=base_clsf.cuts, 
    1187 #                                   betas=betas, 
    1188 #                                   subsets=subsets, 
    1189 #                                   rss_per_subset=rss_per_subset, 
    1190 #                                   gcv_per_subset=gcv_per_subset) 
    1191 #        else: 
    1192 #            return _EarthClassifier(base_clsf, data if self.store_examples else None) 
    1193 #     
    1194 #     
    1195 #    def pruning_pass(self, base_clsf, examples): 
    1196 #        """ Prune the terms constructed in the forward pass. 
    1197 #        (Pure numpy reimplementation) 
    1198 #        """ 
    1199 #        n_terms = numpy.sum(base_clsf.best_set) 
    1200 #        n_eff_params = n_terms + self.penalty * (n_terms - 1) / 2.0 
    1201 #        data, y, _ = examples.to_numpy_MA() 
    1202 #        data = data.filled(0.0) 
    1203 #        best_set = numpy.asarray(base_clsf.best_set, dtype=bool) 
    1204 #         
    1205 #        bx = base_matrix(data, base_clsf.best_set, 
    1206 #                         base_clsf.dirs, base_clsf.cuts, 
    1207 #                         ) 
    1208 #         
    1209 #        bx_used = bx[:, best_set] 
    1210 #        subsets, rss_per_subset = subsets_selection_xtx(bx_used, y) # TODO: Use leaps like library 
    1211 #        gcv_per_subset = [gcv(rss, bx.shape[0], i + self.penalty * (i - 1) / 2.0) \ 
    1212 #                              for i, rss in enumerate(rss_per_subset, 1)] 
    1213 #        gcv_per_subset = numpy.array(gcv_per_subset) 
    1214 #         
    1215 #        best_i = numpy.argmin(gcv_per_subset[1:]) + 1 # Ignore the intercept 
    1216 #        best_ind = subsets[best_i, :best_i + 1] 
    1217 #        bs_i = 0 
    1218 #        for i, b in enumerate(best_set): 
    1219 #            if b: 
    1220 #                best_set[i] = bs_i in best_ind 
    1221 #                bs_i += 1 
    1222 #                 
    1223 #        bx_subset = bx[:, best_set] 
    1224 #        betas, rss, rank, s = numpy.linalg.lstsq(bx_subset, y) 
    1225 #        return best_set, betas, rss, subsets, rss_per_subset, gcv_per_subset 
    1226 #     
    1227 #         
    1228 #class _EarthClassifier(Orange.core.ClassifierFD): 
    1229 #    def __init__(self, base_classifier=None, examples=None, best_set=None, 
    1230 #                 dirs=None, cuts=None, betas=None, subsets=None, 
    1231 #                 rss_per_subset=None, 
    1232 #                 gcv_per_subset=None): 
    1233 #        self._base_classifier = base_classifier 
    1234 #        self.examples = examples 
    1235 #        self.domain = base_classifier.domain 
    1236 #        self.class_var = base_classifier.class_var 
    1237 #         
    1238 #        best_set = base_classifier.best_set if best_set is None else best_set 
    1239 #        dirs = base_classifier.dirs if dirs is None else dirs 
    1240 #        cuts = base_classifier.cuts if cuts is None else cuts 
    1241 #        betas = base_classifier.betas if betas is None else betas 
    1242 #         
    1243 #        self.best_set = numpy.asarray(best_set, dtype=bool) 
    1244 #        self.dirs = numpy.array(dirs, dtype=int) 
    1245 #        self.cuts = numpy.array(cuts) 
    1246 #        self.betas = numpy.array(betas) 
    1247 #         
    1248 #        self.subsets = subsets 
    1249 #        self.rss_per_subset = rss_per_subset 
    1250 #        self.gcv_per_subset = gcv_per_subset 
    1251 #         
    1252 #    @property 
    1253 #    def num_terms(self): 
    1254 #        """ Number of terms in the model (including the intercept). 
    1255 #        """ 
    1256 #        return numpy.sum(numpy.asarray(self.best_set, dtype=int)) 
    1257 #     
    1258 #    @property 
    1259 #    def max_terms(self): 
    1260 #        """ Maximum number of terms (as specified in the learning step). 
    1261 #        """ 
    1262 #        return self.best_set.size 
    1263 #     
    1264 #    @property 
    1265 #    def num_preds(self): 
    1266 #        """ Number of predictors (variables) included in the model. 
    1267 #        """ 
    1268 #        return len(self.used_attributes(term)) 
    1269 #     
    1270 #    def __call__(self, example, what=Orange.core.GetValue): 
    1271 #        value = self.predict(example) 
    1272 #        if isinstance(self.class_var, Orange.feature.Continuous): 
    1273 #            value = self.class_var(value) 
    1274 #        else: 
    1275 #            value = self.class_var(int(round(value))) 
    1276 #             
    1277 #        dist = Orange.statistics.distribution.Distribution(self.class_var) 
    1278 #        dist[value] = 1.0 
    1279 #        if what == Orange.core.GetValue: 
    1280 #            return value 
    1281 #        elif what == Orange.core.GetProbabilities: 
    1282 #            return dist 
    1283 #        else: 
    1284 #            return (value, dist) 
    1285 #     
    1286 #    def base_matrix(self, examples=None): 
    1287 #        """ Return the base matrix (bx) of the Earth model for the table. 
    1288 #        If table is not supplied the base matrix of the training examples  
    1289 #        is returned. 
    1290 #         
    1291 #         
    1292 #        :param examples: Input examples for the base matrix. 
    1293 #        :type examples: Orange.data.Table  
    1294 #         
    1295 #        """ 
    1296 #        if examples is None: 
    1297 #            examples = self.examples 
    1298 #             
    1299 #        if examples is None: 
    1300 #            raise ValueError("base matrix is only available if 'store_examples=True'") 
    1301 #         
    1302 #        if isinstance(examples, Orange.data.Table): 
    1303 #            data, cls, w = examples.to_numpy_MA() 
    1304 #            data = data.filled(0.0) 
    1305 #        else: 
    1306 #            data = numpy.asarray(examples) 
    1307 #             
    1308 #        return base_matrix(data, self.best_set, self.dirs, self.cuts) 
    1309 #     
    1310 #    def _anova_order(self): 
    1311 #        """ Return indices that sort the terms into the 'ANOVA' format. 
    1312 #        """ 
    1313 #        terms = [([], 0)] # intercept 
    1314 #        for i, used in enumerate(self.best_set[1:], 1): 
    1315 #            attrs = sorted(self.used_attributes(i)) 
    1316 #            if used and attrs: 
    1317 #                terms.append((attrs, i)) 
    1318 #        terms = sotred(terms, key=lambda t:(len(t[0]), t[0])) 
    1319 #        return [i for _, i in terms] 
    1320 #     
    1321 #    def format_model(self, percision=3, indent=3): 
    1322 #        return format_model(self, percision, indent) 
    1323 #     
    1324 #    def print_model(self, percision=3, indent=3): 
    1325 #        print self.format_model(percision, indent) 
    1326 #         
    1327 #    def predict(self, ex): 
    1328 #        """ Return the predicted value (float) for example. 
    1329 #        """ 
    1330 #        x = Orange.data.Table(self.domain, [ex]) 
    1331 #        x, c, w = x.to_numpy_MA() 
    1332 #        x = x.filled(0.0)[0] 
    1333 #         
    1334 #        bx = numpy.ones(self.best_set.shape) 
    1335 #        betas = numpy.zeros_like(self.betas) 
    1336 #        betas[0] = self.betas[0] 
    1337 #        beta_i = 0 
    1338 #        for termi in range(1, len(self.best_set)): 
    1339 #            dirs = self.dirs[termi] 
    1340 #            cuts = self.cuts[termi] 
    1341 #            dir_p1 = numpy.where(dirs == 1)[0] 
    1342 #            dir_m1 = numpy.where(dirs == -1)[0] 
    1343 #            dir_2 = numpy.where(dirs == 2)[0] 
    1344 #             
    1345 #            x1 = x[dir_p1] - cuts[dir_p1] 
    1346 #            x2 = cuts[dir_m1] - x[dir_m1] 
    1347 #            x3 = x[dir_2] 
    1348 #             
    1349 #            x1 = numpy.maximum(x1, 0.0) 
    1350 #            x2 = numpy.maximum(x2, 0.0) 
    1351 # 
    1352 #            X = numpy.hstack((x1, x2, x3)) 
    1353 #            X = numpy.cumprod(X) 
    1354 #             
    1355 #            bx[termi] = X[-1] if X.size else 0.0 
    1356 #            if self.best_set[termi] != 0: 
    1357 #                beta_i += 1 
    1358 #                betas[beta_i] = self.betas[beta_i] 
    1359 # 
    1360 #        return numpy.sum(bx[self.best_set] * betas) 
    1361 #             
    1362 #    def used_attributes(self, term=None): 
    1363 #        """ Return a list of used attributes. If term (index) is given 
    1364 #        return only attributes used in that single term. 
    1365 #         
    1366 #        """ 
    1367 #        if term is None: 
    1368 #            terms = numpy.where(self.best_set)[0] 
    1369 #        else: 
    1370 #            terms = [term] 
    1371 #        attrs = set() 
    1372 #        for ti in terms: 
    1373 #            attri = numpy.where(self.dirs[ti] != 0.0)[0] 
    1374 #            attrs.update([self.domain.attributes[i] for i in attri]) 
    1375 #        return attrs 
    1376 #         
    1377 #    def evimp(self, used_only=True): 
    1378 #        """ Return the estimated variable importance. 
    1379 #        """ 
    1380 #        return evimp(self, used_only) 
    1381 #         
    1382 #    def __reduce__(self): 
    1383 #        return (EarthClassifier, (self._base_classifier, self.examples, 
    1384 #                                  self.best_set, self.dirs, self.cuts, 
    1385 #                                  self.betas, self.subsets, 
    1386 #                                  self.rss_per_subset, self.gcv_per_subset), 
    1387 #                {}) 
    1388                                   
Note: See TracChangeset for help on using the changeset viewer.