Example: callbacks_first

Demonstrates the use of generic callbacks

#!/usr/bin/env python
# -*- coding: utf-8 -*-
from amplpy import AMPL
from pprint import pprint

import amplpy_gurobi as ampls
SOLVER = "gurobi"


# Example description
# Demonstrates the use of generic callbacks

MUTE = True # Set to false to enable more chatter
PRINT_CANDO = False

class Callback(ampls.GenericCallback):
    '''Callback class that implements user defined functionalities'''
    def __init__(self):
        self.calls = {w: 0 for w in ampls.Where}
        self.not_mapped = {}
        self.last_bound = 10e10

    def handle_not_mapped(self):
        '''Function to record names and number of occurences of callback
           calls not mapped by ampls' GenericCallback'''
        where_name = self.get_where_string()
        if where_name in self.not_mapped:
            self.not_mapped[where_name]+=1
        else:
            self.not_mapped[where_name]=1
    
    def run(self):
        '''Function executed at each callback call'''
        # Get the generic stage where the callback has been called
        t = self.get_ampl_where()
        # Record it
        self.calls[t] += 1

        if t == ampls.Where.MSG and not MUTE:
            print(self.get_message())
            return 0

        if t == ampls.Where.LPSOLVE and not MUTE:
            print(f"LP solve, {self.get_value(ampls.Value.ITERATIONS)} iterations")
            return 0

        if t == ampls.Where.PRESOLVE and not MUTE:
            print("Presolve, eliminated {} rows and {} columns.".format(
                        self.get_value(ampls.Value.PRE_DELROWS),
                        self.get_value(ampls.Value.PRE_DELCOLS)))
            return 0

        if t in [ampls.Where.MIPNODE, ampls.Where.MIPSOL]:
            obj_bound = self.get_value(ampls.Value.MIP_OBJBOUND)
            if self.last_bound == obj_bound: return 0 # Don't do anything if bound hasn't improved

            self.last_bound = obj_bound
            nonzeroes = sum(x != 0 for x in self.get_solution_vector())
            gap = self.get_value(ampls.Value.MIP_RELATIVEGAP)
            print(f"{t.name}: nonzeroes={nonzeroes} obj={self.get_obj()} bound={obj_bound} gap={gap}")

        if t == ampls.Where.MIP: 
            if not PRINT_CANDO:
                return 0
            # Certain solvers have stop at points during execution of a MIP
            # solve that are neither a new node or a new solution. Check for functionalities
            if self.can_do(ampls.CanDo.GET_MIP_SOLUTION):
                print("In MIP where I can get MIP solution")
            if self.can_do(ampls.CanDo.GET_LP_SOLUTION):
                print("In MIP where I can get LP solution")
            if self.can_do(ampls.CanDo.IMPORT_SOLUTION):
                print("In MIP where I can import an heuristic solution")

        if t == ampls.Where.NOTMAPPED:
            # In case of not mapped, we record the solver-specific stage
            self.handle_not_mapped()
        return 0


# Load a model with amplpy
ampl = AMPL()
ampl.read("models/queens.mod")
ampl.param["size"]=10

# Export to ampls
model = ampl.to_ampls(SOLVER)

# In case of CPLEX, multithreading when using callbacks must
# be handled by the user. See cplex-genericbranch
if SOLVER == "cplex": model.set_option("threads", 1)

# Create a callback and set it in ampls
cb = Callback()
model.set_callback(cb)

# Start the optimization
model.optimize()

# Import the solution back to amplpy
ampl.import_solution(model)
ampls_obj = model.get_obj()
ampl_obj = ampl.get_current_objective().value()
assert(ampls_obj == ampl_obj)

# Show the number of calls to the callack
pprint(cb.calls)
print(cb.not_mapped)

Example: callbacks_stopping

Stopping the solution process given a certain criteria using callbacks

#!/usr/bin/env python
# -*- coding: utf-8 -*-
from math import e
from amplpy import AMPL
import time

import amplpy_gurobi as ampls
SOLVER = "gurobi"

