source: orange/Orange/projection/linear.py @ 10611:f438273040dc

Revision 10611:f438273040dc, 71.7 KB checked in by anze <anze.staric@…>, 2 years ago (diff)

Improved coding style in projection.linear.
Removed dead code and fixed code spacing. Fixed a bunch of PyCharm warnings.

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