KMeans support as range generation strategy for tDL++

ontolearn/learners/tree_learner_refinement_inherit.py

@@ -20,6 +20,7 @@
 import numpy as np
 import pandas as pd
 from sklearn import tree
+from sklearn.cluster import KMeans
 import itertools
 from shap import TreeExplainer
 import matplotlib.pyplot as plt
@@ -32,8 +33,6 @@
 
 logging.getLogger("fitter").setLevel(logging.CRITICAL)
 logger = logging.getLogger(__name__)
-
-
 class TDL_refinement(TDL):
     def __init__(
         self,
@@ -60,6 +59,8 @@ def __init__(
         plot_importance_evo: bool = False,
         refine_iterations: int = 3,
         topk: int = 30,
+        range_strategy: str = "distribution",
+        n_clusters: int = 5,
     ):
         super().__init__(
             knowledge_base,
@@ -86,6 +87,13 @@ def __init__(
         self.plot_importance_evo = plot_importance_evo
         self.refine_iterations = refine_iterations
         self.topk = topk
+        # how numeric data property values are partitioned into ranges: "distribution" or "kmeans"
+        if range_strategy not in ("distribution", "kmeans"):
+            raise ValueError(f"range_strategy must be 'distribution' or 'kmeans', got {range_strategy!r}.")
+        if n_clusters < 1:
+            raise ValueError(f"n_clusters must be >= 1, got {n_clusters}.")
+        self.range_strategy = range_strategy
+        self.n_clusters = n_clusters
         self.concept_per_iteration: List = list()
         self.top_feature_dicts: List = list()
         self.initial_importance_dict: Dict = dict()
@@ -94,31 +102,32 @@ def __init__(
 
     def _pack_data_property_with_range_to_dl_concept(self, property: OWLDataProperty, literal_range: tuple, type) -> str:
         """Repack the ranged data property into an DL Concept String"""
+        facet_restrictions_list: List[OWLFacetRestriction] = []
 
         xsdtype: OWLDatatype = self.data_property_datatype_dict.get(property)
         if xsdtype == OWLDatatype(XSDVocabulary.INTEGER):
-            literal_range = (int(literal_range[0]), int(literal_range[1]))
+            literal_range = tuple(None if b is None else int(b) for b in literal_range)
 
-        min_restr = OWLFacetRestriction(
-            OWLFacet.MIN_INCLUSIVE,
-            OWLLiteral(literal_range[0], xsdtype),
-        )
-        if literal_range[1] is None:
-            max_restr = OWLFacetRestriction(
-                OWLFacet.MAX_INCLUSIVE,
+        if literal_range[0] is not None:
+            min_restr = OWLFacetRestriction(
+                OWLFacet.MIN_INCLUSIVE,
                 OWLLiteral(literal_range[0], xsdtype),
             )
-        else:
+            facet_restrictions_list.append(min_restr)
+
+        if literal_range[1] is not None:
             max_restr = OWLFacetRestriction(
                 OWLFacet.MAX_INCLUSIVE,
                 OWLLiteral(literal_range[1], xsdtype),
             )
-
+            facet_restrictions_list.append(max_restr)
+
+        if not facet_restrictions_list:
+            raise ValueError(f"Cannot create a datatype restriction with no facet restrictions for property {property} and range {literal_range}.")
         dt_range = OWLDatatypeRestriction(
             type_=xsdtype,
-            facet_restrictions=[min_restr, max_restr],
+            facet_restrictions=facet_restrictions_list,
         )
-
         packed_owl_class_expression = OWLDataSomeValuesFrom(property, dt_range)
         return packed_owl_class_expression
 
@@ -133,10 +142,32 @@ def _compute_dt_pdf_ranges(self, best_dist: dict) -> List[tuple]:
         std: np.float64
         mean = frozen_dist.stats(moments="m")
         std = frozen_dist.std()
-        ranges = [(mean - std, mean), (mean, mean + std), (mean - 2 * std, mean - 1 * std), (mean + std, mean + 2 * std), (mean - 3 * std, mean - 2 * std), (mean + 2 * std, mean + 3 * std)]
-
+        #TODO insert also single range <3std and > 3std, if mean unstable return uniform backup ranges
+        ranges = [(mean - std, mean), (mean, mean + std), (mean - 2 * std, mean - 1 * std), (mean + std, mean + 2 * std), (mean - 3 * std, mean - 2 * std), (mean + 2 * std, mean + 3 * std), (None, mean - 3 * std), (mean + 3 * std, None)]
         return ranges
 
+    def _compute_kmeans_ranges(self, dt_values) -> List[tuple]:
+        """Cluster the values with 1D k-means and return the [min, max] interval of each cluster."""
+        values = np.asarray(dt_values, dtype=float)
+        values = values[np.isfinite(values)]
+        unique_values = np.unique(values)
+        if len(unique_values) == 0:
+            return []
+        k = min(self.n_clusters, len(unique_values))
+        labels = KMeans(n_clusters=k, n_init=10, random_state=0).fit_predict(values.reshape(-1, 1))
+        ranges = [(values[labels == c].min(), values[labels == c].max()) for c in range(k)]
+        return sorted(ranges)
+
+    def _compute_ranges(self, dt_values) -> List[tuple]:
+        """Partition data property values into ranges using the configured range_strategy."""
+        if self.range_strategy == "kmeans":
+            return self._compute_kmeans_ranges(dt_values)
+        return self._compute_dt_pdf_ranges(self._find_best_distribution_for_dp(dt_values))
+
+    def _in_interval(self,v, interval) -> bool:
+        lo, hi = interval
+        return (lo is None or lo <= v) and (hi is None or v <= hi)
+
     def _find_best_distribution_for_dp(self, dt_values: List[float], plot: bool = False) -> dict:
         """Iteratively fit distributions on the data and choose the best one based on K-S statistic."""
         assert len(dt_values) > 0, "Data property values list cannot be empty for distribution fitting."
@@ -170,19 +201,21 @@ def _find_best_distribution_for_dp(self, dt_values: List[float], plot: bool = Fa
 
     def _extract_refined_ranges_from_data_properties(self, data_property_values: dict, range_values):
         """Extract refined ranges from data property values based on the original range values."""
+        # one-sided ranges (e.g. [>= v] or [<= v]) are not refined
+        if range_values is None or len(range_values) != 2 or any(b is None for b in range_values):
+            return None
         lb, ub = range_values
         prop_values = np.sort(np.array(data_property_values))
-        # where can maybe be runtime optimized
-        lb_index = np.where(prop_values[np.abs(prop_values - lb).argmin()] == prop_values)
-        ub_index = np.where(prop_values[np.abs(prop_values - ub).argmin()] == prop_values)
-        # if self.verbose>10:
-        logger.debug("lb: %s  ub: %s", lb_index[0][0], ub_index[0][0])
-        if lb_index[0][0] >= ub_index[0][0]:
+        # window = all values v with lb <= v <= ub, including every duplicate of either bound
+        lb_index = np.searchsorted(prop_values, lb, side="left")
+        ub_index = np.searchsorted(prop_values, ub, side="right")
+        window = prop_values[lb_index:ub_index]
+        logger.debug("lb: %s  ub: %s", lb_index, ub_index)
+        if len(np.unique(window)) < 2:
             if self.verbose > 0:
-                logger.warning("Lower bound and upper bound indices are the same. Refinement may not be effective.")
+                logger.warning("Fewer than two distinct values in range [%s, %s]. Refinement may not be effective.", lb, ub)
             return None
-        dist = self._find_best_distribution_for_dp(prop_values[lb_index[0][0] : ub_index[0][0]], plot=False)
-        return self._compute_dt_pdf_ranges(dist)
+        return self._compute_ranges(window)
 
     def _extract_ranges_from_data_properties(
         self,
@@ -196,41 +229,40 @@ def _extract_ranges_from_data_properties(
             return
         ranges_dict: dict[OWLDataProperty, List[Tuple[float, float]]] = dict()
         generated_dt_class_expressions = dict()
-        best_dists = dict()
 
         for prop in data_properties_dict:
-            best_dist = self._find_best_distribution_for_dp(data_properties_dict[prop])
-
-            logger.debug("prop%r", prop)
-            logger.debug(best_dist)
-            best_dists[prop._iri] = best_dist
-            ranges_dict.setdefault(prop, self._compute_dt_pdf_ranges(best_dist))
+            ranges_dict.setdefault(prop, self._compute_ranges(data_properties_dict[prop]))
 
         for r in ranges_dict:
             for interval in range(len(ranges_dict[r])):
-                if (np.isnan(ranges_dict[r][interval][0]) or np.isinf(ranges_dict[r][interval][0])) or (np.isnan(ranges_dict[r][interval][1]) or np.isinf(ranges_dict[r][interval][1])):
+                lo, hi = ranges_dict[r][interval]
+
+                if (lo is  not None and (np.isnan(lo) or np.isinf(lo))) or \
+                   (hi is not None and (np.isnan(hi) or np.isinf(hi))):
                     continue
+                #if (np.isnan(ranges_dict[r][interval][0]) or np.isinf(ranges_dict[r][interval][0])) or (np.isnan(ranges_dict[r][interval][1]) or np.isinf(ranges_dict[r][interval][1])):
+
                 new_class_expression = self._pack_data_property_with_range_to_dl_concept(r, ranges_dict[r][interval], data_properties_dict[r][0])
                 str_dl_concept = owl_expression_to_dl(new_class_expression)
                 # add concept to features
                 if str_dl_concept not in features:
                     features[str_dl_concept] = new_class_expression
+
         # assign concept per individual, if satisfied
-
         for ind, prop in itertools.product(per_individual_data_properties.keys(), data_properties_dict.keys()):
             if prop in per_individual_data_properties.get(ind, {}):
                 for v in per_individual_data_properties[ind][prop]:
                     for interval in ranges_dict[prop]:
-                        if (np.isnan(interval[0]) or np.isinf(interval[0])) or (np.isnan(interval[1]) or np.isinf(interval[1])):
+                        if any(x is not None and (np.isnan(x) or np.isinf(x)) for x in interval):
                             continue
-                        if interval[0] <= v and v <= interval[1]:
+                        if self._in_interval(v, interval):
                             new_class_expression = self._pack_data_property_with_range_to_dl_concept(prop, interval, v)
                             str_dl_concept = owl_expression_to_dl(new_class_expression)
                             individuals_to_feature_mapping[ind.str].add(str_dl_concept)
                             if new_class_expression not in generated_dt_class_expressions:
                                 generated_dt_class_expressions[new_class_expression] = set()
                             generated_dt_class_expressions[new_class_expression].add(ind)
-                            break
+                            #break
         return generated_dt_class_expressions
 
     def _extract_data_property_features(
@@ -435,7 +467,7 @@ def refine_numerical_features(self, topk_expressions: list[OWLClassExpression],
                         continue
 
                     for r in refined_ranges:
-                        if (np.isnan(r[0]) or np.isinf(r[0])) or (np.isnan(r[1]) or np.isinf(r[1])):
+                        if any(x is not None and (np.isnan(x) or np.isinf(x)) for x in r):
                             continue
                         new_class_expression = self._pack_data_property_with_range_to_dl_concept(e.get_property(), r, e.get_filler().get_facet_restrictions()[0].get_facet_value().get_datatype())
                         str_dl_concept = owl_expression_to_dl(new_class_expression)
@@ -452,9 +484,9 @@ def refine_numerical_features(self, topk_expressions: list[OWLClassExpression],
                     for ind in individuals_for_e:
                         for v in self.per_individual_data_properties[ind][e.get_property()]:
                             for interval in refined_ranges:
-                                if (np.isnan(interval[0]) or np.isinf(interval[0])) or (np.isnan(interval[1]) or np.isinf(interval[1])):
+                                if any(x is not None and (np.isnan(x) or np.isinf(x)) for x in interval):
                                     continue
-                                if interval[0] <= v <= interval[1]:
+                                if self._in_interval(v, interval):
                                     new_class_expression = self._pack_data_property_with_range_to_dl_concept(
                                         e.get_property(), interval, e.get_filler().get_facet_restrictions()[0].get_facet_value().get_datatype()
                                     )
@@ -676,13 +708,12 @@ def fit(self, learning_problem: PosNegLPStandard = None, max_runtime: int = None
             X = self._refine_features_with_shap(X, y, learning_problem)
 
         self._report_and_plot(X, y)
-
         perm_feature_topk = False
         if perm_feature_topk:
             # plot permutation feature importance
             res = permutation_importance(self.clf, X.values, y.values, n_repeats=10, random_state=0, n_jobs=12)
             self.plot_perm_feature_importances(res, [owl_expression_to_dl(f) for f in self.features])
-        shap_feature_topk = False
+        shap_feature_topk = True
         # calculate SHAP values, to get topk features
         if shap_feature_topk:
             tree_explainer = TreeExplainer(self.clf)
@@ -692,7 +723,7 @@ def fit(self, learning_problem: PosNegLPStandard = None, max_runtime: int = None
                     global_importance = np.abs(sVal[:, :, 1]).mean(axis=0)
                 else:
                     global_importance = np.abs(sVal).mean(axis=0)
-
+            print([owl_expression_to_dl(f) for f in self.features])
             self.plot_shap_feature_importances(global_importance, [owl_expression_to_dl(f) for f in self.features])
 
         self._finalize_predictions(X, y)

tests/test_tdl_refinement.py

@@ -77,6 +77,29 @@ def test_initialization_custom(self):
         self.assertEqual(model.kwargs_classifier["random_state"], 42)
         self.assertEqual(model.feature_refinement, True)
         self.assertEqual(model.refine_iterations, 3)
+        self.assertEqual(model.range_strategy, "distribution")
+
+    def test_invalid_range_strategy(self):
+        """Test that an unknown range_strategy or non-positive n_clusters is rejected."""
+        with self.assertRaises(ValueError):
+            TDL_refinement(knowledge_base=self.kb_family, range_strategy="quantiles", verbose=0)
+        with self.assertRaises(ValueError):
+            TDL_refinement(knowledge_base=self.kb_family, range_strategy="kmeans", n_clusters=0, verbose=0)
+
+    def test_kmeans_ranges(self):
+        """Test that k-means ranges are sorted, disjoint and cover every value."""
+        model = TDL_refinement(knowledge_base=self.kb_family, range_strategy="kmeans", n_clusters=3, verbose=0)
+        values = [1.0, 1.1, 1.2, 5.0, 5.2, 9.9, 10.0, 10.1]
+        ranges = model._compute_ranges(values)
+        self.assertEqual(ranges, [(1.0, 1.2), (5.0, 5.2), (9.9, 10.1)])
+        for v in values:
+            self.assertEqual(sum(model._in_interval(v, r) for r in ranges), 1)
+
+    def test_kmeans_ranges_fewer_unique_values_than_clusters(self):
+        """Test that k is capped at the number of distinct values."""
+        model = TDL_refinement(knowledge_base=self.kb_family, range_strategy="kmeans", n_clusters=5, verbose=0)
+        self.assertEqual(model._compute_ranges([2.0, 2.0, 7.0]), [(2.0, 2.0), (7.0, 7.0)])
+        self.assertEqual(model._compute_ranges([]), [])
 
 
 class TestTDLFilterLogic(unittest.TestCase):
@@ -763,6 +786,24 @@ def test_inverted_bounds_return_none(self):
             values, (90.0, 10.0))
         self.assertIsNone(result)
 
+    def test_window_includes_upper_bound_and_duplicates(self):
+        """Values equal to ub (incl. duplicates) are refined; values outside [lb, ub] are not."""
+        model = TDL_refinement(knowledge_base=self.kb, range_strategy="kmeans", n_clusters=5, verbose=0)
+        values = [0.0, 1.0, 1.0, 2.0, 3.0, 3.0, 4.0, 4.0, 4.0, 9.0]
+        result = model._extract_refined_ranges_from_data_properties(values, (1.0, 4.0))
+        self.assertEqual(result, [(1.0, 1.0), (2.0, 2.0), (3.0, 3.0), (4.0, 4.0)])
+
+    def test_window_with_two_levels_keeps_both(self):
+        model = TDL_refinement(knowledge_base=self.kb, range_strategy="kmeans", n_clusters=5, verbose=0)
+        result = model._extract_refined_ranges_from_data_properties([2.0, 2.0, 7.0, 7.0], (2.0, 7.0))
+        self.assertEqual(result, [(2.0, 2.0), (7.0, 7.0)])
+
+    def test_bounds_between_data_values_do_not_snap_outside(self):
+        """lb/ub that are not data values must not pull in the nearest value outside the range."""
+        model = TDL_refinement(knowledge_base=self.kb, range_strategy="kmeans", n_clusters=5, verbose=0)
+        result = model._extract_refined_ranges_from_data_properties([1.4, 1.6, 1.8, 2.6], (1.5, 2.5))
+        self.assertEqual(result, [(1.6, 1.6), (1.8, 1.8)])
+
 class TestExtractRangesFromDataProperties(unittest.TestCase):
     """Test range extraction"""
 

tests/test_tdl_refinement_regression.py

@@ -73,6 +73,26 @@ def test_regression_mutagenesis(self):
             print("mutagenesis: ", q)
             assert q >= 0.70
 
+    def test_regression_mutagenesis_kmeans(self):
+        path = "KGs/Mutagenesis/mutagenesis.owl"
+        kb = KnowledgeBase(path=path)
+        with open("LPs/Mutagenesis/lps.json") as json_file:
+            settings = json.load(json_file)
+        model = TDL_refinement(knowledge_base=kb, report_classification=True, kwargs_classifier={"random_state": 1, "max_depth": 3},
+                                     use_inverse= False, use_data_properties_numeric=True, use_data_properties_boolean=False, use_data_properties_string=False, use_data_properties_date=False,
+                                     use_nominals = False, use_card_restrictions = False, feature_refinement=True, refine_iterations=3,
+                                     range_strategy="kmeans", n_clusters=5)
+        for str_target_concept, examples in settings['problems'].items():
+            p = set(examples['positive_examples'])
+            n = set(examples['negative_examples'])
+            typed_pos = set(map(OWLNamedIndividual, map(IRI.create, p)))
+            typed_neg = set(map(OWLNamedIndividual, map(IRI.create, n)))
+            lp = PosNegLPStandard(pos=typed_pos, neg=typed_neg)
+            h = model.fit(learning_problem=lp).best_hypotheses()
+            q = compute_f1_score(individuals=frozenset({i for i in kb.individuals(h)}), pos=lp.pos, neg=lp.neg)
+            print("mutagenesis kmeans: ", q)
+            assert q >= 0.70
+
     def test_regression_carcinogenesis(self):
         path = "KGs/Carcinogenesis/carcinogenesis.owl"
         # (1) Load a knowledge graph.