# Example description
# This example shows how to stop the solution process of a MIP 
# problem with different criteria: one based on MIP gap, 
# the other based on the time since an improvement in the solution
# has been achieved.

# Set to true to activate the MIP gap criteria
STOP_AT_MIPGAP = False
# Set to true to stop when no solution improvement has been detected
# for a certain time
STOP_AT_NOSOLUTIONIMPROVEMENT = True

## Desired MIP gap where to stop execution
DESIREDGAP = 0.01
# Define the maximum time without improvement in 
# the solution (in seconds)
MAX_TIME = 10  
 
MUTE = True


class StoppingCallback(ampls.GenericCallback):
    '''User defined class implementing the stopping functionality.
       Returning -1 from a callback makes the solver interrupt the
       solution process'''
    def __init__(self):
        self.last_bound = 10e10
        self.last_obj = None
        self.last_obj_time = None
        
    def _evaluate_mipgap(self, mipsol : bool):
        # At each solution and each node, we check if the desired gap 
        # abs(bestobjective-objbound)/bestobjective
        # has been achieved. If so, terminate the solution process.
        obj_bound = self.get_value(ampls.Value.MIP_OBJBOUND)
        gap = self.get_value(ampls.Value.MIP_RELATIVEGAP)
        if not mipsol: # print only if bound has improved
            if self.last_bound == obj_bound: return 0
            self.last_bound = obj_bound;
        nnz =  sum(x != 0 for x in self.get_solution_vector())
        print(f"Current objective: {self.getObj()} - bound: {obj_bound} - relgap: {gap} - nonzeroes: {nnz}")
        if gap < DESIREDGAP:
            print("Desired gap reached, terminating")
            return -1
        return 0
    
    def _evaluate_solutionprogress(self):
        # Each time we see a new solultion, check its value;
        # If no improvement in the solution value has been detected
        # for a specified amount of time, interrupt the solution process.
        obj = self.get_obj() 
        print(f"Objective value {obj}")
        if obj != self.last_obj:
            if self.last_obj_time is not None:
                print(f"Improvement in the solution after {time.time() - self.last_obj_time} seconds.")
            self.last_obj_time=time.time()
            self.last_obj = obj
        else:
            if time.time() - self.last_obj_time > MAX_TIME:
                print(f"No improvement in the solution for {time.time() - self.last_obj_time} seconds, terminating.")
                return -1
        return 0

    def run(self):
        returnvalue=0
        t = self.get_ampl_where()
        if t == ampls.Where.MIPNODE:
            if STOP_AT_MIPGAP:
                if not MUTE: print(f"New MIP node, count {self.get_value(ampls.Value.MIP_NODES)}")
                return self._evaluate_mipgap(False)
        if t == ampls.Where.MIPSOL:
            print(f"New MIP solution.")
            if STOP_AT_MIPGAP:
                if self._evaluate_mipgap(True)==-1:
                    return -1
            if STOP_AT_NOSOLUTIONIMPROVEMENT:
                return self._evaluate_solutionprogress()
        return 0


# Load a model with amplpy
ampl = AMPL()
ampl.read("models/queens.mod")

# Export to ampls
model = ampl.to_ampls(SOLVER);
# Specify option to return the mip gap to AMPL after solving,
# see https://mp.ampl.com/features-guide.html, useful if
# the MIP gap value is then needed in amplpy
model.setOption("mip:return_gap", 7)

# Declare a callback
cb = StoppingCallback()

# In case of CPLEX, multithreading when using callbacks must
# be handled by the user. See cplex-genericbranch
if SOLVER == "cplex": model.set_option("threads", 1)

# Set the callback and start the optimization
model.set_callback(cb)
model.optimize()

# When the optimization finishes, show the results
obj = model.get_obj()
if abs(158 - obj) >= DESIREDGAP:
    print("Error")

# Import the solution back to amplpy
ampl.import_solution(model)
ampl_obj = ampl.get_current_objective().value()
if STOP_AT_MIPGAP:
    if abs(158 - ampl_obj) >= DESIREDGAP:
        print("Error")

