kneeliverse.evaluation

  1# coding: utf-8
  2
  3'''
  4The following module provides a set of methods
  5used for evaluation knee detection algorithms.
  6'''
  7
  8__author__ = 'Mário Antunes'
  9__version__ = '1.0'
 10__email__ = 'mario.antunes@ua.pt'
 11__status__ = 'Development'
 12__license__ = 'MIT'
 13__copyright__ = '''
 14Copyright (c) 2021-2023 Stony Brook University
 15Copyright (c) 2021-2023 The Research Foundation of SUNY
 16'''
 17
 18import enum
 19import math
 20import logging
 21import numpy as np
 22import kneeliverse.linear_fit as lf
 23import kneeliverse.metrics as metrics
 24
 25
 26logger = logging.getLogger(__name__)
 27
 28
 29class Strategy(enum.Enum):
 30    """
 31    Enum data type that represents the strategy of MAE, MSE and RMSE
 32    """
 33    knees = 'knees'
 34    expected = 'expected'
 35    best = 'best'
 36    worst = 'worst'
 37
 38    def __str__(self):
 39        return self.value
 40
 41
 42def get_neighbourhood_points(points: np.ndarray, a: int, b: int, t: float) -> tuple:
 43    """
 44    Get the neighbourhood (closest points) from a to b.
 45
 46    The neighbourhood is defined as the longest straitgh line (defined by $R^2$).
 47
 48    Args:
 49        points (np.ndarray): numpy array with the points (x, y)
 50        a (int): the initial point of the search
 51        b (int): the left limit of the search
 52        t (float): $R^2$ threshold
 53
 54    Returns:
 55        tuple: (neighbourhood index, r2, slope)
 56    """
 57
 58    x = points[:, 0]
 59    y = points[:, 1]
 60    return get_neighbourhood(x, y, a, b, t)
 61
 62
 63def get_neighbourhood_fast_points(points: np.ndarray, a: int, b: int, t: float) -> tuple:
 64    """
 65    Get the neighbourhood (closest points) from a to b.
 66
 67    The neighbourhood is defined as the longest straitgh line (defined by $R^2$).
 68    This version uses a inaccurate binary search to speedup the search.
 69
 70    Args:
 71        points (np.ndarray): numpy array with the points (x, y)
 72        a (int): the initial point of the search
 73        b (int): the left limit of the search
 74        t (float): $R^2$ threshold
 75
 76    Returns:
 77        tuple: (neighbourhood index, r2, slope)
 78    """
 79
 80    x = points[:, 0]
 81    y = points[:, 1]
 82    return get_neighbourhood_fast(x, y, a, b, t)
 83
 84
 85def get_neighbourhood_binary(x: np.ndarray, y: np.ndarray, a: int, b: int, t=0.9) -> int:
 86    """
 87    Get the index of the point within the range $[b, a]$ where the $R^2$ is close to the threshold.
 88
 89    This version uses a inaccurate binary search to speedup the search.
 90
 91    Args:
 92        x (np.ndarray): the value of the points in the x axis coordinates
 93        y (np.ndarray): the value of the points in the y axis coordinates
 94        a (int): the initial point of the search
 95        b (int): the left limit of the search
 96        t (float): $R^2$ threshold (default 0.9)
 97
 98    Returns:
 99        int: index of the point
100    """
101
102    i = b
103    right = a
104
105    while abs(i-right) > 1:
106        coef = lf.linear_fit(x[i:a+1], y[i:a+1])
107        r2 = lf.linear_r2(x[i:a+1], y[i:a+1], coef)
108
109        if r2 < t:
110            i = int((i+right)/2.0)
111        else:
112            right = i
113            i = int((b+right)/2.0)
114
115    return i
116
117
118def get_neighbourhood_fast(x: np.ndarray, y: np.ndarray, a: int, b: int, t: float = 0.9) -> tuple:
119    """
120    Get the neighbourhood (closest points) from a to b.
121
122    The neighbourhood is defined as the longest straitgh line (defined by $R^2$).
123    This version uses a inaccurate binary search to speedup the search.
124
125    Args:
126        x (np.ndarray): the value of the points in the x axis coordinates
127        y (np.ndarray): the value of the points in the y axis coordinates
128        a (int): the initial point of the search
129        b (int): the left limit of the search
130        t (float): $R^2$ threshold (default 0.9)
131
132    Returns:
133        tuple: (neighbourhood index, r2, slope)
134    """
135    # speedup when the search using an inaccurate binary search
136    i = get_neighbourhood_binary(x, y, a, b, t)
137    b, slope = lf.linear_fit(x[i:a+1], y[i:a+1])
138    r2 = lf.linear_r2(x[i:a+1], y[i:a+1], (b, slope))
139    previous_res = (i, r2, slope)
140
141    # Linear search to improve accuracy
142    while r2 < t and i < a:
143        i += 1
144        coef = lf.linear_fit(x[i:a+1], y[i:a+1])
145        r2 = lf.linear_r2(x[i:a+1], y[i:a+1], coef)
146        _, slope = coef
147        previous_res = (i, r2, slope)
148
149    return previous_res
150
151
152def get_neighbourhood(x: np.ndarray, y: np.ndarray, a: int, b: int, t: float = 0.9) -> tuple:
153    """
154    Get the neighbourhood (closest points) from a to b.
155
156    The neighbourhood is defined as the longest straitgh line (defined by $R^2$).
157
158    Args:
159        x (np.ndarray): the value of the points in the x axis coordinates
160        y (np.ndarray): the value of the points in the y axis coordinates
161        a (int): the initial point of the search
162        b (int): the left limit of the search
163        t (float): $R^2$ threshold (default 0.9)
164
165    Returns:
166        tuple: (neighbourhood index, r2, slope)
167    """
168
169    r2 = 1.0
170    i = a - 1
171    _, slope = lf.linear_fit(x[i:a+1], y[i:a+1])
172
173    while r2 > t and i > b:
174        # print('.')
175        previous_res = (i, r2, slope)
176        i -= 1
177        coef = lf.linear_fit(x[i:a+1], y[i:a+1])
178        r2 = lf.linear_r2(x[i:a+1], y[i:a+1], coef)
179        _, slope = coef
180        #print(f'{i} -> {r2}')
181
182    if r2 > t:
183        return i, r2, slope
184    else:
185        return previous_res
186
187
188def accuracy_knee(points: np.ndarray, knees: np.ndarray, t: float = 0.9) -> tuple:
189    """Compute the accuracy heuristic for a set of knees.
190
191    The heuristic is based on the average distance of X and Y axis, the slope and the $R^2$.
192    In this version it is used the left neighbourhood of the knee.
193
194    Args:
195        points (np.ndarray): numpy array with the points (x, y)
196        knees (np.ndarray): knees indexes
197        t (float): $R^2$ threshold (default 0.9)
198
199    Returns:
200        tuple: (average_x, average_y, average_slope, average_coeffients, cost)
201    """
202    x = points[:, 0]
203    y = points[:, 1]
204
205    total_x = math.fabs(x[-1] - x[0])
206    total_y = math.fabs(y[-1] - y[0])
207
208    distances_x = []
209    distances_y = []
210    slopes = []
211    coeffients = []
212
213    previous_knee = 0
214    for i in range(len(knees)):
215        idx, r2, slope = get_neighbourhood_fast(x, y, knees[i], previous_knee)
216
217        delta_x = x[idx] - x[knees[i]]
218        delta_y = y[idx] - y[knees[i]]
219
220        distances_x.append(math.fabs(delta_x))
221        distances_y.append(math.fabs(delta_y))
222        slopes.append(math.fabs(slope))
223        coeffients.append(r2)
224
225        previous_knee = knees[i]
226
227    slopes = np.array(slopes)
228    slopes = slopes/slopes.max()
229
230    coeffients = np.array(coeffients)
231    coeffients = coeffients/coeffients.max()
232
233    distances_x = np.array(distances_x)/total_x
234    distances_y = np.array(distances_y)/total_y
235    average_x = np.average(distances_x)
236    average_y = np.average(distances_y)
237
238    average_slope = np.average(slopes)
239    average_coeffients = np.average(coeffients)
240
241    #p = slopes * distances_y * coeffients
242    p = slopes * distances_y
243    #cost = (average_x * average_y) / (average_slope)
244    cost = average_x / np.average(p)
245
246    return average_x, average_y, average_slope, average_coeffients, cost
247
248
249def accuracy_trace(points: np.ndarray, knees: np.ndarray) -> tuple:
250    """Compute the accuracy heuristic for a set of knees.
251
252    The heuristic is based on the average distance of X and Y axis, the slope and the $R^2$.
253    In this version it is used the points from the current knee to the previous.
254
255    Args:
256        points (np.ndarray): numpy array with the points (x, y)
257        knees (np.ndarray): knees indexes
258
259    Returns:
260        tuple: (average_x, average_y, average_slope, average_coeffients, cost)
261    """
262    x = points[:, 0]
263    y = points[:, 1]
264
265    distances_x = []
266    distances_y = []
267    slopes = []
268    coeffients = []
269
270    total_x = math.fabs(x[-1] - x[0])
271    total_y = math.fabs(y[-1] - y[0])
272
273    previous_knee_x = x[knees[0]]
274    previous_knee_y = y[knees[0]]
275
276    delta_x = x[0] - previous_knee_x
277    delta_y = y[0] - previous_knee_y
278    distances_x.append(math.fabs(delta_x))
279    distances_y.append(math.fabs(delta_y))
280
281    coef = lf.linear_fit(x[0:knees[0]+1], y[0:knees[0]+1])
282    r2 = lf.linear_r2(x[0:knees[0]+1], y[0:knees[0]+1], coef)
283    coeffients.append(r2)
284    _, slope = coef
285    slopes.append(math.fabs(slope))
286
287    for i in range(1, len(knees)):
288        knee_x = x[knees[i]]
289        knee_y = y[knees[i]]
290
291        delta_x = previous_knee_x - knee_x
292        delta_y = previous_knee_y - knee_y
293
294        coef = lf.linear_fit(x[knees[i-1]:knees[i]+1],
295                             y[knees[i-1]:knees[i]+1])
296        r2 = lf.linear_r2(x[knees[i-1]:knees[i]+1],
297                          y[knees[i-1]:knees[i]+1], coef)
298
299        distances_x.append(math.fabs(delta_x))
300        distances_y.append(math.fabs(delta_y))
301        _, slope = coef
302        slopes.append(math.fabs(slope))
303        coeffients.append(r2)
304
305        previous_knee_x = knee_x
306        previous_knee_y = knee_y
307
308    distances_x = np.array(distances_x)/total_x
309    distances_y = np.array(distances_y)/total_y
310    slopes = np.array(slopes)
311    slopes = slopes/slopes.max()
312
313    coeffients = np.array(coeffients)
314    coeffients = coeffients/coeffients.max()
315
316    coeffients[coeffients < 0] = 0.0
317    p = slopes * distances_y * coeffients
318    #p = slopes * distances_y
319
320    average_x = np.average(distances_x)
321    average_y = np.average(distances_y)
322    average_slope = np.average(slopes)
323    average_coeffients = np.average(coeffients)
324
325    cost = average_x / np.average(p)
326
327    return average_x, average_y, average_slope, average_coeffients, cost
328
329
330def mae(points: np.ndarray, knees: np.ndarray, expected: np.ndarray, s: Strategy = Strategy.expected) -> float:
331    """
332    Estimates the worst case Mean Absolute Error (MAE) for the given
333    knee and expected points.
334
335    Suppports different size arrays, and estimates the MAE based 
336    on the worst case.
337    It uses the euclidean distance to find the closer points,
338    and computes the error based on the closest point.
339
340    Args:
341        points (np.ndarray): numpy array with the points (x, y)
342        knees (np.ndarray): knees indexes
343        expected (np.ndarray): numpy array with the expected knee points (x, y)
344        s (Strategy): enum that controls the point matching (default Strategy.expected)
345
346    Returns:
347        float: the worst case MAE
348    """
349    # get the knee points
350    knee_points = points[knees]
351
352    error = 0.0
353
354    if s is Strategy.knees:
355        a = knee_points
356        b = expected
357    elif s is Strategy.expected:
358        a = expected
359        b = knee_points
360    elif s is Strategy.best:
361        if len(expected) <= len(knee_points):
362            a = expected
363            b = knee_points
364        else:
365            a = knee_points
366            b = expected
367    else:
368        if len(expected) >= len(knee_points):
369            a = expected
370            b = knee_points
371        else:
372            a = knee_points
373            b = expected
374
375    for p in a:
376        distances = np.linalg.norm(b-p, axis=1)
377        idx = np.argmin(distances)
378        error += np.sum(np.abs(p-b[idx]))
379
380    return error / (len(a)*2.0)
381
382
383def mse(points: np.ndarray, knees: np.ndarray, expected: np.ndarray, s: Strategy = Strategy.expected) -> float:
384    """
385    Estimates the worst case Mean Squared Error (MSE) for the given
386    knee and expected points.
387
388    Suppports different size arrays, and estimates the MSE based 
389    on the worst case.
390    It uses the euclidean distance to find the closer points,
391    and computes the error based on the closest point.
392
393    Args:
394        points (np.ndarray): numpy array with the points (x, y)
395        knees (np.ndarray): knees indexes
396        expected (np.ndarray): numpy array with the expected knee points (x, y)
397        s (Strategy): enum that controls the point matching (default Strategy.expected)
398
399    Returns:
400        float: the worst case MSE
401    """
402    # get the knee points
403    knee_points = points[knees]
404
405    error = 0.0
406
407    if s is Strategy.knees:
408        a = knee_points
409        b = expected
410    elif s is Strategy.expected:
411        a = expected
412        b = knee_points
413    elif s is Strategy.best:
414        if len(expected) <= len(knee_points):
415            a = expected
416            b = knee_points
417        else:
418            a = knee_points
419            b = expected
420    else:
421        if len(expected) >= len(knee_points):
422            a = expected
423            b = knee_points
424        else:
425            a = knee_points
426            b = expected
427
428    for p in a:
429        distances = np.linalg.norm(b-p, axis=1)
430        idx = np.argmin(distances)
431        error += np.sum(np.square(p-b[idx]))
432
433    return error / (len(a)*2.0)
434
435
436def rmse(points: np.ndarray, knees: np.ndarray, expected: np.ndarray, s: Strategy = Strategy.expected) -> float:
437    """
438    Estimates the worst case Root Mean Squared Error (RMSE) for the given
439    knee and expected points.
440
441    Suppports different size arrays, and estimates the RMSE based 
442    on the worst case.
443    It uses the euclidean distance to find the closer points,
444    and computes the error based on the closest point.
445
446    Args:
447        points (np.ndarray): numpy array with the points (x, y)
448        knees (np.ndarray): knees indexes
449        expected (np.ndarray): numpy array with the expected knee points (x, y)
450        s (Strategy): enum that controls the point matching (default Strategy.expected)
451
452    Returns:
453        float: the worst case RMSE
454    """
455    return math.sqrt(mse(points, knees, expected, s))
456
457
458def rmspe(points: np.ndarray, knees: np.ndarray, expected: np.ndarray, s: Strategy = Strategy.expected, eps: float = 1e-16) -> float:
459    """
460    Estimates the worst case Root Mean Squared Percentage Error (RMSPE) for the given knee and expected points.
461
462    Suppports different size arrays, and estimates the RMSPE based on the worst case.
463    It uses the euclidean distance to find the closer points, and computes the error based on the closest point.
464
465    Args:
466        points (np.ndarray): numpy array with the points (x, y)
467        knees (np.ndarray): knees indexes
468        expected (np.ndarray): numpy array with the expected knee points (x, y)
469        s (Strategy): enum that controls the point matching (default Strategy.expected)
470        eps (float): eps value to prevent division by zero (default: 1E-16)
471
472    Returns:
473        float: the worst case RMSPE
474    """
475    # get the knee points
476    knee_points = points[knees]
477
478    if s is Strategy.knees:
479        a = knee_points
480        b = expected
481    elif s is Strategy.expected:
482        a = expected
483        b = knee_points
484    elif s is Strategy.best:
485        if len(expected) <= len(knee_points):
486            a = expected
487            b = knee_points
488        else:
489            a = knee_points
490            b = expected
491    else:
492        if len(expected) >= len(knee_points):
493            a = expected
494            b = knee_points
495        else:
496            a = knee_points
497            b = expected
498
499    errors = []
500
501    for p in a:
502        distances = np.linalg.norm(b-p, axis=1)
503        idx = np.argmin(distances)
504        e = (p - b[idx]) / (p + eps)
505        errors.extend(e)
506    errors = np.array(errors)
507
508    return np.sqrt(np.mean(np.square(errors)))
509
510
511def cm(points: np.ndarray, knees: np.ndarray, expected: np.ndarray, t: float = 0.01) -> np.ndarray:
512    """
513    Computes the Confusion Matrix based on the knees and expected points.
514
515    Args:
516        points (np.ndarray): numpy array with the points (x, y)
517        knees (np.ndarray): knees indexes
518        expected (np.ndarray): numpy array with the expected knee points (x, y)
519        t (float): the maximum allowed distance in percentage (default 0.01)
520
521    Returns:
522        np.ndarray: the confusion matrix
523    """
524
525    #dx = math.fabs(points[-1][0] - points[0][0])
526    max_x, _ = points.max(axis=0)
527    min_x, _ = points.min(axis=0)
528    dx = math.fabs(max_x - min_x)
529    used_knees = []
530    knees_points_x = points[knees][:, 0]
531    tp = fn = fp = tn = 0
532    for px, _ in expected:
533        distances = np.fabs(knees_points_x - px)/dx
534        idx = np.argmin(distances)
535        #logger.info(f"{px} / {idx} / {distances[idx]} = {distances[idx] <= t}")
536        if distances[idx] <= t and idx not in used_knees:
537            tp += 1
538            used_knees.append(idx)
539        else:
540            fn += 1
541    fp = max(len(knees) - tp, 0)
542    tn = len(points) - (tp+fp+fn)
543
544    return np.array([[tp, fp], [fn, tn]])
545
546
547def accuracy(cm: np.ndarray) -> float:
548    """
549    Computes accuracy based on a Confusion Matrix.
550
551    Args:
552        cm (np.ndarray): the confusion matrix
553
554    Returns:
555        float: the accuracy
556    """
557
558    tp, fp = cm[0]
559    fn, tn = cm[1]
560
561    return (tp+tn)/(tp+tn+fp+fn)
562
563
564def f1score(cm: np.ndarray) -> float:
565    """
566    Computes F1-Score based on a Confusion Matrix.
567
568    Args:
569        cm (np.ndarray): the confusion matrix
570
571    Returns:
572        float: the F1-Score
573    """
574
575    tp, fp = cm[0]
576    fn, _ = cm[1]
577
578    return (2.0*tp)/(2*tp+fp+fn)
579
580
581def mcc(cm: np.ndarray) -> float:
582    """
583    Computes Matthews Correlation Coefficient (MCC) based on a Confusion Matrix.
584
585    Args:
586        cm (np.ndarray): the confusion matrix
587
588    Returns:
589        float: the mcc
590    """
591
592    tp, fp = cm[0]
593    fn, tn = cm[1]
594
595    n = tp*tn - fp*fn
596    d = math.sqrt((tp+fp)*(tp+fn)*(tn+fp)*(tn+fn))
597
598    return n/d
599
600
601def compute_global_rmse(points: np.ndarray, reduced: np.ndarray, cache:dict=None) -> float:
602    """
603    Computes the global RMSE for a point reduction set.
604
605    For each segment within the point reduction set, computes a straight line.
606    For each line segment it computes the error with the original set.
607 
608    Args:
609        points (np.ndarray): the original points
610        reduced (np.ndarray): the indexes of the reduced points
611        cache (dict): cache used to speedup the computation (default None)
612
613    Returns:
614        float: the global RMSE
615    """
616    # Setup the cache
617    if cache is None: cache = {}
618
619    segment_errors = np.zeros(len(reduced)-1)
620
621    left = reduced[0]
622    for i in range(1, len(reduced)):
623        right = reduced[i]
624
625        # Get data from cache
626        if (left, right) not in cache:
627            pt = points[left:right+1]
628            coef = lf.linear_fit_points(pt)
629            y_hat = lf.linear_transform_points(pt, coef)
630            y = pt[:,1]
631            segment_error = np.sum(np.square((y-y_hat)))
632            cache[(left, right)] = segment_error
633        
634        segment_errors[i-1] = cache[(left, right)]
635        left = right
636
637    # compute the cost function
638    return math.sqrt(np.sum(segment_errors)/len(points))
639
640
641def mip(points: np.ndarray, reduced: np.ndarray) -> tuple:
642    """
643    Computes the median improvement per point (MIP).
644
645    Improvement is defined as a relation between a final cost against a reference cost.
646    The final cost is the RSME with the complete reduced set.
647    The reference cost id different for each point.
648    It is computed as the RMSE of the reconstruction without that specific point. 
649
650    Args:
651        points (np.ndarray): the original points
652        reduced (np.ndarray): the indexes of the reduced points
653
654    Returns:
655        tuple: the median improvement per point (MIP) and the MAD
656    """
657    # Setup the cache
658    cache = {}
659
660    ip = np.zeros(len(reduced)-2)
661
662    # Compute the final RSME
663    cost_fin = compute_global_rmse(points, reduced, cache)
664
665    for i in range(1, len(reduced)-1):
666        # Compute the reference RMSE
667        cost_ref = compute_global_rmse(points, np.delete(reduced, i), cache)
668        ip[i-1] = cost_ref - cost_fin
669
670    mip = np.median(ip)
671
672    return mip, np.median(np.absolute(ip - mip))
673
674
675def compute_cost(points: np.ndarray, segment_errors:np.ndarray, cost:metrics.Metrics, cache:dict) -> float:
676    """
677    Compute the cost of multi-point fitting using the segment errors.
678    
679    It uses a cache to speedup the compuation of the cost.
680    The cache contains the cost of previously computed segments.
681
682    Args:
683        points (np.ndarray): the original points
684        segment_errors (np.ndarray): errors from each of the segments
685        cost (metrics.Metrics): the metric used for the cost calculation
686        cache (dict): the cache used to store the segments costs
687
688    Returns:
689        float: the global cost of a multi-point fitting
690    """
691    # methods = {metrics.Metrics.r2: metrics.r2,metrics.Metrics.rpd: metrics.rpd,metrics.Metrics.rmsle: metrics.rmsle,metrics.Metrics.rmspe: metrics.rmspe,metrics.Metrics.smape: metrics.smape}
692    # cost = methods[cost](np.array(y), np.array(y_hat))
693
694    # We have to account for the head of each segment
695    total = len(points) + len(segment_errors) - 1
696
697    if cost is metrics.Metrics.r2:
698        # Get tss from the cache
699        if 'tss' not in cache:
700            # Compute TSS and add it to cache (single operation)
701            y = points[:,1]
702            y_mean = np.mean(y)
703            tss = np.sum(np.square(y-y_mean))
704            cache['tss'] = tss
705        tss = cache['tss']
706        rss = np.sum(segment_errors)
707        cost = 1.0 - rss if tss == 0 else 1.0 - (rss/tss)
708        #return np.sum(np.square(y-y_hat))
709    elif cost is metrics.Metrics.rmsle:
710        #return np.sum(np.square((np.log(y+1) - np.log(y_hat+1))))
711        cost = math.sqrt(np.sum(segment_errors)/total)
712    elif cost is metrics.Metrics.rmspe:
713        #return np.sum(np.square((y - y_hat) / (y+eps)))
714        cost = math.sqrt(np.sum(segment_errors)/total)
715    elif cost is metrics.Metrics.rpd:
716        #return np.sum(np.abs((y - y_hat) / (np.maximum(y, y_hat)+eps)))
717        cost = np.sum(segment_errors)/total
718    else:
719        #return np.sum(2.0 * np.abs(y_hat - y) / (np.abs(y) + np.abs(y_hat) + eps))
720        cost = np.sum(segment_errors)/total
721
722    cost = 0 if cost < 0 else cost
723
724    return cost
725
726
727def compute_partial_cost(y:np.ndarray, y_hat:np.ndarray, cost: metrics.Metrics, eps: float = 1e-16) -> float:
728    """
729    Compute the partial cost of a multi-point fitting.
730
731    Args:
732        y (np.ndarray): the y values from the points
733        y_hat (np.ndarray): the y values computed with the linear aproximation
734        cost (metrics.Metrics): the metric used for the cost calculation
735        eps (float): eps value to prevent division by zero (default: 1E-16)
736
737    Returns:
738        float: the global cost of a multi-point fitting
739    """
740    if cost is metrics.Metrics.r2:
741        return np.sum(np.square(y-y_hat))
742    elif cost is metrics.Metrics.rmsle:
743        return np.sum(np.square((np.log(y+1) - np.log(y_hat+1))))
744    elif cost is metrics.Metrics.rmspe:
745        return np.sum(np.square((y - y_hat) / (y+eps)))
746    elif cost is metrics.Metrics.rpd:
747        return np.sum(np.abs((y - y_hat) / (np.maximum(y, y_hat)+eps)))
748    else:
749        return np.sum(2.0 * np.abs(y_hat - y) / (np.abs(y) + np.abs(y_hat) + eps))
750
751
752def compute_global_cost(points: np.ndarray, reduced: np.ndarray, cost: metrics.Metrics = metrics.Metrics.rpd, cache:dict=None) -> float:
753    """
754    Compute the cost of multi-point fitting using the segment errors.
755    
756    It uses a cache to speedup the compuation of the cost.
757    The cache contains the cost of previously computed segments.
758
759    Args:
760        points (np.ndarray): the original points
761        reduced (np.ndarray): the reduced set of points used for the fitting
762        cost (metrics.Metrics): the metric used for the cost calculation (default: metrics.Metrics.rpd)
763        cache (dict): the cache used to store the segments costs (default: None)
764
765    Returns:
766        float: the global cost of a multi-point fitting
767    """
768    # Setup the cache
769    if cache is None: cache = {}
770
771    segment_errors = np.zeros(len(reduced)-1)
772
773    left = reduced[0]
774    for i in range(1, len(reduced)):
775        right = reduced[i]
776
777        # Get data from cache
778        if (left, right) not in cache:
779            pt = points[left:right+1]
780            if len(pt) <= 2:
781                cache[(left, right)] = 0
782            else:
783                y_hat = lf.linear_fit_transform_points(pt)
784                y = pt[:,1]
785                segment_error = compute_partial_cost(y, y_hat, cost)
786                cache[(left, right)] = segment_error
787        
788        segment_errors[i-1] = cache[(left, right)]
789        left = right
790
791    return compute_cost(points, segment_errors, cost, cache)
792
793
794def compute_global_segment_cost(points: np.ndarray, reduced: np.ndarray, cost: metrics.Metrics = metrics.Metrics.rpd) -> tuple:
795    """
796    Compute the cost of multi-point fitting using the segment errors.
797    
798    Legacy function that does not use a cache or combines the segments costs.
799    Used for unit testing.
800
801    Args:
802        points (np.ndarray): the original points
803        reduced (np.ndarray): the reduced set of points used for the fitting
804        cost (metrics.Metrics): the metric used for the cost calculation (default: metrics.Metrics.rpd)
805
806    Returns:
807        tuple: the global cost, and partial costs
808    """
809    y, y_hat = [], []
810
811    cost_segment = []
812
813    left = reduced[0]
814    for i in range(1, len(reduced)):
815        right = reduced[i]
816        pt = points[left:right+1]
817        
818        y_temp, y_hat_temp = lf.linear_fit_transform_points(pt)
819        
820        y_hat.extend(y_hat_temp)
821        y.extend(y_temp)
822
823        # compute the cost function
824        #c = compute_cost(y_temp, y_hat_temp, cost)
825        c = lf.linear_hv_residuals_points(pt) 
826        cost_segment.append(c)
827
828        left = right
829
830    # compute the cost function
831    return compute_cost(y, y_hat, cost), np.array(cost_segment)