Add pounders filtering to tranquilo

.pre-commit-config.yaml

@@ -6,7 +6,7 @@ repos:
     -   id: debug-statements
     -   id: end-of-file-fixer
 -   repo: https://github.com/asottile/reorder_python_imports
-    rev: v3.8.3
+    rev: v3.8.5
     hooks:
     -   id: reorder-python-imports
         types: [python]
@@ -47,7 +47,7 @@ repos:
         additional_dependencies: [black==22.3.0]
         exclude: docs/source/how_to_guides/optimization/how_to_specify_constraints.md
 -   repo: https://github.com/psf/black
-    rev: 22.8.0
+    rev: 22.10.0
     hooks:
       - id: black
         language_version: python3.9
@@ -72,7 +72,7 @@ repos:
     -   id: check-useless-excludes
     # -   id: identity  # Prints all files passed to pre-commits. Debugging.
 -   repo: https://github.com/nbQA-dev/nbQA
-    rev: 1.5.2
+    rev: 1.5.3
     hooks:
     -   id: nbqa-black
         additional_dependencies: [black==22.3.0]

src/estimagic/benchmarking/process_benchmark_results.py

@@ -167,7 +167,7 @@ def _make_history_monotone(df, target_col, direction="minimize"):
 
     """
     sorted_df = df.sort_values(["problem", "algorithm", "n_evaluations"])
-    grouped = sorted_df.groupby(["problem", "algorithm"])[target_col]
+    grouped = sorted_df.groupby(["problem", "algorithm"], group_keys=False)[target_col]
 
     if direction == "minimize":
         out = grouped.apply(np.minimum.accumulate)

src/estimagic/optimization/tranquilo/aggregate_models.py

@@ -46,7 +46,7 @@ def get_aggregator(aggregator, functype, model_info):
 
     # determine if aggregator is compatible with functype and model_info
     aggregator_compatible_with_functype = {
-        "scalar": ("identity", "sum"),
+        "scalar": ("identity", "identity_old", "sum"),
         "least_squares": ("least_squares_linear",),
         "likelihood": (
             "sum",
@@ -57,13 +57,14 @@ def get_aggregator(aggregator, functype, model_info):
     aggregator_compatible_with_model_info = {
         # keys are names of aggregators and values are functions of model_info that
         # return False in case of incompatibility
-        "identity": _is_second_order_model,
+        "identity": lambda x: True,
         "sum": _is_second_order_model,
         "information_equality_linear": lambda model_info: not _is_second_order_model(
             model_info
         ),
-        "least_squares_linear": lambda model_info: model_info.has_intercepts
-        and not _is_second_order_model(model_info),
+        "least_squares_linear": lambda model_info: not _is_second_order_model(
+            model_info
+        ),
     }
 
     if _using_built_in_aggregator:
@@ -73,11 +74,6 @@ def get_aggregator(aggregator, functype, model_info):
                 f"Aggregator {_aggregator_name} is not compatible with functype "
                 f"{functype}. It would not produce a quadratic main model."
             )
-        if functype == "scalar" and not _is_second_order_model(model_info):
-            raise ValueError(
-                f"ModelInfo {model_info} is not compatible with functype scalar. "
-                "It would not produce a quadratic main model."
-            )
         if not aggregator_compatible_with_model_info[_aggregator_name](model_info):
             raise ValueError(
                 f"ModelInfo {model_info} is not compatible with aggregator "
@@ -137,9 +133,12 @@ def aggregator_identity(vector_model, fvec_center, model_info):
     2. ModelInfo: has squares or interactions
 
     """
-    intercept = float(fvec_center)
+    intercept = float(vector_model.intercepts)
     linear_terms = np.squeeze(vector_model.linear_terms)
-    square_terms = np.squeeze(vector_model.square_terms)
+    if model_info.has_squares or model_info.has_interactions:
+        square_terms = np.squeeze(vector_model.square_terms)
+    else:
+        square_terms = np.zeros((len(linear_terms), len(linear_terms)))
     return intercept, linear_terms, square_terms
 
 
@@ -157,11 +156,7 @@ def aggregator_sum(vector_model, fvec_center, model_info):
     2. ModelInfo: has squares or interactions
 
     """
-    if model_info.has_intercepts:
-        vm_intercepts = vector_model.intercepts
-    else:
-        vm_intercepts = fvec_center
-
+    vm_intercepts = vector_model.intercepts
     intercept = vm_intercepts.sum(axis=0)
     linear_terms = vector_model.linear_terms.sum(axis=0)
     square_terms = vector_model.square_terms.sum(axis=0)
@@ -210,10 +205,7 @@ def aggregator_information_equality_linear(vector_model, fvec_center, model_info
 
     """
     vm_linear_terms = vector_model.linear_terms
-    if model_info.has_intercepts:
-        vm_intercepts = vector_model.intercepts
-    else:
-        vm_intercepts = fvec_center
+    vm_intercepts = vector_model.intercepts
 
     fisher_information = vm_linear_terms.T @ vm_linear_terms
 

src/estimagic/optimization/tranquilo/filter_points.py

@@ -1,4 +1,6 @@
 import numpy as np
