Computing constraint violations

[millskyle]

I have a problem where the primal solution found by SCIP (using model.getBestSol()) seems to not obey the constraints produced from an ecole.observation.NodeBipartiteObs. My understanding is the code below should produce an output where every element of V is True, since a primal solution found by SCIP is necessarily feasible and thus does not violate any of the constraints (Ax <= b should be true). In this case, x is coming from SCIP, A and b are coming from ecole.

My observation is a custom function returning a model.as_pyscipopt() (that I will then proceed to run .optimize() on; I know this is not legal for the actual competition, and I should make a copy), as well as an ecole.observation.NodeBipartiteObs object.

Does anyone know where I am going wrong? I have a hunch that perhaps the variables in in the ecole.observation.NodeBipartiteObs are not in the same order produced by model.getVars(), but I cannot find any documentation about the ordering of variables in either.

I appreciate any help I can get here

Here is the minimal-working example code. Note that the line
instances_path = Path("./1_item_placement/train/") is user-specific.

from scipy import sparse
import numpy as np
import ecole
import pyscipopt
from pathlib import Path

import torch
import torch.nn.functional as F

class ObservationFunction():
    def __init__(self, problem):
        # called once for each problem benchmark
        self.problem = problem  # to devise problem-specific observations

    def seed(self, seed):
        # called before each episode
        # use this seed to make your code deterministic
        pass

    def before_reset(self, model):
        # called when a new episode is about to start
        pass

    def extract(self, model, done):
        return model.as_pyscipopt()


if __name__=="__main__":
    instances_path = Path("./1_item_placement/train/")
    instances = iter(list(instances_path.glob('*.mps.gz')))

    scip_parameters = {'separating/maxrounds': 0, 'presolving/maxrestarts': 0, 'limits/time': 300}

    env = ecole.environment.PrimalSearch(trials_per_node=500,
                                        depth_freq=1,
                                        depth_start=0, #start on the root node
                                        depth_stop=0, #do not allow it to proceed deeper
                                        observation_function=(ObservationFunction(problem=None),
                                                              ecole.observation.NodeBipartite()),
                                        scip_params=scip_parameters)

    for seed, instance in enumerate(instances):
        observation, action_set, _, done, _ = env.reset(str(instance))
        while not done:

            model, node_bipartite = observation  

            model.setPresolve(pyscipopt.scip.PY_SCIP_PARAMSETTING.OFF) 
            model.setHeuristics(pyscipopt.scip.PY_SCIP_PARAMSETTING.AGGRESSIVE) 

            model.setParam('limits/bestsol', -1)  # stop as many solutions as possible...
            model.setParam('limits/time', 10)     # ...within a time limit of X seconds
            model.optimize()

            assert model.getNSols() > 0  # make sure we found at least one solution

            sol = model.getBestSol()
            sol_vals = np.asarray([model.getSolVal(sol, var) for var in model.getVars(transformed=False)])
            


            b = node_bipartite.row_features[:,0]                 #The constraint objective of a MILP is Ax <= b. 
                                                                 #   These should be the bs as referenced here: 
                                                                 #      https://doc.ecole.ai/py/en/stable/reference/observations.html#ecole.observation.NodeBipartiteObs.RowFeatures
            x = sol_vals                                         #These are the xs, the solution found by SCIP
            _A = node_bipartite.edge_features.values             #These are the entries in the constraint matrix A
            _i = node_bipartite.edge_features.indices            #And these are the indices of these entries

            #Construct the constraint matrix as a sparse matrix: 
            A = sparse.csr_matrix((_A.squeeze(), (_i.T[:,0].squeeze() ,_i.T[:,1].squeeze())), shape=(len(b),len(x)))
            
            #Compute whether or not the constraint objective is violated:
            V = A@x <= b
            print(V)

            assert V.all()

[gasse]

Hi @millskyle ,

A few things could go wrong. One thing for sure which I can see is that you're collecting solution values for variables in the original space

            sol_vals = np.asarray([model.getSolVal(sol, var) for var in model.getVars(transformed=False)])

The NodeBipartite observations operate in the transformed space, so you'll want transformed=True here.

