fgsea/src/esCalculation.cpp at 9a06694dfc7b54a0a698061a97db15945ede725c · alserglab/fgsea · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
#include "esCalculation.h"
#include <algorithm>

score_t calcES(const vector<int64_t> &ranks, const vector<int> &p, int64_t NS) {
    // p must be sorted
    int n = (int) ranks.size();
    int k = (int) p.size();
    score_t res{NS, 0, n - k, 0};
    score_t cur{NS, 0, n - k, 0};
    int last = -1;
    for (int pos : p) {
        cur.coef_const += pos - last - 1;
        if (res.abs() < cur.abs()) {
            res = cur;
        }
        cur.coef_NS += ranks[pos];
        if (res.abs() < cur.abs()) {
            res = cur;
        }
        last = pos;
    }
    return res;
}

score_t calcES(const vector<int64_t> &ranks, const vector<int> &p) {
    // p must be sorted
    int64_t NS = 0;
    for (int pos : p) {
        NS += ranks[pos];
    }
    return calcES(ranks, p, NS);
}

score_t calcPositiveES(const vector<int64_t> &ranks, const vector<int> &p, int64_t NS) {
    // p must be sorted
    int n = (int) ranks.size();
    int k = (int) p.size();
    score_t res{NS, 0, n - k, 0};
    score_t cur{NS, 0, n - k, 0};
    int last = -1;
    for (int pos : p) {
        cur.coef_NS += ranks[pos];
        cur.coef_const += pos - last - 1;
        res = max(res, cur);
        last = pos;
    }
    return res;
}

score_t calcPositiveES(const vector<int64_t> &ranks, const vector<int> &p) {
    // p must be sorted
    int64_t NS = 0;
    for (int pos : p) {
        NS += ranks[pos];
    }
    return calcPositiveES(ranks, p, NS);
}

bool compareStat(const vector<double> &ranks, const vector<int> &p, double NS, double bound){
    // p must be sorted
    int n = (int) ranks.size();
    int k = (int) p.size();
    double cur = 0.0;
    double q1 = 1.0 / (n - k);
    double q2 = 1.0 / NS;
    int last = -1;
    for (int pos : p) {
        cur += q2 * ranks[pos] - q1 * (pos - last - 1);
        if (cur > bound) {
            return true;
        }
        last = pos;
    }
    return false;
}