diff --git a/Makefile b/Makefile index bdb3c4bed9..84a6ebe253 100644 --- a/Makefile +++ b/Makefile @@ -3,6 +3,16 @@ PYTHON_VER=2.6 GFLAGS_DIR=../gflags-1.4 SWIG_BINARY=swig PROTOBUF_DIR=../protobuf-2.3.0 +# This is the root directory of the CLP installation. Please undefine if CLP is +# not installed. If you have installed CBC, CLP_DIR can have the same value as +# CBC_DIR. +#CLP_DIR=../cbc-2.6.2 +# This is the root directory of the CBC installation. Please undefine if CBC is +# not installed. +#CBC_DIR=../cbc-2.6.2 +# This is the root directory of glpk installation. Please undefine if GLPK is +# not installed. +#GLPK_DIR=../glpk-4.45 # ----- You should not need to modify the following, unless the ----- # ----- configuration is not standard ------ @@ -14,6 +24,25 @@ GFLAGS_INC = -I$(GFLAGS_DIR)/include # This is needed to find protocol buffers. PROTOBUF_INC = -I$(PROTOBUF_DIR)/include +# Define CLP_DIR if unset and if CBC_DIR is set. +ifdef CBC_DIR +ifndef CLP_DIR +CLP_DIR=$(CBC_DIR) +endif +endif +# This is needed to find Coin LP include files. +ifdef CLP_DIR +CLP_INC = -I$(CLP_DIR)/include -DUSE_CLP +endif +# This is needed to find Coin Branch and Cut include files. +ifdef CBC_DIR +CBC_INC = -I$(CBC_DIR)/include -DUSE_CBC +endif +# This is needed to find GLPK include files. +ifdef GLPK_DIR +GLPK_INC = -I$(GLPK_DIR)/include -DUSE_GLPK +endif + # Compilation flags DEBUG=-O3 -DNDEBUG JNIDEBUG=-O1 -DNDEBUG @@ -31,6 +60,15 @@ GFLAGS_LNK = -Wl,-rpath $(GFLAGS_DIR)/lib -L$(GFLAGS_DIR)/lib -lgflags ZLIB_LNK = -lz # This is needed to find libprotobuf.a PROTOBUF_LNK = -Wl,-rpath $(PROTOBUF_DIR)/lib -L$(PROTOBUF_DIR)/lib -lprotobuf -lpthread +ifdef GLPK_DIR +GLPK_LNK = -Wl,-rpath $(GLPK_DIR)/lib -L$(GLPK_DIR)/lib -lglpk +endif +ifdef CLP_DIR +CLP_LNK = -Wl,-rpath $(CLP_DIR)/lib/coin -L$(CLP_DIR)/lib/coin -lClp -lCoinUtils +endif +ifdef CBC_DIR +CBC_LNK = -Wl,-rpath $(CBC_DIR)/lib/coin -L$(CBC_DIR)/lib/coin -lCbcSolver -lCbc -lCgl -lOsi -lOsiCbc -lOsiClp -lOsiVol -lVol +endif # Detect 32 bit or 64 bit OS and define ARCH flags correctly. LBITS := $(shell getconf LONG_BIT) ifeq ($(LBITS),64) @@ -59,24 +97,37 @@ JAVAC_BIN=javac JAVA_BIN=java JNILIBEXT=jnilib FIX_SWIG= + +ifdef GLPK_DIR +GLPK_LNK = -L$(GLPK_DIR)/lib -lglpk +endif +ifdef CLP_DIR +CLP_LNK = -L$(CLP_DIR)/lib/coin -lClp -lOs +endif +ifdef CBC_DIR +CBC_LNK = -L$(CBC_DIR)/lib/coin -lCbcSolver -lCbc -lCgl -lOsi -lOsiCbc -lOsiClp +endif endif CFLAGS= $(SYSCFLAGS) $(DEBUG) -I. $(GFLAGS_INC) $(ARCH) \ - -Wno-deprecated $(PROTOBUF_INC) + -Wno-deprecated $(PROTOBUF_INC) $(CBC_INC) $(CLP_INC) $(GLPK_INC) JNIFLAGS= $(SYSCFLAGS) $(JNIDEBUG) -I. $(GFLAGS_INC) $(ARCH) \ -Wno-deprecated $(PROTOBUF_INC) $(CBC_INC) $(CLP_INC) $(GLPK_INC) LDFLAGS=$(GFLAGS_LNK) $(ZLIB_LNK) $(PROTOBUF_LNK) $(SYS_LNK) +LDLPDEPS= $(GLPK_LNK) $(CBC_LNK) $(CLP_LNK) # Real targets help: @echo Please define target: @echo " - constraint programming: cplibs cpexe pycp javacp" + @echo " - mathematical programming: lplibs lpexe pylp javalp" @echo " - algorithms: algorithmslibs pyalgorithms javaalgorithms" @echo " - graph: graphlibs pygraph javagraph" @echo " - misc: clean" -all: cplibs cpexe pycp javacp algorithmslibs pyalgorithms javaalgorithms graphlibs pygraph javagraph +all: cplibs cpexe pycp javacp algorithmslibs pyalgorithms javaalgorithms graphlibs pygraph javagraph lplibs lpexe pylp javalp + CP_LIBS = \ librouting.a \ libconstraint_solver.a @@ -93,6 +144,11 @@ BASE_LIBS = \ cplibs: $(CP_LIBS) $(BASE_LIBS) +LP_LIBS = \ + liblinear_solver.a + +lplibs: $(LP_LIBS) $(BASE_LIBS) + CPBINARIES = \ costas_array \ cryptarithm \ @@ -106,6 +162,12 @@ CPBINARIES = \ cpexe: $(CPBINARIES) +LPBINARIES = \ + integer_solver_example \ + linear_solver_example + +lpexe: $(LPBINARIES) + ALGORITHMS_LIBS = \ libalgorithms.a @@ -115,6 +177,7 @@ clean: rm -f *.a rm -f objs/*.o rm -f $(CPBINARIES) + rm -f $(LPBINARIES) rm -f */*wrap* rm -f */*.pb.* rm -f *.so @@ -232,6 +295,38 @@ objs/utilities.o:constraint_solver/utilities.cc libconstraint_solver.a: $(CONSTRAINT_SOLVER_LIB_OBJS) ar rv libconstraint_solver.a $(CONSTRAINT_SOLVER_LIB_OBJS) +# Linear Solver Library + +LINEAR_SOLVER_LIB_OBJS = \ + objs/cbc_interface.o \ + objs/clp_interface.o \ + objs/glpk_interface.o \ + objs/linear_solver.o \ + objs/linear_solver.pb.o + +objs/cbc_interface.o:linear_solver/cbc_interface.cc linear_solver/linear_solver.pb.h + $(CCC) $(CFLAGS) -c linear_solver/cbc_interface.cc -o objs/cbc_interface.o + +objs/clp_interface.o:linear_solver/clp_interface.cc linear_solver/linear_solver.pb.h + $(CCC) $(CFLAGS) -c linear_solver/clp_interface.cc -o objs/clp_interface.o + +objs/glpk_interface.o:linear_solver/glpk_interface.cc linear_solver/linear_solver.pb.h + $(CCC) $(CFLAGS) -c linear_solver/glpk_interface.cc -o objs/glpk_interface.o + +objs/linear_solver.o:linear_solver/linear_solver.cc linear_solver/linear_solver.pb.h + $(CCC) $(CFLAGS) -c linear_solver/linear_solver.cc -o objs/linear_solver.o + +objs/linear_solver.pb.o:linear_solver/linear_solver.pb.cc + $(CCC) $(CFLAGS) -c linear_solver/linear_solver.pb.cc -o objs/linear_solver.pb.o + +linear_solver/linear_solver.pb.cc:linear_solver/linear_solver.proto + $(PROTOBUF_DIR)/bin/protoc --proto_path=linear_solver --cpp_out=linear_solver linear_solver/linear_solver.proto + +linear_solver/linear_solver.pb.h:linear_solver/linear_solver.pb.cc + +liblinear_solver.a: $(LINEAR_SOLVER_LIB_OBJS) + ar rv liblinear_solver.a $(LINEAR_SOLVER_LIB_OBJS) + # Util library. UTIL_LIB_OBJS=\ @@ -424,6 +519,20 @@ objs/tsp.o: examples/tsp.cc tsp: $(CP_LIBS) $(BASE_LIBS) objs/tsp.o $(CCC) $(CFLAGS) $(LDFLAGS) objs/tsp.o $(CP_LIBS) $(BASE_LIBS) -o tsp +# Linear Programming Examples + +objs/linear_solver_example.o: examples/linear_solver_example.cc + $(CCC) $(CFLAGS) -c examples/linear_solver_example.cc -o objs/linear_solver_example.o + +linear_solver_example: $(LP_LIBS) $(BASE_LIBS) objs/linear_solver_example.o + $(CCC) $(CFLAGS) $(LDFLAGS) objs/linear_solver_example.o $(LP_LIBS) $(BASE_LIBS) $(LDLPDEPS) -o linear_solver_example + +objs/integer_solver_example.o: examples/integer_solver_example.cc + $(CCC) $(CFLAGS) -c examples/integer_solver_example.cc -o objs/integer_solver_example.o + +integer_solver_example: $(LP_LIBS) $(BASE_LIBS) objs/integer_solver_example.o + $(CCC) $(CFLAGS) $(LDFLAGS) objs/integer_solver_example.o $(LP_LIBS) $(BASE_LIBS) $(LDLPDEPS) -o integer_solver_example + # SWIG # pywrapknapsack_solver @@ -484,6 +593,21 @@ objs/routing_wrap.o: constraint_solver/routing_wrap.cc _pywraprouting.so: objs/routing_wrap.o $(CP_LIBS) $(BASE_LIBS) $(LD) -o _pywraprouting.so objs/routing_wrap.o $(CP_LIBS) $(BASE_LIBS) $(LDFLAGS) +# pywraplp + +pylp: _pywraplp.so linear_solver/pywraplp.py $(LP_LIBS) $(BASE_LIBS) + +linear_solver/pywraplp.py: linear_solver/linear_solver.swig linear_solver/linear_solver.h linear_solver/linear_solver.pb.h base/base.swig + $(SWIG_BINARY) $(CLP_INC) $(CBC_INC) $(GLPK_INC) -c++ -python -o linear_solver/linear_solver_wrap.cc -module pywraplp linear_solver/linear_solver.swig + +linear_solver/linear_solver_wrap.cc: linear_solver/pywraplp.py + +objs/linear_solver_wrap.o: linear_solver/linear_solver_wrap.cc + $(CCC) $(CFLAGS) $(PYTHON_INC) -c linear_solver/linear_solver_wrap.cc -o objs/linear_solver_wrap.o + +_pywraplp.so: objs/linear_solver_wrap.o $(LP_LIBS) $(BASE_LIBS) + $(LD) -o _pywraplp.so objs/linear_solver_wrap.o $(LP_LIBS) $(BASE_LIBS) $(LDLPDEPS) $(LDFLAGS) + # ---------- Java Support ---------- # javacp @@ -746,3 +870,36 @@ com/google/ortools/flow/samples/FlowExample.class: javacp com/google/ortools/flo run_FlowExample: compile_FlowExample $(JAVA_BIN) -Djava.library.path=`pwd` -cp .:com.google.ortools.flow.jar com.google.ortools.flow.samples.FlowExample +# javalp + +javalp: com.google.ortools.linearsolver.jar libjnilinearsolver.$(JNILIBEXT) +linear_solver/linear_solver_java_wrap.cc: linear_solver/linear_solver.swig base/base.swig util/data.swig linear_solver/linear_solver.h + $(SWIG_BINARY) $(CLP_INC) $(CBC_INC) $(GLPK_INC) -c++ -java -o linear_solver/linear_solver_java_wrap.cc -package com.google.ortools.linearsolver -outdir com/google/ortools/linearsolver linear_solver/linear_solver.swig + +objs/linear_solver_java_wrap.o: linear_solver/linear_solver_java_wrap.cc + $(CCC) $(JNIFLAGS) $(JAVA_INC) -c linear_solver/linear_solver_java_wrap.cc -o objs/linear_solver_java_wrap.o + +com.google.ortools.linearsolver.jar: linear_solver/linear_solver_java_wrap.cc + $(JAVAC_BIN) com/google/ortools/linearsolver/*.java + jar cf com.google.ortools.linearsolver.jar com/google/ortools/linearsolver/*.class + +libjnilinearsolver.$(JNILIBEXT): objs/linear_solver_java_wrap.o $(LP_LIBS) $(BASE_LIBS) + $(LD) -o libjnilinearsolver.$(JNILIBEXT) objs/linear_solver_java_wrap.o $(LP_LIBS) $(BASE_LIBS) $(LDLPDEPS) $(LDFLAGS) + +# Java Algorithms Examples + +compile_LinearSolverExample: com/google/ortools/linearsolver/samples/LinearSolverExample.class + +com/google/ortools/linearsolver/samples/LinearSolverExample.class: javacp com/google/ortools/linearsolver/samples/LinearSolverExample.java + $(JAVAC_BIN) -cp com.google.ortools.linearsolver.jar com/google/ortools/linearsolver/samples/LinearSolverExample.java + +run_LinearSolverExample: compile_LinearSolverExample + $(JAVA_BIN) -Djava.library.path=`pwd` -cp .:com.google.ortools.linearsolver.jar com.google.ortools.linearsolver.samples.LinearSolverExample + +compile_IntegerSolverExample: com/google/ortools/linearsolver/samples/IntegerSolverExample.class + +com/google/ortools/linearsolver/samples/IntegerSolverExample.class: javacp com/google/ortools/linearsolver/samples/IntegerSolverExample.java + $(JAVAC_BIN) -cp com.google.ortools.linearsolver.jar com/google/ortools/linearsolver/samples/IntegerSolverExample.java + +run_IntegerSolverExample: compile_IntegerSolverExample + $(JAVA_BIN) -Djava.library.path=`pwd` -cp .:com.google.ortools.linearsolver.jar com.google.ortools.linearsolver.samples.IntegerSolverExample diff --git a/Makefile.msv b/Makefile.msv index 9d4ad057d1..f30eaef9b8 100644 --- a/Makefile.msv +++ b/Makefile.msv @@ -4,6 +4,18 @@ ZLIB_DIR=..\\zlib-1.2.5 ZLIB_NAME=zlib.lib SWIG_BINARY=..\\swigwin-2.0.0\\swig.exe PROTOBUF_DIR=..\\protobuf-2.3.0 +# This is the root directory of the CLP installation. Please undefine if CLP is +# not installed. If you have installed CBC, CLP_DIR can have the same value as +# CBC_DIR. +#CLP_DIR=..\\cbc-2.6.2 +# This is the root directory of the CBC installation. Please undefine if CBC is +# not installed. +#CBC_DIR=..\\cbc-2.6.2 +# This is the root directory of glpk installation. Please undefine if GLPK is +# not installed. +#GLPK_VER=4_45 +#GLPK_DIR=..\\glpk-4.45 + # This describes the python installation. PYTHON_INC=/Ic:\\Python27\\include @@ -17,20 +29,49 @@ CCC=cl /EHsc /MD GFLAGS_INC = /I$(GFLAGS_DIR)\\src\\windows /I$(GFLAGS_DIR)\\src ZLIB_INC = /I$(ZLIB_DIR) PROTOBUF_INC = /I$(PROTOBUF_DIR)\\include + +# Define CLP_DIR if unset and if CBC_DIR is set. +!ifdef CBC_DIR +!ifndef CLP_DIR +CLP_DIR=$(CBC_DIR) +!endif +!endif +# This is needed to find Coin LP include files and libraries. +!ifdef CLP_DIR +CLP_INC = /I$(CLP_DIR)\\include /DUSE_CLP +CLP_SWIG = -DUSE_CLP +CLP_LNK = $(CLP_DIR)\\lib\\coin\libClp.lib $(CLP_DIR)\\lib\\coin\libCoinUtils.lib +!endif +# This is needed to find Coin Branch and Cut include files and libraries. +!ifdef CBC_DIR +CBC_INC = /I$(CBC_DIR)\\include /DUSE_CBC +CBC_SWIG = -DUSE_CBC +CBC_LNK = $(CBC_DIR)\\lib\\coin\\libCbcSolver.lib $(CBC_DIR)\\lib\\coin\\libCbc.lib $(CBC_DIR)\\lib\\coin\\libCgl.lib $(CBC_DIR)\\lib\\coin\\libOsi.lib $(CBC_DIR)\\lib\\coin\\libOsiCbc.lib $(CBC_DIR)\\lib\\coin\\libOsiClp.lib +!endif +# This is needed to find GLPK include files and libraries. +!ifdef GLPK_DIR +GLPK_INC = /I$(GLPK_DIR)\\include /DUSE_GLPK +GLPK_SWIG = -DUSE_GLPK +GLPK_LNK = $(GLPK_DIR)\\lib\\glpk_$(GLPK_VER).lib +!endif + CFLAGS= -nologo $(SYSCFLAGS) $(DEBUG) /I. $(GFLAGS_INC) $(ZLIB_INC)\ - $(PROTOBUF_INC) + $(PROTOBUF_INC) $(CBC_INC) $(CLP_INC) $(GLPK_INC) LD = cl /EHsc GFLAGS_LNK = $(GFLAGS_DIR)\\vsprojects\\libgflags\\Release\\libgflags.lib ZLIB_LNK = $(ZLIB_DIR)/$(ZLIB_NAME) PROTOBUF_LNK = $(PROTOBUF_DIR)\\lib\\libprotobuf.lib LDFLAGS=$(GFLAGS_LNK) $(ZLIB_LNK) $(PROTOBUF_LNK) psapi.lib ws2_32.lib +LDLPDEPS=$(CBC_LNK) $(CLP_LNK) $(GLPK_LNK) + # Real targets all: @echo Please define target: @echo " - constraint programming: cplibs, cpexe, pycp" + @echo " - mathematical programming: lplibs, lpexe, pylp" @echo " - algorithms: algoritmlibs, pyalgorithms" @echo " - graph: graphlibs, pygraph" @echo " - misc: clean" @@ -51,6 +92,10 @@ graphlibs: $(GRAPH_LIBS) $(BASE_LIBS) cplibs: $(CPLIBS) $(BASE_LIBS) +LPLIBS = \ + linear_solver.lib + +lplibs: $(LPLIBS) $(BASE_LIBS) CPBINARIES= \ costas_array.exe \ @@ -65,6 +110,12 @@ CPBINARIES= \ cpexe: $(CPBINARIES) +LPBINARIES = \ + integer_solver_example.exe \ + linear_solver_example.exe + +lpexe: $(LPBINARIES) + ALGORITHM_LIBS= \ algorithms.lib @@ -74,11 +125,15 @@ clean: del *.lib del objs\\*.obj del $(CPBINARIES) + del $(LPBINARIES) del *.dll del *.exp del constraint_solver\constraint_solver_wrap.cc del constraint_solver\_pywrapcp.pyd del constraint_solver\pywrapcp.py + del linear_solver\linear_solver_wrap.cc + del linear_solver\_pywraplp.pyd + del linear_solver\pywraplp.py del algorithms\knapsack_solver_wrap.cc del algorithms\_pywrapknapsack_solver.pyd del algorithms\pywrapknapsack_solver.py @@ -197,7 +252,40 @@ objs/utilities.obj:constraint_solver/utilities.cc constraint_solver.lib: $(CONSTRAINT_SOLVER_LIB_OBJS) lib /OUT:constraint_solver.lib $(CONSTRAINT_SOLVER_LIB_OBJS) +# Linear Solver Library + +LINEAR_SOLVER_LIB_OBJS = \ + objs/cbc_interface.obj \ + objs/clp_interface.obj \ + objs/glpk_interface.obj \ + objs/linear_solver.obj \ + objs/linear_solver.pb.obj + +objs/cbc_interface.obj:linear_solver/cbc_interface.cc linear_solver/linear_solver.pb.h + $(CCC) $(CFLAGS) -c linear_solver\\cbc_interface.cc /Foobjs\\cbc_interface.obj + +objs/clp_interface.obj:linear_solver/clp_interface.cc linear_solver/linear_solver.pb.h + $(CCC) $(CFLAGS) -c linear_solver\\clp_interface.cc /Foobjs/clp_interface.obj + +objs/glpk_interface.obj:linear_solver/glpk_interface.cc linear_solver/linear_solver.pb.h + $(CCC) $(CFLAGS) -c linear_solver\\glpk_interface.cc /Foobjs\\glpk_interface.obj + +objs/linear_solver.obj:linear_solver/linear_solver.cc linear_solver/linear_solver.pb.h + $(CCC) $(CFLAGS) -c linear_solver\\linear_solver.cc /Foobjs\\linear_solver.obj + +objs/linear_solver.pb.obj:linear_solver/linear_solver.pb.cc + $(CCC) $(CFLAGS) -c linear_solver\\linear_solver.pb.cc /Foobjs\\linear_solver.pb.obj + +linear_solver/linear_solver.pb.cc:linear_solver/linear_solver.proto + $(PROTOBUF_DIR)\\bin\\protoc --proto_path=linear_solver --cpp_out=linear_solver linear_solver\\linear_solver.proto + +linear_solver/linear_solver.pb.h:linear_solver/linear_solver.pb.cc + +linear_solver.lib: $(LINEAR_SOLVER_LIB_OBJS) + lib /OUT:linear_solver.lib $(LINEAR_SOLVER_LIB_OBJS) + # GFlags copy library. + libgflags.dll: $(GFLAGS_DIR)\vsprojects\libgflags\Release\libgflags.dll copy $(GFLAGS_DIR)\vsprojects\libgflags\Release\libgflags.dll . @@ -374,6 +462,19 @@ objs/tsp.obj: examples/tsp.cc tsp.exe: $(CPLIBS) $(BASE_LIBS) objs/tsp.obj libgflags.dll $(CCC) $(CFLAGS) $(LDFLAGS) objs/tsp.obj $(CPLIBS) $(BASE_LIBS) /Fetsp.exe +# Linear Programming Examples + +objs/integer_solver_example.obj: examples/integer_solver_example.cc + $(CCC) $(CFLAGS) -c examples/integer_solver_example.cc /Foobjs/integer_solver_example.obj + +integer_solver_example.exe: $(LPLIBS) $(BASE_LIBS) objs/integer_solver_example.obj libgflags.dll + $(CCC) $(CFLAGS) $(LDFLAGS) objs/integer_solver_example.obj $(LPLIBS) $(BASE_LIBS) $(LDLPDEPS) /Feinteger_solver_example.exe + +objs/linear_solver_example.obj: examples/linear_solver_example.cc + $(CCC) $(CFLAGS) -c examples/linear_solver_example.cc /Foobjs/linear_solver_example.obj + +linear_solver_example.exe: $(LPLIBS) $(BASE_LIBS) objs/linear_solver_example.obj libgflags.dll + $(CCC) $(CFLAGS) $(LDFLAGS) objs/linear_solver_example.obj $(LPLIBS) $(BASE_LIBS) $(LDLPDEPS) /Felinear_solver_example.exe # SWIG @@ -447,3 +548,23 @@ _pywraprouting.dll: objs/routing_wrap.obj $(CPLIBS) $(BASE_LIBS) libgflags.dll constraint_solver/_pywraprouting.pyd: _pywraprouting.dll copy _pywraprouting.dll constraint_solver\\_pywraprouting.pyd +# pywraplp + +pylp: linear_solver/_pywraplp.pyd linear_solver/pywraplp.py $(LPLIBS) $(BASE_LIBS) + +linear_solver/pywraplp.py: linear_solver/linear_solver.swig linear_solver/linear_solver.h base/base.swig + $(SWIG_BINARY) $(CLP_SWIG) $(CBC_SWIG) $(GLPK_SWIG) -c++ -python -o linear_solver\\linear_solver_wrap.cc -module pywraplp linear_solver\\linear_solver.swig + +linear_solver/linear_solver_wrap.cc: linear_solver/pywraplp.py + +objs/linear_solver_wrap.obj: linear_solver/linear_solver_wrap.cc + $(CCC) $(CFLAGS) $(PYTHON_INC) /I__WIN32__ -c linear_solver\\linear_solver_wrap.cc /Foobjs\\linear_solver_wrap.obj + +_pywraplp.dll: objs/linear_solver_wrap.obj $(LPLIBS) $(BASE_LIBS) libgflags.dll + link /DLL /OUT:_pywraplp.dll objs\\linear_solver_wrap.obj $(LPLIBS) $(BASE_LIBS) $(LDLPDEPS) $(LDFLAGS) $(PYTHON_LNK) + +linear_solver/_pywraplp.pyd: _pywraplp.dll + copy _pywraplp.dll linear_solver\\_pywraplp.pyd + + + diff --git a/com/google/ortools/linearsolver/samples/IntegerSolverExample.java b/com/google/ortools/linearsolver/samples/IntegerSolverExample.java new file mode 100644 index 0000000000..59b4c351e3 --- /dev/null +++ b/com/google/ortools/linearsolver/samples/IntegerSolverExample.java @@ -0,0 +1,54 @@ +// Copyright 2010 Google +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +package com.google.ortools.linearsolver.samples; + +import com.google.ortools.linearsolver.MPConstraint; +import com.google.ortools.linearsolver.MPSolver; +import com.google.ortools.linearsolver.MPVariable; + +/** + * Example of linear solver usage. + * + */ + +public class IntegerSolverExample { + + static { + System.loadLibrary("jnilinearsolver"); + } + + + private static void runFirstIntegerExample(int solverType) { + MPSolver solver = new MPSolver("runFirstIntegerExample", solverType); + double infinity = solver.infinity(); + MPVariable x1 = solver.makeIntVar(0.0, infinity, "x1"); + MPVariable x2 = solver.makeIntVar(0.0, infinity, "x2"); + + solver.addObjectiveTerm(x1, 1); + solver.addObjectiveTerm(x2, 2); + + MPConstraint ct = solver.makeConstraint(17, infinity); + ct.addTerm(x1, 3); + ct.addTerm(x2, 2); + + // Check the solution. + solver.solve(); + System.out.println("Objective value = " + solver.objectiveValue()); + } + + public static void main(String[] args) throws Exception { + runFirstIntegerExample(MPSolver.CLP_LINEAR_PROGRAMMING); + runFirstIntegerExample(MPSolver.GLPK_LINEAR_PROGRAMMING); + } +} diff --git a/com/google/ortools/linearsolver/samples/LinearSolverExample.java b/com/google/ortools/linearsolver/samples/LinearSolverExample.java new file mode 100644 index 0000000000..0248fcfd09 --- /dev/null +++ b/com/google/ortools/linearsolver/samples/LinearSolverExample.java @@ -0,0 +1,78 @@ +// Copyright 2010 Google +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +package com.google.ortools.linearsolver.samples; + +import com.google.ortools.linearsolver.MPConstraint; +import com.google.ortools.linearsolver.MPSolver; +import com.google.ortools.linearsolver.MPVariable; + +/** + * Example of linear solver usage. + * + */ + +public class LinearSolverExample { + + static { + System.loadLibrary("jnilinearsolver"); + } + + + private static void runFirstLinearExample(int solverType) { + MPSolver solver = new MPSolver("runFirstLinearExample", solverType); + double infinity = solver.infinity(); + MPVariable x1 = solver.makeNumVar(0.0, infinity, "x1"); + MPVariable x2 = solver.makeNumVar(0.0, infinity, "x2"); + MPVariable x3 = solver.makeNumVar(0.0, infinity, "x3"); + System.out.println("Number of variables = " + solver.numVariables()); + + solver.addObjectiveTerm(x1, 10); + solver.addObjectiveTerm(x2, 6); + solver.addObjectiveTerm(x3, 4); + solver.setMaximization(); + + MPConstraint c0 = solver.makeConstraint(-infinity, 100.0); + c0.addTerm(x1, 1); + c0.addTerm(x2, 1); + c0.addTerm(x3, 1); + MPConstraint c1 = solver.makeConstraint(-infinity, 600.0); + c1.addTerm(x1, 10); + c1.addTerm(x2, 4); + c1.addTerm(x3, 5); + MPConstraint c2 = solver.makeConstraint(-infinity, 300.0); + c2.addTerm(x1, 2); + c2.addTerm(x2, 2); + c2.addTerm(x3, 6); + System.out.println("Number of constraints = " + solver.numConstraints()); + + // The problem has an optimal solution. + solver.solve(); + + System.out.println("Optimal value = " + solver.objectiveValue()); + System.out.println("x1 = " + x1.solutionValue() + + " + reduced cost = " + x1.reducedCost()); + System.out.println("x2 = " + x2.solutionValue() + + " + reduced cost = " + x2.reducedCost()); + System.out.println("x3 = " + x3.solutionValue() + + " + reduced cost = " + x3.reducedCost()); + System.out.println("c0 dual value = " + c0.dualValue()); + System.out.println("c1 dual value = " + c1.dualValue()); + System.out.println("c2 dual value = " + c2.dualValue()); + } + + public static void main(String[] args) throws Exception { + runFirstLinearExample(MPSolver.CLP_LINEAR_PROGRAMMING); + runFirstLinearExample(MPSolver.GLPK_LINEAR_PROGRAMMING); + } +} diff --git a/linear_solver/__init__.py b/linear_solver/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/linear_solver/cbc_interface.cc b/linear_solver/cbc_interface.cc new file mode 100644 index 0000000000..29a4f81634 --- /dev/null +++ b/linear_solver/cbc_interface.cc @@ -0,0 +1,510 @@ +// Copyright 2010 Google +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. +// + +#include "base/commandlineflags.h" +#include "base/integral_types.h" +#include "base/logging.h" +#include "base/scoped_ptr.h" +#include "base/stringprintf.h" +#include "base/timer.h" +#include "base/stl_util-inl.h" + +#include "linear_solver/linear_solver.h" + +#if defined(USE_CBC) + +#undef PACKAGE +#undef VERSION +#include "coin/CbcMessage.hpp" +#include "coin/CbcModel.hpp" +#include "coin/CoinModel.hpp" +#include "coin/OsiClpSolverInterface.hpp" + +#include "coin/CglGomory.hpp" +#include "coin/CglProbing.hpp" +#include "coin/CglKnapsackCover.hpp" +#include "coin/CglOddHole.hpp" +#include "coin/CglClique.hpp" +#include "coin/CglFlowCover.hpp" +#include "coin/CglMixedIntegerRounding.hpp" +#include "coin/config_cbc.h" + +// Heuristics + +DECLARE_double(solver_timeout_in_seconds); +DECLARE_string(solver_write_model); + +// Parameters for cuts. ** OBSOLETE ** +// TODO(user): Add doc on parameters. +DEFINE_bool(cbc_probing_use_objective, true, ""); +DEFINE_int32(cbc_probing_max_pass, 3, ""); +DEFINE_int32(cbc_probing_max_probe, 100, ""); +DEFINE_int32(cbc_probing_max_look, 50, ""); +DEFINE_int32(cbc_probing_row_cuts, 3, ""); + +DEFINE_int32(cbc_gomory_limit, 300, ""); + +DEFINE_double(cbc_odd_hole_minimum_violation, 0.005, ""); +DEFINE_double(cbc_odd_hole_minimum_violation_per, 0.00002, ""); +DEFINE_int32(cbc_odd_hole_maximum_entries, 200, ""); + +DEFINE_bool(cbc_clique_start_clique_report, false, ""); +DEFINE_bool(cbc_clique_row_clique_report, false, ""); + +namespace operations_research { + +class CBCInterface : public MPSolverInterface { + public: + // Constructor that takes a name for the underlying glpk solver. + explicit CBCInterface(MPSolver* const solver); + virtual ~CBCInterface(); + + // ----- Reset ----- + virtual void Reset(); + + // Sets the optimization direction (min/max). + virtual void SetOptimizationDirection(bool maximize); + + // ----- Solve ----- + // Solve the problem using the parameter values specified. + virtual MPSolver::ResultStatus Solve(const MPSolverParameters& param); + + // TODO(user): separate the solve from the model extraction. + virtual void ExtractModel() {} + + // Write model + virtual void WriteModel(const string& filename); + + // SuppressOutput. + virtual void SuppressOutput(); + + // Query problem type. + virtual bool IsContinuous() const { return false; } + virtual bool IsLP() const { return false; } + virtual bool IsMIP() const { return true; } + + // Modify bounds. + virtual void SetVariableBounds(int var_index, double lb, double ub); + virtual void SetVariableInteger(int var_index, bool integer); + virtual void SetConstraintBounds(int row_index, double lb, double ub); + + // Add constraint incrementally. + void AddRowConstraint(MPConstraint* const ct); + // Add variable incrementally. + void AddVariable(MPVariable* const var); + // Change a coefficient in a constraint. + virtual void SetCoefficient(MPConstraint* const constraint, + MPVariable* const variable, + double coefficient) { + sync_status_ = MUST_RELOAD; + } + // Clear a constraint from all its terms. + virtual void ClearConstraint(MPConstraint* const constraint) { + sync_status_ = MUST_RELOAD; + } + + // Change a coefficient in the linear objective. + virtual void SetObjectiveCoefficient(MPVariable* const variable, + double coefficient) { + sync_status_ = MUST_RELOAD; + } + // Change the constant term in the linear objective. + virtual void SetObjectiveOffset(double value) { + sync_status_ = MUST_RELOAD; + } + // Clear the objective from all its terms. + virtual void ClearObjective() { + sync_status_ = MUST_RELOAD; + } + + // Number of simplex iterations + virtual int64 iterations() const; + // Number of branch-and-bound nodes. Only available for discrete problems. + virtual int64 nodes() const; + // Best objective bound. Only available for discrete problems. + virtual double best_objective_bound() const; + + virtual void ExtractNewVariables() {} + virtual void ExtractNewConstraints() {} + virtual void ExtractObjective() {} + + virtual string SolverVersion() const { + return PACKAGE_STRING; + } + private: + // Reset best objective bound to +/- infinity depending on the + // objective sense. + void ResetBestObjectiveBound(); + + // Set all parameters in the underlying solver. + virtual void SetParameters(const MPSolverParameters& param); + // Set each parameter in the underlying solver. + virtual void SetRelativeMipGap(double value); + virtual void SetPresolveMode(int value); + virtual void SetLpAlgorithm(int value); + + OsiClpSolverInterface osi_; + bool quiet_; + // TODO(user): remove and query number of iterations directly from CbcModel + int64 iterations_; + int64 nodes_; + double best_objective_bound_; + // Special way to handle the relative MIP gap parameter. + double relative_mip_gap_; +}; + +// ----- Solver ----- + +// Creates a LP/MIP instance with the specified name and minimization objective. +CBCInterface::CBCInterface(MPSolver* const solver) + : MPSolverInterface(solver), + quiet_(false), iterations_(0), nodes_(0), + best_objective_bound_(-std::numeric_limits::infinity()), + relative_mip_gap_(MPSolverParameters::kDefaultRelativeMipGap) { + osi_.setStrParam(OsiProbName, solver_->name_); + osi_.setObjSense(1); +} + +CBCInterface::~CBCInterface() {} + +// Reset the solver. +void CBCInterface::Reset() { + sync_status_ = MODEL_SYNCHRONIZED; + osi_.reset(); + osi_.setObjSense(maximize_ ? -1 : 1); + osi_.setStrParam(OsiProbName, solver_->name_); +} + +void CBCInterface::ResetBestObjectiveBound() { + if (maximize_) { + best_objective_bound_ = std::numeric_limits::infinity(); + } else { + best_objective_bound_ = -std::numeric_limits::infinity(); + } +} + +void CBCInterface::SetOptimizationDirection(bool maximize) { + InvalidateSolutionSynchronization(); + if (sync_status_ == MODEL_SYNCHRONIZED) { + osi_.setObjSense(maximize ? -1 : 1); + } else { + sync_status_ = MUST_RELOAD; + } +} + +void CBCInterface::WriteModel(const string& filename) { + if (solver_->IsLPFormat(filename)) { + osi_.writeLp(filename.c_str(), ""); + } else { + // If filename does not end in ".gz", CBC will + // append ".gz" to the filename. + osi_.writeMps(filename.c_str(), ""); + } +} + +void CBCInterface::SuppressOutput() { + quiet_ = true; +} + +void CBCInterface::SetVariableBounds(int var_index, double lb, double ub) { + InvalidateSolutionSynchronization(); + if (sync_status_ == MODEL_SYNCHRONIZED) { + osi_.setColBounds(var_index, lb, ub); + } else { + sync_status_ = MUST_RELOAD; + } +} + +void CBCInterface::SetVariableInteger(int var_index, bool integer) { + InvalidateSolutionSynchronization(); + // TODO(user) : Check if this is actually a change. + if (sync_status_ == MODEL_SYNCHRONIZED) { + if (integer) { + osi_.setInteger(var_index); + } else { + osi_.setContinuous(var_index); + } + } else { + sync_status_ = MUST_RELOAD; + } +} + +void CBCInterface::SetConstraintBounds(int index, double lb, double ub) { + InvalidateSolutionSynchronization(); + if (sync_status_ == MODEL_SYNCHRONIZED) { + osi_.setRowBounds(index, lb, ub); + } else { + sync_status_ = MUST_RELOAD; + } +} + +void CBCInterface::AddRowConstraint(MPConstraint* const ct) { + sync_status_ = MUST_RELOAD; +} + +void CBCInterface::AddVariable(MPVariable* const var) { + sync_status_ = MUST_RELOAD; +} + +// Solve the LP/MIP. Returns true only if the optimal solution was revealed. +// Returns the status of the search. +MPSolver::ResultStatus CBCInterface::Solve(const MPSolverParameters& param) { + WallTimer timer; + timer.Start(); + + // Special case if the model is empty since CBC is not able to + // handle this special case by itself. + if (solver_->variables_.size() == 0 && solver_->constraints_.size() == 0) { + sync_status_ = SOLUTION_SYNCHRONIZED; + result_status_ = MPSolver::OPTIMAL; + objective_value_ = solver_->linear_objective_.offset_; + best_objective_bound_ = solver_->linear_objective_.offset_; + return result_status_; + } + + // Finish preparing the problem. + // Define variables. + switch (sync_status_) { + case MUST_RELOAD: { + Reset(); + CHECK_EQ(MODEL_SYNCHRONIZED, sync_status_); + CoinModel build; + // Create dummy variable for objective offset. + build.addColumn(0, NULL, NULL, 1.0, 1.0, + solver_->linear_objective_.offset_, "dummy", false); + const int nb_vars = solver_->variables_.size(); + for (int i = 0; i < nb_vars; ++i) { + MPVariable* const var = solver_->variables_[i]; + var->set_index(i + 1); // offset by 1 because of dummy variable. + hash_map::const_iterator it = + solver_->linear_objective_.coefficients_.find(var); + const double obj_coeff = + it == solver_->linear_objective_.coefficients_.end() ? + 0.0 : + it->second; + if (var->name().empty()) { + build.addColumn(0, NULL, NULL, var->lb(), var->ub(), obj_coeff, + NULL, var->integer()); + } else { + build.addColumn(0, NULL, NULL, var->lb(), var->ub(), obj_coeff, + var->name().c_str(), var->integer()); + } + } + + // Define constraints. + int max_row_length = 0; + int constraint_index = 0; + for (int i = 0; i < solver_->constraints_.size(); ++i) { + MPConstraint* const ct = solver_->constraints_[i]; + ct->set_index(constraint_index++); + if (ct->coefficients_.size() > max_row_length) { + max_row_length = ct->coefficients_.size(); + } + } + scoped_array indices(new int[max_row_length]); + scoped_array coefs(new double[max_row_length]); + + for (int i = 0; i < solver_->constraints_.size(); ++i) { + MPConstraint* const ct = solver_->constraints_[i]; + const int size = ct->coefficients_.size(); + int j = 0; + for (hash_map::const_iterator it = + ct->coefficients_.begin(); + it != ct->coefficients_.end(); + ++it) { + const int index = it->first->index(); + DCHECK_NE(kNoIndex, index); + indices[j] = index; + coefs[j] = it->second; + j++; + } + if (ct->name().empty()) { + build.addRow(size, indices.get(), coefs.get(), ct->lb(), ct->ub()); + } else { + build.addRow(size, indices.get(), coefs.get(), ct->lb(), ct->ub(), + ct->name().c_str()); + } + } + osi_.loadFromCoinModel(build); + break; + } + case MODEL_SYNCHRONIZED: { + break; + } + case SOLUTION_SYNCHRONIZED: { + break; + } + } + + // Changing objective sense through OSI so that the model file + // (written through OSI) has the correct objective sense + osi_.setObjSense(maximize_ ? -1 : 1); + + VLOG(1) << StringPrintf("Model built in %.3f seconds.", timer.Get()); + + WriteModelToPredefinedFiles(); + + ResetBestObjectiveBound(); + + // Solve + CbcModel model(osi_); + + if (quiet_) { + model.solver()->setHintParam(OsiDoReducePrint, true, OsiHintTry); + model.setLogLevel(-1); + } else { + model.solver()->setHintParam(OsiDoReducePrint, true, OsiHintTry); + model.setLogLevel(1); + } + + // Time limit. + if (solver_->time_limit()) { + VLOG(1) << "Setting time limit = " << solver_->time_limit() << " ms."; + model.setMaximumSeconds(solver_->time_limit() / 1000.0); + } + + // And solve. + timer.Restart(); + + // Here we use the default function from the command-line CBC solver. + // This enables to activate all the features and get the same performance + // as the CBC stand-alone executable. The syntax is ugly, however. + SetParameters(param); + // Always turn presolve on (it's the CBC default and it consistently + // improves performance). + model.setTypePresolve(0); + // Special way to set the relative MIP gap parameter as it cannot be set + // through callCbc. + model.setAllowableFractionGap(relative_mip_gap_); + int return_status = callCbc("-solve", model); + const int kBadReturnStatus = 777; + CHECK_NE(kBadReturnStatus, return_status); // Should never happen according + // to the CBC source + + VLOG(1) << StringPrintf("Solved in %.3f seconds.", timer.Get()); + + // Get the results + objective_value_ = model.getObjValue(); + VLOG(1) << "objective=" << objective_value_; + const double* const values = model.bestSolution(); + + if (values != NULL) { + // if optimal or feasible solution is found. + for (int i = 0; i < solver_->variables_.size(); ++i) { + MPVariable* const var = solver_->variables_[i]; + const int var_index = var->index(); + const double val = values[var_index]; + var->set_solution_value(val); + VLOG(3) << var->name() << "=" << val; + } + } else { + VLOG(1) << "No feasible solution found."; + } + + // Check the status: optimal, infeasible, etc. + int tmp_status = model.status(); + + VLOG(1) << "cbc result status: " << tmp_status; + /* Final status of problem + (info from cbc/v2_6_2/Cbc/src/CbcSolver.cpp) + Some of these can be found out by is...... functions + -1 before branchAndBound + 0 finished - check isProvenOptimal or isProvenInfeasible to see + if solution found + (or check value of best solution) + 1 stopped - on maxnodes, maxsols, maxtime + 2 difficulties so run was abandoned + (5 event user programmed event occurred) + */ + switch (tmp_status) { + case 0: + // Order of tests counts; if model.isContinuousUnbounded() returns true, + // then so does model.isProvenInfeasible()! + if (model.isProvenOptimal()) { + result_status_ = MPSolver::OPTIMAL; + } else if (model.isContinuousUnbounded()) { + result_status_ = MPSolver::UNBOUNDED; + } else if (model.isProvenInfeasible()) { + result_status_ = MPSolver::INFEASIBLE; + } else { + LOG(FATAL) << "Unknown solver status."; + } + break; + case 1: + result_status_ = MPSolver::FEASIBLE; + break; + default: + result_status_ = MPSolver::ABNORMAL; + break; + } + + iterations_ = model.getIterationCount(); + nodes_ = model.getNodeCount(); + best_objective_bound_ = model.getBestPossibleObjValue(); + VLOG(1) << "best objective bound=" << best_objective_bound_; + + sync_status_ = SOLUTION_SYNCHRONIZED; + + return result_status_; +} + +MPSolverInterface* BuildCBCInterface(MPSolver* const solver) { + return new CBCInterface(solver); +} + +// ------ Query statistics on the solution and the solve ------ + +int64 CBCInterface::iterations() const { + CheckSolutionIsSynchronized(); + return iterations_; +} + +int64 CBCInterface::nodes() const { + CheckSolutionIsSynchronized(); + return nodes_; +} + +double CBCInterface::best_objective_bound() const { + CheckSolutionIsSynchronized(); + CheckBestObjectiveBoundExists(); + return best_objective_bound_; +} + +// ----- Parameters ----- + +// The support for parameters in CBC is intentionally sparse. There is +// a memory leak in callCbc that prevents to pass parameters through +// it, so handling parameters would require an comprehensive rewrite +// of the code. I will improve the parameter support only if there is +// a relevant use case. + +void CBCInterface::SetParameters(const MPSolverParameters& param) { + SetCommonParameters(param); + SetMIPParameters(param); +} + +void CBCInterface::SetRelativeMipGap(double value) { + relative_mip_gap_ = value; +} + +void CBCInterface::SetPresolveMode(int value) { + SetUnsupportedIntegerParam(MPSolverParameters::PRESOLVE); +} + +void CBCInterface::SetLpAlgorithm(int value) { + SetUnsupportedIntegerParam(MPSolverParameters::LP_ALGORITHM); +} + +} // namespace operations_research +#endif // #if defined(USE_CBC) diff --git a/linear_solver/clp_interface.cc b/linear_solver/clp_interface.cc new file mode 100644 index 0000000000..f4ef89aa94 --- /dev/null +++ b/linear_solver/clp_interface.cc @@ -0,0 +1,569 @@ +// Copyright 2010 Google +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. +// + +#include "base/commandlineflags.h" +#include "base/integral_types.h" +#include "base/logging.h" +#include "base/scoped_ptr.h" +#include "base/stringprintf.h" +#include "base/timer.h" +#include "base/strutil.h" +#include "base/concise_iterator.h" +#include "base/stl_util-inl.h" + +#include "linear_solver/linear_solver.h" + +#if defined(USE_CLP) || defined(USE_CBC) + +#undef PACKAGE +#undef VERSION +#include "coin/ClpSimplex.hpp" +#include "coin/CoinBuild.hpp" +#include "coin/ClpMessage.hpp" +#include "coin/config_clp.h" + +DECLARE_double(solver_timeout_in_seconds); +DECLARE_string(solver_write_model); + +namespace operations_research { + +class CLPInterface : public MPSolverInterface { + public: + // Constructor that takes a name for the underlying CLP solver. + explicit CLPInterface(MPSolver* const solver); + ~CLPInterface(); + + // Sets the optimization direction (min/max). + virtual void SetOptimizationDirection(bool maximize); + + // ----- Solve ----- + // Solve the problem using the parameter values specified. + virtual MPSolver::ResultStatus Solve(const MPSolverParameters& param); + + // ----- Model modifications and extraction ----- + // Resets extracted model + virtual void Reset(); + + // Modify bounds. + virtual void SetVariableBounds(int var_index, double lb, double ub); + virtual void SetVariableInteger(int var_index, bool integer); + virtual void SetConstraintBounds(int row_index, double lb, double ub); + + // Add constraint incrementally. + void AddRowConstraint(MPConstraint* const ct); + // Add variable incrementally. + void AddVariable(MPVariable* const var); + // Change a coefficient in a constraint. + virtual void SetCoefficient(MPConstraint* const constraint, + MPVariable* const variable, + double coefficient); + // Clear a constraint from all its terms. + virtual void ClearConstraint(MPConstraint* const constraint); + + // Change a coefficient in the linear objective. + virtual void SetObjectiveCoefficient(MPVariable* const variable, + double coefficient); + // Change the constant term in the linear objective. + virtual void SetObjectiveOffset(double value); + // Clear the objective from all its terms (linear and quadratic) + virtual void ClearObjective(); + + + // ------ Query statistics on the solution and the solve ------ + // Number of simplex iterations + virtual int64 iterations() const; + // Number of branch-and-bound nodes. Only available for discrete problems. + virtual int64 nodes() const; + // Best objective bound. Only available for discrete problems. + virtual double best_objective_bound() const; + // Checks whether a feasible solution exists. + virtual void CheckSolutionExists() const; + + // ----- Misc ----- + // Write model + virtual void WriteModel(const string& filename); + + // SuppressOutput. + virtual void SuppressOutput(); + + // Query problem type. + virtual bool IsContinuous() const { return true; } + virtual bool IsLP() const { + return true; + } + virtual bool IsMIP() const { return false; } + + virtual void ExtractNewVariables(); + virtual void ExtractNewConstraints(); + virtual void ExtractObjective(); + + virtual string SolverVersion() const { + return PACKAGE_STRING; + } + + private: + // Create dummy variable to be able to create empty constraints. + void CreateDummyVariableForEmptyConstraints(); + + // Set all parameters in the underlying solver. + virtual void SetParameters(const MPSolverParameters& param); + // Set each parameter in the underlying solver. + virtual void SetRelativeMipGap(double value); + virtual void SetPresolveMode(int value); + virtual void SetLpAlgorithm(int value); + + scoped_ptr clp_; // TODO(user) : remove pointer. + scoped_ptr options_; // For parameter setting. +}; + +// ----- Solver ----- + +// Creates a LP/MIP instance with the specified name and minimization objective. +CLPInterface::CLPInterface(MPSolver* const solver) + : MPSolverInterface(solver), + clp_(new ClpSimplex), options_(new ClpSolve) { + clp_->setStrParam(ClpProbName, solver_->name_); + clp_->setOptimizationDirection(1); +} + +CLPInterface::~CLPInterface() {} + +void CLPInterface::Reset() { + clp_.reset(new ClpSimplex); + clp_->setOptimizationDirection(maximize_ ? -1 : 1); + ResetExtractionInformation(); +} + +void CLPInterface::WriteModel(const string& filename) { + // CLP does not support the LP format natively. It only supports it + // through OsiClpSolverInterface. + // TODO(user) : Implement support for .lp format. + if (HasSuffixString(filename, ".lp")) { + LOG(WARNING) << "CLP does not support the LP format, " + << "writing in MPS format instead."; + } + clp_->writeMps(filename.c_str()); +} + +void CLPInterface::SuppressOutput() { + clp_->setLogLevel(-1); +} + +// ------ Model modifications and extraction ----- + +// Not cached +void CLPInterface::SetOptimizationDirection(bool maximize) { + InvalidateSolutionSynchronization(); + clp_->setOptimizationDirection(maximize ? -1 : 1); +} + +void CLPInterface::SetVariableBounds(int var_index, double lb, double ub) { + InvalidateSolutionSynchronization(); + if (var_index != kNoIndex) { + // Not cached if the variable has been extracted + DCHECK_LE(var_index, last_variable_index_); + clp_->setColumnBounds(var_index, lb, ub); + } else { + sync_status_ = MUST_RELOAD; + } +} + +// Ignore as CLP does not solve models with integer variables +void CLPInterface::SetVariableInteger(int var_index, bool integer) {} + +void CLPInterface::SetConstraintBounds(int index, double lb, double ub) { + InvalidateSolutionSynchronization(); + if (index != kNoIndex) { + // Not cached if the row has been extracted + DCHECK_LE(index, last_constraint_index_); + clp_->setRowBounds(index, lb, ub); + } else { + sync_status_ = MUST_RELOAD; + } +} + +void CLPInterface::SetCoefficient(MPConstraint* const constraint, + MPVariable* const variable, + double coefficient) { + InvalidateSolutionSynchronization(); + const int constraint_index = constraint->index(); + const int variable_index = variable->index(); + if (constraint_index != kNoIndex && variable_index != kNoIndex) { + // The modification of the coefficient for an extracted row and + // variable is not cached. + DCHECK_LE(constraint_index, last_constraint_index_); + DCHECK_LE(variable_index, last_variable_index_); + clp_->modifyCoefficient(constraint_index, variable_index, coefficient); + } else { + // The modification of an unextracted row or variable is cached + // and handled in ExtractModel. + sync_status_ = MUST_RELOAD; + } +} + +// Not cached +void CLPInterface::ClearConstraint(MPConstraint* const constraint) { + InvalidateSolutionSynchronization(); + const int constraint_index = constraint->index(); + // Constraint may not have been extracted yet. + if (constraint_index != kNoIndex) { + for (ConstIter > + it(constraint->coefficients_); !it.at_end(); ++it) { + const int var_index = it->first->index(); + DCHECK_NE(kNoIndex, var_index); + clp_->modifyCoefficient(constraint_index, var_index, 0.0); + } + } +} + +// Cached +void CLPInterface::SetObjectiveCoefficient(MPVariable* const variable, + double coefficient) { + sync_status_ = MUST_RELOAD; +} + +// Cached +void CLPInterface::SetObjectiveOffset(double value) { + sync_status_ = MUST_RELOAD; +} + +// Clear objective of all its terms (linear and quadratic) +void CLPInterface::ClearObjective() { + InvalidateSolutionSynchronization(); + // Clear linear terms + for (ConstIter > + it(solver_->linear_objective_.coefficients_); + !it.at_end(); ++it) { + const int var_index = it->first->index(); + // Variable may have not been extracted yet. + if (var_index == kNoIndex) { + DCHECK_NE(MODEL_SYNCHRONIZED, sync_status_); + } else { + clp_->setObjectiveCoefficient(var_index, 0.0); + } + } + // Clear constant term. + clp_->setObjectiveOffset(0.0); + +} + + +void CLPInterface::AddRowConstraint(MPConstraint* const ct) { + sync_status_ = MUST_RELOAD; +} + +void CLPInterface::AddVariable(MPVariable* const var) { + sync_status_ = MUST_RELOAD; +} + +void CLPInterface::CreateDummyVariableForEmptyConstraints() { + clp_->setColumnBounds(kDummyVariableIndex, 0.0, 0.0); + clp_->setObjectiveCoefficient(kDummyVariableIndex, 0.0); + // Workaround for peculiar signature of setColumnName. + std::string dummy_name = "dummy"; + int var_index = kDummyVariableIndex; + clp_->setColumnName(var_index, dummy_name); +} + +// Define new variables and add them to existing constraints. +void CLPInterface::ExtractNewVariables() { + // Define new variables + int total_num_vars = solver_->variables_.size(); + if (total_num_vars > last_variable_index_) { + if (last_variable_index_ == 0 && last_constraint_index_ == 0) { + // Faster extraction when nothing has been extracted yet. + clp_->resize(0, total_num_vars + 1); + CreateDummyVariableForEmptyConstraints(); + for (int i = 0; i < total_num_vars; ++i) { + MPVariable* const var = solver_->variables_[i]; + int var_index = i + 1; // offset by 1 because of dummy variable. + var->set_index(var_index); + if (!var->name().empty()) { + std::string std_name = var->name(); + clp_->setColumnName(var_index, std_name); + } + clp_->setColumnBounds(var_index, var->lb(), var->ub()); + } + } else { + // TODO(user): This could perhaps be made slightly faster by + // iterating through old constraints, constructing by hand the + // column-major representation of the addition to them and call + // clp_->addColumns. But this is good enough for now. + // Create new variables. + for (int j = last_variable_index_; j < total_num_vars; ++j) { + MPVariable* const var = solver_->variables_[j]; + DCHECK_EQ(kNoIndex, var->index()); + int var_index = j + 1; // offset by 1 because of dummy variable. + var->set_index(var_index); + // The true objective coefficient will be set later in ExtractObjective. + double tmp_obj_coef = 0.0; + clp_->addColumn(0, NULL, NULL, var->lb(), var->ub(), tmp_obj_coef); + if (!var->name().empty()) { + std::string std_name = var->name(); + clp_->setColumnName(var_index, std_name); + } + } + // Add new variables to existing constraints. + for (int i = 0; i < last_constraint_index_; i++) { + MPConstraint* const ct = solver_->constraints_[i]; + for (ConstIter > it(ct->coefficients_); + !it.at_end(); ++it) { + const int var_index = it->first->index(); + DCHECK_NE(kNoIndex, var_index); + if (var_index >= last_variable_index_) { + clp_->modifyCoefficient(i, var_index, it->second); + } + } + } + } + } +} + +// Define new constraints on old and new variables. +void CLPInterface::ExtractNewConstraints() { + int total_num_rows = solver_->constraints_.size(); + if (last_constraint_index_ < total_num_rows) { + int max_row_length = 0; + for (int i = last_constraint_index_; i < total_num_rows; ++i) { + MPConstraint* const ct = solver_->constraints_[i]; + DCHECK_EQ(kNoIndex, ct->index()); + ct->set_index(i); + if (ct->coefficients_.size() > max_row_length) { + max_row_length = ct->coefficients_.size(); + } + } + // Make space for dummy variable. + max_row_length = max(1, max_row_length); + scoped_array indices(new int[max_row_length]); + scoped_array coefs(new double[max_row_length]); + CoinBuild build_object; + for (int i = last_constraint_index_; + i < solver_->constraints_.size(); ++i) { + MPConstraint* const ct = solver_->constraints_[i]; + DCHECK_NE(kNoIndex, ct->index()); + int size = ct->coefficients_.size(); + if (size == 0) { + // Add dummy variable to be able to build the constraint. + indices[0] = kDummyVariableIndex; + coefs[0] = 1.0; + size = 1; + } + int j = 0; + for (ConstIter > it(ct->coefficients_); + !it.at_end(); ++it) { + const int index = it->first->index(); + DCHECK_NE(kNoIndex, index); + indices[j] = index; + coefs[j] = it->second; + j++; + } + build_object.addRow(size, + indices.get(), + coefs.get(), + ct->lb(), + ct->ub()); + } + // Add and name the rows. + clp_->addRows(build_object); + for (int i = last_constraint_index_; i < total_num_rows; ++i) { + MPConstraint* const ct = solver_->constraints_[i]; + if (!ct->name().empty()) { + std::string std_name = ct->name(); + clp_->setRowName(ct->index(), std_name); + } + } + } +} + +void CLPInterface::ExtractObjective() { + // Linear objective: set objective coefficients for all variables + // (some might have been modified) + for (ConstIter > + it(solver_->linear_objective_.coefficients_); + !it.at_end(); ++it) { + clp_->setObjectiveCoefficient(it->first->index(), it->second); + } + + // Constant term. Use -offset instead of +offset because CLP does + // not follow conventions. + clp_->setObjectiveOffset(-solver_->linear_objective_.offset_); + +} + +// Extracts model and solve the LP/MIP. Returns the status of the search. +MPSolver::ResultStatus CLPInterface::Solve(const MPSolverParameters& param) { + WallTimer timer; + timer.Start(); + + // Special case if the model is empty since CLP is not able to + // handle this special case by itself. + if (solver_->variables_.size() == 0 && solver_->constraints_.size() == 0) { + sync_status_ = SOLUTION_SYNCHRONIZED; + result_status_ = MPSolver::OPTIMAL; + objective_value_ = solver_->linear_objective_.offset_; + return result_status_; + } + + ExtractModel(); + VLOG(1) << StringPrintf("Model built in %.3f seconds.", timer.Get()); + + WriteModelToPredefinedFiles(); + + // Time limit. + if (solver_->time_limit()) { + VLOG(1) << "Setting time limit = " << solver_->time_limit() << " ms."; + clp_->setMaximumSeconds(solver_->time_limit() / 1000.0); + } else { + clp_->setMaximumSeconds(-1.0); + } + + // Start from a fresh set of default parameters and set them to + // specified values. + options_.reset(new ClpSolve); + SetParameters(param); + + + // Solve + timer.Restart(); + clp_->initialSolve(*options_); + VLOG(1) << StringPrintf("Solved in %.3f seconds.", timer.Get()); + + // Get the results + objective_value_ = clp_->objectiveValue(); + VLOG(1) << "objective=" << objective_value_; + const double* const values = clp_->getColSolution(); + const double* const reduced_costs = clp_->getReducedCost(); + for (int i = 0; i < solver_->variables_.size(); ++i) { + MPVariable* const var = solver_->variables_[i]; + int var_index = var->index(); + double val = values[var_index]; + var->set_solution_value(val); + VLOG(3) << var->name() << ": value = " << val; + double reduced_cost = reduced_costs[var_index]; + var->set_reduced_cost(reduced_cost); + VLOG(4) << var->name() << ": reduced cost = " << reduced_cost; + } + const double* const dual_values = clp_->getRowPrice(); + for (int i = 0; i < solver_->constraints_.size(); ++i) { + MPConstraint* const ct = solver_->constraints_[i]; + double dual_value = 0.0; + dual_value = dual_values[i]; + ct->set_dual_value(dual_value); + VLOG(4) << "row " << ct->index() << ": dual value = " << dual_value; + } + + // Check the status: optimal, infeasible, etc. + int tmp_status = clp_->status(); + VLOG(1) << "clp result status: " << tmp_status; + switch (tmp_status) { + case CLP_SIMPLEX_FINISHED: + result_status_ = MPSolver::OPTIMAL; + break; + case CLP_SIMPLEX_INFEASIBLE: + result_status_ = MPSolver::INFEASIBLE; + break; + case CLP_SIMPLEX_UNBOUNDED: + result_status_ = MPSolver::UNBOUNDED; + break; + case CLP_SIMPLEX_STOPPED: + result_status_ = MPSolver::FEASIBLE; + break; + default: + result_status_ = MPSolver::ABNORMAL; + break; + } + + sync_status_ = SOLUTION_SYNCHRONIZED; + return result_status_; +} + +MPSolverInterface* BuildCLPInterface(MPSolver* const solver) { + return new CLPInterface(solver); +} + +// ------ Query statistics on the solution and the solve ------ + +int64 CLPInterface::iterations() const { + CheckSolutionIsSynchronized(); + return clp_->getIterationCount(); +} + +int64 CLPInterface::nodes() const { + LOG(FATAL) << "Number of nodes only available for discrete problems"; + return kUnknownNumberOfNodes; +} + +double CLPInterface::best_objective_bound() const { + LOG(FATAL) << "Best objective bound only available for discrete problems"; + return 0.0; +} + +void CLPInterface::CheckSolutionExists() const { + if (false) { + } else { + // Call default implementation + MPSolverInterface::CheckSolutionExists(); + } +} + +// ------ Parameters ------ + +void CLPInterface::SetParameters(const MPSolverParameters& param) { + SetCommonParameters(param); +} + +void CLPInterface::SetRelativeMipGap(double value) { + LOG(WARNING) << "The relative MIP gap is only available " + << "for discrete problems."; +} + +void CLPInterface::SetPresolveMode(int value) { + switch (value) { + case MPSolverParameters::PRESOLVE_OFF: { + options_->setPresolveType(ClpSolve::presolveOff); + break; + } + case MPSolverParameters::PRESOLVE_ON: { + options_->setPresolveType(ClpSolve::presolveOn); + break; + } + default: { + SetIntegerParamToUnsupportedValue(MPSolverParameters::PRESOLVE, value); + } + } +} + +void CLPInterface::SetLpAlgorithm(int value) { + switch (value) { + case MPSolverParameters::DUAL: { + options_->setSolveType(ClpSolve::useDual); + break; + } + case MPSolverParameters::PRIMAL: { + options_->setSolveType(ClpSolve::usePrimal); + break; + } + case MPSolverParameters::BARRIER: { + options_->setSolveType(ClpSolve::useBarrier); + break; + } + default: { + SetIntegerParamToUnsupportedValue(MPSolverParameters::LP_ALGORITHM, + value); + } + } +} + +} // namespace operations_research +#endif // #if defined(USE_CBC) || defined(USE_CLP) diff --git a/linear_solver/glpk_interface.cc b/linear_solver/glpk_interface.cc new file mode 100644 index 0000000000..eb4c7fe425 --- /dev/null +++ b/linear_solver/glpk_interface.cc @@ -0,0 +1,780 @@ +// Copyright 2010 Google +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. +// + +#include "base/commandlineflags.h" +#include "base/integral_types.h" +#include "base/logging.h" +#include "base/scoped_ptr.h" +#include "base/stringprintf.h" +#include "base/timer.h" +#include "base/concise_iterator.h" +#include "base/stl_util-inl.h" + +#include "linear_solver/linear_solver.h" + +#if defined(USE_GLPK) + +extern "C" { + #include "glpk.h" +} + +DECLARE_double(solver_timeout_in_seconds); +DECLARE_string(solver_write_model); + +namespace operations_research { + +// Class to store information gathered in the callback +class GLPKInformation { + public: + explicit GLPKInformation(bool maximize): num_all_nodes_(0) { + ResetBestObjectiveBound(maximize); + } + void Reset(bool maximize) { + num_all_nodes_ = 0; + ResetBestObjectiveBound(maximize); + } + void ResetBestObjectiveBound(bool maximize) { + if (maximize) { + best_objective_bound_ = std::numeric_limits::infinity(); + } else { + best_objective_bound_ = -std::numeric_limits::infinity(); + } + } + int num_all_nodes_; + double best_objective_bound_; +}; + + + +// Function to be called in the GLPK callback +void GLPKGatherInformationCallback(glp_tree* tree, void* info) { + CHECK_NOTNULL(tree); + CHECK_NOTNULL(info); + GLPKInformation* glpk_info = reinterpret_cast(info); + switch (glp_ios_reason(tree)) { + // The best bound and the number of nodes change only when GLPK + // branches, generates cuts or finds an integer solution. + case GLP_ISELECT: + case GLP_IROWGEN: + case GLP_IBINGO: { + // Get total number of nodes + glp_ios_tree_size(tree, NULL, NULL, &glpk_info->num_all_nodes_); + // Get best bound + int node_id = glp_ios_best_node(tree); + if (node_id > 0) { + glpk_info->best_objective_bound_ = glp_ios_node_bound(tree, node_id); + } + break; + } + default: + break; + } +} + +// ----- GLPK Solver ----- + +class GLPKInterface : public MPSolverInterface { + public: + // Constructor that takes a name for the underlying glpk solver. + GLPKInterface(MPSolver* const solver, bool mip); + ~GLPKInterface(); + + // Sets the optimization direction (min/max). + virtual void SetOptimizationDirection(bool maximize); + + // ----- Solve ----- + // Solve the problem using the parameter values specified. + virtual MPSolver::ResultStatus Solve(const MPSolverParameters& param); + + // ----- Model modifications and extraction ----- + // Resets extracted model + virtual void Reset(); + + // Modify bounds. + virtual void SetVariableBounds(int var_index, double lb, double ub); + virtual void SetVariableInteger(int var_index, bool integer); + virtual void SetConstraintBounds(int row_index, double lb, double ub); + + // Add Constraint incrementally. + void AddRowConstraint(MPConstraint* const ct); + // Add variable incrementally. + void AddVariable(MPVariable* const var); + // Change a coefficient in a constraint. + virtual void SetCoefficient(MPConstraint* const constraint, + MPVariable* const variable, + double coefficient); + // Clear a constraint from all its terms. + virtual void ClearConstraint(MPConstraint* const constraint); + // Change a coefficient in the linear objective + virtual void SetObjectiveCoefficient(MPVariable* const variable, + double coefficient); + // Change the constant term in the linear objective. + virtual void SetObjectiveOffset(double value); + // Clear the objective from all its terms. + virtual void ClearObjective(); + + // ------ Query statistics on the solution and the solve ------ + // Number of simplex iterations + virtual int64 iterations() const; + // Number of branch-and-bound nodes. Only available for discrete problems. + virtual int64 nodes() const; + // Best objective bound. Only available for discrete problems. + virtual double best_objective_bound() const; + // Checks whether a feasible solution exists. + virtual void CheckSolutionExists() const; + // Checks whether information on the best objective bound exists. + virtual void CheckBestObjectiveBoundExists() const; + + // ----- Misc ----- + // Write model + virtual void WriteModel(const string& filename); + + // SuppressOutput. + virtual void SuppressOutput(); + + // Query problem type. + virtual bool IsContinuous() const { return IsLP(); } + virtual bool IsLP() const { return !mip_; } + virtual bool IsMIP() const { return mip_; } + + virtual void ExtractNewVariables(); + virtual void ExtractNewConstraints(); + virtual void ExtractObjective(); + + virtual string SolverVersion() const { + return StringPrintf("GLPK %s", glp_version()); + } + + private: + // Configure the solver's parameters. + void ConfigureGLPKParameters(const MPSolverParameters& param); + + // Set all parameters in the underlying solver. + virtual void SetParameters(const MPSolverParameters& param); + // Set each parameter in the underlying solver. + virtual void SetRelativeMipGap(double value); + virtual void SetPresolveMode(int value); + virtual void SetLpAlgorithm(int value); + + void ExtractOldConstraints(); + void ExtractOneConstraint(MPConstraint* const constraint, + int* const indices, + double* const coefs); + + glp_prob* lp_; + bool mip_; + + // Parameters + glp_smcp lp_param_; + glp_iocp mip_param_; + // For the callback + scoped_ptr mip_callback_info_; +}; + +// Creates a LP/MIP instance with the specified name and minimization objective. +GLPKInterface::GLPKInterface(MPSolver* const solver, bool mip) + : MPSolverInterface(solver), lp_(NULL), mip_(mip), + mip_callback_info_(NULL) { + lp_ = glp_create_prob(); + glp_set_prob_name(lp_, solver_->name_.c_str()); + glp_set_obj_dir(lp_, GLP_MIN); + mip_callback_info_.reset(new GLPKInformation(maximize_)); +} + +// Frees the LP memory allocations. +GLPKInterface::~GLPKInterface() { + CHECK_NOTNULL(lp_); + glp_delete_prob(lp_); + lp_ = NULL; +} + +void GLPKInterface::Reset() { + CHECK_NOTNULL(lp_); + glp_delete_prob(lp_); + lp_ = glp_create_prob(); + glp_set_prob_name(lp_, solver_->name_.c_str()); + glp_set_obj_dir(lp_, maximize_ ? GLP_MAX : GLP_MIN); + ResetExtractionInformation(); +} + +void GLPKInterface::WriteModel(const string& filename) { + if (solver_->IsLPFormat(filename)) { + glp_write_lp(lp_, NULL, filename.c_str()); + } else { + glp_write_mps(lp_, GLP_MPS_DECK, NULL, filename.c_str()); + } +} + +static int no_print_hook(void *info, const char *s) { + return 1; +} + +void GLPKInterface::SuppressOutput() { + glp_term_hook(no_print_hook, NULL); +} + +// ------ Model modifications and extraction ----- + +// Not cached +void GLPKInterface::SetOptimizationDirection(bool maximize) { + InvalidateSolutionSynchronization(); + glp_set_obj_dir(lp_, maximize ? GLP_MAX : GLP_MIN); +} + +void GLPKInterface::SetVariableBounds(int var_index, double lb, double ub) { + InvalidateSolutionSynchronization(); + if (var_index != kNoIndex) { + // Not cached if the variable has been extracted. + DCHECK(lp_ != NULL); + const double infinity = solver_->infinity(); + if (lb != -infinity) { + if (ub != infinity) { + if (lb == ub) { + glp_set_col_bnds(lp_, var_index, GLP_FX, lb, ub); + } else { + glp_set_col_bnds(lp_, var_index, GLP_DB, lb, ub); + } + } else { + glp_set_col_bnds(lp_, var_index, GLP_LO, lb, 0.0); + } + } else if (ub != infinity) { + glp_set_col_bnds(lp_, var_index, GLP_UP, 0.0, ub); + } else { + glp_set_col_bnds(lp_, var_index, GLP_FR, 0.0, 0.0); + } + } else { + sync_status_ = MUST_RELOAD; + } +} + +void GLPKInterface::SetVariableInteger(int var_index, bool integer) { + InvalidateSolutionSynchronization(); + if (mip_) { + if (var_index != kNoIndex) { + // Not cached if the variable has been extracted. + glp_set_col_kind(lp_, var_index, integer ? GLP_IV : GLP_CV); + } else { + sync_status_ = MUST_RELOAD; + } + } +} + +void GLPKInterface::SetConstraintBounds(int index, double lb, double ub) { + InvalidateSolutionSynchronization(); + if (index != kNoIndex) { + // Not cached if the row has been extracted + DCHECK(lp_ != NULL); + const double infinity = solver_->infinity(); + if (lb != -infinity) { + if (ub != infinity) { + if (lb == ub) { + glp_set_row_bnds(lp_, index, GLP_FX, lb, ub); + } else { + glp_set_row_bnds(lp_, index, GLP_DB, lb, ub); + } + } else { + glp_set_row_bnds(lp_, index, GLP_LO, lb, 0.0); + } + } else if (ub != infinity) { + glp_set_row_bnds(lp_, index, GLP_UP, 0.0, ub); + } else { + glp_set_row_bnds(lp_, index, GLP_FR, 0.0, 0.0); + } + } else { + sync_status_ = MUST_RELOAD; + } +} + +void GLPKInterface::SetCoefficient(MPConstraint* const constraint, + MPVariable* const variable, + double coefficient) { + InvalidateSolutionSynchronization(); + // GLPK does not allow to modify one coefficient at a time, so we + // extract the whole constraint again, if it has been extracted + // already and if it does not contain new variables. Otherwise, we + // cache the modification. + if (constraint->index() != kNoIndex && + (sync_status_ == MODEL_SYNCHRONIZED || + !constraint->ContainsNewVariables())) { + const int size = constraint->coefficients_.size(); + scoped_array indices(new int[size + 1]); + scoped_array coefs(new double[size + 1]); + ExtractOneConstraint(constraint, indices.get(), coefs.get()); + } +} + +// Not cached +void GLPKInterface::ClearConstraint(MPConstraint* const constraint) { + InvalidateSolutionSynchronization(); + const int constraint_index = constraint->index(); + // Constraint may have not been extracted yet. + if (constraint_index != kNoIndex) { + glp_set_mat_row(lp_, constraint_index, 0, NULL, NULL); + } +} + +// Cached +void GLPKInterface::SetObjectiveCoefficient(MPVariable* const variable, + double coefficient) { + sync_status_ = MUST_RELOAD; +} + +// Cached +void GLPKInterface::SetObjectiveOffset(double value) { + sync_status_ = MUST_RELOAD; +} + +// Clear objective of all its terms (linear) +void GLPKInterface::ClearObjective() { + InvalidateSolutionSynchronization(); + for (ConstIter > + it(solver_->linear_objective_.coefficients_); + !it.at_end(); ++it) { + const int var_index = it->first->index(); + // Variable may have not been extracted yet. + if (var_index == kNoIndex) { + DCHECK_NE(MODEL_SYNCHRONIZED, sync_status_); + } else { + glp_set_obj_coef(lp_, var_index, 0.0); + } + } + // Constant term. + glp_set_obj_coef(lp_, 0, 0.0); +} + +void GLPKInterface::AddRowConstraint(MPConstraint* const ct) { + sync_status_ = MUST_RELOAD; +} + +void GLPKInterface::AddVariable(MPVariable* const var) { + sync_status_ = MUST_RELOAD; +} + +// Define new variables and add them to existing constraints. +void GLPKInterface::ExtractNewVariables() { + int total_num_vars = solver_->variables_.size(); + if (total_num_vars > last_variable_index_) { + glp_add_cols(lp_, total_num_vars - last_variable_index_); + for (int j = last_variable_index_; j < solver_->variables_.size(); ++j) { + MPVariable* const var = solver_->variables_[j]; + // GLPK convention is to start indexing at 1. + const int var_index = j + 1; + var->set_index(var_index); + if (!var->name().empty()) { + glp_set_col_name(lp_, var_index, var->name().c_str()); + } + SetVariableBounds(var->index(), var->lb(), var->ub()); + SetVariableInteger(var->index(), var->integer()); + + // The true objective coefficient will be set later in ExtractObjective. + double tmp_obj_coef = 0.0; + glp_set_obj_coef(lp_, var->index(), tmp_obj_coef); + } + // Add new variables to the existing constraints. + ExtractOldConstraints(); + } +} + +// Extract again existing constraints if they contain new variables. +void GLPKInterface::ExtractOldConstraints() { + int max_constraint_size = solver_->ComputeMaxConstraintSize( + 0, last_constraint_index_); + // The first entry in the following arrays is dummy, to be + // consistent with glpk API. + scoped_array indices(new int[max_constraint_size + 1]); + scoped_array coefs(new double[max_constraint_size + 1]); + + for (int i = 0; i < last_constraint_index_; ++i) { + MPConstraint* const ct = solver_->constraints_[i]; + DCHECK_NE(kNoIndex, ct->index()); + const int size = ct->coefficients_.size(); + if (size == 0) { + continue; + } + // Update the constraint's coefficients if it contains new variables. + if (ct->ContainsNewVariables()) { + ExtractOneConstraint(ct, indices.get(), coefs.get()); + } + } +} + +// Extract one constraint. Arrays indices and coefs must be +// preallocated to have enough space to contain the constraint's +// coefficients. +void GLPKInterface::ExtractOneConstraint(MPConstraint* const constraint, + int* const indices, + double* const coefs) { + // GLPK convention is to start indexing at 1. + int k = 1; + for (ConstIter > it(constraint->coefficients_); + !it.at_end(); ++it) { + const int var_index = it->first->index(); + DCHECK_NE(kNoIndex, var_index); + indices[k] = var_index; + coefs[k] = it->second; + ++k; + } + glp_set_mat_row(lp_, constraint->index(), k-1, indices, coefs); +} + +// Define new constraints on old and new variables. +void GLPKInterface::ExtractNewConstraints() { + int total_num_rows = solver_->constraints_.size(); + if (last_constraint_index_ < total_num_rows) { + // Define new constraints + glp_add_rows(lp_, total_num_rows - last_constraint_index_); + int num_coefs = 0; + for (int i = last_constraint_index_; i < total_num_rows; ++i) { + // GLPK convention is to start indexing at 1. + const int constraint_index = i + 1; + MPConstraint* ct = solver_->constraints_[i]; + ct->set_index(constraint_index); + if (ct->name().empty()) { + glp_set_row_name(lp_, constraint_index, + StringPrintf("ct_%i", i).c_str()); + } else { + glp_set_row_name(lp_, constraint_index, ct->name().c_str()); + } + // All constraints are set to be of the type <= limit_ . + SetConstraintBounds(constraint_index, ct->lb(), ct->ub()); + num_coefs += ct->coefficients_.size(); + } + + // Fill new constraints with coefficients + if (last_variable_index_ == 0 && last_constraint_index_ == 0) { + // Faster extraction when nothing has been extracted yet: build + // and load whole matrix at once instead of constructing rows + // separately. + + // The first entry in the following arrays is dummy, to be + // consistent with glpk API. + scoped_array variable_indices(new int[num_coefs + 1]); + scoped_array constraint_indices(new int[num_coefs + 1]); + scoped_array coefs(new double[num_coefs + 1]); + int k = 1; + for (int i = 0; i < solver_->constraints_.size(); ++i) { + MPConstraint* ct = solver_->constraints_[i]; + for (hash_map::const_iterator it = + ct->coefficients_.begin(); + it != ct->coefficients_.end(); + ++it) { + DCHECK_NE(kNoIndex, it->first->index()); + constraint_indices[k] = ct->index(); + variable_indices[k] = it->first->index(); + coefs[k] = it->second; + ++k; + } + } + CHECK_EQ(num_coefs + 1, k); + glp_load_matrix(lp_, num_coefs, constraint_indices.get(), + variable_indices.get(), coefs.get()); + } else { + // Build each new row separately. + int max_constraint_size = solver_->ComputeMaxConstraintSize( + last_constraint_index_, total_num_rows); + // The first entry in the following arrays is dummy, to be + // consistent with glpk API. + scoped_array indices(new int[max_constraint_size + 1]); + scoped_array coefs(new double[max_constraint_size + 1]); + for (int i = last_constraint_index_; i < total_num_rows; i++) { + ExtractOneConstraint(solver_->constraints_[i], indices.get(), + coefs.get()); + } + } + } +} + +void GLPKInterface::ExtractObjective() { + // Linear objective: set objective coefficients for all variables + // (some might have been modified). + for (hash_map::const_iterator it = + solver_->linear_objective_.coefficients_.begin(); + it != solver_->linear_objective_.coefficients_.end(); + ++it) { + glp_set_obj_coef(lp_, it->first->index(), it->second); + } + // Constant term. + glp_set_obj_coef(lp_, 0, solver_->linear_objective_.offset_); +} + +// Solve the problem using the parameter values specified. +MPSolver::ResultStatus GLPKInterface::Solve(const MPSolverParameters& param) { + WallTimer timer; + timer.Start(); + + ExtractModel(); + VLOG(1) << StringPrintf("Model built in %.3f seconds.", timer.Get()); + + WriteModelToPredefinedFiles(); + + // Configure parameters at every solve, even when the model has not + // been changed, in case some of the parameters such as the time + // limit have been changed since the last solve. + ConfigureGLPKParameters(param); + + // Solve + timer.Restart(); + if (mip_) { + // glp_intopt requires to solve the root LP separately. + int simplex_status = glp_simplex(lp_, &lp_param_); + // If the root LP was solved successfully, solve the MIP. + if (simplex_status == 0) { + glp_intopt(lp_, &mip_param_); + } else { + // Something abnormal occurred during the root LP solve. It is + // highly unlikely that an integer feasible solution is + // available at this point, so we don't put any effort in trying + // to recover it. + result_status_ = MPSolver::ABNORMAL; + sync_status_ = SOLUTION_SYNCHRONIZED; + return result_status_; + } + } else { + glp_simplex(lp_, &lp_param_); + } + VLOG(1) << StringPrintf("Solved in %.3f seconds.", timer.Get()); + + // Get the results. + if (mip_) { + objective_value_ = glp_mip_obj_val(lp_); + } else { + objective_value_ = glp_get_obj_val(lp_); + } + VLOG(1) << "objective=" << objective_value_; + for (int i = 0; i < solver_->variables_.size(); ++i) { + MPVariable* const var = solver_->variables_[i]; + double val; + if (mip_) { + val = glp_mip_col_val(lp_, var->index()); + } else { + val = glp_get_col_prim(lp_, var->index()); + } + var->set_solution_value(val); + VLOG(3) << var->name() << ": value =" << val; + if (!mip_) { + double reduced_cost; + reduced_cost = glp_get_col_dual(lp_, var->index()); + var->set_reduced_cost(reduced_cost); + VLOG(4) << var->name() << ": reduced cost = " << reduced_cost; + } + } + if (!mip_) { + for (int i = 0; i < solver_->constraints_.size(); ++i) { + MPConstraint* const ct = solver_->constraints_[i]; + double dual_value = glp_get_row_dual(lp_, ct->index()); + ct->set_dual_value(dual_value); + VLOG(4) << "row " << ct->index() << ": dual value = " << dual_value; + } + } + + // Check the status: optimal, infeasible, etc. + if (mip_) { + int tmp_status = glp_mip_status(lp_); + VLOG(1) << "gplk result status: " << tmp_status; + if (tmp_status == GLP_OPT) { + result_status_ = MPSolver::OPTIMAL; + } else if (tmp_status == GLP_FEAS) { + result_status_ = MPSolver::FEASIBLE; + } else if (tmp_status == GLP_NOFEAS) { + // For infeasible problems, GLPK actually seems to return + // GLP_UNDEF. So this is never (?) reached. Return infeasible + // in case GLPK returns a correct status in future versions. + result_status_ = MPSolver::INFEASIBLE; + } else { + result_status_ = MPSolver::ABNORMAL; + // GLPK does not have a status code for unbounded MIP models, so + // we return an abnormal status instead. + } + } else { + int tmp_status = glp_get_status(lp_); + VLOG(1) << "gplk result status: " << tmp_status; + if (tmp_status == GLP_OPT) { + result_status_ = MPSolver::OPTIMAL; + } else if (tmp_status == GLP_FEAS) { + result_status_ = MPSolver::FEASIBLE; + } else if (tmp_status == GLP_NOFEAS || + tmp_status == GLP_INFEAS) { + // For infeasible problems, GLPK actually seems to return + // GLP_UNDEF. So this is never (?) reached. Return infeasible + // in case GLPK returns a correct status in future versions. + result_status_ = MPSolver::INFEASIBLE; + } else if (tmp_status == GLP_UNBND) { + // For unbounded problems, GLPK actually seems to return + // GLP_UNDEF. So this is never (?) reached. Return unbounded + // in case GLPK returns a correct status in future versions. + result_status_ = MPSolver::UNBOUNDED; + } else { + result_status_ = MPSolver::ABNORMAL; + } + } + + sync_status_ = SOLUTION_SYNCHRONIZED; + + return result_status_; +} + +MPSolverInterface* BuildGLPKInterface(MPSolver* const solver, bool mip) { + return new GLPKInterface(solver, mip); +} + +// ------ Query statistics on the solution and the solve ------ + +int64 GLPKInterface::iterations() const { + if (mip_) { + LOG(WARNING) << "Total number of iterations is not available"; + return kUnknownNumberOfIterations; + } else { + CheckSolutionIsSynchronized(); + return lpx_get_int_parm(lp_, LPX_K_ITCNT); + } +} + +int64 GLPKInterface::nodes() const { + if (mip_) { + CheckSolutionIsSynchronized(); + return mip_callback_info_->num_all_nodes_; + } else { + LOG(FATAL) << "Number of nodes only available for discrete problems"; + return kUnknownNumberOfNodes; + } +} + +double GLPKInterface::best_objective_bound() const { + if (mip_) { + CheckSolutionIsSynchronized(); + CheckBestObjectiveBoundExists(); + if (solver_->variables_.size() == 0 && solver_->constraints_.size() == 0) { + // Special case for empty model. + return solver_->linear_objective_.offset_; + } else { + return mip_callback_info_->best_objective_bound_; + } + } else { + LOG(FATAL) << "Best objective bound only available for discrete problems"; + return 0.0; + } +} + +void GLPKInterface::CheckSolutionExists() const { + if (result_status_ == MPSolver::ABNORMAL) { + LOG(WARNING) << "Ignoring ABNORMAL status from GLPK: This status may or may" + << " not indicate that a solution exists."; + } else { + // Call default implementation + MPSolverInterface::CheckSolutionExists(); + } +} + +void GLPKInterface::CheckBestObjectiveBoundExists() const { + if (result_status_ == MPSolver::ABNORMAL) { + LOG(WARNING) << "Ignoring ABNORMAL status from GLPK: This status may or may" + << " not indicate that information is available on the best" + << " objective bound."; + } else { + // Call default implementation + MPSolverInterface::CheckBestObjectiveBoundExists(); + } +} + +// ------ Parameters ------ + +void GLPKInterface::ConfigureGLPKParameters(const MPSolverParameters& param) { + if (mip_) { + glp_init_iocp(&mip_param_); + // Time limit + if (solver_->time_limit()) { + VLOG(1) << "Setting time limit = " << solver_->time_limit() << " ms."; + mip_param_.tm_lim = solver_->time_limit(); + } + // Initialize structures related to the callback. + mip_param_.cb_func = GLPKGatherInformationCallback; + mip_callback_info_->Reset(maximize_); + mip_param_.cb_info = mip_callback_info_.get(); + // TODO(user): switch some cuts on? All cuts are off by default!? + } + + // Configure LP parameters in all cases since they will be used to + // solve the root LP in the MIP case. + glp_init_smcp(&lp_param_); + // Time limit + if (solver_->time_limit()) { + VLOG(1) << "Setting time limit = " << solver_->time_limit() << " ms."; + lp_param_.tm_lim = solver_->time_limit(); + } + + // Should give a numerically better representation of the problem. + glp_scale_prob(lp_, GLP_SF_AUTO); + + // Use advanced initial basis (options: standard / advanced / Bixby's). + glp_adv_basis(lp_, NULL); + + // Set parameters specified by the user. + SetParameters(param); +} + +void GLPKInterface::SetParameters(const MPSolverParameters& param) { + SetCommonParameters(param); + if (mip_) { + SetMIPParameters(param); + } +} + +void GLPKInterface::SetRelativeMipGap(double value) { + if (mip_) { + mip_param_.mip_gap = value; + } else { + LOG(WARNING) << "The relative MIP gap is only available " + << "for discrete problems."; + } +} + +void GLPKInterface::SetPresolveMode(int value) { + switch (value) { + case MPSolverParameters::PRESOLVE_OFF: { + mip_param_.presolve = GLP_OFF; + lp_param_.presolve = GLP_OFF; + break; + } + case MPSolverParameters::PRESOLVE_ON: { + mip_param_.presolve = GLP_ON; + lp_param_.presolve = GLP_ON; + break; + } + default: { + SetIntegerParamToUnsupportedValue(MPSolverParameters::PRESOLVE, value); + } + } +} + +void GLPKInterface::SetLpAlgorithm(int value) { + switch (value) { + case MPSolverParameters::DUAL: { + // Use dual, and if it fails, switch to primal. + lp_param_.meth = GLP_DUALP; + break; + } + case MPSolverParameters::PRIMAL: { + lp_param_.meth = GLP_PRIMAL; + break; + } + case MPSolverParameters::BARRIER: + default: { + SetIntegerParamToUnsupportedValue(MPSolverParameters::LP_ALGORITHM, + value); + } + } +} + +} // namespace operations_research +#endif // #if defined(USE_GLPK) diff --git a/linear_solver/linear_solver.cc b/linear_solver/linear_solver.cc new file mode 100644 index 0000000000..eff4171606 --- /dev/null +++ b/linear_solver/linear_solver.cc @@ -0,0 +1,815 @@ +// Copyright 2010 Google +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. +// +// (Laurent Perron) +// +// A C++ wrapper around the GNU linear programming kit and Coin LP. + +// This is a complete TODO list: (* is higher priority) +// - support primal/dual/barrier for all algorithms that supports them. +// - support write model in different format (at least MPS and LP). +// * Better support for warm restart (store/restore basis?). +// - Compile Coin LP with extra packages +// * Barrier code: see http://www.cise.ufl.edu/research/sparse/amd/ +// * : see http://www.cise.ufl.edu/research/sparse/cholmod/ +// - Include support for more type of constraints: +// * logical ones with linearization. +// * SOS1 and SOS2. +// - Support for semi-continuous variables. +// * Support for absolute and relative optimization steps. +// - Support for generalized convex cost. +// - Support for MIP callbacks if available. +// * Implement MakeXXXVarArray as on Solver CP class. +// * Implement optimized model modification methods (AddVariable) +// to avoid reloading the full model. +// * Make sure the simplex can warmstart on the previous solution if the model +// has been changed incrementally. + +#include "linear_solver/linear_solver.h" + +#include +#include + +#include "base/commandlineflags.h" +#include "base/integral_types.h" +#include "base/logging.h" +#include "base/scoped_ptr.h" +#include "base/stringprintf.h" +#include "base/timer.h" +#include "base/concise_iterator.h" +#include "base/map-util.h" +#include "base/stl_util-inl.h" +#include "linear_solver/linear_solver.pb.h" + +DEFINE_string(solver_write_model, "", "path of the file to write the model to"); + +namespace { + +// Insert name in name_set and check for duplicates. +void CheckDuplicateName(const string& name, + hash_set* name_set) { + if (!name.empty()) { + pair::iterator, bool> result = name_set->insert(name); + if (!result.second) { + LOG(FATAL) << "Duplicate name: " << name; + } + } +} + +} // namespace + +namespace operations_research { + +// ----- MPConstraint ----- + +void MPConstraint::AddTerm(MPVariable* const var, double coeff) { + CHECK_NOTNULL(var); + if (coeff != 0.0) { + double* coefficient = FindOrNull(coefficients_, var); + if (coefficient != NULL) { + (*coefficient) += coeff; + interface_->SetCoefficient(this, var, *coefficient); + } else { + coefficients_[var] = coeff; + interface_->SetCoefficient(this, var, coeff); + } + } +} + +void MPConstraint::AddTerm(MPVariable* const var) { + AddTerm(var, 1.0); +} + +void MPConstraint::SetCoefficient(MPVariable* const var, double coeff) { + CHECK_NOTNULL(var); + coefficients_[var] = coeff; + interface_->SetCoefficient(this, var, coeff); +} + +void MPConstraint::Clear() { + interface_->ClearConstraint(this); + coefficients_.clear(); +} + +void MPConstraint::SetBounds(double lb, double ub) { + const bool change = lb != lb_ || ub != ub_; + lb_ = lb; + ub_ = ub; + if (index_ != MPSolverInterface::kNoIndex && change) { + interface_->SetConstraintBounds(index_, lb_, ub_); + } +} + +double MPConstraint::dual_value() const { + CHECK(interface_->IsContinuous()) << + "Dual value only available for continuous problems"; + interface_->CheckSolutionIsSynchronized(); + interface_->CheckSolutionExists(); + return dual_value_; +} + +bool MPConstraint::ContainsNewVariables() { + const int last_variable_index = interface_->last_variable_index(); + for (ConstIter > it(coefficients_); + !it.at_end(); ++it) { + const int variable_index = it->first->index(); + if (variable_index >= last_variable_index || + variable_index == MPSolverInterface::kNoIndex) { + return true; + } + } + return false; +} + +// ----- MPObjective ----- + +void MPObjective::AddTerm(MPVariable* const var, double coeff) { + CHECK_NOTNULL(var); + if (coeff != 0.0) { + double* coefficient = FindOrNull(coefficients_, var); + if (coefficient != NULL) { + (*coefficient) += coeff; + interface_->SetObjectiveCoefficient(var, *coefficient); + } else { + coefficients_[var] = coeff; + interface_->SetObjectiveCoefficient(var, coeff); + } + } +} + +void MPObjective::AddTerm(MPVariable* const var) { + AddTerm(var, 1.0); +} + +void MPObjective::SetCoefficient(MPVariable* const var, double coeff) { + CHECK_NOTNULL(var); + coefficients_[var] = coeff; + interface_->SetObjectiveCoefficient(var, coeff); +} + +void MPObjective::AddOffset(double value) { + offset_ += value; + interface_->SetObjectiveOffset(offset_); +} + +void MPObjective::SetOffset(double value) { + offset_ = value; + interface_->SetObjectiveOffset(offset_); +} + +void MPObjective::Clear() { + interface_->ClearObjective(); + coefficients_.clear(); + offset_ = 0.0; +} + +// ----- MPVariable ----- + +double MPVariable::solution_value() const { + interface_->CheckSolutionIsSynchronized(); + interface_->CheckSolutionExists(); + return solution_value_; +} + +double MPVariable::reduced_cost() const { + CHECK(interface_->IsContinuous()) << + "Reduced cost only available for continuous problems"; + interface_->CheckSolutionIsSynchronized(); + interface_->CheckSolutionExists(); + return reduced_cost_; +} + +void MPVariable::SetBounds(double lb, double ub) { + const bool change = lb != lb_ || ub != ub_; + lb_ = lb; + ub_ = ub; + if (index_ != MPSolverInterface::kNoIndex && change) { + interface_->SetVariableBounds(index_, lb_, ub_); + } +} + +void MPVariable::SetInteger(bool integer) { + if (integer_ != integer) { + integer_ = integer; + if (index_ != MPSolverInterface::kNoIndex) { + interface_->SetVariableInteger(index_, integer); + } + } +} + + +// ----- Objective ----- + +double MPSolver::objective_value() const { + return interface_->objective_value(); +} + +double MPSolver::best_objective_bound() const { + return interface_->best_objective_bound(); +} + +void MPSolver::ClearObjective() { + linear_objective_.Clear(); +} + +void MPSolver::AddObjectiveTerm(MPVariable* const var, double coeff) { + linear_objective_.AddTerm(var, coeff); +} + +void MPSolver::AddObjectiveTerm(MPVariable* const var) { + linear_objective_.AddTerm(var); +} + +void MPSolver::SetObjectiveCoefficient(MPVariable* const var, double coeff) { + linear_objective_.SetCoefficient(var, coeff); +} + +void MPSolver::AddObjectiveOffset(double value) { + linear_objective_.AddOffset(value); +} + +void MPSolver::SetObjectiveOffset(double value) { + linear_objective_.SetOffset(value); +} + + +void MPSolver::SetOptimizationDirection(bool maximize) { + interface_->maximize_ = maximize; + interface_->SetOptimizationDirection(maximize); +} + // Minimizing or maximizing? +bool MPSolver::Maximization() { + return interface_->maximize_; +} + +bool MPSolver::Minimization() { + return !interface_->maximize_; +} + +// ----- Version ----- + +string MPSolver::SolverVersion() const { + return interface_->SolverVersion(); +} + +// ----- Solver ----- + +#if defined(USE_CLP) || defined(USE_CBC) +extern MPSolverInterface* BuildCLPInterface(MPSolver* const solver); +#endif +#if defined(USE_CBC) +extern MPSolverInterface* BuildCBCInterface(MPSolver* const solver); +#endif +#if defined(USE_GLPK) +extern MPSolverInterface* BuildGLPKInterface(MPSolver* const solver, bool mip); +#endif + +const int64 MPSolverInterface::kUnknownNumberOfIterations = -1; +const int64 MPSolverInterface::kUnknownNumberOfNodes = -1; +const int MPSolverInterface::kNoIndex = -1; + +namespace { +MPSolverInterface* BuildSolverInterface( + MPSolver* const solver, MPSolver::OptimizationProblemType problem_type) { + switch (problem_type) { +#if defined(USE_GLPK) + case MPSolver::GLPK_LINEAR_PROGRAMMING: + return BuildGLPKInterface(solver, false); + case MPSolver::GLPK_MIXED_INTEGER_PROGRAMMING: + return BuildGLPKInterface(solver, true); +#endif +#if defined(USE_CLP) || defined(USE_CBC) + case MPSolver::CLP_LINEAR_PROGRAMMING: + return BuildCLPInterface(solver); +#endif +#if defined(USE_CBC) + case MPSolver::CBC_MIXED_INTEGER_PROGRAMMING: + return BuildCBCInterface(solver); +#endif + default: + LOG(FATAL) << "Linear solver not recognized."; + } + return NULL; +} +} // namespace + +// Creates a LP/MIP instance with the specified name and minimization objective. +MPSolver::MPSolver(const string& name, OptimizationProblemType problem_type) + : name_(name), + interface_(BuildSolverInterface(this, problem_type)), + linear_objective_(interface_.get()), + time_limit_(0.0), + write_model_filename_("") { + timer_.Restart(); +} + +// Frees the LP memory allocations. +MPSolver::~MPSolver() { + Clear(); +} + +// ----- Names management ----- + +bool MPSolver::CheckNameValidity(const string& name) { + // Allow names that conform to the LP and MPS format. + const int kMaxNameLength = 255; + if (name.size() > kMaxNameLength) { + LOG(WARNING) << "Invalid name " << name + << ": length > " << kMaxNameLength << "." + << " Will be unable to write model to file."; + return false; + } + if (name.find_first_of(" +-*<>=:\\") != string::npos) { + LOG(WARNING) << "Invalid name " << name + << ": contains forbidden character: +-*<>=:\\ space." + << " Will be unable to write model to file."; + return false; + } + size_t first_occurrence = name.find_first_of(".0123456789"); + if (first_occurrence != string::npos && first_occurrence == 0) { + LOG(WARNING) << "Invalid name " << name + << ": first character should not be . or a number." + << " Will be unable to write model to file."; + return false; + } + return true; +} + +bool MPSolver::CheckAllNamesValidity() { + for (int i = 0; i < variables_.size(); ++i) { + if (!CheckNameValidity(variables_[i]->name())) { + return false; + } + } + for (int i = 0; i < constraints_.size(); ++i) { + if (!CheckNameValidity(constraints_[i]->name())) { + return false; + } + } + return true; +} + +// ----- Load from protobuf ----- +MPSolver::LoadStatus MPSolver::Load(const MPModelProto& model) { + hash_map variables; + for (int i = 0; i < model.variables_size(); ++i) { + const MPVariableProto& var_proto = model.variables(i); + const string& id = var_proto.id(); + if (!ContainsKey(variables, id)) { + MPVariable* variable = MakeNumVar(var_proto.lb(), var_proto.ub(), id); + variable->SetInteger(var_proto.integer()); + variables[id] = variable; + } else { + return MPSolver::DUPLICATE_VARIABLE_ID; + } + } + for (int i = 0; i < model.constraints_size(); ++i) { + const MPConstraintProto& ct_proto = model.constraints(i); + const string& ct_id = ct_proto.has_id() ? ct_proto.id() : ""; + MPConstraint* const ct = MakeRowConstraint(ct_proto.lb(), + ct_proto.ub(), + ct_id); + for (int j = 0; j < ct_proto.terms_size(); ++j) { + const MPTermProto& term_proto = ct_proto.terms(j); + const string& id = term_proto.variable_id(); + MPVariable* variable = FindPtrOrNull(variables, id); + if (variable != NULL) { + ct->AddTerm(variable, term_proto.coefficient()); + } else { + return MPSolver::UNKNOWN_VARIABLE_ID; + } + } + } + for (int i = 0; i < model.objective_terms_size(); ++i) { + const MPTermProto& term_proto = model.objective_terms(i); + const string& id = term_proto.variable_id(); + MPVariable* variable = FindPtrOrNull(variables, id); + if (variable != NULL) { + AddObjectiveTerm(variable, term_proto.coefficient()); + } else { + return MPSolver::UNKNOWN_VARIABLE_ID; + } + } + SetOptimizationDirection(model.maximize()); + return MPSolver::NO_ERROR; +} + +void MPSolver::Clear() { + ClearObjective(); + STLDeleteElements(&variables_); + STLDeleteElements(&constraints_); + variables_.clear(); + constraints_.clear(); + interface_->Reset(); + SetMinimization(); +} + +void MPSolver::Reset() { + interface_->Reset(); +} + +void MPSolver::SuppressOutput() { + interface_->SuppressOutput(); +} + +MPVariable* MPSolver::MakeVar( + double lb, double ub, bool integer, const string& name) { + CheckNameValidity(name); + CheckDuplicateName(name, &variables_names_); + MPVariable* v = new MPVariable(lb, ub, integer, name, interface_.get()); + variables_.push_back(v); + interface_->AddVariable(v); + return v; +} + +MPVariable* MPSolver::MakeNumVar( + double lb, double ub, const string& name) { + return MakeVar(lb, ub, false, name); +} + +MPVariable* MPSolver::MakeIntVar( + double lb, double ub, const string& name) { + return MakeVar(lb, ub, true, name); +} + +MPVariable* MPSolver::MakeBoolVar(const string& name) { + return MakeVar(0.0, 1.0, true, name); +} + +void MPSolver::MakeVarArray(int nb, + double lb, + double ub, + bool integer, + const string& name, + vector* vars) { + CHECK_GE(nb, 0); + for (int i = 0; i < nb; ++i) { + if (name.empty()) { + vars->push_back(MakeVar(lb, ub, integer, name)); + } else { + string vname = StringPrintf("%s%d", name.c_str(), i); + vars->push_back(MakeVar(lb, ub, integer, vname)); + } + } +} + +void MPSolver::MakeNumVarArray(int nb, + double lb, + double ub, + const string& name, + vector* vars) { + MakeVarArray(nb, lb, ub, false, name, vars); +} + +void MPSolver::MakeIntVarArray(int nb, + double lb, + double ub, + const string& name, + vector* vars) { + MakeVarArray(nb, lb, ub, true, name, vars); +} + +void MPSolver::MakeBoolVarArray(int nb, + const string& name, + vector* vars) { + MakeVarArray(nb, 0.0, 1.0, true, name, vars); +} + +// Creates a new row constraint, adds it to the LP/MIP and returns it. +// MPSolver owns the Constraint. Do not free the memory yourself. +MPConstraint* MPSolver::MakeRowConstraint(double lb, double ub) { + return MakeRowConstraint(lb, ub, ""); +} + +MPConstraint* MPSolver::MakeRowConstraint() { + return MakeRowConstraint(-infinity(), infinity(), ""); +} + +// Creates a new row constraint, adds it to the LP/MIP and returns it. +// MPSolver owns the Constraint. Do not free the memory yourself. +MPConstraint* MPSolver::MakeRowConstraint(double lb, double ub, + const string& name) { + CheckNameValidity(name); + CheckDuplicateName(name, &constraints_names_); + MPConstraint* const constraint = + new MPConstraint(lb, ub, name, interface_.get()); + constraints_.push_back(constraint); + interface_->AddRowConstraint(constraint); + return constraint; +} + +MPConstraint* MPSolver::MakeRowConstraint(const string& name) { + return MakeRowConstraint(-infinity(), infinity(), name); +} + +// Compute the size of the constraint with the largest number of +// coefficients with index in [min_constraint_index, +// max_constraint_index) +int MPSolver::ComputeMaxConstraintSize(int min_constraint_index, + int max_constraint_index) const { + int max_constraint_size = 0; + DCHECK_GE(min_constraint_index, 0); + DCHECK_LE(max_constraint_index, constraints_.size()); + for (int i = min_constraint_index; i < max_constraint_index; ++i) { + MPConstraint* const ct = constraints_[i]; + if (ct->coefficients_.size() > max_constraint_size) { + max_constraint_size = ct->coefficients_.size(); + } + } + return max_constraint_size; +} + +MPSolver::ResultStatus MPSolver::Solve() { + return interface_->Solve(); +} + +MPSolver::ResultStatus MPSolver::Solve(const MPSolverParameters ¶m) { + return interface_->Solve(param); +} + +int64 MPSolver::iterations() const { + return interface_->iterations(); +} + +int64 MPSolver::nodes() const { + return interface_->nodes(); +} + +// ---------- MPSolverInterface ---------- + +const int MPSolverInterface::kDummyVariableIndex = 0; + +MPSolverInterface::MPSolverInterface(MPSolver* const solver) + : solver_(solver), sync_status_(MODEL_SYNCHRONIZED), + result_status_(MPSolver::NOT_SOLVED), maximize_(false), + last_constraint_index_(0), last_variable_index_(0), + objective_value_(0.0) {} + +MPSolverInterface::~MPSolverInterface() {} + +void MPSolverInterface::WriteModelToPredefinedFiles() { + if (!FLAGS_solver_write_model.empty()) { + if (!solver_->CheckAllNamesValidity()) { + LOG(FATAL) << "Invalid name. Unable to write model to file"; + } + WriteModel(FLAGS_solver_write_model); + } + const string filename = solver_->write_model_filename(); + if (!filename.empty()) { + if (!solver_->CheckAllNamesValidity()) { + LOG(FATAL) << "Invalid name. Unable to write model to file"; + } + WriteModel(filename); + } +} + +// Extracts model stored in MPSolver +void MPSolverInterface::ExtractModel() { + + switch (sync_status_) { + case MUST_RELOAD: { + ExtractNewVariables(); + ExtractNewConstraints(); + ExtractObjective(); + + last_constraint_index_ = solver_->constraints_.size(); + last_variable_index_ = solver_->variables_.size(); + sync_status_ = MODEL_SYNCHRONIZED; + break; + } + case MODEL_SYNCHRONIZED: { + // Everything has already been extracted. + CHECK_EQ(last_constraint_index_, solver_->constraints_.size()); + CHECK_EQ(last_variable_index_, solver_->variables_.size()); + break; + } + case SOLUTION_SYNCHRONIZED: { + // Nothing has changed since last solve. + CHECK_EQ(last_constraint_index_, solver_->constraints_.size()); + CHECK_EQ(last_variable_index_, solver_->variables_.size()); + break; + } + } +} + +void MPSolverInterface::ResetExtractionInformation() { + sync_status_ = MUST_RELOAD; + last_constraint_index_ = 0; + last_variable_index_ = 0; + for (int j = 0; j < solver_->variables_.size(); ++j) { + MPVariable* const var = solver_->variables_[j]; + var->set_index(kNoIndex); + } + for (int i = 0; i < solver_->constraints_.size(); ++i) { + MPConstraint* const ct = solver_->constraints_[i]; + ct->set_index(kNoIndex); + } +} + +void MPSolverInterface::CheckSolutionIsSynchronized() const { + CHECK_EQ(SOLUTION_SYNCHRONIZED, sync_status_) << + "The model has been changed since the solution was last computed."; +} + +// Default version that can be overwritten by a solver-specific +// version to accomodate for the quirks of each solver. +void MPSolverInterface::CheckSolutionExists() const { + CHECK(result_status_ == MPSolver::OPTIMAL || + result_status_ == MPSolver::FEASIBLE) << + "No solution exists."; +} + +// Default version that can be overwritten by a solver-specific +// version to accomodate for the quirks of each solver. +void MPSolverInterface::CheckBestObjectiveBoundExists() const { + CHECK(result_status_ == MPSolver::OPTIMAL || + result_status_ == MPSolver::FEASIBLE) + << "No information is available for the best objective bound."; +} + +double MPSolverInterface::objective_value() const { + CheckSolutionIsSynchronized(); + CheckSolutionExists(); + return objective_value_; +} + +void MPSolverInterface::InvalidateSolutionSynchronization() { + if (sync_status_ == SOLUTION_SYNCHRONIZED) { + sync_status_ = MODEL_SYNCHRONIZED; + } +} + +// Solve with default parameters. +MPSolver::ResultStatus MPSolverInterface::Solve() { + MPSolverParameters default_param; + return Solve(default_param); +} + +void MPSolverInterface::SetCommonParameters(const MPSolverParameters& param) { + SetPresolveMode(param.GetIntegerParam(MPSolverParameters::PRESOLVE)); + // TODO(user): In the future, we could distinguish between the + // algorithm to solve the root LP and the algorithm to solve node + // LPs. Not sure if underlying solvers support it. + int value = param.GetIntegerParam(MPSolverParameters::LP_ALGORITHM); + if (value != MPSolverParameters::kDefaultIntegerParamValue) { + SetLpAlgorithm(value); + } +} + +void MPSolverInterface::SetMIPParameters(const MPSolverParameters& param) { + SetRelativeMipGap(param.GetDoubleParam(MPSolverParameters::RELATIVE_MIP_GAP)); +} + +void MPSolverInterface::SetUnsupportedDoubleParam( + MPSolverParameters::DoubleParam param) const { + LOG(WARNING) << "Trying to set an unsupported parameter: " << param << "."; +} +void MPSolverInterface::SetUnsupportedIntegerParam( + MPSolverParameters::IntegerParam param) const { + LOG(WARNING) << "Trying to set an unsupported parameter: " << param << "."; +} +void MPSolverInterface::SetDoubleParamToUnsupportedValue( + MPSolverParameters::DoubleParam param, int value) const { + LOG(WARNING) << "Trying to set a supported parameter: " << param + << " to an unsupported value: " << value; +} +void MPSolverInterface::SetIntegerParamToUnsupportedValue( + MPSolverParameters::IntegerParam param, double value) const { + LOG(WARNING) << "Trying to set a supported parameter: " << param + << " to an unsupported value: " << value; +} + +// ---------- MPSolverParameters ---------- + +const double MPSolverParameters::kDefaultRelativeMipGap = 1e-4; +const MPSolverParameters::PresolveValues MPSolverParameters::kDefaultPresolve = + MPSolverParameters::PRESOLVE_ON; + +const double MPSolverParameters::kDefaultDoubleParamValue = -1.0; +const int MPSolverParameters::kDefaultIntegerParamValue = -1; +const double MPSolverParameters::kUnknownDoubleParamValue = -2.0; +const int MPSolverParameters::kUnknownIntegerParamValue = -2; + +// The constructor sets all parameters to their default value. +MPSolverParameters::MPSolverParameters() + : relative_mip_gap_value_(kDefaultRelativeMipGap), + presolve_value_(kDefaultPresolve), + lp_algorithm_value_(kDefaultIntegerParamValue), + lp_algorithm_is_default_(true) {} + +void MPSolverParameters::SetDoubleParam(MPSolverParameters::DoubleParam param, + double value) { + switch (param) { + case RELATIVE_MIP_GAP: { + relative_mip_gap_value_ = value; + break; + } + default: { + LOG(ERROR) << "Trying to set an unknown parameter: " << param << "."; + } + } +} + +void MPSolverParameters::SetIntegerParam(MPSolverParameters::IntegerParam param, + int value) { + switch (param) { + case PRESOLVE: { + if (value != PRESOLVE_OFF && value != PRESOLVE_ON) { + LOG(ERROR) << "Trying to set a supported parameter: " << param + << " to an unknown value: " << value; + } + presolve_value_ = value; + break; + } + case LP_ALGORITHM: { + if (value != DUAL && value != PRIMAL && value != BARRIER) { + LOG(ERROR) << "Trying to set a supported parameter: " << param + << " to an unknown value: " << value; + } + lp_algorithm_value_ = value; + lp_algorithm_is_default_ = false; + break; + } + default: { + LOG(ERROR) << "Trying to set an unknown parameter: " << param << "."; + } + } +} + +void MPSolverParameters::ResetDoubleParam( + MPSolverParameters::DoubleParam param) { + switch (param) { + case RELATIVE_MIP_GAP: { + relative_mip_gap_value_ = kDefaultRelativeMipGap; + break; + } + default: { + LOG(ERROR) << "Trying to reset an unknown parameter: " << param << "."; + } + } +} + +void MPSolverParameters::ResetIntegerParam( + MPSolverParameters::IntegerParam param) { + switch (param) { + case PRESOLVE: { + presolve_value_ = kDefaultPresolve; + break; + } + case LP_ALGORITHM: { + lp_algorithm_is_default_ = true; + break; + } + default: { + LOG(ERROR) << "Trying to reset an unknown parameter: " << param << "."; + } + } +} + +void MPSolverParameters::Reset() { + ResetDoubleParam(RELATIVE_MIP_GAP); + ResetIntegerParam(PRESOLVE); + ResetIntegerParam(LP_ALGORITHM); +} + +double MPSolverParameters::GetDoubleParam( + MPSolverParameters::DoubleParam param) const { + switch (param) { + case RELATIVE_MIP_GAP: { + return relative_mip_gap_value_; + } + default: { + LOG(ERROR) << "Trying to get an unknown parameter: " << param << "."; + return kUnknownDoubleParamValue; + } + } +} + +int MPSolverParameters::GetIntegerParam( + MPSolverParameters::IntegerParam param) const { + switch (param) { + case PRESOLVE: { + return presolve_value_; + } + case LP_ALGORITHM: { + if (lp_algorithm_is_default_) return kDefaultIntegerParamValue; + return lp_algorithm_value_; + } + default: { + LOG(ERROR) << "Trying to get an unknown parameter: " << param << "."; + return kUnknownIntegerParamValue; + } + } +} + +} // namespace operations_research diff --git a/linear_solver/linear_solver.h b/linear_solver/linear_solver.h new file mode 100644 index 0000000000..56beb87ec9 --- /dev/null +++ b/linear_solver/linear_solver.h @@ -0,0 +1,742 @@ +// Copyright 2010 Google +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. +// +// (Laurent Perron). +// +// A C++ wrapper for the GNU linear programming / mixed integer programming kit, +// Coin LP Solver, and Coin Branch and Cut Solver. + +// ----------------------------------- +// What is Linear Programming ? +// +// In mathematics, linear programming (LP) is a technique for optimization of +// a linear objective function, subject to linear equality and linear +// inequality constraints. Informally, linear programming determines the way +// to achieve the best outcome (such as maximum profit or lowest cost) in a +// given mathematical model and given some list of requirements represented +// as linear equations. The linear programming problem was first shown to be +// solvable in polynomial time by Leonid Khachiyan in 1979, but a larger +// theoretical and practical breakthrough in the field came in 1984 when +// Narendra Karmarkar introduced a new interior point method for solving +// linear programming problems. +// +// ----------------------------------- +// What is Mixed Integer Programming ? +// +// If only some of the unknown variables are required to be integers, then +// the problem is called a mixed integer programming (MIP) problem. These are +// generally also NP-hard. +// +// ----------------------------------- +// Check Wikipedia for more detail: +// +// http://en.wikipedia.org/wiki/Linear_programming +// +// ----------------------------------- +// Example of a Linear Programming: +// +// mimimize: +// f1*x1+f2*x2+...fn*xn +// subject to: +// a1*x1+a2*x2+...an*xn >= k1 +// b1*x1+b2*x2+...bn*xn <= k2 +// c1*x1+c2*x2+...cn*xn = k3 +// ...... +// u1 <= x1 <= v1 +// u2 <= x2 <= v2 +// ..... +// +// As can be seen, Linear Programming has: +// 1) linear objective function +// 2) linear constraint +// +// Note: +// The objective function is linear and convex. +// The constraints form a convex space if feasible. +// +// ----------------------------------- + +#ifndef LINEAR_SOLVER_LINEAR_SOLVER_H_ +#define LINEAR_SOLVER_LINEAR_SOLVER_H_ + + +#include +#include +#include + +#include "base/commandlineflags.h" +#include "base/integral_types.h" +#include "base/logging.h" +#include "base/macros.h" +#include "base/scoped_ptr.h" +#include "base/timer.h" +#include "base/strutil.h" +#include "base/sparsetable.h" +#include "linear_solver/linear_solver.pb.h" + +namespace operations_research { + +class MPModelProto; +class MPSolverInterface; +class MPSolverParameters; + +// A class to express a variable that will appear in a constraint. +class MPVariable { + public: + const string& name() const { return name_; } + + void SetInteger(bool integer); + bool integer() const { return integer_; } + + double solution_value() const; + // Only available for continuous problems. + double reduced_cost() const; + + int index() const { return index_; } + + double lb() const { return lb_; } + double ub() const { return ub_; } + void SetLB(double lb) { SetBounds(lb, ub_); } + void SetUB(double ub) { SetBounds(lb_, ub); } + void SetBounds(double lb, double ub); + protected: + friend class MPSolver; + friend class MPSolverInterface; + friend class CBCInterface; + friend class CLPInterface; + friend class GLPKInterface; + + MPVariable(double lb, double ub, bool integer, const string& name, + MPSolverInterface* const interface) + : lb_(lb), ub_(ub), integer_(integer), name_(name), index_(-1), + solution_value_(0.0), reduced_cost_(0.0), interface_(interface) {} + + void set_index(int index) { index_ = index; } + void set_solution_value(double value) { solution_value_ = value; } + void set_reduced_cost(double reduced_cost) { reduced_cost_ = reduced_cost; } + private: + double lb_; + double ub_; + bool integer_; + const string name_; + int index_; + double solution_value_; + double reduced_cost_; + MPSolverInterface* const interface_; + DISALLOW_COPY_AND_ASSIGN(MPVariable); +}; + +// A class to express constraints for a linear programming problem. A +// constraint is represented as a linear equation/inequality. +class MPConstraint { + public: + const string& name() const { return name_; } + + // Clears all variables and coefficients. + void Clear(); + + // Add (var * coeff) to the current constraint. + void AddTerm(MPVariable* const var, double coeff); + // Add var to the current constraint. + void AddTerm(MPVariable* const var); + // Set the coefficient of the variable on the constraint. + void SetCoefficient(MPVariable* const var, double coeff); + + double lb() const { return lb_; } + double ub() const { return ub_; } + void SetLB(double lb) { SetBounds(lb, ub_); } + void SetUB(double ub) { SetBounds(lb_, ub); } + void SetBounds(double lb, double ub); + + // Only available for continuous problems. + double dual_value() const; + + int index() const { return index_; } + protected: + friend class MPSolver; + friend class MPSolverInterface; + friend class CBCInterface; + friend class CLPInterface; + friend class GLPKInterface; + + // Creates a constraint and updates the pointer to its MPSolverInterface. + MPConstraint(double lb, + double ub, + const string& name, + MPSolverInterface* const interface) + : lb_(lb), ub_(ub), name_(name), index_(-1), dual_value_(0.0), + interface_(interface) {} + + void set_index(int index) { index_ = index; } + void set_dual_value(double dual_value) { dual_value_ = dual_value; } + private: + // Returns true if the constraint contains variables that have not + // been extracted yet. + bool ContainsNewVariables(); + + // Mapping var -> coefficient. + hash_map coefficients_; + + // The lower bound for the linear constraint. + double lb_; + // The upper bound for the linear constraint. + double ub_; + // Name. + const string name_; + int index_; + double dual_value_; + MPSolverInterface* const interface_; + DISALLOW_COPY_AND_ASSIGN(MPConstraint); +}; + +// A class to express a linear objective function +class MPObjective { + public: + // Clears all variables and coefficients. + void Clear(); + + // Add (var * coeff) to the objective + void AddTerm(MPVariable* const var, double coeff); + // Add var to the objective + void AddTerm(MPVariable* const var); + // Set the coefficient of the variable in the objective + void SetCoefficient(MPVariable* const var, double coeff); + + // Add constant term to the objective. + void AddOffset(double value); + // Set constant term in the objective. + void SetOffset(double value); + + private: + friend class MPSolver; + friend class MPSolverInterface; + friend class CBCInterface; + friend class CLPInterface; + friend class GLPKInterface; + + // Creates an objective and updates the pointer to its parent 'MPSolver'. + explicit MPObjective(MPSolverInterface* const interface) + : offset_(0.0), interface_(interface) {} + + // Mapping var -> coefficient. + hash_map coefficients_; + // Constant term. + double offset_; + + MPSolverInterface* const interface_; + DISALLOW_COPY_AND_ASSIGN(MPObjective); +}; + +class MPSolver { + public: + + // The LP/MIP problem type. + enum OptimizationProblemType { +#if defined(USE_GLPK) + GLPK_LINEAR_PROGRAMMING, + GLPK_MIXED_INTEGER_PROGRAMMING, +#endif +#if defined(USE_CLP) + CLP_LINEAR_PROGRAMMING, +#endif +#if defined(USE_CBC) + CBC_MIXED_INTEGER_PROGRAMMING, +#endif + }; + + enum ResultStatus { + OPTIMAL, // optimal + FEASIBLE, // feasible, or stopped by limit. + INFEASIBLE, // proven infeasible + UNBOUNDED, // unbounded + ABNORMAL, // abnormal, i.e., error of some kind. + NOT_SOLVED // not been solved yet. + }; + + enum LoadStatus { + NO_ERROR, + DUPLICATE_VARIABLE_ID, + UNKNOWN_VARIABLE_ID + }; + + // Constructor that takes a name for the underlying solver. + MPSolver(const string& name, OptimizationProblemType problem_type); + virtual ~MPSolver(); + + // ----- Load model from protobuf ----- +#ifndef SWIG + // TODO(user): fix swig support. + LoadStatus Load(const MPModelProto& model); +#endif + + // ----- Init and Clear ----- + void Init() {} // To remove. + void Clear(); + + // ----- Variables ------ + // Returns the number of variables. + int NumVariables() const { return variables_.size(); } + const vector& variables() const { return variables_; } + + // Create a variable with the given bounds. + MPVariable* MakeVar(double lb, double ub, bool integer, const string& name); + MPVariable* MakeNumVar(double lb, double ub, const string& name); + MPVariable* MakeIntVar(double lb, double ub, const string& name); + MPVariable* MakeBoolVar(const string& name); + + void MakeVarArray(int nb, + double lb, + double ub, + bool integer, + const string& name, + vector* vars); + void MakeNumVarArray(int nb, + double lb, + double ub, + const string& name, + vector* vars); + void MakeIntVarArray(int nb, + double lb, + double ub, + const string& name, + vector* vars); + void MakeBoolVarArray(int nb, + const string& name, + vector* vars); + + // ----- Constraints ----- + // Returns the number of constraints. + int NumConstraints() const { return constraints_.size(); } + + // Returns a pointer to a newly created constraint for the linear programming + // problem. The MPSolver class assumes ownership of the constraint. + MPConstraint* MakeRowConstraint(double lb, double ub); + MPConstraint* MakeRowConstraint(); + MPConstraint* MakeRowConstraint(double lb, double ub, const string& name); + MPConstraint* MakeRowConstraint(const string& name); + + // ----- Objective ----- + + // Return the objective value of the best solution found so far. It + // is the optimal objective value if the problem has been solved to + // optimality. + double objective_value() const; + + // Returns the best objective bound. In case of minimization, it is + // a lower bound on the objective value of the optimal integer + // solution. Only available for discrete problems. + double best_objective_bound() const; + + // Clear objective. + void ClearObjective(); + // Add var * coeff to the objective function. + void AddObjectiveTerm(MPVariable* const var, double coeff); + // Add var to the objective function. + void AddObjectiveTerm(MPVariable* const var); + // Set the objective coefficient of a variable in the objective. + void SetObjectiveCoefficient(MPVariable* const var, double coeff); + // Add constant term to the objective. + void AddObjectiveOffset(double value); + // Set constant term in the objective. + void SetObjectiveOffset(double value); + + + // Sets the optimization direction (min/max). + void SetOptimizationDirection(bool maximize); + // Minimizing or maximizing? + bool Maximization(); + // Minimizing or maximizing? + bool Minimization(); + + // Set minimization mode. + void SetMinimization() { SetOptimizationDirection(false); } + // Set maximization mode. + void SetMaximization() { SetOptimizationDirection(true); } + + // ----- Solve ----- + // Solve the problem using default parameter values. + ResultStatus Solve(); + // Solve the problem using the parameter values specified. + ResultStatus Solve(const MPSolverParameters& param); + + // Reset extracted model to solve from scratch + void Reset(); + + // Misc. + static double infinity() { + return std::numeric_limits::infinity(); + } + + // SuppressOutput. + void SuppressOutput(); + + // Set the name of the file where the solver writes out the model + void set_write_model_filename(const string &filename) { + write_model_filename_ = filename; + } + + string write_model_filename() const { + return write_model_filename_; + } + + // Return true if filename ends in ".lp" + bool IsLPFormat(const string &filename) { + return HasSuffixString (filename, ".lp"); + } + + // Set Time limit in ms. (0 = no limit). + void set_time_limit(int64 time_limit) { + DCHECK_GE(time_limit, 0); + time_limit_ = time_limit; + } + + int64 time_limit() const { + return time_limit_; + } + + // wall_time() in ms since the creation of the solver. + int64 wall_time() const { + return timer_.GetInMs(); + } + + // Number of simplex iterations + int64 iterations() const; + // Number of branch-and-bound nodes. Only available for discrete problems. + int64 nodes() const; + + // Check validity of a name. + bool CheckNameValidity(const string& name); + // Check validity of all variables and constraints names. + bool CheckAllNamesValidity(); + + // return a string describing the engine used. + string SolverVersion() const; + + friend class GLPKInterface; + friend class CLPInterface; + friend class CBCInterface; + friend class MPSolverInterface; + + private: + // Compute the size of the constraint with the largest number of + // coefficients with index in [min_constraint_index, + // max_constraint_index) + int ComputeMaxConstraintSize(int min_constraint_index, + int max_constraint_index) const; + + // The name of the linear programming problem. + const string name_; + + // The solver interface. + scoped_ptr interface_; + + // vector of problem variables. + vector variables_; + hash_set variables_names_; + + // The list of constraints for the problem. + vector constraints_; + hash_set constraints_names_; + + // The linear objective function + MPObjective linear_objective_; + + + // Time limit in ms. + int64 time_limit_; + + // Name of the file where the solver writes out the model when Solve + // is called. If empty, no file is written. + string write_model_filename_; + + WallTimer timer_; + + DISALLOW_COPY_AND_ASSIGN(MPSolver); +}; + +// This class stores parameter settings for LP and MIP solvers. +// How to add a new parameter: +// - Add the new Foo parameter in the DoubleParam or IntegerParam enum. +// - If it is a categorical param, add a FooValues enum. +// - Decide if the wrapper should define a default value for it: yes +// if it controls the properties of the solution (example: +// tolerances) or if it consistently improves performance, no +// otherwise. If yes, define kDefaultFoo. +// - Add a foo_value_ member and, if no default value is defined, a +// foo_is_default_ member. +// - Add code to handle Foo in Set...Param, Reset...Param, +// Get...Param, Reset and the constructor. +// - In class MPSolverInterface, add a virtual method SetFoo and +// implement it for each solver. +// - Add a test in linear_solver_test.cc. +class MPSolverParameters { + public: + // Enumeration of parameters that take continuous values. + enum DoubleParam { + RELATIVE_MIP_GAP = 0, // Limit for relative MIP gap. + }; + + // Enumeration of parameters that take integer or categorical values. + enum IntegerParam { + PRESOLVE = 1000, // Presolve mode. + LP_ALGORITHM = 1001 // Algorithm to solve linear programs. + }; + + // For each categorical parameter, enumeration of possible values. + enum PresolveValues { + PRESOLVE_OFF = 0, // Presolve is off. + PRESOLVE_ON = 1 // Presolve is on. + }; + + enum LpAlgorithmValues { + DUAL = 10, // Dual simplex. + PRIMAL = 11, // Primal simplex. + BARRIER = 12 // Barrier algorithm. + }; + + // Values to indicate that a parameter is set to the solver's + // default value. + static const double kDefaultDoubleParamValue; + static const int kDefaultIntegerParamValue; + // Values to indicate that a parameter is unknown. + static const double kUnknownDoubleParamValue; + static const int kUnknownIntegerParamValue; + + // Default values for parameters. Only parameters that define the + // properties of the solution returned need to have a default value + // (that is the same for all solvers). You can also define a default + // value for performance parameters when you are confident it is a + // good choice (example: always turn presolve on). + static const double kDefaultRelativeMipGap; + static const PresolveValues kDefaultPresolve; + + // The constructor sets all parameters to their default value. + MPSolverParameters(); + + // Set parameter to a specific value. + void SetDoubleParam(MPSolverParameters::DoubleParam param, double value); + void SetIntegerParam(MPSolverParameters::IntegerParam param, int value); + // Reset parameter to the default value. + void ResetDoubleParam(MPSolverParameters::DoubleParam param); + void ResetIntegerParam(MPSolverParameters::IntegerParam param); + // Set all parameters to their default value. + void Reset(); + // Get parameter's value. + double GetDoubleParam(MPSolverParameters::DoubleParam param) const; + int GetIntegerParam(MPSolverParameters::IntegerParam param) const; + + private: + // Parameter value for each parameter. + double relative_mip_gap_value_; + int presolve_value_; + int lp_algorithm_value_; + // Boolean value indicating whether each parameter is set to the + // solver's default value. Only parameters for which the wrapper + // does not define a default value need such an indicator. + bool lp_algorithm_is_default_; + + DISALLOW_COPY_AND_ASSIGN(MPSolverParameters); +}; + +// This class serves as a proxy to open sources linear solver. +class MPSolverInterface { + public: + enum SynchronizationStatus { + // The underlying solver (CLP, GLPK, ...) and MPSolver are not in + // sync for the model nor for the solution. + MUST_RELOAD, + // The underlying solver and MPSolver are in sync for the model + // but not for the solution: the model has changed since the + // solution was computed last. + MODEL_SYNCHRONIZED, + // The underlying solver and MPSolver are in sync for the model and + // the solution. + SOLUTION_SYNCHRONIZED + }; + + // When the underlying solver does not provide the number of simplex + // iterations. + static const int64 kUnknownNumberOfIterations; + // When the underlying solver does not provide the number of simplex + // nodes. + static const int64 kUnknownNumberOfNodes; + // When the index of a variable or constraint has not been assigned yet. + static const int kNoIndex; + + // Constructor that takes a name for the underlying glpk solver. + explicit MPSolverInterface(MPSolver* const solver); + virtual ~MPSolverInterface(); + + // ----- Solve ----- + // Solves problem with default parameter values. Returns true if the + // solution is optimal. Calls WriteModelToPredefinedFiles as a + // temporary solution to allow the user to write the model to a + // file. + MPSolver::ResultStatus Solve(); + // Same as Solve(), except it uses the specified parameter values. + virtual MPSolver::ResultStatus Solve(const MPSolverParameters& param) = 0; + + // ----- Model modifications and extraction ----- + // Resets extracted model + virtual void Reset() = 0; + + // Sets the optimization direction (min/max). + virtual void SetOptimizationDirection(bool minimize) = 0; + + // Modify bounds of an extracted variable. + virtual void SetVariableBounds(int index, double lb, double ub) = 0; + + // Modify integrality of an extracted variable. + virtual void SetVariableInteger(int index, bool integer) = 0; + + // Modify bounds of an extracted variable. + virtual void SetConstraintBounds(int index, double lb, double ub) = 0; + + // Add a constraint. + virtual void AddRowConstraint(MPConstraint* const ct) = 0; + + // Add a variable. + virtual void AddVariable(MPVariable* const var) = 0; + + // Change a coefficient in a constraint. + virtual void SetCoefficient(MPConstraint* const constraint, + MPVariable* const variable, + double coefficient) = 0; + + // Clear a constraint from all its terms. + virtual void ClearConstraint(MPConstraint* const constraint) = 0; + + // Change a coefficient in the linear objective. + virtual void SetObjectiveCoefficient(MPVariable* const variable, + double coefficient) = 0; + + // Change the constant term in the linear objective. + virtual void SetObjectiveOffset(double value) = 0; + + // Clear the objective from all its terms. + virtual void ClearObjective() = 0; + + // ------ Query statistics on the solution and the solve ------ + // Number of simplex iterations + virtual int64 iterations() const = 0; + // Number of branch-and-bound nodes + virtual int64 nodes() const = 0; + // Best objective bound. Only available for discrete problems. + virtual double best_objective_bound() const = 0; + // Objective value of the best solution found so far. + double objective_value() const; + + // Checks whether the solution is synchronized with the model, + // i.e. whether the model has changed since the solution was + // computed last. + void CheckSolutionIsSynchronized() const; + // Checks whether a feasible solution exists. + virtual void CheckSolutionExists() const; + // Checks whether information on the best objective bound exists. + virtual void CheckBestObjectiveBoundExists() const; + + // ----- Misc ----- + // Write model to file. + virtual void WriteModel(const string& filename) = 0; + + // SuppressOutput. + virtual void SuppressOutput() = 0; + + // Query problem type. For simplicity, the distinction between + // continuous and discrete is based on the declaration of the user + // when the solver is created (example: GLPK_LINEAR_PROGRAMMING + // vs. GLPK_MIXED_INTEGER_PROGRAMMING), not on the actual content of + // the model. + // Returns true if the problem is continuous. + virtual bool IsContinuous() const = 0; + // Returns true if the problem is continuous and linear. + virtual bool IsLP() const = 0; + // Returns true if the problem is discrete and linear. + virtual bool IsMIP() const = 0; + + int last_variable_index() const { + return last_variable_index_; + } + // Returns the result status of the last solve. + MPSolver::ResultStatus result_status() const { + CheckSolutionIsSynchronized(); + return result_status_; + } + // Returns a string describing the solver. + virtual string SolverVersion() const = 0; + + friend class MPSolver; + protected: + MPSolver* const solver_; + // Indicates whether the model and the solution are synchronized. + SynchronizationStatus sync_status_; + // Indicates whether the solve has reached optimality, + // infeasibility, a limit, etc. + MPSolver::ResultStatus result_status_; + bool maximize_; + + // Index of last constraint extracted + int last_constraint_index_; + // Index of last variable extracted + int last_variable_index_; + + // The value of the objective function. + double objective_value_; + + // Index of dummy variable created for empty constraints or the + // objective offset. + static const int kDummyVariableIndex; + + // Writes out the model to a file specified by the + // --solver_write_model command line argument or + // MPSolver::set_write_model_filename. + // The file is written by each solver interface (CBC, CLP, GLPK) and + // each behaves a little differently. + // If filename ends in ".lp", then the file is written in the + // LP format (except for the CLP solver that does not support the LP + // format). In all other cases it is written in the MPS format. + void WriteModelToPredefinedFiles(); + + // Extracts model stored in MPSolver + void ExtractModel(); + virtual void ExtractNewVariables() = 0; + virtual void ExtractNewConstraints() = 0; + virtual void ExtractObjective() = 0; + void ResetExtractionInformation(); + // Change synchronization status from SOLUTION_SYNCHRONIZED to + // MODEL_SYNCHRONIZED. To be used for model changes. + void InvalidateSolutionSynchronization(); + + // Set parameters common to LP and MIP in the underlying solver. + void SetCommonParameters(const MPSolverParameters& param); + // Set MIP specific parameters in the underlying solver. + void SetMIPParameters(const MPSolverParameters& param); + // Set all parameters in the underlying solver. + virtual void SetParameters(const MPSolverParameters& param) = 0; + // Set an unsupported parameter. + void SetUnsupportedDoubleParam(MPSolverParameters::DoubleParam param) const; + void SetUnsupportedIntegerParam(MPSolverParameters::IntegerParam param) const; + // Set a supported parameter to an unsupported value. + void SetDoubleParamToUnsupportedValue(MPSolverParameters::DoubleParam param, + int value) const; + void SetIntegerParamToUnsupportedValue(MPSolverParameters::IntegerParam param, + double value) const; + // Set each parameter in the underlying solver. + virtual void SetRelativeMipGap(double value) = 0; + virtual void SetPresolveMode(int value) = 0; + virtual void SetLpAlgorithm(int value) = 0; +}; + +} // namespace operations_research + +#endif // LINEAR_SOLVER_LINEAR_SOLVER_H_ diff --git a/linear_solver/linear_solver.proto b/linear_solver/linear_solver.proto new file mode 100644 index 0000000000..5191b3d84b --- /dev/null +++ b/linear_solver/linear_solver.proto @@ -0,0 +1,59 @@ +// Copyright 2010 Google +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +// Linear solver Stubby services + +syntax = "proto2"; + +package operations_research; + +message MPVariableProto { + required string id = 1; // key + required double lb = 2; + required double ub = 3; + required bool integer = 4; +} + +message MPQuadraticTermProto { + required string variable1_id = 1; + required string variable2_id = 2; + required double coefficient = 3; +} + +message MPVariableSampleProto { + required string variable_id = 1; + repeated double sample_values = 2; +} + +message MPTermProto { + required string variable_id = 1; + required double coefficient = 2; +} + +message MPConstraintProto { + required double lb = 1; + required double ub = 2; + repeated MPTermProto terms = 3; + optional string id = 4; +} + +message MPModelProto { + repeated MPVariableProto variables = 1; + required bool maximize = 2; + repeated MPTermProto objective_terms = 3; + repeated MPQuadraticTermProto quadratic_objective_terms = 4; + repeated MPQuadraticTermProto lower_triangular_objective_terms = 5; + repeated MPVariableSampleProto covariance_objective_samples = 6; + repeated MPConstraintProto constraints = 7; + optional string name = 8; +} diff --git a/linear_solver/linear_solver.swig b/linear_solver/linear_solver.swig new file mode 100644 index 0000000000..338f957e21 --- /dev/null +++ b/linear_solver/linear_solver.swig @@ -0,0 +1,493 @@ +// Copyright (C) 2009 and onwards Google +// lperron@google.com (Laurent Perron) + +%include "base/base.swig" + +// Include the file we want to wrap a first time. +%{ +#include "linear_solver/linear_solver.h" +%} + +#ifdef SWIGPYTHON + +// Define the renaming of methods. +%rename (BoolVar) MakeBoolVar; +%rename (CovarianceMatrix) MPCovarianceMatrix; +%rename (IntVar) MakeIntVar; +%rename (LowerTriangularMatrix) MPLowerTriangularMatrix; +%rename (NumVar) MakeNumVar; +%rename (QuadraticMatrix) MPQuadraticMatrix; +%rename (Constraint) MPConstraint; +%rename (Constraint) MakeRowConstraint; +%rename (Solver) MPSolver; +%rename (Variable) MPVariable; +%ignore MakeVarArray; +%ignore MakeNumVarArray; +%ignore MakeIntVarArray; +%ignore MakeBoolVarArray; + +namespace operations_research { +%pythoncode { + +class LinearExpr(object): + def DoVisit(self, coeffs, multiplier): + raise NotImplementedError + + def Visit(self, coeffs): + return self.DoVisit(coeffs, 1.0) + + def __add__(self, expr): + if isinstance(expr, (int, long, float)): + return SumCst(self, expr) + else: + return Sum(self, expr) + + def __radd__(self, cst): + if isinstance(cst, (int, long, float)): + return SumCst(self, cst) + else: + raise TypeError + + def __sub__(self, expr): + if isinstance(expr, (int, long, float)): + return SumCst(self, -expr) + else: + return Sum(self, ProductCst(expr, -1)) + + def __rsub__(self, cst): + if isinstance(cst, (int, long, float)): + return SumCst(ProductCst(self, -1), cst) + else: + raise TypeError + + def __mul__(self, cst): + if isinstance(cst, (int, long, float)): + return ProductCst(self, cst) + else: + raise TypeError + + def __rmul__(self, cst): + if isinstance(cst, (int, long, float)): + return ProductCst(self, cst) + else: + raise TypeError + + def __div__(self, cst): + if isinstance(cst, (int, long, float)): + if cst == 0.0: + raise ZeroDivisionError + else: + return ProductCst(self, 1.0 / cst) + else: + raise TypeError + + def __truediv__(self, cst): + if isinstance(cst, (int, long, float)): + if cst == 0.0: + raise ZeroDivisionError + else: + return ProductCst(self, 1.0 / cst) + else: + raise TypeError + + def __neg__(self): + return ProductCst(self, -1) + + def __eq__(self, arg): + if isinstance(arg, (int, long, float)): + return LinearConstraint(self, arg, arg) + else: + return LinearConstraint(Sum(self, ProductCst(arg, -1)), 0.0, 0.0) + + def __ge__(self, arg): + if isinstance(arg, (int, long, float)): + return LinearConstraint(self, arg, 1e308) + else: + return LinearConstraint(Sum(self, ProductCst(arg, -1)), 0.0, 1e308) + + def __le__(self, arg): + if isinstance(arg, (int, long, float)): + return LinearConstraint(self, -1e308, arg) + else: + return LinearConstraint(Sum(self, ProductCst(arg, -1)), -1e308, 0.0) + + +class ProductCst(LinearExpr): + def __init__(self, expr, coef): + self.__expr = expr + self.__coef = coef + + def __str__(self): + if (self.__coef == -1): + return '-' + str(self.__expr) + else: + return '(' + str(self.__coef) + ' * ' + str(self.__expr) + ')' + + def DoVisit(self, coeffs, multiplier): + current_multiplier = multiplier * self.__coef + if current_multiplier: + return self.__expr.DoVisit(coeffs, current_multiplier) + return 0.0 + +class Sum(LinearExpr): + def __init__(self, left, right): + self.__left = left + self.__right = right + + def __str__(self): + return '(' + str(self.__left) + ' + ' + str(self.__right) + ')' + + def DoVisit(self, coeffs, multiplier): + constant = self.__left.DoVisit(coeffs, multiplier) + constant += self.__right.DoVisit(coeffs, multiplier) + return constant + +class SumArray(LinearExpr): + def __init__(self, array): + self.__array = array + + def __str__(self): + return 'Sum(' + str(self.__array) + ')' + + def DoVisit(self, coeffs, multiplier): + constant = 0.0 + for t in self.__array: + if isinstance(t, (int, long, float)): + constant += t * multiplier + else: + constant += t.DoVisit(coeffs, multiplier) + return constant + + +class SumCst(LinearExpr): + def __init__(self, expr, cst): + self.__expr = expr + self.__cst = cst + + def __str__(self): + return '(' + str(self.__expr) + ' + ' + str(self.__cst) + ')' + + def DoVisit(self, coeffs, multiplier): + constant = self.__expr.DoVisit(coeffs, multiplier) + return constant + self.__cst * multiplier + +class LinearConstraint(object): + def __init__(self, expr, lb, ub): + self.__expr = expr + self.__lb = lb + self.__ub = ub + + def __str__(self): + if self.__lb > -1e308 and self.__ub < 1e308: + if self.__lb == self.__ub: + return str(self.__expr) + ' == ' + str(self.__lb) + else: + return (str(self.__lb) + ' <= ' + str(self.__expr) + + " <= " + str(self.__ub)) + elif self.__lb > -1e308: + return str(self.__expr) + ' >= ' + str(self.__lb) + elif self.__ub < 1e308: + return str(self.__expr) + ' <= ' + str(self.__ub) + else: + return 'true inequality' + + def Extract(self, solver): + coeffs = {} + constant = self.__expr.Visit(coeffs) + lb = -solver.infinity() + ub = solver.infinity() + if self.__lb > -1e308: + lb = self.__lb - constant + if self.__ub < 1e308: + ub = self.__ub - constant + + constraint = solver.RowConstraint(lb, ub) + for v, c, in coeffs.iteritems(): + constraint.AddTerm(v, float(c)) + return constraint +} + +%extend MPVariable { + string __str__() { + return self->name(); + } + string __repr__() { + return self->name(); + } + %pythoncode { + def __add__(self, expr): + if isinstance(expr, (int, long, float)): + return SumCst(self, expr) + else: + return Sum(self, expr) + + def __radd__(self, cst): + if isinstance(cst, (int, long, float)): + return SumCst(self, cst) + else: + raise TypeError + + def __sub__(self, expr): + if isinstance(expr, (int, long, float)): + return SumCst(self, -expr) + else: + return Sum(self, ProductCst(expr, -1)) + + def __rsub__(self, cst): + if isinstance(cst, (int, long, float)): + return SumCst(ProductCst(self, -1), cst) + else: + raise TypeError + + def __mul__(self, cst): + if isinstance(cst, (int, long, float)): + return ProductCst(self, cst) + else: + raise TypeError + + def __rmul__(self, cst): + if isinstance(cst, (int, long, float)): + return ProductCst(self, cst) + else: + raise TypeError + + def __div__(self, cst): + if isinstance(cst, (int, long, float)): + if cst == 0.0: + raise ZeroDivisionError + else: + return ProductCst(self, 1.0 / cst) + else: + raise TypeError + + def __truediv__(self, cst): + if isinstance(cst, (int, long, float)): + if cst == 0.0: + raise ZeroDivisionError + else: + return ProductCst(self, 1.0 / cst) + else: + raise TypeError + + def __neg__(self): + return ProductCst(self, -1) + + def __eq__(self, arg): + if isinstance(arg, (int, long, float)): + return LinearConstraint(self, arg, arg) + else: + return LinearConstraint(Sum(self, ProductCst(arg, -1)), 0.0, 0.0) + + def __ge__(self, arg): + if isinstance(arg, (int, long, float)): + return LinearConstraint(self, arg, 1e308) + else: + return LinearConstraint(Sum(self, ProductCst(arg, -1)), 0.0, 1e308) + + def __le__(self, arg): + if isinstance(arg, (int, long, float)): + return LinearConstraint(self, -1e308, arg) + else: + return LinearConstraint(Sum(self, ProductCst(arg, -1)), -1e308, 0.0) + + def Visit(self, coeffs): + return self.DoVisit(coeffs, 1.0) + + def DoVisit(self, coeffs, multiplier): + if self in coeffs: + coeffs[self] += multiplier + else: + coeffs[self] = multiplier + return 0.0 + } +} + +%extend MPSolver { + %pythoncode { + def Add(self, constraint): + return constraint.Extract(self) + + def Sum(self, expr_array): + result = SumArray(expr_array) + return result + + # Compatibility + def RowConstraint(self, *args): + return self.Constraint(*args) + + def Minimize(self, expr): + self.ClearObjective() + self.SetMinimization() + coeffs = {} + offset = expr.Visit(coeffs) + self.AddObjectiveOffset(offset) + for v, c, in coeffs.iteritems(): + self.AddObjectiveTerm(v, float(c)) + + def Maximize(self, expr): + self.ClearObjective() + self.SetMaximization() + coeffs = {} + offset = expr.Visit(coeffs) + self.AddObjectiveOffset(offset) + for v, c, in coeffs.iteritems(): + self.AddObjectiveTerm(v, float(c)) + } +} + +%apply SWIGTYPE *DISOWN { MPBaseQuadraticMatrix * }; + +} // namespace operations_research + +#endif + +namespace operations_research { + +// Hack to allow us to relinquish ownership of matrix objects after calling +// a method that claims ownership, such as MPSolver::SetQuadraticObjective. +// TODO(user): Find a nicer way to do this. (The SWIG docs are really bad!) +%typemap(javacode) MPBaseQuadraticMatrix %{ + /** + * Asserts that this object is owned by C++, so Java should never delete it. + */ + public void disown() { + swigCMemOwn = false; + } +%} +} // namespace operations_research + +#ifdef SWIGJAVA +%module(directors="1") operations_research; + +namespace operations_research { +// Rename rules on MPVariable. +%rename (reducedCost) MPVariable::reduced_cost; +%rename (setBounds) MPVariable::SetBounds; +%rename (setInteger) MPVariable::SetInteger; +%rename (setLb) MPVariable::SetLB; +%rename (setUb) MPVariable::SetUB; +%rename (solutionValue) MPVariable::solution_value; + +// Rename rules on MPConstraint. +%rename (addTerm) MPConstraint::AddTerm; +%rename (dualValue) MPConstraint::dual_value; +%rename (setBounds) MPConstraint::SetBounds; +%rename (setCoefficient) MPConstraint::SetCoefficient; +%rename (setLb) MPConstraint::SetLB; +%rename (setUb) MPConstraint::SetUB; + +// Rename rules on MPObjective. +%rename (addOffset) MPObjective::AddOffset; +%rename (addTerm) MPObjective::AddObjective; +%rename (clear) MPObjective::Clear; +%rename (setCoefficient) MPObjective::SetCoefficient; +%rename (setOffset) MPObjective::SetOffset; + +// Rename rules on MPSolverParameters. +%rename (getDoubleParam) MPSolverParameters::GetDoubleParam; +%rename (getIntegerParam) MPSolverParameters::GetIntegerParam; +%rename (reset) MPSolverParameters::Reset; +%rename (resetDoubleParam) MPSolverParameters::ResetDoubleParam; +%rename (resetIntegerParam) MPSolverParameters::ResetIntegerParam; +%rename (setDoubleParam) MPSolverParameters::SetDoubleParam; +%rename (setIntegerParam) MPSolverParameters::SetIntegerParam; + +// Rename rules on MPSolver. +%rename (addObjectiveOffset) MPSolver::AddObjectiveOffset; +%rename (addObjectiveTerm) MPSolver::AddObjectiveTerm; +%rename (bestObjectiveBound) MPSolver::best_objective_bound; +%rename (checkAllNamesValidity) MPSolver::CheckAllNamesValidity; +%rename (checkNameValidity) MPSolver::CheckNameValidity; +%rename (clear) MPSolver::Clear; +%rename (clearObjective) MPSolver::ClearObjective; +%rename (init) MPSolver::Init; +%rename (isLpFormat) MPSolver::IsLPFormat; +%rename (makeBoolVar) MPSolver::MakeBoolVar; +%rename (makeIntVar) MPSolver::MakeIntVar; +%rename (makeNumVar) MPSolver::MakeNumVar; +%rename (makeConstraint) MPSolver::MakeRowConstraint; +%rename (makeVar) MPSolver::MakeVar; +%rename (minimization) MPSolver::Minimization; +%rename (numConstraints) MPSolver::NumConstraints; +%rename (numVariables) MPSolver::NumVariables; +%rename (objectiveValue) MPSolver::objective_value; +%rename (reset) MPSolver::Reset; +%rename (setMaximization) MPSolver::SetMaximization; +%rename (setMinimization) MPSolver::SetMinimization; +%rename (setObjectiveCoefficient) MPSolver::SetObjectiveCoefficient; +%rename (setObjectiveOffset) MPSolver::SetObjectiveOffset; +%rename (setOptimizationDirection) MPSolver::SetOptimizationDirection; +%rename (setTimeLimit) MPSolver::set_time_limit; +%rename (setWriteModelFilename) MPSolver::set_write_model_filename; +%rename (solve) MPSolver::Solve; +%rename (solverVersion) MPSolver::SolverVersion; +%rename (suppressOutput) MPSolver::SuppressOutput; +%rename (timeLimit) MPSolver::time_limit; +%rename (wallTime) MPSolver::wall_time; +%rename (writeModelFilename) MPSolver::write_model_filename; +// Ignore code on MPSolver. +%ignore MPSolver::MakeVarArray; +%ignore MPSolver::MakeNumVarArray; +%ignore MPSolver::MakeIntVarArray; +%ignore MPSolver::MakeBoolVarArray; +// Add java code on MPSolver. +%typemap(javacode) MPSolver %{ + public MPVariable[] makeVarArray(int count, + double lb, + double ub, + boolean integer) { + MPVariable[] array = new MPVariable[count]; + for (int i = 0; i < count; ++i) { + array[i] = makeVar(lb, ub, integer, ""); + } + return array; + } + + public MPVariable[] makeVarArray(int count, + double lb, + double ub, + boolean integer, + String var_name) { + MPVariable[] array = new MPVariable[count]; + for (int i = 0; i < count; ++i) { + array[i] = makeVar(lb, ub, integer, var_name + i); + } + return array; + } + + public MPVariable[] makeNumVarArray(int count, double lb, double ub) { + return makeVarArray(count, lb, ub, false); + } + + public MPVariable[] makeNumVarArray(int count, + double lb, + double ub, + String var_name) { + return makeVarArray(count, lb, ub, false, var_name); + } + + public MPVariable[] makeIntVarArray(int count, double lb, double ub) { + return makeVarArray(count, lb, ub, true); + } + + public MPVariable[] makeIntVarArray(int count, + double lb, + double ub, + String var_name) { + return makeVarArray(count, lb, ub, true, var_name); + } + + public MPVariable[] makeBoolVarArray(int count) { + return makeVarArray(count, 0.0, 1.0, true); + } + + public MPVariable[] makeBoolVarArray(int count, String var_name) { + return makeVarArray(count, 0.0, 1.0, true, var_name); + } +%} +} // namespace operations_research + +#endif // SWIGJAVA + +// Wrap linear_solver includes +%include "linear_solver/linear_solver.h" diff --git a/python/integer_solver_example.py b/python/integer_solver_example.py new file mode 100644 index 0000000000..37745bb20c --- /dev/null +++ b/python/integer_solver_example.py @@ -0,0 +1,113 @@ +# Copyright 2010 Google +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""pywraplp example file.""" + + + +from google.apputils import app +import gflags + +from linear_solver import pywraplp + +FLAGS = gflags.FLAGS + + +class PyWrapMIPExamples(object): + """Class that contains a collection of MIP examples.""" + + def RunFirstMIPExample(self, mode): + """Minimal MIP Example.""" + solver = pywraplp.Solver('RunFirstMIPExample', mode) + infinity = solver.infinity() + x1 = solver.IntVar(0.0, infinity, 'x1') + x2 = solver.IntVar(0.0, infinity, 'x2') + + solver.AddObjectiveTerm(x1, 1) + solver.AddObjectiveTerm(x2, 2) + + ct = solver.Constraint(17, infinity) + ct.AddTerm(x1, 3) + ct.AddTerm(x2, 2) + + # Check the solution. + solver.Solve() + print solver.objective_value() + + def RunAllFirstMIPExample(self): + self.RunFirstMIPExample(pywraplp.Solver.GLPK_MIXED_INTEGER_PROGRAMMING) + self.RunFirstMIPExample(pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING) + + def RunFirstMIPExampleNewAPI(self, mode): + """Minimal MIP Example with New API.""" + solver = pywraplp.Solver('RunFirstMIPExample', mode) + infinity = solver.infinity() + x1 = solver.IntVar(0.0, infinity, 'x1') + x2 = solver.IntVar(0.0, infinity, 'x2') + + solver.Minimize(x1 + 2 * x2) + solver.Add(3 * x1 + 2 * x2 >= 17) + + # Check the solution. + solver.Solve() + print solver.objective_value() + + def RunAllFirstMIPExampleNewAPI(self): + self.RunFirstMIPExampleNewAPI( + pywraplp.Solver.GLPK_MIXED_INTEGER_PROGRAMMING) + self.RunFirstMIPExampleNewAPI(pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING) + + def RunSuccessiveObjectives(self, mode): + """Example with succesive objectives.""" + solver = pywraplp.Solver('RunSuccessiveObjectives', mode) + x1 = solver.NumVar(0, 10, 'var1') + x2 = solver.NumVar(0, 10, 'var2') + solver.Add(x1 + 2*x2 <= 10) + solver.Maximize(x1) + + # Check the solution + solver.Solve() + print x1.solution_value() + print x2.solution_value() + + solver.Maximize(x2) + + # Check the solution + solver.Solve() + print x1.solution_value() + print x2.solution_value() + + solver.Minimize(-x1) + + # Check the solution + solver.Solve() + print x1.solution_value() + print x2.solution_value() + + def RunAllSuccessiveObjectives(self): + self.RunSuccessiveObjectives(pywraplp.Solver.GLPK_MIXED_INTEGER_PROGRAMMING) + self.RunSuccessiveObjectives(pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING) + + def RunAllExamples(self): + self.RunAllFirstMIPExample() + self.RunAllFirstMIPExampleNewAPI() + self.RunAllSuccessiveObjectives() + + +def main(unused_argv): + mip_example = PyWrapMIPExamples() + mip_example.RunAllExamples() + + +if __name__ == '__main__': + app.run() diff --git a/python/linear_solver_example.py b/python/linear_solver_example.py new file mode 100644 index 0000000000..f345781559 --- /dev/null +++ b/python/linear_solver_example.py @@ -0,0 +1,147 @@ +# Copyright 2010 Google +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""pywraplp example file.""" + + + +from google.apputils import app +import gflags + +from linear_solver import pywraplp + +FLAGS = gflags.FLAGS + + +class PyWrapLPExamples(object): + """Class that contains a collection of LP examples.""" + + def RunFirstLinearExample(self, mode): + """Minimal Linear Example.""" + solver = pywraplp.Solver('RunFirstLinearExample', mode) + infinity = solver.infinity() + x1 = solver.NumVar(0.0, infinity, 'x1') + x2 = solver.NumVar(0.0, infinity, 'x2') + x3 = solver.NumVar(0.0, infinity, 'x3') + print solver.NumVariables() + + solver.AddObjectiveTerm(x1, 10) + solver.AddObjectiveTerm(x2, 6) + solver.AddObjectiveTerm(x3, 4) + solver.SetMaximization() + + c0 = solver.Constraint(-infinity, 100.0) + c0.AddTerm(x1, 1) + c0.AddTerm(x2, 1) + c0.AddTerm(x3, 1) + c1 = solver.Constraint(-infinity, 600.0) + c1.AddTerm(x1, 10) + c1.AddTerm(x2, 4) + c1.AddTerm(x3, 5) + c2 = solver.Constraint(-infinity, 300.0) + c2.AddTerm(x1, 2) + c2.AddTerm(x2, 2) + c2.AddTerm(x3, 6) + print solver.NumConstraints() + + # The problem has an optimal solution. + solver.Solve() + + print solver.objective_value() + print 'x1 = ', x1.solution_value(), ', reduced_cost = ', x1.reduced_cost() + print 'x2 = ', x2.solution_value(), ', reduced_cost = ', x2.reduced_cost() + print 'x3 = ', x3.solution_value(), ', reduced_cost = ', x3.reduced_cost() + print 'c0 dual value = ', c0.dual_value() + print 'c1 dual value = ', c1.dual_value() + print 'c2 dual value = ', c2.dual_value() + + def RunAllFirstLinearExample(self): + self.RunFirstLinearExample(pywraplp.Solver.GLPK_LINEAR_PROGRAMMING) + self.RunFirstLinearExample(pywraplp.Solver.CLP_LINEAR_PROGRAMMING) + + def RunFirstLinearExampleNewAPI(self, mode): + """Minimal LP Example with New API.""" + solver = pywraplp.Solver('RunFirstLinearExample', mode) + infinity = solver.infinity() + x1 = solver.NumVar(0.0, infinity, 'x1') + x2 = solver.NumVar(0.0, infinity, 'x2') + x3 = solver.NumVar(0.0, infinity, 'x3') + print solver.NumVariables() + + solver.Maximize(10 * x1 + 6 * x2 + 4 * x3) + + c0 = solver.Add(x1 + x2 + x3 <= 100.0) + c1 = solver.Add(10 * x1 + 4 * x2 + 5 * x3 <= 600) + c2 = solver.Add(2 * x1 + 2 * x2 + 6 * x3 <= 300) + + print solver.NumConstraints() + + # The problem has an optimal solution. + solver.Solve() + + print solver.objective_value() + print 'x1 = ', x1.solution_value(), ', reduced_cost = ', x1.reduced_cost() + print 'x2 = ', x2.solution_value(), ', reduced_cost = ', x2.reduced_cost() + print 'x3 = ', x3.solution_value(), ', reduced_cost = ', x3.reduced_cost() + print 'c0 dual value = ', c0.dual_value() + print 'c1 dual value = ', c1.dual_value() + print 'c2 dual value = ', c2.dual_value() + + def RunAllFirstLinearExampleNewAPI(self): + self.RunFirstLinearExampleNewAPI(pywraplp.Solver.GLPK_LINEAR_PROGRAMMING) + self.RunFirstLinearExampleNewAPI(pywraplp.Solver.CLP_LINEAR_PROGRAMMING) + + def RunSuccessiveObjectives(self, mode): + """Example with succesive objectives.""" + solver = pywraplp.Solver('RunSuccessiveObjectives', mode) + x1 = solver.NumVar(0, 10, 'var1') + x2 = solver.NumVar(0, 10, 'var2') + solver.Add(x1 + 2*x2 <= 10) + solver.Maximize(x1) + + # Check the solution + solver.Solve() + print x1.solution_value() + print x2.solution_value() + + solver.Maximize(x2) + + # Check the solution + solver.Solve() + print x1.solution_value() + print x2.solution_value() + + solver.Minimize(-x1) + + # Check the solution + solver.Solve() + print x1.solution_value() + print x2.solution_value() + + def RunAllSuccessiveObjectives(self): + self.RunSuccessiveObjectives(pywraplp.Solver.GLPK_LINEAR_PROGRAMMING) + self.RunSuccessiveObjectives(pywraplp.Solver.CLP_LINEAR_PROGRAMMING) + + def RunAllExamples(self): + self.RunAllFirstLinearExample() + self.RunAllFirstLinearExampleNewAPI() + self.RunAllSuccessiveObjectives() + + +def main(unused_argv): + lp_example = PyWrapLPExamples() + lp_example.RunAllExamples() + + +if __name__ == '__main__': + app.run()