source: orange/orange/orngProjectionPursuit.py @ 6390:b3b1f405d7eb

Revision 6390:b3b1f405d7eb, 7.8 KB checked in by gregor <gregor@…>, 4 years ago (diff)
  • projection pursuit indices
Line 
1import orange
2import numpy
3
4import scipy.special
5import scipy.optimize
6import scipy.stats
7
8from pylab import *
9
10def sqrtm(mat):
11    """ Retruns the square root of the matrix mat """
12    U, S, V = numpy.linalg.svd(mat)
13    D = numpy.diag(numpy.sqrt(S))
14    return numpy.dot(numpy.dot(U,D),V)
15
16def standardize(mat):
17    """ Subtracts means and multiplies by diagonal elements of inverse
18        square root of covariance matrix.
19    """
20    av = numpy.average(mat, axis=0)
21    sigma = numpy.corrcoef(mat, rowvar=0)
22    srSigma = sqrtm(sigma)
23    isrSigma = numpy.linalg.inv(srSigma)
24    return (mat-av) * numpy.diag(isrSigma)
25
26
27def friedman_tmp_func(alpha, Z=numpy.zeros((1,1)), J=5, n=1):
28    alpha = numpy.array(alpha)
29    pols = [scipy.special.legendre(j) for j in range(0,J+1)]
30    vals0 = [numpy.dot(alpha.T, Z[i,:]) for i in range(n)]
31    def f_tmp(x): return 2*x-1
32    vals = map(f_tmp, map(scipy.stats.zprob, vals0))
33    val = [1./n*sum(map(p, vals))**2 for p in pols]
34    return vals, pols, - 0.5 * sum([(2*j+1)*v for j, v in enumerate(val)])
35
36
37class ProjectionPursuit:
38    FRIEDMAN = 0
39    MOMENT = 1
40    SILHUETTE = 2
41    HARTINGAN = 3
42   
43    def __init__(self, data, index = FRIEDMAN, dim=2, maxiter=10):
44        self.dim = dim
45        if type(data) == orange.ExampleTable:
46            self.dataNP = data.toNumpy()[0]         # TODO: check if conversion of discrete values works ok
47        else:
48            self.dataNP = data
49        self.Z = standardize(self.dataNP)
50        self.totalSize, self.nVars = numpy.shape(self.Z)
51        self.maxiter = maxiter
52        self.currentOptimum = None
53        self.index = index
54
55    def optimize(self, maxiter = 5, opt_method=scipy.optimize.fmin):
56        func = self.getIndex()
57        if self.currentOptimum != None:
58            x = self.currentOptimum
59        else:
60            x = numpy.random.rand(self.dim * self.nVars)
61        alpha = opt_method(func, x, maxiter=maxiter).reshape(self.dim * self.nVars,1)
62        self.currentOptimum = alpha
63        print alpha, len(alpha)
64        optValue = func(alpha)
65        if self.dim == 2:
66            alpha1 = alpha[:self.nVars]
67            alpha2 = alpha[self.nVars:]
68            alpha = numpy.append(alpha1, alpha2, axis=1)
69        projectedData = numpy.dot(self.Z, alpha)
70        return alpha, optValue, projectedData
71
72    def find_optimum(self, opt_method=scipy.optimize.fmin):
73        func = self.getIndex()
74        alpha = opt_method(func, \
75                           numpy.random.rand(self.dim * self.nVars),\
76                           maxiter=self.maxiter).reshape(self.dim * self.nVars,1)
77        print alpha, len(alpha)
78        optValue = func(alpha)
79        if self.dim == 2:
80            alpha1 = alpha[:self.nVars]
81            alpha2 = alpha[self.nVars:]
82            alpha = numpy.append(alpha1, alpha2, axis=1)
83        projectedData = numpy.dot(self.Z, alpha)
84        return alpha, optValue, projectedData       
85
86    def getIndex(self):
87        if self.index == self.FRIEDMAN:
88            return self.getFriedmanIndex()
89        elif self.index == self.MOMENT:
90            return self.getMomentIndex()
91        elif self.index == self.SILHUETTE:
92            return self.getSilhouetteBasedIndex()
93        elif self.index == self.HARTINGAN:
94            return self.getHartinganBasedIndex()
95       
96
97    def getFriedmanIndex(self, J=5):
98        if self.dim == 1:
99            def func(alpha, Z=self.Z, J=J, n=self.totalSize):
100                vals, pols, val = friedman_tmp_func(alpha, Z=Z, J=J, n=n)
101                return val
102        elif self.dim == 2:
103            def func(alpha, Z=self.Z, J=J, n=self.totalSize):
104                alpha1, alpha2 = alpha[:self.nVars], alpha[self.nVars:]
105                vals1, pols, val1 = friedman_tmp_func(alpha1, Z=Z, J=J, n=n)
106                vals2, pols, val2 = friedman_tmp_func(alpha2, Z=Z, J=J, n=n)
107                val12 = - 0.5 * sum([sum([(2*j+1)*(2*k+1)*vals1[j]*vals2[k] for k in range(0, J+1-j)]) \
108                             for j in range(0,J+1)])
109##                print val1, val2
110                return 0.5 * (val1 + val2 + val12)
111        return func
112       
113
114    def getMomentIndex(self): # lahko dodas faktor 1./12
115        if self.dim == 1:
116            def func(alpha):
117                smpl = numpy.dot(self.Z, alpha)
118                return scipy.stats.kstat(smpl, n=3) ** 2 + 0.25 * scipy.stats.kstat(smpl, n=4)
119        else:
120            print "To do."
121        return func
122
123    def getSilhouetteBasedIndex(self, nClusters=5):
124        import orngClustering
125        def func(alpha, nClusters=nClusters):
126            alpha1, alpha2 = alpha[:self.nVars], alpha[self.nVars:]
127            alpha1 = alpha1.reshape((self.nVars,1))
128            alpha2 = alpha2.reshape(self.nVars,1)
129            alpha = numpy.append(alpha1, alpha2, axis=1)
130            smpl = numpy.dot(self.Z, alpha)
131            smpl = orange.ExampleTable(smpl)
132            km = orngClustering.KMeans(smpl, centroids=nClusters)
133            score = orngClustering.score_silhouette(km)
134            return -score
135        import functools
136        silhIndex = functools.partial(func, nClusters=nClusters)
137        return silhIndex
138       
139
140    def getHartinganBasedIndex(self, nClusters=5):
141        import orngClustering
142        def func(alpha, nClusters=nClusters):
143            alpha1, alpha2 = alpha[:self.nVars], alpha[self.nVars:]
144            alpha1 = alpha1.reshape((self.nVars,1))
145            alpha2 = alpha2.reshape(self.nVars,1)
146            alpha = numpy.append(alpha1, alpha2, axis=1)
147            smpl = numpy.dot(self.Z, alpha)
148            smpl = orange.ExampleTable(smpl)
149            km1 = orngClustering.KMeans(smpl, centroids=nClusters)
150            km2 = orngClustering.KMeans(smpl, centroids=nClusters)
151           
152            score = (self.totalSize - nClusters - 1) * (km1.score-km2.score) / (km2.score)
153            return -score
154        import functools
155        hartinganIndex = functools.partial(func, nClusters=nClusters)
156        return hartinganIndex
157               
158
159
160
161def draw_scatter_hist(x,y, fileName="lala.png"):
162    from matplotlib.ticker import NullFormatter
163    nullfmt   = NullFormatter()         # no labels
164
165    clf()
166
167    # definitions for the axes
168    left, width = 0.1, 0.65
169    bottom, height = 0.1, 0.65
170    bottom_h = left_h = left+width+0.02
171
172    rect_scatter = [left, bottom, width, height]
173    rect_histx = [left, bottom_h, width, 0.2]
174    rect_histy = [left_h, bottom, 0.2, height]
175
176    # start with a rectangular Figure
177    figure(1, figsize=(8,8))
178
179    axScatter = axes(rect_scatter)
180    axHistx = axes(rect_histx)
181    axHisty = axes(rect_histy)
182
183    # no labels
184    axHistx.xaxis.set_major_formatter(nullfmt)
185    axHisty.yaxis.set_major_formatter(nullfmt)
186
187    # the scatter plot:
188    axScatter.scatter(x, y)
189
190    # now determine nice limits by hand:
191    binwidth = 0.25
192    xymax = numpy.max([numpy.max(np.fabs(x)), numpy.max(np.fabs(y))])
193    lim = (int(xymax/binwidth) + 1) * binwidth
194
195    axScatter.set_xlim( (-lim, lim) )
196    axScatter.set_ylim( (-lim, lim) )
197
198    bins = numpy.arange(-lim, lim + binwidth, binwidth)
199    axHistx.hist(x, bins=bins)
200    axHisty.hist(y, bins=bins, orientation='horizontal')
201
202    axHistx.set_xlim(axScatter.get_xlim())
203    axHisty.set_ylim(axScatter.get_ylim())
204
205    savefig(fileName)
206
207
208if __name__=="__main__":
209##    data = orange.ExampleTable("c:\\Work\\Subgroup discovery\\iris.tab")
210    data = orange.ExampleTable(r"E:\Development\Orange Datasets\UCI\iris.tab")
211         
212    data = data.select(data.domain.attributes)
213
214    impmin = orange.ImputerConstructor_minimal(data) 
215    data = impmin(data) 
216
217    ppy = ProjectionPursuit(data, dim=2, maxiter=100)
218    #ppy.friedman_index(J=5)
219    #ppy.silhouette_based_index(nClusters=2)
220
221##    import os
222##    os.chdir("C:\\Work\\Subgroup discovery")
223    #draw_scatter_hist(ppy.friedmanProjData[:,0], ppy.friedmanProjData[:,1])
224    #draw_scatter_hist(ppy.silhouetteProjData[:,0], ppy.silhouetteProjData[:,1])
225
226    print ppy.optimize()
Note: See TracBrowser for help on using the repository browser.