# Print MIP gap
ampl_rel_gap = ampl.get_value("Initial.relmipgap")
print(f"From AMPL we see a relative gap of {ampl_rel_gap}")

Example: load_model

Shows how to export a model from amplpy to ampls, how to solve it and how to get the solution as a vector or as a dictionary

#!/usr/bin/env python
# -*- coding: utf-8 -*-
from amplpy import AMPL

import amplpy_gurobi as ampls 
SOLVER = "gurobi"

# Example description
# Shows how to export a model from amplpy to ampls,
# how to solve it and how to get the solution as a vector or
# as a dictionary

def doStuff(a):
    '''Generic function doing the optimization and reading the results'''

    # Optimize with default settings
    a.optimize()
    print("Model status:", a.get_status())
    # Print the objective function
    print("Objective:", a.get_obj())
    
    # Get the solution as vector and count the nonzeroes
    sol = a.get_solution_vector()
    countnz = sum(x != 0 for x in sol)

    # Get the solution as dictionary (name : value)
    sol = a.get_solution_dict()
    nonzeroes = {name : value for name, value in sol.items() if value != 0}
    for (name, value) in nonzeroes.items():
        print("{} = {}".format(name, value))

    print(f"Non zeroes: vector = {countnz}, dict = {len(nonzeroes)}")

# Using amplpy
ampl = AMPL()
ampl.read("models/queens.mod")
ampl.param["size"]=10
# Export to specified solver
ampls_model = ampl.to_ampls(SOLVER)
# Call generic function
doStuff(ampls_model)

Example: multiple_models

Shows how to create multiple models with AMPL, export an AMPLs instance for a specific solver, set some solver options and get status/solution, and get the results back in AMPL. It also shows how to use the native C API of various solvers

#!/usr/bin/env python
# -*- coding: utf-8 -*-
import amplpy
import pandas as pd
from typing import Any, Dict

import amplpy_cplex as ampls
SOLVER = "cplex"

# Example description
# Shows how to create multiple models with AMPL, export an AMPLs instance
# for a specific solver, set some solver options and get status/solution, and get the results
# back in AMPL. It also shows how to use the native C API of various solvers


def makeAmplModel(num_vars :int=10, flip_objective:bool = False, infeasible:bool = False, 
                  presolve_level : bool = 10):
    """ Create an instance of AMPL and a model"""
    ampl = amplpy.AMPL()
    # Note that if a model is detected as infeasible by the AMPL presolver,
    # it will not be exported to the solver, and an exception will be thrown.
    ampl.set_option("presolve", presolve_level)
    
    ampl.eval("""
    set varIndexes within Integers;
    set constraintIndexes within Integers;

    param varCosts { varIndexes };
    param rightHandSides { constraintIndexes };

    var x {i in varIndexes} >=0;

    maximize objectiveFunction:
        sum { i in varIndexes } varCosts[i] * x[i];

    subject to mainConstraints { j in constraintIndexes } :
        x[j] + x[j+1] <= rightHandSides[j];"""
    )

    varIndexDf = pd.DataFrame(
        data = [(i,i if not flip_objective else -i) for i in range(num_vars)],
        columns = ["varIndexes","varCosts"]
    ).set_index("varIndexes")
    ampl.setData(varIndexDf,set_name = "varIndexes")
    
    constraintIndexesDf = pd.DataFrame(
        data = [(i,1 if not infeasible else -1) for i in range(num_vars-1)],
        columns = ["constraintIndexes","rightHandSides"]
    ).set_index("constraintIndexes")
    ampl.setData(constraintIndexesDf,set_name = "constraintIndexes")
    return ampl

