source: orange/Orange/evaluation/scoring.py @ 10422:fa72dfc419f9

Revision 10422:fa72dfc419f9, 98.3 KB checked in by anzeh <anze.staric@…>, 2 years ago (diff)

Emmit a warning when AUC could not be computed.

Line 
1import math
2import functools
3from operator import add
4
5import numpy
6
7import Orange
8from Orange import statc, corn
9from Orange.misc import deprecated_keywords, deprecated_function_name, \
10    deprecation_warning, environ
11from Orange.evaluation import testing
12
13try:
14    import matplotlib
15    HAS_MATPLOTLIB = True
16except ImportError:
17    matplotlib = None
18    HAS_MATPLOTLIB = False
19
20#### Private stuff
21
22def log2(x):
23    """Calculate logarithm in base 2."""
24    return math.log(x) / math.log(2)
25
26def check_non_zero(x):
27    """Throw Value Error when x = 0."""
28    if x == 0.:
29        raise ValueError, "Cannot compute the score: no examples or sum of weights is 0."
30
31def gettotweight(res):
32    """Sum all the weights."""
33    totweight = reduce(lambda x, y: x + y.weight, res.results, 0)
34    if totweight == 0.:
35        raise ValueError, "Cannot compute the score: sum of weights is 0."
36    return totweight
37
38def gettotsize(res):
39    """Get number of result instances."""
40    if len(res.results):
41        return len(res.results)
42    else:
43        raise ValueError, "Cannot compute the score: no examples."
44
45
46def split_by_iterations(res):
47    """Split ExperimentResults of a multiple iteratation test into a list
48    of ExperimentResults, one for each iteration.
49    """
50    if res.number_of_iterations < 2:
51        return [res]
52       
53    ress = [Orange.evaluation.testing.ExperimentResults(
54                1, res.classifier_names, res.class_values,
55                res.weights, classifiers=res.classifiers,
56                loaded=res.loaded, test_type=res.test_type, labels=res.labels)
57            for _ in range(res.number_of_iterations)]
58    for te in res.results:
59        ress[te.iteration_number].results.append(te)
60    return ress
61
62def split_by_classifiers(res):
63    """Split an instance of :obj:`ExperimentResults` into a list of
64    :obj:`ExperimentResults`, one for each classifier.
65    """
66    split_res = []
67    for i in range(len(res.classifierNames)):
68        r = Orange.evaluation.testing.ExperimentResults(res.numberOfIterations,
69                    [res.classifierNames[i]], res.classValues,
70                    weights=res.weights, baseClass=res.baseClass,
71                    classifiers=[res.classifiers[i]] if res.classifiers else [],
72                    test_type=res.test_type, labels=res.labels)
73        r.results = []
74        for te in res.results:
75            r.results.append(Orange.evaluation.testing.TestedExample(
76                te.iterationNumber, te.actualClass, n=1, weight=te.weight))
77            r.results[-1].classes = [te.classes[i]]
78            r.results[-1].probabilities = [te.probabilities[i]]
79        split_res.append(r)
80    return split_res
81
82
83def class_probabilities_from_res(res, **argkw):
84    """Calculate class probabilities."""
85    probs = [0.] * len(res.class_values)
86    if argkw.get("unweighted", 0) or not res.weights:
87        for tex in res.results:
88            probs[int(tex.actual_class)] += 1.
89        totweight = gettotsize(res)
90    else:
91        totweight = 0.
92        for tex in res.results:
93            probs[tex.actual_class] += tex.weight
94            totweight += tex.weight
95        check_non_zero(totweight)
96    return [prob / totweight for prob in probs]
97
98
99@deprecated_keywords({
100    "foldN": "fold_n",
101    "reportSE": "report_se",
102    "iterationIsOuter": "iteration_is_outer"})
103def statistics_by_folds(stats, fold_n, report_se, iteration_is_outer):
104    # remove empty folds, turn the matrix so that learner is outer
105    if iteration_is_outer:
106        if not stats:
107            raise ValueError, "Cannot compute the score: no examples or sum of weights is 0."
108        number_of_learners = len(stats[0])
109        stats = filter(lambda (x, fN): fN > 0, zip(stats, fold_n))
110        stats = [[x[lrn] / fN for x, fN in stats]
111                 for lrn in range(number_of_learners)]
112    else:
113        stats = [[x / Fn for x, Fn in filter(lambda (x, Fn): Fn > 0,
114                 zip(lrnD, fold_n))] for lrnD in stats]
115
116    if not stats:
117        raise ValueError, "Cannot compute the score: no classifiers"
118    if not stats[0]:
119        raise ValueError, "Cannot compute the score: no examples or sum of weights is 0."
120   
121    if report_se:
122        return [(statc.mean(x), statc.sterr(x)) for x in stats]
123    else:
124        return [statc.mean(x) for x in stats]
125   
126def ME(res, **argkw):
127    MEs = [0.] * res.number_of_learners
128
129    if argkw.get("unweighted", 0) or not res.weights:
130        for tex in res.results:
131            MEs = map(lambda res, cls, ac = float(tex.actual_class):
132                      res + abs(float(cls) - ac), MEs, tex.classes)
133        totweight = gettotsize(res)
134    else:
135        for tex in res.results:
136            MEs = map(lambda res, cls, ac=float(tex.actual_class), tw=tex.weight:
137                       res + tw * abs(float(cls) - ac), MEs, tex.classes)
138        totweight = gettotweight(res)
139
140    return [x / totweight for x in MEs]
141
142MAE = ME
143
144
145class ConfusionMatrix:
146    """
147    Classification result summary.
148    """
149    #: True Positive predictions
150    TP = 0.
151    #:True Negative predictions
152    TN = 0.
153    #:False Positive predictions
154    FP = 0.
155    #: False Negative predictions
156    FN = 0.
157
158    @deprecated_keywords({"predictedPositive": "predicted_positive",
159                          "isPositive": "is_positive"})
160    def addTFPosNeg(self, predicted_positive, is_positive, weight=1.0):
161        """
162        Update confusion matrix with result of a single classification.
163
164        :param predicted_positive: positive class value was predicted
165        :param is_positive: correct class value is positive
166        :param weight: weight of the selected instance
167         """
168        if predicted_positive:
169            if is_positive:
170                self.TP += weight
171            else:
172                self.FP += weight
173        else:
174            if is_positive:
175                self.FN += weight
176            else:
177                self.TN += weight
178
179
180
181#########################################################################
182# PERFORMANCE MEASURES:
183# Scores for evaluation of numeric predictions
184
185def check_argkw(dct, lst):
186    """Return True if any item from lst has a non-zero value in dct."""
187    return reduce(lambda x, y: x or y, [dct.get(k, 0) for k in lst])
188
189def regression_error(res, **argkw):
190    """Return the regression error (default: MSE)."""
191    if argkw.get("SE", 0) and res.number_of_iterations > 1:
192        # computes the scores for each iteration, then averages
193        scores = [[0.] * res.number_of_iterations
194                  for _ in range(res.number_of_learners)]
195        norm = None
196        if argkw.get("norm-abs", 0) or argkw.get("norm-sqr", 0):
197            norm = [0.] * res.number_of_iterations
198
199        # counts examples in each iteration
200        nIter = [0] * res.number_of_iterations
201        # average class in each iteration
202        a = [0] * res.number_of_iterations
203        for tex in res.results:
204            nIter[tex.iteration_number] += 1
205            a[tex.iteration_number] += float(tex.actual_class)
206        a = [a[i] / nIter[i] for i in range(res.number_of_iterations)]
207
208        if argkw.get("unweighted", 0) or not res.weights:
209            # iterate accross test cases
210            for tex in res.results:
211                ai = float(tex.actual_class)
212                nIter[tex.iteration_number] += 1
213
214                # compute normalization, if required
215                if argkw.get("norm-abs", 0):
216                    norm[tex.iteration_number] += abs(ai - a[tex.iteration_number])
217                elif argkw.get("norm-sqr", 0):
218                    norm[tex.iteration_number] += (ai - a[tex.iteration_number])**2
219
220                # iterate accross results of different regressors
221                for i, cls in enumerate(tex.classes):
222                    if argkw.get("abs", 0):
223                        scores[i][tex.iteration_number] += abs(float(cls) - ai)
224                    else:
225                        scores[i][tex.iteration_number] += (float(cls) - ai)**2
226        else: # unweighted != 0
227            raise NotImplementedError, "weighted error scores with SE not implemented yet"
228
229        if argkw.get("norm-abs") or argkw.get("norm-sqr"):
230            scores = [[x / n for x, n in zip(y, norm)] for y in scores]
231        else:
232            scores = [[x / ni for x, ni in zip(y, nIter)] for y in scores]
233
234        if argkw.get("R2"):
235            scores = [[1.0 - x for x in y] for y in scores]
236
237        if argkw.get("sqrt", 0):
238            scores = [[math.sqrt(x) for x in y] for y in scores]
239
240        return [(statc.mean(x), statc.std(x)) for x in scores]
241       
242    else: # single iteration (testing on a single test set)
243        scores = [0.] * res.number_of_learners
244        norm = 0.
245
246        if argkw.get("unweighted", 0) or not res.weights:
247            a = sum([tex.actual_class for tex in res.results]) \
248                / len(res.results)
249            for tex in res.results:
250                if argkw.get("abs", 0):
251                    scores = map(lambda res, cls, ac=float(tex.actual_class):
252                                 res + abs(float(cls) - ac), scores, tex.classes)
253                else:
254                    scores = map(lambda res, cls, ac=float(tex.actual_class):
255                                 res + (float(cls) - ac)**2, scores, tex.classes)
256
257                if argkw.get("norm-abs", 0):
258                    norm += abs(tex.actual_class - a)
259                elif argkw.get("norm-sqr", 0):
260                    norm += (tex.actual_class - a)**2
261            totweight = gettotsize(res)
262        else:
263            # UNFINISHED
264            MSEs = [0.] * res.number_of_learners
265            for tex in res.results:
266                MSEs = map(lambda res, cls, ac=float(tex.actual_class),
267                           tw=tex.weight:
268                           res + tw * (float(cls) - ac)**2, MSEs, tex.classes)
269            totweight = gettotweight(res)
270
271        if argkw.get("norm-abs", 0) or argkw.get("norm-sqr", 0):
272            scores = [s / norm for s in scores]
273        else: # normalize by number of instances (or sum of weights)
274            scores = [s / totweight for s in scores]
275
276        if argkw.get("R2"):
277            scores = [1. - s for s in scores]
278
279        if argkw.get("sqrt", 0):
280            scores = [math.sqrt(x) for x in scores]
281
282        return scores
283
284def MSE(res, **argkw):
285    """Compute mean-squared error."""
286    return regression_error(res, **argkw)
287   
288def RMSE(res, **argkw):
289    """Compute root mean-squared error."""
290    argkw.setdefault("sqrt", True)
291    return regression_error(res, **argkw)
292
293def MAE(res, **argkw):
294    """Compute mean absolute error."""
295    argkw.setdefault("abs", True)
296    return regression_error(res, **argkw)
297
298def RSE(res, **argkw):
299    """Compute relative squared error."""
300    argkw.setdefault("norm-sqr", True)
301    return regression_error(res, **argkw)
302
303def RRSE(res, **argkw):
304    """Compute relative squared error."""
305    argkw.setdefault("norm-sqr", True)
306    argkw.setdefault("sqrt", True)
307    return regression_error(res, **argkw)
308
309def RAE(res, **argkw):
310    """Compute relative absolute error."""
311    argkw.setdefault("abs", True)
312    argkw.setdefault("norm-abs", True)
313    return regression_error(res, **argkw)
314
315def R2(res, **argkw):
316    """Compute the coefficient of determination, R-squared."""
317    argkw.setdefault("norm-sqr", True)
318    argkw.setdefault("R2", True)
319    return regression_error(res, **argkw)
320
321def MSE_old(res, **argkw):
322    """Compute mean-squared error."""
323    if argkw.get("SE", 0) and res.number_of_iterations > 1:
324        MSEs = [[0.] * res.number_of_iterations
325                for _ in range(res.number_of_learners)]
326        nIter = [0] * res.number_of_iterations
327        if argkw.get("unweighted", 0) or not res.weights:
328            for tex in res.results:
329                ac = float(tex.actual_class)
330                nIter[tex.iteration_number] += 1
331                for i, cls in enumerate(tex.classes):
332                    MSEs[i][tex.iteration_number] += (float(cls) - ac)**2
333        else:
334            raise ValueError, "weighted RMSE with SE not implemented yet"
335        MSEs = [[x / ni for x, ni in zip(y, nIter)] for y in MSEs]
336        if argkw.get("sqrt", 0):
337            MSEs = [[math.sqrt(x) for x in y] for y in MSEs]
338        return [(statc.mean(x), statc.std(x)) for x in MSEs]
339       
340    else:
341        MSEs = [0.] * res.number_of_learners
342        if argkw.get("unweighted", 0) or not res.weights:
343            for tex in res.results:
344                MSEs = map(lambda res, cls, ac=float(tex.actual_class):
345                           res + (float(cls) - ac)**2, MSEs, tex.classes)
346            totweight = gettotsize(res)
347        else:
348            for tex in res.results:
349                MSEs = map(lambda res, cls, ac=float(tex.actual_class),
350                           tw=tex.weight: res + tw * (float(cls) - ac)**2,
351                           MSEs, tex.classes)
352            totweight = gettotweight(res)
353
354        if argkw.get("sqrt", 0):
355            MSEs = [math.sqrt(x) for x in MSEs]
356        return [x / totweight for x in MSEs]
357
358def RMSE_old(res, **argkw):
359    """Compute root mean-squared error."""
360    argkw.setdefault("sqrt", 1)
361    return MSE_old(res, **argkw)
362
363#########################################################################
364# PERFORMANCE MEASURES:
365# Scores for evaluation of classifiers
366
367class CA(list):
368    """
369    Compute percentage of matches between predicted and actual class.
370
371    :param test_results: :obj:`~Orange.evaluation.testing.ExperimentResults`
372                         or list of :obj:`ConfusionMatrix`.
373    :param report_se: include standard error in result.
374    :param ignore_weights: ignore instance weights.
375    :rtype: list of scores, one for each learner.
376
377    Standard errors are estimated from deviation of CAs across folds (if
378    test_results were produced by cross_validation) or approximated under
379    the assumption of normal distribution otherwise.
380    """
381    CONFUSION_MATRIX = 0
382    CONFUSION_MATRIX_LIST = 1
383    CLASSIFICATION = 2
384    CROSS_VALIDATION = 3
385
386    @deprecated_keywords({"reportSE": "report_se",
387                          "unweighted": "ignore_weights"})
388    def __init__(self, test_results, report_se = False, ignore_weights=False):
389        super(CA, self).__init__()
390        self.report_se = report_se
391        self.ignore_weights = ignore_weights
392
393        input_type = self.get_input_type(test_results)
394        if input_type == self.CONFUSION_MATRIX:
395            self[:] = [self.from_confusion_matrix(test_results)]
396        elif input_type == self.CONFUSION_MATRIX_LIST:
397            self[:] = self.from_confusion_matrix_list(test_results)
398        elif input_type == self.CLASSIFICATION:
399            self[:] = self.from_classification_results(test_results)
400        elif input_type == self.CROSS_VALIDATION:
401            self[:] = self.from_crossvalidation_results(test_results)
402
403    def from_confusion_matrix(self, cm):
404        all_predictions = 0.
405        correct_predictions = 0.
406        if isinstance(cm, ConfusionMatrix):
407            all_predictions += cm.TP + cm.FN + cm.FP + cm.TN
408            correct_predictions += cm.TP + cm.TN
409        else:
410            for r, row in enumerate(cm):
411                for c, column in enumerate(row):
412                    if r == c:
413                        correct_predictions += column
414                    all_predictions += column
415
416        check_non_zero(all_predictions)
417        ca = correct_predictions / all_predictions
418
419        if self.report_se:
420            return ca, ca * (1 - ca) / math.sqrt(all_predictions)
421        else:
422            return ca
423
424    def from_confusion_matrix_list(self, confusion_matrices):
425        return [self.from_confusion_matrix(cm) for cm in confusion_matrices]
426
427    def from_classification_results(self, test_results):
428        CAs = [0.] * test_results.number_of_learners
429        totweight = 0.
430        for tex in test_results.results:
431            w = 1. if self.ignore_weights else tex.weight
432            CAs = map(lambda res, cls: res + (cls == tex.actual_class and w),
433                      CAs, tex.classes)
434            totweight += w
435        check_non_zero(totweight)
436        ca = [x / totweight for x in CAs]
437
438        if self.report_se:
439            return [(x, x * (1 - x) / math.sqrt(totweight)) for x in ca]
440        else:
441            return ca
442
443    def from_crossvalidation_results(self, test_results):
444        CAsByFold = [[0.] * test_results.number_of_iterations
445                     for _ in range(test_results.number_of_learners)]
446        foldN = [0.] * test_results.number_of_iterations
447
448        for tex in test_results.results:
449            w = 1. if self.ignore_weights else tex.weight
450            for lrn in range(test_results.number_of_learners):
451                CAsByFold[lrn][tex.iteration_number] += (tex.classes[lrn] == 
452                    tex.actual_class) and w
453            foldN[tex.iteration_number] += w
454
455        return statistics_by_folds(CAsByFold, foldN, self.report_se, False)
456
457    def get_input_type(self, test_results):
458        if isinstance(test_results, ConfusionMatrix):
459            return self.CONFUSION_MATRIX
460        elif isinstance(test_results, testing.ExperimentResults):
461            if test_results.number_of_iterations == 1:
462                return self.CLASSIFICATION
463            else:
464                return self.CROSS_VALIDATION
465        elif isinstance(test_results, list):
466            return self.CONFUSION_MATRIX_LIST
467
468
469@deprecated_keywords({"reportSE": "report_se",
470                      "unweighted": "ignore_weights"})
471def AP(res, report_se=False, ignore_weights=False, **argkw):
472    """Compute the average probability assigned to the correct class."""
473    if res.number_of_iterations == 1:
474        APs=[0.] * res.number_of_learners
475        if ignore_weights or not res.weights:
476            for tex in res.results:
477                APs = map(lambda res, probs: res + probs[tex.actual_class],
478                          APs, tex.probabilities)
479            totweight = gettotsize(res)
480        else:
481            totweight = 0.
482            for tex in res.results:
483                APs = map(lambda res, probs: res + probs[tex.actual_class] *
484                          tex.weight, APs, tex.probabilities)
485                totweight += tex.weight
486        check_non_zero(totweight)
487        return [AP / totweight for AP in APs]
488
489    APsByFold = [[0.] * res.number_of_learners
490                 for _ in range(res.number_of_iterations)]
491    foldN = [0.] * res.number_of_iterations
492    if ignore_weights or not res.weights:
493        for tex in res.results:
494            APsByFold[tex.iteration_number] = map(lambda res, probs:
495                res + probs[tex.actual_class],
496                APsByFold[tex.iteration_number], tex.probabilities)
497            foldN[tex.iteration_number] += 1
498    else:
499        for tex in res.results:
500            APsByFold[tex.iteration_number] = map(lambda res, probs:
501                res + probs[tex.actual_class] * tex.weight,
502                APsByFold[tex.iteration_number], tex.probabilities)
503            foldN[tex.iteration_number] += tex.weight
504
505    return statistics_by_folds(APsByFold, foldN, report_se, True)
506
507
508@deprecated_keywords({"reportSE": "report_se",
509                      "unweighted": "ignore_weights"})
510def Brier_score(res, report_se=False, ignore_weights=False, **argkw):
511    """Compute the Brier score, defined as the average (over test instances)
512    of :math:`\sum_x(t(x) - p(x))^2`, where x is a class value, t(x) is 1 for
513    the actual class value and 0 otherwise, and p(x) is the predicted
514    probability of x.
515    """
516    # Computes an average (over examples) of sum_x(t(x) - p(x))^2, where
517    #    x is class,
518    #    t(x) is 0 for 'wrong' and 1 for 'correct' class
519    #    p(x) is predicted probabilty.
520    # There's a trick: since t(x) is zero for all classes but the
521    # correct one (c), we compute the sum as sum_x(p(x)^2) - 2*p(c) + 1
522    # Since +1 is there for each example, it adds 1 to the average
523    # We skip the +1 inside the sum and add it just at the end of the function
524    # We take max(result, 0) to avoid -0.0000x due to rounding errors
525
526    if res.number_of_iterations == 1:
527        MSEs=[0.] * res.number_of_learners
528        if ignore_weights or not res.weights:
529            totweight = 0.
530            for tex in res.results:
531                MSEs = map(lambda res, probs: res + reduce(
532                    lambda s, pi: s + pi**2, probs, 0) - 
533                    2 * probs[tex.actual_class], MSEs, tex.probabilities)
534                totweight += tex.weight
535        else:
536            for tex in res.results:
537                MSEs = map(lambda res, probs: res + tex.weight * reduce(
538                    lambda s, pi: s + pi**2, probs, 0) - 
539                    2 * probs[tex.actual_class], MSEs, tex.probabilities)
540            totweight = gettotweight(res)
541        check_non_zero(totweight)
542        if report_se:
543            ## change this, not zero!!!
544            return [(max(x / totweight + 1., 0), 0) for x in MSEs]
545        else:
546            return [max(x / totweight + 1., 0) for x in MSEs]
547
548    BSs = [[0.] * res.number_of_learners
549           for _ in range(res.number_of_iterations)]
550    foldN = [0.] * res.number_of_iterations
551
552    if ignore_weights or not res.weights:
553        for tex in res.results:
554            BSs[tex.iteration_number] = map(lambda rr, probs: rr + reduce(
555                lambda s, pi: s + pi**2, probs, 0) -
556                2 * probs[tex.actual_class], BSs[tex.iteration_number],
557                tex.probabilities)
558            foldN[tex.iteration_number] += 1
559    else:
560        for tex in res.results:
561            BSs[tex.iteration_number] = map(lambda res, probs:
562                res + tex.weight * reduce(lambda s, pi: s + pi**2, probs, 0) -
563                2 * probs[tex. actual_class], BSs[tex.iteration_number],
564                tex.probabilities)
565            foldN[tex.iteration_number] += tex.weight
566
567    stats = statistics_by_folds(BSs, foldN, report_se, True)
568    if report_se:
569        return [(x + 1., y) for x, y in stats]
570    else:
571        return [x + 1. for x in stats]
572
573def BSS(res, **argkw):
574    return [1 - x / 2 for x in apply(Brier_score, (res, ), argkw)]
575
576def IS_ex(Pc, P):
577    """Pc aposterior probability, P aprior"""
578    if Pc >= P:
579        return -log2(P) + log2(Pc)
580    else:
581        return -(-log2(1 - P) + log2(1 - Pc))
582
583
584@deprecated_keywords({"reportSE": "report_se"})
585def IS(res, apriori=None, report_se=False, **argkw):
586    """Compute the information score as defined by
587    `Kononenko and Bratko (1991) \
588    <http://www.springerlink.com/content/g5p7473160476612/>`_.
589    Argument :obj:`apriori` gives the apriori class
590    distribution; if it is omitted, the class distribution is computed from
591    the actual classes of examples in :obj:`res`.
592    """
593    if not apriori:
594        apriori = class_probabilities_from_res(res)
595
596    if res.number_of_iterations == 1:
597        ISs = [0.] * res.number_of_learners
598        if argkw.get("unweighted", 0) or not res.weights:
599            for tex in res.results:
600              for i in range(len(tex.probabilities)):
601                    cls = tex.actual_class
602                    ISs[i] += IS_ex(tex.probabilities[i][cls], apriori[cls])
603            totweight = gettotsize(res)
604        else:
605            for tex in res.results:
606              for i in range(len(tex.probabilities)):
607                    cls = tex.actual_class
608                    ISs[i] += (IS_ex(tex.probabilities[i][cls], apriori[cls]) *
609                               tex.weight)
610            totweight = gettotweight(res)
611        if report_se:
612            return [(IS / totweight, 0) for IS in ISs]
613        else:
614            return [IS / totweight for IS in ISs]
615
616       
617    ISs = [[0.] * res.number_of_iterations
618           for _ in range(res.number_of_learners)]
619    foldN = [0.] * res.number_of_iterations
620
621    # compute info scores for each fold   
622    if argkw.get("unweighted", 0) or not res.weights:
623        for tex in res.results:
624            for i in range(len(tex.probabilities)):
625                cls = tex.actual_class
626                ISs[i][tex.iteration_number] += IS_ex(tex.probabilities[i][cls],
627                                                apriori[cls])
628            foldN[tex.iteration_number] += 1
629    else:
630        for tex in res.results:
631            for i in range(len(tex.probabilities)):
632                cls = tex.actual_class
633                ISs[i][tex.iteration_number] += IS_ex(tex.probabilities[i][cls],
634                                                apriori[cls]) * tex.weight
635            foldN[tex.iteration_number] += tex.weight
636
637    return statistics_by_folds(ISs, foldN, report_se, False)
638   
639
640def Wilcoxon(res, statistics, **argkw):
641    res1, res2 = [], []
642    for ri in split_by_iterations(res):
643        stats = statistics(ri, **argkw)
644        if len(stats) != 2:
645            raise TypeError, "Wilcoxon compares two classifiers, no more, no less"
646        res1.append(stats[0])
647        res2.append(stats[1])
648    return statc.wilcoxont(res1, res2)
649
650def rank_difference(res, statistics, **argkw):
651    if not res.results:
652        raise TypeError, "no experiments"
653
654    k = len(res.results[0].classes)
655    if k < 2:
656        raise TypeError, "nothing to compare (less than two classifiers given)"
657    if k == 2:
658        return apply(Wilcoxon, (res, statistics), argkw)
659    else:
660        return apply(Friedman, (res, statistics), argkw)
661
662
663@deprecated_keywords({"res": "test_results",
664                      "classIndex": "class_index",
665                      "unweighted": "ignore_weights"})
666def confusion_matrices(test_results, class_index=-1,
667                       ignore_weights=False, cutoff=.5):
668    """
669    Return confusion matrices for test_results.
670
671    :param test_results: test results
672    :param class_index: index of class value for which the confusion matrices
673                        are to be computed.
674    :param ignore_weights: ignore instance weights.
675    :params cutoff: cutoff for probability
676
677    :rtype: list of :obj:`ConfusionMatrix`
678    """
679    tfpns = [ConfusionMatrix() for _ in range(test_results.number_of_learners)]
680   
681    if class_index < 0:
682        numberOfClasses = len(test_results.class_values)
683        if class_index < -1 or numberOfClasses > 2:
684            cm = [[[0.] * numberOfClasses for _ in range(numberOfClasses)]
685                  for _ in range(test_results.number_of_learners)]
686            if ignore_weights or not test_results.weights:
687                for tex in test_results.results:
688                    trueClass = int(tex.actual_class)
689                    for li, pred in enumerate(tex.classes):
690                        predClass = int(pred)
691                        if predClass < numberOfClasses:
692                            cm[li][trueClass][predClass] += 1
693            else:
694                for tex in test_results.results:
695                    trueClass = int(tex.actual_class)
696                    for li, pred in tex.classes:
697                        predClass = int(pred)
698                        if predClass < numberOfClasses:
699                            cm[li][trueClass][predClass] += tex.weight
700            return cm
701           
702        elif test_results.baseClass >= 0:
703            class_index = test_results.baseClass
704        else:
705            class_index = 1
706
707    if cutoff != .5:
708        if ignore_weights or not test_results.weights:
709            for lr in test_results.results:
710                isPositive=(lr.actual_class==class_index)
711                for i in range(test_results.number_of_learners):
712                    tfpns[i].addTFPosNeg(lr.probabilities[i][class_index] >
713                                         cutoff, isPositive)
714        else:
715            for lr in test_results.results:
716                isPositive=(lr.actual_class == class_index)
717                for i in range(test_results.number_of_learners):
718                    tfpns[i].addTFPosNeg(lr.probabilities[i][class_index] > 
719                                         cutoff, isPositive, lr.weight)
720    else:
721        if ignore_weights or not test_results.weights:
722            for lr in test_results.results:
723                isPositive = (lr.actual_class == class_index)
724                for i in range(test_results.number_of_learners):
725                    tfpns[i].addTFPosNeg(lr.classes[i] == class_index,
726                                         isPositive)
727        else:
728            for lr in test_results.results:
729                isPositive = (lr.actual_class == class_index)
730                for i in range(test_results.number_of_learners):
731                    tfpns[i].addTFPosNeg(lr.classes[i] == class_index,
732                                         isPositive, lr.weight)
733    return tfpns
734
735
736# obsolete (renamed)
737compute_confusion_matrices = confusion_matrices
738
739
740@deprecated_keywords({"confusionMatrix": "confusion_matrix"})
741def confusion_chi_square(confusion_matrix):
742    """
743    Return chi square statistic of the confusion matrix
744    (higher value indicates that prediction is not by chance).
745    """
746    if isinstance(confusion_matrix, ConfusionMatrix) or \
747       not isinstance(confusion_matrix[1], list):
748        return _confusion_chi_square(confusion_matrix)
749    else:
750        return map(_confusion_chi_square, confusion_matrix)
751
752def _confusion_chi_square(confusion_matrix):
753    if isinstance(confusion_matrix, ConfusionMatrix):
754        c = confusion_matrix
755        confusion_matrix = [[c.TP, c.FN], [c.FP, c.TN]]
756    dim = len(confusion_matrix)
757    rowPriors = [sum(r) for r in confusion_matrix]
758    colPriors = [sum(r[i] for r in confusion_matrix) for i in range(dim)]
759    total = sum(rowPriors)
760    rowPriors = [r / total for r in rowPriors]
761    colPriors = [r / total for r in colPriors]
762    ss = 0
763    for ri, row in enumerate(confusion_matrix):
764        for ci, o in enumerate(row):
765            e = total * rowPriors[ri] * colPriors[ci]
766            if not e:
767                return -1, -1, -1
768            ss += (o - e)**2 / e
769    df = (dim - 1)**2
770    return ss, df, statc.chisqprob(ss, df)
771
772@deprecated_keywords({"confm": "confusion_matrix"})
773def sens(confusion_matrix):
774    """
775    Return `sensitivity
776    <http://en.wikipedia.org/wiki/Sensitivity_and_specificity>`_ (proportion
777    of actual positives which are correctly identified as such).
778    """
779    if type(confusion_matrix) == list:
780        return [sens(cm) for cm in confusion_matrix]
781    else:
782        tot = confusion_matrix.TP+confusion_matrix.FN
783        if tot < 1e-6:
784            import warnings
785            warnings.warn("Can't compute sensitivity: one or both classes have no instances")
786            return -1
787
788        return confusion_matrix.TP / tot
789
790
791@deprecated_keywords({"confm": "confusion_matrix"})
792def recall(confusion_matrix):
793    """
794    Return `recall <http://en.wikipedia.org/wiki/Precision_and_recall>`_
795    (fraction of relevant instances that are retrieved).
796    """
797    return sens(confusion_matrix)
798
799
800@deprecated_keywords({"confm": "confusion_matrix"})
801def spec(confusion_matrix):
802    """
803    Return `specificity
804    <http://en.wikipedia.org/wiki/Sensitivity_and_specificity>`_
805    (proportion of negatives which are correctly identified).
806    """
807    if type(confusion_matrix) == list:
808        return [spec(cm) for cm in confusion_matrix]
809    else:
810        tot = confusion_matrix.FP+confusion_matrix.TN
811        if tot < 1e-6:
812            import warnings
813            warnings.warn("Can't compute specificity: one or both classes have no instances")
814            return -1
815        return confusion_matrix.TN / tot
816
817
818@deprecated_keywords({"confm": "confusion_matrix"})
819def PPV(confusion_matrix):
820    """
821    Return `positive predictive value
822    <http://en.wikipedia.org/wiki/Positive_predictive_value>`_ (proportion of
823    subjects with positive test results who are correctly diagnosed)."""
824    if type(confusion_matrix) == list:
825        return [PPV(cm) for cm in confusion_matrix]
826    else:
827        tot = confusion_matrix.TP + confusion_matrix.FP
828        if tot < 1e-6:
829            import warnings
830            warnings.warn("Can't compute PPV: one or both classes have no instances")
831            return -1
832        return confusion_matrix.TP/tot
833
834
835@deprecated_keywords({"confm": "confusion_matrix"})
836def precision(confusion_matrix):
837    """
838    Return `precision <http://en.wikipedia.org/wiki/Precision_and_recall>`_
839    (retrieved instances that are relevant).
840    """
841    return PPV(confusion_matrix)
842
843@deprecated_keywords({"confm": "confusion_matrix"})
844def NPV(confusion_matrix):
845    """
846    Return `negative predictive value
847    <http://en.wikipedia.org/wiki/Negative_predictive_value>`_ (proportion of
848    subjects with a negative test result who are correctly diagnosed).
849     """
850    if type(confusion_matrix) == list:
851        return [NPV(cm) for cm in confusion_matrix]
852    else:
853        tot = confusion_matrix.FN + confusion_matrix.TN
854        if tot < 1e-6:
855            import warnings
856            warnings.warn("Can't compute NPV: one or both classes have no instances")
857            return -1
858        return confusion_matrix.TN / tot
859
860@deprecated_keywords({"confm": "confusion_matrix"})
861def F1(confusion_matrix):
862    """
863    Return `F1 score <http://en.wikipedia.org/wiki/F1_score>`_
864    (harmonic mean of precision and recall).
865    """
866    if type(confusion_matrix) == list:
867        return [F1(cm) for cm in confusion_matrix]
868    else:
869        p = precision(confusion_matrix)
870        r = recall(confusion_matrix)
871        if p + r > 0:
872            return 2. * p * r / (p + r)
873        else:
874            import warnings
875            warnings.warn("Can't compute F1: P + R is zero or not defined")
876            return -1
877
878
879@deprecated_keywords({"confm": "confusion_matrix"})
880def Falpha(confusion_matrix, alpha=1.0):
881    """
882    Return the alpha-mean of precision and recall over the given confusion
883    matrix.
884    """
885    if type(confusion_matrix) == list:
886        return [Falpha(cm, alpha=alpha) for cm in confusion_matrix]
887    else:
888        p = precision(confusion_matrix)
889        r = recall(confusion_matrix)
890        return (1. + alpha) * p * r / (alpha * p + r)
891
892
893@deprecated_keywords({"confm": "confusion_matrix"})
894def MCC(confusion_matrix):
895    """
896    Return `Matthew correlation coefficient
897    <http://en.wikipedia.org/wiki/Matthews_correlation_coefficient>`_
898    (correlation coefficient between the observed and predicted binary
899    classifications).
900    """
901    # code by Boris Gorelik
902    if type(confusion_matrix) == list:
903        return [MCC(cm) for cm in confusion_matrix]
904    else:
905        truePositive = confusion_matrix.TP
906        trueNegative = confusion_matrix.TN
907        falsePositive = confusion_matrix.FP
908        falseNegative = confusion_matrix.FN
909         
910        try:   
911            r = (((truePositive * trueNegative) -
912                  (falsePositive * falseNegative)) /
913                 math.sqrt((truePositive + falsePositive) *
914                           (truePositive + falseNegative) * 
915                           (trueNegative + falsePositive) * 
916                           (trueNegative + falseNegative)))
917        except ZeroDivisionError:
918            # Zero difision occurs when there is either no true positives
919            # or no true negatives i.e. the problem contains only one
920            # type of classes.
921            import warnings
922            warnings.warn("Can't compute MCC: TP or TN is zero or not defined")
923            r = None
924
925    return r
926
927
928@deprecated_keywords({"bIsListOfMatrices": "b_is_list_of_matrices"})
929def scotts_pi(confusion_matrix, b_is_list_of_matrices=True):
930   """Compute Scott's Pi for measuring inter-rater agreement for nominal data
931
932   http://en.wikipedia.org/wiki/Scott%27s_Pi
933   Scott's Pi is a statistic for measuring inter-rater reliability for nominal
934   raters.
935
936   @param confusion_matrix: confusion matrix, or list of confusion matrices.
937                            To obtain non-binary confusion matrix, call
938                            Orange.evaluation.scoring.confusion_matrices
939                            and set the class_index parameter to -2.
940   @param b_is_list_of_matrices: specifies whether confm is list of matrices.
941                            This function needs to operate on non-binary
942                            confusion matrices, which are represented by python
943                            lists, therefore one needs a way to distinguish
944                            between a single matrix and list of matrices
945   """
946
947   if b_is_list_of_matrices:
948       try:
949           return [scotts_pi(cm, b_is_list_of_matrices=False)
950                   for cm in confusion_matrix]
951       except TypeError:
952           # Nevermind the parameter, maybe this is a "conventional" binary
953           # confusion matrix and bIsListOfMatrices was specified by mistake
954           return scottsPiSingle(confusion_matrix, bIsListOfMatrices=False)
955   else:
956       if isinstance(confusion_matrix, ConfusionMatrix):
957           confusion_matrix = numpy.array([[confusion_matrix.TP,
958               confusion_matrix.FN], [confusion_matrix.FP,
959               confusion_matrix.TN]], dtype=float)
960       else:
961           confusion_matrix = numpy.array(confusion_matrix, dtype=float)
962
963       marginalSumOfRows = numpy.sum(confusion_matrix, axis=0)
964       marginalSumOfColumns = numpy.sum(confusion_matrix, axis=1)
965       jointProportion = (marginalSumOfColumns + marginalSumOfRows) / \
966                           (2. * numpy.sum(confusion_matrix))
967       # In the eq. above, 2. is what the Wikipedia page calls
968       # the number of annotators. Here we have two annotators:
969       # the observed (true) labels (annotations) and the predicted by
970       # the learners.
971
972       prExpected = numpy.sum(jointProportion ** 2)
973       prActual = numpy.sum(numpy.diag(confusion_matrix)) / \
974                  numpy.sum(confusion_matrix)
975
976       ret = (prActual - prExpected) / (1.0 - prExpected)
977       return ret
978
979@deprecated_keywords({"classIndex": "class_index",
980                      "unweighted": "ignore_weights"})
981def AUCWilcoxon(res, class_index=-1, ignore_weights=False, **argkw):
982    """Compute the area under ROC (AUC) and its standard error using
983    Wilcoxon's approach proposed by Hanley and McNeal (1982). If
984    :obj:`class_index` is not specified, the first class is used as
985    "the positive" and others are negative. The result is a list of
986    tuples (aROC, standard error).
987
988    If test results consist of multiple folds, you need to split them using
989    :obj:`split_by_iterations` and perform this test on each fold separately.
990    """
991    useweights = res.weights and not ignore_weights
992    problists, tots = corn.computeROCCumulative(res, class_index, useweights)
993
994    results = []
995
996    totPos, totNeg = tots[1], tots[0]
997    N = totPos + totNeg
998    for plist in problists:
999        highPos, lowNeg = totPos, 0.
1000        W, Q1, Q2 = 0., 0., 0.
1001        for prob in plist:
1002            thisPos, thisNeg = prob[1][1], prob[1][0]
1003            highPos -= thisPos
1004            W += thisNeg * (highPos + thisPos / 2.)
1005            Q2 += thisPos * (lowNeg**2  + lowNeg*thisNeg  + thisNeg**2 / 3.)
1006            Q1 += thisNeg * (highPos**2 + highPos*thisPos + thisPos**2 / 3.)
1007
1008            lowNeg += thisNeg
1009
1010        W  /= (totPos*totNeg)
1011        Q1 /= (totNeg*totPos**2)
1012        Q2 /= (totPos*totNeg**2)
1013
1014        SE = math.sqrt((W * (1 - W) + (totPos - 1) * (Q1 - W**2) +
1015                       (totNeg - 1) * (Q2 - W**2)) / (totPos * totNeg))
1016        results.append((W, SE))
1017    return results
1018
1019AROC = AUCWilcoxon # for backward compatibility, AROC is obsolete
1020
1021
1022@deprecated_keywords({"classIndex": "class_index",
1023                      "unweighted": "ignore_weights"})
1024def compare_2_AUCs(res, lrn1, lrn2, class_index=-1,
1025                   ignore_weights=False, **argkw):
1026    return corn.compare2ROCs(res, lrn1, lrn2, class_index,
1027                             res.weights and not ignore_weights)
1028
1029# for backward compatibility, compare_2_AROCs is obsolete
1030compare_2_AROCs = compare_2_AUCs
1031
1032
1033@deprecated_keywords({"classIndex": "class_index"})
1034def compute_ROC(res, class_index=-1):
1035    """Compute a ROC curve as a list of (x, y) tuples, where x is
1036    1-specificity and y is sensitivity.
1037    """
1038    problists, tots = corn.computeROCCumulative(res, class_index)
1039
1040    results = []
1041    totPos, totNeg = tots[1], tots[0]
1042
1043    for plist in problists:
1044        curve=[(1., 1.)]
1045        TP, TN = totPos, 0.
1046        FN, FP = 0., totNeg
1047        for prob in plist:
1048            thisPos, thisNeg = prob[1][1], prob[1][0]
1049            # thisPos go from TP to FN
1050            TP -= thisPos
1051            FN += thisPos
1052            # thisNeg go from FP to TN
1053            TN += thisNeg
1054            FP -= thisNeg
1055
1056            sens = TP / (TP + FN)
1057            spec = TN / (FP + TN)
1058            curve.append((1 - spec, sens))
1059        results.append(curve)
1060
1061    return results   
1062
1063## TC's implementation of algorithms, taken from:
1064## T Fawcett: ROC Graphs: Notes and Practical Considerations for
1065## Data Mining Researchers, submitted to KDD Journal.
1066def ROC_slope((P1x, P1y, P1fscore), (P2x, P2y, P2fscore)):
1067    if P1x == P2x:
1068        return 1e300
1069    return (P1y - P2y) / (P1x - P2x)
1070
1071
1072@deprecated_keywords({"keepConcavities": "keep_concavities"})
1073def ROC_add_point(P, R, keep_concavities=1):
1074    if keep_concavities:
1075        R.append(P)
1076    else:
1077        while True:
1078            if len(R) < 2:
1079                R.append(P)
1080                return R
1081            else:
1082                T = R.pop()
1083                T2 = R[-1]
1084                if ROC_slope(T2, T) > ROC_slope(T, P):
1085                    R.append(T)
1086                    R.append(P)
1087                    return R
1088    return R
1089
1090
1091@deprecated_keywords({"classIndex": "class_index",
1092                      "keepConcavities": "keep_concavities"})
1093def TC_compute_ROC(res, class_index=-1, keep_concavities=1):
1094    problists, tots = corn.computeROCCumulative(res, class_index)
1095
1096    results = []
1097    P, N = tots[1], tots[0]
1098
1099    for plist in problists:
1100        ## corn gives an increasing by scores list, we need a decreasing
1101        plist.reverse()
1102        TP = 0.
1103        FP = 0.
1104        curve=[]
1105        fPrev = 10e300 # "infinity" score at 0., 0.
1106        for prob in plist:
1107            f = prob[0]
1108            if f != fPrev:
1109                if P:
1110                    tpr = TP / P
1111                else:
1112                    tpr = 0.
1113                if N:
1114                    fpr = FP / N
1115                else:
1116                    fpr = 0.
1117                curve = ROC_add_point((fpr, tpr, fPrev), curve,
1118                                      keep_concavities)
1119                fPrev = f
1120            thisPos, thisNeg = prob[1][1], prob[1][0]
1121            TP += thisPos
1122            FP += thisNeg
1123        if P:
1124            tpr = TP / P
1125        else:
1126            tpr = 0.
1127        if N:
1128            fpr = FP / N
1129        else:
1130            fpr = 0.
1131        curve = ROC_add_point((fpr, tpr, f), curve, keep_concavities) ## ugly
1132        results.append(curve)
1133
1134    return results
1135
1136## returns a list of points at the intersection of the tangential
1137## iso-performance line and the given ROC curve
1138## for given values of FPcost, FNcost and pval
1139def TC_best_thresholds_on_ROC_curve(FPcost, FNcost, pval, curve):
1140    m = (FPcost * (1. - pval)) / (FNcost * pval)
1141
1142    ## put the iso-performance line in point (0., 1.)
1143    x0, y0 = (0., 1.)
1144    x1, y1 = (1., 1. + m)
1145    d01 = math.sqrt((x1 - x0) * (x1 - x0) + (y1 - y0) * (y1 - y0))
1146
1147    ## calculate and find the closest point to the line
1148    firstp = 1
1149    mind = 0.
1150    a = x0 * y1 - x1 * y0
1151    closestPoints = []
1152    for (x, y, fscore) in curve:
1153        d = ((y0 - y1) * x + (x1 - x0) * y + a) / d01
1154        d = abs(d)
1155        if firstp or d < mind:
1156            mind, firstp = d, 0
1157            closestPoints = [(x, y, fscore)]
1158        else:
1159            if abs(d - mind) <= 0.0001: ## close enough
1160                closestPoints.append( (x, y, fscore) )
1161    return closestPoints         
1162
1163def frange(start, end=None, inc=None):
1164    """A range function, that does accept float increments..."""
1165
1166    if end is None:
1167        end = start + 0.
1168        start = 0.
1169
1170    if inc is None or inc == 0:
1171        inc = 1.
1172
1173    L = [start]
1174    while 1:
1175        next = start + len(L) * inc
1176        if inc > 0 and next >= end:
1177            L.append(end)
1178            break
1179        elif inc < 0 and next <= end:
1180            L.append(end)
1181            break
1182        L.append(next)
1183       
1184    return L
1185
1186## input ROCcurves are of form [ROCcurves1, ROCcurves2, ... ROCcurvesN],
1187## where ROCcurvesX is a set of ROC curves,
1188## where a (one) ROC curve is a set of (FP, TP) points
1189##
1190## for each (sub)set of input ROC curves
1191## returns the average ROC curve and an array of (vertical) standard deviations
1192@deprecated_keywords({"ROCcurves": "roc_curves"})
1193def TC_vertical_average_ROC(roc_curves, samples = 10):
1194    def INTERPOLATE((P1x, P1y, P1fscore), (P2x, P2y, P2fscore), X):
1195        if (P1x == P2x) or P1x < X > P2x or P1x > X < P2x:
1196            raise ValueError, "assumptions for interpolation are not met: P1 = %f,%f P2 = %f,%f X = %f" % (P1x, P1y, P2x, P2y, X)
1197        dx = float(P2x) - float(P1x)
1198        dy = float(P2y) - float(P1y)
1199        m = dy / dx
1200        return P1y + m * (X - P1x)
1201
1202    def TP_FOR_FP(FPsample, ROC, npts):
1203        i = 0
1204        while i < npts - 1:
1205            (fp, _, _) = ROC[i + 1]
1206            if fp <= FPsample:
1207                i += 1
1208            else:
1209                break
1210        (fp, tp, _) = ROC[i]
1211        if fp == FPsample:
1212            return tp
1213        elif fp < FPsample and i + 1 < len(ROC):
1214            return INTERPOLATE(ROC[i], ROC[i + 1], FPsample)
1215        elif fp < FPsample and i + 1 == len(ROC): # return the last
1216            return ROC[i][1]
1217        raise ValueError, "cannot compute: TP_FOR_FP in TC_vertical_average_ROC"
1218        #return 0.
1219
1220    average = []
1221    stdev = []
1222    for ROCS in roc_curves:
1223        npts = []
1224        for c in ROCS:
1225            npts.append(len(c))
1226        nrocs = len(ROCS)
1227
1228        TPavg = []
1229        TPstd = []
1230        for FPsample in frange(0., 1., 1. / samples):
1231            TPsum = []
1232            for i in range(nrocs):
1233                ##TPsum = TPsum + TP_FOR_FP(FPsample, ROCS[i], npts[i])
1234                TPsum.append(TP_FOR_FP(FPsample, ROCS[i], npts[i]))
1235            TPavg.append((FPsample, statc.mean(TPsum)))
1236            if len(TPsum) > 1:
1237                stdv = statc.std(TPsum)
1238            else:
1239                stdv = 0.
1240            TPstd.append(stdv)
1241
1242        average.append(TPavg)
1243        stdev.append(TPstd)
1244
1245    return average, stdev
1246
1247## input ROCcurves are of form [ROCcurves1, ROCcurves2, ... ROCcurvesN],
1248## where ROCcurvesX is a set of ROC curves,
1249## where a (one) ROC curve is a set of (FP, TP) points
1250##
1251## for each (sub)set of input ROC curves
1252## returns the average ROC curve, an array of vertical standard deviations and an array of horizontal standard deviations
1253@deprecated_keywords({"ROCcurves": "roc_curves"})
1254def TC_threshold_average_ROC(roc_curves, samples=10):
1255    def POINT_AT_THRESH(ROC, npts, thresh):
1256        i = 0
1257        while i < npts - 1:
1258            (px, py, pfscore) = ROC[i]
1259            if pfscore > thresh:
1260                i += 1
1261            else:
1262                break
1263        return ROC[i]
1264
1265    average = []
1266    stdevV = []
1267    stdevH = []
1268    for ROCS in roc_curves:
1269        npts = []
1270        for c in ROCS:
1271            npts.append(len(c))
1272        nrocs = len(ROCS)
1273
1274        T = []
1275        for c in ROCS:
1276            for (px, py, pfscore) in c:
1277##                try:
1278##                    T.index(pfscore)
1279##                except:
1280                T.append(pfscore)
1281        T.sort()
1282        T.reverse() ## ugly
1283
1284        TPavg = []
1285        TPstdV = []
1286        TPstdH = []
1287        for tidx in frange(0, (len(T) - 1.), float(len(T)) / samples):
1288            FPsum = []
1289            TPsum = []
1290            for i in range(nrocs):
1291                (fp, tp, _) = POINT_AT_THRESH(ROCS[i], npts[i], T[int(tidx)])
1292                FPsum.append(fp)
1293                TPsum.append(tp)
1294            TPavg.append((statc.mean(FPsum), statc.mean(TPsum)))
1295            ## vertical standard deviation
1296            if len(TPsum) > 1:
1297                stdv = statc.std(TPsum)
1298            else:
1299                stdv = 0.
1300            TPstdV.append( stdv )
1301            ## horizontal standard deviation
1302            if len(FPsum) > 1:
1303                stdh = statc.std(FPsum)
1304            else:
1305                stdh = 0.
1306            TPstdH.append( stdh )
1307
1308        average.append(TPavg)
1309        stdevV.append(TPstdV)
1310        stdevH.append(TPstdH)
1311
1312    return average, stdevV, stdevH
1313
1314## Calibration Curve
1315## returns an array of (curve, yesClassPredictions, noClassPredictions)
1316## elements, where:
1317##  - curve is an array of points (x, y) on the calibration curve
1318##  - yesClassRugPoints is an array of (x, 1) points
1319##  - noClassRugPoints is an array of (x, 0) points
1320@deprecated_keywords({"classIndex": "class_index"})
1321def compute_calibration_curve(res, class_index=-1):
1322    ## merge multiple iterations into one
1323    mres = Orange.evaluation.testing.ExperimentResults(1, res.classifier_names,
1324        res.class_values, res.weights, classifiers=res.classifiers,
1325        loaded=res.loaded, test_type=res.test_type, labels=res.labels)
1326    for te in res.results:
1327        mres.results.append(te)
1328
1329    problists, tots = corn.computeROCCumulative(mres, class_index)
1330
1331    results = []
1332
1333    bins = 10 ## divide interval between 0. and 1. into N bins
1334
1335    for plist in problists:
1336        yesClassRugPoints = [] 
1337        noClassRugPoints = []
1338
1339        yesBinsVals = [0] * bins
1340        noBinsVals = [0] * bins
1341        for (f, (thisNeg, thisPos)) in plist:
1342            yesClassRugPoints.append((f, thisPos)) # 1.
1343            noClassRugPoints.append((f, thisNeg)) # 1.
1344
1345            index = int(f * bins )
1346            index = min(index, bins - 1) ## just in case for value 1.
1347            yesBinsVals[index] += thisPos
1348            noBinsVals[index] += thisNeg
1349
1350        curve = []
1351        for cn in range(bins):
1352            f = float(cn * 1. / bins) + (1. / 2. / bins)
1353            yesVal = yesBinsVals[cn]
1354            noVal = noBinsVals[cn]
1355            allVal = yesVal + noVal
1356            if allVal == 0.: continue
1357            y = float(yesVal) / float(allVal)
1358            curve.append((f,  y))
1359
1360        ## smooth the curve
1361        maxnPoints = 100
1362        if len(curve) >= 3:
1363#            loessCurve = statc.loess(curve, -3, 0.6)
1364            loessCurve = statc.loess(curve, maxnPoints, 0.5, 3)
1365        else:
1366            loessCurve = curve
1367        clen = len(loessCurve)
1368        if clen > maxnPoints:
1369            df = clen / maxnPoints
1370            if df < 1: df = 1
1371            curve = [loessCurve[i]  for i in range(0, clen, df)]
1372        else:
1373            curve = loessCurve
1374        ## remove the third value (variance of epsilon?) that suddenly
1375        ## appeared in the output of the statc.loess function
1376        curve = [c[:2] for c in curve]
1377        results.append((curve, yesClassRugPoints, noClassRugPoints))
1378
1379    return results
1380
1381
1382## Lift Curve
1383## returns an array of curve elements, where:
1384##  - curve is an array of points ((TP + FP) / (P + N), TP / P, (th, FP / N))
1385##    on the Lift Curve
1386@deprecated_keywords({"classIndex": "class_index"})
1387def compute_lift_curve(res, class_index=-1):
1388    ## merge multiple iterations into one
1389    mres = Orange.evaluation.testing.ExperimentResults(1, res.classifier_names,
1390        res.class_values, res.weights, classifiers=res.classifiers,
1391        loaded=res.loaded, test_type=res.test_type, labels=res.labels)
1392    for te in res.results:
1393        mres.results.append(te)
1394
1395    problists, tots = corn.computeROCCumulative(mres, class_index)
1396
1397    results = []
1398    P, N = tots[1], tots[0]
1399    for plist in problists:
1400        ## corn gives an increasing by scores list, we need a decreasing
1401        plist.reverse()
1402        TP = 0.
1403        FP = 0.
1404        curve = [(0., 0., (10e300, 0.))]
1405        for (f, (thisNeg, thisPos)) in plist:
1406            TP += thisPos
1407            FP += thisNeg
1408            curve.append(((TP + FP) / (P + N), TP, (f, FP / (N or 1))))
1409        results.append(curve)
1410
1411    return P, N, results
1412
1413
1414class CDT:
1415  """Store the number of concordant (C), discordant (D) and tied (T) pairs."""
1416  def __init__(self, C=0., D=0., T=0.):
1417    self.C, self.D, self.T = C, D, T
1418   
1419def is_CDT_empty(cdt):
1420    return cdt.C + cdt.D + cdt.T < 1e-20
1421
1422
1423@deprecated_keywords({"classIndex": "class_index",
1424                      "unweighted": "ignore_weights"})
1425def compute_CDT(res, class_index=-1, ignore_weights=False, **argkw):
1426    """Obsolete, don't use."""
1427    if class_index < 0:
1428        if res.baseClass >= 0:
1429            class_index = res.baseClass
1430        else:
1431            class_index = 1
1432           
1433    useweights = res.weights and not ignore_weights
1434    weightByClasses = argkw.get("weightByClasses", True)
1435
1436    if res.number_of_iterations>1:
1437        CDTs = [CDT() for _ in range(res.number_of_learners)]
1438        iterationExperiments = split_by_iterations(res)
1439        for exp in iterationExperiments:
1440            expCDTs = corn.computeCDT(exp, class_index, useweights)
1441            for i in range(len(CDTs)):
1442                CDTs[i].C += expCDTs[i].C
1443                CDTs[i].D += expCDTs[i].D
1444                CDTs[i].T += expCDTs[i].T
1445        for i in range(res.number_of_learners):
1446            if is_CDT_empty(CDTs[0]):
1447                return corn.computeCDT(res, class_index, useweights)
1448       
1449        return CDTs
1450    else:
1451        return corn.computeCDT(res, class_index, useweights)
1452
1453# Backward compatibility
1454def replace_use_weights(fun):
1455    if environ.orange_no_deprecated_members:
1456        return fun
1457
1458    @functools.wraps(fun)
1459    def wrapped(*args, **kwargs):
1460        use_weights = kwargs.pop("useWeights", None)
1461        if use_weights is not None:
1462            deprecation_warning("useWeights", "ignore_weights")
1463            kwargs["ignore_weights"] = not use_weights
1464        return fun(*args, **kwargs)
1465    return wrapped
1466
1467class AUC(list):
1468    """
1469    Compute the area under ROC curve given a set of experimental results.
1470    For multivalued class problems, return the result of
1471    :obj:`by_weighted_pairs`.
1472    If testing consisted of multiple folds, each fold is scored and the
1473    average score is returned. If a fold contains only instances with the
1474    same class value, folds will be merged.
1475
1476    :param test_results: test results to score
1477    :param ignore_weights: ignore instance weights when calculating score
1478    :param method: DEPRECATED, call the appropriate method directly.
1479    """
1480
1481    @replace_use_weights
1482    def __init__(self, test_results=None, method=0, ignore_weights=False):
1483
1484        super(AUC, self).__init__()
1485
1486        self.ignore_weights=ignore_weights
1487        self.method=method
1488
1489        if test_results is not None:
1490            self.__call__(test_results)
1491
1492    def __call__(self, test_results):
1493        if len(test_results.class_values) < 2:
1494            raise ValueError("Cannot compute AUC on a single-class problem")
1495        elif len(test_results.class_values) == 2:
1496            self._compute_for_binary_class(test_results)
1497        else:
1498            self._compute_for_multi_value_class(test_results, self.method)
1499
1500    @classmethod
1501    def by_weighted_pairs(cls, res, ignore_weights=False):
1502        """
1503        Compute AUC for each pair of classes (ignoring instances of all other
1504        classes) and average the results, weighting them by the number of
1505        pairs of instances from these two classes (e.g. by the product of
1506        probabilities of the two classes). AUC computed in this way still
1507        behaves as the concordance index, e.g., gives the probability that two
1508        randomly chosen instances from different classes will be correctly
1509        recognized (if the classifier knows from which two classes the
1510        instances came).
1511        """
1512        auc = AUC(ignore_weights=ignore_weights)
1513        auc._compute_for_multi_value_class(res, method=cls.ByWeightedPairs)
1514        return auc
1515
1516    @classmethod
1517    def by_pairs(cls, res, ignore_weights=False):
1518        """
1519        Similar to by_weighted_pairs, except that the average over class pairs
1520        is not weighted. This AUC is, like the binary version, independent of
1521        class distributions, but it is not related to the concordance index
1522        any more.
1523        """
1524        auc = AUC(ignore_weights=ignore_weights)
1525        auc._compute_for_multi_value_class(res, method=cls.ByPairs)
1526        return auc
1527
1528    @classmethod
1529    def weighted_one_against_all(cls, res, ignore_weights=False):
1530        """
1531        For each class, it computes AUC for this class against all others (that
1532        is, treating other classes as one class). The AUCs are then averaged by
1533        the class probabilities. This is related to the concordance index in
1534        which we test the classifier's (average) capability of distinguishing
1535        the instances from a specified class from those that come from other
1536        classes.
1537        Unlike the binary AUC, the measure is not independent of class
1538        distributions.
1539        """
1540        auc = AUC(ignore_weights=ignore_weights)
1541        auc._compute_for_multi_value_class(res,
1542            method=cls.WeightedOneAgainstAll)
1543        return auc
1544
1545    @classmethod
1546    def one_against_all(cls, res, ignore_weights=False):
1547        """
1548        Similar to weighted_one_against_all, except that the average
1549        is not weighted.
1550        """
1551        auc = AUC(ignore_weights=ignore_weights)
1552        auc._compute_for_multi_value_class(res, method=cls.OneAgainstAll)
1553        return auc
1554
1555    @classmethod
1556    def single_class(cls, res, class_index=-1, ignore_weights=False):
1557        """
1558        Compute AUC where the class with the given class_index is singled
1559        out and all other classes are treated as a single class.
1560        """
1561        if class_index < 0:
1562            if res.base_class >= 0:
1563                class_index = res.base_class
1564            else:
1565                class_index = 1
1566
1567        auc = AUC(ignore_weights=ignore_weights)
1568        auc._compute_for_single_class(res, class_index)
1569        return auc
1570
1571    @classmethod
1572    def pair(cls, res, class_index1, class_index2, ignore_weights=False):
1573        """
1574        Computes AUC between a pair of classes, ignoring instances from all
1575        other classes.
1576        """
1577        auc = AUC(ignore_weights=ignore_weights)
1578        auc._compute_for_pair_of_classes(res, class_index1, class_index2)
1579        return auc
1580
1581    @classmethod
1582    def matrix(cls, res, ignore_weights=False):
1583        """
1584        Compute a (lower diagonal) matrix with AUCs for all pairs of classes.
1585        If there are empty classes, the corresponding elements in the matrix
1586        are -1.
1587        """
1588        auc = AUC(ignore_weights=ignore_weights)
1589        auc._compute_matrix(res)
1590        return auc
1591
1592    def _compute_for_binary_class(self, res):
1593        """AUC for binary classification problems"""
1594        if res.number_of_iterations > 1:
1595            return self._compute_for_multiple_folds(
1596                self._compute_one_class_against_all,
1597                split_by_iterations(res),
1598                (-1, res, res.number_of_iterations))
1599        else:
1600            auc, _ = self._compute_one_class_against_all(res, -1)
1601            self[:] = auc
1602            return self
1603
1604    def _compute_for_multi_value_class(self, res, method=0):
1605        """AUC for multiclass classification problems"""
1606        numberOfClasses = len(res.class_values)
1607
1608        if res.number_of_iterations > 1:
1609            iterations = split_by_iterations(res)
1610            all_ite = res
1611        else:
1612            iterations = [res]
1613            all_ite = None
1614
1615        # by pairs
1616        sum_aucs = [0.] * res.number_of_learners
1617        usefulClassPairs = 0.
1618
1619        prob = None
1620        if method in [self.ByWeightedPairs, self.WeightedOneAgainstAll]:
1621            prob = class_probabilities_from_res(res)
1622
1623        if method in [self.ByWeightedPairs, self.ByPairs]:
1624            for classIndex1 in range(numberOfClasses):
1625                for classIndex2 in range(classIndex1):
1626                    subsum_aucs = self._compute_for_multiple_folds(
1627                        self._compute_one_class_against_another, iterations,
1628                        (classIndex1, classIndex2, all_ite,
1629                        res.number_of_iterations))
1630                    if subsum_aucs:
1631                        if method == self.ByWeightedPairs:
1632                            p_ij = prob[classIndex1] * prob[classIndex2]
1633                            subsum_aucs = [x * p_ij  for x in subsum_aucs]
1634                            usefulClassPairs += p_ij
1635                        else:
1636                            usefulClassPairs += 1
1637                        sum_aucs = map(add, sum_aucs, subsum_aucs)
1638        else:
1639            for classIndex in range(numberOfClasses):
1640                subsum_aucs = self._compute_for_multiple_folds(
1641                    self._compute_one_class_against_all,
1642                    iterations, (classIndex, all_ite,
1643                                 res.number_of_iterations))
1644                if subsum_aucs:
1645                    if method == self.ByWeightedPairs:
1646                        p_i = prob[classIndex]
1647                        subsum_aucs = [x * p_i  for x in subsum_aucs]
1648                        usefulClassPairs += p_i
1649                    else:
1650                        usefulClassPairs += 1
1651                    sum_aucs = map(add, sum_aucs, subsum_aucs)
1652
1653        if usefulClassPairs > 0:
1654            sum_aucs = [x/usefulClassPairs for x in sum_aucs]
1655
1656        self[:] = sum_aucs
1657        return self
1658
1659    # computes the average AUC over folds using "AUCcomputer" (AUC_i or AUC_ij)
1660    # it returns the sum of what is returned by the computer,
1661    # unless at a certain fold the computer has to resort to computing
1662    # over all folds or even this failed;
1663    # in these cases the result is returned immediately
1664    def _compute_for_multiple_folds(self, auc_computer, iterations,
1665                                 computer_args):
1666        """Compute the average AUC over folds using :obj:`auc_computer`."""
1667        subsum_aucs = [0.] * iterations[0].number_of_learners
1668        for ite in iterations:
1669            aucs, foldsUsed = auc_computer(*(ite, ) + computer_args)
1670            if not aucs:
1671                import warnings
1672                warnings.warn("AUC cannot be computed (all instances belong to the same class).")
1673                return
1674            if not foldsUsed:
1675                self[:] = aucs
1676                return aucs
1677            subsum_aucs = map(add, subsum_aucs, aucs)
1678        self[:] = subsum_aucs
1679        return self
1680
1681    # Computes AUC
1682    # in multivalued class problem, AUC is computed as one against all
1683    # results over folds are averages
1684    # if some folds examples from one class only, the folds are merged
1685    def _compute_for_single_class(self, res, class_index):
1686        if res.number_of_iterations > 1:
1687            self._compute_for_multiple_folds(
1688                self._compute_one_class_against_all, split_by_iterations(res),
1689                (class_index, res, res.number_of_iterations))
1690        else:
1691            self._compute_one_class_against_all(res, class_index)
1692
1693    # Computes AUC for a pair of classes (as if there were no other classes)
1694    # results over folds are averages
1695    # if some folds have examples from one class only, the folds are merged
1696    def _compute_for_pair_of_classes(self, res, class_index1, class_index2):
1697        if res.number_of_iterations > 1:
1698            self._compute_for_multiple_folds(
1699                self._compute_one_class_against_another,
1700                split_by_iterations(res),
1701                (class_index1, class_index2, res, res.number_of_iterations))
1702        else:
1703            self._compute_one_class_against_another(res, class_index1,
1704                                                    class_index2)
1705
1706    # computes AUC between class i and the other classes
1707    # (treating them as the same class)
1708    @deprecated_keywords({"classIndex": "class_index",
1709                          "divideByIfIte": "divide_by_if_ite"})
1710    def _compute_one_class_against_all(self, ite, class_index, all_ite=None,
1711                                      divide_by_if_ite=1.0):
1712        """Compute AUC between class i and all the other classes)"""
1713        return self._compute_auc(corn.computeCDT, ite, all_ite,
1714            divide_by_if_ite, (class_index, not self.ignore_weights))
1715
1716
1717    # computes AUC between classes i and j as if there are no other classes
1718    def _compute_one_class_against_another(
1719        self, ite, class_index1, class_index2,
1720        all_ite=None, divide_by_if_ite=1.):
1721        """
1722        Compute AUC between classes i and j as if there are no other classes.
1723        """
1724        return self._compute_auc(corn.computeCDTPair, ite,
1725            all_ite, divide_by_if_ite,
1726            (class_index1, class_index2, not self.ignore_weights))
1727
1728    # computes AUC using a specified 'cdtComputer' function
1729    # It tries to compute AUCs from 'ite' (examples from a single iteration)
1730    # and, if C+D+T=0, from 'all_ite' (entire test set). In the former case,
1731    # the AUCs are divided by 'divideByIfIte'.
1732    # Additional flag is returned which is True in the former case,
1733    # or False in the latter.
1734    @deprecated_keywords({"cdt_computer": "cdtComputer",
1735                          "divideByIfIte": "divide_by_if_ite",
1736                          "computerArgs": "computer_args"})
1737    def _compute_auc(self, cdt_computer, ite, all_ite, divide_by_if_ite,
1738                     computer_args):
1739        """
1740        Compute AUC using a :obj:`cdt_computer`.
1741        """
1742        cdts = cdt_computer(*(ite, ) + computer_args)
1743        if not is_CDT_empty(cdts[0]):
1744            return [(cdt.C + cdt.T / 2) / (cdt.C + cdt.D + cdt.T) /
1745                    divide_by_if_ite for cdt in cdts], True
1746
1747        if all_ite:
1748            cdts = cdt_computer(*(all_ite, ) + computer_args)
1749            if not is_CDT_empty(cdts[0]):
1750                return [(cdt.C + cdt.T / 2) / (cdt.C + cdt.D + cdt.T)
1751                        for cdt in cdts], False
1752
1753        return False, False
1754
1755    def _compute_matrix(self, res):
1756        numberOfClasses = len(res.class_values)
1757        number_of_learners = res.number_of_learners
1758        if res.number_of_iterations > 1:
1759            iterations, all_ite = split_by_iterations(res), res
1760        else:
1761            iterations, all_ite = [res], None
1762        aucs = [[[] for _ in range(numberOfClasses)]
1763                for _ in range(number_of_learners)]
1764        for classIndex1 in range(numberOfClasses):
1765            for classIndex2 in range(classIndex1):
1766                pair_aucs = self._compute_for_multiple_folds(
1767                    self._compute_one_class_against_another, iterations,
1768                    (classIndex1, classIndex2, all_ite,
1769                     res.number_of_iterations))
1770                if pair_aucs:
1771                    for lrn in range(number_of_learners):
1772                        aucs[lrn][classIndex1].append(pair_aucs[lrn])
1773                else:
1774                    for lrn in range(number_of_learners):
1775                        aucs[lrn][classIndex1].append(-1)
1776        self[:] = aucs
1777        return aucs
1778
1779#Backward compatibility
1780AUC.ByWeightedPairs = 0
1781AUC.ByPairs = 1
1782AUC.WeightedOneAgainstAll = 2
1783AUC.OneAgainstAll = 3
1784
1785@replace_use_weights
1786def AUC_binary(res, ignore_weights=False):
1787    auc = deprecated_function_name(AUC)(ignore_weights=ignore_weights)
1788    auc._compute_for_binary_class(res)
1789    return auc
1790
1791@replace_use_weights
1792def AUC_multi(res, ignore_weights=False, method=0):
1793    auc = deprecated_function_name(AUC)(ignore_weights=ignore_weights,
1794        method=method)
1795    auc._compute_for_multi_value_class(res)
1796    return auc
1797
1798
1799@deprecated_keywords({"AUCcomputer": "auc_computer",
1800                      "computerArgs": "computer_args"})
1801def AUC_iterations(auc_computer, iterations, computer_args):
1802    auc = deprecated_function_name(AUC)()
1803    auc._compute_for_multiple_folds(auc_computer, iterations, computer_args)
1804    return auc
1805
1806def AUC_x(cdtComputer, ite, all_ite, divide_by_if_ite, computer_args):
1807    auc = deprecated_function_name(AUC)()
1808    result = auc._compute_auc(cdtComputer, ite, all_ite, divide_by_if_ite,
1809                              computer_args)
1810    return result
1811
1812@replace_use_weights
1813def AUC_i(ite, class_index, ignore_weights=False, all_ite=None,
1814          divide_by_if_ite=1.):
1815    auc = deprecated_function_name(AUC)()
1816    result = auc._compute_one_class_against_another(ite, class_index,
1817        all_ite=None, divide_by_if_ite=1.)
1818    return result
1819
1820
1821@replace_use_weights
1822def AUC_ij(ite, class_index1, class_index2, ignore_weights=False,
1823           all_ite=None, divide_by_if_ite=1.):
1824    auc = deprecated_function_name(AUC)(ignore_weights=ignore_weights)
1825    result = auc._compute_one_class_against_another(
1826        ite, class_index1, class_index2, all_ite=None, divide_by_if_ite=1.)
1827    return result
1828
1829
1830
1831
1832#AUC_binary = replace_use_weights(deprecated_function_name(AUC()._compute_for_binary_class))
1833#AUC_multi = replace_use_weights(deprecated_function_name(AUC._compute_for_multi_value_class))
1834#AUC_iterations = replace_use_weights(deprecated_function_name(AUC._compute_for_multiple_folds))
1835#AUC_x = replace_use_weights(deprecated_function_name(AUC._compute_auc))
1836#AUC_i = replace_use_weights(deprecated_function_name(AUC._compute_one_class_against_all))
1837#AUC_ij = replace_use_weights(deprecated_function_name(AUC._compute_one_class_against_another))
1838
1839AUC_single = replace_use_weights(
1840             deprecated_keywords({"classIndex": "class_index"})(
1841             deprecated_function_name(AUC.single_class)))
1842AUC_pair = replace_use_weights(
1843           deprecated_keywords({"classIndex1": "class_index1",
1844                                "classIndex2": "class_index2"})(
1845           deprecated_function_name(AUC.pair)))
1846AUC_matrix = replace_use_weights(deprecated_function_name(AUC.matrix))
1847
1848
1849@deprecated_keywords({"unweighted": "ignore_weights"})
1850def McNemar(res, ignore_weights=False, **argkw):
1851    """
1852    Compute a triangular matrix with McNemar statistics for each pair of
1853    classifiers. The statistics is distributed by chi-square distribution with
1854    one degree of freedom; critical value for 5% significance is around 3.84.
1855    """
1856    nLearners = res.number_of_learners
1857    mcm = []
1858    for i in range(nLearners):
1859       mcm.append([0.] * res.number_of_learners)
1860
1861    if not res.weights or ignore_weights:
1862        for i in res.results:
1863            actual = i.actual_class
1864            classes = i.classes
1865            for l1 in range(nLearners):
1866                for l2 in range(l1, nLearners):
1867                    if classes[l1] == actual:
1868                        if classes[l2] != actual:
1869                            mcm[l1][l2] += 1
1870                    elif classes[l2] == actual:
1871                        mcm[l2][l1] += 1
1872    else:
1873        for i in res.results:
1874            actual = i.actual_class
1875            classes = i.classes
1876            for l1 in range(nLearners):
1877                for l2 in range(l1, nLearners):
1878                    if classes[l1] == actual:
1879                        if classes[l2] != actual:
1880                            mcm[l1][l2] += i.weight
1881                    elif classes[l2] == actual:
1882                        mcm[l2][l1] += i.weight
1883
1884    for l1 in range(nLearners):
1885        for l2 in range(l1, nLearners):
1886            su = mcm[l1][l2] + mcm[l2][l1]
1887            if su:
1888                mcm[l2][l1] = (abs(mcm[l1][l2] - mcm[l2][l1]) - 1)**2 / su
1889            else:
1890                mcm[l2][l1] = 0
1891
1892    for l1 in range(nLearners):
1893        mcm[l1]=mcm[l1][:l1]
1894
1895    return mcm
1896
1897
1898def McNemar_of_two(res, lrn1, lrn2, ignore_weights=False):
1899    """
1900    McNemar_of_two computes a McNemar statistics for a pair of classifier,
1901    specified by indices learner1 and learner2.
1902    """
1903    tf = ft = 0.
1904    if not res.weights or ignore_weights:
1905        for i in res.results:
1906            actual = i.actual_class
1907            if i.classes[lrn1] == actual:
1908                if i.classes[lrn2] != actual:
1909                    tf += i.weight
1910            elif i.classes[lrn2] == actual:
1911                    ft += i.weight
1912    else:
1913        for i in res.results:
1914            actual = i.actual_class
1915            if i.classes[lrn1] == actual:
1916                if i.classes[lrn2] != actual:
1917                    tf += 1.
1918            elif i.classes[lrn2] == actual:
1919                    ft += 1.
1920
1921    su = tf + ft
1922    if su:
1923        return (abs(tf - ft) - 1)**2 / su
1924    else:
1925        return 0
1926
1927
1928def Friedman(res, stat=CA):
1929    """
1930    Compare classifiers with Friedman test, treating folds as different examles.
1931    Returns F, p and average ranks.
1932    """
1933    res_split = split_by_iterations(res)
1934    res = [stat(r) for r in res_split]
1935   
1936    N = len(res)
1937    k = len(res[0])
1938    sums = [0.] * k
1939    for r in res:
1940        ranks = [k - x + 1 for x in statc.rankdata(r)]
1941        if stat == Brier_score: # reverse ranks for Brier_score (lower better)
1942            ranks = [k + 1 - x for x in ranks]
1943        sums = map(add, ranks, sums)
1944
1945    T = sum(x * x for x in sums)
1946    sums = [x / N for x in sums]
1947
1948    F = 12. / (N * k * (k + 1)) * T  - 3 * N * (k + 1)
1949
1950    return F, statc.chisqprob(F, k - 1), sums
1951
1952
1953def Wilcoxon_pairs(res, avgranks, stat=CA):
1954    """
1955    Return a triangular matrix, where element[i][j] stores significance of
1956    difference between the i-th and the j-th classifier, as computed by the
1957    Wilcoxon test. The element is positive if the i-th is better than the j-th,
1958    negative if it is worse, and 1 if they are equal.
1959    Arguments are ExperimentResults, average ranks (as returned by Friedman)
1960    and, optionally, a statistics; greater values should mean better results.
1961    """
1962    res_split = split_by_iterations(res)
1963    res = [stat(r) for r in res_split]
1964
1965    k = len(res[0])
1966    bt = []
1967    for m1 in range(k):
1968        nl = []
1969        for m2 in range(m1 + 1, k):
1970            t, p = statc.wilcoxont([r[m1] for r in res], [r[m2] for r in res])
1971            if avgranks[m1]<avgranks[m2]:
1972                nl.append(p)
1973            elif avgranks[m2]<avgranks[m1]:
1974                nl.append(-p)
1975            else:
1976                nl.append(1)
1977        bt.append(nl)
1978    return bt
1979
1980
1981@deprecated_keywords({"allResults": "all_results",
1982                      "noConfidence": "no_confidence"})
1983def plot_learning_curve_learners(file, all_results, proportions, learners,
1984                                 no_confidence=0):
1985    plot_learning_curve(file, all_results, proportions,
1986        [Orange.misc.getobjectname(learners[i], "Learner %i" % i)
1987        for i in range(len(learners))], no_confidence)
1988
1989
1990@deprecated_keywords({"allResults": "all_results",
1991                      "noConfidence": "no_confidence"})
1992def plot_learning_curve(file, all_results, proportions, legend,
1993                        no_confidence=0):
1994    import types
1995    fopened = 0
1996    if type(file) == types.StringType:
1997        file = open(file, "wt")
1998        fopened = 1
1999       
2000    file.write("set yrange [0:1]\n")
2001    file.write("set xrange [%f:%f]\n" % (proportions[0], proportions[-1]))
2002    file.write("set multiplot\n\n")
2003    CAs = [CA(x, report_se=True) for x in all_results]
2004
2005    file.write("plot \\\n")
2006    for i in range(len(legend) - 1):
2007        if not no_confidence:
2008            file.write("'-' title '' with yerrorbars pointtype %i,\\\n" % (i+1))
2009        file.write("'-' title '%s' with linespoints pointtype %i,\\\n" % (legend[i], i+1))
2010    if not no_confidence:
2011        file.write("'-' title '' with yerrorbars pointtype %i,\\\n" % (len(legend)))
2012    file.write("'-' title '%s' with linespoints pointtype %i\n" % (legend[-1], len(legend)))
2013
2014    for i in range(len(legend)):
2015        if not no_confidence:
2016            for p in range(len(proportions)):
2017                file.write("%f\t%f\t%f\n" % (proportions[p], CAs[p][i][0], 1.96*CAs[p][i][1]))
2018            file.write("e\n\n")
2019
2020        for p in range(len(proportions)):
2021            file.write("%f\t%f\n" % (proportions[p], CAs[p][i][0]))
2022        file.write("e\n\n")
2023
2024    if fopened:
2025        file.close()
2026
2027
2028def print_single_ROC_curve_coordinates(file, curve):
2029    import types
2030    fopened=0
2031    if type(file)==types.StringType:
2032        file=open(file, "wt")
2033        fopened=1
2034
2035    for coord in curve:
2036        file.write("%5.3f\t%5.3f\n" % tuple(coord))
2037
2038    if fopened:
2039        file.close()
2040
2041
2042def plot_ROC_learners(file, curves, learners):
2043    plot_ROC(file, curves, [Orange.misc.getobjectname(learners[i], "Learner %i" % i) for i in range(len(learners))])
2044   
2045def plot_ROC(file, curves, legend):
2046    import types
2047    fopened=0
2048    if type(file)==types.StringType:
2049        file=open(file, "wt")
2050        fopened=1
2051
2052    file.write("set yrange [0:1]\n")
2053    file.write("set xrange [0:1]\n")
2054    file.write("set multiplot\n\n")
2055
2056    file.write("plot \\\n")
2057    for leg in legend:
2058        file.write("'-' title '%s' with lines,\\\n" % leg)
2059    file.write("'-' title '' with lines\n")
2060
2061    for curve in curves:
2062        for coord in curve:
2063            file.write("%5.3f\t%5.3f\n" % tuple(coord))
2064        file.write("e\n\n")
2065
2066    file.write("1.0\t1.0\n0.0\t0.0e\n\n")         
2067
2068    if fopened:
2069        file.close()
2070
2071
2072@deprecated_keywords({"allResults": "all_results"})
2073def plot_McNemar_curve_learners(file, all_results, proportions, learners, reference=-1):
2074    plot_McNemar_curve(file, all_results, proportions, [Orange.misc.getobjectname(learners[i], "Learner %i" % i) for i in range(len(learners))], reference)
2075
2076
2077@deprecated_keywords({"allResults": "all_results"})
2078def plot_McNemar_curve(file, all_results, proportions, legend, reference=-1):
2079    if reference<0:
2080        reference=len(legend)-1
2081       
2082    import types
2083    fopened=0
2084    if type(file)==types.StringType:
2085        file=open(file, "wt")
2086        fopened=1
2087       
2088    #file.write("set yrange [0:1]\n")
2089    #file.write("set xrange [%f:%f]\n" % (proportions[0], proportions[-1]))
2090    file.write("set multiplot\n\n")
2091    file.write("plot \\\n")
2092    tmap=range(reference)+range(reference+1, len(legend))
2093    for i in tmap[:-1]:
2094        file.write("'-' title '%s' with linespoints pointtype %i,\\\n" % (legend[i], i+1))
2095    file.write("'-' title '%s' with linespoints pointtype %i\n" % (legend[tmap[-1]], tmap[-1]))
2096    file.write("\n")
2097
2098    for i in tmap:
2099        for p in range(len(proportions)):
2100            file.write("%f\t%f\n" % (proportions[p], McNemar_of_two(all_results[p], i, reference)))
2101        file.write("e\n\n")
2102
2103    if fopened:
2104        file.close()
2105
2106default_point_types=("{$\\circ$}", "{$\\diamond$}", "{$+$}", "{$\\times$}", "{$|$}")+tuple([chr(x) for x in range(97, 122)])
2107default_line_types=("\\setsolid", "\\setdashpattern <4pt, 2pt>", "\\setdashpattern <8pt, 2pt>", "\\setdashes", "\\setdots")
2108
2109@deprecated_keywords({"allResults": "all_results"})
2110def learning_curve_learners_to_PiCTeX(file, all_results, proportions, **options):
2111    return apply(learning_curve_to_PiCTeX, (file, all_results, proportions), options)
2112
2113
2114@deprecated_keywords({"allResults": "all_results"})
2115def learning_curve_to_PiCTeX(file, all_results, proportions, **options):
2116    import types
2117    fopened=0
2118    if type(file)==types.StringType:
2119        file=open(file, "wt")
2120        fopened=1
2121
2122    nexamples=len(all_results[0].results)
2123    CAs = [CA(x, report_se=True) for x in all_results]
2124
2125    graphsize=float(options.get("graphsize", 10.0)) #cm
2126    difprop=proportions[-1]-proportions[0]
2127    ntestexamples=nexamples*proportions[-1]
2128    xunit=graphsize/ntestexamples
2129
2130    yshift=float(options.get("yshift", -ntestexamples/20.))
2131   
2132    pointtypes=options.get("pointtypes", default_point_types)
2133    linetypes=options.get("linetypes", default_line_types)
2134
2135    if options.has_key("numberedx"):
2136        numberedx=options["numberedx"]
2137        if type(numberedx)==types.IntType:
2138            if numberedx>0:
2139                numberedx=[nexamples*proportions[int(i/float(numberedx)*len(proportions))] for i in range(numberedx)]+[proportions[-1]*nexamples]
2140            elif numberedx<0:
2141                numberedx = -numberedx
2142                newn=[]
2143                for i in range(numberedx+1):
2144                    wanted=proportions[0]+float(i)/numberedx*difprop
2145                    best=(10, 0)
2146                    for t in proportions:
2147                        td=abs(wanted-t)
2148                        if td<best[0]:
2149                            best=(td, t)
2150                    if not best[1] in newn:
2151                        newn.append(best[1])
2152                newn.sort()
2153                numberedx=[nexamples*x for x in newn]
2154        elif type(numberedx[0])==types.FloatType:
2155            numberedx=[nexamples*x for x in numberedx]
2156    else:
2157        numberedx=[nexamples*x for x in proportions]
2158
2159    file.write("\\mbox{\n")
2160    file.write(\\beginpicture\n")
2161    file.write(\\setcoordinatesystem units <%10.8fcm, %5.3fcm>\n\n" % (xunit, graphsize))   
2162    file.write(\\setplotarea x from %5.3f to %5.3f, y from 0 to 1\n" % (0, ntestexamples))   
2163    file.write(\\axis bottom invisible\n")# label {#examples}\n")
2164    file.write("      ticks short at %s /\n" % reduce(lambda x,y:x+" "+y, ["%i"%(x*nexamples+0.5) for x in proportions]))
2165    if numberedx:
2166        file.write("            long numbered at %s /\n" % reduce(lambda x,y:x+y, ["%i " % int(x+0.5) for x in numberedx]))
2167    file.write("  /\n")
2168    file.write(\\axis left invisible\n")# label {classification accuracy}\n")
2169    file.write("      shiftedto y=%5.3f\n" % yshift)
2170    file.write("      ticks short from 0.0 to 1.0 by 0.05\n")
2171    file.write("            long numbered from 0.0 to 1.0 by 0.25\n")
2172    file.write("  /\n")
2173    if options.has_key("default"):
2174        file.write(\\setdashpattern<1pt, 1pt>\n")
2175        file.write(\\plot %5.3f %5.3f %5.3f %5.3f /\n" % (0., options["default"], ntestexamples, options["default"]))
2176   
2177    for i in range(len(CAs[0])):
2178        coordinates=reduce(lambda x,y:x+" "+y, ["%i %5.3f" % (proportions[p]*nexamples, CAs[p][i][0]) for p in range(len(proportions))])
2179        if linetypes:
2180            file.write(%s\n" % linetypes[i])
2181            file.write(\\plot %s /\n" % coordinates)
2182        if pointtypes:
2183            file.write(\\multiput %s at %s /\n" % (pointtypes[i], coordinates))
2184
2185    file.write(\\endpicture\n")
2186    file.write("}\n")
2187    if fopened:
2188        file.close()
2189    file.close()
2190    del file
2191
2192def legend_learners_to_PiCTeX(file, learners, **options):
2193  return apply(legend_to_PiCTeX, (file, [Orange.misc.getobjectname(learners[i], "Learner %i" % i) for i in range(len(learners))]), options)
2194   
2195def legend_to_PiCTeX(file, legend, **options):
2196    import types
2197    fopened=0
2198    if type(file)==types.StringType:
2199        file=open(file, "wt")
2200        fopened=1
2201
2202    pointtypes=options.get("pointtypes", default_point_types)
2203    linetypes=options.get("linetypes", default_line_types)
2204
2205    file.write("\\mbox{\n")
2206    file.write(\\beginpicture\n")
2207    file.write(\\setcoordinatesystem units <5cm, 1pt>\n\n")
2208    file.write(\\setplotarea x from 0.000 to %5.3f, y from 0 to 12\n" % len(legend))
2209
2210    for i in range(len(legend)):
2211        if linetypes:
2212            file.write(%s\n" % linetypes[i])
2213            file.write(\\plot %5.3f 6 %5.3f 6 /\n" % (i, i+0.2))
2214        if pointtypes:
2215            file.write(\\put {%s} at %5.3f 6\n" % (pointtypes[i], i+0.1))
2216        file.write(\\put {%s} [lb] at %5.3f 0\n" % (legend[i], i+0.25))
2217
2218    file.write(\\endpicture\n")
2219    file.write("}\n")
2220    if fopened:
2221        file.close()
2222    file.close()
2223    del file
2224
2225
2226def compute_friedman(avranks, N):
2227    """ Returns a tuple composed of (friedman statistic, degrees of freedom)
2228    and (Iman statistic - F-distribution, degrees of freedoma) given average
2229    ranks and a number of tested data sets N.
2230    """
2231
2232    k = len(avranks)
2233
2234    def friedman(N, k, ranks):
2235        return 12*N*(sum([rank**2.0 for rank in ranks]) - (k*(k+1)*(k+1)/4.0) )/(k*(k+1))
2236
2237    def iman(fried, N, k):
2238        return (N-1)*fried/(N*(k-1) - fried)
2239
2240    f = friedman(N, k, avranks)
2241    im = iman(f, N, k)
2242    fdistdof = (k-1, (k-1)*(N-1))
2243
2244    return (f, k-1), (im, fdistdof)
2245
2246def compute_CD(avranks, N, alpha="0.05", type="nemenyi"):
2247    """ Returns critical difference for Nemenyi or Bonferroni-Dunn test
2248    according to given alpha (either alpha="0.05" or alpha="0.1") for average
2249    ranks and number of tested data sets N. Type can be either "nemenyi" for
2250    for Nemenyi two tailed test or "bonferroni-dunn" for Bonferroni-Dunn test.
2251    """
2252
2253    k = len(avranks)
2254   
2255    d = {("nemenyi", "0.05"): [0, 0, 1.959964, 2.343701, 2.569032, 2.727774,
2256                               2.849705, 2.94832, 3.030879, 3.101730, 3.163684,
2257                               3.218654, 3.268004, 3.312739, 3.353618, 3.39123,
2258                               3.426041, 3.458425, 3.488685, 3.517073, 3.543799]
2259        , ("nemenyi", "0.1"): [0, 0, 1.644854, 2.052293, 2.291341, 2.459516,
2260                               2.588521, 2.692732, 2.779884, 2.854606, 2.919889,
2261                               2.977768, 3.029694, 3.076733, 3.119693, 3.159199,
2262                               3.195743, 3.229723, 3.261461, 3.291224, 3.319233]
2263        , ("bonferroni-dunn", "0.05"): [0, 0, 1.960, 2.241, 2.394, 2.498, 2.576,
2264                                        2.638, 2.690, 2.724, 2.773],
2265         ("bonferroni-dunn", "0.1"): [0, 0, 1.645, 1.960, 2.128, 2.241, 2.326,
2266                                      2.394, 2.450, 2.498, 2.539]}
2267
2268    #can be computed in R as qtukey(0.95, n, Inf)**0.5
2269    #for (x in c(2:20)) print(qtukey(0.95, x, Inf)/(2**0.5)
2270
2271    q = d[(type, alpha)]
2272
2273    cd = q[k]*(k*(k+1)/(6.0*N))**0.5
2274
2275    return cd
2276 
2277
2278def graph_ranks(filename, avranks, names, cd=None, cdmethod=None, lowv=None, highv=None, width=6, textspace=1, reverse=False, **kwargs):
2279    """
2280    Draws a CD graph, which is used to display  the differences in methods'
2281    performance.
2282    See Janez Demsar, Statistical Comparisons of Classifiers over
2283    Multiple Data Sets, 7(Jan):1--30, 2006.
2284
2285    Needs matplotlib to work.
2286
2287    :param filename: Output file name (with extension). Formats supported
2288                     by matplotlib can be used.
2289    :param avranks: List of average methods' ranks.
2290    :param names: List of methods' names.
2291
2292    :param cd: Critical difference. Used for marking methods that whose
2293               difference is not statistically significant.
2294    :param lowv: The lowest shown rank, if None, use 1.
2295    :param highv: The highest shown rank, if None, use len(avranks).
2296    :param width: Width of the drawn figure in inches, default 6 in.
2297    :param textspace: Space on figure sides left for the description
2298                      of methods, default 1 in.
2299    :param reverse:  If True, the lowest rank is on the right. Default\: False.
2300    :param cdmethod: None by default. It can be an index of element in avranks
2301                     or or names which specifies the method which should be
2302                     marked with an interval.
2303    """
2304    if not HAS_MATPLOTLIB:
2305        import sys
2306        print >> sys.stderr, "Function requires matplotlib. Please install it."
2307        return
2308
2309    width = float(width)
2310    textspace = float(textspace)
2311
2312    def nth(l,n):
2313        """
2314        Returns only nth elemnt in a list.
2315        """
2316        n = lloc(l,n)
2317        return [ a[n] for a in l ]
2318
2319    def lloc(l,n):
2320        """
2321        List location in list of list structure.
2322        Enable the use of negative locations:
2323        -1 is the last element, -2 second last...
2324        """
2325        if n < 0:
2326            return len(l[0])+n
2327        else:
2328            return n
2329
2330    def mxrange(lr):
2331        """
2332        Multiple xranges. Can be used to traverse matrices.
2333        This function is very slow due to unknown number of
2334        parameters.
2335
2336        >>> mxrange([3,5])
2337        [(0, 0), (0, 1), (0, 2), (1, 0), (1, 1), (1, 2)]
2338
2339        >>> mxrange([[3,5,1],[9,0,-3]])
2340        [(3, 9), (3, 6), (3, 3), (4, 9), (4, 6), (4, 3)]
2341
2342        """
2343        if not len(lr):
2344            yield ()
2345        else:
2346            #it can work with single numbers
2347            index = lr[0]
2348            if type(1) == type(index):
2349                index = [ index ]
2350            for a in range(*index):
2351                for b in mxrange(lr[1:]):
2352                    yield tuple([a] + list(b))
2353
2354
2355    from matplotlib.figure import Figure
2356    from matplotlib.backends.backend_agg import FigureCanvasAgg
2357
2358    def print_figure(fig, *args, **kwargs):
2359        canvas = FigureCanvasAgg(fig)
2360        canvas.print_figure(*args, **kwargs)
2361
2362    sums = avranks
2363
2364    tempsort =  sorted([ (a,i) for i,a in  enumerate(sums) ], reverse=reverse)
2365    ssums = nth(tempsort, 0)
2366    sortidx = nth(tempsort, 1)
2367    nnames = [ names[x] for x in sortidx ]
2368   
2369    if lowv is None:
2370        lowv = min(1, int(math.floor(min(ssums))))
2371    if highv is None:
2372        highv = max(len(avranks), int(math.ceil(max(ssums))))
2373
2374    cline = 0.4
2375
2376    k = len(sums)
2377
2378    lines = None
2379
2380    linesblank = 0
2381    scalewidth = width - 2*textspace
2382
2383    def rankpos(rank):
2384        if not reverse:
2385            a = rank - lowv
2386        else:
2387            a = highv - rank
2388        return textspace+scalewidth/(highv-lowv)*a
2389
2390    distanceh = 0.25
2391
2392    if cd and cdmethod is None:
2393   
2394        #get pairs of non significant methods
2395
2396        def get_lines(sums, hsd):
2397
2398            #get all pairs
2399            lsums = len(sums)
2400            allpairs = [ (i,j) for i,j in mxrange([[lsums], [lsums]]) if j > i ]
2401            #remove not significant
2402            notSig = [ (i,j) for i,j in allpairs if abs(sums[i]-sums[j]) <= hsd ]
2403            #keep only longest
2404           
2405            def no_longer((i,j), notSig):
2406                for i1,j1 in notSig:
2407                    if (i1 <= i and j1 > j) or (i1 < i and j1 >= j):
2408                        return False
2409                return True
2410
2411            longest = [ (i,j) for i,j in notSig if no_longer((i,j),notSig) ]
2412           
2413            return longest
2414
2415        lines = get_lines(ssums, cd)
2416        linesblank = 0.2 + 0.2 + (len(lines)-1)*0.1
2417
2418        #add scale
2419        distanceh = 0.25
2420        cline += distanceh
2421
2422    #calculate height needed height of an image
2423    minnotsignificant = max(2*0.2, linesblank)
2424    height = cline + ((k+1)/2)*0.2 + minnotsignificant
2425
2426    fig = Figure(figsize=(width, height))
2427    ax = fig.add_axes([0,0,1,1]) #reverse y axis
2428    ax.set_axis_off()
2429
2430    hf = 1./height # height factor
2431    wf = 1./width
2432
2433    def hfl(l): 
2434        return [ a*hf for a in l ]
2435
2436    def wfl(l): 
2437        return [ a*wf for a in l ]
2438
2439
2440    # Upper left corner is (0,0).
2441
2442    ax.plot([0,1], [0,1], c="w")
2443    ax.set_xlim(0, 1)
2444    ax.set_ylim(1, 0)
2445
2446    def line(l, color='k', **kwargs):
2447        """
2448        Input is a list of pairs of points.
2449        """
2450        ax.plot(wfl(nth(l,0)), hfl(nth(l,1)), color=color, **kwargs)
2451
2452    def text(x, y, s, *args, **kwargs):
2453        ax.text(wf*x, hf*y, s, *args, **kwargs)
2454
2455    line([(textspace, cline), (width-textspace, cline)], linewidth=0.7)
2456   
2457    bigtick = 0.1
2458    smalltick = 0.05
2459
2460
2461    import numpy
2462    tick = None
2463    for a in list(numpy.arange(lowv, highv, 0.5)) + [highv]:
2464        tick = smalltick
2465        if a == int(a): tick = bigtick
2466        line([(rankpos(a), cline-tick/2),(rankpos(a), cline)], linewidth=0.7)
2467
2468    for a in range(lowv, highv+1):
2469        text(rankpos(a), cline-tick/2-0.05, str(a), ha="center", va="bottom")
2470
2471    k = len(ssums)
2472
2473    for i in range((k+1)/2):
2474        chei = cline+ minnotsignificant + i *0.2
2475        line([(rankpos(ssums[i]), cline), (rankpos(ssums[i]), chei), (textspace-0.1, chei)], linewidth=0.7)
2476        text(textspace-0.2, chei, nnames[i], ha="right", va="center")
2477
2478    for i in range((k+1)/2, k):
2479        chei = cline + minnotsignificant + (k-i-1)*0.2
2480        line([(rankpos(ssums[i]), cline), (rankpos(ssums[i]), chei), (textspace+scalewidth+0.1, chei)], linewidth=0.7)
2481        text(textspace+scalewidth+0.2, chei, nnames[i], ha="left", va="center")
2482
2483    if cd and cdmethod is None:
2484
2485        #upper scale
2486        if not reverse:
2487            begin, end = rankpos(lowv), rankpos(lowv+cd)
2488        else:
2489            begin, end = rankpos(highv), rankpos(highv - cd)
2490           
2491        line([(begin, distanceh), (end, distanceh)], linewidth=0.7)
2492        line([(begin, distanceh + bigtick/2), (begin, distanceh - bigtick/2)], linewidth=0.7)
2493        line([(end, distanceh + bigtick/2), (end, distanceh - bigtick/2)], linewidth=0.7)
2494        text((begin+end)/2, distanceh - 0.05, "CD", ha="center", va="bottom")
2495
2496        #non significance lines   
2497        def draw_lines(lines, side=0.05, height=0.1):
2498            start = cline + 0.2
2499            for l,r in lines: 
2500                line([(rankpos(ssums[l])-side, start), (rankpos(ssums[r])+side, start)], linewidth=2.5) 
2501                start += height
2502
2503        draw_lines(lines)
2504
2505    elif cd:
2506        begin = rankpos(avranks[cdmethod]-cd)
2507        end = rankpos(avranks[cdmethod]+cd)
2508        line([(begin, cline), (end, cline)], linewidth=2.5) 
2509        line([(begin, cline + bigtick/2), (begin, cline - bigtick/2)], linewidth=2.5)
2510        line([(end, cline + bigtick/2), (end, cline - bigtick/2)], linewidth=2.5)
2511 
2512    print_figure(fig, filename, **kwargs)
2513
2514def mlc_hamming_loss(res):
2515    """
2516    Schapire and Singer (2000) presented Hamming Loss, which id defined as:
2517   
2518    :math:`HammingLoss(H,D)=\\frac{1}{|D|} \\sum_{i=1}^{|D|} \\frac{Y_i \\vartriangle Z_i}{|L|}`
2519    """
2520    losses = [0]*res.number_of_learners
2521    label_num = len(res.labels)
2522    example_num = gettotsize(res)
2523
2524    for e in res.results:
2525        aclass = e.actual_class
2526        for i, labels in enumerate(e.classes):
2527            labels = map(int, labels)
2528            if len(labels) <> len(aclass):
2529                raise ValueError, "The dimensions of the classified output and the actual class array do not match."
2530            for j in range(label_num):
2531                if labels[j] != aclass[j]:
2532                    losses[i] += 1
2533           
2534    return [float(x)/(label_num*example_num) for x in losses]
2535
2536def mlc_accuracy(res, forgiveness_rate = 1.0):
2537    """
2538    Godbole & Sarawagi, 2004 uses the metrics accuracy, precision, recall as follows:
2539     
2540    :math:`Accuracy(H,D)=\\frac{1}{|D|} \\sum_{i=1}^{|D|} \\frac{|Y_i \\cap Z_i|}{|Y_i \\cup Z_i|}`
2541   
2542    Boutell et al. (2004) give a more generalized version using a parameter :math:`\\alpha \\ge 0`,
2543    called forgiveness rate:
2544   
2545    :math:`Accuracy(H,D)=\\frac{1}{|D|} \\sum_{i=1}^{|D|} (\\frac{|Y_i \\cap Z_i|}{|Y_i \\cup Z_i|})^{\\alpha}`
2546    """
2547    accuracies = [0.0]*res.number_of_learners
2548    example_num = gettotsize(res)
2549   
2550    for e in res.results:
2551        aclass = e.actual_class
2552        for i, labels in enumerate(e.classes):
2553            labels = map(int, labels)
2554            if len(labels) <> len(aclass):
2555                raise ValueError, "The dimensions of the classified output and the actual class array do not match."
2556           
2557            intersection = 0.0
2558            union = 0.0
2559            for real, pred in zip(labels, aclass):
2560                if real and pred:
2561                    intersection += 1
2562                if real or pred:
2563                    union += 1
2564
2565            if union:
2566                accuracies[i] += intersection / union
2567           
2568    return [math.pow(x/example_num,forgiveness_rate) for x in accuracies]
2569
2570def mlc_precision(res):
2571    """
2572    :math:`Precision(H,D)=\\frac{1}{|D|} \\sum_{i=1}^{|D|} \\frac{|Y_i \\cap Z_i|}{|Z_i|}`
2573    """
2574    precisions = [0.0]*res.number_of_learners
2575    example_num = gettotsize(res)
2576   
2577    for e in res.results:
2578        aclass = e.actual_class
2579        for i, labels in enumerate(e.classes):
2580            labels = map(int, labels)
2581            if len(labels) <> len(aclass):
2582                raise ValueError, "The dimensions of the classified output and the actual class array do not match."
2583           
2584            intersection = 0.0
2585            predicted = 0.0
2586            for real, pred in zip(labels, aclass):
2587                if real and pred:
2588                    intersection += 1
2589                if real:
2590                    predicted += 1
2591            if predicted:
2592                precisions[i] += intersection / predicted
2593           
2594    return [x/example_num for x in precisions]
2595
2596def mlc_recall(res):
2597    """
2598    :math:`Recall(H,D)=\\frac{1}{|D|} \\sum_{i=1}^{|D|} \\frac{|Y_i \\cap Z_i|}{|Y_i|}`
2599    """
2600    recalls = [0.0]*res.number_of_learners
2601    example_num = gettotsize(res)
2602   
2603    for e in res.results:
2604        aclass = e.actual_class
2605        for i, labels in enumerate(e.classes):
2606            labels = map(int, labels)
2607            if len(labels) <> len(aclass):
2608                raise ValueError, "The dimensions of the classified output and the actual class array do not match."
2609           
2610            intersection = 0.0
2611            actual = 0.0
2612            for real, pred in zip(labels, aclass):
2613                if real and pred:
2614                    intersection += 1
2615                if pred:
2616                    actual += 1
2617            if actual:
2618                recalls[i] += intersection / actual
2619           
2620    return [x/example_num for x in recalls]
2621
2622#def mlc_ranking_loss(res):
2623#    pass
2624#
2625#def mlc_average_precision(res):
2626#    pass
2627#
2628#def mlc_hierarchical_loss(res):
2629#    pass
2630
2631
2632def mt_average_score(res, score, weights=None):
2633    """
2634    Compute individual scores for each target and return the (weighted) average.
2635
2636    One method can be used to compute scores for all targets or a list of
2637    scoring methods can be passed to use different methods for different
2638    targets. In the latter case, care has to be taken if the ranges of scoring
2639    methods differ.
2640    For example, when the first target is scored from -1 to 1 (1 best) and the
2641    second from 0 to 1 (0 best), using `weights=[0.5,-1]` would scale both
2642    to a span of 1, and invert the second so that higher scores are better.
2643
2644    :param score: Single-target scoring method or a list of such methods
2645                  (one for each target).
2646    :param weights: List of real weights, one for each target,
2647                    for a weighted average.
2648
2649    """
2650    if not len(res.results):
2651        raise ValueError, "Cannot compute the score: no examples."
2652    if res.number_of_learners < 1:
2653        return []
2654    n_classes = len(res.results[0].actual_class)
2655    if weights is None:
2656        weights = [1.] * n_classes
2657    if not hasattr(score, '__len__'):
2658        score = [score] * n_classes
2659    elif len(score) != n_classes:
2660        raise ValueError, "Number of scoring methods and targets do not match."
2661    # save original classes
2662    clsss = [te.classes for te in res.results]
2663    aclsss = [te.actual_class for te in res.results]
2664    # compute single target scores
2665    single_scores = []
2666    for i in range(n_classes):
2667        for te, clss, aclss in zip(res.results, clsss, aclsss):
2668            te.classes = [cls[i] for cls in clss]
2669            te.actual_class = aclss[i]
2670        single_scores.append(score[i](res))
2671    # restore original classes
2672    for te, clss, aclss in zip(res.results, clsss, aclsss):
2673        te.classes = clss
2674        te.actual_class = aclss
2675    return [sum(w * s for w, s in zip(weights, scores)) / sum(weights)
2676        for scores in zip(*single_scores)]
2677
2678def mt_flattened_score(res, score):
2679    """
2680    Flatten (concatenate into a single list) the predictions of multiple
2681    targets and compute a single-target score.
2682   
2683    :param score: Single-target scoring method.
2684    """
2685    res2 = Orange.evaluation.testing.ExperimentResults(res.number_of_iterations,
2686        res.classifier_names, class_values=res.class_values,
2687        weights=res.weights, classifiers=res.classifiers, loaded=res.loaded,
2688        test_type=Orange.evaluation.testing.TEST_TYPE_SINGLE, labels=res.labels)
2689    for te in res.results:
2690        for i, ac in enumerate(te.actual_class):
2691            te2 = Orange.evaluation.testing.TestedExample(
2692                iteration_number=te.iteration_number, actual_class=ac)
2693            for c, p in zip(te.classes, te.probabilities):
2694                te2.add_result(c[i], p[i])
2695            res2.results.append(te2)
2696    return score(res2)
2697
2698
2699################################################################################
2700if __name__ == "__main__":
2701    avranks =  [3.143, 2.000, 2.893, 1.964]
2702    names = ["prva", "druga", "tretja", "cetrta" ]
2703    cd = compute_CD(avranks, 14)
2704    #cd = compute_CD(avranks, 10, type="bonferroni-dunn")
2705    print cd
2706
2707    print compute_friedman(avranks, 14)
2708
2709    #graph_ranks("test.eps", avranks, names, cd=cd, cdmethod=0, width=6, textspace=1.5)
Note: See TracBrowser for help on using the repository browser.