source: orange/Orange/projection/linear.py @ 10715:94cb2d76e3f1

Revision 10715:94cb2d76e3f1, 71.8 KB checked in by anze <anze.staric@…>, 2 years ago (diff)

Merged in patch from #1161.

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                         center=Xm,
166                         scale=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, ddof=0)
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.variances, 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, ddof=0)
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.variances, 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    #: Domain of the data set that was used to construct principal component subspace.
1202    input_domain = None
1203
1204    #: Domain used in returned data sets. This domain has a continuous
1205    #: variable for each axis in the projected space,
1206    #: and no class variable(s).
1207    output_domain = None
1208
1209    #: Array containing means of each variable in the data set that was used
1210    #: to construct the projection.
1211    center = numpy.array(())
1212
1213    #: An array containing standard deviations of each variable in the data
1214    #: set that was used to construct the projection.
1215    scale = numpy.array(())
1216
1217    #: True, if standardization was used when constructing the projection. If
1218    #: set, instances will be standardized before being projected.
1219    standardize = True
1220
1221    #: Array containing projection (vectors that describe the
1222    #: transformation from input to output domain).
1223    projection = numpy.array(()).reshape(0, 0)
1224
1225    def __init__(self, **kwds):
1226        self.__dict__.update(kwds)
1227
1228        features = []
1229        for i in range(len(self.projection)):
1230            f = feature.Continuous("Comp.%d" % (i + 1))
1231            f.get_value_from = lambda ex, w: self._project_single(ex, w, f, i)
1232            features.append(f)
1233
1234        self.output_domain = Orange.data.Domain(features,
1235                                                self.input_domain.class_var,
1236                                                class_vars=self.input_domain.class_vars)
1237
1238    def _project_single(self, example, return_what, new_feature, feature_idx):
1239        ex = Orange.data.Table([example]).to_numpy("a")[0]
1240        ex -= self.center
1241        if self.standardize:
1242            ex /= self.scale
1243        return new_feature(numpy.dot(self.projection[feature_idx, :], ex.T)[0])
1244
1245    def __call__(self, dataset):
1246        """
1247        Project data.
1248
1249        :param dataset: input data set
1250        :type dataset: :class:`Orange.data.Table`
1251
1252        :rtype: :class:`Orange.data.Table`
1253        """
1254        if not isinstance(dataset, data.Table):
1255            dataset = data.Table([dataset])
1256        if len(self.projection.T) != len(dataset.domain.features):
1257            dataset = data.Table(self.input_domain, dataset)
1258
1259        X, = dataset.to_numpy_MA("a")
1260        Xm, U = self.center, self.projection
1261        n, m = X.shape
1262
1263        if m != len(self.projection.T):
1264            raise ValueError, "Invalid number of features"
1265
1266        Xd = X - Xm
1267
1268        if self.standardize:
1269            Xd /= self.scale
1270
1271        self.A = numpy.dot(Xd, U.T)
1272
1273        class_, classes = dataset.to_numpy("c")[0], dataset.to_numpy("m")[0]
1274        return data.Table(self.output_domain, numpy.hstack((self.A, class_, classes)))
1275
1276
1277#color table for biplot
1278Colors = ['bo', 'go', 'yo', 'co', 'mo']
1279
1280class PCA(object):
1281    """
1282    Orthogonal transformation of data into a set of uncorrelated variables called
1283    principal components. This transformation is defined in such a way that the
1284    first variable has as high variance as possible.
1285
1286    :param standardize: perform standardization of the data set.
1287    :type standardize: boolean
1288    :param max_components: maximum number of retained components.
1289    :type max_components: int
1290    :param variance_covered: percent of the variance to cover with components.
1291    :type variance_covered: float
1292    :param use_generalized_eigenvectors: use generalized eigenvectors (ie.
1293        multiply data matrix with inverse of its covariance matrix).
1294    :type use_generalized_eigenvectors: boolean
1295    """
1296
1297    #: Delta degrees of freedom used for numpy operations.
1298    #: 1 means normalization with (N-1) in cov and std operations
1299    ddof = 1
1300
1301    def __new__(cls, dataset=None, **kwds):
1302        optimizer = object.__new__(cls)
1303
1304        if dataset:
1305            optimizer.__init__(**kwds)
1306            return optimizer(dataset)
1307        else:
1308            return optimizer
1309
1310    def __init__(self, standardize=True, max_components=0, variance_covered=1,
1311                 use_generalized_eigenvectors=0, ddof=1):
1312        self.standardize = standardize
1313        self.max_components = max_components
1314        self.variance_covered = min(1, variance_covered)
1315        self.use_generalized_eigenvectors = use_generalized_eigenvectors
1316        self.ddof=ddof
1317
1318    def __call__(self, dataset):
1319        """
1320        Perform a PCA analysis on a data set and return a linear projector
1321        that maps data into principal component subspace.
1322
1323        :param dataset: input data set.
1324        :type dataset: :class:`Orange.data.Table`
1325
1326        :rtype: :class:`~Orange.projection.linear.PcaProjector`
1327        """
1328        Xd, stdev, mean, relevant_features = self._normalize_data(dataset)
1329
1330        #use generalized eigenvectors
1331        if self.use_generalized_eigenvectors:
1332            inv_covar = numpy.linalg.inv(numpy.dot(Xd.T, Xd))
1333            Xg = numpy.dot(inv_covar, Xd.T)
1334        else:
1335            Xg = Xd.T
1336
1337        components, variances = self._perform_pca(dataset, Xd, Xg)
1338        components = self._insert_zeros_for_constant_features(len(dataset.domain.features),
1339                                                              components,
1340                                                              relevant_features)
1341
1342        variances, components, variance_sum = self._select_components(variances, components)
1343
1344        n, m = components.shape
1345
1346        return PcaProjector(input_domain=dataset.domain,
1347                            center=mean,
1348                            scale=stdev,
1349                            standardize=self.standardize,
1350                            projection=components,
1351                            variances=variances,
1352                            variance_sum=variance_sum)
1353
1354    def _normalize_data(self, dataset):
1355        if not len(dataset) or not len(dataset.domain.features):
1356            raise ValueError("Empty dataset")
1357        X = dataset.to_numpy_MA("a")[0]
1358
1359        Xm = numpy.mean(X, axis=0)
1360        Xd = X - Xm
1361
1362        if not Xd.any():
1363            raise ValueError("All features are constant")
1364
1365        #take care of the constant features
1366        stdev = numpy.std(Xd, axis=0, ddof=self.ddof)
1367        stdev[stdev == 0] = 1. # to prevent division by zero
1368        relevant_features = stdev != 0
1369        Xd = Xd[:, relevant_features]
1370        if self.standardize:
1371            Xd /= stdev[relevant_features]
1372        return Xd, stdev, Xm, relevant_features
1373
1374    def _perform_pca(self, dataset, Xd, Xg):
1375        n, m = Xd.shape
1376        if n < m:
1377            _, D, U = numpy.linalg.svd(Xd, full_matrices=0)
1378            D *= D / (n - self.ddof)
1379        else:
1380            C = numpy.dot(Xg, Xd) / (n - self.ddof)
1381            U, D, T = numpy.linalg.svd(C)
1382            U = U.# eigenvectors are now in rows
1383        return U, D
1384
1385    def _select_components(self, D, U):
1386        variance_sum = D.sum()
1387        #select eigen vectors
1388        if self.variance_covered != 1:
1389            nfeatures = numpy.searchsorted(numpy.cumsum(D) / variance_sum,
1390                                           self.variance_covered) + 1
1391            U = U[:nfeatures, :]
1392            D = D[:nfeatures]
1393        if self.max_components > 0:
1394            U = U[:self.max_components, :]
1395            D = D[:self.max_components]
1396        return D, U, variance_sum
1397
1398    def _insert_zeros_for_constant_features(self, original_dimension, U, relevant_features):
1399        n, m = U.shape
1400        if m != original_dimension:
1401            U_ = numpy.zeros((n, original_dimension))
1402            U_[:, relevant_features] = U
1403            U = U_
1404        return U
1405
1406Pca = PCA
1407
1408
1409class Spca(PCA):
1410    def _perform_pca(self, dataset, Xd, Xg):
1411        # define the Laplacian matrix
1412        c = dataset.to_numpy("c")[0]
1413        l = -numpy.array(numpy.hstack([(c != v) for v in c]), dtype='f')
1414        l -= numpy.diagflat(numpy.sum(l, axis=0))
1415
1416        Xg = numpy.dot(Xg, l)
1417
1418        return Pca._perform_pca(self, dataset, Xd, Xg)
1419
1420
1421@deprecated_members({"pc_domain": "output_domain"})
1422class PcaProjector(Projector):
1423    #: Array containing variances of principal components.
1424    variances = numpy.array(())
1425
1426    #: Sum of all variances in the data set that was used to construct the PCA space.
1427    variance_sum = 0.
1428
1429    def __str__(self):
1430        ncomponents = 10
1431        s = self.variance_sum
1432        cs = numpy.cumsum(self.variances) / s
1433        stdev = numpy.sqrt(self.variances)
1434        return "\n".join([
1435            "PCA SUMMARY",
1436            "",
1437            "Std. deviation of components:",
1438            " ".join(["              "] +
1439                     ["%10s" % a.name for a in self.output_domain.attributes]),
1440            " ".join(["Std. deviation"] +
1441                     ["%10.3f" % a for a in stdev]),
1442            " ".join(["Proportion Var"] +
1443                     ["%10.3f" % a for a in self.variances / s * 100]),
1444            " ".join(["Cumulative Var"] +
1445                     ["%10.3f" % a for a in cs * 100]),
1446            "",
1447            ]) if len(self.output_domain) <= ncomponents else\
1448        "\n".join([
1449            "PCA SUMMARY",
1450            "",
1451            "Std. deviation of components:",
1452            " ".join(["              "] +
1453                     ["%10s" % a.name for a in self.output_domain.attributes[:ncomponents]] +
1454                     ["%10s" % "..."] +
1455                     ["%10s" % self.output_domain.attributes[-1].name]),
1456            " ".join(["Std. deviation"] +
1457                     ["%10.3f" % a for a in stdev[:ncomponents]] +
1458                     ["%10s" % ""] +
1459                     ["%10.3f" % stdev[-1]]),
1460            " ".join(["Proportion Var"] +
1461                     ["%10.3f" % a for a in self.variances[:ncomponents] / s * 100] +
1462                     ["%10s" % ""] +
1463                     ["%10.3f" % (self.variances[-1] / s * 100)]),
1464            " ".join(["Cumulative Var"] +
1465                     ["%10.3f" % a for a in cs[:ncomponents] * 100] +
1466                     ["%10s" % ""] +
1467                     ["%10.3f" % (cs[-1] * 100)]),
1468            "",
1469            ])
1470
1471
1472    def scree_plot(self, filename=None, title='Scree Plot'):
1473        """
1474        Draw a scree plot of principal components
1475
1476        :param filename: Name of the file to which the plot will be saved.
1477                         If None, plot will be displayed instead.
1478        :type filename: str
1479        :param title: Plot title
1480        :type title: str
1481        """
1482        import pylab as plt
1483
1484        s = self.variance_sum
1485        vc = self.variances / s
1486        cs = numpy.cumsum(self.variances) / s
1487
1488        fig = plt.figure()
1489        ax = fig.add_subplot(111)
1490
1491        x_axis = range(len(self.variances))
1492        plt.grid(True)
1493
1494        ax.set_xlabel('Principal Component Number')
1495        ax.set_ylabel('Proportion of Variance')
1496        ax.set_title(title + "\n")
1497        ax.plot(x_axis, vc, color="red")
1498        ax.scatter(x_axis, vc, color="red", label="Variance")
1499
1500        ax.plot(x_axis, cs, color="orange")
1501        ax.scatter(x_axis, cs, color="orange", label="Cumulative Variance")
1502        ax.legend(loc=0)
1503
1504        ax.axis([-0.5, len(self.variances) - 0.5, 0, 1])
1505
1506        if filename:
1507            plt.savefig(filename)
1508        else:
1509            plt.show()
1510
1511    def biplot(self, filename=None, components=(0, 1), title='Biplot'):
1512        """
1513        Draw biplot for PCA. Actual projection must be performed via pca(data)
1514        before bipot can be used.
1515
1516        :param filename: Name of the file to which the plot will be saved.
1517                         If None, plot will be displayed instead.
1518        :type filename: str
1519        :param components: List of two components to plot.
1520        :type components: list
1521        :param title: Plot title
1522        :type title: str
1523        """
1524        import pylab as plt
1525
1526        if len(components) < 2:
1527            raise ValueError, 'Two components are needed for biplot'
1528
1529        if not (0 <= min(components) <= max(components) < len(self.variances)):
1530            raise ValueError, 'Invalid components'
1531
1532        X = self.A[:, components[0]]
1533        Y = self.A[:, components[1]]
1534
1535        vectorsX = self.variances[:, components[0]]
1536        vectorsY = self.variances[:, components[1]]
1537
1538        max_load_value = self.variances.max() * 1.5
1539
1540        fig = plt.figure()
1541        ax1 = fig.add_subplot(111)
1542        ax1.set_title(title + "\n")
1543        ax1.set_xlabel("PC%s (%d%%)" % (components[0], self.variances[components[0]] / self.variance_sum * 100))
1544        ax1.set_ylabel("PC%s (%d%%)" % (components[1], self.variances[components[1]] / self.variance_sum * 100))
1545        ax1.xaxis.set_label_position('bottom')
1546        ax1.xaxis.set_ticks_position('bottom')
1547        ax1.yaxis.set_label_position('left')
1548        ax1.yaxis.set_ticks_position('left')
1549
1550        ax1.plot(X, Y, Colors[0])
1551
1552
1553        #eliminate double axis on right
1554        ax0 = ax1.twinx()
1555        ax0.yaxis.set_visible(False)
1556
1557        ax2 = ax0.twiny()
1558        ax2.xaxis.set_label_position('top')
1559        ax2.xaxis.set_ticks_position('top')
1560        ax2.yaxis.set_label_position('right')
1561        ax2.yaxis.set_ticks_position('right')
1562        for tl in ax2.get_xticklabels():
1563            tl.set_color('r')
1564        for tl in ax2.get_yticklabels():
1565            tl.set_color('r')
1566
1567        arrowprops = dict(facecolor='red', edgecolor='red', width=1, headwidth=4)
1568
1569        for (x, y, a) in zip(vectorsX, vectorsY, self.input_domain.attributes):
1570            if max(x, y) < 0.1:
1571                continue
1572            print x, y, a
1573            ax2.annotate('', (x, y), (0, 0), arrowprops=arrowprops)
1574            ax2.text(x * 1.1, y * 1.2, a.name, color='red')
1575
1576        ax2.set_xlim(-max_load_value, max_load_value)
1577        ax2.set_ylim(-max_load_value, max_load_value)
1578
1579        if filename:
1580            plt.savefig(filename)
1581        else:
1582            plt.show()
1583
1584
1585class Fda(object):
1586    """
1587    Construct a linear projection of data using FDA. When using this projection optimization method, data is always
1588    standardized prior to being projected.
1589
1590    If data instances are provided to the constructor,
1591    the optimization algorithm is called and the resulting projector
1592    (:class:`~Orange.projection.linear.FdaProjector`) is
1593    returned instead of the optimizer (instance of this class).
1594
1595    :rtype: :class:`~Orange.projection.linear.Fda` or
1596            :class:`~Orange.projection.linear.FdaProjector`
1597    """
1598
1599    def __new__(cls, dataset=None):
1600        self = object.__new__(cls)
1601        if dataset:
1602            self.__init__()
1603            return self.__call__(dataset)
1604        else:
1605            return self
1606
1607    def __call__(self, dataset):
1608        """
1609        Perform a FDA analysis on a data set and return a linear projector
1610        that maps data into another vector space.
1611
1612        :param dataset: input data set.
1613        :type dataset: :class:`Orange.data.Table`
1614
1615        :rtype: :class:`~Orange.projection.linear.FdaProjector`
1616        """
1617        X, Y = dataset.to_numpy_MA("a/c")
1618
1619        Xm = numpy.mean(X, axis=0)
1620        X -= Xm
1621
1622        #take care of the constant features
1623        stdev = numpy.std(X, axis=0)
1624        relevant_features = stdev != 0
1625        stdev[stdev == 0] = 1.
1626        X /= stdev
1627        X = X[:, relevant_features]
1628
1629        instances, features = X.shape
1630        class_count = len(set(Y))
1631        # special case when we have two classes
1632        if class_count == 2:
1633            data1 = MA.take(X, numpy.argwhere(Y == 0).flatten(), axis=0)
1634            data2 = MA.take(X, numpy.argwhere(Y != 0).flatten(), axis=0)
1635            miDiff = MA.average(data1, axis=1) - MA.average(data2, axis=1)
1636            covMatrix = (MA.dot(data1.T, data1) + MA.dot(data2.T, data2)) / instances
1637            U = numpy.linalg.inv(covMatrix) * miDiff
1638            D = numpy.array([1])
1639        else:
1640            # compute means and average covariances of examples in each class group
1641            Sw = MA.zeros([features, features])
1642            for v in set(Y):
1643                d = MA.take(X, numpy.argwhere(Y == v).flatten(), axis=0)
1644                d -= numpy.mean(d, axis=0)
1645                Sw += MA.dot(d.T, d)
1646            Sw /= instances
1647            total = MA.dot(X.T, X) / float(instances)
1648            Sb = total - Sw
1649
1650            matrix = numpy.linalg.inv(Sw) * Sb
1651            D, U = numpy.linalg.eigh(matrix)
1652
1653        sorted_indices = [i for _, i in sorted([(ev, i)
1654        for i, ev in enumerate(D)], reverse=True)]
1655        U = numpy.take(U, sorted_indices, axis=1)
1656        D = numpy.take(D, sorted_indices)
1657
1658        #insert zeros for constant features
1659        n, m = U.shape
1660        if m != M:
1661            U_ = numpy.zeros((n, M))
1662            U_[:, relevant_features] = U
1663            U = U_
1664
1665        out_domain = data.Domain([feature.Continuous("Comp.%d" %
1666                                                     (i + 1)) for
1667                                  i in range(len(D))], False)
1668
1669        return FdaProjector(input_domain=dataset.domain,
1670                            output_domain=out_domain,
1671                            center=Xm,
1672                            scale=stdev,
1673                            standardize=True,
1674                            eigen_vectors=U,
1675                            projection=U,
1676                            eigen_values=D)
1677
1678
1679class FdaProjector(Projector):
1680    """
1681    .. attribute:: eigen_vectors
1682
1683        Synonymous for :obj:`~Orange.projection.linear.Projector.projection`.
1684
1685    .. attribute:: eigen_values
1686
1687        Array containing eigenvalues corresponding to eigenvectors.
1688
1689    """
1690
1691
1692@deprecated_keywords({"dataMatrix": "data_matrix",
1693                      "classArray": "class_array",
1694                      "NComps": "ncomps",
1695                      "useGeneralizedEigenvectors": "use_generalized_eigenvectors"})
1696def create_pca_projection(data_matrix, class_array=None, ncomps=-1, use_generalized_eigenvectors=True):
1697    import warnings
1698
1699    warnings.warn("Deprecated in favour of Orange"
1700                  ".projection.linear.Pca.",
1701                  DeprecationWarning)
1702    if isinstance(data_matrix, numpy.ma.core.MaskedArray):
1703        data_matrix = numpy.array(data_matrix)
1704    if isinstance(class_array, numpy.ma.core.MaskedArray):
1705        class_array = numpy.array(class_array)
1706
1707    data_matrix = numpy.transpose(data_matrix)
1708
1709    s = numpy.sum(data_matrix, axis=0) / float(len(data_matrix))
1710    data_matrix -= s       # substract average value to get zero mean
1711
1712    if class_array is not None and use_generalized_eigenvectors:
1713        covarMatrix = numpy.dot(numpy.transpose(data_matrix), data_matrix)
1714        try:
1715            matrix = numpy.linalg.inv(covarMatrix)
1716        except numpy.linalg.LinAlgError:
1717            return None, None
1718        matrix = numpy.dot(matrix, numpy.transpose(data_matrix))
1719    else:
1720        matrix = numpy.transpose(data_matrix)
1721
1722    # compute dataMatrixT * L * dataMatrix
1723    if class_array is not None:
1724        # define the Laplacian matrix
1725        l = numpy.zeros((len(data_matrix), len(data_matrix)))
1726        for i in range(len(data_matrix)):
1727            for j in range(i + 1, len(data_matrix)):
1728                l[i, j] = -int(class_array[i] != class_array[j])
1729                l[j, i] = -int(class_array[i] != class_array[j])
1730
1731        s = numpy.sum(l, axis=0)      # doesn't matter which axis since the matrix l is symmetrical
1732        for i in range(len(data_matrix)):
1733            l[i, i] = -s[i]
1734
1735        matrix = numpy.dot(matrix, l)
1736
1737    matrix = numpy.dot(matrix, data_matrix)
1738
1739    vals, vectors = numpy.linalg.eig(matrix)
1740    if vals.dtype.kind == "c":       # if eigenvalues are complex numbers then do nothing
1741        return None, None
1742    vals = list(vals)
1743
1744    if ncomps == -1:
1745        ncomps = len(vals)
1746    ncomps = min(ncomps, len(vals))
1747
1748    ret_vals = []
1749    ret_indices = []
1750    for i in range(ncomps):
1751        ret_vals.append(max(vals))
1752        bestind = vals.index(max(vals))
1753        ret_indices.append(bestind)
1754        vals[bestind] = -1
1755
1756    return ret_vals, numpy.take(vectors.T, ret_indices,
1757                                axis=0)         # i-th eigenvector is the i-th column in vectors so we have to transpose the array
1758
1759createPCAProjection = create_pca_projection
Note: See TracBrowser for help on using the repository browser.