def solveModel(ampl:amplpy.AMPL, options:Dict[str,Any] = {}, useNativeCall: bool=True):
    """Solves a model by exporting it to ampls and using either the generic optimize()
       function or the solver native C API"""
    model= ampl.to_ampls(SOLVER)
    for o, v in options.items():
        model.set_option(o,v) # Set options from
       # ampl driver (see for example https://dev.ampl.com/solvers/gurobi/options.html,  
       #  https://dev.ampl.com/solvers/xpress/options.html# )
       # Most of them will work the same way across solvers
    
    # For specific solvers there are also wrappers for the 
    # gurobi parameter accessors. See https://www.gurobi.com/documentation/10.0/refman/c_parameter_examples.html
    if SOLVER == "gurobi":
        model.set_param("Threads", 4) 

    # A third option is to use the solver C API directy, by getting access to the underlying 
    # gurobi model/environment pointer. This is suitabe to use functionalities provided
    # by the specific solver C API
    # The two statmeents below are therefore equivalent, the first uses the Gurobi C API directly,
    # the latter uses the (generic) function provided by ampls
    if useNativeCall:
        if SOLVER == "gurobi":
            grbnative = model.get_grb_model()
            ampls.GRBoptimize(grbnative)
        if SOLVER == "cplex":
            cplexnativeenv = model.get_cplex_env()
            cplexnativemodel = model.get_cplex_lp()
            ampls.CPXlpopt(cplexnativeenv, cplexnativemodel)
        if SOLVER == "xpress":
            xpressnative = model.getXPRSprob()
            ampls.XPRSlpoptimize(xpressnative, None)
    else:
        model.optimize()
    ampl.import_solution(model)
    return model



def createAndSolveSimpleModel(useNativeCall=True):
    num_vars = 10
    ampl = makeAmplModel(num_vars=num_vars)
    model = solveModel(ampl, useNativeCall=useNativeCall)
    
    print(f"Status is {model.get_status().name}")
    assert model.get_status() == ampls.Status.OPTIMAL

    # Num Vars is even -> last var idx is odd -> expect odd idx vars to be 1
    # Num Vars is odd -> last var idx is even -> expect even idx vars to be 1
    expectedSolution = [abs(i%2-num_vars %2) for i in range(num_vars )]    
    expectedObjective = sum( [i*expectedSolution[i] for i in range(num_vars )])    
    solverObjective = model.getObj()
    assert solverObjective == expectedObjective
    # Check that AMPL has received the correct objective
    assert solverObjective == ampl.getCurrentObjective().value()
    print("Completed Simple Model Test.")

def createAndSolveInfeasibleModel(presolve_level : bool = 10, useNativeCall=True):
    ampl = makeAmplModel(infeasible=True, presolve_level=presolve_level)
    # Set some solver options
    # Setting iisfind to 1 to find the Irreducible Infeasibility Subset
    model = solveModel(ampl, {"outlev" : 1, "iisfind" : 1}, useNativeCall=useNativeCall)
    print(f"Status is {model.get_status().name}")
    assert model.get_status() == ampls.Status.INFEASIBLE
    
    # Display IIS for infeasible model by converting it to a pandas dataframe
    print(ampl.get_data("_varname, _var.iis").to_pandas().set_index('_varname'))
    print(ampl.get_data("_conname, _con.iis").to_pandas().set_index('_conname'))
    print("Completed Infeasible Model Test.")


if __name__ == "__main__":
    try:
        # With presolve set to a high level, AMPL detects the infeasibilitty
        # and does not export the model. Throws an explainatory runtime error
        createAndSolveInfeasibleModel(10)
    except ampls.AMPLSolverException as e:
        print(e)
    # Turning off presolve makes ampl actually export the model
    # In the fucntion we'll set some options to find the source of the infeasibility
    createAndSolveInfeasibleModel(0)

    # Repeated solves below for feasible model:
    createAndSolveSimpleModel()
    createAndSolveSimpleModel()
    

Example: options

How to set options in the solver driver using: 1) solver driver options 2) solver specific parameters

#!/usr/bin/env python
# -*- coding: utf-8 -*-
from amplpy import AMPL

import amplpy_cplex as ampls
SOLVER = "cplex"

# Example description
# How to set options in the solver driver using:
#   1) solver driver options 
#   2) solver specific parameters

def define_model():
    ampl = AMPL()
    ampl.eval('var x binary; var y binary; var z binary;'
              'minimize TotalSum: z + 1;'
              'subj to C1: x+y >= 1;')
    return ampl

