Aluode/PerceptionLabPortable

1# Authors: The scikit-learn developers
2# SPDX-License-Identifier: BSD-3-Clause
3
4# Uses the pool adjacent violators algorithm (PAVA), with the
5# enhancement of searching for the longest decreasing subsequence to
6# pool at each step.
7
8import numpy as np
9from cython cimport floating
10
11
12def _inplace_contiguous_isotonic_regression(floating[::1] y, floating[::1] w):
13    cdef:
14        Py_ssize_t n = y.shape[0], i, k
15        floating prev_y, sum_wy, sum_w
16        Py_ssize_t[::1] target = np.arange(n, dtype=np.intp)
17
18    # target describes a list of blocks.  At any time, if [i..j] (inclusive) is
19    # an active block, then target[i] := j and target[j] := i.
20
21    # For "active" indices (block starts):
22    # w[i] := sum{w_orig[j], j=[i..target[i]]}
23    # y[i] := sum{y_orig[j]*w_orig[j], j=[i..target[i]]} / w[i]
24
25    with nogil:
26        i = 0
27        while i < n:
28            k = target[i] + 1
29            if k == n:
30                break
31            if y[i] < y[k]:
32                i = k
33                continue
34            sum_wy = w[i] * y[i]
35            sum_w = w[i]
36            while True:
37                # We are within a decreasing subsequence.
38                prev_y = y[k]
39                sum_wy += w[k] * y[k]
40                sum_w += w[k]
41                k = target[k] + 1
42                if k == n or prev_y < y[k]:
43                    # Non-singleton decreasing subsequence is finished,
44                    # update first entry.
45                    y[i] = sum_wy / sum_w
46                    w[i] = sum_w
47                    target[i] = k - 1
48                    target[k - 1] = i
49                    if i > 0:
50                        # Backtrack if we can.  This makes the algorithm
51                        # single-pass and ensures O(n) complexity.
52                        i = target[i - 1]
53                    # Otherwise, restart from the same point.
54                    break
55        # Reconstruct the solution.
56        i = 0
57        while i < n:
58            k = target[i] + 1
59            y[i + 1 : k] = y[i]
60            i = k
61
62
63def _make_unique(const floating[::1] X,
64                 const floating[::1] y,
65                 const floating[::1] sample_weights):
66    """Average targets for duplicate X, drop duplicates.
67
68    Aggregates duplicate X values into a single X value where
69    the target y is a (sample_weighted) average of the individual
70    targets.
71
72    Assumes that X is ordered, so that all duplicates follow each other.
73    """
74    unique_values = len(np.unique(X))
75
76    if floating is float:
77        dtype = np.float32
78    else:
79        dtype = np.float64
80
81    cdef floating[::1] y_out = np.empty(unique_values, dtype=dtype)
82    cdef floating[::1] x_out = np.empty_like(y_out)
83    cdef floating[::1] weights_out = np.empty_like(y_out)
84
85    cdef floating current_x = X[0]
86    cdef floating current_y = 0
87    cdef floating current_weight = 0
88    cdef int i = 0
89    cdef int j
90    cdef floating x
91    cdef int n_samples = len(X)
92    cdef floating eps = np.finfo(dtype).resolution
93
94    for j in range(n_samples):
95        x = X[j]
96        if x - current_x >= eps:
97            # next unique value
98            x_out[i] = current_x
99            weights_out[i] = current_weight
100            y_out[i] = current_y / current_weight
101            i += 1
102            current_x = x
103            current_weight = sample_weights[j]
104            current_y = y[j] * sample_weights[j]
105        else:
106            current_weight += sample_weights[j]
107            current_y += y[j] * sample_weights[j]
108
109    x_out[i] = current_x
110    weights_out[i] = current_weight
111    y_out[i] = current_y / current_weight
112    return(
113        np.asarray(x_out[:i+1]),
114        np.asarray(y_out[:i+1]),
115        np.asarray(weights_out[:i+1]),
116    )
117