# source:orange/Orange/projection/linear.py@10490:a1479e654132

Revision 10490:a1479e654132, 76.9 KB checked in by Lan Zagar <lan.zagar@…>, 2 years ago (diff)

Minor improvements.

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
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
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]
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
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
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):
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)
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)
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)
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)
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
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:
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
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:
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 is 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",
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
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()
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, max_components=0, variance_covered=1,
1352                 use_generalized_eigenvectors=0):
1353        self.standardize = standardize
1354        self.max_components = max_components
1355        self.variance_covered = min(1, variance_covered)
1356        self.use_generalized_eigenvectors = use_generalized_eigenvectors
1357
1358    def _pca(self, dataset, Xd, Xg):
1359        n,m = Xd.shape
1360        if n < m:
1361            C = numpy.ma.dot(Xg.T, Xd.T)
1362            V, D, T = numpy.linalg.svd(C)
1363            U = numpy.ma.dot(V.T, Xd) / numpy.sqrt(D.reshape(-1, 1))
1364        else:
1365            C = numpy.ma.dot(Xg, Xd)
1366            U, D, T = numpy.linalg.svd(C)
1367
1368        U = U.# eigenvectors are now in rows
1369        return U, D
1370
1371    def __call__(self, dataset):
1372        """
1373        Perform a PCA analysis on a data set and return a linear projector
1374        that maps data into principal component subspace.
1375
1376        :param dataset: input data set.
1377        :type dataset: :class:`Orange.data.Table`
1378
1379        :rtype: :class:`~Orange.projection.linear.PcaProjector`
1380        """
1381
1382        X = dataset.to_numpy_MA("a")[0]
1383        N,M = X.shape
1384        Xm = numpy.mean(X, axis=0)
1385        Xd = X - Xm
1386
1387        #take care of the constant features
1388        stdev = numpy.std(Xd, axis=0)
1389        relevant_features = stdev != 0
1390        Xd = Xd[:, relevant_features]
1391        if self.standardize:
1392            Xd /= stdev[relevant_features]
1393
1394        #use generalized eigenvectors
1395        if self.use_generalized_eigenvectors:
1396            inv_covar = numpy.linalg.inv(numpy.dot(Xd.T, Xd))
1397            Xg = numpy.dot(inv_covar, Xd.T)
1398        else:
1399            Xg = Xd.T
1400
1401        #actual pca
1402        n, m = Xd.shape
1403        U, D = self._pca(dataset, Xd, Xg)
1404
1405        #insert zeros for constant features
1406        n, m = U.shape
1407        if m != M:
1408            U_ = numpy.zeros((n, M))
1409            U_[:, relevant_features] = U
1410            U = U_
1411
1412        variance_sum = D.sum()
1413
1414        #select eigen vectors
1415        if self.variance_covered != 1:
1416            nfeatures = numpy.searchsorted(numpy.cumsum(D) / variance_sum,
1417                                           self.variance_covered) + 1
1418            U = U[:nfeatures, :]
1419            D = D[:nfeatures]
1420
1421        if self.max_components > 0:
1422            U = U[:self.max_components, :]
1423            D = D[:self.max_components]
1424
1425        n, m = U.shape
1426        pc_domain = Orange.data.Domain([Orange.feature.Continuous("Comp.%d"%
1427            (i + 1)) for i in range(n)], False)
1428
1429        return PcaProjector(input_domain=dataset.domain,
1430            output_domain = pc_domain,
1431            pc_domain = pc_domain,
1432            mean = Xm,
1433            stdev = stdev,
1434            standardize = self.standardize,
1435            eigen_vectors = U,
1436            projection = U,
1437            eigen_values = D,
1438            variance_sum = variance_sum)
1439
1440
1441class Spca(Pca):
1442    def _pca(self, dataset, Xd, Xg):
1443        # define the Laplacian matrix
1444        c = dataset.to_numpy("c")[0]
1445        l = -numpy.array(numpy.hstack([(c != v) for v in c]), dtype='f')
1446        l -= numpy.diagflat(numpy.sum(l, axis=0))
1447
1448        Xg = numpy.dot(Xg, l)
1449
1450        return Pca._pca(self, dataset, Xd, Xg)
1451
1452class PcaProjector(Projector):
1453    """
1454    .. attribute:: pc_domain
1455
1456        Synonymous for :obj:`~Orange.projection.linear.Projector.output_domain`.
1457
1458    .. attribute:: eigen_vectors
1459
1460        Synonymous for :obj:`~Orange.projection.linear.Projector.projection`.
1461
1462    .. attribute:: eigen_values
1463
1464        Array containing standard deviations of principal components.
1465
1466    .. attribute:: variance_sum
1467
1468        Sum of all variances in the data set that was used to construct the PCA
1469        space.
1470
1471    """
1472
1473    def __init__(self, **kwds):
1474        self.__dict__.update(kwds)
1475
1476    def __str__(self):
1477        ncomponents = 10
1478        s = self.variance_sum
1479        cs = numpy.cumsum(self.eigen_values) / s
1480        return "\n".join([
1481            "PCA SUMMARY",
1482            "",
1483            "Std. deviation of components:",
1484            " ".join(["              "] +
1485                     ["%10s" % a.name for a in self.pc_domain.attributes]),
1486            " ".join(["Std. deviation"] +
1487                     ["%10.3f" % a for a in self.eigen_values]),
1488            " ".join(["Proportion Var"] +
1489                     ["%10.3f" % a for a in  self.eigen_values / s * 100]),
1490            " ".join(["Cumulative Var"] +
1491                     ["%10.3f" % a for a in cs * 100]),
1492            "",
1494            #" ".join(["%10s"%""] + ["%10s" % a.name for a in self.pc_domain]),
1495            #"\n".join([
1496            #    " ".join([a.name] + ["%10.3f" % b for b in self.eigen_vectors.T[i]])
1497            #          for i, a in enumerate(self.input_domain.attributes)
1498            #          ])
1499        ]) if len(self.pc_domain) <= ncomponents else\
1500        "\n".join([
1501            "PCA SUMMARY",
1502            "",
1503            "Std. deviation of components:",
1504            " ".join(["              "] +
1505                     ["%10s" % a.name for a in self.pc_domain.attributes[:ncomponents]] +
1506                     ["%10s" % "..."] +
1507                     ["%10s" % self.pc_domain.attributes[-1].name]),
1508            " ".join(["Std. deviation"] +
1509                     ["%10.3f" % a for a in self.eigen_values[:ncomponents]] +
1510                     ["%10s" % ""] +
1511                     ["%10.3f" % self.eigen_values[-1]]),
1512            " ".join(["Proportion Var"] +
1513                     ["%10.3f" % a for a in self.eigen_values[:ncomponents] / s * 100] +
1514                     ["%10s" % ""] +
1515                     ["%10.3f" % (self.eigen_values[-1] / s * 100)]),
1516            " ".join(["Cumulative Var"] +
1517                     ["%10.3f" % a for a in cs[:ncomponents] * 100] +
1518                     ["%10s" % ""] +
1519                     ["%10.3f" % (cs[-1] * 100)]),
1520            "",
1522            #" ".join(["%16s" % ""] +
1523            #         ["%8s" % a.name for a in self.pc_domain.attributes[:ncomponents]] +
1524            #         ["%8s" % "..."] +
1525            #         ["%8s" % self.pc_domain.attributes[-1].name]),
1526            #"\n".join([
1527            #    " ".join(["%16.16s" %a.name] +
1528            #             ["%8.3f" % b for b in self.eigen_vectors.T[i, :ncomponents]] +
1529            #             ["%8s" % ""] +
1530            #             ["%8.3f" % self.eigen_vectors.T[i, -1]])
1531            #          for i, a in enumerate(self.input_domain.attributes)
1532            #          ])
1533        ])
1534
1535
1536
1537    ################ Plotting functions ###################
1538
1539    def scree_plot(self, filename = None, title = 'Scree plot'):
1540        """
1541        Draw a scree plot of principal components
1542
1543        :param filename: Name of the file to which the plot will be saved. \
1544        If None, plot will be displayed instead.
1545        :type filename: str
1546        :param title: Plot title
1547        :type title: str
1548        """
1549        import pylab as plt
1550
1551        s = self.variance_sum
1552        vc = self.eigen_values / s
1553        cs = numpy.cumsum(self.eigen_values) / s
1554
1555        fig = plt.figure()
1557
1558        x_axis = range(len(self.eigen_values))
1559        x_labels = ["PC%d" % (i + 1, ) for i in x_axis]
1560
1561        ax.set_xticks(x_axis)
1562        ax.set_xticklabels(x_labels)
1563        plt.setp(ax.get_xticklabels(), "rotation", 90)
1564        plt.grid(True)
1565
1566        ax.set_xlabel('Principal components')
1567        ax.set_ylabel('Proportion of Variance')
1568        ax.set_title(title + "\n")
1569        ax.plot(x_axis, vc, color="red")
1570        ax.scatter(x_axis, vc, color="red", label="Variance")
1571
1572        ax.plot(x_axis, cs, color="orange")
1573        ax.scatter(x_axis, cs, color="orange", label="Cumulative Variance")
1574        ax.legend(loc=0)
1575
1576        ax.axis([-0.5, len(self.eigen_values) - 0.5, 0, 1])
1577
1578        if filename:
1579            plt.savefig(filename)
1580        else:
1581            plt.show()
1582
1583    def biplot(self, filename = None, components = [0,1], title = 'Biplot'):
1584        """
1585        Draw biplot for PCA. Actual projection must be performed via pca(data)
1586        before bipot can be used.
1587
1588        :param filename: Name of the file to which the plot will be saved. \
1589        If None, plot will be displayed instead.
1590        :type plot: str
1591        :param components: List of two components to plot.
1592        :type components: list
1593        :param title: Plot title
1594        :type title: str
1595        """
1596        import pylab as plt
1597
1598        if len(components) < 2:
1599            raise orange.KernelException, 'Two components are needed for biplot'
1600
1601        if not (0 <= min(components) <= max(components) < len(self.eigen_values)):
1602            raise orange.KernelException, 'Invalid components'
1603
1604        X = self.A[:,components[0]]
1605        Y = self.A[:,components[1]]
1606
1607        vectorsX = self.eigen_vectors[:,components[0]]
1608        vectorsY = self.eigen_vectors[:,components[1]]
1609
1610
1611        #TO DO -> pc.biplot (maybe)
1612        #trDataMatrix = dataMatrix / lam
1614
1615        #max_data_value = numpy.max(abs(trDataMatrix)) * 1.05
1616        max_load_value = self.eigen_vectors.max() * 1.5
1617
1618        #plt.clf()
1619        fig = plt.figure()
1621        ax1.set_title(title + "\n")
1622        ax1.set_xlabel("PC%s (%d%%)" % (components[0], self.eigen_values[components[0]] / self.variance_sum * 100))
1623        ax1.set_ylabel("PC%s (%d%%)" % (components[1], self.eigen_values[components[1]] / self.variance_sum * 100))
1624        ax1.xaxis.set_label_position('bottom')
1625        ax1.xaxis.set_ticks_position('bottom')
1626        ax1.yaxis.set_label_position('left')
1627        ax1.yaxis.set_ticks_position('left')
1628
1629        #if self._classArray == None:
1630        #trDataMatrix = transpose(trDataMatrix)
1631        ax1.plot(X, Y, Colors[0])
1632        #else:
1633        #suboptimal
1634        #    classValues = []
1635        #    for classValue in self._classArray:
1636        #        if classValue not in classValues:
1637        #            classValues.append(classValue)
1638        #    for i in range(len(classValues)):
1639        #        choice = numpy.array([classValues[i] == cv for cv in self._classArray])
1640        #        partialDataMatrix = transpose(trDataMatrix[choice])
1641        #        ax1.plot(partialDataMatrix[0], partialDataMatrix[1],
1642        #                 Colors[i % len(Colors)], label = str(classValues[i]))
1643        #    ax1.legend()
1644
1645        #ax1.set_xlim(-max_data_value, max_data_value)
1646        #ax1.set_ylim(-max_data_value, max_data_value)
1647
1648        #eliminate double axis on right
1649        ax0 = ax1.twinx()
1650        ax0.yaxis.set_visible(False)
1651
1652        ax2 = ax0.twiny()
1653        ax2.xaxis.set_label_position('top')
1654        ax2.xaxis.set_ticks_position('top')
1655        ax2.yaxis.set_label_position('right')
1656        ax2.yaxis.set_ticks_position('right')
1657        for tl in ax2.get_xticklabels():
1658            tl.set_color('r')
1659        for tl in ax2.get_yticklabels():
1660            tl.set_color('r')
1661
1662        arrowprops = dict(facecolor = 'red', edgecolor = 'red', width = 1, headwidth = 4)
1663
1664        for (x, y, a) in zip(vectorsX, vectorsY,self.input_domain.attributes):
1665            if max(x, y) < 0.1:
1666                continue
1667            print x, y, a
1668            ax2.annotate('', (x, y), (0, 0), arrowprops = arrowprops)
1669            ax2.text(x * 1.1, y * 1.2, a.name, color = 'red')
1670
1673
1674        if filename:
1675            plt.savefig(filename)
1676        else:
1677            plt.show()
1678
1679
1680class Fda(object):
1681    """
1682    Construct a linear projection of data using FDA. When using this projection optimization method, data is always
1683    standardized prior to being projected.
1684
1685    If data instances are provided to the constructor,
1686    the optimization algorithm is called and the resulting projector
1687    (:class:`~Orange.projection.linear.FdaProjector`) is
1688    returned instead of the optimizer (instance of this class).
1689
1690    :rtype: :class:`~Orange.projection.linear.Fda` or
1691            :class:`~Orange.projection.linear.FdaProjector`
1692    """
1693
1694    def __new__(cls, data = None):
1695        self = object.__new__(cls)
1696        if data:
1697            self.__init__()
1698            return self.__call__(data)
1699        else:
1700            return self
1701
1702    def __call__(self, dataset):
1703        """
1704        Perform a FDA analysis on a data set and return a linear projector
1705        that maps data into another vector space.
1706
1707        :param dataset: input data set.
1708        :type dataset: :class:`Orange.data.Table`
1709
1710        :rtype: :class:`~Orange.projection.linear.FdaProjector`
1711        """
1712        X, Y = dataset.to_numpy_MA("a/c")
1713
1714        Xm = numpy.mean(X, axis=0)
1715        X = X - Xm
1716
1717        #take care of the constant features
1718        stdev = numpy.std(X, axis=0)
1719        relevant_features = stdev != 0
1720        stdev[stdev == 0] = 1.
1721        X /= stdev
1722        X = X[:,relevant_features]
1723
1724        instances, features = X.shape
1725        class_count = len(set(Y))
1726        # special case when we have two classes
1727        if class_count == 2:
1728            data1 = MA.take(X, numpy.argwhere(Y == 0).flatten(), axis=0)
1729            data2 = MA.take(X, numpy.argwhere(Y != 0).flatten(), axis=0)
1730            miDiff = MA.average(data1, axis=1) - MA.average(data2, axis=1)
1731            covMatrix = (MA.dot(data1.T, data1) + MA.dot(data2.T, data2)) / instances
1732            U = numpy.linalg.inv(covMatrix) * miDiff
1733            D = numpy.array([1])
1734        else:
1735            # compute means and average covariances of examples in each class group
1736            Sw = MA.zeros([features, features])
1737            for v in set(Y):
1738                d = MA.take(X, numpy.argwhere(Y == v).flatten(), axis=0)
1739                d = d - numpy.mean(d, axis=0)
1740                Sw += MA.dot(d.T, d)
1741            Sw /= instances
1742            total = MA.dot(X.T, X)/float(instances)
1743            Sb = total - Sw
1744
1745            matrix = numpy.linalg.inv(Sw)*Sb
1746            D, U = numpy.linalg.eigh(matrix)
1747
1748        sorted_indices = [i for _,i in sorted([(ev, i)
1749                          for i, ev in enumerate(D)], reverse=True)]
1750        U = numpy.take(U, sorted_indices, axis = 1)
1751        D = numpy.take(D, sorted_indices)
1752
1753        #insert zeros for constant features
1754        n, m = U.shape
1755        if m != M:
1756            U_ = numpy.zeros((n,M))
1757            U_[:,relevant_features] = U
1758            U = U_
1759
1760        out_domain = Orange.data.Domain([Orange.feature.Continuous("Comp.%d"%
1761                                                                  (i+1)) for
1762                                         i in range(len(D))], False)
1763
1764        return FdaProjector(input_domain = dataset.domain,
1765            output_domain = out_domain,
1766            mean = Xm,
1767            stdev = stdev,
1768            standardize = True,
1769            eigen_vectors = U,
1770            projection = U,
1771            eigen_values = D)
1772
1773class FdaProjector(Projector):
1774    """
1775    .. attribute:: eigen_vectors
1776
1777        Synonymous for :obj:`~Orange.projection.linear.Projector.projection`.
1778
1779    .. attribute:: eigen_values
1780
1781        Array containing eigenvalues corresponding to eigenvectors.
1782
1783    """
1784
1785    def __init__(self, **kwds):
1786        self.__dict__.update(kwds)
1787
1788
1789
1790@deprecated_keywords({"dataMatrix": "data_matrix",
1791                      "classArray": "class_array",
1792                      "NComps": "ncomps",
1793                      "useGeneralizedEigenvectors": "use_generalized_eigenvectors"})
1794def create_pca_projection(data_matrix, class_array = None, ncomps = -1, use_generalized_eigenvectors = 1):
1795    import warnings
1796    warnings.warn("Deprecated in favour of Orange"
1797                  ".projection.linear.Pca.",
1798        DeprecationWarning)
1800        data_matrix = numpy.array(data_matrix)
1801    if class_array != None and type(class_array) == numpy.ma.core.MaskedArray:
1802        class_array = numpy.array(class_array)
1803
1804    data_matrix = numpy.transpose(data_matrix)
1805
1806    s = numpy.sum(data_matrix, axis=0)/float(len(data_matrix))
1807    data_matrix -= s       # substract average value to get zero mean
1808
1809    if class_array != None and use_generalized_eigenvectors:
1810        covarMatrix = numpy.dot(numpy.transpose(data_matrix), data_matrix)
1811        try:
1812            matrix = inv(covarMatrix)
1813        except:
1814            return None, None
1815        matrix = numpy.dot(matrix, numpy.transpose(data_matrix))
1816    else:
1817        matrix = numpy.transpose(data_matrix)
1818
1819    # compute dataMatrixT * L * dataMatrix
1820    if class_array != None:
1821        # define the Laplacian matrix
1822        l = numpy.zeros((len(data_matrix), len(data_matrix)))
1823        for i in range(len(data_matrix)):
1824            for j in range(i+1, len(data_matrix)):
1825                l[i,j] = -int(class_array[i] != class_array[j])
1826                l[j,i] = -int(class_array[i] != class_array[j])
1827
1828        s = numpy.sum(l, axis=0)      # doesn't matter which axis since the matrix l is symmetrical
1829        for i in range(len(data_matrix)):
1830            l[i,i] = -s[i]
1831
1832        matrix = numpy.dot(matrix, l)
1833
1834    matrix = numpy.dot(matrix, data_matrix)
1835
1836    vals, vectors = numpy.linalg.eig(matrix)
1837    if vals.dtype.kind == "c":       # if eigenvalues are complex numbers then do nothing
1838        return None, None
1839    vals = list(vals)
1840
1841    if ncomps == -1:
1842        ncomps = len(vals)
1843    ncomps = min(ncomps, len(vals))
1844
1845    ret_vals = []
1846    ret_indices = []
1847    for i in range(ncomps):
1848        ret_vals.append(max(vals))
1849        bestind = vals.index(max(vals))
1850        ret_indices.append(bestind)
1851        vals[bestind] = -1
1852
1853    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
1854
1855createPCAProjection = create_pca_projection
Note: See TracBrowser for help on using the repository browser.