Changeset 9676:895804462430 in orange for Orange/orng/orngPCA.py


Ignore:
Timestamp:
02/06/12 11:47:54 (2 years ago)
Author:
Miha Stajdohar <miha.stajdohar@…>
Branch:
default
Parents:
9675:aa9d846ae025 (diff), 9672:3b64ad491c7e (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merged.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • Orange/orng/orngPCA.py

    r9671 r9676  
    1717    warnings.warn("Importing of matplotlib has failed.\nPlotting functions will be unavailable.") 
    1818    mathlib_import = False 
    19      
     19 
    2020#color table for biplot 
    21 Colors = ['bo','go','yo','co','mo'] 
     21Colors = ['bo', 'go', 'yo', 'co', 'mo'] 
    2222 
    2323def defaultImputer(dataset): 
     
    4141def centerData(dataMatrix): 
    4242    """Perfomrs centering od data along rows, returns center and centered data.""" 
    43     center = numpy.sum(dataMatrix, axis = 0)/float(len(dataMatrix)) 
     43    center = numpy.sum(dataMatrix, axis=0) / float(len(dataMatrix)) 
    4444    return center, (dataMatrix - center) 
    4545 
     
    4747    """Performs standardization of data along rows. Throws error if constant 
    4848    variable is present.""" 
    49     scale = numpy.std(dataMatrix, axis = 0) 
     49    scale = numpy.std(dataMatrix, axis=0) 
    5050    if 0. in scale: 
    5151        raise orange.KernelException, "Constant variable, cannot standardize!" 
    52     return scale, dataMatrix * 1./scale 
     52    return scale, dataMatrix * 1. / scale 
    5353 
    5454class PCA(object): 
     
    8888    PCAClassifier object if projection was successful or None if projection has failed 
    8989    """ 
    90      
    91     def __new__(cls, dataset = None, attributes = None, rows = None, standardize = 0, 
    92                  imputer = defaultImputer, continuizer = defaultContinuizer, 
    93                  maxNumberOfComponents = 10, varianceCovered = 0.95, 
    94                  useGeneralizedVectors = 0): 
    95          
     90 
     91    def __new__(cls, dataset=None, attributes=None, rows=None, standardize=0, 
     92                 imputer=defaultImputer, continuizer=defaultContinuizer, 
     93                 maxNumberOfComponents=10, varianceCovered=0.95, 
     94                 useGeneralizedVectors=0): 
     95 
    9696        learner = object.__new__(cls) 
    97          
     97 
    9898        if dataset: 
    9999            learner.__init__(attributes, rows, standardize, imputer, continuizer, 
     
    102102        else: 
    103103            return learner 
    104      
     104 
    105105    #decide what should be inmutable and what not 
    106106    #imutable for sure: imputer and continuizer 
    107     def __init__(self, attributes = None, rows = None, standardize = 0, 
    108                  imputer = defaultImputer, continuizer = defaultContinuizer, 
    109                  maxNumberOfComponents = 10, varianceCovered = 0.95, 
    110                  useGeneralizedVectors = 0): 
    111          
     107    def __init__(self, attributes=None, rows=None, standardize=0, 
     108                 imputer=defaultImputer, continuizer=defaultContinuizer, 
     109                 maxNumberOfComponents=10, varianceCovered=0.95, 
     110                 useGeneralizedVectors=0): 
     111 
    112112        self.attributes = attributes 
    113113        self.rows = rows 
     
    118118        self.varianceCovered = varianceCovered 
    119119        self.useGeneralizedVectors = useGeneralizedVectors 
    120      
     120 
    121121    def __call__(self, dataset): 
    122          
     122 
    123123        #Modify dataset 
    124124        dataset = self._datasetPreprocessing(dataset) 
    125125        dataMatrix, classArray, center, deviation = self._dataMatrixPreprocessing(dataset) 
    126          
     126 
    127127        #Perform projection 
    128128        evalues, evectors = self._createPCAProjection(dataMatrix, classArray) 
    129          
     129 
    130130        #check if return status is None, None 
    131131        if (evalues, evectors) == (None, None): 
    132132            print "Principal component could not be performed (complex eigenvalues or singular matrix if generalized eigenvectors were used)" 
    133133            return None 
    134          
     134 
    135135        #return PCAClassifier 
    136         return PCAClassifier(domain = self.attributes, 
    137                              imputer = self.imputer, 
    138                              continuizer = self.continuizer, 
    139                              center = center, 
    140                              deviation = deviation, 
    141                              evalues = evalues, 
    142                              loadings = evectors) 
     136        return PCAClassifier(domain=self.attributes, 
     137                             imputer=self.imputer, 
     138                             continuizer=self.continuizer, 
     139                             center=center, 
     140                             deviation=deviation, 
     141                             evalues=evalues, 
     142                             loadings=evectors) 
    143143 
    144144    def _datasetPreprocessing(self, dataset): 
     
    151151        if self.attributes: 
    152152            dataset = dataset.select(self.attributes) 
    153          
     153 
    154154        #we need to retain only selected attributes without class attribute     
    155155        self.attributes = [att.name for att in dataset.domain.attributes] 
    156156 
    157157        imputer = self.imputer(dataset) 
    158              
     158 
    159159        if self.rows: 
    160160            dataset = dataset.select(self.rows) 
     
    165165        dataset = dataset.translate(domain) 
    166166 
    167         return dataset         
     167        return dataset 
    168168 
    169169    def _dataMatrixPreprocessing(self, dataset): 
     
    175175 
    176176        dataMatrix, classArray, x = dataset.toNumpy() 
    177          
     177 
    178178        center, dataMatrix = centerData(dataMatrix) 
    179          
     179 
    180180        deviation = None 
    181181        if self.standardize: 
     
    183183 
    184184        return dataMatrix, classArray, center, deviation 
    185      
     185 
    186186    def _createPCAProjection(self, dataMatrix, classArray): 
    187187        """ 
     
    195195        Snapshot with generalized vectors: M * (t(M) * M)^(-1) * t(M) * L * M * x = lambda * M * x 
    196196        """ 
    197          
     197 
    198198        n, d = numpy.shape(dataMatrix) 
    199          
     199 
    200200        if classArray != None: 
    201201            L = numpy.zeros((len(dataMatrix), len(dataMatrix))) 
    202202            for i in range(len(dataMatrix)): 
    203                 for j in range(i+1, len(dataMatrix)): 
    204                     L[i,j] = -int(classArray[i] != classArray[j]) 
    205                     L[j,i] = -int(classArray[i] != classArray[j]) 
     203                for j in range(i + 1, len(dataMatrix)): 
     204                    L[i, j] = -int(classArray[i] != classArray[j]) 
     205                    L[j, i] = -int(classArray[i] != classArray[j]) 
    206206 
    207207            s = numpy.sum(L, axis=0)      # doesn't matter which axis since the matrix L is symmetrical 
    208208            for i in range(len(dataMatrix)): 
    209                 L[i,i] = -s[i] 
    210                  
     209                L[i, i] = -s[i] 
     210 
    211211            matrix = dot(transpose(dataMatrix), L) 
    212212        else: 
    213213            matrix = transpose(dataMatrix) 
    214          
     214 
    215215        if self.useGeneralizedVectors: 
    216216            temp_matrix = dot(transpose(dataMatrix), dataMatrix) 
     
    220220                return None, None 
    221221            matrix = dot(temp_matrix, matrix) 
    222          
     222 
    223223        if n < d: 
    224224            #snapshot method 
     
    232232            evalues = evalues[positiveArray] 
    233233            evectors = evectors[positiveArray] 
    234             evectors = transpose(1./sqrt(evalues)) * transpose(dot(evectors, dataMatrix)) 
     234            evectors = transpose(1. / sqrt(evalues)) * transpose(dot(evectors, dataMatrix)) 
    235235        else: 
    236236            covMatrix = dot(matrix, dataMatrix) 
     
    243243            evalues = evalues[positiveArray] 
    244244            evectors = evectors[positiveArray] 
    245                          
     245 
    246246        order = (numpy.argsort(evalues)[::-1]) 
    247247        N = len(evalues) 
    248248        maxComp = self.maxNumberOfComponents 
    249249        variance = self.varianceCovered 
    250          
     250 
    251251        if maxComp == -1: 
    252252            maxComp = N 
     
    259259        #                   -> value < 1 : explains less than previous attribute 
    260260        for i in range(len(order)): 
    261             variance -= evalues[i] / trace  
     261            variance -= evalues[i] / trace 
    262262            if variance < 0: 
    263263                return evalues[: i + 1] * scale, numpy.take(evectors, order[: i + 1], 1) 
    264          
    265         return evalues * scale, numpy.take(evectors, order, 1)  
     264 
     265        return evalues * scale, numpy.take(evectors, order, 1) 
    266266 
    267267class PCAClassifier(object): 
    268      
     268 
    269269    def __init__(self, domain, imputer, continuizer, center, deviation, evalues, loadings): 
    270270        #data checking and modifying 
     
    277277        self.evalues = evalues 
    278278        self.loadings = loadings 
    279          
     279 
    280280        #last predicition performed -> used for biplot 
    281281        self._dataMatrix = None 
    282282        self._classArray = None 
    283      
     283 
    284284    def __call__(self, dataset): 
    285                  
     285 
    286286        try: 
    287287            #retain class attribute 
     
    298298        dataMatrix -= self.center 
    299299        if self.deviation != None: 
    300             dataMatrix *= 1./self.deviation 
    301              
     300            dataMatrix *= 1. / self.deviation 
     301 
    302302        #save transformed data 
    303303        self._dataMatrix = numpy.dot(dataMatrix, self.loadings) 
    304304 
    305         attributes = [orange.FloatVariable("PC%d" % (i + 1, )) for i in range(len(self.evalues))] 
     305        attributes = [orange.FloatVariable("PC%d" % (i + 1,)) for i in range(len(self.evalues))] 
    306306        new_domain = orange.Domain(attributes) 
    307307        new_table = orange.ExampleTable(new_domain, self._dataMatrix) 
     
    312312            self._classArray = numpy.array([row.getclass() for row in classTable]) 
    313313            new_table = orange.ExampleTable([new_table, classTable]) 
    314          
     314 
    315315        return new_table 
    316      
     316 
    317317    def __str__(self): 
    318          
     318 
    319319        n, d = numpy.shape(self.loadings) 
    320320        comulative = 0.0 
    321          
     321 
    322322        summary = "PCA SUMMARY\n\nCenter:\n\n" 
    323323        summary += " %15s  " * len(self.domain) % tuple(self.domain) + "\n" 
    324324        summary += " %15.4f  " * len(self.center) % tuple(self.center) + "\n\n" 
    325      
     325 
    326326        if self.deviation != None: 
    327327            summary += "Deviation:\n\n" 
    328328            summary += " %15s  " * len(self.domain) % tuple(self.domain) + "\n" 
    329329            summary += " %15.4f  " * len(self.deviation) % tuple(self.deviation) + "\n\n" 
    330     
     330 
    331331        summary += "Importance of components:\n\n %12s   %12s   %12s\n" % ("eigenvalues", "proportion", "cumulative") 
    332332        for evalue in self.evalues: 
    333             comulative += evalue/n 
    334             summary += " %12.4f   %12.4f   %12.4f\n" % (evalue, evalue/n, comulative) 
    335          
    336         summary += "\nLoadings:\n\n" 
    337         summary += "      PC%d" * d % tuple(range(1, d + 1)) + "\n" 
    338         for attr_num in range(len(self.domain)): 
    339             summary += " % 8.4f" * d % tuple(self.loadings[attr_num]) 
    340             summary += "   %-30s\n" % self.domain[attr_num] 
     333            comulative += evalue / n 
     334            summary += " %12.4f   %12.4f   %12.4f\n" % (evalue, evalue / n, comulative) 
     335 
     336        #summary += "\nLoadings:\n\n" 
     337        #summary += "      PC%d" * d % tuple(range(1, d + 1)) + "\n" 
     338        #for attr_num in range(len(self.domain)): 
     339        #    summary += " % 8.4f" * d % tuple(self.loadings[attr_num]) 
     340        #    summary += "   %-30s\n" % self.domain[attr_num] 
    341341 
    342342        return summary 
    343      
     343 
    344344    ################ Ploting functions ################### 
    345      
    346     def plot(self, title = 'Scree plot', filename = 'scree_plot.png'): 
     345 
     346    def plot(self, title='Scree plot', filename='scree_plot.png'): 
    347347        """ 
    348348        Draws a scree plot. Matplotlib is needed. 
     
    354354            If None, plot will be shown 
    355355        """ 
    356          
     356 
    357357        if not mathlib_import: 
    358358            raise orange.KernelException, "Matplotlib was not imported!" 
    359          
     359 
    360360        #plt.clf() -> opens two windows 
    361361        fig = plt.figure() 
    362362        ax = fig.add_subplot(111) 
    363          
     363 
    364364        x_axis = range(len(self.evalues)) 
    365         x_labels = ["PC%d" % (i + 1, ) for i in x_axis] 
    366          
     365        x_labels = ["PC%d" % (i + 1,) for i in x_axis] 
     366 
    367367        ax.set_xticks(x_axis) 
    368368        ax.set_xticklabels(x_labels) 
     
    370370        ax.set_ylabel('Variance') 
    371371        ax.set_title(title + "\n") 
    372         ax.bar(left = x_axis, height = self.evalues, align = 'center') 
    373         ax.axis([-0.5, len(self.evalues) - 0.5, 0, self.evalues[0]*1.05]) 
     372        ax.bar(left=x_axis, height=self.evalues, align='center') 
     373        ax.axis([-0.5, len(self.evalues) - 0.5, 0, self.evalues[0] * 1.05]) 
    374374        if filename: 
    375375            plt.savefig(filename) 
    376376        else: 
    377377            plt.show() 
    378              
    379     def biplot(self, choices = [1,2], scale = 1., title = 'Biplot', 
    380                filename = 'biplot.png'): 
     378 
     379    def biplot(self, choices=[1, 2], scale=1., title='Biplot', 
     380               filename='biplot.png'): 
    381381        """ 
    382382        Draws biplot for PCA. Matplotlib is needed. Actual projection must be 
     
    392392            If None, plot will be shown 
    393393        """ 
    394          
     394 
    395395        if not mathlib_import: 
    396396            raise orange.KernelException, "Matplotlib was not imported!" 
    397          
     397 
    398398        if self._dataMatrix == None: 
    399399            raise orange.KernelException, "No data available for biplot!" 
    400          
     400 
    401401        if len(choices) != 2: 
    402402            raise orange.KernelException, 'You have to choose exactly two components' 
    403          
     403 
    404404        if max(choices[0], choices[1]) > len(self.evalues) or min(choices[0], choices[1]) < 1: 
    405405            raise orange.KernelException, 'Invalid choices' 
    406          
     406 
    407407        choice = numpy.array([i == choices[0] - 1 or i == choices[1] - 1 for i in range(len(self.evalues))]) 
    408          
    409         dataMatrix = self._dataMatrix[:,choice] 
    410         loadings = self.loadings[:,choice] 
     408 
     409        dataMatrix = self._dataMatrix[:, choice] 
     410        loadings = self.loadings[:, choice] 
    411411        lam = (numpy.array(self.evalues)[choice]) 
    412412        lam *= sqrt(len(self._dataMatrix)) 
    413          
     413 
    414414        if scale < 0. or scale > 1.: 
    415415            print "Warning: 'scale' is outside [0, 1]" 
    416         lam = lam**scale 
    417          
     416        lam = lam ** scale 
     417 
    418418        #TO DO -> pc.biplot (maybe) 
    419419        trDataMatrix = dataMatrix / lam 
    420420        trLoadings = loadings * lam 
    421          
     421 
    422422        max_data_value = numpy.max(abs(trDataMatrix)) * 1.05 
    423423        max_load_value = numpy.max(abs(trLoadings)) * 1.5 
    424          
     424 
    425425        #plt.clf() 
    426426        fig = plt.figure() 
     
    432432        ax1.xaxis.set_ticks_position('bottom') 
    433433        ax1.yaxis.set_label_position('left') 
    434         ax1.yaxis.set_ticks_position('left')         
    435          
     434        ax1.yaxis.set_ticks_position('left') 
     435 
    436436        if self._classArray == None: 
    437437            trDataMatrix = transpose(trDataMatrix) 
     
    447447                partialDataMatrix = transpose(trDataMatrix[choice]) 
    448448                ax1.plot(partialDataMatrix[0], partialDataMatrix[1], 
    449                          Colors[i % len(Colors)], label = str(classValues[i])) 
     449                         Colors[i % len(Colors)], label=str(classValues[i])) 
    450450            ax1.legend() 
    451          
     451 
    452452        ax1.set_xlim(-max_data_value, max_data_value) 
    453453        ax1.set_ylim(-max_data_value, max_data_value) 
    454          
     454 
    455455        #eliminate double axis on right 
    456456        ax0 = ax1.twinx() 
    457457        ax0.yaxis.set_visible(False) 
    458                  
     458 
    459459        ax2 = ax0.twiny() 
    460460        ax2.xaxis.set_label_position('top') 
     
    466466        for tl in ax2.get_yticklabels(): 
    467467            tl.set_color('r') 
    468          
    469         arrowprops = dict(facecolor = 'red', edgecolor = 'red', width = 1, headwidth = 4) 
     468 
     469        arrowprops = dict(facecolor='red', edgecolor='red', width=1, headwidth=4) 
    470470        #using annotations instead of arrows because there is a strange implementation 
    471471        #of arrows in matplotlib version 0.99 
    472472        for i in range(len(trLoadings)): 
    473473            x, y = trLoadings[i] 
    474             ax2.annotate('', (x, y), (0, 0), arrowprops = arrowprops) 
    475             ax2.text(x * 1.1, y * 1.2, self.domain[i], color = 'red') 
    476              
     474            ax2.annotate('', (x, y), (0, 0), arrowprops=arrowprops) 
     475            ax2.text(x * 1.1, y * 1.2, self.domain[i], color='red') 
     476 
    477477        ax2.set_xlim(-max_load_value, max_load_value) 
    478478        ax2.set_ylim(-max_load_value, max_load_value) 
    479          
     479 
    480480        if filename: 
    481481            plt.savefig(filename) 
Note: See TracChangeset for help on using the changeset viewer.