def solve(ampl: AMPL):
    mod = ampl.to_ampls(SOLVER, ["sol:stub=test_multi"])

    # Use AMPL driver options
    # See https://mp.ampl.com/features-guide.html
    # and https://dev.ampl.com/solvers/index.html
    mod.set_option("outlev", 1)

    # Use solver specific parameters
    if SOLVER=="gurobi":
        mod.set_param(ampls.GRB_INT_PAR_SOLUTIONLIMIT, 5)
    if SOLVER=="cplex":
        mod.set_param(ampls.CPXPARAM_MIP_Pool_Capacity, 5)
    if SOLVER=="xpress":
        mod.set_param(ampls.XPRS_MSP_SOLUTIONS, 5)

    # Optimize
    mod.optimize()

    # Import into AMPL
    ampl.import_solution(mod)
    print(f"Found {ampl.get_value('Initial.nsol')} solutions")
    print(f"Objective value: {mod.get_obj()}")


if __name__ == "__main__":
    a = define_model()
    solve(a)

Example: pattern_generation

Solves a cutting stock problem using two AMPL problems for pattern generation and shows how to terminate the solution of the final MIP with a callback by increasing the allowed MIP gap depending on how many solutions have been found

#!/usr/bin/env python
# -*- coding: utf-8 -*-
from math import floor
from amplpy import AMPL

import amplpy_gurobi as ampls
SOLVER = "gurobi"


# Example description
# Solves a cutting stock problem using two AMPL problems for pattern generation
# and shows how to terminate the solution of the final MIP with a callback by
# increasing the allowed MIP gap depending on how many solutions have been found

def create_models() -> tuple:
    master="""param nPatterns integer > 0;

    set PATTERNS = 1..nPatterns;  # patterns
    set WIDTHS;                   # finished widths
    param order {WIDTHS} >= 0;    # rolls of width j ordered
    param overrun;                # permitted overrun on any width
    param rawWidth;               # width of raw rolls to be cut
    param rolls {WIDTHS,PATTERNS} >= 0, default 0;   
                                  # rolls of width i in pattern j
    var Cut {PATTERNS} integer >= 0;  # raw rolls to cut in each pattern

    minimize TotalRawRolls: sum {p in PATTERNS} Cut[p];
    subject to OrderLimits {w in WIDTHS}:
      order[w] <= sum {p in PATTERNS} rolls[w,p] * Cut[p] <= order[w] + overrun;"""

    subproblem = """set SIZES;          # ordered widths
    param cap >= 0;     # width of raw rolls 
    param val {SIZES};  # shadow prices of ordered widths (from Master)
    var Qty {SIZES} integer >= 0;  # number of each width in generated pattern
    maximize TotVal: sum {s in SIZES} val[s] * Qty[s];
    subject to Cap: sum {s in SIZES} s * Qty[s] <= cap;"""

    # Data
    roll_width = 6000
    overrun = 30
    orders = {
    1630: 172, 1625: 714,  1620: 110, 1617: 262, 1540: 32, 1529: 100, 1528: 76,
    1505: 110, 1504: 20, 1484: 58, 1466: 15, 1450: 10, 1283: 40, 1017: 50,
    970: 70, 930: 8, 916: 210, 898: 395, 894: 49, 881: 17, 855: 20, 844: 10,
    805: 718, 787: 17, 786: 710, 780: 150, 754: 34, 746: 15, 707: 122, 698: 7,
    651: 10, 644: 15, 638: 10, 605: 10, 477: 4, 473: 34, 471: 25, 468: 10,
    460: 908, 458: 161, 453: 765, 447: 21, 441: 20, 422: 318, 421: 22,
    419: 382, 396: 22,  309: 123,266: 35
    }
    widths = list(sorted(orders.keys(), reverse=True))

    # Create master
    Master = AMPL()
    Master.eval(master)
    # Send scalar values
    Master.param["nPatterns"] = len(widths)
    Master.param["overrun"] = overrun
    Master.param["rawWidth"] = roll_width
    # Send order vector
    Master.set["WIDTHS"] = widths
    Master.param["order"] = orders
    # Generate and send initial pattern matrix
    Master.param["rolls"] = {
        (widths[i], 1+i): int(floor(roll_width/widths[i]))
        for i in range(len(widths))
    }
    # Define a param for sending new patterns to master problem
    Master.eval("param newPat {WIDTHS} integer >= 0;")
    # Set solve options
    Master.option["solver"] = SOLVER
    Master.option["relax_integrality"] =  1

     
    # Create subproblem
    Sub = AMPL()
    Sub.eval(subproblem)

    Sub.set["SIZES"] = widths
    Sub.param["cap"] = roll_width
    Sub.option["solver"] = SOLVER
    return Master,Sub