+from numba import njit
+from scipy.linalg import qr_multiply
 
 
 def get_sample_filter(sample_filter="keep_all"):
@@ -18,9 +20,10 @@ def get_sample_filter(sample_filter="keep_all"):
 
     """
     built_in_filters = {
-        "discard_all": _discard_all,
-        "keep_all": _keep_all,
-        "drop_collinear": _drop_collinear,
+        "discard_all": discard_all,
+        "keep_all": keep_all,
+        "drop_collinear": drop_collinear,
+        "drop_pounders": drop_collinear_pounders,
     }
 
     if isinstance(sample_filter, str) and sample_filter in built_in_filters:
@@ -33,14 +36,129 @@ def get_sample_filter(sample_filter="keep_all"):
     return out
 
 
-def _discard_all(xs, indices, state):
+def discard_all(xs, indices, state):
     return state.x.reshape(1, -1), np.array([state.index])
 
 
-def _keep_all(xs, indices, state):
+def keep_all(xs, indices, state):
     return xs, indices
 
 
-def _drop_collinear(xs, indices, state):
+def drop_collinear(xs, indices, state):
     """Make sure that the points that are kept are linearly independent."""
     raise NotImplementedError()
+
+
+def drop_collinear_pounders(xs, indices, state):
+    """Drop collinear points using pounders filtering."""
+    if xs.shape[0] <= xs.shape[1] + 1:
+        filtered_xs, filtered_indices = xs, indices
+    else:
+        filtered_xs, filtered_indices = _drop_collinear_pounders(xs, indices, state)
+
+    return filtered_xs, filtered_indices
+
+
+def _drop_collinear_pounders(xs, indices, state):
+    theta2 = 1e-4
+    n_samples, n_params = xs.shape
+    n_poly_terms = n_params * (n_params + 1) // 2
+
+    indices_reverse = indices[::-1]
+    indexer_reverse = np.arange(n_samples)[::-1]
+
+    radius = state.radius
+    center = state.x
+    index_center = int(np.where(indices_reverse == state.index)[0])
+    centered_xs = (xs - center) / radius
+
+    (
+        linear_features,
+        square_features,
+        idx_list_n_plus_1,
+        index,
+    ) = _get_polynomial_feature_matrices(
+        centered_xs,
+        indexer_reverse,
+        index_center,
+        n_params,
+        n_samples,
+        n_poly_terms,
+    )
+
+    indexer_filtered = indexer_reverse[idx_list_n_plus_1].tolist()
+    _index_center = indexer_reverse[index_center]
+
+    counter = n_params + 1
+
+    while (counter < n_samples) and (index >= 0):
+        if index == _index_center:
+            index -= 1
+            continue
+
+        linear_features[counter, 1:] = centered_xs[index]
+        square_features[counter, :] = _square_features(linear_features[counter, 1:])
+
+        linear_features_pad = np.zeros((n_samples, n_samples))
+        linear_features_pad[:n_samples, : n_params + 1] = linear_features
+
+        n_z_mat, _ = qr_multiply(
+            linear_features_pad[: counter + 1, :],
+            square_features.T[:n_poly_terms, : counter + 1],
+        )
+        beta = np.linalg.svd(n_z_mat.T[n_params + 1 :], compute_uv=False)
+
+        if beta[min(counter - n_params, n_poly_terms) - 1] > theta2:
+            indexer_filtered += [index]
+            counter += 1
+
+        index -= 1
+
+    filtered_indices = indices[indexer_filtered]
+    filtered_xs = xs[indexer_filtered]
+
+    return filtered_xs, filtered_indices
+
+
+def _get_polynomial_feature_matrices(
+    centered_xs, indexer, index_center, n_params, n_samples, n_poly_terms
+):
+    linear_features = np.zeros((n_samples, n_params + 1))
+    square_features = np.zeros((n_samples, n_poly_terms))
+
+    linear_features[0, 1:] = centered_xs[indexer[index_center]]
+    square_features[0, :] = _square_features(linear_features[0, 1:]).flatten()
+
+    _is_center_in_head = index_center < n_params
+    idx_list_n = [i for i in range(n_params + _is_center_in_head) if i != index_center]
+    idx_list_n_plus_1 = [index_center] + idx_list_n
+
+    linear_features[:, 0] = 1
+    linear_features[: n_params + 1, 1:] = centered_xs[indexer[idx_list_n_plus_1]]
+    square_features[: n_params + 1, :] = _square_features(
+        linear_features[: n_params + 1, 1:]
+    )
+
+    idx = n_samples - _is_center_in_head - len(idx_list_n) - 1
+
+    return linear_features, square_features, idx_list_n_plus_1, idx
+
+
+@njit
+def _square_features(x):
+    n_samples, n_params = np.atleast_2d(x).shape
+    n_poly_terms = n_params * (n_params + 1) // 2
+
+    poly_terms = np.empty((n_poly_terms, n_samples), x.dtype)
+    xt = x.T
+
+    idx = 0
+    for i in range(n_params):
+        poly_terms[idx] = 0.5 * xt[i] ** 2
+        idx += 1
+
+        for j in range(i + 1, n_params):
+            poly_terms[idx] = xt[i] * xt[j] / np.sqrt(2)
+            idx += 1
+
+    return poly_terms.T