A second thing which could cause trouble is that NodeBipartite extracts information about the current node's LP, and thus it concerns only the variables which are present in the current LP (LP columns, see here). I think in general all transformed variables are in the LP, but sometimes it might not be the case (fixed variables maybe ?). And most importantly, the ordering of the variables in the LP might be different from the ordering of variables in the transformed problem formulation. To figure things out, I invite you to print the variable indexes in the transformed problem, and also their indexes in the LP (if they are indeed part of the LP). You can use the following API:

var.getIndex()
var.isInLP()
var.getCol().getLPPos()

See also related discussions here and here.

Unfortunately the Ecole documentation is not very clear about that.

I hope that helps !

Best,
Maxime

[millskyle]

Hi Maxime, thank you for the reply.

Even if I index it according to the LP, the constraints still seem to be violated according to my calculation. See the two lines demarcated with #!!!!!!!!!!!!.

I'm at a loss as to why this isn't aligning. If constraints (i.e. Rows) are represented in NodeBipartite as Ax <= b, where b is NodeBipartiteObs.RowFeatures.bias and A is the reconstructed matrix from NodeBipartiteObs.edge_features, am I correct in asserting that A@x <= b should be true for any feasible solution x?

Kyle

from scipy import sparse
import numpy as np
import ecole
import ecole as ec
import pyscipopt
from pathlib import Path

import torch
import torch.nn.functional as F

class ObservationFunction():
    def __init__(self, problem):
        # called once for each problem benchmark
        self.problem = problem  # to devise problem-specific observations
        self.nbp = ec.observation.NodeBipartite()

    def seed(self, seed):
        # called before each episode
        # use this seed to make your code deterministic
        pass

    def before_reset(self, model):
        self.nbp.before_reset(model)
        # called when a new episode is about to start
        pass

    def extract(self, model, done):
        return model.as_pyscipopt(), self.nbp.extract(model, done)


if __name__=="__main__":
    instances_path = Path("/home/kmills/git/ml4co-competition/instances/1_item_placement/train/")
    instances = iter(list(instances_path.glob('*.mps.gz')))
    scip_parameters = {'separating/maxrounds': 0, 'presolving/maxrestarts': 0, 'limits/time': 300}

    env = ecole.environment.PrimalSearch(trials_per_node=500,
                                         depth_freq=1,
                                         depth_start=0, #start on the root node
                                         depth_stop=0, #do not allow it to proceed deeper
                                         observation_function=ObservationFunction(problem=None),
                                         scip_params=scip_parameters)

    for seed, instance in enumerate(instances):
        observation, action_set, _, done, _ = env.reset(str(instance))
        while not done:

            model, node_bipartite = observation  

            best_sol = model.getBestSol()
            
            # make two arrays of None
            sol_vals = [None for _ in model.getVars(transformed=True)] 
            nbp_vals = [None for _ in sol_vals]

            for var in model.getVars(transformed=True):

                #get the position in the LP problem; this should
                #be the index in the NodeBipartiteObs
                assert var.isInLP(), "Var must be in LP to proceed." 
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                lp_index = var.getCol().getLPPos()
                sol_vals[lp_index] =model.getSolVal(best_sol, var)
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

                # THis ends up printing something like:
                # Variable t_place_70_1 has index 472 in the LP problem. It has value 1.0 from pyscipopt and value -0.0 in the nbp observation.
                # This isn't necessarily problematic since NBP contains the LP solution.
                print(f"Variable {var} has index {lp_index} in the LP problem. It has value {model.getSolVal(best_sol, var)} from pyscipopt and value {node_bipartite.column_features[lp_index, int(node_bipartite.ColumnFeatures.solution_value)]} in the nbp observation.")


            x = np.asarray(sol_vals)

            b = node_bipartite.row_features[:, int(node_bipartite.RowFeatures.bias)]
            _A = node_bipartite.edge_features.values             
            _i = node_bipartite.edge_features.indices            

            #Construct the constraint matrix as a sparse matrix: 
            A = sparse.csr_matrix((_A.squeeze(), (_i.T[:,0].squeeze() ,_i.T[:,1].squeeze())), shape=(len(b),len(x)))
            
            #Compute whether or not the constraint objective is violated:
            V = A@x <= b
            print(V)

            assert V.all()