Master,Sub = create_models()
# Main cycle, solves the relaxed master and use the duals
# to populate the subproblem.
# Solve the subproblem to generate new patterns and add them
# to the master
while True:
    Master.solve()

    Sub.param["val"].setValues(Master.con["OrderLimits"].getValues())
    Sub.solve()
    if Sub.obj["TotVal"].value() <= 1.00001:
        break

    Master.param["newPat"].setValues(Sub.var["Qty"].getValues())
    Master.eval("let nPatterns := nPatterns + 1;")
    Master.eval("let {w in WIDTHS} rolls[w, nPatterns] := newPat[w];")



# At this point, we export the non-relaxed master problem
# and solve with ampls
Master.option["relax_integrality"] = 0
# Export model to ampls
# If we plan to import the results back to AMPL, we have to explicitly set what additional
# suffixes we want returned at this stage
# In this case, we want to return the mip gap as a suffix
ampls_model = Master.to_ampls(SOLVER, ["return_mipgap=5"])

class MyCallback(ampls.GenericCallback):
    """Callback implementing the stopping rule"""
    def __init__(self, stoprule):
      super(MyCallback, self).__init__()
      self._stoprule = stoprule
      self._current = 0
      self._continueOpt = True

    def setCurrentGap(self):
      print("Increasing gap tolerance to %.2f%%" % \
                    (100*self._stoprule["gaptol"][self._current]))
      ampls_model.setAMPLParameter(ampls.SolverParams.DBL_MIPGap,
                             self._stoprule["gaptol"][self._current])
      self._current += 1
    def run(self):
        where = self.getAMPLWhere()
        if where == ampls.Where.MIPNODE:
            runtime = self.getValue(ampls.Value.RUNTIME)
            if runtime >= self._stoprule["time"][self._current]:
                print(f"Current is: {self._stoprule['time'][self._current]}")
                print(f"Stopping optimization at {runtime} seconds")
                self._continueOpt = True
                return -1
        return 0



# Callback"s stopping rule is created here...
stopdict = { "time"   : (  1,    2,   3, 60 ),
             "gaptol" : ( .0002, .002, .02, .1 )
}
# ...and initialized in the constructor
callback = MyCallback(stopdict)
ampls_model.setCallback(callback)

# Invoke solver
# Most solvers (e.g. Gurobi https://support.gurobi.com/hc/en-us/articles/360039233191-How-do-I-change-parameters-in-a-callback-)
# do not support changing optimization parameters from a callback
# Instead we have to stop the optimization, change the desired parameters
# and start it again
callback._continueOpt = True
while callback._continueOpt:
  callback._continueOpt = False
  ampls_model.optimize()
  if callback._continueOpt:
    callback.setCurrentGap()

# Import solution from the solver
Master.importSolution(ampls_model)

print(f"From AMPL MIPGap={100*Master.getValue('TotalRawRolls.relmipgap'):.3f}%")
print(f"Objective value: {Master.getValue('TotalRawRolls')}")

Example: tsp_callback

This example uses generic callbacks to solve a travelling salesperson problem with MTZ cuts or with sub-tour elimination constraints.

#!/usr/bin/env python
# -*- coding: utf-8 -*-
from tsp_helpers import tsp_model, ford_fulkerson, UnionFind

import amplpy_gurobi as ampls
SOLVER = "gurobi"

# Example description
# This example uses generic callbacks to solve a travelling salesperson problem 
# with MTZ cuts or with sub-tour elimination constraints.

