improved numerical stability with exp

Author: Jouni Helske

hdbscan/_prediction_utils.pyx

@@ -225,8 +225,10 @@ cpdef np.ndarray[np.float64_t, ndim=1] outlier_membership_vector(neighbor,
     if softmax:
         result = per_cluster_scores(neighbor, lambda_, clusters, tree,
                                     max_lambda_dict, cluster_tree)
-        result = np.exp(result)
-        result[~np.isfinite(result)] = np.finfo(np.double).max
+        # Scale for numerical stability, mathematically equivalent with old
+        # version due to the scaling with the sum in below.
+        result = np.exp(result - np.nanmax(result))
+        #result[~np.isfinite(result)] = np.finfo(np.double).max
     else:
         result = per_cluster_scores(neighbor, lambda_, clusters, tree,
                                     max_lambda_dict, cluster_tree)
@@ -310,8 +312,10 @@ cpdef np.ndarray[np.float64_t, ndim=2] all_points_outlier_membership_vector(
                                 max_lambda_dict,
                                 cluster_tree)
     if softmax:
-        result = np.exp(per_cluster_scores)
-        result[~np.isfinite(result)] = np.finfo(np.double).max
+        # Scale for numerical stability, mathematically equivalent with old
+        # version due to the scaling with the sum in below.
+        result = np.exp(per_cluster_scores - np.nanmax(per_cluster_scores))
+        #result[~np.isfinite(result)] = np.finfo(np.double).max
     else:
         result = per_cluster_scores
 
@@ -354,4 +358,3 @@ cpdef all_points_prob_in_some_cluster(
         result[point] = (heights.max() / max_lambda)
 
     return result
-