[KyleM4t1qbit]

Just following up if you have any more ideas @gasse. Thank you!

[gasse]

Hi @KyleM4t1qbit . I am a bit under water at the moment with a conference deadline, but should become clear again in a couple of days. I'll have a look at this issue ASAP then.

Best,
Maxime

[gasse]

Hi @millskyle ,

So I've looked into it and it appears the issue comes from Ecole. In particular, in Ecole's NodeBipartite observation function we had decided at some point to normalize each row by its L2 norm (computed only on the coefficients), as the observations's purpose was to be used by ML algorithms, which tend to like normalized inputs better. This way, a constraint such as 10x5 <= 3 becomes 1x5 <= 0.3. However, it looks like this was implemented only for the row right-hand-sides (b in your example), but not for the row coefficients (matrix A in your example). As a result, the constraints returned by NodeBipartite are not valid any more in the Ax <= b form. I'll open an issue in Ecole, meanwhile here is a quick fix:

from scipy import sparse
import scipy.sparse.linalg
import numpy as np
import ecole
import ecole as ec
import pyscipopt
from pathlib import Path


class ObservationFunction():
    def __init__(self, problem):
        # called once for each problem benchmark
        self.problem = problem  # to devise problem-specific observations
        self.nbp = ec.observation.NodeBipartite()

    def seed(self, seed):
        # called before each episode
        # use this seed to make your code deterministic
        pass

    def before_reset(self, model):
        self.nbp.before_reset(model)
        # called when a new episode is about to start
        pass

    def extract(self, model, done):
        return model.as_pyscipopt(), self.nbp.extract(model, done)


if __name__=="__main__":
    instances_path = Path("instances/1_item_placement/train/")
    instances = iter(list(instances_path.glob('*.mps.gz')))
    scip_parameters = {'separating/maxrounds': 0, 'presolving/maxrestarts': 0, 'limits/time': 300}

    env = ecole.environment.PrimalSearch(trials_per_node=500,
                                         depth_freq=1,
                                         depth_start=0, #start on the root node
                                         depth_stop=0, #do not allow it to proceed deeper
                                         observation_function=ObservationFunction(problem=None),
                                         scip_params=scip_parameters)

    for seed, instance in enumerate(instances):
        observation, action_set, _, done, _ = env.reset(str(instance))
        while not done:

            model, node_bipartite = observation  

            best_sol = model.getBestSol()

            # get variables in transformed problem
            trans_vars = model.getVars(transformed=True)

            # get variables in LP (ordered by LP position)
            lp_vars = [var for var in trans_vars if var.isInLP()]
            lp_vars = sorted(lp_vars, key=lambda var: var.getCol().getLPPos())

            # get LP variable values in the best MILP solution so far 
            lp_bestsolval = [model.getSolVal(best_sol, var) for var in lp_vars]
            print("Best solution values")
            print(lp_bestsolval)

            # get LP variable values in the current node's LP relaxation
            lp_lpsolval = [model.getSolVal(None, var) for var in lp_vars]
            print("LP solution values")
            print(lp_lpsolval)

            b = node_bipartite.row_features[:, int(node_bipartite.RowFeatures.bias)]
            x = np.asarray(lp_lpsolval)
            _A = node_bipartite.edge_features.values             
            _i = node_bipartite.edge_features.indices            
            A = sparse.csr_matrix((_A, (_i[0] ,_i[1])), shape=(len(b), len(x)))

            # bugfix
            row_norms = sparse.linalg.norm(A, axis=1)
            row_norms[row_norms == 0] = 1
            for i in range(A.shape[0]):
                A[i] /= row_norms[i] if row_norms[i] > 0 else 1

            for i, (aix, bi) in enumerate(zip(A@x, b)):
                print(f"c{i:03d}: {aix} <= {bi} is {aix <= bi}")

            assert all(A@x - b <= 1e-5)

[dchetelat]

Hi, just a little update. Maxime's suggestion to row-normalize the A coefficients in Ecole have been implemented and will be effective in the next Ecole release (pull request).