source: orange/Orange/projection/linear.py @ 10477:82d4a84af229

Revision 10477:82d4a84af229, 76.9 KB checked in by Matija Polajnar <matija.polajnar@…>, 2 years ago (diff)

Bugfix: missing import in linear projection.

Line 
1#TODO: eliminate create_pls_projection (transform into a class)
2#TODO: Projector as a preprocessor
3
4import Orange
5from Orange import orangeom
6import math
7import random
8import numpy
9
10from Orange.preprocess.scaling import ScaleLinProjData
11from Orange.orng import orngVisFuncts as visfuncts
12from Orange.misc import deprecated_keywords
13from Orange.misc import deprecated_members
14
15try:
16    import numpy.ma as MA
17except:
18    import numpy.core.ma as MA
19
20#implementation
21FAST_IMPLEMENTATION = 0
22SLOW_IMPLEMENTATION = 1
23LDA_IMPLEMENTATION = 2
24
25LAW_LINEAR = 0
26LAW_SQUARE = 1
27LAW_GAUSSIAN = 2
28LAW_KNN = 3
29LAW_LINEAR_PLUS = 4
30
31DR_PCA = 0
32DR_SPCA = 1
33DR_PLS = 2
34
35def normalize(x):
36    return x / numpy.linalg.norm(x)
37
38def center(matrix):
39    '''centers all variables, i.e. subtracts averages in colomns
40    and divides them by their standard deviations'''
41    n,m = numpy.shape(matrix)
42    return (matrix - numpy.multiply(matrix.mean(axis = 0),
43                                    numpy.ones((n,m))))/numpy.std(matrix,
44                                                                  axis = 0)
45
46
47class FreeViz:
48    """
49    Contains an easy-to-use interface to the core of the method, which is
50    written in C++. Differs from other linear projection optimizers in that it itself can store the data
51    to make iterative optimization and visualization possible. It can, however, still be used as any other
52    projection optimizer by calling (:obj:`~Orange.projection.linear.FreeViz.__call__`) it.
53   
54    .. attribute:: attract_g
55   
56    Coefficient for the attractive forces. By increasing or decreasing the ratio
57    between :obj:`attract_g` and :obj:`repel_g`, you can make one kind of the
58    forces stronger. Default: 1.
59   
60    .. attribute:: repel_g
61   
62    Coefficient for the repulsive forces. By increasing or decreasing the ratio
63    between :obj:`attract_g` and :obj:`repel_g`, you can make one kind of the
64    forces stronger. Default: 1.
65   
66    .. attribute:: force_balancing
67   
68    If set (default is False), the forces are balanced so that the total sum of
69    the attractive equals the total of repulsive, before they are multiplied by
70    the above factors. (By our experience, this gives bad results so you may
71    want to leave this alone.)
72   
73    .. attribute:: law
74   
75    Can be LAW_LINEAR, LAW_SQUARE, LAW_GAUSSIAN, LAW_KNN or LAW_LINEAR_PLUS.
76    Default is LAW_LINEAR, which means that the attractive forces increase
77    linearly by the distance and the repulsive forces are inversely
78    proportional to the distance. LAW_SQUARE would make them rise or fall with
79    the square of the distance, LAW_GAUSSIAN is based on a kind of
80    log-likelihood estimation, LAW_KNN tries to directly optimize the
81    classification accuracy of the kNN classifier in the projection space, and
82    in LAW_LINEAR_PLUS both forces rise with the square of the distance,
83    yielding a method that is somewhat similar to PCA. We found the first law
84    perform the best, with the second to not far behind.
85   
86    .. attribute:: force_sigma
87   
88    The sigma to be used in LAW_GAUSSIAN and LAW_KNN.
89   
90    .. attribute:: mirror_symmetry
91   
92    If enabled, it keeps the projection of the second attribute on the upper
93    side of the graph (the first is always on the right-hand x-axis). This is
94    useful when comparing whether two projections are the same, but has no
95    effect on the projection's clarity or its classification accuracy.
96
97    There are some more, undescribed, methods of a more internal nature.
98   
99    """
100   
101    def __init__(self, graph = None):
102        if not graph:
103            graph = ScaleLinProjData()
104        self.graph = graph
105
106        self.implementation = 0
107        self.attract_g = 1.0
108        self.repel_g = 1.0
109        self.law = LAW_LINEAR
110        self.restrain = 0
111        self.force_balancing = 0
112        self.force_sigma = 1.0
113        self.mirror_symmetry = 1
114        self.use_generalized_eigenvectors = 1
115
116        # s2n heuristics parameters
117        self.steps_before_update = 10
118        self.s2n_spread = 5
119        self.s2n_place_attributes = 50
120        self.s2n_mix_data = None
121        self.auto_set_parameters = 1
122        self.class_permutation_list = None
123        self.attrs_num = [5, 10, 20, 30, 50, 70, 100, 150, 200, 300, 500, 750,
124                          1000]
125
126    def __call__(self, dataset=None):
127        """
128        Perform FreeViz optimization on the dataset, if given, and return a resulting
129        linear :class:`~Orange.projection.linear.Projector`. If no dataset is given,
130        the projection currently stored within the FreeViz object is returned as
131        a :class:`~Orange.projection.linear.Projector`.
132
133        :param dataset: input data set.
134        :type dataset: :class:`Orange.data.Table`
135
136        :rtype: :class:`~Orange.projection.linear.Projector`
137        """
138        if dataset:
139            self.graph.setData(dataset)
140            self.show_all_attributes()
141
142            self.radial_anchors()
143            self.optimize_separation()
144
145            X = dataset.to_numpy_MA("a")[0]
146            Xm = numpy.mean(X, axis=0)
147            Xd = X - Xm
148            stdev = numpy.std(Xd, axis=0)
149        else:
150            Xm = numpy.zeros(len(self.graph.anchor_data))
151            stdev = None
152
153        graph = self.graph
154
155        U = numpy.array([val[:2] for val in self.graph.anchor_data]).T
156
157        domain = graph.data_domain
158        if len(domain) > len(self.graph.anchor_data):
159            domain = Orange.data.Domain([graph.data_domain[a]
160                                         for _,_,a in self.graph.anchor_data],
161                graph.data_domain.class_var)
162
163        return Orange.projection.linear.Projector(input_domain = domain,
164            mean = Xm,
165            stdev = stdev,
166            standardize = False,
167            projection = U)
168
169
170    def clear_data(self):
171        self.s2n_mix_data = None
172        self.class_permutation_list = None
173       
174    clearData = clear_data
175
176    def set_statusbar_text(self, *args):
177        pass
178   
179    setStatusBarText = set_statusbar_text
180
181    def show_all_attributes(self):
182        self.graph.anchor_data = [(0,0, a.name)
183                                 for a in self.graph.data_domain.attributes]
184        self.radial_anchors()
185       
186    showAllAttributes = show_all_attributes
187
188    def get_shown_attribute_list(self):
189        return [anchor[2] for anchor in self.graph.anchor_data]
190
191    getShownAttributeList = get_shown_attribute_list
192
193    def radial_anchors(self):
194        """
195        Reset the projection so that the anchors (projections of attributes)
196        are placed evenly around the circle.
197       
198        """
199        attr_list = self.get_shown_attribute_list()
200        if not attr_list:
201            return
202        if hasattr(self, "parentName") and "3d" in self.parentName.lower():
203            self.graph.anchor_data = self.graph.create_anchors(len(attr_list), attr_list)
204            return
205        phi = 2*math.pi/len(attr_list)
206        self.graph.anchor_data = [(math.cos(i*phi), math.sin(i*phi), a)
207                                 for i, a in enumerate(attr_list)]
208
209    radialAnchors = radial_anchors
210
211    def random_anchors(self):
212        """
213        Set the projection to a random one.
214       
215        """
216        if not self.graph.have_data:
217            return
218        attr_list = self.get_shown_attribute_list()
219        if not attr_list:
220            return
221        if "3d" in self.parentName.lower():
222            if self.restrain == 0:
223                def ranch(i, label):
224                    r = 0.3+0.7*random.random()
225                    phi = 2*math.pi*random.random()
226                    theta = math.pi*random.random()
227                    return (r*math.sin(theta)*math.cos(phi),
228                            r*math.sin(theta)*math.sin(phi),
229                            r*math.cos(theta),
230                            label)
231            elif self.restrain == 1:
232                def ranch(i, label):
233                    phi = 2*math.pi*random.random()
234                    theta = math.pi*random.random()
235                    r = 1.
236                    return (r*math.sin(theta)*math.cos(phi),
237                            r*math.sin(theta)*math.sin(phi),
238                            r*math.cos(theta),
239                            label)
240            else:
241                self.graph.anchor_data = self.graph.create_anchors(len(attr_list), attr_list)
242                def ranch(i, label):
243                    r = 0.3+0.7*random.random()
244                    return (r*self.graph.anchor_data[i][0],
245                            r*self.graph.anchor_data[i][1],
246                            r*self.graph.anchor_data[i][2],
247                            label)
248
249            anchors = [ranch(*a) for a in enumerate(attr_list)]
250
251            if not self.restrain == 1:
252                maxdist = math.sqrt(max([x[0]**2+x[1]**2+x[2]**2 for x in anchors]))
253                anchors = [(x[0]/maxdist, x[1]/maxdist, x[2]/maxdist, x[3]) for x in anchors]
254
255            self.graph.anchor_data = anchors
256            return
257
258        if self.restrain == 0:
259            def ranch(i, label):
260                r = 0.3+0.7*random.random()
261                phi = 2*math.pi*random.random()
262                return (r*math.cos(phi), r*math.sin(phi), label)
263
264        elif self.restrain == 1:
265            def ranch(i, label):
266                phi = 2*math.pi*random.random()
267                return (math.cos(phi), math.sin(phi), label)
268
269        else:
270            def ranch(i, label):
271                r = 0.3+0.7*random.random()
272                phi = 2*math.pi * i / max(1, len(attr_list))
273                return (r*math.cos(phi), r*math.sin(phi), label)
274
275        anchors = [ranch(*a) for a in enumerate(attr_list)]
276
277        if not self.restrain == 1:
278            maxdist = math.sqrt(max([x[0]**2+x[1]**2 for x in anchors]))
279            anchors = [(x[0]/maxdist, x[1]/maxdist, x[2]) for x in anchors]
280
281        if not self.restrain == 2 and self.mirror_symmetry:
282            #### Need to rotate and mirror here
283            pass
284
285        self.graph.anchor_data = anchors
286
287    randomAnchors = random_anchors
288
289    @deprecated_keywords({"singleStep": "single_step"})
290    def optimize_separation(self, steps = 10, single_step = False, distances=None):
291        """
292        Optimize the class separation. If you did not change any of the settings
293        which are not documented above, it will call a fast C++ routine which
294        will make :obj:`steps` optimization steps at a time, after which the
295        graph (if one is given) is updated. If :obj:`single_step` is True, it
296        will do that only once,
297        otherwise it calls it on and on, and compares the current positions of
298        the anchors with those 50 calls ago. If no anchor moved for more than
299        1e-3, it stops. In Orange Canvas the optimization is also stopped if
300        someone outside (namely, the stop button) manages to set the FreeViz's
301        flag attribute
302        :obj:`Orange.projection.linear.FreeViz.cancel_optimization`.
303        """
304        # check if we have data and a discrete class
305        if (not self.graph.have_data or len(self.graph.raw_data) == 0
306            or not (self.graph.data_has_class or distances)):
307            return
308        ai = self.graph.attribute_name_index
309        attr_indices = [ai[label] for label in self.get_shown_attribute_list()]
310        if not attr_indices: return
311
312        if self.implementation == FAST_IMPLEMENTATION and not hasattr(self, '_use_3D'): # TODO
313            return self.optimize_fast_separation(steps, single_step, distances)
314
315        if self.__class__ != FreeViz: from PyQt4.QtGui import qApp
316        if single_step: steps = 1
317        if self.implementation == SLOW_IMPLEMENTATION:
318            impl = self.optimize_slow_separation
319        elif self.implementation == LDA_IMPLEMENTATION:
320            impl = self.optimize_lda_separation
321        xanchors = None
322        yanchors = None
323        zanchors = None
324
325        if hasattr(self, '_use_3D'):
326            if self.implementation == SLOW_IMPLEMENTATION:
327                impl = self.optimize_slow_separation_3D
328            elif self.implementation == LDA_IMPLEMENTATION:
329                impl = self.optimize_lda_separation_3D
330            else:
331                print('Unimplemented method!')
332                return
333
334            for c in range((single_step and 1) or 50):
335                for i in range(steps):
336                    if self.__class__ != FreeViz and self.cancel_optimization == 1:
337                        return
338                    self.graph.anchor_data, (xanchors, yanchors, zanchors) = impl(attr_indices,
339                                                                                 self.graph.anchor_data,
340                                                                                 xanchors,
341                                                                                 yanchors,
342                                                                                 zanchors)
343                if self.__class__ != FreeViz: qApp.processEvents()
344                if hasattr(self.graph, "updateGraph"): self.graph.updateData()
345        else:
346            for c in range((single_step and 1) or 50):
347                for i in range(steps):
348                    if self.__class__ != FreeViz and self.cancel_optimization == 1:
349                        return
350                    self.graph.anchor_data, (xanchors, yanchors) = impl(attr_indices,
351                                                                       self.graph.anchor_data,
352                                                                       xanchors,
353                                                                       yanchors)
354                if self.__class__ != FreeViz: qApp.processEvents()
355                if hasattr(self.graph, "updateGraph"): self.graph.updateData()
356
357    optimizeSeparation = optimize_separation
358
359    @deprecated_keywords({"singleStep": "single_step"})
360    def optimize_fast_separation(self, steps = 10, single_step = False, distances=None):
361        optimizer = [orangeom.optimizeAnchors, orangeom.optimizeAnchorsRadial,
362                     orangeom.optimizeAnchorsR][self.restrain]
363        ai = self.graph.attribute_name_index
364        attr_indices = [ai[label] for label in self.get_shown_attribute_list()]
365        if not attr_indices: return
366
367        # repeat until less than 1% energy decrease in 5 consecutive iterations*steps steps
368        positions = [numpy.array([x[:2] for x in self.graph.anchor_data])]
369        needed_steps = 0
370
371        valid_data = self.graph.get_valid_list(attr_indices)
372        n_valid = sum(valid_data) 
373        if not n_valid:
374            return 0
375
376        data = numpy.compress(valid_data, self.graph.no_jittering_scaled_data,
377                              axis=1)
378        data = numpy.transpose(data).tolist()
379        if self.__class__ != FreeViz: from PyQt4.QtGui import qApp
380
381        if distances:
382            if n_valid != len(valid_data):
383                classes = Orange.misc.SymMatrix(n_valid)
384                r = 0
385                for ro, vr in enumerate(valid_data):
386                    if not vr:
387                        continue
388                    c = 0
389                    for co, vr in enumerate(valid_data):
390                        if vr:
391                            classes[r, c] = distances[ro, co]
392                            c += 1
393                    r += 1 
394            else:
395                classes = distances
396        else:
397            classes = numpy.compress(valid_data,
398                                     self.graph.original_data[self.graph.data_class_index]).tolist()
399        while 1:
400            self.graph.anchor_data = optimizer(data, classes,
401                                              self.graph.anchor_data,
402                                              attr_indices,
403                                              attractG = self.attract_g,
404                                              repelG = self.repel_g,
405                                              law = self.law,
406                                              sigma2 = self.force_sigma,
407                                              dynamicBalancing = self.force_balancing,
408                                              steps = steps,
409                                              normalizeExamples = self.graph.normalize_examples,
410                                              contClass = 2 if distances
411                                              else self.graph.data_has_continuous_class,
412                                              mirrorSymmetry = self.mirror_symmetry)
413            needed_steps += steps
414
415            if self.__class__ != FreeViz:
416                qApp.processEvents()
417
418            if hasattr(self.graph, "updateData"):
419                self.graph.potentials_emp = None
420                self.graph.updateData()
421
422            positions = positions[-49:]+[numpy.array([x[:2] for x
423                                                      in self.graph.anchor_data])]
424            if len(positions)==50:
425                m = max(numpy.sum((positions[0]-positions[49])**2), 0)
426                if m < 1e-3: break
427            if single_step or (self.__class__ != FreeViz
428                               and self.cancel_optimization):
429                break
430        return needed_steps
431
432    optimize_FAST_Separation = optimize_fast_separation
433
434    @deprecated_keywords({"attrIndices": "attr_indices",
435                          "anchorData": "anchor_data",
436                          "XAnchors": "xanchors",
437                          "YAnchors": "yanchors"})
438    def optimize_lda_separation(self, attr_indices, anchor_data, xanchors = None, yanchors = None):
439        if (not self.graph.have_data or len(self.graph.raw_data) == 0
440            or not self.graph.data_has_discrete_class): 
441            return anchor_data, (xanchors, yanchors)
442        class_count = len(self.graph.data_domain.classVar.values)
443        valid_data = self.graph.get_valid_list(attr_indices)
444        selected_data = numpy.compress(valid_data,
445                                       numpy.take(self.graph.no_jittering_scaled_data,
446                                                  attr_indices, axis = 0),
447                                       axis = 1)
448
449        if xanchors == None:
450            xanchors = numpy.array([a[0] for a in anchor_data], numpy.float)
451        if yanchors == None:
452            yanchors = numpy.array([a[1] for a in anchor_data], numpy.float)
453
454        trans_proj_data = self.graph.create_projection_as_numeric_array(attr_indices,
455                                                                    validData = valid_data,
456                                                                    xanchors = xanchors,
457                                                                    yanchors = yanchors,
458                                                                    scaleFactor = self.graph.scale_factor,
459                                                                    normalize = self.graph.normalize_examples,
460                                                                    useAnchorData = 1)
461        if trans_proj_data == None:
462            return anchor_data, (xanchors, yanchors)
463
464        proj_data = numpy.transpose(trans_proj_data)
465        x_positions, y_positions, classData = (proj_data[0], proj_data[1],
466                                               proj_data[2])
467
468        averages = []
469        for i in range(class_count):
470            ind = classData == i
471            xpos = numpy.compress(ind, x_positions)
472            ypos = numpy.compress(ind, y_positions)
473            xave = numpy.sum(xpos)/len(xpos)
474            yave = numpy.sum(ypos)/len(ypos)
475            averages.append((xave, yave))
476
477        # compute the positions of all the points. we will try to move all points so that the center will be in the (0,0)
478        x_center_vector = -numpy.sum(x_positions) / len(x_positions)
479        y_center_vector = -numpy.sum(y_positions) / len(y_positions)
480        center_vector_length = math.sqrt(x_center_vector*x_center_vector +
481                                         y_center_vector*y_center_vector)
482
483        mean_destination_vectors = []
484
485        for i in range(class_count):
486            xdir = 0.0; ydir = 0.0; rs = 0.0
487            for j in range(class_count):
488                if i==j: continue
489                r = math.sqrt((averages[i][0] - averages[j][0])**2 +
490                              (averages[i][1] - averages[j][1])**2)
491                if r == 0.0:
492                    xdir += math.cos((i/float(class_count))*2*math.pi)
493                    ydir += math.sin((i/float(class_count))*2*math.pi)
494                    r = 0.0001
495                else:
496                    xdir += (1/r**3) * ((averages[i][0] - averages[j][0]))
497                    ydir += (1/r**3) * ((averages[i][1] - averages[j][1]))
498                #rs += 1/r
499            #actualDirAmpl = math.sqrt(xDir**2 + yDir**2)
500            #s = abs(xDir)+abs(yDir)
501            #xDir = rs * (xDir/s)
502            #yDir = rs * (yDir/s)
503            mean_destination_vectors.append((xdir, ydir))
504
505
506        maxlength = math.sqrt(max([x**2 + y**2 for (x,y)
507                                   in mean_destination_vectors]))
508        mean_destination_vectors = [(x/(2*maxlength), y/(2*maxlength)) for (x,y)
509                                    in mean_destination_vectors]     # normalize destination vectors to some normal values
510        mean_destination_vectors = [(mean_destination_vectors[i][0]+averages[i][0],
511                                     mean_destination_vectors[i][1]+averages[i][1])
512                                    for i in range(len(mean_destination_vectors))]    # add destination vectors to the class averages
513        #mean_destination_vectors = [(x + x_center_vector/5, y + y_center_vector/5) for (x,y) in mean_destination_vectors]   # center mean values
514        mean_destination_vectors = [(x + x_center_vector, y + y_center_vector)
515                                    for (x,y) in mean_destination_vectors]   # center mean values
516
517        fxs = numpy.zeros(len(x_positions), numpy.float)        # forces
518        fys = numpy.zeros(len(x_positions), numpy.float)
519
520        for c in range(class_count):
521            ind = (classData == c)
522            numpy.putmask(fxs, ind, mean_destination_vectors[c][0]-x_positions)
523            numpy.putmask(fys, ind, mean_destination_vectors[c][1]-y_positions)
524
525        # compute gradient for all anchors
526        gxs = numpy.array([sum(fxs * selected_data[i])
527                           for i in range(len(anchor_data))], numpy.float)
528        gys = numpy.array([sum(fys * selected_data[i])
529                           for i in range(len(anchor_data))], numpy.float)
530
531        m = max(max(abs(gxs)), max(abs(gys)))
532        gxs /= (20*m); gys /= (20*m)
533
534        newxanchors = xanchors + gxs
535        newyanchors = yanchors + gys
536
537        # normalize so that the anchor most far away will lie on the circle
538        m = math.sqrt(max(newxanchors**2 + newyanchors**2))
539        newxanchors /= m
540        newyanchors /= m
541
542        #self.parentWidget.updateGraph()
543
544        """
545        for a in range(len(anchor_data)):
546            x = anchor_data[a][0]; y = anchor_data[a][1];
547            self.parentWidget.graph.addCurve("lll%i" % i, QColor(0, 0, 0), QColor(0, 0, 0), 10, style = QwtPlotCurve.Lines, symbol = QwtSymbol.NoSymbol, xData = [x, x+gxs[a]], yData = [y, y+gys[a]], forceFilledSymbols = 1, lineWidth=3)
548
549        for i in range(class_count):
550            self.parentWidget.graph.addCurve("lll%i" % i, QColor(0, 0, 0), QColor(0, 0, 0), 10, style = QwtPlotCurve.Lines, symbol = QwtSymbol.NoSymbol, xData = [averages[i][0], mean_destination_vectors[i][0]], yData = [averages[i][1], mean_destination_vectors[i][1]], forceFilledSymbols = 1, lineWidth=3)
551            self.parentWidget.graph.addCurve("lll%i" % i, QColor(0, 0, 0), QColor(0, 0, 0), 10, style = QwtPlotCurve.Lines, xData = [averages[i][0], averages[i][0]], yData = [averages[i][1], averages[i][1]], forceFilledSymbols = 1, lineWidth=5)
552        """
553        #self.parentWidget.graph.repaint()
554        #self.graph.anchor_data = [(newxanchors[i], newyanchors[i], anchor_data[i][2]) for i in range(len(anchor_data))]
555        #self.graph.updateData(attrs, 0)
556        return [(newxanchors[i], newyanchors[i], anchor_data[i][2])
557                for i in range(len(anchor_data))], (newxanchors, newyanchors)
558
559    optimize_LDA_Separation = optimize_lda_separation
560
561    @deprecated_keywords({"attrIndices": "attr_indices",
562                          "anchorData": "anchor_data",
563                          "XAnchors": "xanchors",
564                          "YAnchors": "yanchors"})
565    def optimize_slow_separation(self, attr_indices, anchor_data, xanchors = None, yanchors = None):
566        if (not self.graph.have_data or len(self.graph.raw_data) == 0
567            or not self.graph.data_has_discrete_class): 
568            return anchor_data, (xanchors, yanchors)
569        valid_data = self.graph.get_valid_list(attr_indices)
570        selected_data = numpy.compress(valid_data, numpy.take(self.graph.no_jittering_scaled_data,
571                                                              attr_indices,
572                                                              axis = 0),
573                                       axis = 1)
574
575        if xanchors == None:
576            xanchors = numpy.array([a[0] for a in anchor_data], numpy.float)
577        if yanchors == None:
578            yanchors = numpy.array([a[1] for a in anchor_data], numpy.float)
579
580        trans_proj_data = self.graph.create_projection_as_numeric_array(attr_indices,
581                                                                    validData = valid_data,
582                                                                    xanchors = xanchors,
583                                                                    yanchors = yanchors,
584                                                                    scaleFactor = self.graph.scale_factor,
585                                                                    normalize = self.graph.normalize_examples,
586                                                                    useAnchorData = 1)
587        if trans_proj_data == None:
588            return anchor_data, (xanchors, yanchors)
589
590        proj_data = numpy.transpose(trans_proj_data)
591        x_positions = proj_data[0]; x_positions2 = numpy.array(x_positions)
592        y_positions = proj_data[1]; y_positions2 = numpy.array(y_positions)
593        class_data = proj_data[2]  ; class_data2 = numpy.array(class_data)
594
595        fxs = numpy.zeros(len(x_positions), numpy.float)        # forces
596        fys = numpy.zeros(len(x_positions), numpy.float)
597        gxs = numpy.zeros(len(anchor_data), numpy.float)        # gradients
598        gys = numpy.zeros(len(anchor_data), numpy.float)
599
600        rotate_array = range(len(x_positions))
601        rotate_array = rotate_array[1:] + [0]
602        for i in range(len(x_positions)-1):
603            x_positions2 = numpy.take(x_positions2, rotate_array)
604            y_positions2 = numpy.take(y_positions2, rotate_array)
605            class_data2 = numpy.take(class_data2, rotate_array)
606            dx = x_positions2 - x_positions
607            dy = y_positions2 - y_positions
608            rs2 = dx**2 + dy**2
609            rs2 += numpy.where(rs2 == 0.0, 0.0001, 0.0)    # replace zeros to avoid divisions by zero
610            rs = numpy.sqrt(rs2)
611
612            F = numpy.zeros(len(x_positions), numpy.float)
613            classDiff = numpy.where(class_data == class_data2, 1, 0)
614            numpy.putmask(F, classDiff, 150*self.attract_g*rs2)
615            numpy.putmask(F, 1-classDiff, -self.repel_g/rs2)
616            fxs += F * dx / rs
617            fys += F * dy / rs
618
619        # compute gradient for all anchors
620        gxs = numpy.array([sum(fxs * selected_data[i])
621                           for i in range(len(anchor_data))], numpy.float)
622        gys = numpy.array([sum(fys * selected_data[i])
623                           for i in range(len(anchor_data))], numpy.float)
624
625        m = max(max(abs(gxs)), max(abs(gys)))
626        gxs /= (20*m); gys /= (20*m)
627
628        newxanchors = xanchors + gxs
629        newyanchors = yanchors + gys
630
631        # normalize so that the anchor most far away will lie on the circle
632        m = math.sqrt(max(newxanchors**2 + newyanchors**2))
633        newxanchors /= m
634        newyanchors /= m
635        return [(newxanchors[i], newyanchors[i], anchor_data[i][2])
636                for i in range(len(anchor_data))], (newxanchors, newyanchors)
637
638    optimize_SLOW_Separation = optimize_slow_separation
639
640
641    @deprecated_keywords({"attrIndices": "attr_indices",
642                          "anchorData": "anchor_data",
643                          "XAnchors": "xanchors",
644                          "YAnchors": "yanchors"})
645    def optimize_lda_separation_3D(self, attr_indices, anchor_data, xanchors = None, yanchors = None, zanchors = None):
646        if (not self.graph.have_data or len(self.graph.raw_data) == 0
647            or not self.graph.data_has_discrete_class): 
648            return anchor_data, (xanchors, yanchors, zanchors)
649        class_count = len(self.graph.data_domain.classVar.values)
650        valid_data = self.graph.get_valid_list(attr_indices)
651        selected_data = numpy.compress(valid_data,
652                                       numpy.take(self.graph.no_jittering_scaled_data,
653                                                  attr_indices, axis = 0),
654                                       axis = 1)
655
656        if xanchors == None:
657            xanchors = numpy.array([a[0] for a in anchor_data], numpy.float)
658        if yanchors == None:
659            yanchors = numpy.array([a[1] for a in anchor_data], numpy.float)
660        if zanchors == None:
661            zanchors = numpy.array([a[2] for a in anchor_data], numpy.float)
662
663        trans_proj_data = self.graph.create_projection_as_numeric_array(attr_indices,
664                                                                    validData = valid_data,
665                                                                    xanchors = xanchors,
666                                                                    yanchors = yanchors,
667                                                                    zanchors = zanchors,
668                                                                    scaleFactor = self.graph.scale_factor,
669                                                                    normalize = self.graph.normalize_examples,
670                                                                    useAnchorData = 1)
671        if trans_proj_data == None:
672            return anchor_data, (xanchors, yanchors, zanchors)
673
674        proj_data = numpy.transpose(trans_proj_data)
675        x_positions, y_positions, z_positions, classData = (proj_data[0],
676                                                            proj_data[1],
677                                                            proj_data[2],
678                                                            proj_data[3])
679
680        averages = []
681        for i in range(class_count):
682            ind = classData == i
683            xpos = numpy.compress(ind, x_positions)
684            ypos = numpy.compress(ind, y_positions)
685            zpos = numpy.compress(ind, z_positions)
686            xave = numpy.sum(xpos)/len(xpos)
687            yave = numpy.sum(ypos)/len(ypos)
688            zave = numpy.sum(zpos)/len(zpos)
689            averages.append((xave, yave, zave))
690
691        # compute the positions of all the points. we will try to move all points so that the center will be in the (0,0)
692        x_center_vector = -numpy.sum(x_positions) / len(x_positions)
693        y_center_vector = -numpy.sum(y_positions) / len(y_positions)
694        z_center_vector = -numpy.sum(z_positions) / len(z_positions)
695        center_vector_length = math.sqrt(x_center_vector*x_center_vector +
696                                         y_center_vector*y_center_vector +
697                                         z_center_vector*z_center_vector)
698
699        mean_destination_vectors = []
700
701        for i in range(class_count):
702            xdir = 0.0; ydir = 0.0; zdir = 0.0; rs = 0.0
703            for j in range(class_count):
704                if i==j: continue
705                r = math.sqrt((averages[i][0] - averages[j][0])**2 +
706                              (averages[i][1] - averages[j][1])**2)
707                if r == 0.0:
708                    xdir += math.cos((i/float(class_count))*2*math.pi)
709                    ydir += math.sin((i/float(class_count))*2*math.pi)
710                    r = 0.0001
711                else:
712                    xdir += (1/r**3) * ((averages[i][0] - averages[j][0]))
713                    ydir += (1/r**3) * ((averages[i][1] - averages[j][1]))
714                #rs += 1/r
715            #actualDirAmpl = math.sqrt(xDir**2 + yDir**2)
716            #s = abs(xDir)+abs(yDir)
717            #xDir = rs * (xDir/s)
718            #yDir = rs * (yDir/s)
719            mean_destination_vectors.append((xdir, ydir))
720
721
722        maxlength = math.sqrt(max([x**2 + y**2 for (x,y)
723                                   in mean_destination_vectors]))
724        mean_destination_vectors = [(x/(2*maxlength), y/(2*maxlength)) for (x,y)
725                                    in mean_destination_vectors]     # normalize destination vectors to some normal values
726        mean_destination_vectors = [(mean_destination_vectors[i][0]+averages[i][0],
727                                     mean_destination_vectors[i][1]+averages[i][1])
728                                    for i in range(len(mean_destination_vectors))]    # add destination vectors to the class averages
729        #mean_destination_vectors = [(x + x_center_vector/5, y + y_center_vector/5) for (x,y) in mean_destination_vectors]   # center mean values
730        mean_destination_vectors = [(x + x_center_vector, y + y_center_vector)
731                                    for (x,y) in mean_destination_vectors]   # center mean values
732
733        fxs = numpy.zeros(len(x_positions), numpy.float)        # forces
734        fys = numpy.zeros(len(x_positions), numpy.float)
735
736        for c in range(class_count):
737            ind = (classData == c)
738            numpy.putmask(fxs, ind, mean_destination_vectors[c][0]-x_positions)
739            numpy.putmask(fys, ind, mean_destination_vectors[c][1]-y_positions)
740
741        # compute gradient for all anchors
742        gxs = numpy.array([sum(fxs * selected_data[i])
743                           for i in range(len(anchor_data))], numpy.float)
744        gys = numpy.array([sum(fys * selected_data[i])
745                           for i in range(len(anchor_data))], numpy.float)
746
747        m = max(max(abs(gxs)), max(abs(gys)))
748        gxs /= (20*m); gys /= (20*m)
749
750        newxanchors = xanchors + gxs
751        newyanchors = yanchors + gys
752
753        # normalize so that the anchor most far away will lie on the circle
754        m = math.sqrt(max(newxanchors**2 + newyanchors**2))
755        newxanchors /= m
756        newyanchors /= m
757
758        #self.parentWidget.updateGraph()
759
760        """
761        for a in range(len(anchor_data)):
762            x = anchor_data[a][0]; y = anchor_data[a][1];
763            self.parentWidget.graph.addCurve("lll%i" % i, QColor(0, 0, 0), QColor(0, 0, 0), 10, style = QwtPlotCurve.Lines, symbol = QwtSymbol.NoSymbol, xData = [x, x+gxs[a]], yData = [y, y+gys[a]], forceFilledSymbols = 1, lineWidth=3)
764
765        for i in range(class_count):
766            self.parentWidget.graph.addCurve("lll%i" % i, QColor(0, 0, 0), QColor(0, 0, 0), 10, style = QwtPlotCurve.Lines, symbol = QwtSymbol.NoSymbol, xData = [averages[i][0], mean_destination_vectors[i][0]], yData = [averages[i][1], mean_destination_vectors[i][1]], forceFilledSymbols = 1, lineWidth=3)
767            self.parentWidget.graph.addCurve("lll%i" % i, QColor(0, 0, 0), QColor(0, 0, 0), 10, style = QwtPlotCurve.Lines, xData = [averages[i][0], averages[i][0]], yData = [averages[i][1], averages[i][1]], forceFilledSymbols = 1, lineWidth=5)
768        """
769        #self.parentWidget.graph.repaint()
770        #self.graph.anchor_data = [(newxanchors[i], newyanchors[i], anchor_data[i][2]) for i in range(len(anchor_data))]
771        #self.graph.updateData(attrs, 0)
772        return [(newxanchors[i], newyanchors[i], anchor_data[i][2])
773                for i in range(len(anchor_data))], (newxanchors, newyanchors)
774
775    optimize_LDA_Separation_3D = optimize_lda_separation_3D
776
777    @deprecated_keywords({"attrIndices": "attr_indices",
778                          "anchorData": "anchor_data",
779                          "XAnchors": "xanchors",
780                          "YAnchors": "yanchors"})
781    def optimize_slow_separation_3D(self, attr_indices, anchor_data, xanchors = None, yanchors = None, zanchors = None):
782        if (not self.graph.have_data or len(self.graph.raw_data) == 0
783            or not self.graph.data_has_discrete_class): 
784            return anchor_data, (xanchors, yanchors, zanchors)
785        valid_data = self.graph.get_valid_list(attr_indices)
786        selected_data = numpy.compress(valid_data, numpy.take(self.graph.no_jittering_scaled_data,
787                                                              attr_indices,
788                                                              axis = 0),
789                                       axis = 1)
790
791        if xanchors == None:
792            xanchors = numpy.array([a[0] for a in anchor_data], numpy.float)
793        if yanchors == None:
794            yanchors = numpy.array([a[1] for a in anchor_data], numpy.float)
795        if zanchors == None:
796            zanchors = numpy.array([a[2] for a in anchor_data], numpy.float)
797
798        trans_proj_data = self.graph.create_projection_as_numeric_array(attr_indices,
799                                                                    validData = valid_data,
800                                                                    XAnchors = xanchors,
801                                                                    YAnchors = yanchors,
802                                                                    ZAnchors = zanchors,
803                                                                    scaleFactor = self.graph.scale_factor,
804                                                                    normalize = self.graph.normalize_examples,
805                                                                    useAnchorData = 1)
806        if trans_proj_data == None:
807            return anchor_data, (xanchors, yanchors, zanchors)
808
809        proj_data = numpy.transpose(trans_proj_data)
810        x_positions = proj_data[0]; x_positions2 = numpy.array(x_positions)
811        y_positions = proj_data[1]; y_positions2 = numpy.array(y_positions)
812        z_positions = proj_data[2]; z_positions2 = numpy.array(z_positions)
813        class_data = proj_data[3];  class_data2 = numpy.array(class_data)
814
815        fxs = numpy.zeros(len(x_positions), numpy.float)        # forces
816        fys = numpy.zeros(len(x_positions), numpy.float)
817        fzs = numpy.zeros(len(x_positions), numpy.float)
818        gxs = numpy.zeros(len(anchor_data), numpy.float)        # gradients
819        gys = numpy.zeros(len(anchor_data), numpy.float)
820        gzs = numpy.zeros(len(anchor_data), numpy.float)
821
822        rotate_array = range(len(x_positions))
823        rotate_array = rotate_array[1:] + [0]
824        for i in range(len(x_positions)-1):
825            x_positions2 = numpy.take(x_positions2, rotate_array)
826            y_positions2 = numpy.take(y_positions2, rotate_array)
827            z_positions2 = numpy.take(z_positions2, rotate_array)
828            class_data2 = numpy.take(class_data2, rotate_array)
829            dx = x_positions2 - x_positions
830            dy = y_positions2 - y_positions
831            dz = z_positions2 - z_positions
832            rs2 = dx**2 + dy**2 + dz**2
833            rs2 += numpy.where(rs2 == 0.0, 0.0001, 0.0)    # replace zeros to avoid divisions by zero
834            rs = numpy.sqrt(rs2)
835
836            F = numpy.zeros(len(x_positions), numpy.float)
837            classDiff = numpy.where(class_data == class_data2, 1, 0)
838            numpy.putmask(F, classDiff, 150*self.attract_g*rs2)
839            numpy.putmask(F, 1-classDiff, -self.repel_g/rs2)
840            fxs += F * dx / rs
841            fys += F * dy / rs
842            fzs += F * dz / rs
843
844        # compute gradient for all anchors
845        gxs = numpy.array([sum(fxs * selected_data[i])
846                           for i in range(len(anchor_data))], numpy.float)
847        gys = numpy.array([sum(fys * selected_data[i])
848                           for i in range(len(anchor_data))], numpy.float)
849        gzs = numpy.array([sum(fzs * selected_data[i])
850                           for i in range(len(anchor_data))], numpy.float)
851
852        m = max(max(abs(gxs)), max(abs(gys)), max(abs(gzs)))
853        gxs /= (20*m)
854        gys /= (20*m)
855        gzs /= (20*m)
856
857        newxanchors = xanchors + gxs
858        newyanchors = yanchors + gys
859        newzanchors = zanchors + gzs
860
861        # normalize so that the anchor most far away will lie on the circle
862        m = math.sqrt(max(newxanchors**2 + newyanchors**2 + newzanchors**2))
863        newxanchors /= m
864        newyanchors /= m
865        newzanchors /= m
866        return [(newxanchors[i], newyanchors[i], newzanchors[i], anchor_data[i][3])
867                for i in range(len(anchor_data))], (newxanchors, newyanchors, newzanchors)
868
869    optimize_SLOW_Separation_3D = optimize_slow_separation_3D
870
871
872
873    # ###############################################################
874    # S2N HEURISTIC FUNCTIONS
875    # ###############################################################
876
877
878
879    # place a subset of attributes around the circle. this subset must contain "good" attributes for each of the class values
880    @deprecated_keywords({"setAttributeListInRadviz":
881                          "set_attribute_list_in_radviz"})
882    def s2n_mix_anchors(self, set_attribute_list_in_radviz = 1):
883        # check if we have data and a discrete class
884        if (not self.graph.have_data or len(self.graph.raw_data) == 0
885            or not self.graph.data_has_discrete_class): 
886            self.set_statusbar_text("S2N only works on data with a discrete class value")
887            return
888
889        # compute the quality of attributes only once
890        if self.s2n_mix_data == None:
891            ranked_attrs, ranked_attrs_by_class = visfuncts.findAttributeGroupsForRadviz(self.graph.raw_data,
892                                                                                         visfuncts.S2NMeasureMix())
893            self.s2n_mix_data = (ranked_attrs, ranked_attrs_by_class)
894            class_count = len(ranked_attrs_by_class)
895            attrs = ranked_attrs[:(self.s2n_place_attributes/class_count)*
896                                 class_count]    # select appropriate number of attributes
897        else:
898            class_count = len(self.s2n_mix_data[1])
899            attrs = self.s2n_mix_data[0][:(self.s2n_place_attributes/class_count)*
900                                         class_count]
901
902        if len(attrs) == 0:
903            self.set_statusbar_text("No discrete attributes found")
904            return 0
905
906        arr = [0]       # array that will tell where to put the next attribute
907        for i in range(1,len(attrs)/2): arr += [i,-i]
908
909        phi = (2*math.pi*self.s2n_spread)/(len(attrs)*10.0)
910        anchor_data = []; start = []
911        arr2 = arr[:(len(attrs)/class_count)+1]
912        for cls in range(class_count):
913            start_pos = (2*math.pi*cls)/class_count
914            if self.class_permutation_list: cls = self.class_permutation_list[cls]
915            attrs_cls = attrs[cls::class_count]
916            temp_data = [(arr2[i], math.cos(start_pos + arr2[i]*phi),
917                          math.sin(start_pos + arr2[i]*phi),
918                          attrs_cls[i]) for i in
919                          range(min(len(arr2), len(attrs_cls)))]
920            start.append(len(anchor_data) + len(arr2)/2) # starting indices for each class value
921            temp_data.sort()
922            anchor_data += [(x, y, name) for (i, x, y, name) in temp_data]
923
924        anchor_data = anchor_data[(len(attrs)/(2*class_count)):] + anchor_data[:(len(attrs)/(2*class_count))]
925        self.graph.anchor_data = anchor_data
926        attrNames = [anchor[2] for anchor in anchor_data]
927
928        if self.__class__ != FreeViz:
929            if set_attribute_list_in_radviz:
930                self.parentWidget.setShownAttributeList(attrNames)
931            self.graph.updateData(attrNames)
932            self.graph.repaint()
933        return 1
934
935    s2nMixAnchors = s2n_mix_anchors
936
937    # find interesting linear projection using PCA, SPCA, or PLS
938    @deprecated_keywords({"attrIndices": "attr_indices",
939                          "setAnchors": "set_anchors",
940                          "percentDataUsed": "percent_data_used"})
941    def find_projection(self, method, attr_indices = None, set_anchors = 0, percent_data_used = 100):
942        if not self.graph.have_data: return
943        ai = self.graph.attribute_name_index
944        if attr_indices == None:
945            attributes = self.get_shown_attribute_list()
946            attr_indices = [ai[label] for label in attributes]
947        if len(attr_indices) == 0: return None
948
949        valid_data = self.graph.get_valid_list(attr_indices)
950        if sum(valid_data) == 0: return None
951
952        data_matrix = numpy.compress(valid_data, numpy.take(self.graph.no_jittering_scaled_data,
953                                                            attr_indices,
954                                                            axis = 0),
955                                     axis = 1)
956        if self.graph.data_has_class:
957            class_array = numpy.compress(valid_data,
958                                         self.graph.no_jittering_scaled_data[self.graph.data_class_index])
959
960        if percent_data_used != 100:
961            indices = Orange.data.sample.SubsetIndices2(self.graph.raw_data,
962                                                1.0-(float(percent_data_used)/100.0))
963            try:
964                data_matrix = numpy.compress(indices, data_matrix, axis = 1)
965            except:
966                pass
967            if self.graph.data_has_class:
968                class_array = numpy.compress(indices, class_array)
969
970        ncomps = 3 if hasattr(self, '_use_3D') else 2
971        vectors = None
972        if method == DR_PCA:
973            pca = Pca(standardize=False, max_components=ncomps,
974                use_generalized_eigenvectors=0)
975            domain = Orange.data.Domain([Orange.feature.Continuous("g%d"%i) for i
976                                         in xrange(len(data_matrix))], False)
977            pca = pca(Orange.data.Table(domain, data_matrix.T))
978            vals, vectors = pca.eigen_values, pca.projection
979        elif method == DR_SPCA and self.graph.data_has_class:
980            pca = Spca(standardize=False, max_components=ncomps,
981                use_generalized_eigenvectors=self.use_generalized_eigenvectors)
982            domain = Orange.data.Domain([Orange.feature.Continuous("g%d"%i) for i
983                                         in xrange(len(data_matrix))], Orange.feature.Continuous("c"))
984            pca = pca(Orange.data.Table(domain,
985                numpy.hstack([data_matrix.T, numpy.array(class_array, ndmin=2).T])))
986            vals, vectors = pca.eigen_values, pca.projection
987        elif method == DR_PLS and self.graph.data_has_class:
988            data_matrix = data_matrix.transpose()
989            class_matrix = numpy.transpose(numpy.matrix(class_array))
990            vectors = create_pls_projection(data_matrix, class_matrix, ncomps)
991            vectors = vectors.T
992
993        # test if all values are 0, if there is an invalid number in the array and if there are complex numbers in the array
994        if (vectors == None or not vectors.any() or
995            False in numpy.isfinite(vectors) or False in numpy.isreal(vectors)):
996            self.set_statusbar_text("Unable to compute anchor positions for the selected attributes") 
997            return None
998
999        xanchors = vectors[0]
1000        yanchors = vectors[1]
1001       
1002        if ncomps == 3:
1003            zanchors = vectors[2]
1004            m = math.sqrt(max(xanchors**2 + yanchors**2 + zanchors**2))
1005            zanchors /= m
1006        else:
1007            m = math.sqrt(max(xanchors**2 + yanchors**2))
1008
1009        xanchors /= m
1010        yanchors /= m
1011        names = self.graph.attribute_names
1012        attributes = [names[attr_indices[i]] for i in range(len(attr_indices))]
1013
1014        if set_anchors:
1015            if ncomps == 3:
1016                self.graph.set_anchors(list(xanchors), list(yanchors), list(zanchors), attributes)
1017            else:
1018                self.graph.set_anchors(list(xanchors), list(yanchors), attributes)
1019            if hasattr(self.graph, "updateData"):
1020                self.graph.updateData()
1021            if hasattr(self.graph, "repaint"):
1022                self.graph.repaint()
1023
1024        if ncomps == 3:
1025            return xanchors, yanchors, zanchors, (attributes, attr_indices)
1026        else:
1027            return xanchors, yanchors, (attributes, attr_indices)
1028
1029    findProjection = find_projection
1030
1031
1032FreeViz = deprecated_members({"attractG": "attract_g",
1033                              "repelG": "repel_g",
1034                              "forceBalancing": "force_balancing",
1035                              "forceSigma": "force_sigma",
1036                              "mirrorSymmetry": "mirror_symmetry",
1037                              "useGeneralizedEigenvectors": "use_generalized_eigenvectors",
1038                              "stepsBeforeUpdate": "steps_before_update",
1039                              "s2nSpread": "s2n_spread",
1040                              "s2nPlaceAttributes": "s2n_place_attributes",
1041                              "s2nMixData": "s2n_mix_data",
1042                              "autoSetParameters": "auto_set_parameters",
1043                              "classPermutationList": "class_permutation_list",
1044                              "attrsNum": "attrs_num",
1045                              "cancelOptimization": "cancel_optimization"})(FreeViz)
1046
1047
1048@deprecated_keywords({"X": "x", "Y": "y", "Ncomp": "ncomp"})
1049def create_pls_projection(x,y, ncomp = 2):
1050    '''Predict y from x using first ncomp principal components'''
1051
1052    # data dimensions
1053    n, mx = numpy.shape(x)
1054    my = numpy.shape(y)[1]
1055
1056    # Z-scores of original matrices
1057    ymean = y.mean()
1058    x,y = center(x), center(y)
1059
1060    p = numpy.empty((mx,ncomp))
1061    w = numpy.empty((mx,ncomp))
1062    c = numpy.empty((my,ncomp))
1063    t = numpy.empty((n,ncomp))
1064    u = numpy.empty((n,ncomp))
1065    b = numpy.zeros((ncomp,ncomp))
1066
1067    e,f = x,y
1068
1069    # main algorithm
1070    for i in range(ncomp):
1071
1072        u = numpy.random.random_sample((n,1))
1073        w = normalize(numpy.dot(e.T,u))
1074        t = normalize(numpy.dot(e,w))
1075        c = normalize(numpy.dot(f.T,t))
1076
1077        dif = t
1078        # iterations for loading vector t
1079        while numpy.linalg.norm(dif) > 10e-16:
1080            c = normalize(numpy.dot(f.T,t))
1081            u = numpy.dot(f,c)
1082            w = normalize(numpy.dot(e.T,u))
1083            t0 = normalize(numpy.dot(e,w))
1084            dif = t - t0
1085            t = t0
1086
1087        t[:,i] = t.T
1088        u[:,i] = u.T
1089        c[:,i] = c.T
1090        w[:,i] = w.T
1091
1092        b = numpy.dot(t.T,u)[0,0]
1093        b[i][i] = b
1094        p = numpy.dot(e.T,t)
1095        p[:,i] = p.T
1096        e = e - numpy.dot(t,p.T)
1097        xx = b * numpy.dot(t,c.T)
1098        f = f - xx
1099
1100    # esimated y
1101    #YE = numpy.dot(numpy.dot(t,b),c.t)*numpy.std(y, axis = 0) + ymean
1102    #y = y*numpy.std(y, axis = 0)+ ymean
1103    #BPls = numpy.dot(numpy.dot(numpy.linalg.pinv(p.t),b),c.t)
1104
1105    return w
1106
1107createPLSProjection = create_pls_projection
1108
1109# #############################################################################
1110# class that represents freeviz classifier
1111class FreeVizClassifier(Orange.classification.Classifier):
1112    """
1113    A kNN classifier on the 2D projection of the data, optimized by FreeViz.
1114   
1115    Usually the learner
1116    (:class:`Orange.projection.linear.FreeVizLearner`) is used to construct the
1117    classifier.
1118   
1119    When constructing the classifier manually, the following parameters can
1120    be passed:
1121   
1122    :param data: table of data instances to project to a 2D plane and use for
1123        classification.
1124    :type data: :class:`Orange.data.Table`
1125   
1126    :param freeviz: the FreeViz algorithm instance to use to optimize the 2D
1127        projection.
1128    :type freeviz: :class:`Orange.projection.linear.FreeViz`
1129   
1130    """
1131   
1132    def __init__(self, data, freeviz):
1133        self.freeviz = freeviz
1134
1135        if self.freeviz.__class__ != FreeViz:
1136            self.freeviz.parentWidget.setData(data)
1137            self.freeviz.parentWidget.showAllAttributes = 1
1138        else:
1139            self.freeviz.graph.set_data(data)
1140            self.freeviz.show_all_attributes()
1141
1142        #self.FreeViz.randomAnchors()
1143        self.freeviz.radial_anchors()
1144        self.freeviz.optimize_separation()
1145
1146        graph = self.freeviz.graph
1147        ai = graph.attribute_name_index
1148        labels = [a[2] for a in graph.anchor_data]
1149        indices = [ai[label] for label in labels]
1150
1151        valid_data = graph.get_valid_list(indices)
1152        domain = Orange.data.Domain([graph.data_domain[i].name for i in indices]+
1153                               [graph.data_domain.classVar.name],
1154                               graph.data_domain)
1155        offsets = [graph.attr_values[graph.attribute_names[i]][0]
1156                   for i in indices]
1157        normalizers = [graph.get_min_max_val(i) for i in indices]
1158        selected_data = numpy.take(graph.original_data, indices, axis = 0)
1159        averages = numpy.average(numpy.compress(valid_data, selected_data,
1160                                                axis=1), 1)
1161        class_data = numpy.compress(valid_data,
1162                                    graph.original_data[graph.data_class_index])
1163
1164        graph.create_projection_as_numeric_array(indices, use_anchor_data = 1,
1165                                             remove_missing_data = 0,
1166                                             valid_data = valid_data,
1167                                             jitter_size = -1)
1168        self.classifier = Orange.classification.knn.P2NN(domain,
1169                                      numpy.transpose(numpy.array([numpy.compress(valid_data,
1170                                                                                  graph.unscaled_x_positions),
1171                                                                   numpy.compress(valid_data,
1172                                                                                  graph.unscaled_y_positions),
1173                                                                   class_data])),
1174                                      graph.anchor_data, offsets, normalizers,
1175                                      averages, graph.normalize_examples, law=1)
1176
1177    # for a given instance run argumentation and find out to which class it most often fall
1178    @deprecated_keywords({"example": "instance", "returnType": "return_type"})
1179    def __call__(self, instance, return_type=Orange.classification.Classifier.GetValue):
1180        #instance.setclass(0)
1181        return self.classifier(instance, return_type)
1182
1183FreeVizClassifier = deprecated_members({"FreeViz":"freeviz"})(FreeVizClassifier)
1184
1185class FreeVizLearner(Orange.classification.Learner):
1186    """
1187    A learner that builds a :class:`FreeVizClassifier` on given data. An
1188    instance of :class:`FreeViz` can be passed to the constructor as a
1189    keyword argument :obj:`freeviz`.   
1190
1191    If data instances are provided to the constructor, the learning algorithm
1192    is called and the resulting classifier is returned instead of the learner.
1193   
1194    """
1195    def __new__(cls, freeviz = None, instances = None, weight_id = 0, **argkw):
1196        self = Orange.classification.Learner.__new__(cls, **argkw)
1197        if instances:
1198            self.__init__(freeviz, **argkw)
1199            return self.__call__(instances, weight_id)
1200        else:
1201            return self
1202
1203    def __init__(self, freeviz = None):
1204        if not freeviz:
1205            freeviz = FreeViz()
1206        self.freeviz = freeviz
1207        self.name = "freeviz Learner"
1208
1209    @deprecated_keywords({"examples": "instances", "weightID": "weight_id"})
1210    def __call__(self, instances, weight_id = 0):
1211        return FreeVizClassifier(instances, self.freeviz)
1212
1213FreeVizLearner = deprecated_members({"FreeViz":"freeviz"})(FreeVizLearner)
1214
1215
1216class S2NHeuristicLearner(Orange.classification.Learner):
1217    """
1218    This class is not documented yet.
1219   
1220    """
1221    def __new__(cls, freeviz = None, instances = None, weight_id = 0, **argkw):
1222        self = Orange.classification.Learner.__new__(cls, **argkw)
1223        if instances:
1224            self.__init__(freeviz, **argkw)
1225            return self.__call__(instances, weight_id)
1226        else:
1227            return self
1228
1229    def __init__(self, freeviz = None):
1230        if not freeviz:
1231            freeviz = FreeViz()
1232        self.freeviz = freeviz
1233        self.name = "S2N Feature Selection Learner"
1234
1235    @deprecated_keywords({"examples": "instances", "weightID": "weight_id"})
1236    def __call__(self, instances, weight_id = 0):
1237        return S2NHeuristicClassifier(instances, self.freeviz)
1238
1239S2NHeuristicLearner = deprecated_members({"FreeViz":
1240                                          "freeviz"})(S2NHeuristicLearner)
1241
1242class Projector(object):
1243    """
1244    Stores a linear projection of data and uses it to transform any given data with matching input domain.
1245
1246    .. attribute:: input_domain
1247
1248        Domain of the data set that was used to construct principal component
1249        subspace.
1250
1251    .. attribute:: output_domain
1252
1253        Domain used in returned data sets. This domain has a continuous
1254        variable for each axis in the projected space,
1255        and no class variable(s).
1256
1257    .. attribute:: mean
1258
1259        Array containing means of each variable in the data set that was used
1260        to construct the projection.
1261
1262    .. attribute:: stdev
1263
1264        An array containing standard deviations of each variable in the data
1265        set that was used to construct the projection.
1266
1267    .. attribute:: standardize
1268
1269        True, if standardization was used when constructing the projection. If
1270        set, instances will be standardized before being projected.
1271
1272    .. attribute:: projection
1273
1274        Array containing projection (vectors that describe the
1275        transformation from input to output domain).
1276
1277    """
1278    def __init__(self, **kwds):
1279        self.__dict__.update(kwds)
1280        if not hasattr(self, "output_domain"):
1281            self.output_domain = Orange.data.Domain([Orange.feature.Continuous("a.%d"%(i+1)) for i in range(len(self.projection))], False)
1282
1283
1284    def __call__(self, data):
1285        """
1286        Project data.
1287
1288        :param data: input data set
1289        :type data: :class:`Orange.data.Table`
1290
1291        :rtype: :class:`Orange.data.Table`
1292        """
1293        if type(data) != Orange.data.Table:
1294            data = Orange.data.Table([data])
1295        if len(self.projection.T) != len(data.domain.features):
1296            data = Orange.data.Table(self.input_domain, data)
1297
1298        X = data.to_numpy_MA("a")[0]
1299        Xm, U = self.mean, self.projection
1300        n, m = X.shape
1301
1302        if m != len(self.projection.T):
1303            raise Orange.core.KernelException, "Invalid number of features"
1304
1305        Xd = X - Xm
1306
1307        if self.standardize:
1308            Xd /= self.stdev
1309
1310        self.A = numpy.ma.dot(Xd, U.T)
1311
1312        return Orange.data.Table(self.output_domain, self.A.tolist())
1313
1314#color table for biplot
1315Colors = ['bo','go','yo','co','mo']
1316
1317class Pca(object):
1318    """
1319    Orthogonal transformation of data into a set of uncorrelated variables called
1320    principal components. This transformation is defined in such a way that the
1321    first variable has as high variance as possible.
1322
1323    If data instances are provided to the constructor,
1324    the optimization algorithm is called and the resulting projector
1325    (:class:`~Orange.projection.linear.PcaProjector`) is
1326    returned instead of the optimizer (instance of this class).
1327
1328    :param standardize: perform standardization of the data set.
1329    :type standardize: boolean
1330    :param max_components: maximum number of retained components.
1331    :type max_components: int
1332    :param variance_covered: percent of the variance to cover with components.
1333    :type variance_covered: float
1334    :param use_generalized_eigenvectors: use generalized eigenvectors (ie.
1335        multiply data matrix with inverse of its covariance matrix).
1336    :type use_generalized_eigenvectors: boolean
1337
1338    :rtype: :class:`~Orange.projection.linear.Pca` or
1339            :class:`~Orange.projection.linear.PcaProjector`
1340    """
1341
1342    def __new__(cls, dataset = None, **kwds):
1343        optimizer = object.__new__(cls)
1344        optimizer.__init__(**kwds)
1345
1346        if dataset:
1347            return optimizer(dataset)
1348        else:
1349            return optimizer
1350
1351    def __init__(self, standardize = True,
1352                 max_components = 0, variance_covered = 1,
1353                 use_generalized_eigenvectors = 0):
1354        self.standardize = standardize
1355        self.max_components = max_components
1356        self.variance_covered = variance_covered if variance_covered < 1. else 1
1357        self.use_generalized_eigenvectors = use_generalized_eigenvectors
1358
1359    def _pca(self, dataset, Xd, Xg):
1360        n,m = Xd.shape
1361        if n < m:
1362            C = numpy.ma.dot(Xg.T, Xd.T)
1363            V, D, T = numpy.linalg.svd(C)
1364            U = numpy.ma.dot(V.T, Xd) / numpy.sqrt(D.reshape(-1,1))
1365        else:
1366            C = numpy.ma.dot(Xg, Xd)
1367            U, D, T = numpy.linalg.svd(C)
1368
1369        U = U.# eigenvectors are now in rows
1370        return U, D
1371
1372    def __call__(self, dataset):
1373        """
1374        Perform a PCA analysis on a data set and return a linear projector
1375        that maps data into principal component subspace.
1376
1377        :param dataset: input data set.
1378        :type dataset: :class:`Orange.data.Table`
1379
1380        :rtype: :class:`~Orange.projection.linear.PcaProjector`
1381        """
1382
1383        X = dataset.to_numpy_MA("a")[0]
1384        N,M = X.shape
1385        Xm = numpy.mean(X, axis=0)
1386        Xd = X - Xm
1387
1388        #take care of the constant features
1389        stdev = numpy.std(Xd, axis=0)
1390        relevant_features = stdev != 0
1391        if self.standardize:
1392            stdev[stdev == 0] = 1.
1393            Xd /= stdev
1394        Xd = Xd[:,relevant_features]
1395
1396        #use generalized eigenvectors
1397        if self.use_generalized_eigenvectors:
1398            inv_covar = numpy.linalg.inv(numpy.dot(Xd.T, Xd))
1399            Xg = numpy.dot(inv_covar, Xd.T)
1400        else:
1401            Xg = Xd.T
1402
1403        #actual pca
1404        n,m = Xd.shape
1405        U, D = self._pca(dataset, Xd, Xg)
1406
1407        #insert zeros for constant features
1408        n, m = U.shape
1409        if m != M:
1410            U_ = numpy.zeros((n,M))
1411            U_[:,relevant_features] = U
1412            U = U_
1413
1414        variance_sum = D.sum()
1415
1416        #select eigen vectors
1417        if self.variance_covered != 1:
1418            nfeatures = numpy.nonzero(numpy.cumsum(D) / sum(D) >= self.variance_covered)[0][0] + 1
1419            U = U[:nfeatures, :]
1420            D = D[:nfeatures]
1421
1422        if self.max_components > 0:
1423            U = U[:self.max_components, :]
1424            D = D[:self.max_components]
1425
1426        n, m = U.shape
1427        pc_domain = Orange.data.Domain([Orange.feature.Continuous("Comp.%d"%
1428                                                                  (i+1)) for i in range(n)], False)
1429
1430        return PcaProjector(input_domain = dataset.domain,
1431            output_domain = pc_domain,
1432            pc_domain = pc_domain,
1433            mean = Xm,
1434            stdev = stdev,
1435            standardize = self.standardize,
1436            eigen_vectors = U,
1437            projection = U,
1438            eigen_values = D,
1439            variance_sum = variance_sum)
1440
1441
1442class Spca(Pca):
1443    def _pca(self, dataset, Xd, Xg):
1444        # define the Laplacian matrix
1445        c = dataset.to_numpy("c")[0]
1446        l = -numpy.array(numpy.hstack( [(c != v) for v in c]), dtype='f')
1447        l -= numpy.diagflat(numpy.sum(l, axis=0))
1448
1449        Xg = numpy.dot(Xg, l)
1450
1451        return Pca._pca(self, dataset, Xd, Xg)
1452
1453class PcaProjector(Projector):
1454    """
1455    .. attribute:: pc_domain
1456
1457        Synonymous for :obj:`~Orange.projection.linear.Projector.output_domain`.
1458
1459    .. attribute:: eigen_vectors
1460
1461        Synonymous for :obj:`~Orange.projection.linear.Projector.projection`.
1462
1463    .. attribute:: eigen_values
1464
1465        Array containing standard deviations of principal components.
1466
1467    .. attribute:: variance_sum
1468
1469        Sum of all variances in the data set that was used to construct the PCA
1470        space.
1471
1472    """
1473
1474    def __init__(self, **kwds):
1475        self.__dict__.update(kwds)
1476
1477    def __str__(self):
1478        ncomponents = 10
1479        s = self.variance_sum
1480        cs = numpy.cumsum(self.eigen_values) / s
1481        return "\n".join([
1482            "PCA SUMMARY",
1483            "",
1484            "Std. deviation of components:",
1485            " ".join(["              "] +
1486                     ["%10s" % a.name for a in self.pc_domain.attributes]),
1487            " ".join(["Std. deviation"] +
1488                     ["%10.3f" % a for a in self.eigen_values]),
1489            " ".join(["Proportion Var"] +
1490                     ["%10.3f" % a for a in  self.eigen_values / s * 100]),
1491            " ".join(["Cumulative Var"] +
1492                     ["%10.3f" % a for a in cs * 100]),
1493            "",
1494            #"Loadings:",
1495            #" ".join(["%10s"%""] + ["%10s" % a.name for a in self.pc_domain]),
1496            #"\n".join([
1497            #    " ".join([a.name] + ["%10.3f" % b for b in self.eigen_vectors.T[i]])
1498            #          for i, a in enumerate(self.input_domain.attributes)
1499            #          ])
1500        ]) if len(self.pc_domain) <= ncomponents else\
1501        "\n".join([
1502            "PCA SUMMARY",
1503            "",
1504            "Std. deviation of components:",
1505            " ".join(["              "] +
1506                     ["%10s" % a.name for a in self.pc_domain.attributes[:ncomponents]] +
1507                     ["%10s" % "..."] +
1508                     ["%10s" % self.pc_domain.attributes[-1].name]),
1509            " ".join(["Std. deviation"] +
1510                     ["%10.3f" % a for a in self.eigen_values[:ncomponents]] +
1511                     ["%10s" % ""] +
1512                     ["%10.3f" % self.eigen_values[-1]]),
1513            " ".join(["Proportion Var"] +
1514                     ["%10.3f" % a for a in self.eigen_values[:ncomponents] / s * 100] +
1515                     ["%10s" % ""] +
1516                     ["%10.3f" % (self.eigen_values[-1] / s * 100)]),
1517            " ".join(["Cumulative Var"] +
1518                     ["%10.3f" % a for a in cs[:ncomponents] * 100] +
1519                     ["%10s" % ""] +
1520                     ["%10.3f" % (cs[-1] * 100)]),
1521            "",
1522            #"Loadings:",
1523            #" ".join(["%16s" % ""] +
1524            #         ["%8s" % a.name for a in self.pc_domain.attributes[:ncomponents]] +
1525            #         ["%8s" % "..."] +
1526            #         ["%8s" % self.pc_domain.attributes[-1].name]),
1527            #"\n".join([
1528            #    " ".join(["%16.16s" %a.name] +
1529            #             ["%8.3f" % b for b in self.eigen_vectors.T[i, :ncomponents]] +
1530            #             ["%8s" % ""] +
1531            #             ["%8.3f" % self.eigen_vectors.T[i, -1]])
1532            #          for i, a in enumerate(self.input_domain.attributes)
1533            #          ])
1534        ])
1535
1536
1537
1538    ################ Plotting functions ###################
1539
1540    def scree_plot(self, filename = None, title = 'Scree plot'):
1541        """
1542        Draw a scree plot of principal components
1543
1544        :param filename: Name of the file to which the plot will be saved. \
1545        If None, plot will be displayed instead.
1546        :type filename: str
1547        :param title: Plot title
1548        :type title: str
1549        """
1550        import pylab as plt
1551
1552        s = self.variance_sum
1553        vc = self.eigen_values / s
1554        cs = numpy.cumsum(self.eigen_values) / s
1555
1556        fig = plt.figure()
1557        ax = fig.add_subplot(111)
1558
1559        x_axis = range(len(self.eigen_values))
1560        x_labels = ["PC%d" % (i + 1, ) for i in x_axis]
1561
1562        ax.set_xticks(x_axis)
1563        ax.set_xticklabels(x_labels)
1564        plt.setp(ax.get_xticklabels(), "rotation", 90)
1565        plt.grid(True)
1566
1567        ax.set_xlabel('Principal components')
1568        ax.set_ylabel('Proportion of Variance')
1569        ax.set_title(title + "\n")
1570        ax.plot(x_axis, vc, color="red")
1571        ax.scatter(x_axis, vc, color="red", label="Variance")
1572
1573        ax.plot(x_axis, cs, color="orange")
1574        ax.scatter(x_axis, cs, color="orange", label="Cumulative Variance")
1575        ax.legend(loc=0)
1576
1577        ax.axis([-0.5, len(self.eigen_values) - 0.5, 0, 1])
1578
1579        if filename:
1580            plt.savefig(filename)
1581        else:
1582            plt.show()
1583
1584    def biplot(self, filename = None, components = [0,1], title = 'Biplot'):
1585        """
1586        Draw biplot for PCA. Actual projection must be performed via pca(data)
1587        before bipot can be used.
1588
1589        :param filename: Name of the file to which the plot will be saved. \
1590        If None, plot will be displayed instead.
1591        :type plot: str
1592        :param components: List of two components to plot.
1593        :type components: list
1594        :param title: Plot title
1595        :type title: str
1596        """
1597        import pylab as plt
1598
1599        if len(components) < 2:
1600            raise orange.KernelException, 'Two components are needed for biplot'
1601
1602        if not (0 <= min(components) <= max(components) < len(self.eigen_values)):
1603            raise orange.KernelException, 'Invalid components'
1604
1605        X = self.A[:,components[0]]
1606        Y = self.A[:,components[1]]
1607
1608        vectorsX = self.eigen_vectors[:,components[0]]
1609        vectorsY = self.eigen_vectors[:,components[1]]
1610
1611
1612        #TO DO -> pc.biplot (maybe)
1613        #trDataMatrix = dataMatrix / lam
1614        #trLoadings = loadings * lam
1615
1616        #max_data_value = numpy.max(abs(trDataMatrix)) * 1.05
1617        max_load_value = self.eigen_vectors.max() * 1.5
1618
1619        #plt.clf()
1620        fig = plt.figure()
1621        ax1 = fig.add_subplot(111)
1622        ax1.set_title(title + "\n")
1623        ax1.set_xlabel("PC%s (%d%%)" % (components[0], self.eigen_values[components[0]] / self.variance_sum * 100))
1624        ax1.set_ylabel("PC%s (%d%%)" % (components[1], self.eigen_values[components[1]] / self.variance_sum * 100))
1625        ax1.xaxis.set_label_position('bottom')
1626        ax1.xaxis.set_ticks_position('bottom')
1627        ax1.yaxis.set_label_position('left')
1628        ax1.yaxis.set_ticks_position('left')
1629
1630        #if self._classArray == None:
1631        #trDataMatrix = transpose(trDataMatrix)
1632        ax1.plot(X, Y, Colors[0])
1633        #else:
1634        #suboptimal
1635        #    classValues = []
1636        #    for classValue in self._classArray:
1637        #        if classValue not in classValues:
1638        #            classValues.append(classValue)
1639        #    for i in range(len(classValues)):
1640        #        choice = numpy.array([classValues[i] == cv for cv in self._classArray])
1641        #        partialDataMatrix = transpose(trDataMatrix[choice])
1642        #        ax1.plot(partialDataMatrix[0], partialDataMatrix[1],
1643        #                 Colors[i % len(Colors)], label = str(classValues[i]))
1644        #    ax1.legend()
1645
1646        #ax1.set_xlim(-max_data_value, max_data_value)
1647        #ax1.set_ylim(-max_data_value, max_data_value)
1648
1649        #eliminate double axis on right
1650        ax0 = ax1.twinx()
1651        ax0.yaxis.set_visible(False)
1652
1653        ax2 = ax0.twiny()
1654        ax2.xaxis.set_label_position('top')
1655        ax2.xaxis.set_ticks_position('top')
1656        ax2.yaxis.set_label_position('right')
1657        ax2.yaxis.set_ticks_position('right')
1658        for tl in ax2.get_xticklabels():
1659            tl.set_color('r')
1660        for tl in ax2.get_yticklabels():
1661            tl.set_color('r')
1662
1663        arrowprops = dict(facecolor = 'red', edgecolor = 'red', width = 1, headwidth = 4)
1664
1665        for (x, y, a) in zip(vectorsX, vectorsY,self.input_domain.attributes):
1666            if max(x, y) < 0.1:
1667                continue
1668            print x, y, a
1669            ax2.annotate('', (x, y), (0, 0), arrowprops = arrowprops)
1670            ax2.text(x * 1.1, y * 1.2, a.name, color = 'red')
1671
1672        ax2.set_xlim(-max_load_value, max_load_value)
1673        ax2.set_ylim(-max_load_value, max_load_value)
1674
1675        if filename:
1676            plt.savefig(filename)
1677        else:
1678            plt.show()
1679
1680
1681class Fda(object):
1682    """
1683    Construct a linear projection of data using FDA. When using this projection optimization method, data is always
1684    standardized prior to being projected.
1685
1686    If data instances are provided to the constructor,
1687    the optimization algorithm is called and the resulting projector
1688    (:class:`~Orange.projection.linear.FdaProjector`) is
1689    returned instead of the optimizer (instance of this class).
1690
1691    :rtype: :class:`~Orange.projection.linear.Fda` or
1692            :class:`~Orange.projection.linear.FdaProjector`
1693    """
1694
1695    def __new__(cls, data = None):
1696        self = object.__new__(cls)
1697        if data:
1698            self.__init__()
1699            return self.__call__(data)
1700        else:
1701            return self
1702
1703    def __call__(self, dataset):
1704        """
1705        Perform a FDA analysis on a data set and return a linear projector
1706        that maps data into another vector space.
1707
1708        :param dataset: input data set.
1709        :type dataset: :class:`Orange.data.Table`
1710
1711        :rtype: :class:`~Orange.projection.linear.FdaProjector`
1712        """
1713        X, Y = dataset.to_numpy_MA("a/c")
1714
1715        Xm = numpy.mean(X, axis=0)
1716        X = X - Xm
1717
1718        #take care of the constant features
1719        stdev = numpy.std(X, axis=0)
1720        relevant_features = stdev != 0
1721        stdev[stdev == 0] = 1.
1722        X /= stdev
1723        X = X[:,relevant_features]
1724
1725        instances, features = X.shape
1726        class_count = len(set(Y))
1727        # special case when we have two classes
1728        if class_count == 2:
1729            data1 = MA.take(X, numpy.argwhere(Y == 0).flatten(), axis=0)
1730            data2 = MA.take(X, numpy.argwhere(Y != 0).flatten(), axis=0)
1731            miDiff = MA.average(data1, axis=1) - MA.average(data2, axis=1)
1732            covMatrix = (MA.dot(data1.T, data1) + MA.dot(data2.T, data2)) / instances
1733            U = numpy.linalg.inv(covMatrix) * miDiff
1734            D = numpy.array([1])
1735        else:
1736            # compute means and average covariances of examples in each class group
1737            Sw = MA.zeros([features, features])
1738            for v in set(Y):
1739                d = MA.take(X, numpy.argwhere(Y == v).flatten(), axis=0)
1740                d = d - numpy.mean(d, axis=0)
1741                Sw += MA.dot(d.T, d)
1742            Sw /= instances
1743            total = MA.dot(X.T, X)/float(instances)
1744            Sb = total - Sw
1745
1746            matrix = numpy.linalg.inv(Sw)*Sb
1747            D, U = numpy.linalg.eigh(matrix)
1748
1749        sorted_indices = [i for _,i in sorted([(ev, i)
1750                          for i, ev in enumerate(D)], reverse=True)]
1751        U = numpy.take(U, sorted_indices, axis = 1)
1752        D = numpy.take(D, sorted_indices)
1753
1754        #insert zeros for constant features
1755        n, m = U.shape
1756        if m != M:
1757            U_ = numpy.zeros((n,M))
1758            U_[:,relevant_features] = U
1759            U = U_
1760
1761        out_domain = Orange.data.Domain([Orange.feature.Continuous("Comp.%d"%
1762                                                                  (i+1)) for
1763                                         i in range(len(D))], False)
1764
1765        return FdaProjector(input_domain = dataset.domain,
1766            output_domain = out_domain,
1767            mean = Xm,
1768            stdev = stdev,
1769            standardize = True,
1770            eigen_vectors = U,
1771            projection = U,
1772            eigen_values = D)
1773
1774class FdaProjector(Projector):
1775    """
1776    .. attribute:: eigen_vectors
1777
1778        Synonymous for :obj:`~Orange.projection.linear.Projector.projection`.
1779
1780    .. attribute:: eigen_values
1781
1782        Array containing eigenvalues corresponding to eigenvectors.
1783
1784    """
1785
1786    def __init__(self, **kwds):
1787        self.__dict__.update(kwds)
1788
1789
1790
1791@deprecated_keywords({"dataMatrix": "data_matrix",
1792                      "classArray": "class_array",
1793                      "NComps": "ncomps",
1794                      "useGeneralizedEigenvectors": "use_generalized_eigenvectors"})
1795def create_pca_projection(data_matrix, class_array = None, ncomps = -1, use_generalized_eigenvectors = 1):
1796    import warnings
1797    warnings.warn("Deprecated in favour of Orange"
1798                  ".projection.linear.Pca.",
1799        DeprecationWarning)
1800    if type(data_matrix) == numpy.ma.core.MaskedArray:
1801        data_matrix = numpy.array(data_matrix)
1802    if class_array != None and type(class_array) == numpy.ma.core.MaskedArray:
1803        class_array = numpy.array(class_array)
1804
1805    data_matrix = numpy.transpose(data_matrix)
1806
1807    s = numpy.sum(data_matrix, axis=0)/float(len(data_matrix))
1808    data_matrix -= s       # substract average value to get zero mean
1809
1810    if class_array != None and use_generalized_eigenvectors:
1811        covarMatrix = numpy.dot(numpy.transpose(data_matrix), data_matrix)
1812        try:
1813            matrix = inv(covarMatrix)
1814        except:
1815            return None, None
1816        matrix = numpy.dot(matrix, numpy.transpose(data_matrix))
1817    else:
1818        matrix = numpy.transpose(data_matrix)
1819
1820    # compute dataMatrixT * L * dataMatrix
1821    if class_array != None:
1822        # define the Laplacian matrix
1823        l = numpy.zeros((len(data_matrix), len(data_matrix)))
1824        for i in range(len(data_matrix)):
1825            for j in range(i+1, len(data_matrix)):
1826                l[i,j] = -int(class_array[i] != class_array[j])
1827                l[j,i] = -int(class_array[i] != class_array[j])
1828
1829        s = numpy.sum(l, axis=0)      # doesn't matter which axis since the matrix l is symmetrical
1830        for i in range(len(data_matrix)):
1831            l[i,i] = -s[i]
1832
1833        matrix = numpy.dot(matrix, l)
1834
1835    matrix = numpy.dot(matrix, data_matrix)
1836
1837    vals, vectors = numpy.linalg.eig(matrix)
1838    if vals.dtype.kind == "c":       # if eigenvalues are complex numbers then do nothing
1839        return None, None
1840    vals = list(vals)
1841
1842    if ncomps == -1:
1843        ncomps = len(vals)
1844    ncomps = min(ncomps, len(vals))
1845
1846    ret_vals = []
1847    ret_indices = []
1848    for i in range(ncomps):
1849        ret_vals.append(max(vals))
1850        bestind = vals.index(max(vals))
1851        ret_indices.append(bestind)
1852        vals[bestind] = -1
1853
1854    return ret_vals, numpy.take(vectors.T, ret_indices, axis = 0)         # i-th eigenvector is the i-th column in vectors so we have to transpose the array
1855
1856createPCAProjection = create_pca_projection
Note: See TracBrowser for help on using the repository browser.