VERBOSE = True
ENABLE_MTZ = False


if ENABLE_MTZ:
    ENABLE_CB_MIPSOL = False  # not necessary if MTZ cuts are in place
else:
    ENABLE_CB_MIPSOL = True  # must add sub-tour elimination constraints

ENABLE_CB_MIPNODE = not ENABLE_CB_MIPSOL

class MyCallback(ampls.GenericCallback):
    CALL_COUNT_MIPSOL = 0
    CALL_COUNT_MIPNODE = 0

    def mipsol(self):
        self.CALL_COUNT_MIPSOL += 1
        sol = self.get_solution_vector()
        nv = sum(abs(x) > 1e-5 for x in sol)
        if VERBOSE:
            print("MIPSOL #{}, nnz={}".format(self.CALL_COUNT_MIPSOL, nv))

        values = {
            xvars[i]: sol[i]
            for i in xvars
            if abs(sol[i]) > 1e-5
        }
        uf = UnionFind()
        for (u, v) in values:
            uf.link(u, v)
        groups = uf.groups()
        if len(groups) == 1:
            print('Valid solution!')
        else:
            for grp in groups:
                print('> sub-tour: ', grp)
                cutvarnames = [ampls.tuple2var('x', i, j)
                               for i in grp for j in grp if i != j]
                coeffs = [1 for i in range(len(cutvarnames))]
                self.add_lazy(cutvarnames, coeffs,
                             ampls.CutDirection.LE, len(grp)-1)
        return 0

    def mipnode(self):
        self.CALL_COUNT_MIPNODE += 1
        if VERBOSE:
            print("MIPNODE #{}!".format(self.CALL_COUNT_MIPNODE))
        if self.CALL_COUNT_MIPNODE >= 1000:
            return 1
        sol = self.getSolutionVector()
        capacities = {
            xvars[i]: sol[i]
            for i in xvars
            if abs(sol[i]) > 1e-5
        }
        for i in range(1, len(vertices)):
            max_flow, partition = ford_fulkerson(
                capacities, vertices[0], vertices[i]
            )
            if max_flow <= 1-1e-5:
                p1 = set(partition)
                p2 = set(vertices) - p1
                min_cut = sum(
                    capacities.get((i, j), 0)
                    for i in p1 for j in p2
                )
                cutvarnames = [ampls.tuple2var('x', i, j) for i in p1 for j in p2]
                coeffs = [1 for i in range(len(cutvarnames))]
                self.add_cut(cutvarnames, coeffs, ampls.CutDirection.GE, 1)
                print('> max-flow: {}, min-cut: {}, must be == 1'.format(
                    max_flow, min_cut))
                return 0
        return 0

    def run(self):
        try:
            if ENABLE_CB_MIPSOL and self.getAMPLWhere() == ampls.Where.MIPSOL:
                return self.mipsol()
            elif ENABLE_CB_MIPNODE and self.getAMPLWhere() == ampls.Where.MIPNODE:
                return self.mipnode()
            else:
                return 0
        except Exception as e:
            print('Error:', e)
            return 1

# Create the ampls model instance
ampl = tsp_model('tsp/tsp_51_1.txt', ENABLE_MTZ)
m = ampl.to_ampls(SOLVER)
print("Model loaded, nvars=", m.getNumVars())

if ENABLE_CB_MIPSOL:  # needs lazy constraints
    m.enableLazyConstraints()

# Initialize global structures
var_map = dict(m.getVarMapFiltered('x'))
xvars = {
    index: ampls.var2tuple(var)[1:]
    for var, index in var_map.items()}
vertices = list(sorted(set(
    [x[0] for x in xvars.values()] + [x[1] for x in xvars.values()])))

# Set the callback
cb = MyCallback()
m.setCallback(cb)

# Start the optimzation
m.optimize()

# Get the results
obj = m.getObj()
nvars = m.getNumVars()
print("Solved for {} variables, objective {}".format(nvars, obj))
if m.getStatus() == ampls.Status.OPTIMAL:
    ampl.importSolution(m)
    ampl.display('total')

previous

Architecture

next

Reference