Changeset 8920:00df413757bc in orange


Ignore:
Timestamp:
09/07/11 11:55:03 (3 years ago)
Author:
ales_erjavec <ales.erjavec@…>
Branch:
default
Convert:
8b58013f6c2881d912ff0170b967320a35f023a6
Message:

Changed EvalSubsetsUsingXtx in earth package to return an error code if there are lin. dep. terms in bx.
Fixes #932 (not completely - there are still alot of places where exit() can get called but are nested more deeply (most importantly malloc errors) and would require to much work, but this one is the only one that fails under normal circumstances).

Files:
3 edited

Legend:

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

    r8884 r8920  
    426426    """ 
    427427    return  rss / (n * (1 - n_effective_params / n) ** 2) 
    428      
    429  
    430 def subsets_selection_xtx_numpy(X, Y): 
    431     """ A numpy implementation of EvalSubsetsUsingXtx in the Earth package.  
    432     """ 
    433     X = numpy.asarray(X) 
    434     Y = numpy.asarray(Y) 
    435      
    436     var_count= X.shape[1] 
    437     rss_vec = numpy.zeros(var_count) 
    438     working_set = range(var_count) 
    439     subsets = numpy.zeros((var_count, var_count), dtype=int) 
    440      
    441     for subset_size in reversed(range(var_count)): 
    442         subsets[subset_size, :subset_size + 1] = working_set 
    443         X_work = X[:, working_set] 
    444         b, res, rank, s = numpy.linalg.lstsq(X_work, Y) 
    445         if res.size > 0: 
    446             rss_vec[subset_size] = numpy.sum(res) 
    447         else: 
    448             rss_vec[subset_size] = numpy.sum((Y - numpy.dot(X_work, b)) ** 2) 
    449              
    450         XtX = numpy.dot(X_work.T, X_work) 
    451         iXtX = numpy.linalg.pinv(XtX) 
    452         diag = numpy.diag(iXtX) 
    453          
    454         if subset_size == 0: 
    455             break 
    456          
    457         delta_rss = b ** 2 / diag 
    458         delete_i = numpy.argmin(delta_rss[1:]) + 1 # Keep the intercept 
    459         del working_set[delete_i] 
    460     return subsets, rss_vec 
    461  
    462428 
    463429""" 
     
    603569     ctypeslib.ndpointer(dtype=ctypes.c_double, ndim=1)  #WeightsArg 
    604570     ] 
     571     
     572_c_eval_subsets_xtx.restype = ctypes.c_int 
    605573 
    606574def subset_selection_xtx(X, Y): 
     
    623591    weights = numpy.ones((cases,), dtype=ctypes.c_double, order="F") 
    624592     
    625     _c_eval_subsets_xtx(subsets, rss_vec, cases, resp_count, var_count, 
     593    rval = _c_eval_subsets_xtx(subsets, rss_vec, cases, resp_count, var_count, 
    626594                        X, Y, weights) 
     595    if rval != 0: 
     596        raise numpy.linalg.LinAlgError("Lin. dep. terms in X") 
    627597     
    628598    subsets_ind = numpy.zeros((var_count, var_count), dtype=int) 
     
    632602    return subsets_ind, rss_vec 
    633603     
     604def subset_selection_xtx_numpy(X, Y): 
     605    """ A numpy implementation of EvalSubsetsUsingXtx in the Earth package. 
     606    Using numpy.linalg.lstsq 
     607     
     608    """ 
     609    X = numpy.asarray(X) 
     610    Y = numpy.asarray(Y) 
     611     
     612    var_count= X.shape[1] 
     613    rss_vec = numpy.zeros(var_count) 
     614    working_set = range(var_count) 
     615    subsets = numpy.zeros((var_count, var_count), dtype=int) 
     616     
     617    for subset_size in reversed(range(var_count)): 
     618        subsets[subset_size, :subset_size + 1] = working_set 
     619        X_work = X[:, working_set] 
     620        b, res, rank, s = numpy.linalg.lstsq(X_work, Y) 
     621        if res.size > 0: 
     622            rss_vec[subset_size] = numpy.sum(res) 
     623        else: 
     624            rss_vec[subset_size] = numpy.sum((Y - numpy.dot(X_work, b)) ** 2) 
     625             
     626        XtX = numpy.dot(X_work.T, X_work) 
     627        iXtX = numpy.linalg.pinv(XtX) 
     628        diag = numpy.diag(iXtX) 
     629         
     630        if subset_size == 0: 
     631            break 
     632         
     633        delta_rss = b ** 2 / diag 
     634        delete_i = numpy.argmin(delta_rss[1:]) + 1 # Keep the intercept 
     635        del working_set[delete_i] 
     636    return subsets, rss_vec 
     637     
    634638     
    635639def pruning_pass(bx, y, penalty, pruned_terms=-1): 
     
    639643     
    640644    """ 
    641     subsets, rss_vec = subset_selection_xtx(bx, y) 
     645    try: 
     646        subsets, rss_vec = subset_selection_xtx(bx, y) 
     647    except numpy.linalg.LinAlgError: 
     648        subsets, rss_vec = subset_selection_xtx_numpy(bx, y) 
    642649     
    643650    cases, terms = bx.shape 
     
    986993            return score[self.score_what] 
    987994     
    988 #class ScoreRSS(scoring.Score): 
    989 #    __new__ = _orange__new__(scoring.Score) 
    990 #    def __init__(self): 
    991 #        self._cache_data = None 
    992 #        self._cache_rss = None 
    993 #         
    994 #    def __call__(self, attr, data, weight_id=None): 
    995 #        ref = self._cache_data 
    996 #        if ref is not None and ref is data: 
    997 #            rss = self._cache_rss 
    998 #        else: 
    999 #            x, y = data.to_numpy_MA("1A/c") 
    1000 #            subsets, rss = subsets_selection_xtx_numpy(x, y) 
    1001 #            rss_diff = -numpy.diff(rss) 
    1002 #            rss = numpy.zeros_like(rss) 
    1003 #            for s_size in range(1, subsets.shape[0]): 
    1004 #                subset = subsets[s_size, :s_size + 1] 
    1005 #                rss[subset] += rss_diff[s_size - 1] 
    1006 #            rss = rss[1:] #Drop the intercept 
    1007 #            self._cache_data = data 
    1008 #            self._cache_rss = rss 
    1009 ##        print sorted(zip(rss, data.domain.attributes)) 
    1010 #        index = list(data.domain.attributes).index(attr) 
    1011 #        return rss[index] 
     995class ScoreRSS(scoring.Score): 
     996     
     997    handles_discrete = False 
     998    handles_continuous = True 
     999    computes_thresholds = False 
     1000     
     1001    def __new__(cls, attr=None, data=None, weight_id=None, **kwargs): 
     1002        self = scoring.Score.__new__(cls) 
     1003        if attr is not None and data is not None: 
     1004            self.__init__(**kwargs) 
     1005            # TODO: Should raise a warning, about caching 
     1006            return self.__call__(attr, data, weight_id) 
     1007        elif not attr and not data: 
     1008            return self 
     1009        else: 
     1010            raise ValueError("Both 'attr' and 'data' arguments expected.") 
     1011         
     1012    def __init__(self): 
     1013        self._cache_data = None 
     1014        self._cache_rss = None 
     1015         
     1016    def __call__(self, attr, data, weight_id=None): 
     1017        ref = self._cache_data 
     1018        if ref is not None and ref is data: 
     1019            rss = self._cache_rss 
     1020        else: 
     1021            x, y = data.to_numpy_MA("1A/c") 
     1022            try: 
     1023                subsets, rss = subset_selection_xtx(x, y) 
     1024            except numpy.linalg.LinAlgError: 
     1025                subsets, rss = subset_selection_xtx_numpy(x, y) 
     1026            rss_diff = -numpy.diff(rss) 
     1027            rss = numpy.zeros_like(rss) 
     1028            for s_size in range(1, subsets.shape[0]): 
     1029                subset = subsets[s_size, :s_size + 1] 
     1030                rss[subset] += rss_diff[s_size - 1] 
     1031            rss = rss[1:] #Drop the intercept 
     1032            self._cache_data = data 
     1033            self._cache_rss = rss 
     1034#        print sorted(zip(rss, data.domain.attributes)) 
     1035        index = list(data.domain.attributes).index(attr) 
     1036        return rss[index] 
    10121037         
    10131038 
  • source/orange/earth.cpp

    r8735 r8920  
    7373     - Added #include <limits> 
    7474     - Changed include of earth.h to earth.ppp and moved it before the module level defines 
     75     - Changed EvalSubsetsUsingXtX to return an error code if lin. dep. terms in bx 
    7576 
    7677    - TODO: Move global vars inside the functions using them (most are local) 
     
    24942495// The "Xtx" in the name refers to the X'X matrix. 
    24952496 
    2496 static void EvalSubsetsUsingXtx( 
     2497static int EvalSubsetsUsingXtx( 
    24972498    bool   PruneTerms[],    // out: nMaxTerms x nMaxTerms 
    24982499    double RssVec[],        // out: nMaxTerms x 1, RSS of each subset 
     
    25212522    for (int i = 0; i < nMaxTerms; i++) 
    25222523        WorkingSet[i] = true; 
     2524 
     2525    int error_code = 0; 
    25232526    for (int nUsedCols = nMaxTerms; nUsedCols > 0; nUsedCols--) { 
    25242527        int nRank; 
     
    25282531 
    25292532        if(nRank != nUsedCols) 
    2530             error("nRank %d != nUsedCols %d " 
    2531                 "(probably because of lin dep terms in bx)\n", 
    2532                 nRank, nUsedCols); 
     2533        { 
     2534            error_code = 1; 
     2535            break; 
     2536        } 
     2537//            error("nRank %d != nUsedCols %d " 
     2538//                "(probably because of lin dep terms in bx)\n", 
     2539//                nRank, nUsedCols); 
    25332540 
    25342541        RssVec[nUsedCols-1] = Rss; 
     
    25642571    free1(Diags); 
    25652572    free1(Betas); 
     2573 
     2574    return error_code; 
    25662575} 
    25672576 
     
    30193028} 
    30203029 
    3021 extern "C" void EarthEvalSubsetsUsingXtx( 
     3030extern "C" int EarthEvalSubsetsUsingXtx( 
    30223031    bool   PruneTerms[],    // out: nMaxTerms x nMaxTerms 
    30233032    double RssVec[],        // out: nMaxTerms x 1, RSS of each subset 
     
    30293038    const double WeightsArg[]) // in: nCases x 1, can be NULL 
    30303039{ 
    3031     EvalSubsetsUsingXtx(PruneTerms, RssVec, nCases, nResp, nMaxTerms, bx, y, WeightsArg); 
     3040    return EvalSubsetsUsingXtx(PruneTerms, RssVec, nCases, nResp, nMaxTerms, bx, y, WeightsArg); 
    30323041} 
    30333042 
  • source/orange/earth.hpp

    r8735 r8920  
    3838     - Removed definition for bool 
    3939     - Added extern "C" definitions for ForwardPass and EvalSubsetsUsingXtX 
     40     - Changed EvalSubsetsUsingXtX to return an error code if lin. dep. terms in bx 
    4041 
    4142 */ 
     
    204205    const char *sPredNames[]);   // in: predictor names, can be NULL 
    205206 
    206 EARTH_EXPORT void EarthEvalSubsetsUsingXtx( 
     207EARTH_EXPORT int EarthEvalSubsetsUsingXtx( 
    207208    bool   PruneTerms[],    // out: nMaxTerms x nMaxTerms 
    208209    double RssVec[],        // out: nMaxTerms x 1, RSS of each subset 
Note: See TracChangeset for help on using the changeset viewer.