source: orange/orange/orngDimReduction.py @ 6538:a5f65d7f0b2c

Revision 6538:a5f65d7f0b2c, 9.6 KB checked in by Mitar <Mitar@…>, 4 years ago (diff)

Made XPM version of the icon 32x32.

Line 
1import orange
2import numpy
3import numpy.ma as MA
4import numpy.linalg as linalg
5
6
7class Base:
8    def getSorted(self, values, vectors):
9        pairs = [(val, i) for i, val in enumerate(values)]
10        pairs.sort()
11        pairs.reverse()
12        indices = [pair[1] for pair in pairs]
13        newValues = [values[i] for i in indices]
14        newVectors = MA.take(vectors, indices, axis = 1)
15        return newValues, newVectors
16       
17    # substract average value to get zero mean
18    def center(self, data):
19        suma = data.sum(axis=0)/float(len(data))
20        data -= suma                           
21        return data
22   
23    # normalize variance to 1 for each attribute
24    def normalize(self, data):
25        data /= MA.std(data, axis=0)      # divide with standard deviation to get deviation = 1 for all attrs
26        return data
27
28class PCA(Base):
29    def __init__(self, data = None):
30        self.setData(data)
31       
32    def setData(self, data):
33        self.data = data
34        self.normalizedData = None
35        self.eigVectors = None
36        self.eigValues = None
37               
38           
39    def compute(self):
40        if not self.data:
41            return
42        if type(self.eigVectors) == MA.MaskedArray and type(self.eigValues) == MA.MaskedArray:
43            return
44       
45        if type(self.data) == orange.ExampleTable:
46            data = self.data.toNumpyMA("a")[0]
47        else:
48            data = self.data
49           
50        data = self.center(data)
51        data = self.normalize(data)
52        self.normalizedData = data      #
53       
54        covMatrix = MA.dot(data.T, data) / len(data)
55        eigVals, eigVectors = linalg.eigh(covMatrix)
56        self.eigValues, self.eigVectors = self.getSorted(eigVals, eigVectors)
57               
58       
59    def getCount(self, nPCs = None, varianceExplained = None):
60        if nPCs == None and varianceExplained == None:
61            nPCs = len(self.eigValues)
62        if varianceExplained == None:
63            nPCs = min(nPCs, len(self.eigVectors))
64        else:
65            total = 0; nPCs = 0
66            while total < varianceExplained and ind < len(self.eigValues):
67                total += self.eigValues[ind]
68                nPCs += 1
69        return nPCs
70       
71    def getEigValues(self, nPCs = None, varianceExplained = None):
72        if not self.data: return None
73        self.compute()
74       
75        nPCs = self.getCount(nPCs, varianceExplained)
76        return self.eigValues[:nPCs]
77
78
79    def getEigVectors(self, nPCs = None, varianceExplained = None):
80        if not self.data: return None
81        self.compute()
82       
83        nPCs = self.getCount(nPCs, varianceExplained)
84        return MA.take(self.eigVectors, range(nPCs), axis=1)
85       
86       
87    def getProjectedData(self, nPCs = None, varianceExplained = None):
88        if not self.data: return None
89        self.compute()
90       
91        vectors = self.getEigVectors(nPCs, varianceExplained)
92        return MA.dot(self.normalizedData, vectors)
93   
94    def getExampleTable(self, nPCs = None, varianceExplained = None, keepOriginal = 1):
95        if not self.data or type(self.data) != orange.ExampleTable: return None
96        count = self.getCount(nPCs, varianceExplained)
97        projData = self.getProjectedData(nPCs = count)
98        newDomain = orange.Domain([orange.FloatVariable("PC %d" % (d+1)) for d in range(count)], 0)
99        newTable = orange.ExampleTable(newDomain, projData.data)
100        if keepOriginal:
101            return orange.ExampleTable([self.data, newTable])
102        else:
103            if self.data.domain.classVar:
104                return orange.ExampleTable([self.data.select([self.data.domain.classVar]), newTable])
105            else:
106                return newTable
107       
108# ###########################################################
109# Fisher discriminant analysis
110class FDA(Base):
111    def __init__(self, data = None):
112        self.setData(data)
113       
114    def setData(self, data):
115        if data != None and type(data) not in [orange.ExampleTable, tuple]:
116            print "invalid data type. Only ExampleTable and tuple types are allowed"
117            self.data = None
118            return
119        self.data = data
120        self.normalizedData = None
121        self.eigVectors = None
122        self.eigValues = None
123       
124    def compute(self):
125        if self.data == None:
126            return
127        if type(self.eigVectors) == MA.MaskedArray and type(self.eigValues) == MA.MaskedArray:
128            return
129       
130        if type(self.data) == orange.ExampleTable:
131            data, classes = self.data.toNumpyMA("a/c")
132        elif type(self.data) == tuple:
133            data, classes = self.data
134
135        data = self.center(data)
136        data = self.normalize(data)
137        self.normalizedData = data
138        exampleCount, attrCount = data.shape
139        classCount = len(set(classes))
140        # special case when we have two classes
141        if classCount == 2:
142            data1 = MA.take(data, numpy.argwhere(classes == 0).flatten(), axis=0)
143            data2 = MA.take(data, numpy.argwhere(classes != 0).flatten(), axis=0)
144            miDiff = MA.average(data1, axis=1) - MA.average(data2, axis=1)
145            covMatrix = (MA.dot(data1.T, data1) + MA.dot(data2.T, data2)) / exampleCount
146            self.eigVectors = linalg.inv(covMatrix) * miDiff
147            self.eigValues = numpy.array([1])
148        else:
149            # compute means and average covariances of examples in each class group
150            Sw = MA.zeros([attrCount, attrCount])
151            for v in set(classes):
152                d = MA.take(data, numpy.argwhere(classes == v).flatten(), axis=0)
153                d = self.center(d)
154                Sw += MA.dot(d.T, d)
155            Sw /= exampleCount
156            total = MA.dot(data.T, data)/float(exampleCount)
157            Sb = total - Sw
158                       
159            matrix = linalg.inv(Sw)*Sb
160            eigVals, eigVectors = linalg.eigh(matrix)
161            self.eigValues, self.eigVectors = self.getSorted(eigVals, eigVectors)
162           
163    def getCount(self, nPCs = None, varianceExplained = None):
164        if nPCs == None and varianceExplained == None:
165            nPCs = len(self.eigValues)
166        if varianceExplained == None:
167            nPCs = min(nPCs, len(self.eigVectors))
168        else:
169            total = 0; nPCs = 0
170            while total < varianceExplained and ind < len(self.eigValues):
171                total += self.eigValues[ind]
172                nPCs += 1
173        return nPCs
174       
175    def getEigValues(self, nPCs = None, varianceExplained = None):
176        if not self.data: return None
177        self.compute()
178       
179        nPCs = self.getCount(nPCs, varianceExplained)
180        return self.eigValues[:nPCs]
181
182
183    def getEigVectors(self, nPCs = None, varianceExplained = None):
184        if not self.data: return None
185        self.compute()
186       
187        nPCs = self.getCount(nPCs, varianceExplained)
188        return MA.take(self.eigVectors, range(nPCs), axis=1)
189
190           
191    def getProjectedData(self, nPCs = None, varianceExplained = None):
192        if not self.data: return None
193        self.compute()
194       
195        vectors = self.getEigVectors(nPCs, varianceExplained)
196        return MA.dot(self.normalizedData, vectors)
197   
198    def getExampleTable(self, nPCs = None, varianceExplained = None, keepOriginal = 1):
199        if not self.data or type(self.data) != orange.ExampleTable: return None
200        count = self.getCount(nPCs, varianceExplained)
201        projData = self.getProjectedData(nPCs = count)
202        newDomain = orange.Domain([orange.FloatVariable("Component %d" % (d+1)) for d in range(count)], 0)
203        newTable = orange.ExampleTable(newDomain, projData.data)
204        if keepOriginal:
205            return orange.ExampleTable([self.data, newTable])
206        else:
207            if self.data.domain.classVar:
208                return orange.ExampleTable([self.data.select([self.data.domain.classVar]), newTable])
209            else:
210                return newTable
211
212       
213def PCAOnExampleTable(table, keepOriginal = 1, nPCs = -1):
214    data = table.toNumpyMA("a")[0]
215    projData, vectors, values = pca(data, nPCs)
216    newDomain = orange.Domain([orange.FloatVariable("PC %d" % (d+1)) for d in range(len(vectors))], 0)
217    newTable = orange.ExampleTable(newDomain, projData.data)
218    if keepOriginal:
219        return orange.ExampleTable([table, newTable])
220    else:
221        return newTable
222   
223
224def pca(data, nPCs = -1):
225    domain = None
226   
227    suma = data.sum(axis=0)/float(len(data))
228    data -= suma       # substract average value to get zero mean
229    data /= MA.std(data, axis=0)
230    covMatrix = MA.dot(data.T, data) / len(data)
231
232    eigVals, eigVectors = linalg.eigh(covMatrix)
233    eigVals = list(eigVals)
234   
235    if nPCs == -1:
236        nPCs = len(eigVals)
237    nPCs = min(nPCs, len(eigVals))
238   
239    pairs = [(val, i) for i, val in enumerate(eigVals)]
240    pairs.sort()
241    pairs.reverse()
242    indices = [pair[1] for pair in pairs[:nPCs]]  # take indices of the wanted number of principal components
243
244    vectors = MA.take(eigVectors, indices, axis = 1)
245    values = [eigVals[i] for i in indices]
246    projectedData = MA.dot(data, vectors)
247   
248    return projectedData, vectors, values
249
250
251if __name__ == "__main__":
252    #d = orange.Domain([orange.FloatVariable("X1"), orange.FloatVariable("X2")], 0)
253    #data = orange.ExampleTable(d, [[x,x] for x in range(100)])
254    data = orange.ExampleTable(r"E:\Development\Orange Datasets\UCI\iris.tab")
255    x = PCAOnExampleTable(data)
256
257    data2 = data.toNumpyMA("a")[0]
258    projData, vectors, values = pca(data2)
259    print vectors
260    print values
261   
262    instance = PCA(data)
263    print instance.getEigValues()
264    print instance.getEigVectors()
265    print instance.getProjectedData()
266   
267    fda = FDA(data)
268    print fda.getProjectedData()
Note: See TracBrowser for help on using the repository browser.