Feature/visualize wc

PEPit/pep.py

@@ -45,6 +45,14 @@ class PEP(object):
 
         G_value (ndarray): the value of the Gram matrix G that the solver found.
         F_value (ndarray): the value of the vector of :class:`Expression`s F that the solver found.
+
+        trim_dim (bool): trims the apperent useless dimensions of the found worst-case function
+                         (best used with dimension_reduction_heuristic)
+
+        toggled_dimensions (int): number of dimensions output when using eval() on :class:`Point` objects
+                                 (after solving the PEP)
+
+        estimated_dimension (int): estimated dimension of computed worst-case instance.
 
         residual (ndarray): the dual value found by the solver to the lmi constraints G >> 0.
 
@@ -101,6 +109,10 @@ def __init__(self):
         # The constraint G >= 0 is the only constraint that is not defined from the class Constraint.
         # Its dual value, called residual, is then stored in the following attribute.
         self.residual = None
+
+        self.trim_dim = False
+        self.toggled_dimensions = None
+        self.estimated_dimension = None
 
     @staticmethod
     def _reset_classes():
@@ -575,6 +587,7 @@ def _solve_with_wrapper(self, wrapper, verbose=1, return_primal_or_dual="dual",
         # leading to different dual values. The ones we store here provide the proof of the obtained guarantee.
         self.residual = wrapper.assign_dual_values()
         G_value, F_value = wrapper.get_primal_variables()
+        nb_eigenvalues = G_value.shape[0]
 
         # Perform a dimension reduction if required
         if dimension_reduction_heuristic:
@@ -637,6 +650,7 @@ def _solve_with_wrapper(self, wrapper, verbose=1, return_primal_or_dual="dual",
         self.F_value = F_value
         self._eval_points_and_function_values(F_value, G_value, verbose=verbose)
         dual_objective = self.check_feasibility(wc_value, verbose=verbose)
+        self.estimated_dimension = nb_eigenvalues
 
         # Return the value of the minimal performance metric
         if return_primal_or_dual == "dual":
@@ -646,6 +660,25 @@ def _solve_with_wrapper(self, wrapper, verbose=1, return_primal_or_dual="dual",
         else:
             raise ValueError("The argument \'return_primal_or_dual\' must be \'dual\' or \`primal\`."
                              "Got {}".format(return_primal_or_dual))
+
+    def trim_dimension(self, trim_dim=False, toggled_dimensions=None):
+        """
+        Activate the dimension trimmer (convenient for plotting worst-case trajectories).
+	`toggled_dimensions` dimensions are output when using `Point.eval`.
+
+        Args:
+            trim_dim (bool): activates/deactivates the dimension trimmer.
+            toggled_dimensions (int, optional): number of dimensions that are kept when evaluating.
+                                                default to the estimated number of nonzero eigenvalues.
+
+        """
+
+        self.trim_dim = trim_dim
+        if toggled_dimensions is not None:
+            self.toggled_dimensions = min(toggled_dimensions, Point.counter)
+        elif self.estimated_dimension is not None:
+            self.toggled_dimensions = min(self.estimated_dimension, Point.counter)
+        Point._toggled_dimensions = self.toggled_dimensions
 
     def check_feasibility(self, wc_value, verbose=1):
         """
@@ -920,3 +953,12 @@ def _eval_points_and_function_values(self, F_value, G_value, verbose=1):
                                 raise TypeError(
                                     "Expressions are made of function values, inner products and constants only!"
                                     "Got {}".format(type(sub_expression)))
+
+        # Adapts to the trim_dimensions and toggled_dimensions
+        if not self.trim_dim:
+            toggled_dimensions = Point.counter
+        else:
+            if self.toggled_dimensions is not None:
+                Point._toggled_dimensions = min(self.toggled_dimensions, nb_points)
+            else:
+                Point._toggled_dimensions = min(self.estimated_dimension, nb_points)

PEPit/point.py

@@ -20,6 +20,7 @@ class Point(object):
                                    Keys are :class:`Point` objects.
                                    And values are their associated coefficients.
         counter (int): counts the number of leaf :class:`Point` objects.
+        _toggled_dimensions (int): shows only the first `toggled_dimensions` dimensions when using `eval`.
 
     :class:`Point` objects can be added or subtracted together.
     They can also be multiplied and divided by a scalar value.
@@ -52,6 +53,9 @@ class Point(object):
     # namely the number of points needed to linearly generate the others.
     counter = 0
     list_of_leaf_points = list()
+    # Number of dimensions to be output when evaluating (None corresponds to
+    # all dimensions being output
+    _toggled_dimensions = None
 
     def __init__(self,
                  is_leaf=True,
@@ -282,7 +286,6 @@ def eval(self):
             ValueError("The PEP must be solved to evaluate Points!") if the PEP has not been solved yet.
 
         """
-
         # If the attribute value is not None, then simply return it.
         # Otherwise, compute it and return it.
         if self._value is None:
@@ -298,6 +301,23 @@ def eval(self):
 
         return self._value
 
+    def eval_ld(self):
+        """
+        Evaluation in lower dimension.
+
+        Returns:
+            value (np.array): The trimmed value of this :class:`Point` after the corresponding PEP was solved numerically.
+
+        Raises:
+            ValueError("The PEP must be solved to evaluate Points!") if the PEP has not been solved yet.
+
+        """
+        if Point._toggled_dimensions is not None:
+            toggled_dimensions = Point._toggled_dimensions
+        else:
+            toggled_dimensions = Point.counter
+
+        return self.eval()[:toggled_dimensions]
 
 # Define a null Point initialized to 0.
 null_point = Point(is_leaf=False, decomposition_dict=dict())

PEPit/tools/interpolator.py

@@ -0,0 +1,59 @@
+import cvxpy as cp
+import numpy as np
+
+class Interpolator(object):
+    """
+    The class :class:`Interpolator` is designed to help identifying worst-case examples.
+
+    Attributes:
+        list_of_triplets (list): list of triplets
+        
+        L (float): curvature of the quadratic upper bound on the function (possibly np.inf)
+        m (float): curvature of the quadratic lower bound on the function (possibly np.inf)
+        d (int): dimension
+        
+    """
+    def __init__(self, list_of_triplets, mu, L, d, options='lowest'):
+        self.x_list, self.g_list, self.f_list =  list_of_triplets
+        self.nb_eval = len(self.x_list)
+        self.L = L
+        self.mu = mu
+        self.d = d
+        self.options = options
+
+    def __set_constraint__(self,
+                           xi, gi, fi,
+                           xj, gj, fj,
+                           ):
+        if self.L is not np.inf:
+            constraint = (0 >= fj - fi + gj @ (xi - xj)
+                          + 1 / 2 / (self.L-self.mu) * cp.norm(gi - gj,2) ** 2
+                          + self.mu * self.L / 2 / (self.L-self.mu) *  cp.norm(xi - xj,2) ** 2
+                          - self.mu / (self.L - self.mu) * (gi-gj)@(xi-xj) )
+        else:
+            constraint = (0 >= fj - fi + gj @ (xi - xj)
+                          + self.mu / 2 *  cp.norm(xi - xj,2) ** 2 )
+
+        return constraint
+
+
+    def evaluate(self, x):
+        fx = cp.Variable(1)
+        gx = cp.Variable((self.d,))
+        cons = []
+        for i in range(self.nb_eval):
+            fi = self.f_list[i]
+            gi = self.g_list[i]
+            xi = self.x_list[i]
+            cons.append(self.__set_constraint__(xi,gi,fi,x,gx,fx))
+            cons.append(self.__set_constraint__(x,gx,fx,xi,gi,fi))
+
+        if self.options == 'highest':
+            prob = cp.Problem(cp.Maximize(fx), cons)
+        if self.options == 'lowest':
+            prob = cp.Problem(cp.Minimize(fx), cons)
+        prob.solve(verbose=False)
+        return fx.value
+
+
+