source: orange/orange/orngPade.py @ 8042:ffcb93bc9028

Revision 8042:ffcb93bc9028, 34.2 KB checked in by markotoplak, 3 years ago (diff)

Hierarchical clustering: also catch RuntimeError when importing matplotlib (or the documentation could not be built on server).

Line 
1import orange, statc
2import math, numpy, string, os
3
4from orangeom import star, dist
5from copy import deepcopy
6
7#pathQHULL = r"c:\qhull"
8#pathQHULL = r"c:\D\ai\Orange\test\squin\qhull"
9if __name__ != "__main__":
10    pathQHULL = os.path.split(__file__)[0]
11
12class Cache:
13    pass
14
15
16def clearCache(cache):
17    cache.points = cache.tri = cache.stars = cache.dts = cache.deltas = cache.findNearest = cache.attrStat = None
18
19
20def makeEmptyCache(cache = None):
21    if not cache:
22        cache = Cache()
23    cache.data = None
24    cache.attributes = []
25    cache.contAttributes = []
26    cache.npoints = 0
27    clearCache(cache)
28    return cache
29
30
31def makeBasicCache(data, cache = None):
32    if not cache:
33        cache = Cache()
34
35    cache.data = data             
36    cache.npoints = len(data)
37
38    attributes = cache.data.domain.attributes           
39    cache.contIndices = [i for i, attr in enumerate(attributes) if attr.varType == orange.VarTypes.Continuous]
40    cache.contAttributes = [attributes[i] for i in cache.contIndices]
41    cache.attributes = [(attr.name, attr.varType) for attr in cache.contAttributes]
42
43    clearCache(cache)
44    return cache
45
46
47
48def triangulate(cache, points):
49    import orangeom
50    return orangeom.qhull(cache.data.toNumpy("a")[0][:, cache.contIndices]).tolist()
51#    return [map(int, string.split(l[i],' ')) for i in xrange(num_of_triangles+1) if l[i]!='']
52
53#
54def simplex_with_xnDBG(cache, xp, Star, d):
55    zm=0.0001 # zm ... zelo malo
56    r=[]
57    xnd = [xp[i] for i in xrange(len(xp)-1) if i!=d] # vse koordinate, razen d-te
58    xndz = xp[-1]
59    for s in Star:
60        for (p,j) in [(cache.points[i],i) for i in s]:
61            pd = [p[i] for i in xrange(len(p)-1) if i!=d] # vse koordinate, razen d-te
62            if reduce(lambda x,y: x+y,[math.fabs(i) for i in list(numpy.array(xnd)-numpy.array(pd))])<=zm and math.fabs(xp[d]-p[d])>zm:
63                r += [(math.fabs(xp[d]-p[d]),j)]
64    if len(r)>0:
65        r.sort()
66        xdt = cache.points[r[0][1]]
67        dt = xdt[d]-xp[d]
68        D = (xdt[-1]-xndz)/r[0][0]
69        if dt<0:
70            D = -D
71        return D
72    return None
73
74def simplex_with_xn(cache, xn, Star):
75    for simplex in Star:
76        bl = [numpy.linalg.det(a) for a in inside(cache, xn, simplex)]
77        if reduce(lambda x,y: x and y, [i<0 for i in bl]) or reduce(lambda x,y: x and y, [i>0 for i in bl]):
78            return simplex
79    return None
80
81def change(cache, i,j,n):
82    if i==j:
83        return n+[1]
84    return cache.points[j][:-1]+[1]
85
86def inside(cache, vertex,simplex):
87    return [numpy.array([change(cache, i,j,vertex) for j in simplex]) for i in simplex]
88
89
90def firstTriangle(cache, dimensions, progressCallback = None, **args):
91    if len(cache.contAttributes) == 1:
92        return triangles1D(cache, False)
93
94    if not cache.points:       
95        cache.points = orange.ExampleTable(orange.Domain(cache.contAttributes, cache.data.domain.classVar), cache.data).native(0)
96    points = cache.points
97    npoints = len(points)
98
99    if not cache.tri:
100        cache.tri = triangulate(cache, points)
101    tri = cache.tri
102       
103    if not cache.stars:
104        cache.stars = [star(x, tri) for x in xrange(npoints)]
105    S = cache.stars
106
107    if not cache.dts:       
108        cache.dts = [min([ min([ dist(points[x][:-1],points[v][:-1]) for v in simplex if v!=x]) for simplex in S[x]])*.1 for x in xrange(npoints)]
109
110
111    if progressCallback:
112        nPoints = 100.0/npoints
113       
114    for x, (S, xp, dt, deltas) in enumerate(zip(cache.stars, points, cache.dts, cache.deltas)):
115        for d in dimensions:
116            xn = xp[:-1]
117            DBG=0
118#            if xn[0]>16 and xn[1]<44.5 and xn[1]>43:
119#            #if xn[0]>4.7 and xn[0]<4.9 and xn[1]<24.5 and xn[1]>23.5:
120#                DBG=1
121#                #print "DBG"
122           
123            O = numpy.array(xp[:-1])
124
125            xn[d] += dt
126            swx = simplex_with_xn(cache, xn, S)
127            if swx:
128#                if DBG:
129#                    print "iskanje cudnih trikotnikov"
130#                    print swx
131#                    print [points[k] for k in swx]
132                obrni = 1
133#                if DBG:
134#                    print xp
135#                    print swx
136#                    print [points[k][-1] for k in swx]
137#                    vecs = numpy.array([numpy.array(points[p][:-1])-O for p in swx if p!=x])
138#                    vecs = vecs.transpose()
139#                    XN = numpy.array(xn)-O
140#                    coef = numpy.linalg.solve(vecs,XN)
141#                    xnz = sum(coef*[numpy.array(points[p][-1]-xp[-1])for p in swx if p!=x])+xp[-1]
142#                    dd = obrni * (xnz-xp[-1]) / dt
143#                    print "xnz= ",xnz,"\tdd=",dd,"\trazlika (xnz-xp[-1])=",(xnz-xp[-1])
144#                    print "-----------------"
145            else:
146                xn[d] = xp[d]-dt
147                swx = simplex_with_xn(cache, xn, S)
148                if swx:
149                    obrni = -1
150#                    if DBG:
151#                        print xp
152#                        print swx
153#                        print "v levo\t",[points[k] for k in swx]
154#                        print "----------------"
155                else:
156#                    if DBG:
157#                        print xp
158#                        do = simplex_with_xnDBG(cache, xp, S, d)
159#                        print do
160                    do = simplex_with_xnDBG(cache, xp, S, d)
161                    if do:
162                        deltas[d] = do
163                    else:
164                        deltas[d] = "?"
165                    continue
166           
167            vecs = numpy.array([numpy.array(points[p][:-1])-O for p in swx if p!=x])
168            vecs = vecs.transpose()
169            XN = numpy.array(xn)-O
170            coef = numpy.linalg.solve(vecs,XN)
171            xnz = sum(coef*[numpy.array(points[p][-1]-xp[-1])for p in swx if p!=x])+xp[-1]
172
173            deltas[d] = obrni * (xnz-xp[-1]) / dt
174            #print "DELTAS = ",deltas[d]
175
176        if progressCallback:
177            progressCallback(x*nPoints)
178
179
180# calculates a linear regression on the star
181def starRegression(cache, dimensions, progressCallback=None, **args):
182    if len(cache.contAttributes) == 1:
183        return triangles1D(cache, True)
184
185    if not cache.points:       
186        cache.points = orange.ExampleTable(orange.Domain(cache.contAttributes, cache.data.domain.classVar), cache.data).native(0)
187    points = cache.points
188    npoints = len(points)
189
190    if not cache.tri:
191        cache.tri = triangulate(cache, points)
192    tri = cache.tri
193       
194    if not cache.stars:
195        cache.stars = [star(x, tri) for x in xrange(npoints)]
196    S = cache.stars
197
198    points = cache.points
199
200    if progressCallback:
201        nPoints = 100.0/len(points)
202       
203    for x,(S,p) in enumerate(zip(cache.stars,points)):
204        if S==[]:
205            cache.deltas[x] = ['?' for i in dimensions]
206            continue
207        st  =list(set(reduce(lambda x,y: x+y, S)))
208        A = [points[i][:-1] for i in st]
209        b = [[points[i][-1]] for i in st]
210        cache.deltas[x] = [i[0] for i in numpy.linalg.lstsq(A, b)[0]]
211
212        if progressCallback:
213            progressCallback(x*nPoints)
214
215       
216# calculates a univariate linear regression on the star
217def starUnivariateRegression(cache, dimensions, progressCallback = None, **args):
218    if len(cache.contAttributes) == 1:
219        return triangles1D(cache, True)
220       
221    if not cache.points:       
222        cache.points = orange.ExampleTable(orange.Domain(cache.contAttributes, cache.data.domain.classVar), cache.data).native(0)
223    points = cache.points
224    npoints = len(points)
225
226    if not cache.tri:
227        cache.tri = triangulate(cache, points)
228    tri = cache.tri
229       
230    if not cache.stars:
231        cache.stars = [star(x, tri) for x in xrange(npoints)]
232    S = cache.stars
233
234    points = cache.points
235
236    if progressCallback:
237        nPoints = 100.0/len(points)
238       
239    for x,(S,p) in enumerate(zip(cache.stars,points)):
240        if S==[]:
241            cache.deltas[x] = ['?' for i in dimensions]
242            continue
243        st = list(set(reduce(lambda x,y: x+y, S)))
244        lenst = len(st)
245        avgy = sum([points[i][-1] for i in st])/lenst
246        for d in dimensions:
247            avgx = sum([points[i][d] for i in st])/lenst
248            sxx2 = sum([(points[i][d]-avgx)**2 for i in st])
249            if sxx2:
250                sxx = sum([(points[i][d]-avgx)*(points[i][-1]-avgy) for i in st])
251                b = sxx/sxx2
252                cache.deltas[x][d] = b
253            else:
254                cache.deltas[x][d] = '?'
255
256        if progressCallback:
257            progressCallback(x*nPoints)
258
259
260def triangles1D(cache, bothWays):
261    data = cache.data
262    attrIdx = cache.contIndices[0]
263    points = [(ex[attrIdx], float(ex.getclass()), [i]) for i, ex in enumerate(data) if not ex.getclass().isSpecial()]
264    points.sort()
265
266    i, e = 0, len(points)
267    while i < e:
268        ex = points[i]
269        ni = i+1
270        while ni < e and points[ni][0] == ex[0]:
271            ex[2] += points[ni][2]
272            ex[1] += points[ni][1]
273            ni += 1
274        i = ni
275    points = [(a, c/len(i), i) for a, c, i in points]       
276
277    for delta in cache.deltas:
278        delta[0] = "?"
279
280    if len(points) < 2:
281        return
282
283    der = (points[-1][1] - points[-2][1]) / (points[-1][0] - points[-2][0])
284    for pidx in points[-1][2]:
285        cache.deltas[pidx][0] = der
286
287    if bothWays:
288        der = (points[1][1] - points[0][1]) / (points[1][0] - points[0][0])
289        for pidx in points[0][2]:
290            cache.deltas[pidx][0] = der
291
292        for i in range(1, len(cache.deltas)-1):
293            x1, x2, x3 = points[i-1][0], points[i][0], points[i+1][0]
294            sx = x1+x2+x3
295            div = 3*(x1**2+x2**2+x3**2) - sx**2
296            if div > 1e-6:
297                der = (3*(x1*points[i-1][1] + x2*points[i][1] + x3*points[i+1][1]) - sx * (points[i-1][1] + points[i][1] + points[i+1][1])) / div
298            for pidx in points[i][2]:
299                cache.deltas[pidx][0] = der
300
301    else:
302        for i in range(len(cache.deltas)-1):
303            der = (points[i][1] - points[i+1][1]) / (points[i][0] - points[i+1][0])
304            for pidx in points[i][2]:
305                cache.deltas[pidx][0] = der
306       
307   
308# PCA
309def pca(self,PCAthreshold=.8):
310    if not self.points:       
311        self.points = orange.ExampleTable(orange.Domain(self.contAttributes, self.data.domain.classVar), self.data).native(0)
312    #x = numpy.array(list(self.points))
313    x = numpy.array([[5.6359, 1.8749, 0.2587, 0.9281],[-9.0487, 7.8200, 0.6616, 0.7096],[-7.3699, -3.3882, 0.2115, 0.9430],[-2.1107, 6.4463, 0.6226, 0.1166],[-9.2594, -3.0736, 0.2573, 0.0401],[5.4889, 6.1418, 0.8147, 0.4558],[1.0844, -2.8861, 0.0941, 0.2960],[-8.0125, -7.9141, 0.0397, 0.5051],[-7.2509, -4.6146, 0.2589, 0.3370],[4.3882, 2.4636, 0.9816, 0.9167]])
314    xlen = len(x)
315    #print x
316    #print xlen
317    meanx = numpy.mean(x,0)
318    #print meanx
319    x = numpy.matrix([i-meanx for i in x])
320    #print x
321    Cov = (numpy.transpose(x)*x)/xlen
322    #print "Cov = ",Cov
323    evals,evecs = numpy.linalg.eig(Cov)
324    sum_evals = sum(list(evals))
325    eigv = [(evals[i]/sum_evals,i) for i in xrange(len(evals))]
326    eigv.sort()
327    eigv.reverse()
328    total = 0
329    for i in eigv:
330        total += i[0]
331        if total > PCAthreshold:
332            limit = i[1]
333            break
334    Ev = evecs[:,0:limit+1]
335    print eigv
336    print Ev
337    #pca_data = (Ev' * data')';
338    #pca_example = (Ev' * example')'
339
340
341def Dpca(self):
342    self.pca()
343    return
344    if not self.deltas:
345        self.deltas = [[None] * len(self.contAttributes) for x in xrange(len(self.data))]
346
347    dimensions = [d for d in self.dimensions if not self.deltas[0][d]]
348
349    if not dimensions:
350        return
351
352    if not self.points:       
353        self.points = orange.ExampleTable(orange.Domain(self.contAttributes, self.data.domain.classVar), self.data).native(0)
354    points = self.points
355    npoints = len(points)
356
357    if not self.tri:
358        print self.dimension
359        self.tri = self.triangulate(points)
360    tri = self.tri
361
362    if not self.stars:
363        self.stars = [star(x, tri) for x in xrange(npoints)]
364    S = self.stars
365
366    points = import85
367    nPoints = 100.0/len(points)
368
369    self.progressBarInit()
370    for x,(S,p) in enumerate(zip(self.stars,points)):
371        if S==[]:
372            self.deltas[x] = ['?' for i in dimensions]
373            continue
374        st = list(set(reduce(lambda x,y: x+y, S)))
375        lenst = len(st)
376        avgy = sum([points[i][-1] for i in st])/lenst
377        for di, d in enumerate(dimensions):
378            avgx = sum([points[i][di] for i in st])/lenst
379            sxx2 = sum([(points[i][di]-avgx)**2 for i in st])
380            if sxx2:
381                sxx = sum([(points[i][di]-avgx)*(points[i][-1]-avgy) for i in st])
382                b = sxx/sxx2
383                self.deltas[x][di] = b
384            else:
385                self.deltas[x][di] = '?'
386        self.progressBarSet(x*nPoints)
387
388    self.progressBarFinished()
389
390def star1D(i,p):
391    lenp = len(p)
392    if lenp==1:
393        return []
394    if i==0:
395        return [p[1]]
396    elif i==lenp-1:
397        return [p[lenp-2]]
398    else:
399        return [p[i-1],p[i+1]]
400
401#def qing1D(cache, dimensions, progressCallback = None, persistence=0.4):
402#    Dim = len(dimensions)
403#    print "dimension = ",Dim
404#    if not cache.points:       
405#        cache.points = orange.ExampleTable(orange.Domain(cache.contAttributes, cache.data.domain.classVar), cache.data).native(0)
406#    p = cache.points
407#    lenp = len(p)
408#    p.sort()
409#    # in > 1D we have to look for the points in the cylinder, project them to 1D and continue as before
410#    for d in dimensions:
411#        for x, (pts,deltas) in enumerate(zip(cache.points,cache.deltas)):
412#            px = pts[0:d]+pts[d+1:]
413#            px = px[:-1]
414#            print x,px
415#            if x>10: break
416#        print
417
418
419def qing1D(cache, dimensions, progressCallback = None, **args):
420    persistence = args.get("persistence", 0.4)
421    Dim = len(dimensions)
422    #print "dimension = ",Dim
423    if not cache.points:       
424        cache.points = orange.ExampleTable(orange.Domain(cache.contAttributes, cache.data.domain.classVar), cache.data).native(0)
425    p = cache.points
426    lenp = len(p)
427    p.sort()
428   
429    # in > 1D we have to look for the points in the cylinder, project them to 1D and continue as before
430    if Dim>1:
431        pass
432       
433    sosedi = [(p[i],star1D(i,p)) for i in xrange(lenp)]
434    mini = []
435    maxi = []
436    for i,s in enumerate(sosedi):
437        if s[0][1] < min([j[1] for j in s[1]]):
438            mini += [i]
439        elif s[0][1] > max([j[1] for j in s[1]]):
440            maxi += [i]
441           
442    #print "vsi minimumi: ",mini
443    #print "vsi maximumi: ",maxi
444    mini1 = deepcopy(mini)
445
446    for i in mini:
447        fv = p[i][1]
448        left = [j for j in maxi if j<i]
449        if left:
450            left = max(left)
451            left = (math.fabs(fv-p[left][1]),left)
452        right = [j for j in maxi if j>i]
453        if right:
454            right = min(right)
455            right = (math.fabs(fv-p[right][1]),right)
456        t = [(i,fv)]
457        ti = []
458        if left:
459            ti += [left]
460        if right:
461            ti += [right]
462        mti = min(ti)
463        if mti[0] <= persistence:
464            maxi.remove(mti[1])
465            mini1.remove(i)
466    #print "min. po krajsanju: ",mini1
467    #print "max. po krajsanju: ",maxi
468    ext = [(p[m],-1) for m in mini1]+[(p[m],1) for m in maxi]
469    ext.sort()
470    elist = zip(mini1,[-1]*len(mini1))+zip(maxi,[1]*len(maxi))
471    elist.sort()
472    ekstremi = []
473    if elist[0][1]==-1: # ce je najprej minimum
474        ekstrem = min([(p[k][-1],k,-1) for k in xrange(maxi[0])]) # poiscemo pravi min med tockami do prvega maxa
475    else:
476        ekstrem = max([(p[k][-1],k,1) for k in xrange(mini1[0])]) # sicer poiscemo prvi max med tockami do prvega mina
477    ekstremi += [ekstrem]
478    for i in xrange(1,len(elist)-1):
479        if elist[i][-1]==1:
480            ekstrem = max([(p[k][-1],k,1) for k in xrange(elist[i-1][0],elist[i+1][0]+1)])
481        else:
482            ekstrem = min([(p[k][-1],k,-1) for k in xrange(elist[i-1][0],elist[i+1][0]+1)])
483        ekstremi += [ekstrem]       
484
485    if elist[-1][1]==-1: # ce je zadnji ekstrem minimum
486        ekstrem = min([(p[k][-1],k,-1) for k in xrange(maxi[-1],lenp)]) # poiscemo pravi min med tockami do prvega maxa
487    else:
488        ekstrem = max([(p[k][-1],k,1) for k in xrange(mini1[-1],lenp)]) # sicer poiscemo prvi max med tockami do prvega mina
489    ekstremi += [ekstrem]
490    lookup = [k[1:] for k in ekstremi]
491    if progressCallback:
492        nPoints = 100.0/lenp
493    for x, deltas in enumerate(cache.deltas):
494        for d in dimensions:
495            #print x,p[x]
496            # ce je x ekstrem, ima odvod 0, oz. gledamo levo/desno limito, ce je na robu
497            if x in [j[0] for j in lookup]:
498                if x == lookup[0][0]:
499                    #deltas[d] = -lookup[0][1]
500                    deltas[d] = (p[lookup[1][0]][1] - p[x][1])/(p[lookup[1][0]][0] - p[x][0])
501                elif x == lookup[-1][0]:
502                    #deltas[d] = lookup[-1][1]
503                    deltas[d] = (p[lookup[-2][0]][1] - p[x][1])/(p[lookup[-2][0]][0] - p[x][0])
504                else:
505                    deltas[d] = 0.
506                continue
507            j=0
508            while lookup[j][0] < x:
509                j+=1
510            if lookup[j][1] == 1:
511                #deltas[d] = 1. #obrni * (xnz-xp[-1]) / dt
512                deltas[d] = (p[x][1] - p[lookup[j][0]][1])/(p[x][0] - p[lookup[j][0]][0])
513            else:
514                #deltas[d] = -1.
515                deltas[d] = (p[x][1] - p[lookup[j][0]][1])/(p[x][0] - p[lookup[j][0]][0])
516        if progressCallback:
517            progressCallback(x*nPoints)
518
519    #print [k[-1] for k in ekstremi]
520
521
522# regression in a tube
523def tubedRegression(cache, dimensions, progressCallback = None, **args):
524    if not cache.findNearest:
525        cache.findNearest = orange.FindNearestConstructor_BruteForce(cache.data, distanceConstructor=orange.ExamplesDistanceConstructor_Euclidean(), includeSame=True)
526       
527    if not cache.attrStat:
528        cache.attrStat = orange.DomainBasicAttrStat(cache.data)
529
530    normalizers = cache.findNearest.distance.normalizers
531
532    if progressCallback:
533        nExamples = len(cache.data)
534        nPoints = 100.0/nExamples/len(dimensions)
535
536    effNeighbours = len(cache.contAttributes) > 1 and cache.nNeighbours or len(cache.deltas)
537   
538    for di, d in enumerate(dimensions):
539        contIdx = cache.contIndices[d]
540
541        minV, maxV = cache.attrStat[contIdx].min, cache.attrStat[contIdx].max
542        if minV == maxV:
543            continue
544       
545        oldNormalizer = normalizers[cache.contIndices[d]]
546        normalizers[cache.contIndices[d]] = 0
547
548        for exi, ref_example in enumerate(cache.data):
549            if ref_example[contIdx].isSpecial():
550                cache.deltas[exi][d] = "?"
551                continue
552
553            ref_x = float(ref_example[contIdx])
554
555            Sx = Sy = Sxx = Syy = Sxy = n = 0.0
556
557            nn = cache.findNearest(ref_example, 0, True)
558            nn = [ex for ex in nn if not ex[contIdx].isSpecial()][:effNeighbours]
559            mx = [abs(ex[contIdx] - ref_x) for ex in nn]
560            if not mx:
561                cache.deltas[exi][d] = "?"
562                continue
563            if max(mx) < 1e-10:
564                kw = math.log(.001)
565            else:
566                kw = math.log(.001) / max(mx)**2
567            for ex in nn[:effNeighbours]:
568                ex_x = float(ex[contIdx])
569                ex_y = float(ex.getclass())
570                w = math.exp(kw*(ex_x-ref_x)**2)
571                Sx += w * ex_x
572                Sy += w * ex_y
573                Sxx += w * ex_x**2
574                Syy += w * ex_y**2
575                Sxy += w * ex_x * ex_y
576                n += w
577
578            div = n*Sxx-Sx**2
579            if div:# and i<40:
580                b = (Sxy*n - Sx*Sy) / div
581               
582#                div = Sx*Sy/n - Sxy
583#                if abs(div) < 1e-10:
584#                    cache.errors[exi][d] = 1
585#                else:
586#                    B = ((Syy - Sy**2/n) - (Sxx - Sx**2/n)) / 2 / div
587#
588#                    b_p = -B + math.sqrt(B**2+1)
589#                    a = Sy/n - b_p * Sx/n
590#                    error1 = 1/(1+b_p**2) * (Syy + a**2 + b_p**2*Sxx - 2*a*Sy + 2*a*b_p*Sx - 2*b_p*Sxy)
591#
592#                    b_2 = -B - math.sqrt(B**2+1)
593#                    a = Sy/n - b_p * Sx/n
594#                    error2 = 1/(1+b_p**2) * (Syy + a**2 + b_p**2*Sxx - 2*a*Sy + 2*a*b_p*Sx - 2*b_p*Sxy)
595#                   
596#                    if error1 < error2 and error1 >= 0:
597#                        cache.errors[exi][d] = error1
598#                    elif error2 >= 0:
599#                        cache.errors[exi][d] = error2
600#                    else:
601#                        cache.errors[exi][d] = 42
602#                        print error1, error2
603                           
604                a = (Sy - b*Sx)/n
605                err = (n * a**2 + b**2 * Sxx + Syy + 2*a*b*Sx - 2*a*Sy - 2*b*Sxy)
606                tot = Syy - Sy**2/n
607                mod = tot - err
608                merr = err/(n-2)
609                if merr < 1e-10:
610                    F = 0
611                    Fprob = 1
612                else:
613                    F = mod/merr
614                    Fprob = statc.fprob(F, 1, int(n-2))
615                cache.errors[exi][d] = Fprob
616#                        print "%.4f" % Fprob,
617                #print ("%.3f\t" + "%.0f\t"*6 + "%f\t%f") % (w, ref_x, ex_x, n, a, b, merr, F, Fprob)
618                cache.deltas[exi][d] = b
619            else:
620                cache.deltas[exi][d] = "?"
621
622            if progressCallback:
623                progressCallback((nExamples*di+exi)*nPoints)
624
625        normalizers[cache.contIndices[d]] = oldNormalizer
626
627
628def createClassVar(attributes, MQCNotation = False):
629    import orngMisc
630    if MQCNotation:
631        return orange.EnumVariable("Q", values = ["%s(%s)" % ("".join(["+-"[x] for x in v if x<2]), ", ".join([attr for attr,x in zip(attributes, v) if x<2])) for v in orngMisc.LimitedCounter([3]*len(attributes))])
632    else:
633        return orange.EnumVariable("Q", values = ["Q(%s)" % ", ".join(["+-"[x]+attr for attr, x in zip(attributes, v) if x<2]) for v in orngMisc.LimitedCounter([3]*len(attributes))])
634
635   
636def createQTable(cache, data, dimensions, outputAttr = -1, threshold = 0, MQCNotation = False, derivativeAsMeta = False, differencesAsMeta = False, correlationsAsMeta = False, originalAsMeta = False):
637    nDimensions = len(dimensions)
638   
639    needQ = outputAttr < 0 or derivativeAsMeta
640    if needQ:
641        qVar = createClassVar([cache.attributes[i][0] for i in dimensions], MQCNotation)
642       
643    if outputAttr >= 0:
644        classVar = orange.FloatVariable("df/d"+cache.attributes[outputAttr][0])
645    else:
646        classVar = qVar
647           
648    dom = orange.Domain(data.domain.attributes, classVar)
649    dom.addmetas(data.domain.getmetas())
650    setattr(dom, "constraintAttributes", [cache.contAttributes[i] for i in dimensions])
651
652    if derivativeAsMeta:
653        derivativeID = orange.newmetaid()
654        dom.addmeta(derivativeID, qVar)
655    else:
656        derivativeID = 0
657           
658    metaIDs = []       
659    if differencesAsMeta:
660        for dim in dimensions:
661            metaVar = orange.FloatVariable("df/d"+cache.attributes[dim][0])
662            metaID = orange.newmetaid()
663            dom.addmeta(metaID, metaVar)
664            metaIDs.append(metaID)
665
666    corMetaIDs = []
667    if correlationsAsMeta:
668        for dim in dimensions:
669            metaVar = orange.FloatVariable("corr(%s)" % cache.attributes[dim][0])
670            metaID = orange.newmetaid()
671            dom.addmeta(metaID, metaVar)
672            corMetaIDs.append(metaID)
673        metaVar = orange.FloatVariable("corr")
674        metaID = orange.newmetaid()
675        dom.addmeta(metaID, metaVar)
676        corMetaIDs.append(metaID)
677
678    if originalAsMeta:
679        originalID = orange.newmetaid()
680        dom.addmeta(originalID, data.domain.classVar)
681    else:
682        originalID = 0
683
684
685    paded = orange.ExampleTable(dom, data)
686
687    for i, (pad, alldeltas) in enumerate(zip(paded, cache.deltas)):
688        deltas = [alldeltas[d] for d in dimensions]
689       
690        if needQ:
691            qs = "".join([(delta > threshold and "0") or (delta < -threshold and "1") or (delta == "?" and "?") or "2" for delta in deltas])
692            q = ("?" in qs and "?") or int(qs, 3)
693               
694        if outputAttr >= 0:
695            pad.setclass(alldeltas[outputAttr])
696        else:
697            pad.setclass(q)
698
699        if derivativeAsMeta:
700            pad.setmeta(derivativeID, q)
701
702        if differencesAsMeta:
703            for a in zip(metaIDs, deltas):
704                pad.setmeta(*a)
705
706        if correlationsAsMeta:
707            if hasattr(cache, "errors"):
708                maxerr = -1e20
709                for id, val in zip(corMetaIDs, [cache.errors[i][d] for d in dimensions]):
710                    if val == None:
711                        pad.setmeta(id, "?")
712                    else:
713                        pad.setmeta(id, val)
714                        maxerr = max(maxerr, val)
715                pad.setmeta(corMetaIDs[-1], maxerr)
716            else:
717                minder = 0
718                for id, val in zip(corMetaIDs[:-1], deltas):
719                    if type(val) == str:
720                        pad.setmeta(id, "?")
721                    else:
722                        pad.setmeta(id, abs(val))
723                        minder = min(minder, abs(val))
724                pad.setmeta(corMetaIDs[-1], minder)
725               
726    return paded, derivativeID, metaIDs, corMetaIDs, originalID
727
728
729def pade(data, attributes = None, method = tubedRegression, outputAttr = -1, threshold = 0, MQCNotation = False, derivativeAsMeta = False, differencesAsMeta = False, correlationsAsMeta = False, originalAsMeta = False):
730    cache = makeBasicCache(data)
731    cache.deltas = [[None] * len(cache.contAttributes) for x in xrange(len(data))]
732    if method == tubedRegression:
733        cache.errors = [[None] * len(cache.contAttributes) for x in xrange(len(data))]
734
735    cache.nNeighbours = 30
736
737    if not attributes:
738        attributes = range(len(cache.contAttributes))
739       
740    dimensions = [data.domain.index(attr) for attr in attributes]
741
742    if outputAttr != -1:
743        outputAttr = data.domain.index(outputAttr)
744        if outputAttr not in dimensions:
745            dimensions.append(outputAttr)
746
747    method(cache, dimensions)
748    return createQTable(cache, data, dimensions, outputAttr, threshold, MQCNotation, derivativeAsMeta, differencesAsMeta, correlationsAsMeta, originalAsMeta)
749
750
751### Quin-like measurement of quality
752#
753# We consider these measures completely inappropriate for various reasons.
754# We nevertheless implemented them for the sake of comparison.
755
756# puts examples from 'data' into the corresponding leaves of the subtree
757def splitDataOntoNodes(node, data):
758    if not node.branches:
759        node.setattr("leafExamples", data)
760    else:
761        goeswhere = [(ex, int(node.branchSelector(ex))) for ex in data]
762        for bi, branch in enumerate(node.branches):
763            splitDataOntoNodes(branch, [ex for ex, ii in goeswhere if ii == bi])
764
765
766# checks whether 'ex1' and 'ex2' match the given constraint
767# 'reversedAttributes' is a reversed list of the arguments of the QMC
768# 'constraint' is the index of the class
769# 'clsID' is the id of the meta attribute with the original class value
770#
771# Result: -2 classes or all attribute values are unknown or same for both
772#            examples. This not considered a model's fault and doesn't count
773#         -1 Ambiguity: one attribute value goes along the constraint and
774#            another goes the opposite. This is ambiguity due to the model.
775#          0 False prediction
776#          1 Correct prediction
777#
778# Note: Quin does not distinguish between ambiguity due to the data and
779#   due to the model. We understand the ambiguity in the data as
780#   unavoidable and don't count the corresponding example pairs, while
781#   ambiguity due to the model is model's fault and could even be
782#   counted as misclassification
783
784def checkMatch(ex1, ex2, reversedAttributes, constraint, clsID):
785    cls1, cls2 = ex1[clsID], ex2[clsID]
786    if cls1.isSpecial() or cls2.isSpecial():
787        return -2  # unknown classes - useless example
788    clsc = cmp(cls1, cls2)
789    if not clsc:
790        return -2  # ambiguity due to same class
791   
792    cs = 0
793    for attr in reversedAttributes:
794        constraint, tc = constraint/3, constraint%3
795        if tc==2:
796            continue
797        v1, v2 = ex1[attr], ex2[attr]
798        if v1.isSpecial() or v2.isSpecial():
799            return -2   # unknowns - useless example
800        if tc:
801            c = -cmp(v1, v2)
802        else:
803            c = cmp(v1, v2)
804        if not cs:
805            cs = c
806        elif c != cs:
807            return -1   # ambiguity due to example pair
808    if not cs:
809        return -2       # ambiguity due to same example values
810
811    return cs == clsc
812
813
814def computeAmbiguityAccuracyNode(node, reversedAttributes, clsID):
815    samb = sacc = spairs = 0.0
816    if node.branches:
817        for branch in node.branches:
818            if branch:
819                amb, acc, pairs = computeAmbiguityAccuracyNode(branch, reversedAttributes, clsID)
820                samb += amb
821                sacc += acc
822                spairs += pairs
823    else:
824        constraint = int(node.nodeClassifier.defaultVal)
825        for i, ex1 in enumerate(node.leafExamples):
826            for j in range(i):
827                ma = checkMatch(ex1, node.leafExamples[j], reversedAttributes, constraint, clsID)
828                if ma == -2:
829                    continue
830                if ma == -1:
831                    samb += 1
832                elif ma:
833                    sacc += 1
834                else:
835                    ex2 = node.leafExamples[j]
836                spairs += 1
837    return samb, sacc, spairs
838
839
840def computeAmbiguityAccuracy(tree, data, clsID):
841    splitDataOntoNodes(tree.tree, data)
842    l = tree.domain.constraintAttributes[:]
843    l.reverse()
844    amb, acc, pairs = computeAmbiguityAccuracyNode(tree.tree, l, clsID)
845    return amb/pairs, acc/(pairs-amb)
846
847
848def CVByNodes(data, dimensions = None, method = None, **dic):
849    import orngTree
850    cv = orange.MakeRandomIndicesCV(data, 10)
851    for fold in range(10):
852        train = data.select(cv, fold, negate=1)
853        test = data.select(cv, fold)
854        pa, qid, did, cid = pade(train, dimensions, method, originalAsMeta=True, **dic)
855        tree = orngTree.TreeLearner(pa, maxDepth=4)
856
857        mb, cc = computeAmbiguityAccuracy(tree, test, -1)
858        amb += mb
859        acc += cc
860    return amb/10, acc/10
861
862
863### Better measures of quality (as in better-than-Quin)
864#
865def checkDirectionAccuracyForPair(model, reference, direction, clsID, reversedAttributes):
866    constraint = int(model(reference))
867    prediction = 0
868
869    for attr in reversedAttributes:
870        constraint, tc = constraint/3, constraint%3
871        if tc==2:
872            continue
873
874        v1, v2 = reference[attr], direction[attr]
875        if v1.isSpecial() or v2.isSpecial():
876            continue
877
878        if tc:
879            c = -cmp(v1, v2)
880        else:
881            c = cmp(v1, v2)
882
883        if not prediction:
884            prediction = c
885        elif prediction != c:
886            return -1
887
888    if not prediction:
889        return -3
890   
891    return c == cmp(reference.getclass(), direction.getclass())
892
893
894def computeDirectionAccuracyForPairs(model, data, meter, weightK, clsID, nTests = 0):
895    nTests = nTests or 10*len(data)
896
897    reversedAttrs = model.domain.constraintAttributes[:]
898    reversedAttrs.reverse()
899
900    actTests = acc = amb = unre = 0
901    for i in range(10*len(data)):
902        distance = 0
903        while not distance:
904            ref, dir = data.randomexample(), data.randomexample()
905            distance = meter(ref, dir)
906        weight = math.exp(weightK * distance**2)
907        actTests += weight
908        diracc = checkDirectionAccuracyForPair(model, ref, dir, -1, reversedAttrs)
909        if diracc == -1:
910            amb += weight
911        elif diracc == -3:
912            unre += weight
913        elif diracc:
914            acc += weight
915
916    return acc/actTests, amb/actTests, unre/actTests
917       
918
919def CVByPairs(data, dimensions = None, method = None, **dic):
920    import orngTree
921    cv = orange.MakeRandomIndicesCV(data, 10)
922    meter = orange.ExamplesDistanceConstructor_Euclidean(data)
923
924    maxDist = 0
925    for i in range(100):
926        maxDist = max(maxDist, meter(data.randomexample(), data.randomexample()))
927    weightK = 10.0 / maxDist
928
929    acc = amb = unre = 0
930    for fold in range(10):
931        train = data.select(cv, fold, negate=1)
932        test = data.select(cv, fold)
933        pa, qid, did, cid = pade(train, dimensions, method, originalAsMeta=True, **dic)
934        tree = orngTree.TreeLearner(pa, maxDepth=4)
935
936        tacc, tamb, tunre = computeDirectionAccuracyForPairs(tree, data, meter, weightK, -1)
937        acc += tacc
938        amb += tamb
939        unre += tunre
940       
941    return acc/10, amb/10, unre/10
942
943
944import re
945sidx = {"+": "0", "-": "1", None: "2"}
946numb = r"[+-]?\d*(\.\d*)?(e[+-]\d+)?"
947constraint = r"M(?P<signs>[+-](,[+-])*)\((?P<attrs>[^ ,)]*(,[^ ,)]*)*)\)\s*\(Cov=(?P<cov>(\*\*)|("+numb+r"%="+numb+r"))\(\?(?P<amb>"+numb+r")%\)/(?P<xmpls>"+numb+r")\)"
948re_constraint = re.compile(constraint)
949re_node = re.compile(r"\s*(?P<splitattr>[^ ]*)\s*<=\s*(?P<thresh>[\-0-9.]*)\s*\[\s*"+constraint+r"\]")
950
951def readQuinTreeRec(fle, domain):
952    node = orange.TreeNode()
953    line = fle.readline()
954
955    match = re_node.search(line)
956    if match:
957        splitattr, thresh = match.group("splitattr", "thresh")
958        node.branchSelector = orange.Classifier(lambda ex,rw, attr=domain[splitattr], thr=float(thresh): ex[attr] > thr)
959        node.branchDescriptions = ["<= "+thresh, "> "+thresh]
960        node.branches = [readQuinTreeRec(fle, domain), readQuinTreeRec(fle, domain)]
961        node.branchSizes = [node.branches[0].nExamples, node.branches[1].nExamples]
962    else:
963        match = re_constraint.search(line)
964        if not match:
965            raise "Cannot read line '%s'" % line
966
967    attrs, signs, cov, amb, xmpls = match.group("attrs", "signs", "cov", "amb", "xmpls")
968
969    mqc = dict(zip(attrs.split(","), signs.split(",")))
970    node.nodeClassifier = orange.DefaultClassifier(defaultVal = int("".join([sidx[mqc.get(attr.name, None)] for attr in domain.attributes]), 3))
971    node.setattr("coverage", cov == "**" and -1 or float(cov[:cov.index("%")]))
972    node.setattr("ambiguity", float(amb))
973    node.setattr("nExamples", float(xmpls))
974    return node
975
976def readQuinTree(fname, domain):
977    qVar = createClassVar([attr.name for attr in domain.attributes])
978    domain.setattr("constraintAttributes", domain.attributes)
979    return orange.TreeClassifier(domain = domain, classVar = domain.classVar, tree = readQuinTreeRec(open(fname), domain))
980
Note: See TracBrowser for help on using the repository browser.