mirrors/z3
3
0
Fork 0
mirror of https://github.com/Z3Prover/z3 synced 2025-04-28 03:15:50 +00:00
Code Activity

Merge branch 'opt' of https://github.com/NikolajBjorner/z3 into opt

Browse source
This commit is contained in:
Nikolaj Bjorner 2017-05-23 10:52:27 -07:00
commit db4bbecb27
174 changed files with 30440 additions and 922 deletions
Download patch file Download diff file Expand all files Collapse all files

44
CMakeLists.txt
View file

@ -188,6 +188,7 @@ endif()
# Note for some reason we have to leave off ``-D`` here otherwise
# we get ``-D-DZ3DEBUG`` passed to the compiler
list(APPEND Z3_COMPONENT_CXX_DEFINES $<$<CONFIG:Debug>:Z3DEBUG>)
list(APPEND Z3_COMPONENT_CXX_DEFINES $<$<CONFIG:Debug>:LEAN_DEBUG>)
list(APPEND Z3_COMPONENT_CXX_DEFINES $<$<CONFIG:Release>:_EXTERNAL_RELEASE>)
list(APPEND Z3_COMPONENT_CXX_DEFINES $<$<CONFIG:RelWithDebInfo>:_EXTERNAL_RELEASE>)
@ -213,6 +214,16 @@ endif()
################################################################################
include(${CMAKE_SOURCE_DIR}/cmake/z3_add_cxx_flag.cmake)
################################################################################
# C++ language version
################################################################################
# FIXME: Use CMake's own mechanism for selecting language version
if (("${CMAKE_CXX_COMPILER_ID}" MATCHES "GNU") OR ("${CMAKE_CXX_COMPILER_ID}" MATCHES "Clang"))
z3_add_cxx_flag("-std=c++11" REQUIRED)
else()
message(AUTHOR_WARNING "Not setting C++ language version for compiler")
endif()
################################################################################
# Platform detection
################################################################################
@ -245,6 +256,7 @@ else()
message(FATAL_ERROR "Platform \"${CMAKE_SYSTEM_NAME}\" not recognised")
endif()
list(APPEND Z3_COMPONENT_EXTRA_INCLUDE_DIRS "${CMAKE_SOURCE_DIR}/src")
################################################################################
# GNU multiple precision library support
@ -263,18 +275,6 @@ else()
message(STATUS "Not using libgmp")
endif()
################################################################################
# FOCI2 support
################################################################################
# FIXME: What is this?
option(USE_FOCI2 "Use FOCI2" OFF)
if (USE_FOCI2)
message(FATAL_ERROR "TODO")
message(STATUS "Using FOCI2")
else()
message(STATUS "Not using FOCI2")
endif()
################################################################################
# OpenMP support
################################################################################
@ -288,6 +288,7 @@ if (USE_OPENMP)
message(WARNING "OpenMP support was requested but your compiler doesn't support it")
endif()
endif()
if (OPENMP_FOUND)
list(APPEND Z3_COMPONENT_CXX_FLAGS ${OpenMP_CXX_FLAGS})
# GCC and Clang need to have additional flags passed to the linker.
@ -306,6 +307,23 @@ else()
set(USE_OPENMP OFF CACHE BOOL "Use OpenMP" FORCE)
endif()
################################################################################
# API Log sync
################################################################################
option(API_LOG_SYNC
"Use locking when logging Z3 API calls (experimental)"
OFF
)
if (API_LOG_SYNC)
if (NOT USE_OPENMP)
message(FATAL_ERROR "API_LOG_SYNC feature requires OpenMP")
endif()
list(APPEND Z3_COMPONENT_CXX_DEFINES "-DZ3_LOG_SYNC")
message(STATUS "Using API_LOG_SYNC")
else()
message(STATUS "Not using API_LOG_SYNC")
endif()
################################################################################
# FP math
################################################################################
@ -326,6 +344,8 @@ if (("${TARGET_ARCHITECTURE}" STREQUAL "x86_64") OR ("${TARGET_ARCHITECTURE}" ST
unset(SSE_FLAGS)
endif()
# FIXME: Remove "x.." when CMP0054 is set to NEW
if ("x${CMAKE_CXX_COMPILER_ID}" STREQUAL "xMSVC")
# This is the default for MSVC already but to replicate the

3
README-CMake.md
View file

@ -1,4 +1,4 @@
# Z3's CMake build system
# Z3's CMake build system
[CMake](https://cmake.org/) is a "meta build system" that reads a description
of the project written in the ``CMakeLists.txt`` files and emits a build
@ -293,6 +293,7 @@ The following useful options can be passed to CMake whilst configuring.
* ``ALWAYS_BUILD_DOCS`` - BOOL. If set to ``TRUE`` and ``BUILD_DOCUMENTATION`` is ``TRUE`` then documentation for API bindings will always be built.
Disabling this is useful for faster incremental builds. The documentation can be manually built by invoking the ``api_docs`` target.
* ``LINK_TIME_OPTIMIZATION`` - BOOL. If set to ``TRUE`` link time optimization will be enabled.
* ``API_LOG_SYNC`` - BOOL. If set to ``TRUE`` will enable experimental API log sync feature.
On the command line these can be passed to ``cmake`` using the ``-D`` option. In ``ccmake`` and ``cmake-gui`` these can be set in the user interface.

1
contrib/cmake/src/CMakeLists.txt
View file

@ -35,6 +35,7 @@ endforeach()
# raised if you try to declare a component is dependent on another component
# that has not yet been declared.
add_subdirectory(util)
add_subdirectory(util/lp)
add_subdirectory(math/polynomial)
add_subdirectory(sat)
add_subdirectory(nlsat)

1
contrib/cmake/src/interp/CMakeLists.txt
View file

@ -2,7 +2,6 @@ z3_add_component(interp
SOURCES
iz3base.cpp
iz3checker.cpp
iz3foci.cpp
iz3interp.cpp
iz3mgr.cpp
iz3pp.cpp

1
contrib/cmake/src/shell/CMakeLists.txt
View file

@ -27,6 +27,7 @@ add_executable(shell
opt_frontend.cpp
smtlib_frontend.cpp
z3_log_frontend.cpp
lp_frontend.cpp
# FIXME: shell should really link against libz3 but it can't due to requiring
# use of some hidden symbols. Also libz3 has the ``api_dll`` component which
# we don't want (I think).

2
contrib/cmake/src/smt/CMakeLists.txt
View file

@ -55,6 +55,7 @@ z3_add_component(smt
theory_dl.cpp
theory_dummy.cpp
theory_fpa.cpp
theory_lra.cpp
theory_opt.cpp
theory_pb.cpp
theory_seq.cpp
@ -69,6 +70,7 @@ z3_add_component(smt
euclid
fpa
grobner
lp
macros
normal_forms
parser_util

4
contrib/cmake/src/test/CMakeLists.txt
View file

@ -119,6 +119,7 @@ add_executable(test-z3
upolynomial.cpp
var_subst.cpp
vector.cpp
lp.cpp
${z3_test_extra_object_files}
)
z3_add_install_tactic_rule(${z3_test_deps})
@ -130,3 +131,6 @@ target_link_libraries(test-z3 PRIVATE ${Z3_DEPENDENT_LIBS})
target_include_directories(test-z3 PRIVATE ${Z3_COMPONENT_EXTRA_INCLUDE_DIRS})
z3_append_linker_flag_list_to_target(test-z3 ${Z3_DEPENDENT_EXTRA_CXX_LINK_FLAGS})
z3_add_component_dependencies_to_target(test-z3 ${z3_test_expanded_deps})

35
contrib/cmake/src/util/lp/CMakeLists.txt Normal file
View file

@ -0,0 +1,35 @@
z3_add_component(lp
SOURCES
lp_utils.cpp
binary_heap_priority_queue_instances.cpp
binary_heap_upair_queue_instances.cpp
bound_propagator.cpp
core_solver_pretty_printer_instances.cpp
dense_matrix_instances.cpp
eta_matrix_instances.cpp
indexed_vector_instances.cpp
lar_core_solver_instances.cpp
lp_core_solver_base_instances.cpp
lp_dual_core_solver_instances.cpp
lp_dual_simplex_instances.cpp
lp_primal_core_solver_instances.cpp
lp_primal_simplex_instances.cpp
lp_settings_instances.cpp
lp_solver_instances.cpp
lu_instances.cpp
matrix_instances.cpp
permutation_matrix_instances.cpp
quick_xplain.cpp
row_eta_matrix_instances.cpp
scaler_instances.cpp
sparse_matrix_instances.cpp
square_dense_submatrix_instances.cpp
static_matrix_instances.cpp
random_updater_instances.cpp
COMPONENT_DEPENDENCIES
util
PYG_FILES
lp_params.pyg
)
include_directories(${src_SOURCE_DIR})

3
scripts/mk_project.py
View file

@ -11,6 +11,7 @@ from mk_util import *
def init_project_def():
set_version(4, 5, 1, 0)
add_lib('util', [])
add_lib('lp', ['util'], 'util/lp')
add_lib('polynomial', ['util'], 'math/polynomial')
add_lib('sat', ['util'])
add_lib('nlsat', ['polynomial', 'sat'])
@ -52,7 +53,7 @@ def init_project_def():
add_lib('smt_params', ['ast', 'simplifier', 'pattern', 'bit_blaster'], 'smt/params')
add_lib('proto_model', ['model', 'simplifier', 'smt_params'], 'smt/proto_model')
add_lib('smt', ['bit_blaster', 'macros', 'normal_forms', 'cmd_context', 'proto_model',
'substitution', 'grobner', 'euclid', 'simplex', 'proof_checker', 'pattern', 'parser_util', 'fpa'])
'substitution', 'grobner', 'euclid', 'simplex', 'proof_checker', 'pattern', 'parser_util', 'fpa', 'lp'])
add_lib('bv_tactics', ['tactic', 'bit_blaster', 'core_tactics'], 'tactic/bv')
add_lib('fuzzing', ['ast'], 'test/fuzzing')
add_lib('smt_tactic', ['smt'], 'smt/tactic')

45
scripts/mk_util.py
View file

@ -96,8 +96,6 @@ VER_BUILD=None
VER_REVISION=None
PREFIX=sys.prefix
GMP=False
FOCI2=False
FOCI2LIB=''
VS_PAR=False
VS_PAR_NUM=8
GPROF=False
@ -257,13 +255,6 @@ def test_gmp(cc):
t.commit()
return exec_compiler_cmd([cc, CPPFLAGS, 'tstgmp.cpp', LDFLAGS, '-lgmp']) == 0
def test_foci2(cc,foci2lib):
if is_verbose():
print("Testing FOCI2...")
t = TempFile('tstfoci2.cpp')
t.add('#include<foci2.h>\nint main() { foci2 *f = foci2::create("lia"); return 0; }\n')
t.commit()
return exec_compiler_cmd([cc, CPPFLAGS, '-Isrc/interp', 'tstfoci2.cpp', LDFLAGS, foci2lib]) == 0
def test_openmp(cc):
if not USE_OMP:
@ -650,7 +641,6 @@ def display_help(exit_code):
if not IS_WINDOWS:
print(" -g, --gmp use GMP.")
print(" --gprof enable gprof")
print(" -f <path> --foci2=<path> use foci2 library at path")
print(" --noomp disable OpenMP and all features that require it.")
print(" --log-sync synchronize access to API log files to enable multi-thread API logging.")
print("")
@ -678,13 +668,13 @@ def display_help(exit_code):
# Parse configuration option for mk_make script
def parse_options():
global VERBOSE, DEBUG_MODE, IS_WINDOWS, VS_X64, ONLY_MAKEFILES, SHOW_CPPS, VS_PROJ, TRACE, VS_PAR, VS_PAR_NUM
global DOTNET_ENABLED, DOTNET_KEY_FILE, JAVA_ENABLED, ML_ENABLED, STATIC_LIB, STATIC_BIN, PREFIX, GMP, FOCI2, FOCI2LIB, PYTHON_PACKAGE_DIR, GPROF, GIT_HASH, GIT_DESCRIBE, PYTHON_INSTALL_ENABLED, PYTHON_ENABLED
global DOTNET_ENABLED, DOTNET_KEY_FILE, JAVA_ENABLED, ML_ENABLED, STATIC_LIB, STATIC_BIN, PREFIX, GMP, PYTHON_PACKAGE_DIR, GPROF, GIT_HASH, GIT_DESCRIBE, PYTHON_INSTALL_ENABLED, PYTHON_ENABLED
global LINUX_X64, SLOW_OPTIMIZE, USE_OMP, LOG_SYNC
try:
options, remainder = getopt.gnu_getopt(sys.argv[1:],
'b:df:sxhmcvtnp:gj',
['build=', 'debug', 'silent', 'x64', 'help', 'makefiles', 'showcpp', 'vsproj',
'trace', 'dotnet', 'dotnet-key=', 'staticlib', 'prefix=', 'gmp', 'foci2=', 'java', 'parallel=', 'gprof',
'trace', 'dotnet', 'dotnet-key=', 'staticlib', 'prefix=', 'gmp', 'java', 'parallel=', 'gprof',
'githash=', 'git-describe', 'x86', 'ml', 'optimize', 'noomp', 'pypkgdir=', 'python', 'staticbin', 'log-sync'])
except:
print("ERROR: Invalid command line option")
@ -735,9 +725,6 @@ def parse_options():
VS_PAR_NUM = int(arg)
elif opt in ('-g', '--gmp'):
GMP = True
elif opt in ('-f', '--foci2'):
FOCI2 = True
FOCI2LIB = arg
elif opt in ('-j', '--java'):
JAVA_ENABLED = True
elif opt == '--gprof':
@ -774,8 +761,13 @@ def extract_c_includes(fname):
linenum = 1
for line in f:
m1 = std_inc_pat.match(line)
if m1:
result.append(m1.group(1))
if m1:
root_file_name = m1.group(1)
slash_pos = root_file_name.rfind('/')
if slash_pos >= 0 and root_file_name.find("..") < 0 : #it is a hack for lp include files that behave as continued from "src"
# print(root_file_name)
root_file_name = root_file_name[slash_pos+1:]
result.append(root_file_name)
elif not system_inc_pat.match(line) and non_std_inc_pat.match(line):
raise MKException("Invalid #include directive at '%s':%s" % (fname, line))
linenum = linenum + 1
@ -999,6 +991,7 @@ class Component:
out.write('%s =' % include_defs)
for dep in self.deps:
out.write(' -I%s' % get_component(dep).to_src_dir)
out.write(' -I%s' % os.path.join(REV_BUILD_DIR,"src"))
out.write('\n')
mk_dir(os.path.join(BUILD_DIR, self.build_dir))
if VS_PAR and IS_WINDOWS:
@ -1175,7 +1168,6 @@ class ExeComponent(Component):
for dep in deps:
c_dep = get_component(dep)
out.write(' ' + c_dep.get_link_name())
out.write(' ' + FOCI2LIB)
out.write(' $(LINK_EXTRA_FLAGS)\n')
out.write('%s: %s\n\n' % (self.name, exefile))
@ -1301,7 +1293,6 @@ class DLLComponent(Component):
if dep not in self.reexports:
c_dep = get_component(dep)
out.write(' ' + c_dep.get_link_name())
out.write(' ' + FOCI2LIB)
out.write(' $(SLINK_EXTRA_FLAGS)')
if IS_WINDOWS:
out.write(' /DEF:%s.def' % os.path.join(self.to_src_dir, self.name))
@ -2301,7 +2292,7 @@ def mk_config():
if ONLY_MAKEFILES:
return
config = open(os.path.join(BUILD_DIR, 'config.mk'), 'w')
global CXX, CC, GMP, FOCI2, CPPFLAGS, CXXFLAGS, LDFLAGS, EXAMP_DEBUG_FLAG, FPMATH_FLAGS, HAS_OMP, LOG_SYNC
global CXX, CC, GMP, CPPFLAGS, CXXFLAGS, LDFLAGS, EXAMP_DEBUG_FLAG, FPMATH_FLAGS, HAS_OMP, LOG_SYNC
if IS_WINDOWS:
config.write(
'CC=cl\n'
@ -2411,14 +2402,6 @@ def mk_config():
SLIBEXTRAFLAGS = '%s -lgmp' % SLIBEXTRAFLAGS
else:
CPPFLAGS = '%s -D_MP_INTERNAL' % CPPFLAGS
if FOCI2:
if test_foci2(CXX,FOCI2LIB):
LDFLAGS = '%s %s' % (LDFLAGS,FOCI2LIB)
SLIBEXTRAFLAGS = '%s %s' % (SLIBEXTRAFLAGS,FOCI2LIB)
CPPFLAGS = '%s -D_FOCI2' % CPPFLAGS
else:
print("FAILED\n")
FOCI2 = False
if GIT_HASH:
CPPFLAGS = '%s -DZ3GITHASH=%s' % (CPPFLAGS, GIT_HASH)
CXXFLAGS = '%s -std=c++11' % CXXFLAGS
@ -2953,7 +2936,11 @@ def get_platform_toolset_str():
if len(tokens) < 2:
return default
else:
return 'v' + tokens[0] + tokens[1]
if tokens[0] == "15":
# Visual Studio 2017 reports 15.* but the PlatformToolsetVersion is 141
return "v141"
else:
return 'v' + tokens[0] + tokens[1]
def mk_vs_proj_property_groups(f, name, target_ext, type):
f.write(' <ItemGroup Label="ProjectConfigurations">\n')

2
src/api/dotnet/core/DummyContracts.cs
View file

@ -62,4 +62,4 @@ namespace System.Diagnostics.Contracts
public static void EndContractBlock() { }
public static T ValueAtReturn<T>(out T v) { T[] t = new T[1]; v = t[0]; return v; }
}
}
}

2
src/api/dotnet/dotnet35/Example/App.config
View file

@ -3,4 +3,4 @@
<startup>
<supportedRuntime version="v4.0" sku=".NETFramework,Version=v4.5.2" />
</startup>
</configuration>
</configuration>

2
src/api/dotnet/dotnet35/Example/Example.csproj
View file

@ -75,4 +75,4 @@
<Target Name="AfterBuild">
</Target>
-->
</Project>
</Project>

2
src/api/dotnet/dotnet35/Microsoft.Z3.NET35.csproj
View file

@ -344,4 +344,4 @@
<None Include="packages.config" />
</ItemGroup>
<Import Project="$(MSBuildToolsPath)\Microsoft.CSharp.targets" />
</Project>
</Project>

2
src/api/dotnet/dotnet35/packages.config
View file

@ -1,4 +1,4 @@
<?xml version="1.0" encoding="utf-8"?>
<packages>
<package id="Code.Contract" version="1.0.0" targetFramework="net35" />
</packages>
</packages>

2
src/api/python/setup.py
View file

@ -164,6 +164,6 @@ setup(
package_data={
'z3': [os.path.join('lib', '*'), os.path.join('include', '*.h'), os.path.join('include', 'c++', '*.h')]
},
scripts=[os.path.join('bin', EXECUTABLE_FILE)],
data_files=[('bin',[os.path.join('bin',EXECUTABLE_FILE)])],
cmdclass={'build': build, 'develop': develop, 'sdist': sdist, 'bdist_egg': bdist_egg},
)

6
src/api/python/z3/z3.py
View file

@ -1300,8 +1300,10 @@ class BoolSortRef(SortRef):
if isinstance(val, bool):
return BoolVal(val, self.ctx)
if __debug__:
_z3_assert(is_expr(val), "True, False or Z3 Boolean expression expected. Received %s" % val)
_z3_assert(self.eq(val.sort()), "Value cannot be converted into a Z3 Boolean value")
if not is_expr(val):
_z3_assert(is_expr(val), "True, False or Z3 Boolean expression expected. Received %s" % val)
if not self.eq(val.sort()):
_z3_assert(self.eq(val.sort()), "Value cannot be converted into a Z3 Boolean value")
return val
def subsort(self, other):

4
src/ast/arith_decl_plugin.h
View file

@ -305,6 +305,7 @@ public:
MATCH_UNARY(is_uminus);
MATCH_UNARY(is_to_real);
MATCH_UNARY(is_to_int);
MATCH_UNARY(is_is_int);
MATCH_BINARY(is_sub);
MATCH_BINARY(is_add);
MATCH_BINARY(is_mul);
@ -377,6 +378,9 @@ public:
app * mk_real(int i) {
return mk_numeral(rational(i), false);
}
app * mk_real(rational const& r) {
return mk_numeral(r, false);
}
app * mk_le(expr * arg1, expr * arg2) const { return m_manager.mk_app(m_afid, OP_LE, arg1, arg2); }
app * mk_ge(expr * arg1, expr * arg2) const { return m_manager.mk_app(m_afid, OP_GE, arg1, arg2); }
app * mk_lt(expr * arg1, expr * arg2) const { return m_manager.mk_app(m_afid, OP_LT, arg1, arg2); }

7
src/ast/ast_pp_util.cpp
View file

@ -44,8 +44,11 @@ void ast_pp_util::display_decls(std::ostream& out) {
}
n = coll.get_num_decls();
for (unsigned i = 0; i < n; ++i) {
ast_smt2_pp(out, coll.get_func_decls()[i], env);
out << "\n";
func_decl* f = coll.get_func_decls()[i];
if (f->get_family_id() == null_family_id) {
ast_smt2_pp(out, f, env);
out << "\n";
}
}
}

117
src/ast/rewriter/pb2bv_rewriter.cpp
View file

@ -26,6 +26,7 @@ Notes:
#include"ast_pp.h"
#include"lbool.h"
#include"uint_set.h"
#include"gparams.h"
const unsigned g_primes[7] = { 2, 3, 5, 7, 11, 13, 17};
@ -54,6 +55,8 @@ struct pb2bv_rewriter::imp {
vector<rational> m_coeffs;
bool m_keep_cardinality_constraints;
bool m_keep_pb_constraints;
bool m_pb_num_system;
bool m_pb_totalizer;
unsigned m_min_arity;
template<lbool is_le>
@ -114,23 +117,24 @@ struct pb2bv_rewriter::imp {
return expr_ref((is_le == l_false)?m.mk_true():m.mk_false(), m);
}
#if 0
expr_ref result(m);
switch (is_le) {
case l_true: if (mk_le_tot(sz, args, k, result)) return result; else break;
case l_false: if (mk_ge_tot(sz, args, k, result)) return result; else break;
case l_undef: break;
if (m_pb_totalizer) {
expr_ref result(m);
switch (is_le) {
case l_true: if (mk_le_tot(sz, args, k, result)) return result; else break;
case l_false: if (mk_ge_tot(sz, args, k, result)) return result; else break;
case l_undef: break;
}
}
if (m_pb_num_system) {
expr_ref result(m);
switch (is_le) {
case l_true: if (mk_le(sz, args, k, result)) return result; else break;
case l_false: if (mk_ge(sz, args, k, result)) return result; else break;
case l_undef: if (mk_eq(sz, args, k, result)) return result; else break;
}
}
#endif
#if 0
expr_ref result(m);
switch (is_le) {
case l_true: if (mk_le(sz, args, k, result)) return result; else break;
case l_false: if (mk_ge(sz, args, k, result)) return result; else break;
case l_undef: if (mk_eq(sz, args, k, result)) return result; else break;
}
#endif
// fall back to divide and conquer encoding.
SASSERT(k.is_pos());
expr_ref zero(m), bound(m);
@ -397,22 +401,31 @@ struct pb2bv_rewriter::imp {
for (unsigned j = 0; j + lim - 1 < out.size(); j += n) {
expr_ref tmp(m);
tmp = out[j + lim - 1];
if (j + n < out.size()) {
tmp = m.mk_and(tmp, m.mk_not(out[j + n]));
if (j + n - 1 < out.size()) {
tmp = m.mk_and(tmp, m.mk_not(out[j + n - 1]));
}
ors.push_back(tmp);
}
return ::mk_or(ors);
}
// x0 + 5x1 + 3x2 >= k
// x0 x1 x1 -> s0 s1 s2
// s2 x1 x2 -> s3 s4 s5
// k = 7: s5 or (s4 & not s2 & s0)
// k = 6: s4
// k = 5: s4 or (s3 & not s2 & s1)
// k = 4: s4 or (s3 & not s2 & s0)
// k = 3: s3
//
bool mk_ge(unsigned sz, expr * const* args, rational bound, expr_ref& result) {
if (!create_basis()) return false;
if (!bound.is_unsigned()) return false;
vector<rational> coeffs(m_coeffs);
result = m.mk_true();
expr_ref_vector carry(m), new_carry(m);
for (unsigned i = 0; i < m_base.size(); ++i) {
rational b_i = m_base[i];
m_base.push_back(bound + rational::one());
for (rational b_i : m_base) {
unsigned B = b_i.get_unsigned();
unsigned d_i = (bound % b_i).get_unsigned();
bound = div(bound, b_i);
@ -425,16 +438,16 @@ struct pb2bv_rewriter::imp {
coeffs[j] = div(coeffs[j], b_i);
}
TRACE("pb", tout << "Carry: " << carry << "\n";
for (unsigned j = 0; j < coeffs.size(); ++j) tout << coeffs[j] << " ";
for (auto c : coeffs) tout << c << " ";
tout << "\n";
);
ptr_vector<expr> out;
m_sort.sorting(carry.size(), carry.c_ptr(), out);
expr_ref gt = mod_ge(out, B, d_i + 1);
expr_ref ge = mod_ge(out, B, d_i);
result = mk_or(gt, mk_and(ge, result));
TRACE("pb", tout << result << "\n";);
TRACE("pb", tout << "b: " << b_i << " d: " << d_i << " gt: " << gt << " ge: " << ge << " " << result << "\n";);
new_carry.reset();
for (unsigned j = B - 1; j < out.size(); j += B) {
@ -443,11 +456,11 @@ struct pb2bv_rewriter::imp {
carry.reset();
carry.append(new_carry);
}
TRACE("pb", tout << result << "\n";);
TRACE("pb", tout << "bound: " << bound << " Carry: " << carry << " result: " << result << "\n";);
return true;
}
expr_ref mk_and(expr_ref& a, expr_ref& b) {
expr_ref mk_and(expr_ref& a, expr_ref& b) {
if (m.is_true(a)) return b;
if (m.is_true(b)) return a;
if (m.is_false(a)) return a;
@ -491,6 +504,7 @@ struct pb2bv_rewriter::imp {
for (unsigned i = 0; i < sz; ++i) {
m_coeffs.push_back(pb.get_coeff(f, i));
}
CTRACE("pb", k.is_neg(), tout << expr_ref(m.mk_app(f, sz, args), m) << "\n";);
SASSERT(!k.is_neg());
switch (kind) {
case OP_PB_GE:
@ -567,6 +581,8 @@ struct pb2bv_rewriter::imp {
m_args(m),
m_keep_cardinality_constraints(false),
m_keep_pb_constraints(false),
m_pb_num_system(false),
m_pb_totalizer(false),
m_min_arity(2)
{}
@ -761,6 +777,14 @@ struct pb2bv_rewriter::imp {
void keep_pb_constraints(bool f) {
m_keep_pb_constraints = f;
}
void pb_num_system(bool f) {
m_pb_num_system = f;
}
void pb_totalizer(bool f) {
m_pb_totalizer = f;
}
};
struct card2bv_rewriter_cfg : public default_rewriter_cfg {
@ -774,6 +798,8 @@ struct pb2bv_rewriter::imp {
card2bv_rewriter_cfg(imp& i, ast_manager & m):m_r(i, m) {}
void keep_cardinality_constraints(bool f) { m_r.keep_cardinality_constraints(f); }
void keep_pb_constraints(bool f) { m_r.keep_pb_constraints(f); }
void pb_num_system(bool f) { m_r.pb_num_system(f); }
void pb_totalizer(bool f) { m_r.pb_totalizer(f); }
};
class card_pb_rewriter : public rewriter_tpl<card2bv_rewriter_cfg> {
@ -784,27 +810,62 @@ struct pb2bv_rewriter::imp {
m_cfg(i, m) {}
void keep_cardinality_constraints(bool f) { m_cfg.keep_cardinality_constraints(f); }
void keep_pb_constraints(bool f) { m_cfg.keep_pb_constraints(f); }
void pb_num_system(bool f) { m_cfg.pb_num_system(f); }
void pb_totalizer(bool f) { m_cfg.pb_totalizer(f); }
};
card_pb_rewriter m_rw;
bool keep_cardinality() const {
params_ref const& p = m_params;
return
p.get_bool("keep_cardinality_constraints", false) ||
p.get_bool("sat.cardinality.solver", false) ||
p.get_bool("cardinality.solver", false) ||
gparams::get_module("sat").get_bool("cardinality.solver", false) ||
keep_pb();
}
bool keep_pb() const {
params_ref const& p = m_params;
return
p.get_bool("keep_pb_constraints", false) ||
p.get_bool("sat.pb.solver", false) ||
p.get_bool("pb.solver", false) ||
gparams::get_module("sat").get_bool("pb.solver", false);
}
bool pb_num_system() const {
return m_params.get_bool("pb_num_system", false);
}
bool pb_totalizer() const {
return m_params.get_bool("pb_totalizer", false);
}
imp(ast_manager& m, params_ref const& p):
m(m), m_params(p), m_lemmas(m),
m_fresh(m),
m_num_translated(0),
m_rw(*this, m) {
m_rw.keep_cardinality_constraints(p.get_bool("keep_cardinality_constraints", false));
m_rw.keep_pb_constraints(p.get_bool("keep_pb_constraints", false));
m_rw.keep_cardinality_constraints(keep_cardinality());
m_rw.keep_pb_constraints(keep_pb());
m_rw.pb_num_system(pb_num_system());
m_rw.pb_totalizer(pb_totalizer());
}
void updt_params(params_ref const & p) {
m_params.append(p);
m_rw.keep_cardinality_constraints(m_params.get_bool("keep_cardinality_constraints", false));
m_rw.keep_pb_constraints(m_params.get_bool("keep_pb_constraints", false));
m_rw.keep_cardinality_constraints(keep_cardinality());
m_rw.keep_pb_constraints(keep_pb());
m_rw.pb_num_system(pb_num_system());
m_rw.pb_totalizer(pb_totalizer());
}
void collect_param_descrs(param_descrs& r) const {
r.insert("keep_cardinality_constraints", CPK_BOOL, "(default: true) retain cardinality constraints (don't bit-blast them) and use built-in cardinality solver");
r.insert("keep_pb_constraints", CPK_BOOL, "(default: true) retain pb constraints (don't bit-blast them) and use built-in pb solver");
r.insert("pb_num_system", CPK_BOOL, "(default: false) use pb number system encoding");
r.insert("pb_totalizer", CPK_BOOL, "(default: false) use pb totalizer encoding");
}
unsigned get_num_steps() const { return m_rw.get_num_steps(); }

43
src/ast/rewriter/pb_rewriter.cpp
View file

@ -232,6 +232,7 @@ br_status pb_rewriter::mk_app_core(func_decl * f, unsigned num_args, expr * cons
}
bool is_eq = f->get_decl_kind() == OP_PB_EQ;
br_status st = BR_DONE;
pb_ast_rewriter_util pbu(m);
pb_rewriter_util<pb_ast_rewriter_util> util(pbu);
@ -249,11 +250,14 @@ br_status pb_rewriter::mk_app_core(func_decl * f, unsigned num_args, expr * cons
default: {
bool all_unit = true;
unsigned sz = vec.size();
rational slack(0);
m_args.reset();
m_coeffs.reset();
for (unsigned i = 0; i < sz; ++i) {
m_args.push_back(vec[i].first);
m_coeffs.push_back(vec[i].second);
SASSERT(vec[i].second.is_pos());
slack += vec[i].second;
all_unit &= m_coeffs.back().is_one();
}
if (is_eq) {
@ -271,7 +275,42 @@ br_status pb_rewriter::mk_app_core(func_decl * f, unsigned num_args, expr * cons
result = mk_and(m, sz, m_args.c_ptr());
}
else {
result = m_util.mk_ge(sz, m_coeffs.c_ptr(), m_args.c_ptr(), k);
expr_ref_vector conj(m), disj(m);
for (unsigned i = 0; i < m_args.size(); ++i) {
rational& c = m_coeffs[i];
if (slack < c + k) {
conj.push_back(m_args[i]);
slack -= c;
k -= c;
}
else if (c >= k) {
slack -= c;
disj.push_back(m_args[i]);
}
else {
continue;
}
m_args[i] = m_args.back();
m_coeffs[i] = m_coeffs.back();
--i;
m_args.pop_back();
m_coeffs.pop_back();
}
sz = m_coeffs.size();
if (slack < k) {
conj.push_back(m.mk_false());
}
else if (k.is_pos() && sz > 0) {
conj.push_back(m_util.mk_ge(sz, m_coeffs.c_ptr(), m_args.c_ptr(), k));
}
result = mk_and(conj);
if (!disj.empty()) {
disj.push_back(result);
result = mk_or(disj);
}
if (!disj.empty() || conj.size() > 1) {
st = BR_REWRITE3;
}
}
break;
}
@ -286,7 +325,7 @@ br_status pb_rewriter::mk_app_core(func_decl * f, unsigned num_args, expr * cons
TRACE("pb_validate",
validate_rewrite(f, num_args, args, result););
return BR_DONE;
return st;
}

56
src/ast/rewriter/seq_rewriter.cpp
View file

@ -440,24 +440,16 @@ br_status seq_rewriter::mk_app_core(func_decl * f, unsigned num_args, expr * con
* (seq.unit (_ BitVector 8)) ==> String constant
*/
br_status seq_rewriter::mk_seq_unit(expr* e, expr_ref& result) {
sort * s = m().get_sort(e);
bv_util bvu(m());
if (is_sort_of(s, bvu.get_family_id(), BV_SORT)) {
// specifically we want (_ BitVector 8)
rational n_val;
unsigned int n_size;
if (bvu.is_numeral(e, n_val, n_size)) {
if (n_size == 8) {
// convert to string constant
char ch = (char)n_val.get_int32();
TRACE("seq", tout << "rewrite seq.unit of 8-bit value " << n_val.to_string() << " to string constant \"" << ch << "\"" << std::endl;);
char s_tmp[2] = {ch, '\0'};
symbol s(s_tmp);
result = m_util.str.mk_string(s);
return BR_DONE;
}
}
rational n_val;
unsigned int n_size;
// specifically we want (_ BitVector 8)
if (bvu.is_bv(e) && bvu.is_numeral(e, n_val, n_size) && n_size == 8) {
// convert to string constant
zstring str(n_val.get_unsigned());
TRACE("seq", tout << "rewrite seq.unit of 8-bit value " << n_val.to_string() << " to string constant \"" << str<< "\"" << std::endl;);
result = m_util.str.mk_string(str);
return BR_DONE;
}
return BR_FAILED;
@ -1431,37 +1423,10 @@ br_status seq_rewriter::mk_re_star(expr* a, expr_ref& result) {
}
/*
* (re.range c_1 c_n) = (re.union (str.to.re c1) (str.to.re c2) ... (str.to.re cn))
* (re.range c_1 c_n)
*/
br_status seq_rewriter::mk_re_range(expr* lo, expr* hi, expr_ref& result) {
return BR_FAILED;
TRACE("seq", tout << "rewrite re.range [" << mk_pp(lo, m()) << " " << mk_pp(hi, m()) << "]\n";);
zstring str_lo, str_hi;
if (m_util.str.is_string(lo, str_lo) && m_util.str.is_string(hi, str_hi)) {
if (str_lo.length() == 1 && str_hi.length() == 1) {
unsigned int c1 = str_lo[0];
unsigned int c2 = str_hi[0];
if (c1 > c2) {
// exchange c1 and c2
unsigned int tmp = c1;
c2 = c1;
c1 = tmp;
}
zstring s(c1);
expr_ref acc(m_util.re.mk_to_re(m_util.str.mk_string(s)), m());
for (unsigned int ch = c1 + 1; ch <= c2; ++ch) {
zstring s_ch(ch);
expr_ref acc2(m_util.re.mk_to_re(m_util.str.mk_string(s_ch)), m());
acc = m_util.re.mk_union(acc, acc2);
}
result = acc;
return BR_REWRITE2;
} else {
m().raise_exception("string constants in re.range must have length 1");
}
}
return BR_FAILED;
}
/*
@ -1506,6 +1471,7 @@ br_status seq_rewriter::mk_re_opt(expr* a, expr_ref& result) {
}
br_status seq_rewriter::mk_eq_core(expr * l, expr * r, expr_ref & result) {
TRACE("seq", tout << mk_pp(l, m()) << " = " << mk_pp(r, m()) << "\n";);
expr_ref_vector lhs(m()), rhs(m()), res(m());
bool changed = false;
if (!reduce_eq(l, r, lhs, rhs, changed)) {

30
src/ast/seq_decl_plugin.cpp
View file

@ -876,6 +876,36 @@ bool seq_decl_plugin::is_value(app* e) const {
}
}
bool seq_decl_plugin::are_equal(app* a, app* b) const {
if (a == b) return true;
// handle concatenations
return false;
}
bool seq_decl_plugin::are_distinct(app* a, app* b) const {
if (a == b) {
return false;
}
if (is_app_of(a, m_family_id, OP_STRING_CONST) &&
is_app_of(b, m_family_id, OP_STRING_CONST)) {
return true;
}
if (is_app_of(a, m_family_id, OP_SEQ_UNIT) &&
is_app_of(b, m_family_id, OP_SEQ_UNIT)) {
return true;
}
if (is_app_of(a, m_family_id, OP_SEQ_EMPTY) &&
is_app_of(b, m_family_id, OP_SEQ_UNIT)) {
return true;
}
if (is_app_of(b, m_family_id, OP_SEQ_EMPTY) &&
is_app_of(a, m_family_id, OP_SEQ_UNIT)) {
return true;
}
return false;
}
expr* seq_decl_plugin::get_some_value(sort* s) {
seq_util util(*m_manager);
if (util.is_seq(s)) {

6
src/ast/seq_decl_plugin.h
View file

@ -182,7 +182,11 @@ public:
virtual bool is_value(app * e) const;
virtual bool is_unique_value(app * e) const { return is_value(e); }
virtual bool is_unique_value(app * e) const { return false; }
virtual bool are_equal(app* a, app* b) const;
virtual bool are_distinct(app* a, app* b) const;
virtual expr * get_some_value(sort * s);

75
src/interp/foci2.h
View file

@ -1,75 +0,0 @@
/*++
Copyright (c) 2011 Microsoft Corporation
Module Name:
foci2.h
Abstract:
An interface class for foci2.
Author:
Ken McMillan (kenmcmil)
Revision History:
--*/
#ifndef FOCI2_H
#define FOCI2_H
#include <vector>
#include <string>
#ifdef _WINDOWS
#define FOCI2_EXPORT __declspec(dllexport)
#else
#define FOCI2_EXPORT __attribute__ ((visibility ("default")))
#endif
class foci2 {
public:
virtual ~foci2(){}
typedef int ast;
typedef int symb;
/** Built-in operators */
enum ops {
And = 0, Or, Not, Iff, Ite, Equal, Plus, Times, Floor, Leq, Div, Bool, Int, Array, Tsym, Fsym, Forall, Exists, Distinct, LastOp
};
virtual symb mk_func(const std::string &s) = 0;
virtual symb mk_pred(const std::string &s) = 0;
virtual ast mk_op(ops op, const std::vector<ast> args) = 0;
virtual ast mk_op(ops op, ast) = 0;
virtual ast mk_op(ops op, ast, ast) = 0;
virtual ast mk_op(ops op, ast, ast, ast) = 0;
virtual ast mk_int(const std::string &) = 0;
virtual ast mk_rat(const std::string &) = 0;
virtual ast mk_true() = 0;
virtual ast mk_false() = 0;
virtual ast mk_app(symb,const std::vector<ast> args) = 0;
virtual bool get_func(ast, symb &) = 0;
virtual bool get_pred(ast, symb &) = 0;
virtual bool get_op(ast, ops &) = 0;
virtual bool get_true(ast id) = 0;
virtual bool get_false(ast id) = 0;
virtual bool get_int(ast id, std::string &res) = 0;
virtual bool get_rat(ast id, std::string &res) = 0;
virtual const std::string &get_symb(symb) = 0;
virtual int get_num_args(ast) = 0;
virtual ast get_arg(ast, int) = 0;
virtual void show_ast(ast) = 0;
virtual bool interpolate(const std::vector<ast> &frames, std::vector<ast> &itps, std::vector<int> parents) = 0;
FOCI2_EXPORT static foci2 *create(const std::string &);
};
#endif

25
src/interp/foci2stub/foci2.cpp
View file

@ -1,25 +0,0 @@
/*++
Copyright (c) 2011 Microsoft Corporation
Module Name:
foci2.cpp
Abstract:
Fake foci2, to be replaced
Author:
Ken McMillan (kenmcmil)
Revision History:
--*/
#include "foci2.h"
FOCI2_EXPORT foci2 *foci2::create(const std::string &){
return 0;
}

75
src/interp/foci2stub/foci2.h
View file

@ -1,75 +0,0 @@
/*++
Copyright (c) 2011 Microsoft Corporation
Module Name:
foci2.h
Abstract:
An interface class for foci2.
Author:
Ken McMillan (kenmcmil)
Revision History:
--*/
#ifndef FOCI2_H
#define FOCI2_H
#include <vector>
#include <string>
#ifdef WIN32
#define FOCI2_EXPORT __declspec(dllexport)
#else
#define FOCI2_EXPORT __attribute__ ((visibility ("default")))
#endif
class foci2 {
public:
virtual ~foci2(){}
typedef int ast;
typedef int symb;
/** Built-in operators */
enum ops {
And = 0, Or, Not, Iff, Ite, Equal, Plus, Times, Floor, Leq, Div, Bool, Int, Array, Tsym, Fsym, Forall, Exists, Distinct, LastOp
};
virtual symb mk_func(const std::string &s) = 0;
virtual symb mk_pred(const std::string &s) = 0;
virtual ast mk_op(ops op, const std::vector<ast> args) = 0;
virtual ast mk_op(ops op, ast) = 0;
virtual ast mk_op(ops op, ast, ast) = 0;
virtual ast mk_op(ops op, ast, ast, ast) = 0;
virtual ast mk_int(const std::string &) = 0;
virtual ast mk_rat(const std::string &) = 0;
virtual ast mk_true() = 0;
virtual ast mk_false() = 0;
virtual ast mk_app(symb,const std::vector<ast> args) = 0;
virtual bool get_func(ast, symb &) = 0;
virtual bool get_pred(ast, symb &) = 0;
virtual bool get_op(ast, ops &) = 0;
virtual bool get_true(ast id) = 0;
virtual bool get_false(ast id) = 0;
virtual bool get_int(ast id, std::string &res) = 0;
virtual bool get_rat(ast id, std::string &res) = 0;
virtual const std::string &get_symb(symb) = 0;
virtual int get_num_args(ast) = 0;
virtual ast get_arg(ast, int) = 0;
virtual void show_ast(ast) = 0;
virtual bool interpolate(const std::vector<ast> &frames, std::vector<ast> &itps, std::vector<int> parents) = 0;
FOCI2_EXPORT static foci2 *create(const std::string &);
};
#endif

356
src/interp/iz3foci.cpp
View file

@ -1,356 +0,0 @@
/*++
Copyright (c) 2011 Microsoft Corporation
Module Name:
iz3foci.cpp
Abstract:
Implements a secondary solver using foci2.
Author:
Ken McMillan (kenmcmil)
Revision History:
--*/
#include <sstream>
#include <iostream>
#include "iz3hash.h"
#include "foci2.h"
#include "iz3foci.h"
using namespace stl_ext;
class iz3foci_impl : public iz3secondary {
int frames;
int *parents;
foci2 *foci;
foci2::symb select_op;
foci2::symb store_op;
foci2::symb mod_op;
public:
iz3foci_impl(iz3mgr *mgr, int _frames, int *_parents) : iz3secondary(*mgr) {
frames = _frames;
parents = _parents;
foci = 0;
}
typedef hash_map<ast,foci2::ast> AstToNode;
AstToNode ast_to_node; // maps Z3 ast's to foci expressions
typedef hash_map<foci2::ast,ast> NodeToAst;
NodeToAst node_to_ast; // maps Z3 ast's to foci expressions
// We only use this for FuncDeclToSymbol, which has no range destructor
struct symb_hash {
size_t operator()(const symb &s) const {
return (size_t) s;
}
};
typedef hash_map<symb,foci2::symb> FuncDeclToSymbol;
FuncDeclToSymbol func_decl_to_symbol; // maps Z3 func decls to symbols
typedef hash_map<foci2::symb,symb> SymbolToFuncDecl;
SymbolToFuncDecl symbol_to_func_decl; // maps symbols to Z3 func decls
int from_symb(symb func){
std::string name = string_of_symbol(func);
bool is_bool = is_bool_type(get_range_type(func));
foci2::symb f;
if(is_bool)
f = foci->mk_pred(name);
else
f = foci->mk_func(name);
symbol_to_func_decl[f] = func;
func_decl_to_symbol[func] = f;
return f;
}
// create a symbol corresponding to a DeBruijn index of a particular type
// the type has to be encoded into the name because the same index can
// occur with different types
foci2::symb make_deBruijn_symbol(int index, type ty){
std::ostringstream s;
// s << "#" << index << "#" << type;
return foci->mk_func(s.str());
}
int from_Z3_ast(ast t){
std::pair<ast,foci2::ast> foo(t,0);
std::pair<AstToNode::iterator, bool> bar = ast_to_node.insert(foo);
int &res = bar.first->second;
if(!bar.second) return res;
int nargs = num_args(t);
std::vector<foci2::ast> args(nargs);
for(int i = 0; i < nargs; i++)
args[i] = from_Z3_ast(arg(t,i));
switch(op(t)){
case True:
res = foci->mk_true(); break;
case False:
res = foci->mk_false(); break;
case And:
res = foci->mk_op(foci2::And,args); break;
case Or:
res = foci->mk_op(foci2::Or,args); break;
case Not:
res = foci->mk_op(foci2::Not,args[0]); break;
case Iff:
res = foci->mk_op(foci2::Iff,args); break;
case OP_OEQ: // bit of a mystery, this one...
if(args[0] == args[1]) res = foci->mk_true();
else res = foci->mk_op(foci2::Iff,args);
break;
case Ite:
if(is_bool_type(get_type(t)))
res = foci->mk_op(foci2::And,foci->mk_op(foci2::Or,foci->mk_op(foci2::Not,args[0]),args[1]),foci->mk_op(foci2::Or,args[0],args[2]));
else
res = foci->mk_op(foci2::Ite,args);
break;
case Equal:
res = foci->mk_op(foci2::Equal,args); break;
case Implies:
args[0] = foci->mk_op(foci2::Not,args[0]); res = foci->mk_op(foci2::Or,args); break;
case Xor:
res = foci->mk_op(foci2::Not,foci->mk_op(foci2::Iff,args)); break;
case Leq:
res = foci->mk_op(foci2::Leq,args); break;
case Geq:
std::swap(args[0],args[1]); res = foci->mk_op(foci2::Leq,args); break;
case Gt:
res = foci->mk_op(foci2::Not,foci->mk_op(foci2::Leq,args)); break;
case Lt:
std::swap(args[0],args[1]); res = foci->mk_op(foci2::Not,foci->mk_op(foci2::Leq,args)); break;
case Plus:
res = foci->mk_op(foci2::Plus,args); break;
case Sub:
args[1] = foci->mk_op(foci2::Times,foci->mk_int("-1"),args[1]); res = foci->mk_op(foci2::Plus,args); break;
case Uminus:
res = foci->mk_op(foci2::Times,foci->mk_int("-1"),args[0]); break;
case Times:
res = foci->mk_op(foci2::Times,args); break;
case Idiv:
res = foci->mk_op(foci2::Div,args); break;
case Mod:
res = foci->mk_app(mod_op,args); break;
case Select:
res = foci->mk_app(select_op,args); break;
case Store:
res = foci->mk_app(store_op,args); break;
case Distinct:
res = foci->mk_op(foci2::Distinct,args); break;
case Uninterpreted: {
symb func = sym(t);
FuncDeclToSymbol::iterator it = func_decl_to_symbol.find(func);
foci2::symb f = (it == func_decl_to_symbol.end()) ? from_symb(func) : it->second;
if(foci->get_symb(f).substr(0,3) == "lbl" && args.size()==1) // HACK to handle Z3 labels
res = args[0];
else if(foci->get_symb(f).substr(0,3) == "lbl" && args.size()==0) // HACK to handle Z3 labels
res = foci->mk_true();
else res = foci->mk_app(f,args);
break;
}
case Numeral: {
std::string s = string_of_numeral(t);
res = foci->mk_int(s);
break;
}
case Forall:
case Exists: {
bool is_forall = op(t) == Forall;
foci2::ops qop = is_forall ? foci2::Forall : foci2::Exists;
int bvs = get_quantifier_num_bound(t);
std::vector<int> foci_bvs(bvs);
for(int i = 0; i < bvs; i++){
std::string name = get_quantifier_bound_name(t,i);
//Z3_string name = Z3_get_symbol_string(ctx,sym);
// type ty = get_quantifier_bound_type(t,i);
foci2::symb f = foci->mk_func(name);
foci2::ast v = foci->mk_app(f,std::vector<foci2::ast>());
foci_bvs[i] = v;
}
foci2::ast body = from_Z3_ast(get_quantifier_body(t));
foci_bvs.push_back(body);
res = foci->mk_op(qop,foci_bvs);
node_to_ast[res] = t; // desperate
break;
}
case Variable: { // a deBruijn index
int index = get_variable_index_value(t);
type ty = get_type(t);
foci2::symb symbol = make_deBruijn_symbol(index,ty);
res = foci->mk_app(symbol,std::vector<foci2::ast>());
break;
}
default:
{
std::cerr << "iZ3: unsupported Z3 operator in expression\n ";
print_expr(std::cerr,t);
std::cerr << "\n";
SASSERT(0 && "iZ3: unsupported Z3 operator");
}
}
return res;
}
// convert an expr to Z3 ast
ast to_Z3_ast(foci2::ast i){
std::pair<foci2::ast,ast> foo(i,ast());
std::pair<NodeToAst::iterator,bool> bar = node_to_ast.insert(foo);
if(!bar.second) return bar.first->second;
ast &res = bar.first->second;
if(i < 0){
res = mk_not(to_Z3_ast(-i));
return res;
}
// get the arguments
unsigned n = foci->get_num_args(i);
std::vector<ast> args(n);
for(unsigned j = 0; j < n; j++)
args[j] = to_Z3_ast(foci->get_arg(i,j));
// handle operators
foci2::ops o;
foci2::symb f;
std::string nval;
if(foci->get_true(i))
res = mk_true();
else if(foci->get_false(i))
res = mk_false();
else if(foci->get_op(i,o)){
switch(o){
case foci2::And:
res = make(And,args); break;
case foci2::Or:
res = make(Or,args); break;
case foci2::Not:
res = mk_not(args[0]); break;
case foci2::Iff:
res = make(Iff,args[0],args[1]); break;
case foci2::Ite:
res = make(Ite,args[0],args[1],args[2]); break;
case foci2::Equal:
res = make(Equal,args[0],args[1]); break;
case foci2::Plus:
res = make(Plus,args); break;
case foci2::Times:
res = make(Times,args); break;
case foci2::Div:
res = make(Idiv,args[0],args[1]); break;
case foci2::Leq:
res = make(Leq,args[0],args[1]); break;
case foci2::Distinct:
res = make(Distinct,args);
break;
case foci2::Tsym:
res = mk_true();
break;
case foci2::Fsym:
res = mk_false();
break;
case foci2::Forall:
case foci2::Exists:
{
int nargs = n;
std::vector<ast> bounds(nargs-1);
for(int i = 0; i < nargs-1; i++)
bounds[i] = args[i];
opr oz = o == foci2::Forall ? Forall : Exists;
res = make_quant(oz,bounds,args[nargs-1]);
}
break;
default:
SASSERT(false && "unknown built-in op");
}
}
else if(foci->get_int(i,nval)){
res = make_int(nval);
}
else if(foci->get_func(i,f)){
if(f == select_op){
SASSERT(n == 2);
res = make(Select,args[0],args[1]);
}
else if(f == store_op){
SASSERT(n == 3);
res = make(Store,args[0],args[1],args[2]);
}
else if(f == mod_op){
SASSERT(n == 2);
res = make(Mod,args[0],args[1]);
}
else {
std::pair<int,symb> foo(f,(symb)0);
std::pair<SymbolToFuncDecl::iterator,bool> bar = symbol_to_func_decl.insert(foo);
symb &func_decl = bar.first->second;
if(bar.second){
std::cout << "unknown function symbol:\n";
foci->show_ast(i);
SASSERT(0);
}
res = make(func_decl,args);
}
}
else {
std::cerr << "iZ3: unknown FOCI expression kind\n";
SASSERT(0 && "iZ3: unknown FOCI expression kind");
}
return res;
}
int interpolate(const std::vector<ast> &cnsts, std::vector<ast> &itps){
SASSERT((int)cnsts.size() == frames);
std::string lia("lia");
#ifdef _FOCI2
foci = foci2::create(lia);
#else
foci = 0;
#endif
if(!foci){
std::cerr << "iZ3: cannot find foci lia solver.\n";
SASSERT(0);
}
select_op = foci->mk_func("select");
store_op = foci->mk_func("store");
mod_op = foci->mk_func("mod");
std::vector<foci2::ast> foci_cnsts(frames), foci_itps(frames-1), foci_parents;
if(parents)
foci_parents.resize(frames);
for(int i = 0; i < frames; i++){
foci_cnsts[i] = from_Z3_ast(cnsts[i]);
if(parents)
foci_parents[i] = parents[i];
}
int res = foci->interpolate(foci_cnsts, foci_itps, foci_parents);
if(res == 0){
SASSERT((int)foci_itps.size() == frames-1);
itps.resize(frames-1);
for(int i = 0; i < frames-1; i++){
// foci->show_ast(foci_itps[i]);
itps[i] = to_Z3_ast(foci_itps[i]);
}
}
ast_to_node.clear();
node_to_ast.clear();
func_decl_to_symbol.clear();
symbol_to_func_decl.clear();
delete foci;
return res;
}
};
iz3secondary *iz3foci::create(iz3mgr *mgr, int num, int *parents){
return new iz3foci_impl(mgr,num,parents);
}

32
src/interp/iz3foci.h
View file

@ -1,32 +0,0 @@
/*++
Copyright (c) 2011 Microsoft Corporation
Module Name:
iz3foci.h
Abstract:
Implements a secondary solver using foci2.
Author:
Ken McMillan (kenmcmil)
Revision History:
--*/
#ifndef IZ3FOCI_H
#define IZ3FOCI_H
#include "iz3secondary.h"
/** Secondary prover based on Cadence FOCI. */
class iz3foci {
public:
static iz3secondary *create(iz3mgr *mgr, int num, int *parents);
};
#endif

30
src/interp/iz3interp.cpp
View file

@ -35,7 +35,6 @@
#include "iz3profiling.h"
#include "iz3translate.h"
#include "iz3foci.h"
#include "iz3proof.h"
#include "iz3hash.h"
#include "iz3interp.h"
@ -167,22 +166,6 @@ struct frame_reducer {
#endif
#if 0
static lbool test_secondary(context ctx,
int num,
ast *cnsts,
ast *interps,
int *parents = 0
){
iz3secondary *sp = iz3foci::create(ctx,num,parents);
std::vector<ast> frames(num), interpolants(num-1);
std::copy(cnsts,cnsts+num,frames.begin());
int res = sp->interpolate(frames,interpolants);
if(res == 0)
std::copy(interpolants.begin(),interpolants.end(),interps);
return res ? L_TRUE : L_FALSE;
}
#endif
template<class T>
struct killme {
@ -213,11 +196,7 @@ public:
const std::vector<int> &parents,
std::vector<ast> &interps
){
int num = cnsts.size();
iz3secondary *sp = iz3foci::create(this,num,(int *)(parents.empty()?0:&parents[0]));
int res = sp->interpolate(cnsts, interps);
if(res != 0)
throw iz3_exception("secondary failed");
throw iz3_exception("secondary interpolating prover not supported");
}
void proof_to_interpolant(z3pf proof,
@ -248,10 +227,9 @@ public:
if(is_linear(parents_vec))
parents_vec.clear();
// create a secondary prover
iz3secondary *sp = iz3foci::create(this,num,parents_vec.empty()?0:&parents_vec[0]);
sp_killer.set(sp); // kill this on exit
// secondary prover no longer supported
iz3secondary *sp = NULL;
#define BINARY_INTERPOLATION
#ifndef BINARY_INTERPOLATION
// create a translator

6
src/interp/iz3translate.h
View file

@ -53,12 +53,8 @@ class iz3translation : public iz3base {
: iz3base(mgr,_cnsts,_parents,_theory) {}
};
//#define IZ3_TRANSLATE_DIRECT2
#ifdef _FOCI2
#define IZ3_TRANSLATE_DIRECT
#else
// To use a secondary prover, define IZ3_TRANSLATE_DIRECT instead of this
#define IZ3_TRANSLATE_FULL
#endif
#endif

8
src/interp/iz3translate_direct.cpp
View file

@ -332,7 +332,7 @@ public:
}
}
// get the lits of a Z3 clause as foci terms
// get the lits of a Z3 clause as secondary prover terms
void get_Z3_lits(ast t, std::vector<ast> &lits){
opr dk = op(t);
if(dk == False)
@ -666,9 +666,9 @@ public:
#endif
// interpolate using secondary prover
profiling::timer_start("foci");
profiling::timer_start("secondary prover");
int sat = secondary->interpolate(preds,itps);
profiling::timer_stop("foci");
profiling::timer_stop("secondary prover");
std::cout << "lemma done" << std::endl;
@ -1495,7 +1495,7 @@ public:
return find_nll(new_proofs);
}
// translate a Z3 proof term into a foci proof term
// translate a Z3 proof term into a secondary prover proof term
Iproof::node translate_main(ast proof, non_local_lits *nll, bool expect_clause = true){
non_local_lits *old_nll = nll;

12
src/opt/opt_solver.cpp
View file

@ -26,6 +26,7 @@ Notes:
#include "theory_diff_logic.h"
#include "theory_dense_diff_logic.h"
#include "theory_pb.h"
#include "theory_lra.h"
#include "ast_pp.h"
#include "ast_smt_pp.h"
#include "pp_params.hpp"
@ -143,6 +144,9 @@ namespace opt {
else if (typeid(smt::theory_dense_si&) == typeid(*arith_theory)) {
return dynamic_cast<smt::theory_dense_si&>(*arith_theory);
}
else if (typeid(smt::theory_lra&) == typeid(*arith_theory)) {
return dynamic_cast<smt::theory_lra&>(*arith_theory);
}
else {
UNREACHABLE();
return dynamic_cast<smt::theory_mi_arith&>(*arith_theory);
@ -401,6 +405,14 @@ namespace opt {
return th.mk_ge(m_fm, v, val);
}
if (typeid(smt::theory_lra) == typeid(opt)) {
smt::theory_lra& th = dynamic_cast<smt::theory_lra&>(opt);
SASSERT(val.is_finite());
return th.mk_ge(m_fm, v, val.get_numeral());
}
// difference logic?
if (typeid(smt::theory_dense_si) == typeid(opt) &&
val.get_infinitesimal().is_zero()) {
smt::theory_dense_si& th = dynamic_cast<smt::theory_dense_si&>(opt);

244
src/sat/card_extension.cpp
View file

@ -46,6 +46,8 @@ namespace sat {
m_lit(lit),
m_k(k),
m_size(wlits.size()),
m_slack(0),
m_num_watch(0),
m_max_sum(0) {
for (unsigned i = 0; i < wlits.size(); ++i) {
m_wlits[i] = wlits[i];
@ -237,7 +239,6 @@ namespace sat {
p.negate();
}
TRACE("sat", display(tout << "init watch: ", p, true););
SASSERT(p.lit() == null_literal || value(p.lit()) == l_true);
unsigned sz = p.size(), bound = p.k();
@ -249,7 +250,7 @@ namespace sat {
p.swap(i, j);
}
if (slack < bound) {
slack += p[i].first;
slack += p[j].first;
++num_watch;
}
++j;
@ -266,7 +267,7 @@ namespace sat {
if (slack < bound) {
literal lit = p[j].second;
SASSERT(value(p[j].second) == l_false);
for (unsigned i = j + 1; j < sz; ++i) {
for (unsigned i = j + 1; i < sz; ++i) {
if (lvl(lit) < lvl(p[i].second)) {
lit = p[i].second;
}
@ -280,6 +281,7 @@ namespace sat {
p.set_slack(slack);
p.set_num_watch(num_watch);
}
TRACE("sat", display(tout << "init watch: ", p, true););
}
/*
@ -296,7 +298,7 @@ namespace sat {
Lw = Lw u {l_s}
Lu = Lu \ {l_s}
}
if (Sw < bound) conflict
if (Sw < k) conflict
while (Sw < k + a_max) {
assign (l_max)
a_max = max { ai | li in Lw, li = undef }
@ -308,7 +310,19 @@ namespace sat {
*/
void card_extension::add_index(pb& p, unsigned index, literal lit) {
if (value(lit) == l_undef) {
m_pb_undef.push_back(index);
if (p[index].first > m_a_max) {
m_a_max = p[index].first;
}
}
}
lbool card_extension::add_assign(pb& p, literal alit) {
TRACE("sat", display(tout << "assign: " << alit << "\n", p, true););
SASSERT(!s().inconsistent());
unsigned sz = p.size();
unsigned bound = p.k();
unsigned num_watch = p.num_watch();
@ -316,63 +330,43 @@ namespace sat {
SASSERT(value(alit) == l_false);
SASSERT(p.lit() == null_literal || value(p.lit()) == l_true);
SASSERT(num_watch <= sz);
SASSERT(num_watch > 0);
unsigned index = 0;
unsigned a_max = 0;
unsigned max_index = 0;
m_a_max = 0;
m_pb_undef.reset();
for (; index < num_watch; ++index) {
literal lit = p[index].second;
if (lit == alit) {
break;
}
if (value(lit) == l_undef) {
m_pb_undef.push_back(index);
if (p[index].first > a_max) {
a_max = p[index].first;
max_index = index;
}
}
add_index(p, index, lit);
}
for (unsigned j = index + 1; a_max == 0 && j < num_watch; ++j) {
literal lit = p[j].second;
if (value(lit) == l_undef) {
m_pb_undef.push_back(j);
a_max = p[j].first;
max_index = j;
}
}
for (unsigned j = num_watch; a_max == 0 && j < sz; ++j) {
literal lit = p[j].second;
if (value(lit) == l_undef) {
p.swap(j, num_watch);
m_pb_undef.push_back(num_watch);
a_max = p[num_watch].first;
max_index = num_watch;
}
SASSERT(index < num_watch);
unsigned index1 = index + 1;
for (; m_a_max == 0 && index1 < num_watch; ++index1) {
add_index(p, index1, p[index1].second);
}
unsigned val = p[index].first;
SASSERT(num_watch > 0);
SASSERT(index < num_watch);
SASSERT(value(p[index].second) == l_false);
SASSERT(val <= slack);
slack -= val;
// find literals to swap with:
for (unsigned j = num_watch; j < sz && slack < bound + a_max; ++j) {
if (value(p[j].second) != l_false) {
for (unsigned j = num_watch; j < sz && slack < bound + m_a_max; ++j) {
literal lit = p[j].second;
if (value(lit) != l_false) {
slack += p[j].first;
watch_literal(p, p[j]);
p.swap(num_watch, j);
if (value(p[num_watch].second) == l_undef && a_max < p[num_watch].first) {
m_pb_undef.push_back(num_watch);
a_max = p[num_watch].first;
max_index = num_watch;
}
add_index(p, num_watch, lit);
++num_watch;
}
}
SASSERT(!s().inconsistent());
DEBUG_CODE(for (auto idx : m_pb_undef) { SASSERT(value(p[idx].second) == l_undef); });
if (slack < bound) {
// maintain watching the literal
slack += val;
@ -385,34 +379,38 @@ namespace sat {
}
// swap out the watched literal.
p.set_slack(slack);
--num_watch;
SASSERT(num_watch > 0);
p.set_slack(slack);
p.set_num_watch(num_watch);
p.swap(num_watch, index);
if (num_watch == max_index) {
max_index = index;
}
SASSERT(max_index < sz);
while (slack < bound + a_max && !s().inconsistent()) {
// variable at max-index must be assigned to true.
assign(p, p[max_index].second);
a_max = 0;
// find the next a_max among m_pb_undef
while (!m_pb_undef.empty() && l_undef != value(p[m_pb_undef.back()].second)) {
if (slack < bound + m_a_max) {
TRACE("sat", display(tout, p, false); for(auto j : m_pb_undef) tout << j << "\n";);
literal_vector to_assign;
while (!m_pb_undef.empty()) {
index1 = m_pb_undef.back();
literal lit = p[index1].second;
TRACE("sat", tout << index1 << " " << lit << "\n";);
if (slack >= bound + p[index1].first) {
break;
}
m_pb_undef.pop_back();
to_assign.push_back(lit);
}
if (m_pb_undef.empty()) {
break;
for (literal lit : to_assign) {
if (s().inconsistent()) break;
if (value(lit) == l_undef) {
assign(p, lit);
}
}
max_index = m_pb_undef.back();
a_max = p[max_index].first;
m_pb_undef.pop_back();
}
return s().inconsistent() ? l_false : l_true;
TRACE("sat", display(tout << "assign: " << alit << "\n", p, true););
return l_undef;
}
void card_extension::watch_literal(pb& p, wliteral l) {
@ -467,6 +465,7 @@ namespace sat {
set_conflict(p, lit);
break;
default:
SASSERT(validate_unit_propagation(p, lit));
m_stats.m_num_pb_propagations++;
m_num_propagations_since_pop++;
if (s().m_config.m_drat) {
@ -483,7 +482,8 @@ namespace sat {
}
void card_extension::display(std::ostream& out, pb const& p, bool values) const {
out << p.lit() << "[" << p.size() << "]";
if (p.lit() != null_literal) out << p.lit();
out << "[watch: " << p.num_watch() << ", slack: " << p.slack() << "]";
if (p.lit() != null_literal && values) {
out << "@(" << value(p.lit());
if (value(p.lit()) != l_undef) {
@ -514,7 +514,7 @@ namespace sat {
}
void card_extension::asserted_pb(literal l, ptr_vector<pb>* pbs, pb* p0) {
TRACE("sat", tout << l << " " << !is_tag_empty(pbs) << " " << (p0 != 0) << "\n";);
TRACE("sat", tout << "literal: " << l << " has pb: " << !is_tag_empty(pbs) << " p0 != 0: " << (p0 != 0) << "\n";);
if (!is_tag_empty(pbs)) {
ptr_vector<pb>::iterator begin = pbs->begin();
ptr_vector<pb>::iterator it = begin, it2 = it, end = pbs->end();
@ -524,20 +524,28 @@ namespace sat {
continue;
}
switch (add_assign(p, ~l)) {
case l_false: // conflict
for (; it != end; ++it, ++it2) {
*it2 = *it;
}
SASSERT(s().inconsistent());
pbs->set_end(it2);
return;
case l_true: // unit propagation, keep watching the literal
if (it2 != it) {
*it2 = *it;
}
++it2;
break;
case l_false: // conflict.
SASSERT(s().inconsistent());
for (; it != end; ++it, ++it2) {
*it2 = *it;
}
pbs->set_end(it2);
return;
case l_undef: // watch literal was swapped
if (s().inconsistent()) {
++it;
for (; it != end; ++it, ++it2) {
*it2 = *it;
}
pbs->set_end(it2);
return;
}
break;
}
}
@ -844,6 +852,7 @@ namespace sat {
literal consequent = s().m_not_l;
justification js = s().m_conflict;
TRACE("sat", tout << consequent << " " << js << "\n";);
TRACE("sat", s().display(tout););
m_conflict_lvl = s().get_max_lvl(consequent, js);
if (consequent != null_literal) {
consequent.neg();
@ -942,7 +951,7 @@ namespace sat {
m_bound += offset;
inc_coeff(consequent, offset);
get_pb_antecedents(consequent, p, m_lemma);
TRACE("sat", tout << m_lemma << "\n";);
TRACE("sat", display(tout, p, true); tout << m_lemma << "\n";);
for (unsigned i = 0; i < m_lemma.size(); ++i) {
process_antecedent(~m_lemma[i], offset);
}
@ -989,6 +998,7 @@ namespace sat {
SASSERT(lvl(v) == m_conflict_lvl);
s().reset_mark(v);
--idx;
TRACE("sat", tout << "Unmark: v" << v << "\n";);
--m_num_marks;
js = s().m_justification[v];
offset = get_abs_coeff(v);
@ -1015,7 +1025,6 @@ namespace sat {
m_lemma.reset();
#if 1
m_lemma.push_back(null_literal);
for (unsigned i = 0; 0 <= slack && i < m_active_vars.size(); ++i) {
bool_var v = m_active_vars[i];
@ -1038,6 +1047,7 @@ namespace sat {
}
}
#if 0
if (jus.size() > 1) {
std::cout << jus.size() << "\n";
for (unsigned i = 0; i < jus.size(); ++i) {
@ -1047,6 +1057,7 @@ namespace sat {
active2pb(m_A);
display(std::cout, m_A);
}
#endif
if (slack >= 0) {
@ -1057,7 +1068,7 @@ namespace sat {
if (m_lemma[0] == null_literal) {
m_lemma[0] = m_lemma.back();
m_lemma.pop_back();
unsigned level = lvl(m_lemma[0]);
unsigned level = m_lemma.empty() ? 0 : lvl(m_lemma[0]);
for (unsigned i = 1; i < m_lemma.size(); ++i) {
if (lvl(m_lemma[i]) > level) {
level = lvl(m_lemma[i]);
@ -1066,37 +1077,6 @@ namespace sat {
}
IF_VERBOSE(2, verbose_stream() << "(sat.card set level to " << level << " < " << m_conflict_lvl << ")\n";);
}
#else
++idx;
while (0 <= slack) {
literal lit = lits[idx];
bool_var v = lit.var();
if (m_active_var_set.contains(v)) {
int coeff = get_coeff(v);
if (coeff < 0 && !lit.sign()) {
slack += coeff;
m_lemma.push_back(~lit);
}
else if (coeff > 0 && lit.sign()) {
slack -= coeff;
m_lemma.push_back(~lit);
}
}
if (idx == 0 && slack >= 0) {
IF_VERBOSE(2, verbose_stream() << "(sat.card non-asserting)\n";);
goto bail_out;
}
SASSERT(idx > 0 || slack < 0);
--idx;
}
if (m_lemma.size() >= 2 && lvl(m_lemma[1]) == m_conflict_lvl) {
// TRACE("sat", tout << "Bail out on no progress " << lit << "\n";);
IF_VERBOSE(2, verbose_stream() << "(sat.card bail non-asserting resolvent)\n";);
goto bail_out;
}
#endif
SASSERT(slack < 0);
@ -1108,7 +1088,7 @@ namespace sat {
svector<drat::premise> ps; // TBD fill in
s().m_drat.add(m_lemma, ps);
}
s().m_lemma.reset();
s().m_lemma.append(m_lemma);
for (unsigned i = 1; i < m_lemma.size(); ++i) {
@ -1153,6 +1133,7 @@ namespace sat {
if (level > 0 && !s().is_marked(v) && level == m_conflict_lvl) {
s().mark(v);
TRACE("sat", tout << "Mark: v" << v << "\n";);
++m_num_marks;
}
inc_coeff(l, offset);
@ -1347,20 +1328,30 @@ namespace sat {
void card_extension::get_pb_antecedents(literal l, pb const& p, literal_vector& r) {
if (p.lit() != null_literal) r.push_back(p.lit());
SASSERT(p.lit() == null_literal || value(p.lit()) == l_true);
unsigned k = p.k();
unsigned max_sum = p.max_sum();
for (unsigned i = p.size(); i > 0 && max_sum >= k; ) {
--i;
literal lit = p[i].second;
if (lit == l) {
max_sum -= p[i].first;
}
else if (value(lit) == l_false) {
r.push_back(~p[i].second);
max_sum -= p[i].first;
TRACE("sat", display(tout, p, true););
// unsigned coeff = get_coeff(p, l);
unsigned coeff = 0;
for (unsigned j = 0; j < p.num_watch(); ++j) {
if (p[j].second == l) {
coeff = p[j].first;
break;
}
}
SASSERT(max_sum < k);
unsigned slack = p.slack() - coeff;
unsigned i = p.num_watch();
// skip entries that are not required for unit propagation.
// slack - coeff + w_head < k
unsigned h = 0;
for (; i < p.size() && p[i].first + h + slack < p.k(); ++i) {
h += p[i].first;
}
for (; i < p.size(); ++i) {
literal lit = p[i].second;
SASSERT(l_false == value(lit));
r.push_back(~lit);
}
}
void card_extension::get_card_antecedents(literal l, card const& c, literal_vector& r) {
@ -1720,7 +1711,9 @@ namespace sat {
std::ostream& card_extension::display_justification(std::ostream& out, ext_justification_idx idx) const {
if (is_card_index(idx)) {
card& c = index2card(idx);
out << "bound " << c.lit() << ": ";
out << "bound ";
if (c.lit() != null_literal) out << c.lit();
out << ": ";
for (unsigned i = 0; i < c.size(); ++i) {
out << c[i] << " ";
}
@ -1728,14 +1721,18 @@ namespace sat {
}
else if (is_xor_index(idx)) {
xor& x = index2xor(idx);
out << "xor " << x.lit() << ": ";
out << "xor ";
if (x.lit() != null_literal) out << x.lit();
out << ": ";
for (unsigned i = 0; i < x.size(); ++i) {
out << x[i] << " ";
}
}
else if (is_pb_index(idx)) {
pb& p = index2pb(idx);
out << "pb " << p.lit() << ": ";
out << "pb ";
if (p.lit() != null_literal) out << p.lit();
out << ": ";
for (unsigned i = 0; i < p.size(); ++i) {
out << p[i].first << "*" << p[i].second << " ";
}
@ -1776,6 +1773,21 @@ namespace sat {
}
return true;
}
bool card_extension::validate_unit_propagation(pb const& p, literal alit) {
if (p.lit() != null_literal && value(p.lit()) != l_true) return false;
unsigned sum = 0;
TRACE("sat", display(tout << "validate: " << alit << "\n", p, true););
for (unsigned i = 0; i < p.size(); ++i) {
literal lit = p[i].second;
lbool val = value(lit);
if (val != l_false && lit != alit) {
sum += p[i].first;
}
}
return sum < p.k();
}
bool card_extension::validate_lemma() {
int val = -m_bound;
normalize_active_coeffs();

3
src/sat/card_extension.h
View file

@ -232,10 +232,12 @@ namespace sat {
// pb functionality
unsigned m_a_max;
void copy_pb(card_extension& result);
void asserted_pb(literal l, ptr_vector<pb>* pbs, pb* p);
void init_watch(pb& p, bool is_true);
lbool add_assign(pb& p, literal alit);
void add_index(pb& p, unsigned index, literal lit);
void watch_literal(pb& p, wliteral lit);
void clear_watch(pb& p);
void set_conflict(pb& p, literal lit);
@ -280,6 +282,7 @@ namespace sat {
bool validate_assign(literal_vector const& lits, literal lit);
bool validate_lemma();
bool validate_unit_propagation(card const& c);
bool validate_unit_propagation(pb const& p, literal lit);
bool validate_conflict(literal_vector const& lits, ineq& p);
ineq m_A, m_B, m_C;

35
src/sat/sat_solver.cpp
View file

@ -1938,23 +1938,25 @@ namespace sat {
void solver::learn_lemma_and_backjump() {
TRACE("sat_lemma", tout << "new lemma size: " << m_lemma.size() << "\n" << m_lemma << "\n";);
if (m_config.m_minimize_lemmas) {
minimize_lemma();
reset_lemma_var_marks();
if (m_config.m_dyn_sub_res)
dyn_sub_res();
TRACE("sat_lemma", tout << "new lemma (after minimization) size: " << m_lemma.size() << "\n" << m_lemma << "\n";);
}
else
reset_lemma_var_marks();
literal_vector::iterator it = m_lemma.begin();
literal_vector::iterator end = m_lemma.end();
unsigned new_scope_lvl = 0;
++it;
for(; it != end; ++it) {
bool_var var = (*it).var();
new_scope_lvl = std::max(new_scope_lvl, lvl(var));
if (!m_lemma.empty()) {
if (m_config.m_minimize_lemmas) {
minimize_lemma();
reset_lemma_var_marks();
if (m_config.m_dyn_sub_res)
dyn_sub_res();
TRACE("sat_lemma", tout << "new lemma (after minimization) size: " << m_lemma.size() << "\n" << m_lemma << "\n";);
}
else
reset_lemma_var_marks();
literal_vector::iterator it = m_lemma.begin();
literal_vector::iterator end = m_lemma.end();
++it;
for(; it != end; ++it) {
bool_var var = (*it).var();
new_scope_lvl = std::max(new_scope_lvl, lvl(var));
}
}
unsigned glue = num_diff_levels(m_lemma.size(), m_lemma.c_ptr());
@ -2390,6 +2392,7 @@ namespace sat {
assigned at level 0.
*/
void solver::minimize_lemma() {
SASSERT(!m_lemma.empty());
SASSERT(m_unmark.empty());
//m_unmark.reset();
updt_lemma_lvl_set();

110
src/shell/lp_frontend.cpp Normal file
View file

@ -0,0 +1,110 @@
/*++
Copyright (c) 2016 Microsoft Corporation
Author:
Lev Nachmanson 2016-10-27
--*/
#include "lp_params.hpp"
#include "util/lp/lp_settings.h"
#include "util/lp/mps_reader.h"
#include "timeout.h"
#include "cancel_eh.h"
#include "scoped_timer.h"
#include "rlimit.h"
#include "gparams.h"
#include <signal.h>
static lean::lp_solver<double, double>* g_solver = 0;
static void display_statistics() {
if (g_solver && g_solver->settings().print_statistics) {
// TBD display relevant information about statistics
}
}
static void STD_CALL on_ctrl_c(int) {
signal (SIGINT, SIG_DFL);
#pragma omp critical (g_display_stats)
{
display_statistics();
}
raise(SIGINT);
}
static void on_timeout() {
#pragma omp critical (g_display_stats)
{
display_statistics();
exit(0);
}
}
struct front_end_resource_limit : public lean::lp_resource_limit {
reslimit& m_reslim;
front_end_resource_limit(reslimit& lim):
m_reslim(lim)
{}
virtual bool get_cancel_flag() { return !m_reslim.inc(); }
};
void run_solver(lp_params & params, char const * mps_file_name) {
reslimit rlim;
unsigned timeout = gparams::get().get_uint("timeout", 0);
unsigned rlimit = gparams::get().get_uint("rlimit", 0);
front_end_resource_limit lp_limit(rlim);
scoped_rlimit _rlimit(rlim, rlimit);
cancel_eh<reslimit> eh(rlim);
scoped_timer timer(timeout, &eh);
std::string fn(mps_file_name);
lean::mps_reader<double, double> reader(fn);
reader.set_message_stream(&std::cout); // can be redirected
reader.read();
if (!reader.is_ok()) {
std::cerr << "cannot process " << mps_file_name << std::endl;
return;
}
lean::lp_solver<double, double> * solver = reader.create_solver(false); // false - to create the primal solver
solver->settings().set_resource_limit(lp_limit);
g_solver = solver;
if (params.min()) {
solver->flip_costs();
}
solver->settings().set_message_ostream(&std::cout);
solver->settings().report_frequency = params.rep_freq();
solver->settings().print_statistics = params.print_stats();
solver->find_maximal_solution();
*(solver->settings().get_message_ostream()) << "status is " << lp_status_to_string(solver->get_status()) << std::endl;
if (solver->get_status() == lean::OPTIMAL) {
if (params.min()) {
solver->flip_costs();
}
solver->print_model(std::cout);
}
// #pragma omp critical (g_display_stats)
{
display_statistics();
register_on_timeout_proc(0);
g_solver = 0;
}
delete solver;
}
unsigned read_mps_file(char const * mps_file_name) {
signal(SIGINT, on_ctrl_c);
register_on_timeout_proc(on_timeout);
lp_params p;
param_descrs r;
p.collect_param_descrs(r);
run_solver(p, mps_file_name);
return 0;
}

7
src/shell/lp_frontend.h Normal file
View file

@ -0,0 +1,7 @@
/*
Copyright (c) 2013 Microsoft Corporation. All rights reserved.
Author: Lev Nachmanson
*/
#pragma once
unsigned read_mps_file(char const * mps_file_name);

10
src/shell/main.cpp
View file

@ -35,8 +35,9 @@ Revision History:
#include"error_codes.h"
#include"gparams.h"
#include"env_params.h"
#include "lp_frontend.h"
typedef enum { IN_UNSPECIFIED, IN_SMTLIB, IN_SMTLIB_2, IN_DATALOG, IN_DIMACS, IN_WCNF, IN_OPB, IN_Z3_LOG } input_kind;
typedef enum { IN_UNSPECIFIED, IN_SMTLIB, IN_SMTLIB_2, IN_DATALOG, IN_DIMACS, IN_WCNF, IN_OPB, IN_Z3_LOG, IN_MPS } input_kind;
std::string g_aux_input_file;
char const * g_input_file = 0;
@ -342,6 +343,10 @@ int STD_CALL main(int argc, char ** argv) {
else if (strcmp(ext, "smt") == 0) {
g_input_kind = IN_SMTLIB;
}
else if (strcmp(ext, "mps") == 0 || strcmp(ext, "sif") == 0 ||
strcmp(ext, "MPS") == 0 || strcmp(ext, "SIF") == 0) {
g_input_kind = IN_MPS;
}
}
}
switch (g_input_kind) {
@ -367,6 +372,9 @@ int STD_CALL main(int argc, char ** argv) {
case IN_Z3_LOG:
replay_z3_log(g_input_file);
break;
case IN_MPS:
return_value = read_mps_file(g_input_file);
break;
default:
UNREACHABLE();
}

2
src/smt/smt2_extra_cmds.cpp
View file

@ -44,4 +44,4 @@ public:
void install_smt2_extra_cmds(cmd_context & ctx) {
ctx.insert(alloc(include_cmd));
}
}

21
src/smt/smt_setup.cpp
View file

@ -20,6 +20,7 @@ Revision History:
#include"smt_setup.h"
#include"static_features.h"
#include"theory_arith.h"
#include"theory_lra.h"
#include"theory_dense_diff_logic.h"
#include"theory_diff_logic.h"
#include"theory_utvpi.h"
@ -442,7 +443,7 @@ namespace smt {
m_params.m_arith_propagate_eqs = false;
m_params.m_eliminate_term_ite = true;
m_params.m_nnf_cnf = false;
setup_mi_arith();
setup_r_arith();
}
void setup::setup_QF_LRA(static_features const & st) {
@ -467,6 +468,10 @@ namespace smt {
m_params.m_restart_adaptive = false;
}
m_params.m_arith_small_lemma_size = 32;
setup_r_arith();
}
void setup::setup_QF_LIRA(static_features const& st) {
setup_mi_arith();
}
@ -475,7 +480,7 @@ namespace smt {
m_params.m_relevancy_lvl = 0;
m_params.m_arith_expand_eqs = true;
m_params.m_arith_reflect = false;
m_params.m_arith_propagate_eqs = false;
m_params.m_arith_propagate_eqs = false;
m_params.m_nnf_cnf = false;
setup_i_arith();
}
@ -539,7 +544,7 @@ namespace smt {
m_params.m_relevancy_lvl = 0;
m_params.m_arith_reflect = false;
m_params.m_nnf_cnf = false;
setup_mi_arith();
setup_r_arith();
}
void setup::setup_QF_BV() {
@ -718,6 +723,12 @@ namespace smt {
m_context.register_plugin(alloc(smt::theory_i_arith, m_manager, m_params));
}
void setup::setup_r_arith() {
// to disable theory lra
// m_context.register_plugin(alloc(smt::theory_mi_arith, m_manager, m_params));
m_context.register_plugin(alloc(smt::theory_lra, m_manager, m_params));
}
void setup::setup_mi_arith() {
if (m_params.m_arith_mode == AS_OPTINF) {
m_context.register_plugin(alloc(smt::theory_inf_arith, m_manager, m_params));
@ -929,7 +940,9 @@ namespace smt {
}
if (st.num_theories() == 1 && is_arith(st)) {
if (st.m_has_real)
if ((st.m_has_int && st.m_has_real) || (st.m_num_non_linear != 0))
setup_QF_LIRA(st);
else if (st.m_has_real)
setup_QF_LRA(st);
else
setup_QF_LIA(st);

2
src/smt/smt_setup.h
View file

@ -65,6 +65,7 @@ namespace smt {
void setup_QF_LRA();
void setup_QF_LRA(static_features const & st);
void setup_QF_LIA();
void setup_QF_LIRA(static_features const& st);
void setup_QF_LIA(static_features const & st);
void setup_QF_UFLIA();
void setup_QF_UFLIA(static_features & st);
@ -99,6 +100,7 @@ namespace smt {
void setup_card();
void setup_i_arith();
void setup_mi_arith();
void setup_r_arith();
void setup_fpa();
void setup_str();

6
src/smt/smt_theory.cpp
View file

@ -54,7 +54,7 @@ namespace smt {
}
}
void theory::display_app(std::ostream & out, app * n) const {
std::ostream& theory::display_app(std::ostream & out, app * n) const {
func_decl * d = n->get_decl();
if (n->get_num_args() == 0) {
out << d->get_name();
@ -73,9 +73,10 @@ namespace smt {
else {
out << "#" << n->get_id();
}
return out;
}
void theory::display_flat_app(std::ostream & out, app * n) const {
std::ostream& theory::display_flat_app(std::ostream & out, app * n) const {
func_decl * d = n->get_decl();
if (n->get_num_args() == 0) {
out << d->get_name();
@ -106,6 +107,7 @@ namespace smt {
else {
out << "#" << n->get_id();
}
return out;
}
bool theory::is_relevant_and_shared(enode * n) const {

12
src/smt/smt_theory.h
View file

@ -337,14 +337,14 @@ namespace smt {
virtual void collect_statistics(::statistics & st) const {
}
void display_app(std::ostream & out, app * n) const;
void display_flat_app(std::ostream & out, app * n) const;
void display_var_def(std::ostream & out, theory_var v) const { return display_app(out, get_enode(v)->get_owner()); }
std::ostream& display_app(std::ostream & out, app * n) const;
void display_var_flat_def(std::ostream & out, theory_var v) const { return display_flat_app(out, get_enode(v)->get_owner()); }
std::ostream& display_flat_app(std::ostream & out, app * n) const;
std::ostream& display_var_def(std::ostream & out, theory_var v) const { return display_app(out, get_enode(v)->get_owner()); }
std::ostream& display_var_flat_def(std::ostream & out, theory_var v) const { return display_flat_app(out, get_enode(v)->get_owner()); }
/**
\brief Assume eqs between variable that are equal with respect to the given table.

11
src/smt/tactic/smt_tactic.cpp
View file

@ -21,6 +21,7 @@ Notes:
#include"smt_kernel.h"
#include"smt_params.h"
#include"smt_params_helper.hpp"
#include"lp_params.hpp"
#include"rewriter_types.h"
#include"filter_model_converter.h"
#include"ast_util.h"
@ -64,6 +65,10 @@ public:
return m_params;
}
params_ref & params() {
return m_params_ref;
}
void updt_params_core(params_ref const & p) {
m_candidate_models = p.get_bool("candidate_models", false);
m_fail_if_inconclusive = p.get_bool("fail_if_inconclusive", true);
@ -73,6 +78,7 @@ public:
TRACE("smt_tactic", tout << "updt_params: " << p << "\n";);
updt_params_core(p);
fparams().updt_params(p);
m_params_ref.copy(p);
m_logic = p.get_sym(symbol("logic"), m_logic);
if (m_logic != symbol::null && m_ctx) {
m_ctx->set_logic(m_logic);
@ -84,6 +90,7 @@ public:
r.insert("candidate_models", CPK_BOOL, "(default: false) create candidate models even when quantifier or theory reasoning is incomplete.");
r.insert("fail_if_inconclusive", CPK_BOOL, "(default: true) fail if found unsat (sat) for under (over) approximated goal.");
smt_params_helper::collect_param_descrs(r);
lp_params::collect_param_descrs(r);
}
@ -112,10 +119,12 @@ public:
struct scoped_init_ctx {
smt_tactic & m_owner;
smt_params m_params; // smt-setup overwrites parameters depending on the current assertions.
params_ref m_params_ref;
scoped_init_ctx(smt_tactic & o, ast_manager & m):m_owner(o) {
m_params = o.fparams();
smt::kernel * new_ctx = alloc(smt::kernel, m, m_params);
m_params_ref = o.params();
smt::kernel * new_ctx = alloc(smt::kernel, m, m_params, m_params_ref);
TRACE("smt_tactic", tout << "logic: " << o.m_logic << "\n";);
new_ctx->set_logic(o.m_logic);
if (o.m_callback) {

2615
src/smt/theory_lra.cpp Normal file
View file

File diff suppressed because it is too large Load diff

97
src/smt/theory_lra.h Normal file
View file

@ -0,0 +1,97 @@
/*++
Copyright (c) 2016 Microsoft Corporation
Module Name:
theory_lra.h
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach) 2016-25-3
Nikolaj Bjorner (nbjorner)
Revision History:
--*/
#pragma once
#include "theory_opt.h"
namespace smt {
class theory_lra : public theory, public theory_opt {
class imp;
imp* m_imp;
public:
theory_lra(ast_manager& m, theory_arith_params& ap);
virtual ~theory_lra();
virtual theory* mk_fresh(context* new_ctx);
virtual char const* get_name() const { return "lra"; }
virtual void init(context * ctx);
virtual bool internalize_atom(app * atom, bool gate_ctx);
virtual bool internalize_term(app * term);
virtual void internalize_eq_eh(app * atom, bool_var v);
virtual void assign_eh(bool_var v, bool is_true);
virtual void new_eq_eh(theory_var v1, theory_var v2);
virtual bool use_diseqs() const;
virtual void new_diseq_eh(theory_var v1, theory_var v2);
virtual void push_scope_eh();
virtual void pop_scope_eh(unsigned num_scopes);
virtual void restart_eh();
virtual void relevant_eh(app* e);
virtual void init_search_eh();
virtual final_check_status final_check_eh();
virtual bool is_shared(theory_var v) const;
virtual bool can_propagate();
virtual void propagate();
virtual justification * why_is_diseq(theory_var v1, theory_var v2);
// virtual void flush_eh();
virtual void reset_eh();
virtual void init_model(model_generator & m);
virtual model_value_proc * mk_value(enode * n, model_generator & mg);
virtual bool get_value(enode* n, expr_ref& r);
virtual bool validate_eq_in_model(theory_var v1, theory_var v2, bool is_true) const;
virtual void display(std::ostream & out) const;
virtual void collect_statistics(::statistics & st) const;
// optimization
virtual inf_eps value(theory_var);
virtual inf_eps maximize(theory_var v, expr_ref& blocker, bool& has_shared);
virtual theory_var add_objective(app* term);
virtual expr_ref mk_ge(filter_model_converter& fm, theory_var v, inf_rational const& val);
};
}

2
src/smt/theory_seq.cpp
View file

@ -2176,6 +2176,8 @@ bool theory_seq::simplify_and_solve_eqs() {
return m_new_propagation || ctx.inconsistent();
}
void theory_seq::internalize_eq_eh(app * atom, bool_var v) {}
bool theory_seq::internalize_term(app* term) {
context & ctx = get_context();

2
src/smt/theory_seq.h
View file

@ -343,7 +343,7 @@ namespace smt {
virtual final_check_status final_check_eh();
virtual bool internalize_atom(app* atom, bool) { return internalize_term(atom); }
virtual bool internalize_term(app*);
virtual void internalize_eq_eh(app * atom, bool_var v) {}
virtual void internalize_eq_eh(app * atom, bool_var v);
virtual void new_eq_eh(theory_var, theory_var);
virtual void new_diseq_eh(theory_var, theory_var);
virtual void assign_eh(bool_var v, bool is_true);

46
src/smt/theory_str.cpp
View file

@ -60,6 +60,7 @@ namespace smt {
totalCacheAccessCount(0),
cacheHitCount(0),
cacheMissCount(0),
m_fresh_id(0),
m_find(*this),
m_trail_stack(*this)
{
@ -441,7 +442,12 @@ namespace smt {
}
app * theory_str::mk_fresh_const(char const* name, sort* s) {
return u.mk_skolem(symbol(name), 0, 0, s);
string_buffer<64> buffer;
buffer << name;
buffer << "!tmp";
buffer << m_fresh_id;
m_fresh_id++;
return u.mk_skolem(symbol(buffer.c_str()), 0, 0, s);
}
@ -5552,8 +5558,8 @@ namespace smt {
for (; arg1_grdItor != groundedMap[arg1DeAlias].end(); arg1_grdItor++) {
std::vector<expr*> ndVec;
ndVec.insert(ndVec.end(), arg0_grdItor->first.begin(), arg0_grdItor->first.end());
int arg0VecSize = arg0_grdItor->first.size();
int arg1VecSize = arg1_grdItor->first.size();
size_t arg0VecSize = arg0_grdItor->first.size();
size_t arg1VecSize = arg1_grdItor->first.size();
if (arg0VecSize > 0 && arg1VecSize > 0 && u.str.is_string(arg0_grdItor->first[arg0VecSize - 1]) && u.str.is_string(arg1_grdItor->first[0])) {
ndVec.pop_back();
ndVec.push_back(mk_concat(arg0_grdItor->first[arg0VecSize - 1], arg1_grdItor->first[0]));
@ -5645,8 +5651,8 @@ namespace smt {
}
bool theory_str::is_partial_in_grounded_concat(const std::vector<expr*> & strVec, const std::vector<expr*> & subStrVec) {
int strCnt = strVec.size();
int subStrCnt = subStrVec.size();
size_t strCnt = strVec.size();
size_t subStrCnt = subStrVec.size();
if (strCnt == 0 || subStrCnt == 0) {
return false;
@ -5717,7 +5723,7 @@ namespace smt {
}
// tail nodes
int tailIdx = i + subStrCnt - 1;
size_t tailIdx = i + subStrCnt - 1;
zstring subStrTailVal;
if (u.str.is_string(subStrVec[subStrCnt - 1], subStrTailVal)) {
zstring strTailVal;
@ -5908,8 +5914,14 @@ namespace smt {
app * n2_curr = to_app(n2);
// case 0: n1_curr is const string, n2_curr is const string
if (u.str.is_string(n1_curr) && u.str.is_string(n2_curr)) {
if (n1_curr != n2_curr) {
zstring n1_curr_str, n2_curr_str;
if (u.str.is_string(n1_curr, n1_curr_str) && u.str.is_string(n2_curr, n2_curr_str)) {
TRACE("str", tout << "checking string constants: n1=" << n1_curr_str << ", n2=" << n2_curr_str << std::endl;);
if (n1_curr_str == n2_curr_str) {
// TODO(mtrberzi) potential correction: if n1_curr != n2_curr,
// assert that these two terms are in fact equal, because they ought to be
return true;
} else {
return false;
}
}
@ -6864,7 +6876,7 @@ namespace smt {
// easiest case. we will search within these bounds
} else if (upper_bound_exists && !lower_bound_exists) {
// search between 0 and the upper bound
v_lower_bound == rational::zero();
v_lower_bound = rational::zero();
} else if (lower_bound_exists && !upper_bound_exists) {
// check some finite portion of the search space
v_upper_bound = v_lower_bound + rational(10);
@ -8251,7 +8263,7 @@ namespace smt {
}
}
if (lId == -1)
lId = mLMap.size();
lId = static_cast<int>(mLMap.size());
for (std::set<expr*>::iterator itor2 = nSet.begin(); itor2 != nSet.end(); itor2++) {
bool itorHasEqcValue = false;
get_eqc_value(*itor2, itorHasEqcValue);
@ -8290,7 +8302,7 @@ namespace smt {
}
}
if (rId == -1)
rId = mRMap.size();
rId = static_cast<int>(mRMap.size());
for (itor2 = nSet.begin(); itor2 != nSet.end(); itor2++) {
bool rHasEqcValue = false;
get_eqc_value(*itor2, rHasEqcValue);
@ -9182,7 +9194,7 @@ namespace smt {
// ----------------------------------------------------------------------------------------
int len = atoi(lenStr.encode().c_str());
bool coverAll = false;
svector<int_vector> options;
vector<int_vector, true, long long> options;
int_vector base;
TRACE("str", tout
@ -9227,16 +9239,16 @@ namespace smt {
);
// ----------------------------------------------------------------------------------------
ptr_vector<expr> orList;
ptr_vector<expr> andList;
for (long long i = l; i < h; i++) {
orList.push_back(m.mk_eq(val_indicator, mk_string(longlong_to_string(i).c_str()) ));
if (m_params.m_AggressiveValueTesting) {
literal l = mk_eq(val_indicator, mk_string(longlong_to_string(i).c_str()), false);
ctx.mark_as_relevant(l);
ctx.force_phase(l);
literal lit = mk_eq(val_indicator, mk_string(longlong_to_string(i).c_str()), false);
ctx.mark_as_relevant(lit);
ctx.force_phase(lit);
}
zstring aStr = gen_val_string(len, options[i - l]);
@ -10059,7 +10071,7 @@ namespace smt {
}
TRACE("str", tout << "last bounds are [" << lastBounds.lowerBound << " | " << lastBounds.midPoint << " | " << lastBounds.upperBound << "]!" << lastBounds.windowSize << std::endl;);
binary_search_info newBounds;
expr * newTester;
expr * newTester = 0;
if (lastTesterConstant == "more") {
// special case: if the midpoint, upper bound, and window size are all equal,
// we double the window size and adjust the bounds

2
src/smt/theory_str.h
View file

@ -350,6 +350,8 @@ protected:
unsigned long cacheHitCount;
unsigned long cacheMissCount;
unsigned m_fresh_id;
// cache mapping each string S to Length(S)
obj_map<expr, app*> length_ast_map;

8
src/tactic/portfolio/pb2bv_solver.cpp
View file

@ -20,6 +20,7 @@ Notes:
#include "solver_na2as.h"
#include "tactic.h"
#include "pb2bv_rewriter.h"
#include "th_rewriter.h"
#include "filter_model_converter.h"
#include "ast_pp.h"
#include "model_smt2_pp.h"
@ -29,6 +30,7 @@ class pb2bv_solver : public solver_na2as {
params_ref m_params;
mutable expr_ref_vector m_assertions;
mutable ref<solver> m_solver;
mutable th_rewriter m_th_rewriter;
mutable pb2bv_rewriter m_rewriter;
public:
@ -39,6 +41,7 @@ public:
m_params(p),
m_assertions(m),
m_solver(s),
m_th_rewriter(m, p),
m_rewriter(m, p)
{
}
@ -121,10 +124,11 @@ public:
private:
void flush_assertions() const {
proof_ref proof(m);
expr_ref fml(m);
expr_ref fml1(m), fml(m);
expr_ref_vector fmls(m);
for (unsigned i = 0; i < m_assertions.size(); ++i) {
m_rewriter(m_assertions[i].get(), fml, proof);
m_th_rewriter(m_assertions[i].get(), fml1, proof);
m_rewriter(fml1, fml, proof);
m_solver->assert_expr(fml);
}
m_rewriter.flush_side_constraints(fmls);

14
src/tactic/smtlogics/qflra_tactic.cpp
View file

@ -22,7 +22,6 @@ Notes:
#include"solve_eqs_tactic.h"
#include"elim_uncnstr_tactic.h"
#include"smt_tactic.h"
// include"mip_tactic.h"
#include"recover_01_tactic.h"
#include"ctx_simplify_tactic.h"
#include"probe_arith.h"
@ -72,5 +71,18 @@ tactic * mk_qflra_tactic(ast_manager & m, params_ref const & p) {
// using_params(mk_smt_tactic(), pivot_p)),
// p);
#if 0
params_ref simplex_0, simplex_1, simplex_2;
simplex_0.set_uint("lp.simplex_strategy", 0);
simplex_1.set_uint("lp.simplex_strategy", 1);
simplex_2.set_uint("lp.simplex_strategy", 2);
return par(using_params(mk_smt_tactic(), simplex_0),
using_params(mk_smt_tactic(), simplex_1),
using_params(mk_smt_tactic(), simplex_2));
#else
return using_params(using_params(mk_smt_tactic(), pivot_p), p);
#endif
}

144
src/test/argument_parser.h Normal file
View file

@ -0,0 +1,144 @@
/*
Copyright (c) 2013 Microsoft Corporation. All rights reserved.
Released under Apache 2.0 license as described in the file LICENSE.
Author: Lev Nachmanson
*/
#include <unordered_map>
#include <vector>
#include <string>
#include <set>
#include <iostream>
namespace lean {
class argument_parser {
std::unordered_map<std::string, std::string> m_options;
std::unordered_map<std::string, std::string> m_options_with_after_string;
std::set<std::string> m_used_options;
std::unordered_map<std::string, std::string> m_used_options_with_after_string;
std::vector<std::string> m_free_args;
std::vector<std::string> m_args;
public:
std::string m_error_message;
argument_parser(unsigned argn, char * const* args) {
for (unsigned i = 0; i < argn; i++) {
m_args.push_back(std::string(args[i]));
}
}
void add_option(std::string s) {
add_option_with_help_string(s, "");
}
void add_option_with_help_string(std::string s, std::string help_string) {
m_options[s]=help_string;
}
void add_option_with_after_string(std::string s) {
add_option_with_after_string_with_help(s, "");
}
void add_option_with_after_string_with_help(std::string s, std::string help_string) {
m_options_with_after_string[s]=help_string;
}
bool parse() {
bool status_is_ok = true;
for (unsigned i = 0; i < m_args.size(); i++) {
std::string ar = m_args[i];
if (m_options.find(ar) != m_options.end() )
m_used_options.insert(ar);
else if (m_options_with_after_string.find(ar) != m_options_with_after_string.end()) {
if (i == m_args.size() - 1) {
m_error_message = "Argument is missing after "+ar;
return false;
}
i++;
m_used_options_with_after_string[ar] = m_args[i];
} else {
if (starts_with(ar, "-") || starts_with(ar, "//"))
status_is_ok = false;
m_free_args.push_back(ar);
}
}
return status_is_ok;
}
bool contains(std::unordered_map<std::string, std::string> & m, std::string s) {
return m.find(s) != m.end();
}
bool contains(std::set<std::string> & m, std::string s) {
return m.find(s) != m.end();
}
bool option_is_used(std::string option) {
return contains(m_used_options, option) || contains(m_used_options_with_after_string, option);
}
std::string get_option_value(std::string option) {
auto t = m_used_options_with_after_string.find(option);
if (t != m_used_options_with_after_string.end()){
return t->second;
}
return std::string();
}
bool starts_with(std::string s, char const * prefix) {
return starts_with(s, std::string(prefix));
}
bool starts_with(std::string s, std::string prefix) {
return s.substr(0, prefix.size()) == prefix;
}
std::string usage_string() {
std::string ret = "";
std::vector<std::string> unknown_options;
for (auto t : m_free_args) {
if (starts_with(t, "-") || starts_with(t, "\\")) {
unknown_options.push_back(t);
}
}
if (unknown_options.size()) {
ret = "Unknown options:";
}
for (auto unknownOption : unknown_options) {
ret += unknownOption;
ret += ",";
}
ret += "\n";
ret += "Usage:\n";
for (auto allowed_option : m_options)
ret += allowed_option.first + " " + (allowed_option.second.size() == 0 ? std::string("") : std::string("/") + allowed_option.second) + std::string("\n");
for (auto s : m_options_with_after_string) {
ret += s.first + " " + (s.second.size() == 0? " \"option value\"":("\""+ s.second+"\"")) + "\n";
}
return ret;
}
void print() {
if (m_used_options.size() == 0 && m_used_options_with_after_string.size() == 0 && m_free_args.size() == 0) {
std::cout << "no options are given" << std::endl;
return;
}
std::cout << "options are: " << std::endl;
for (std::string s : m_used_options) {
std::cout << s << std::endl;
}
for (auto & t : m_used_options_with_after_string) {
std::cout << t.first << " " << t.second << std::endl;
}
if (m_free_args.size() > 0) {
std::cout << "free arguments are: " << std::endl;
for (auto & t : m_free_args) {
std::cout << t << " " << std::endl;
}
}
}
};
}

3236
src/test/lp.cpp Normal file
View file

File diff suppressed because it is too large Load diff

1
src/test/main.cpp
View file

@ -240,6 +240,7 @@ int main(int argc, char ** argv) {
TST(pdr);
TST_ARGV(ddnf);
TST(model_evaluator);
TST_ARGV(lp);
TST(get_consequences);
TST(pb2bv);
TST_ARGV(sat_lookahead);

396
src/test/smt_reader.h Normal file
View file

@ -0,0 +1,396 @@
/*
Copyright (c) 2013 Microsoft Corporation. All rights reserved.
Released under Apache 2.0 license as described in the file LICENSE.
Author: Lev Nachmanson
*/
#pragma once
// reads an MPS file reperesenting a Mixed Integer Program
#include <string>
#include <vector>
#include <unordered_map>
#include "util/lp/lp_primal_simplex.h"
#include "util/lp/lp_dual_simplex.h"
#include "util/lp/lar_solver.h"
#include <iostream>
#include <fstream>
#include <functional>
#include <algorithm>
#include "util/lp/mps_reader.h"
#include "util/lp/ul_pair.h"
#include "util/lp/lar_constraints.h"
#include <sstream>
#include <cstdlib>
namespace lean {
template<typename T>
T from_string(const std::string& str) {
std::istringstream ss(str);
T ret;
ss >> ret;
return ret;
}
class smt_reader {
public:
struct lisp_elem {
std::string m_head;
std::vector<lisp_elem> m_elems;
void print() {
if (m_elems.size()) {
std::cout << '(';
std::cout << m_head << ' ';
for (auto & el : m_elems)
el.print();
std::cout << ')';
} else {
std::cout << " " << m_head;
}
}
unsigned size() const { return static_cast<unsigned>(m_elems.size()); }
bool is_simple() const { return size() == 0; }
};
struct formula_constraint {
lconstraint_kind m_kind;
std::vector<std::pair<mpq, std::string>> m_coeffs;
mpq m_right_side;
void add_pair(mpq c, std::string name) {
m_coeffs.push_back(make_pair(c, name));
}
formula_constraint() : m_right_side(numeric_traits<mpq>::zero()) {}
};
lisp_elem m_formula_lisp_elem;
std::unordered_map<std::string, unsigned> m_name_to_var_index;
std::vector<formula_constraint> m_constraints;
bool m_is_OK;
unsigned m_line_number;
std::string m_file_name;
std::ifstream m_file_stream;
std::string m_line;
smt_reader(std::string file_name):
m_is_OK(true),
m_line_number(0),
m_file_name(file_name),
m_file_stream(file_name) {
}
void set_error() {
std::cout << "setting error" << std::endl;
m_is_OK = false;
}
bool is_ok() {
return m_is_OK;
}
bool prefix(const char * pr) {
return m_line.find(pr) == 0;
}
int first_separator() {
unsigned blank_pos = static_cast<unsigned>(m_line.find(' '));
unsigned br_pos = static_cast<unsigned>(m_line.find('('));
unsigned reverse_br_pos = static_cast<unsigned>(m_line.find(')'));
return std::min(blank_pos, std::min(br_pos, reverse_br_pos));
}
void fill_lisp_elem(lisp_elem & lm) {
if (m_line[0] == '(')
fill_nested_elem(lm);
else
fill_simple_elem(lm);
}
void fill_simple_elem(lisp_elem & lm) {
int separator = first_separator();
lean_assert(-1 != separator && separator != 0);
lm.m_head = m_line.substr(0, separator);
m_line = m_line.substr(separator);
}
void fill_nested_elem(lisp_elem & lm) {
lean_assert(m_line[0] == '(');
m_line = m_line.substr(1);
int separator = first_separator();
lm.m_head = m_line.substr(0, separator);
m_line = m_line.substr(lm.m_head.size());
eat_blanks();
while (m_line.size()) {
if (m_line[0] == '(') {
lisp_elem el;
fill_nested_elem(el);
lm.m_elems.push_back(el);
} else {
if (m_line[0] == ')') {
m_line = m_line.substr(1);
break;
}
lisp_elem el;
fill_simple_elem(el);
lm.m_elems.push_back(el);
}
eat_blanks();
}
}
void eat_blanks() {
while (m_line.size()) {
if (m_line[0] == ' ')
m_line = m_line.substr(1);
else
break;
}
}
void fill_formula_elem() {
fill_lisp_elem(m_formula_lisp_elem);
}
void parse_line() {
if (m_line.find(":formula") == 0) {
int first_br = static_cast<int>(m_line.find('('));
if (first_br == -1) {
std::cout << "empty formula" << std::endl;
return;
}
m_line = m_line.substr(first_br);
fill_formula_elem();
}
}
void set_constraint_kind(formula_constraint & c, lisp_elem & el) {
if (el.m_head == "=") {
c.m_kind = EQ;
} else if (el.m_head == ">=") {
c.m_kind = GE;
} else if (el.m_head == "<=") {
c.m_kind = LE;
} else if (el.m_head == ">") {
c.m_kind = GT;
} else if (el.m_head == "<") {
c.m_kind = LT;
} else {
std::cout << "kind " << el.m_head << " is not supported " << std::endl;
set_error();
}
}
void adjust_rigth_side(formula_constraint & /* c*/, lisp_elem & /*el*/) {
// lean_assert(el.m_head == "0"); // do nothing for the time being
}
void set_constraint_coeffs(formula_constraint & c, lisp_elem & el) {
lean_assert(el.m_elems.size() == 2);
set_constraint_coeffs_on_coeff_element(c, el.m_elems[0]);
adjust_rigth_side(c, el.m_elems[1]);
}
bool is_integer(std::string & s) {
if (s.size() == 0) return false;
return atoi(s.c_str()) != 0 || isdigit(s.c_str()[0]);
}
void add_complex_sum_elem(formula_constraint & c, lisp_elem & el) {
if (el.m_head == "*") {
add_mult_elem(c, el.m_elems);
} else if (el.m_head == "~") {
lisp_elem & minel = el.m_elems[0];
lean_assert(minel.is_simple());
c.m_right_side += mpq(str_to_int(minel.m_head));
} else {
std::cout << "unexpected input " << el.m_head << std::endl;
set_error();
return;
}
}
std::string get_name(lisp_elem & name) {
lean_assert(name.is_simple());
lean_assert(!is_integer(name.m_head));
return name.m_head;
}
void add_mult_elem(formula_constraint & c, std::vector<lisp_elem> & els) {
lean_assert(els.size() == 2);
mpq coeff = get_coeff(els[0]);
std::string col_name = get_name(els[1]);
c.add_pair(coeff, col_name);
}
mpq get_coeff(lisp_elem & le) {
if (le.is_simple()) {
return mpq(str_to_int(le.m_head));
} else {
lean_assert(le.m_head == "~");
lean_assert(le.size() == 1);
lisp_elem & el = le.m_elems[0];
lean_assert(el.is_simple());
return -mpq(str_to_int(el.m_head));
}
}
int str_to_int(std::string & s) {
lean_assert(is_integer(s));
return atoi(s.c_str());
}
void add_sum_elem(formula_constraint & c, lisp_elem & el) {
if (el.size()) {
add_complex_sum_elem(c, el);
} else {
lean_assert(is_integer(el.m_head));
int v = atoi(el.m_head.c_str());
mpq vr(v);
c.m_right_side -= vr;
}
}
void add_sum(formula_constraint & c, std::vector<lisp_elem> & sum_els) {
for (auto & el : sum_els)
add_sum_elem(c, el);
}
void set_constraint_coeffs_on_coeff_element(formula_constraint & c, lisp_elem & el) {
if (el.m_head == "*") {
add_mult_elem(c, el.m_elems);
} else if (el.m_head == "+") {
add_sum(c, el.m_elems);
} else {
lean_assert(false); // unexpected input
}
}
void create_constraint(lisp_elem & el) {
formula_constraint c;
set_constraint_kind(c, el);
set_constraint_coeffs(c, el);
m_constraints.push_back(c);
}
void fill_constraints() {
if (m_formula_lisp_elem.m_head != "and") {
std::cout << "unexpected top element " << m_formula_lisp_elem.m_head << std::endl;
set_error();
return;
}
for (auto & el : m_formula_lisp_elem.m_elems)
create_constraint(el);
}
void read() {
if (!m_file_stream.is_open()){
std::cout << "cannot open file " << m_file_name << std::endl;
set_error();
return;
}
while (m_is_OK && getline(m_file_stream, m_line)) {
parse_line();
m_line_number++;
}
m_file_stream.close();
fill_constraints();
}
/*
void fill_lar_solver_on_row(row * row, lar_solver *solver) {
if (row->m_name != m_cost_row_name) {
lar_constraint c(get_lar_relation_from_row(row->m_type), row->m_right_side);
for (auto s : row->m_row_columns) {
var_index i = solver->add_var(s.first);
c.add_variable_to_constraint(i, s.second);
}
solver->add_constraint(&c);
} else {
// ignore the cost row
}
}
void fill_lar_solver_on_rows(lar_solver * solver) {
for (auto row_it : m_rows) {
fill_lar_solver_on_row(row_it.second, solver);
}
}
void create_low_constraint_for_var(column* col, bound * b, lar_solver *solver) {
lar_constraint c(GE, b->m_low);
var_index i = solver->add_var(col->m_name);
c.add_variable_to_constraint(i, numeric_traits<T>::one());
solver->add_constraint(&c);
}
void create_upper_constraint_for_var(column* col, bound * b, lar_solver *solver) {
lar_constraint c(LE, b->m_upper);
var_index i = solver->add_var(col->m_name);
c.add_variable_to_constraint(i, numeric_traits<T>::one());
solver->add_constraint(&c);
}
void create_equality_contraint_for_var(column* col, bound * b, lar_solver *solver) {
lar_constraint c(EQ, b->m_fixed_value);
var_index i = solver->add_var(col->m_name);
c.add_variable_to_constraint(i, numeric_traits<T>::one());
solver->add_constraint(&c);
}
void fill_lar_solver_on_columns(lar_solver * solver) {
for (auto s : m_columns) {
mps_reader::column * col = s.second;
solver->add_var(col->m_name);
auto b = col->m_bound;
if (b == nullptr) return;
if (b->m_free) continue;
if (b->m_low_is_set) {
create_low_constraint_for_var(col, b, solver);
}
if (b->m_upper_is_set) {
create_upper_constraint_for_var(col, b, solver);
}
if (b->m_value_is_fixed) {
create_equality_contraint_for_var(col, b, solver);
}
}
}
*/
unsigned register_name(std::string s) {
auto it = m_name_to_var_index.find(s);
if (it!= m_name_to_var_index.end())
return it->second;
unsigned ret = static_cast<unsigned>(m_name_to_var_index.size());
m_name_to_var_index[s] = ret;
return ret;
}
void add_constraint_to_solver(lar_solver * solver, formula_constraint & fc) {
vector<std::pair<mpq, var_index>> ls;
for (auto & it : fc.m_coeffs) {
ls.push_back(std::make_pair(it.first, solver->add_var(register_name(it.second))));
}
solver->add_constraint(ls, fc.m_kind, fc.m_right_side);
}
void fill_lar_solver(lar_solver * solver) {
for (formula_constraint & fc : m_constraints)
add_constraint_to_solver(solver, fc);
}
lar_solver * create_lar_solver() {
lar_solver * ls = new lar_solver();
fill_lar_solver(ls);
return ls;
}
};
}

73
src/test/test_file_reader.h Normal file
View file

@ -0,0 +1,73 @@
/*
Copyright (c) 2013 Microsoft Corporation. All rights reserved.
Released under Apache 2.0 license as described in the file LICENSE.
Author: Lev Nachmanson
*/
#pragma once
// reads a text file
#include <string>
#include <vector>
#include <unordered_map>
#include <iostream>
#include <fstream>
#include "util/lp/lp_utils.h"
#include "util/lp/lp_solver.h"
namespace lean {
template <typename T>
struct test_result {
lp_status m_status;
T m_cost;
std::unordered_map<std::string, T> column_values;
};
template <typename T>
class test_file_reader {
struct raw_blob {
std::vector<std::string> m_unparsed_strings;
std::vector<raw_blob> m_blobs;
};
struct test_file_blob {
std::string m_name;
std::string m_content;
std::unordered_map<std::string, std::string> m_table;
std::unordered_map<std::string, test_file_blob> m_blobs;
test_result<T> * get_test_result() {
test_result<T> * tr = new test_result<T>();
throw "not impl";
return tr;
}
};
std::ifstream m_file_stream;
public:
// constructor
test_file_reader(std::string file_name) : m_file_stream(file_name) {
if (!m_file_stream.is_open()) {
std::cout << "cannot open file " << "\'" << file_name << "\'" << std::endl;
}
}
raw_blob scan_to_row_blob() {
}
test_file_blob scan_row_blob_to_test_file_blob(raw_blob /* rblob */) {
}
test_result<T> * get_test_result() {
if (!m_file_stream.is_open()) {
return nullptr;
}
raw_blob rblob = scan_to_row_blob();
test_file_blob tblob = scan_row_blob_to_test_file_blob(rblob);
return tblob.get_test_result();
}
};
}

2
src/util/hash.cpp
View file

@ -83,8 +83,8 @@ unsigned string_hash(const char * str, unsigned length, unsigned init_value) {
Z3_fallthrough;
case 1 :
a+=str[0];
Z3_fallthrough;
/* case 0: nothing left to add */
break;
}
mix(a,b,c);
/*-------------------------------------------- report the result */

70
src/util/lp/binary_heap_priority_queue.h Normal file
View file

@ -0,0 +1,70 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#pragma once
#include "util/vector.h"
#include "util/debug.h"
#include "util/lp/lp_utils.h"
namespace lean {
// the elements with the smallest priority are dequeued first
template <typename T>
class binary_heap_priority_queue {
vector<T> m_priorities;
// indexing for A starts from 1
vector<unsigned> m_heap; // keeps the elements of the queue
vector<int> m_heap_inverse; // o = m_heap[m_heap_inverse[o]]
unsigned m_heap_size;
// is is the child place in heap
void swap_with_parent(unsigned i);
void put_at(unsigned i, unsigned h);
void decrease_priority(unsigned o, T newPriority);
public:
#ifdef LEAN_DEBUG
bool is_consistent() const;
#endif
public:
void remove(unsigned o);
unsigned size() const { return m_heap_size; }
binary_heap_priority_queue(): m_heap(1), m_heap_size(0) {} // the empty constructror
// n is the initial queue capacity.
// The capacity will be enlarged two times automatically if needed
binary_heap_priority_queue(unsigned n);
void clear() {
for (unsigned i = 0; i < m_heap_size; i++) {
unsigned o = m_heap[i+1];
m_heap_inverse[o] = -1;
}
m_heap_size = 0;
}
void resize(unsigned n);
void put_to_heap(unsigned i, unsigned o);
void enqueue_new(unsigned o, const T& priority);
// This method can work with an element that is already in the queue.
// In this case the priority will be changed and the queue adjusted.
void enqueue(unsigned o, const T & priority);
void change_priority_for_existing(unsigned o, const T & priority);
T get_priority(unsigned o) const { return m_priorities[o]; }
bool is_empty() const { return m_heap_size == 0; }
/// return the first element of the queue and removes it from the queue
unsigned dequeue_and_get_priority(T & priority);
void fix_heap_under(unsigned i);
void put_the_last_at_the_top_and_fix_the_heap();
/// return the first element of the queue and removes it from the queue
unsigned dequeue();
unsigned peek() const {
lean_assert(m_heap_size > 0);
return m_heap[1];
}
#ifdef LEAN_DEBUG
void print(std::ostream & out);
#endif
};
}

194
src/util/lp/binary_heap_priority_queue.hpp Normal file
View file

@ -0,0 +1,194 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#include "util/vector.h"
#include "util/lp/binary_heap_priority_queue.h"
namespace lean {
// is is the child place in heap
template <typename T> void binary_heap_priority_queue<T>::swap_with_parent(unsigned i) {
unsigned parent = m_heap[i >> 1];
put_at(i >> 1, m_heap[i]);
put_at(i, parent);
}
template <typename T> void binary_heap_priority_queue<T>::put_at(unsigned i, unsigned h) {
m_heap[i] = h;
m_heap_inverse[h] = i;
}
template <typename T> void binary_heap_priority_queue<T>::decrease_priority(unsigned o, T newPriority) {
m_priorities[o] = newPriority;
int i = m_heap_inverse[o];
while (i > 1) {
if (m_priorities[m_heap[i]] < m_priorities[m_heap[i >> 1]])
swap_with_parent(i);
else
break;
i >>= 1;
}
}
#ifdef LEAN_DEBUG
template <typename T> bool binary_heap_priority_queue<T>::is_consistent() const {
for (int i = 0; i < m_heap_inverse.size(); i++) {
int i_index = m_heap_inverse[i];
lean_assert(i_index <= static_cast<int>(m_heap_size));
lean_assert(i_index == -1 || m_heap[i_index] == i);
}
for (unsigned i = 1; i < m_heap_size; i++) {
unsigned ch = i << 1;
for (int k = 0; k < 2; k++) {
if (ch > m_heap_size) break;
if (!(m_priorities[m_heap[i]] <= m_priorities[m_heap[ch]])){
return false;
}
ch++;
}
}
return true;
}
#endif
template <typename T> void binary_heap_priority_queue<T>::remove(unsigned o) {
T priority_of_o = m_priorities[o];
int o_in_heap = m_heap_inverse[o];
if (o_in_heap == -1) {
return; // nothing to do
}
lean_assert(static_cast<unsigned>(o_in_heap) <= m_heap_size);
if (static_cast<unsigned>(o_in_heap) < m_heap_size) {
put_at(o_in_heap, m_heap[m_heap_size--]);
if (m_priorities[m_heap[o_in_heap]] > priority_of_o) {
fix_heap_under(o_in_heap);
} else { // we need to propogate the m_heap[o_in_heap] up
unsigned i = o_in_heap;
while (i > 1) {
unsigned ip = i >> 1;
if (m_priorities[m_heap[i]] < m_priorities[m_heap[ip]])
swap_with_parent(i);
else
break;
i = ip;
}
}
} else {
lean_assert(static_cast<unsigned>(o_in_heap) == m_heap_size);
m_heap_size--;
}
m_heap_inverse[o] = -1;
// lean_assert(is_consistent());
}
// n is the initial queue capacity.
// The capacity will be enlarged two times automatically if needed
template <typename T> binary_heap_priority_queue<T>::binary_heap_priority_queue(unsigned n) :
m_priorities(n),
m_heap(n + 1), // because the indexing for A starts from 1
m_heap_inverse(n, -1),
m_heap_size(0)
{ }
template <typename T> void binary_heap_priority_queue<T>::resize(unsigned n) {
m_priorities.resize(n);
m_heap.resize(n + 1);
m_heap_inverse.resize(n, -1);
}
template <typename T> void binary_heap_priority_queue<T>::put_to_heap(unsigned i, unsigned o) {
m_heap[i] = o;
m_heap_inverse[o] = i;
}
template <typename T> void binary_heap_priority_queue<T>::enqueue_new(unsigned o, const T& priority) {
m_heap_size++;
int i = m_heap_size;
lean_assert(o < m_priorities.size());
m_priorities[o] = priority;
put_at(i, o);
while (i > 1 && m_priorities[m_heap[i >> 1]] > priority) {
swap_with_parent(i);
i >>= 1;
}
}
// This method can work with an element that is already in the queue.
// In this case the priority will be changed and the queue adjusted.
template <typename T> void binary_heap_priority_queue<T>::enqueue(unsigned o, const T & priority) {
if (o >= m_priorities.size()) {
resize(o << 1); // make the size twice larger
}
if (m_heap_inverse[o] == -1)
enqueue_new(o, priority);
else
change_priority_for_existing(o, priority);
}
template <typename T> void binary_heap_priority_queue<T>::change_priority_for_existing(unsigned o, const T & priority) {
if (m_priorities[o] > priority) {
decrease_priority(o, priority);
} else {
m_priorities[o] = priority;
fix_heap_under(m_heap_inverse[o]);
}
}
/// return the first element of the queue and removes it from the queue
template <typename T> unsigned binary_heap_priority_queue<T>::dequeue_and_get_priority(T & priority) {
lean_assert(m_heap_size != 0);
int ret = m_heap[1];
priority = m_priorities[ret];
put_the_last_at_the_top_and_fix_the_heap();
return ret;
}
template <typename T> void binary_heap_priority_queue<T>::fix_heap_under(unsigned i) {
while (true) {
unsigned smallest = i;
unsigned l = i << 1;
if (l <= m_heap_size && m_priorities[m_heap[l]] < m_priorities[m_heap[i]])
smallest = l;
unsigned r = l + 1;
if (r <= m_heap_size && m_priorities[m_heap[r]] < m_priorities[m_heap[smallest]])
smallest = r;
if (smallest != i)
swap_with_parent(smallest);
else
break;
i = smallest;
}
}
template <typename T> void binary_heap_priority_queue<T>::put_the_last_at_the_top_and_fix_the_heap() {
if (m_heap_size > 1) {
put_at(1, m_heap[m_heap_size--]);
fix_heap_under(1);
} else {
m_heap_size--;
}
}
/// return the first element of the queue and removes it from the queue
template <typename T> unsigned binary_heap_priority_queue<T>::dequeue() {
lean_assert(m_heap_size > 0);
int ret = m_heap[1];
put_the_last_at_the_top_and_fix_the_heap();
m_heap_inverse[ret] = -1;
return ret;
}
#ifdef LEAN_DEBUG
template <typename T> void binary_heap_priority_queue<T>::print(std::ostream & out) {
vector<int> index;
vector<T> prs;
while (size()) {
T prior;
int j = dequeue_and_get_priority(prior);
index.push_back(j);
prs.push_back(prior);
out << "(" << j << ", " << prior << ")";
}
out << std::endl;
// restore the queue
for (int i = 0; i < index.size(); i++)
enqueue(index[i], prs[i]);
}
#endif
}

26
src/util/lp/binary_heap_priority_queue_instances.cpp Normal file
View file

@ -0,0 +1,26 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#include "util/lp/numeric_pair.h"
#include "util/lp/binary_heap_priority_queue.hpp"
namespace lean {
template binary_heap_priority_queue<int>::binary_heap_priority_queue(unsigned int);
template unsigned binary_heap_priority_queue<int>::dequeue();
template void binary_heap_priority_queue<int>::enqueue(unsigned int, int const&);
template void binary_heap_priority_queue<double>::enqueue(unsigned int, double const&);
template void binary_heap_priority_queue<mpq>::enqueue(unsigned int, mpq const&);
template void binary_heap_priority_queue<int>::remove(unsigned int);
template unsigned binary_heap_priority_queue<numeric_pair<mpq> >::dequeue();
template unsigned binary_heap_priority_queue<double>::dequeue();
template unsigned binary_heap_priority_queue<mpq>::dequeue();
template void binary_heap_priority_queue<numeric_pair<mpq> >::enqueue(unsigned int, numeric_pair<mpq> const&);
template void binary_heap_priority_queue<numeric_pair<mpq> >::resize(unsigned int);
template void lean::binary_heap_priority_queue<double>::resize(unsigned int);
template binary_heap_priority_queue<unsigned int>::binary_heap_priority_queue(unsigned int);
template void binary_heap_priority_queue<unsigned>::resize(unsigned int);
template unsigned binary_heap_priority_queue<unsigned int>::dequeue();
template void binary_heap_priority_queue<unsigned int>::enqueue(unsigned int, unsigned int const&);
template void binary_heap_priority_queue<unsigned int>::remove(unsigned int);
template void lean::binary_heap_priority_queue<mpq>::resize(unsigned int);
}

50
src/util/lp/binary_heap_upair_queue.h Normal file
View file

@ -0,0 +1,50 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#pragma once
#include <unordered_set>
#include <unordered_map>
#include <queue>
#include "util/vector.h"
#include <set>
#include <utility>
#include "util/lp/binary_heap_priority_queue.h"
typedef std::pair<unsigned, unsigned> upair;
namespace lean {
template <typename T>
class binary_heap_upair_queue {
binary_heap_priority_queue<T> m_q;
std::unordered_map<upair, unsigned> m_pairs_to_index;
svector<upair> m_pairs; // inverse to index
svector<unsigned> m_available_spots;
public:
binary_heap_upair_queue(unsigned size);
unsigned dequeue_available_spot();
bool is_empty() const { return m_q.is_empty(); }
unsigned size() const {return m_q.size(); }
bool contains(unsigned i, unsigned j) const { return m_pairs_to_index.find(std::make_pair(i, j)) != m_pairs_to_index.end();
}
void remove(unsigned i, unsigned j);
bool ij_index_is_new(unsigned ij_index) const;
void enqueue(unsigned i, unsigned j, const T & priority);
void dequeue(unsigned & i, unsigned &j);
T get_priority(unsigned i, unsigned j) const;
#ifdef LEAN_DEBUG
bool pair_to_index_is_a_bijection() const;
bool available_spots_are_correct() const;
bool is_correct() const {
return m_q.is_consistent() && pair_to_index_is_a_bijection() && available_spots_are_correct();
}
#endif
void resize(unsigned size) { m_q.resize(size); }
};
}

110
src/util/lp/binary_heap_upair_queue.hpp Normal file
View file

@ -0,0 +1,110 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#include <set>
#include "util/lp/lp_utils.h"
#include "util/lp/binary_heap_upair_queue.h"
namespace lean {
template <typename T> binary_heap_upair_queue<T>::binary_heap_upair_queue(unsigned size) : m_q(size), m_pairs(size) {
for (unsigned i = 0; i < size; i++)
m_available_spots.push_back(i);
}
template <typename T> unsigned
binary_heap_upair_queue<T>::dequeue_available_spot() {
lean_assert(m_available_spots.empty() == false);
unsigned ret = m_available_spots.back();
m_available_spots.pop_back();
return ret;
}
template <typename T> void binary_heap_upair_queue<T>::remove(unsigned i, unsigned j) {
upair p(i, j);
auto it = m_pairs_to_index.find(p);
if (it == m_pairs_to_index.end())
return; // nothing to do
m_q.remove(it->second);
m_available_spots.push_back(it->second);
m_pairs_to_index.erase(it);
}
template <typename T> bool binary_heap_upair_queue<T>::ij_index_is_new(unsigned ij_index) const {
for (auto it : m_pairs_to_index) {
if (it.second == ij_index)
return false;
}
return true;
}
template <typename T> void binary_heap_upair_queue<T>::enqueue(unsigned i, unsigned j, const T & priority) {
upair p(i, j);
auto it = m_pairs_to_index.find(p);
unsigned ij_index;
if (it == m_pairs_to_index.end()) {
// it is a new pair, let us find a spot for it
if (m_available_spots.empty()) {
// we ran out of empty spots
unsigned size_was = static_cast<unsigned>(m_pairs.size());
unsigned new_size = size_was << 1;
for (unsigned i = size_was; i < new_size; i++)
m_available_spots.push_back(i);
m_pairs.resize(new_size);
}
ij_index = dequeue_available_spot();
// lean_assert(ij_index<m_pairs.size() && ij_index_is_new(ij_index));
m_pairs[ij_index] = p;
m_pairs_to_index[p] = ij_index;
} else {
ij_index = it->second;
}
m_q.enqueue(ij_index, priority);
}
template <typename T> void binary_heap_upair_queue<T>::dequeue(unsigned & i, unsigned &j) {
lean_assert(!m_q.is_empty());
unsigned ij_index = m_q.dequeue();
upair & p = m_pairs[ij_index];
i = p.first;
j = p.second;
m_available_spots.push_back(ij_index);
m_pairs_to_index.erase(p);
}
template <typename T> T binary_heap_upair_queue<T>::get_priority(unsigned i, unsigned j) const {
auto it = m_pairs_to_index.find(std::make_pair(i, j));
if (it == m_pairs_to_index.end())
return T(0xFFFFFF); // big number
return m_q.get_priority(it->second);
}
#ifdef LEAN_DEBUG
template <typename T> bool binary_heap_upair_queue<T>::pair_to_index_is_a_bijection() const {
std::set<int> tmp;
for (auto p : m_pairs_to_index) {
unsigned j = p.second;
unsigned size = tmp.size();
tmp.insert(j);
if (tmp.size() == size)
return false;
}
return true;
}
template <typename T> bool binary_heap_upair_queue<T>::available_spots_are_correct() const {
std::set<int> tmp;
for (auto p : m_available_spots){
tmp.insert(p);
}
if (tmp.size() != m_available_spots.size())
return false;
for (auto it : m_pairs_to_index)
if (tmp.find(it.second) != tmp.end())
return false;
return true;
}
#endif
}

17
src/util/lp/binary_heap_upair_queue_instances.cpp Normal file
View file

@ -0,0 +1,17 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#include "util/lp/binary_heap_upair_queue.hpp"
namespace lean {
template binary_heap_upair_queue<int>::binary_heap_upair_queue(unsigned int);
template binary_heap_upair_queue<unsigned int>::binary_heap_upair_queue(unsigned int);
template unsigned binary_heap_upair_queue<int>::dequeue_available_spot();
template unsigned binary_heap_upair_queue<unsigned int>::dequeue_available_spot();
template void binary_heap_upair_queue<int>::enqueue(unsigned int, unsigned int, int const&);
template void binary_heap_upair_queue<int>::remove(unsigned int, unsigned int);
template void binary_heap_upair_queue<unsigned int>::remove(unsigned int, unsigned int);
template void binary_heap_upair_queue<int>::dequeue(unsigned int&, unsigned int&);
template void binary_heap_upair_queue<unsigned int>::enqueue(unsigned int, unsigned int, unsigned int const&);
template void binary_heap_upair_queue<unsigned int>::dequeue(unsigned int&, unsigned int&);
}

335
src/util/lp/bound_analyzer_on_row.h Normal file
View file

@ -0,0 +1,335 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#pragma once
#include "util/vector.h"
#include "util/lp/linear_combination_iterator.h"
#include "implied_bound.h"
#include "test_bound_analyzer.h"
#include <functional>
#include "util/lp/bound_propagator.h"
// We have an equality : sum by j of row[j]*x[j] = rs
// We try to pin a var by pushing the total by using the variable bounds
// In a loop we drive the partial sum down, denoting the variables of this process by _u.
// In the same loop trying to pin variables by pushing the partial sum up, denoting the variable related to it by _l
namespace lean {
class bound_analyzer_on_row {
linear_combination_iterator<mpq> & m_it;
bound_propagator & m_bp;
unsigned m_row_or_term_index;
int m_column_of_u; // index of an unlimited from above monoid
// -1 means that such a value is not found, -2 means that at least two of such monoids were found
int m_column_of_l; // index of an unlimited from below monoid
impq m_rs;
public :
// constructor
bound_analyzer_on_row(
linear_combination_iterator<mpq> &it,
const numeric_pair<mpq>& rs,
unsigned row_or_term_index,
bound_propagator & bp
)
:
m_it(it),
m_bp(bp),
m_row_or_term_index(row_or_term_index),
m_column_of_u(-1),
m_column_of_l(-1),
m_rs(rs)
{}
unsigned j;
void analyze() {
mpq a; unsigned j;
while (((m_column_of_l != -2) || (m_column_of_u != -2)) && m_it.next(a, j))
analyze_bound_on_var_on_coeff(j, a);
if (m_column_of_u >= 0)
limit_monoid_u_from_below();
else if (m_column_of_u == -1)
limit_all_monoids_from_below();
if (m_column_of_l >= 0)
limit_monoid_l_from_above();
else if (m_column_of_l == -1)
limit_all_monoids_from_above();
}
bool bound_is_available(unsigned j, bool low_bound) {
return (low_bound && low_bound_is_available(j)) ||
(!low_bound && upper_bound_is_available(j));
}
bool upper_bound_is_available(unsigned j) const {
switch (m_bp.get_column_type(j))
{
case column_type::fixed:
case column_type::boxed:
case column_type::upper_bound:
return true;
default:
return false;
}
}
bool low_bound_is_available(unsigned j) const {
switch (m_bp.get_column_type(j))
{
case column_type::fixed:
case column_type::boxed:
case column_type::low_bound:
return true;
default:
return false;
}
}
const impq & ub(unsigned j) const {
lean_assert(upper_bound_is_available(j));
return m_bp.get_upper_bound(j);
}
const impq & lb(unsigned j) const {
lean_assert(low_bound_is_available(j));
return m_bp.get_low_bound(j);
}
const mpq & monoid_max_no_mult(bool a_is_pos, unsigned j, bool & strict) const {
if (a_is_pos) {
strict = !is_zero(ub(j).y);
return ub(j).x;
}
strict = !is_zero(lb(j).y);
return lb(j).x;
}
mpq monoid_max(const mpq & a, unsigned j) const {
if (is_pos(a)) {
return a * ub(j).x;
}
return a * lb(j).x;
}
mpq monoid_max(const mpq & a, unsigned j, bool & strict) const {
if (is_pos(a)) {
strict = !is_zero(ub(j).y);
return a * ub(j).x;
}
strict = !is_zero(lb(j).y);
return a * lb(j).x;
}
const mpq & monoid_min_no_mult(bool a_is_pos, unsigned j, bool & strict) const {
if (!a_is_pos) {
strict = !is_zero(ub(j).y);
return ub(j).x;
}
strict = !is_zero(lb(j).y);
return lb(j).x;
}
mpq monoid_min(const mpq & a, unsigned j, bool& strict) const {
if (is_neg(a)) {
strict = !is_zero(ub(j).y);
return a * ub(j).x;
}
strict = !is_zero(lb(j).y);
return a * lb(j).x;
}
mpq monoid_min(const mpq & a, unsigned j) const {
if (is_neg(a)) {
return a * ub(j).x;
}
return a * lb(j).x;
}
void limit_all_monoids_from_above() {
int strict = 0;
mpq total;
lean_assert(is_zero(total));
m_it.reset();
mpq a; unsigned j;
while (m_it.next(a, j)) {
bool str;
total -= monoid_min(a, j, str);
if (str)
strict++;
}
m_it.reset();
while (m_it.next(a, j)) {
bool str;
bool a_is_pos = is_pos(a);
mpq bound = total / a + monoid_min_no_mult(a_is_pos, j, str);
if (a_is_pos) {
limit_j(j, bound, true, false, strict - static_cast<int>(str) > 0);
}
else {
limit_j(j, bound, false, true, strict - static_cast<int>(str) > 0);
}
}
}
void limit_all_monoids_from_below() {
int strict = 0;
mpq total;
lean_assert(is_zero(total));
m_it.reset();
mpq a; unsigned j;
while (m_it.next(a, j)) {
bool str;
total -= monoid_max(a, j, str);
if (str)
strict++;
}
m_it.reset();
while (m_it.next(a, j)) {
bool str;
bool a_is_pos = is_pos(a);
mpq bound = total / a + monoid_max_no_mult(a_is_pos, j, str);
bool astrict = strict - static_cast<int>(str) > 0;
if (a_is_pos) {
limit_j(j, bound, true, true, astrict);
}
else {
limit_j(j, bound, false, false, astrict);
}
}
}
void limit_monoid_u_from_below() {
// we are going to limit from below the monoid m_column_of_u,
// every other monoid is impossible to limit from below
mpq u_coeff, a;
unsigned j;
mpq bound = -m_rs.x;
m_it.reset();
bool strict = false;
while (m_it.next(a, j)) {
if (j == static_cast<unsigned>(m_column_of_u)) {
u_coeff = a;
continue;
}
bool str;
bound -= monoid_max(a, j, str);
if (str)
strict = true;
}
bound /= u_coeff;
if (numeric_traits<impq>::is_pos(u_coeff)) {
limit_j(m_column_of_u, bound, true, true, strict);
} else {
limit_j(m_column_of_u, bound, false, false, strict);
}
}
void limit_monoid_l_from_above() {
// we are going to limit from above the monoid m_column_of_l,
// every other monoid is impossible to limit from above
mpq l_coeff, a;
unsigned j;
mpq bound = -m_rs.x;
bool strict = false;
m_it.reset();
while (m_it.next(a, j)) {
if (j == static_cast<unsigned>(m_column_of_l)) {
l_coeff = a;
continue;
}
bool str;
bound -= monoid_min(a, j, str);
if (str)
strict = true;
}
bound /= l_coeff;
if (is_pos(l_coeff)) {
limit_j(m_column_of_l, bound, true, false, strict);
} else {
limit_j(m_column_of_l, bound, false, true, strict);
}
}
// // it is the coefficent before the bounded column
// void provide_evidence(bool coeff_is_pos) {
// /*
// auto & be = m_ibounds.back();
// bool low_bound = be.m_low_bound;
// if (!coeff_is_pos)
// low_bound = !low_bound;
// auto it = m_it.clone();
// mpq a; unsigned j;
// while (it->next(a, j)) {
// if (be.m_j == j) continue;
// lean_assert(bound_is_available(j, is_neg(a) ? low_bound : !low_bound));
// be.m_vector_of_bound_signatures.emplace_back(a, j, numeric_traits<impq>::
// is_neg(a)? low_bound: !low_bound);
// }
// delete it;
// */
// }
void limit_j(unsigned j, const mpq& u, bool coeff_before_j_is_pos, bool is_low_bound, bool strict){
m_bp.try_add_bound(u, j, is_low_bound, coeff_before_j_is_pos, m_row_or_term_index, strict);
}
void advance_u(unsigned j) {
if (m_column_of_u == -1)
m_column_of_u = j;
else
m_column_of_u = -2;
}
void advance_l(unsigned j) {
if (m_column_of_l == -1)
m_column_of_l = j;
else
m_column_of_l = -2;
}
void analyze_bound_on_var_on_coeff(int j, const mpq &a) {
switch (m_bp.get_column_type(j)) {
case column_type::low_bound:
if (numeric_traits<mpq>::is_pos(a))
advance_u(j);
else
advance_l(j);
break;
case column_type::upper_bound:
if(numeric_traits<mpq>::is_neg(a))
advance_u(j);
else
advance_l(j);
break;
case column_type::free_column:
advance_u(j);
advance_l(j);
break;
default:
break;
}
}
static void analyze_row(linear_combination_iterator<mpq> &it,
const numeric_pair<mpq>& rs,
unsigned row_or_term_index,
bound_propagator & bp
) {
bound_analyzer_on_row a(it, rs, row_or_term_index, bp);
a.analyze();
}
};
}

47
src/util/lp/bound_propagator.cpp Normal file
View file

@ -0,0 +1,47 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#include "util/lp/lar_solver.h"
namespace lean {
bound_propagator::bound_propagator(lar_solver & ls):
m_lar_solver(ls) {}
column_type bound_propagator::get_column_type(unsigned j) const {
return m_lar_solver.m_mpq_lar_core_solver.m_column_types()[j];
}
const impq & bound_propagator::get_low_bound(unsigned j) const {
return m_lar_solver.m_mpq_lar_core_solver.m_r_low_bounds()[j];
}
const impq & bound_propagator::get_upper_bound(unsigned j) const {
return m_lar_solver.m_mpq_lar_core_solver.m_r_upper_bounds()[j];
}
void bound_propagator::try_add_bound(const mpq & v, unsigned j, bool is_low, bool coeff_before_j_is_pos, unsigned row_or_term_index, bool strict) {
j = m_lar_solver.adjust_column_index_to_term_index(j);
lconstraint_kind kind = is_low? GE : LE;
if (strict)
kind = static_cast<lconstraint_kind>(kind / 2);
if (!bound_is_interesting(j, kind, v))
return;
unsigned k; // index to ibounds
if (is_low) {
if (try_get_val(m_improved_low_bounds, j, k)) {
auto & found_bound = m_ibounds[k];
if (v > found_bound.m_bound || (v == found_bound.m_bound && found_bound.m_strict == false && strict))
found_bound = implied_bound(v, j, is_low, coeff_before_j_is_pos, row_or_term_index, strict);
} else {
m_improved_low_bounds[j] = m_ibounds.size();
m_ibounds.push_back(implied_bound(v, j, is_low, coeff_before_j_is_pos, row_or_term_index, strict));
}
} else { // the upper bound case
if (try_get_val(m_improved_upper_bounds, j, k)) {
auto & found_bound = m_ibounds[k];
if (v < found_bound.m_bound || (v == found_bound.m_bound && found_bound.m_strict == false && strict))
found_bound = implied_bound(v, j, is_low, coeff_before_j_is_pos, row_or_term_index, strict);
} else {
m_improved_upper_bounds[j] = m_ibounds.size();
m_ibounds.push_back(implied_bound(v, j, is_low, coeff_before_j_is_pos, row_or_term_index, strict));
}
}
}
}

27
src/util/lp/bound_propagator.h Normal file
View file

@ -0,0 +1,27 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#pragma once
#include "util/lp/lp_settings.h"
namespace lean {
class lar_solver;
class bound_propagator {
std::unordered_map<unsigned, unsigned> m_improved_low_bounds; // these maps map a column index to the corresponding index in ibounds
std::unordered_map<unsigned, unsigned> m_improved_upper_bounds;
lar_solver & m_lar_solver;
public:
vector<implied_bound> m_ibounds;
public:
bound_propagator(lar_solver & ls);
column_type get_column_type(unsigned) const;
const impq & get_low_bound(unsigned) const;
const impq & get_upper_bound(unsigned) const;
void try_add_bound(const mpq & v, unsigned j, bool is_low, bool coeff_before_j_is_pos, unsigned row_or_term_index, bool strict);
virtual bool bound_is_interesting(unsigned vi,
lean::lconstraint_kind kind,
const rational & bval) {return true;}
unsigned number_of_found_bounds() const { return m_ibounds.size(); }
virtual void consume(mpq const& v, unsigned j) { std::cout << "doh\n"; }
};
}

20
src/util/lp/breakpoint.h Normal file
View file

@ -0,0 +1,20 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#pragma once
namespace lean {
enum breakpoint_type {
low_break, upper_break, fixed_break
};
template <typename X>
struct breakpoint {
unsigned m_j; // the basic column
breakpoint_type m_type;
X m_delta;
breakpoint(){}
breakpoint(unsigned j, X delta, breakpoint_type type):m_j(j), m_type(type), m_delta(delta) {}
};
}

250
src/util/lp/column_info.h Normal file
View file

@ -0,0 +1,250 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#pragma once
#include "util/vector.h"
#include <unordered_map>
#include <string>
#include <algorithm>
#include "util/lp/lp_settings.h"
namespace lean {
inline bool is_valid(unsigned j) { return static_cast<int>(j) >= 0;}
template <typename T>
class column_info {
std::string m_name;
bool m_low_bound_is_set;
bool m_low_bound_is_strict;
bool m_upper_bound_is_set;
bool m_upper_bound_is_strict;
T m_low_bound;
T m_upper_bound;
T m_fixed_value;
bool m_is_fixed;
T m_cost;
unsigned m_column_index;
public:
bool operator==(const column_info & c) const {
return m_name == c.m_name &&
m_low_bound_is_set == c.m_low_bound_is_set &&
m_low_bound_is_strict == c.m_low_bound_is_strict &&
m_upper_bound_is_set == c.m_upper_bound_is_set&&
m_upper_bound_is_strict == c.m_upper_bound_is_strict&&
(!m_low_bound_is_set || m_low_bound == c.m_low_bound) &&
(!m_upper_bound_is_set || m_upper_bound == c.m_upper_bound) &&
m_cost == c.m_cost &&
m_is_fixed == c.m_is_fixed &&
(!m_is_fixed || m_fixed_value == c.m_fixed_value) &&
m_column_index == c.m_column_index;
}
bool operator!=(const column_info & c) const { return !((*this) == c); }
void set_column_index(unsigned j) {
m_column_index = j;
}
// the default constructor
column_info():
m_low_bound_is_set(false),
m_low_bound_is_strict(false),
m_upper_bound_is_set (false),
m_upper_bound_is_strict (false),
m_is_fixed(false),
m_cost(numeric_traits<T>::zero()),
m_column_index(static_cast<unsigned>(-1))
{}
column_info(unsigned column_index) :
m_low_bound_is_set(false),
m_low_bound_is_strict(false),
m_upper_bound_is_set (false),
m_upper_bound_is_strict (false),
m_is_fixed(false),
m_cost(numeric_traits<T>::zero()),
m_column_index(column_index) {
}
column_info(const column_info & ci) {
m_name = ci.m_name;
m_low_bound_is_set = ci.m_low_bound_is_set;
m_low_bound_is_strict = ci.m_low_bound_is_strict;
m_upper_bound_is_set = ci.m_upper_bound_is_set;
m_upper_bound_is_strict = ci.m_upper_bound_is_strict;
m_low_bound = ci.m_low_bound;
m_upper_bound = ci.m_upper_bound;
m_cost = ci.m_cost;
m_fixed_value = ci.m_fixed_value;
m_is_fixed = ci.m_is_fixed;
m_column_index = ci.m_column_index;
}
unsigned get_column_index() const {
return m_column_index;
}
column_type get_column_type() const {
return m_is_fixed? column_type::fixed : (m_low_bound_is_set? (m_upper_bound_is_set? column_type::boxed : column_type::low_bound) : (m_upper_bound_is_set? column_type::upper_bound: column_type::free_column));
}
column_type get_column_type_no_flipping() const {
if (m_is_fixed) {
return column_type::fixed;
}
if (m_low_bound_is_set) {
return m_upper_bound_is_set? column_type::boxed: column_type::low_bound;
}
// we are flipping the bounds!
return m_upper_bound_is_set? column_type::upper_bound
: column_type::free_column;
}
T get_low_bound() const {
lean_assert(m_low_bound_is_set);
return m_low_bound;
}
T get_upper_bound() const {
lean_assert(m_upper_bound_is_set);
return m_upper_bound;
}
bool low_bound_is_set() const {
return m_low_bound_is_set;
}
bool upper_bound_is_set() const {
return m_upper_bound_is_set;
}
T get_shift() {
if (is_fixed()) {
return m_fixed_value;
}
if (is_flipped()){
return m_upper_bound;
}
return m_low_bound_is_set? m_low_bound : numeric_traits<T>::zero();
}
bool is_flipped() {
return m_upper_bound_is_set && !m_low_bound_is_set;
}
bool adjusted_low_bound_is_set() {
return !is_flipped()? low_bound_is_set(): upper_bound_is_set();
}
bool adjusted_upper_bound_is_set() {
return !is_flipped()? upper_bound_is_set(): low_bound_is_set();
}
T get_adjusted_upper_bound() {
return get_upper_bound() - get_low_bound();
}
bool is_fixed() const {
return m_is_fixed;
}
bool is_free() {
return !m_low_bound_is_set && !m_upper_bound_is_set;
}
void set_fixed_value(T v) {
m_is_fixed = true;
m_fixed_value = v;
}
T get_fixed_value() const {
lean_assert(m_is_fixed);
return m_fixed_value;
}
T get_cost() const {
return m_cost;
}
void set_cost(T const & cost) {
m_cost = cost;
}
void set_name(std::string const & s) {
m_name = s;
}
std::string get_name() const {
return m_name;
}
void set_low_bound(T const & l) {
m_low_bound = l;
m_low_bound_is_set = true;
}
void set_upper_bound(T const & l) {
m_upper_bound = l;
m_upper_bound_is_set = true;
}
void unset_low_bound() {
m_low_bound_is_set = false;
}
void unset_upper_bound() {
m_upper_bound_is_set = false;
}
void unset_fixed() {
m_is_fixed = false;
}
bool low_bound_holds(T v) {
return !low_bound_is_set() || v >= m_low_bound -T(0.0000001);
}
bool upper_bound_holds(T v) {
return !upper_bound_is_set() || v <= m_upper_bound + T(0.000001);
}
bool bounds_hold(T v) {
return low_bound_holds(v) && upper_bound_holds(v);
}
bool adjusted_bounds_hold(T v) {
return adjusted_low_bound_holds(v) && adjusted_upper_bound_holds(v);
}
bool adjusted_low_bound_holds(T v) {
return !adjusted_low_bound_is_set() || v >= -T(0.0000001);
}
bool adjusted_upper_bound_holds(T v) {
return !adjusted_upper_bound_is_set() || v <= get_adjusted_upper_bound() + T(0.000001);
}
bool is_infeasible() {
if ((!upper_bound_is_set()) || (!low_bound_is_set()))
return false;
// ok, both bounds are set
bool at_least_one_is_strict = upper_bound_is_strict() || low_bound_is_strict();
if (!at_least_one_is_strict)
return get_upper_bound() < get_low_bound();
// at least on bound is strict
return get_upper_bound() <= get_low_bound(); // the equality is impossible
}
bool low_bound_is_strict() const {
return m_low_bound_is_strict;
}
void set_low_bound_strict(bool val) {
m_low_bound_is_strict = val;
}
bool upper_bound_is_strict() const {
return m_upper_bound_is_strict;
}
void set_upper_bound_strict(bool val) {
m_upper_bound_is_strict = val;
}
};
}

82
src/util/lp/column_namer.h Normal file
View file

@ -0,0 +1,82 @@
#pragma once
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#include <string>
#include "util/lp/linear_combination_iterator.h"
namespace lean {
class column_namer {
public:
virtual std::string get_column_name(unsigned j) const = 0;
template <typename T>
void print_linear_iterator(linear_combination_iterator<T>* it, std::ostream & out) const {
vector<std::pair<T, unsigned>> coeff;
T a;
unsigned i;
while (it->next(a, i)) {
coeff.emplace_back(a, i);
}
print_linear_combination_of_column_indices(coeff, out);
}
template <typename T>
void print_linear_iterator_indices_only(linear_combination_iterator<T>* it, std::ostream & out) const {
vector<std::pair<T, unsigned>> coeff;
T a;
unsigned i;
while (it->next(a, i)) {
coeff.emplace_back(a, i);
}
print_linear_combination_of_column_indices_only(coeff, out);
}
template <typename T>
void print_linear_combination_of_column_indices_only(const vector<std::pair<T, unsigned>> & coeffs, std::ostream & out) const {
bool first = true;
for (const auto & it : coeffs) {
auto val = it.first;
if (first) {
first = false;
} else {
if (numeric_traits<T>::is_pos(val)) {
out << " + ";
} else {
out << " - ";
val = -val;
}
}
if (val == -numeric_traits<T>::one())
out << " - ";
else if (val != numeric_traits<T>::one())
out << T_to_string(val);
out << "_" << it.second;
}
}
template <typename T>
void print_linear_combination_of_column_indices(const vector<std::pair<T, unsigned>> & coeffs, std::ostream & out) const {
bool first = true;
for (const auto & it : coeffs) {
auto val = it.first;
if (first) {
first = false;
} else {
if (numeric_traits<T>::is_pos(val)) {
out << " + ";
} else {
out << " - ";
val = -val;
}
}
if (val == -numeric_traits<T>::one())
out << " - ";
else if (val != numeric_traits<T>::one())
out << val;
out << get_column_name(it.second);
}
}
};
}

43
src/util/lp/conversion_helper.h Normal file
View file

@ -0,0 +1,43 @@
/*
Copyright (c) 2013 Microsoft Corporation. All rights reserved.
Author: Lev Nachmanson
*/
#pragma once
namespace lean {
template <typename V>
struct conversion_helper {
static V get_low_bound(const column_info<mpq> & ci) {
return V(ci.get_low_bound(), ci.low_bound_is_strict()? 1 : 0);
}
static V get_upper_bound(const column_info<mpq> & ci) {
return V(ci.get_upper_bound(), ci.upper_bound_is_strict()? -1 : 0);
}
};
template<>
struct conversion_helper <double> {
static double get_upper_bound(const column_info<mpq> & ci) {
if (!ci.upper_bound_is_strict())
return ci.get_upper_bound().get_double();
double eps = 0.00001;
if (!ci.low_bound_is_set())
return ci.get_upper_bound().get_double() - eps;
eps = std::min((ci.get_upper_bound() - ci.get_low_bound()).get_double() / 1000, eps);
return ci.get_upper_bound().get_double() - eps;
}
static double get_low_bound(const column_info<mpq> & ci) {
if (!ci.low_bound_is_strict())
return ci.get_low_bound().get_double();
double eps = 0.00001;
if (!ci.upper_bound_is_set())
return ci.get_low_bound().get_double() + eps;
eps = std::min((ci.get_upper_bound() - ci.get_low_bound()).get_double() / 1000, eps);
return ci.get_low_bound().get_double() + eps;
}
};
}

117
src/util/lp/core_solver_pretty_printer.h Normal file
View file

@ -0,0 +1,117 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#pragma once
#include <limits>
#include <string>
#include <algorithm>
#include "util/vector.h"
#include <ostream>
#include "util/lp/lp_settings.h"
#include "util/lp/indexed_vector.h"
namespace lean {
template <typename T, typename X> class lp_core_solver_base; // forward definition
template <typename T, typename X>
class core_solver_pretty_printer {
std::ostream & m_out;
typedef std::string string;
lp_core_solver_base<T, X> & m_core_solver;
vector<unsigned> m_column_widths;
vector<vector<string>> m_A;
vector<vector<string>> m_signs;
vector<string> m_costs;
vector<string> m_cost_signs;
vector<string> m_lows; // low bounds
vector<string> m_upps; // upper bounds
vector<string> m_lows_signs;
vector<string> m_upps_signs;
unsigned m_rs_width;
vector<X> m_rs;
unsigned m_title_width;
std::string m_cost_title;
std::string m_basis_heading_title;
std::string m_x_title;
std::string m_low_bounds_title;
std::string m_upp_bounds_title;
std::string m_exact_norm_title;
std::string m_approx_norm_title;
unsigned ncols() { return m_core_solver.m_A.column_count(); }
unsigned nrows() { return m_core_solver.m_A.row_count(); }
unsigned m_artificial_start;
indexed_vector<T> m_w_buff;
indexed_vector<T> m_ed_buff;
vector<T> m_exact_column_norms;
public:
core_solver_pretty_printer(lp_core_solver_base<T, X > & core_solver, std::ostream & out);
void init_costs();
~core_solver_pretty_printer();
void init_rs_width();
T current_column_norm();
void init_m_A_and_signs();
void init_column_widths();
void adjust_width_with_low_bound(unsigned column, unsigned & w);
void adjust_width_with_upper_bound(unsigned column, unsigned & w);
void adjust_width_with_bounds(unsigned column, unsigned & w);
void adjust_width_with_basis_heading(unsigned column, unsigned & w) {
w = std::max(w, (unsigned)T_to_string(m_core_solver.m_basis_heading[column]).size());
}
unsigned get_column_width(unsigned column);
unsigned regular_cell_width(unsigned row, unsigned column, std::string name) {
return regular_cell_string(row, column, name).size();
}
std::string regular_cell_string(unsigned row, unsigned column, std::string name);
void set_coeff(vector<string>& row, vector<string> & row_signs, unsigned col, const T & t, string name);
void print_x();
std::string get_low_bound_string(unsigned j);
std::string get_upp_bound_string(unsigned j);
void print_lows();
void print_upps();
string get_exact_column_norm_string(unsigned col) {
return T_to_string(m_exact_column_norms[col]);
}
void print_exact_norms();
void print_approx_norms();
void print();
void print_basis_heading();
void print_bottom_line() {
m_out << "----------------------" << std::endl;
}
void print_cost();
void print_given_rows(vector<string> & row, vector<string> & signs, X rst);
void print_row(unsigned i);
};
}

383
src/util/lp/core_solver_pretty_printer.hpp Normal file
View file

@ -0,0 +1,383 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#include <limits>
#include <string>
#include <algorithm>
#include "util/lp/lp_utils.h"
#include "util/lp/lp_core_solver_base.h"
#include "util/lp/core_solver_pretty_printer.h"
#include "util/lp/numeric_pair.h"
namespace lean {
template <typename T, typename X>
core_solver_pretty_printer<T, X>::core_solver_pretty_printer(lp_core_solver_base<T, X > & core_solver, std::ostream & out):
m_out(out),
m_core_solver(core_solver),
m_A(core_solver.m_A.row_count(), vector<string>(core_solver.m_A.column_count(), "")),
m_signs(core_solver.m_A.row_count(), vector<string>(core_solver.m_A.column_count(), " ")),
m_costs(ncols(), ""),
m_cost_signs(ncols(), " "),
m_rs(ncols(), zero_of_type<X>()),
m_w_buff(core_solver.m_w),
m_ed_buff(core_solver.m_ed) {
m_low_bounds_title = "low";
m_upp_bounds_title = "upp";
m_exact_norm_title = "exact cn";
m_approx_norm_title = "approx cn";
m_artificial_start = std::numeric_limits<unsigned>::max();
m_column_widths.resize(core_solver.m_A.column_count(), 0),
init_m_A_and_signs();
init_costs();
init_column_widths();
init_rs_width();
m_cost_title = "costs";
m_basis_heading_title = "heading";
m_x_title = "x*";
m_title_width = static_cast<unsigned>(std::max(std::max(m_cost_title.size(), std::max(m_basis_heading_title.size(), m_x_title.size())), m_approx_norm_title.size()));
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::init_costs() {
if (!m_core_solver.use_tableau()) {
vector<T> local_y(m_core_solver.m_m());
m_core_solver.solve_yB(local_y);
for (unsigned i = 0; i < ncols(); i++) {
if (m_core_solver.m_basis_heading[i] < 0) {
T t = m_core_solver.m_costs[i] - m_core_solver.m_A.dot_product_with_column(local_y, i);
set_coeff(m_costs, m_cost_signs, i, t, m_core_solver.column_name(i));
}
}
} else {
for (unsigned i = 0; i < ncols(); i++) {
if (m_core_solver.m_basis_heading[i] < 0) {
set_coeff(m_costs, m_cost_signs, i, m_core_solver.m_d[i], m_core_solver.column_name(i));
}
}
}
}
template <typename T, typename X> core_solver_pretty_printer<T, X>::~core_solver_pretty_printer() {
m_core_solver.m_w = m_w_buff;
m_core_solver.m_ed = m_ed_buff;
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::init_rs_width() {
m_rs_width = static_cast<unsigned>(T_to_string(m_core_solver.get_cost()).size());
for (unsigned i = 0; i < nrows(); i++) {
unsigned wt = static_cast<unsigned>(T_to_string(m_rs[i]).size());
if (wt > m_rs_width) {
m_rs_width = wt;
}
}
}
template <typename T, typename X> T core_solver_pretty_printer<T, X>::current_column_norm() {
T ret = zero_of_type<T>();
for (auto i : m_core_solver.m_ed.m_index)
ret += m_core_solver.m_ed[i] * m_core_solver.m_ed[i];
return ret;
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::init_m_A_and_signs() {
if (numeric_traits<T>::precise() && m_core_solver.m_settings.use_tableau()) {
for (unsigned column = 0; column < ncols(); column++) {
vector<T> t(nrows(), zero_of_type<T>());
for (const auto & c : m_core_solver.m_A.m_columns[column]){
t[c.m_i] = m_core_solver.m_A.get_val(c);
}
string name = m_core_solver.column_name(column);
for (unsigned row = 0; row < nrows(); row ++) {
set_coeff(
m_A[row],
m_signs[row],
column,
t[row],
name);
m_rs[row] += t[row] * m_core_solver.m_x[column];
}
}
} else {
for (unsigned column = 0; column < ncols(); column++) {
m_core_solver.solve_Bd(column); // puts the result into m_core_solver.m_ed
string name = m_core_solver.column_name(column);
for (unsigned row = 0; row < nrows(); row ++) {
set_coeff(
m_A[row],
m_signs[row],
column,
m_core_solver.m_ed[row],
name);
m_rs[row] += m_core_solver.m_ed[row] * m_core_solver.m_x[column];
}
if (!m_core_solver.use_tableau())
m_exact_column_norms.push_back(current_column_norm() + T(1)); // a conversion missing 1 -> T
}
}
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::init_column_widths() {
for (unsigned i = 0; i < ncols(); i++) {
m_column_widths[i] = get_column_width(i);
}
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::adjust_width_with_low_bound(unsigned column, unsigned & w) {
if (!m_core_solver.low_bounds_are_set()) return;
w = std::max(w, (unsigned)T_to_string(m_core_solver.low_bound_value(column)).size());
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::adjust_width_with_upper_bound(unsigned column, unsigned & w) {
w = std::max(w, (unsigned)T_to_string(m_core_solver.upper_bound_value(column)).size());
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::adjust_width_with_bounds(unsigned column, unsigned & w) {
switch (m_core_solver.get_column_type(column)) {
case column_type::fixed:
case column_type::boxed:
adjust_width_with_low_bound(column, w);
adjust_width_with_upper_bound(column, w);
break;
case column_type::low_bound:
adjust_width_with_low_bound(column, w);
break;
case column_type::upper_bound:
adjust_width_with_upper_bound(column, w);
break;
case column_type::free_column:
break;
default:
lean_assert(false);
break;
}
}
template <typename T, typename X> unsigned core_solver_pretty_printer<T, X>:: get_column_width(unsigned column) {
unsigned w = static_cast<unsigned>(std::max((size_t)m_costs[column].size(), T_to_string(m_core_solver.m_x[column]).size()));
adjust_width_with_bounds(column, w);
adjust_width_with_basis_heading(column, w);
for (unsigned i = 0; i < nrows(); i++) {
unsigned cellw = static_cast<unsigned>(m_A[i][column].size());
if (cellw > w) {
w = cellw;
}
}
if (!m_core_solver.use_tableau()) {
w = std::max(w, (unsigned)T_to_string(m_exact_column_norms[column]).size());
if (m_core_solver.m_column_norms.size() > 0)
w = std::max(w, (unsigned)T_to_string(m_core_solver.m_column_norms[column]).size());
}
return w;
}
template <typename T, typename X> std::string core_solver_pretty_printer<T, X>::regular_cell_string(unsigned row, unsigned /* column */, std::string name) {
T t = fabs(m_core_solver.m_ed[row]);
if ( t == 1) return name;
return T_to_string(t) + name;
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::set_coeff(vector<string>& row, vector<string> & row_signs, unsigned col, const T & t, string name) {
if (numeric_traits<T>::is_zero(t)) {
return;
}
if (col > 0) {
if (t > 0) {
row_signs[col] = "+";
row[col] = t != 1? T_to_string(t) + name : name;
} else {
row_signs[col] = "-";
row[col] = t != -1? T_to_string(-t) + name: name;
}
} else { // col == 0
if (t == -1) {
row[col] = "-" + name;
} else if (t == 1) {
row[col] = name;
} else {
row[col] = T_to_string(t) + name;
}
}
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::print_x() {
if (ncols() == 0) {
return;
}
int blanks = m_title_width + 1 - static_cast<int>(m_x_title.size());
m_out << m_x_title;
print_blanks(blanks, m_out);
auto bh = m_core_solver.m_x;
for (unsigned i = 0; i < ncols(); i++) {
string s = T_to_string(bh[i]);
int blanks = m_column_widths[i] - static_cast<int>(s.size());
print_blanks(blanks, m_out);
m_out << s << " "; // the column interval
}
m_out << std::endl;
}
template <typename T, typename X> std::string core_solver_pretty_printer<T, X>::get_low_bound_string(unsigned j) {
switch (m_core_solver.get_column_type(j)){
case column_type::boxed:
case column_type::low_bound:
case column_type::fixed:
if (m_core_solver.low_bounds_are_set())
return T_to_string(m_core_solver.low_bound_value(j));
else
return std::string("0");
break;
default:
return std::string();
}
}
template <typename T, typename X> std::string core_solver_pretty_printer<T, X>::get_upp_bound_string(unsigned j) {
switch (m_core_solver.get_column_type(j)){
case column_type::boxed:
case column_type::upper_bound:
case column_type::fixed:
return T_to_string(m_core_solver.upper_bound_value(j));
break;
default:
return std::string();
}
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::print_lows() {
if (ncols() == 0) {
return;
}
int blanks = m_title_width + 1 - static_cast<unsigned>(m_low_bounds_title.size());
m_out << m_low_bounds_title;
print_blanks(blanks, m_out);
for (unsigned i = 0; i < ncols(); i++) {
string s = get_low_bound_string(i);
int blanks = m_column_widths[i] - static_cast<unsigned>(s.size());
print_blanks(blanks, m_out);
m_out << s << " "; // the column interval
}
m_out << std::endl;
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::print_upps() {
if (ncols() == 0) {
return;
}
int blanks = m_title_width + 1 - static_cast<unsigned>(m_upp_bounds_title.size());
m_out << m_upp_bounds_title;
print_blanks(blanks, m_out);
for (unsigned i = 0; i < ncols(); i++) {
string s = get_upp_bound_string(i);
int blanks = m_column_widths[i] - static_cast<unsigned>(s.size());
print_blanks(blanks, m_out);
m_out << s << " "; // the column interval
}
m_out << std::endl;
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::print_exact_norms() {
if (m_core_solver.use_tableau()) return;
int blanks = m_title_width + 1 - static_cast<int>(m_exact_norm_title.size());
m_out << m_exact_norm_title;
print_blanks(blanks, m_out);
for (unsigned i = 0; i < ncols(); i++) {
string s = get_exact_column_norm_string(i);
int blanks = m_column_widths[i] - static_cast<int>(s.size());
print_blanks(blanks, m_out);
m_out << s << " ";
}
m_out << std::endl;
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::print_approx_norms() {
if (m_core_solver.use_tableau()) return;
int blanks = m_title_width + 1 - static_cast<int>(m_approx_norm_title.size());
m_out << m_approx_norm_title;
print_blanks(blanks, m_out);
for (unsigned i = 0; i < ncols(); i++) {
string s = T_to_string(m_core_solver.m_column_norms[i]);
int blanks = m_column_widths[i] - static_cast<int>(s.size());
print_blanks(blanks, m_out);
m_out << s << " ";
}
m_out << std::endl;
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::print() {
for (unsigned i = 0; i < nrows(); i++) {
print_row(i);
}
print_bottom_line();
print_cost();
print_x();
print_basis_heading();
print_lows();
print_upps();
print_exact_norms();
if (m_core_solver.m_column_norms.size() > 0)
print_approx_norms();
m_out << std::endl;
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::print_basis_heading() {
int blanks = m_title_width + 1 - static_cast<int>(m_basis_heading_title.size());
m_out << m_basis_heading_title;
print_blanks(blanks, m_out);
if (ncols() == 0) {
return;
}
auto bh = m_core_solver.m_basis_heading;
for (unsigned i = 0; i < ncols(); i++) {
string s = T_to_string(bh[i]);
int blanks = m_column_widths[i] - static_cast<unsigned>(s.size());
print_blanks(blanks, m_out);
m_out << s << " "; // the column interval
}
m_out << std::endl;
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::print_cost() {
int blanks = m_title_width + 1 - static_cast<int>(m_cost_title.size());
m_out << m_cost_title;
print_blanks(blanks, m_out);
print_given_rows(m_costs, m_cost_signs, m_core_solver.get_cost());
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::print_given_rows(vector<string> & row, vector<string> & signs, X rst) {
for (unsigned col = 0; col < row.size(); col++) {
unsigned width = m_column_widths[col];
string s = row[col];
int number_of_blanks = width - static_cast<unsigned>(s.size());
lean_assert(number_of_blanks >= 0);
print_blanks(number_of_blanks, m_out);
m_out << s << ' ';
if (col < row.size() - 1) {
m_out << signs[col + 1] << ' ';
}
}
m_out << '=';
string rs = T_to_string(rst);
int nb = m_rs_width - static_cast<int>(rs.size());
lean_assert(nb >= 0);
print_blanks(nb + 1, m_out);
m_out << rs << std::endl;
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::print_row(unsigned i){
print_blanks(m_title_width + 1, m_out);
auto row = m_A[i];
auto sign_row = m_signs[i];
auto rs = m_rs[i];
print_given_rows(row, sign_row, rs);
}
}

15
src/util/lp/core_solver_pretty_printer_instances.cpp Normal file
View file

@ -0,0 +1,15 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#include "util/lp/numeric_pair.h"
#include "util/lp/core_solver_pretty_printer.hpp"
template lean::core_solver_pretty_printer<double, double>::core_solver_pretty_printer(lean::lp_core_solver_base<double, double> &, std::ostream & out);
template void lean::core_solver_pretty_printer<double, double>::print();
template lean::core_solver_pretty_printer<double, double>::~core_solver_pretty_printer();
template lean::core_solver_pretty_printer<lean::mpq, lean::mpq>::core_solver_pretty_printer(lean::lp_core_solver_base<lean::mpq, lean::mpq> &, std::ostream & out);
template void lean::core_solver_pretty_printer<lean::mpq, lean::mpq>::print();
template lean::core_solver_pretty_printer<lean::mpq, lean::mpq>::~core_solver_pretty_printer();
template lean::core_solver_pretty_printer<lean::mpq, lean::numeric_pair<lean::mpq> >::core_solver_pretty_printer(lean::lp_core_solver_base<lean::mpq, lean::numeric_pair<lean::mpq> > &, std::ostream & out);
template lean::core_solver_pretty_printer<lean::mpq, lean::numeric_pair<lean::mpq> >::~core_solver_pretty_printer();
template void lean::core_solver_pretty_printer<lean::mpq, lean::numeric_pair<lean::mpq> >::print();

92
src/util/lp/dense_matrix.h Normal file
View file

@ -0,0 +1,92 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#pragma once
#ifdef LEAN_DEBUG
#include "util/vector.h"
#include "util/lp/matrix.h"
namespace lean {
// used for debugging purposes only
template <typename T, typename X>
class dense_matrix: public matrix<T, X> {
public:
struct ref {
unsigned m_i;
dense_matrix & m_s;
ref(unsigned i, dense_matrix & s) :m_i(i * s.m_n), m_s(s){}
T & operator[] (unsigned j) {
return m_s.m_values[m_i + j];
}
const T & operator[] (unsigned j) const {
return m_s.m_v[m_i + j];
}
};
ref operator[] (unsigned i) {
return ref(i, *this);
}
unsigned m_m; // number of rows
unsigned m_n; // number of const
vector<T> m_values;
dense_matrix(unsigned m, unsigned n);
dense_matrix operator*=(matrix<T, X> const & a) {
lean_assert(column_count() == a.row_count());
dense_matrix c(row_count(), a.column_count());
for (unsigned i = 0; i < row_count(); i++) {
for (unsigned j = 0; j < a.column_count(); j++) {
T v = numeric_traits<T>::zero();
for (unsigned k = 0; k < a.column_count(); k++) {
v += get_elem(i, k) * a(k, j);
}
c.set_elem(i, j, v);
}
}
*this = c;
return *this;
}
dense_matrix & operator=(matrix<T, X> const & other);
dense_matrix & operator=(dense_matrix const & other);
dense_matrix(matrix<T, X> const * other);
void apply_from_right(T * w);
void apply_from_right(vector <T> & w);
T * apply_from_left_with_different_dims(vector<T> & w);
void apply_from_left(vector<T> & w , lp_settings & ) { apply_from_left(w); }
void apply_from_left(vector<T> & w);
void apply_from_left(X * w, lp_settings & );
void apply_from_left_to_X(vector<X> & w, lp_settings & );
virtual void set_number_of_rows(unsigned /*m*/) {}
virtual void set_number_of_columns(unsigned /*n*/) { }
T get_elem(unsigned i, unsigned j) const { return m_values[i * m_n + j]; }
unsigned row_count() const { return m_m; }
unsigned column_count() const { return m_n; }
void set_elem(unsigned i, unsigned j, const T& val) { m_values[i * m_n + j] = val; }
// This method pivots row i to row i0 by muliplying row i by
// alpha and adding it to row i0.
void pivot_row_to_row(unsigned i, const T& alpha, unsigned i0,
const double & pivot_epsilon);
void swap_columns(unsigned a, unsigned b);
void swap_rows(unsigned a, unsigned b);
void multiply_row_by_constant(unsigned row, T & t);
};
template <typename T, typename X>
dense_matrix<T, X> operator* (matrix<T, X> & a, matrix<T, X> & b);
}
#endif

186
src/util/lp/dense_matrix.hpp Normal file
View file

@ -0,0 +1,186 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#include "util/lp/lp_settings.h"
#ifdef LEAN_DEBUG
#include "util/vector.h"
#include "util/lp/numeric_pair.h"
#include "util/lp/dense_matrix.h"
namespace lean {
template <typename T> void print_vector(const vector<T> & t, std::ostream & out);
template <typename T, typename X> dense_matrix<T, X>::dense_matrix(unsigned m, unsigned n) : m_m(m), m_n(n), m_values(m * n, numeric_traits<T>::zero()) {
}
template <typename T, typename X> dense_matrix<T, X>&
dense_matrix<T, X>::operator=(matrix<T, X> const & other){
if ( this == & other)
return *this;
m_values = new T[m_m * m_n];
for (unsigned i = 0; i < m_m; i ++)
for (unsigned j = 0; j < m_n; j++)
m_values[i * m_n + j] = other.get_elem(i, j);
return *this;
}
template <typename T, typename X> dense_matrix<T, X>&
dense_matrix<T, X>::operator=(dense_matrix const & other){
if ( this == & other)
return *this;
m_m = other.m_m;
m_n = other.m_n;
m_values.resize(m_m * m_n);
for (unsigned i = 0; i < m_m; i ++)
for (unsigned j = 0; j < m_n; j++)
m_values[i * m_n + j] = other.get_elem(i, j);
return *this;
}
template <typename T, typename X> dense_matrix<T, X>::dense_matrix(matrix<T, X> const * other) :
m_m(other->row_count()),
m_n(other->column_count()) {
m_values.resize(m_m*m_n);
for (unsigned i = 0; i < m_m; i++)
for (unsigned j = 0; j < m_n; j++)
m_values[i * m_n + j] = other->get_elem(i, j);
}
template <typename T, typename X> void dense_matrix<T, X>::apply_from_right(T * w) {
T * t = new T[m_m];
for (int i = 0; i < m_m; i ++) {
T v = numeric_traits<T>::zero();
for (int j = 0; j < m_m; j++) {
v += w[j]* get_elem(j, i);
}
t[i] = v;
}
for (int i = 0; i < m_m; i++) {
w[i] = t[i];
}
delete [] t;
}
template <typename T, typename X> void dense_matrix<T, X>::apply_from_right(vector <T> & w) {
vector<T> t(m_m, numeric_traits<T>::zero());
for (unsigned i = 0; i < m_m; i ++) {
auto & v = t[i];
for (unsigned j = 0; j < m_m; j++)
v += w[j]* get_elem(j, i);
}
for (unsigned i = 0; i < m_m; i++)
w[i] = t[i];
}
template <typename T, typename X> T* dense_matrix<T, X>::
apply_from_left_with_different_dims(vector<T> & w) {
T * t = new T[m_m];
for (int i = 0; i < m_m; i ++) {
T v = numeric_traits<T>::zero();
for (int j = 0; j < m_n; j++) {
v += w[j]* get_elem(i, j);
}
t[i] = v;
}
return t;
}
template <typename T, typename X> void dense_matrix<T, X>::apply_from_left(vector<T> & w) {
T * t = new T[m_m];
for (unsigned i = 0; i < m_m; i ++) {
T v = numeric_traits<T>::zero();
for (unsigned j = 0; j < m_m; j++) {
v += w[j]* get_elem(i, j);
}
t[i] = v;
}
for (unsigned i = 0; i < m_m; i ++) {
w[i] = t[i];
}
delete [] t;
}
template <typename T, typename X> void dense_matrix<T, X>::apply_from_left(X * w, lp_settings & ) {
T * t = new T[m_m];
for (int i = 0; i < m_m; i ++) {
T v = numeric_traits<T>::zero();
for (int j = 0; j < m_m; j++) {
v += w[j]* get_elem(i, j);
}
t[i] = v;
}
for (int i = 0; i < m_m; i ++) {
w[i] = t[i];
}
delete [] t;
}
template <typename T, typename X> void dense_matrix<T, X>::apply_from_left_to_X(vector<X> & w, lp_settings & ) {
vector<X> t(m_m);
for (int i = 0; i < m_m; i ++) {
X v = zero_of_type<X>();
for (int j = 0; j < m_m; j++) {
v += w[j]* get_elem(i, j);
}
t[i] = v;
}
for (int i = 0; i < m_m; i ++) {
w[i] = t[i];
}
}
// This method pivots row i to row i0 by muliplying row i by
// alpha and adding it to row i0.
template <typename T, typename X> void dense_matrix<T, X>::pivot_row_to_row(unsigned i, const T& alpha, unsigned i0,
const double & pivot_epsilon) {
for (unsigned j = 0; j < m_n; j++) {
m_values[i0 * m_n + j] += m_values[i * m_n + j] * alpha;
if (fabs(m_values[i0 + m_n + j]) < pivot_epsilon) {
m_values[i0 + m_n + j] = numeric_traits<T>::zero();;
}
}
}
template <typename T, typename X> void dense_matrix<T, X>::swap_columns(unsigned a, unsigned b) {
for (unsigned i = 0; i < m_m; i++) {
T t = get_elem(i, a);
set_elem(i, a, get_elem(i, b));
set_elem(i, b, t);
}
}
template <typename T, typename X> void dense_matrix<T, X>::swap_rows(unsigned a, unsigned b) {
for (unsigned i = 0; i < m_n; i++) {
T t = get_elem(a, i);
set_elem(a, i, get_elem(b, i));
set_elem(b, i, t);
}
}
template <typename T, typename X> void dense_matrix<T, X>::multiply_row_by_constant(unsigned row, T & t) {
for (unsigned i = 0; i < m_n; i++) {
set_elem(row, i, t * get_elem(row, i));
}
}
template <typename T, typename X>
dense_matrix<T, X> operator* (matrix<T, X> & a, matrix<T, X> & b){
lean_assert(a.column_count() == b.row_count());
dense_matrix<T, X> ret(a.row_count(), b.column_count());
for (unsigned i = 0; i < ret.m_m; i++)
for (unsigned j = 0; j< ret.m_n; j++) {
T v = numeric_traits<T>::zero();
for (unsigned k = 0; k < a.column_count(); k ++){
v += (a.get_elem(i, k) * b.get_elem(k, j));
}
ret.set_elem(i, j, v);
}
return ret;
}
}
#endif

24
src/util/lp/dense_matrix_instances.cpp Normal file
View file

@ -0,0 +1,24 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#include "util/lp/lp_settings.h"
#include "util/lp/dense_matrix.hpp"
#ifdef LEAN_DEBUG
#include "util/vector.h"
template lean::dense_matrix<double, double> lean::operator*<double, double>(lean::matrix<double, double>&, lean::matrix<double, double>&);
template void lean::dense_matrix<double, double>::apply_from_left(vector<double> &);
template lean::dense_matrix<double, double>::dense_matrix(lean::matrix<double, double> const*);
template lean::dense_matrix<double, double>::dense_matrix(unsigned int, unsigned int);
template lean::dense_matrix<double, double>& lean::dense_matrix<double, double>::operator=(lean::dense_matrix<double, double> const&);
template lean::dense_matrix<lean::mpq, lean::numeric_pair<lean::mpq> >::dense_matrix(lean::matrix<lean::mpq, lean::numeric_pair<lean::mpq> > const*);
template void lean::dense_matrix<lean::mpq, lean::numeric_pair<lean::mpq> >::apply_from_left(vector<lean::mpq>&);
template lean::dense_matrix<lean::mpq, lean::mpq> lean::operator*<lean::mpq, lean::mpq>(lean::matrix<lean::mpq, lean::mpq>&, lean::matrix<lean::mpq, lean::mpq>&);
template lean::dense_matrix<lean::mpq, lean::mpq> & lean::dense_matrix<lean::mpq, lean::mpq>::operator=(lean::dense_matrix<lean::mpq, lean::mpq> const&);
template lean::dense_matrix<lean::mpq, lean::numeric_pair<lean::mpq> >::dense_matrix(unsigned int, unsigned int);
template lean::dense_matrix<lean::mpq, lean::numeric_pair<lean::mpq> >& lean::dense_matrix<lean::mpq, lean::numeric_pair<lean::mpq> >::operator=(lean::dense_matrix<lean::mpq, lean::numeric_pair<lean::mpq> > const&);
template lean::dense_matrix<lean::mpq, lean::numeric_pair<lean::mpq> > lean::operator*<lean::mpq, lean::numeric_pair<lean::mpq> >(lean::matrix<lean::mpq, lean::numeric_pair<lean::mpq> >&, lean::matrix<lean::mpq, lean::numeric_pair<lean::mpq> >&);
template void lean::dense_matrix<lean::mpq, lean::numeric_pair< lean::mpq> >::apply_from_right( vector< lean::mpq> &);
template void lean::dense_matrix<double,double>::apply_from_right(class vector<double> &);
template void lean::dense_matrix<lean::mpq, lean::mpq>::apply_from_left(vector<lean::mpq>&);
#endif

83
src/util/lp/eta_matrix.h Normal file
View file

@ -0,0 +1,83 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#pragma once
#include "util/vector.h"
#include "util/lp/tail_matrix.h"
#include "util/lp/permutation_matrix.h"
namespace lean {
// This is the sum of a unit matrix and a one-column matrix
template <typename T, typename X>
class eta_matrix
: public tail_matrix<T, X> {
#ifdef LEAN_DEBUG
unsigned m_length;
#endif
unsigned m_column_index;
public:
sparse_vector<T> m_column_vector;
T m_diagonal_element;
#ifdef LEAN_DEBUG
eta_matrix(unsigned column_index, unsigned length):
#else
eta_matrix(unsigned column_index):
#endif
#ifdef LEAN_DEBUG
m_length(length),
#endif
m_column_index(column_index) {}
bool is_dense() const { return false; }
void print(std::ostream & out) {
print_matrix(*this, out);
}
bool is_unit() {
return m_column_vector.size() == 0 && m_diagonal_element == 1;
}
bool set_diagonal_element(T const & diagonal_element) {
m_diagonal_element = diagonal_element;
return !lp_settings::is_eps_small_general(diagonal_element, 1e-12);
}
const T & get_diagonal_element() const {
return m_diagonal_element;
}
void apply_from_left(vector<X> & w, lp_settings & );
template <typename L>
void apply_from_left_local(indexed_vector<L> & w, lp_settings & settings);
void apply_from_left_to_T(indexed_vector<T> & w, lp_settings & settings) {
apply_from_left_local(w, settings);
}
void push_back(unsigned row_index, T val ) {
lean_assert(row_index != m_column_index);
m_column_vector.push_back(row_index, val);
}
void apply_from_right(vector<T> & w);
void apply_from_right(indexed_vector<T> & w);
T get_elem(unsigned i, unsigned j) const;
#ifdef LEAN_DEBUG
unsigned row_count() const { return m_length; }
unsigned column_count() const { return m_length; }
void set_number_of_rows(unsigned m) { m_length = m; }
void set_number_of_columns(unsigned n) { m_length = n; }
#endif
void divide_by_diagonal_element() {
m_column_vector.divide(m_diagonal_element);
}
void conjugate_by_permutation(permutation_matrix<T, X> & p);
};
}

136
src/util/lp/eta_matrix.hpp Normal file
View file

@ -0,0 +1,136 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#pragma once
#include "util/vector.h"
#include "util/lp/eta_matrix.h"
namespace lean {
// This is the sum of a unit matrix and a one-column matrix
template <typename T, typename X>
void eta_matrix<T, X>::apply_from_left(vector<X> & w, lp_settings & ) {
auto & w_at_column_index = w[m_column_index];
for (auto & it : m_column_vector.m_data) {
w[it.first] += w_at_column_index * it.second;
}
w_at_column_index /= m_diagonal_element;
}
template <typename T, typename X>
template <typename L>
void eta_matrix<T, X>::
apply_from_left_local(indexed_vector<L> & w, lp_settings & settings) {
const L w_at_column_index = w[m_column_index];
if (is_zero(w_at_column_index)) return;
if (settings.abs_val_is_smaller_than_drop_tolerance(w[m_column_index] /= m_diagonal_element)) {
w[m_column_index] = zero_of_type<L>();
w.erase_from_index(m_column_index);
}
for (auto & it : m_column_vector.m_data) {
unsigned i = it.first;
if (is_zero(w[i])) {
L v = w[i] = w_at_column_index * it.second;
if (settings.abs_val_is_smaller_than_drop_tolerance(v)) {
w[i] = zero_of_type<L>();
continue;
}
w.m_index.push_back(i);
} else {
L v = w[i] += w_at_column_index * it.second;
if (settings.abs_val_is_smaller_than_drop_tolerance(v)) {
w[i] = zero_of_type<L>();
w.erase_from_index(i);
}
}
}
}
template <typename T, typename X>
void eta_matrix<T, X>::apply_from_right(vector<T> & w) {
#ifdef LEAN_DEBUG
// dense_matrix<T, X> deb(*this);
// auto clone_w = clone_vector<T>(w, get_number_of_rows());
// deb.apply_from_right(clone_w);
#endif
T t = w[m_column_index] / m_diagonal_element;
for (auto & it : m_column_vector.m_data) {
t += w[it.first] * it.second;
}
w[m_column_index] = t;
#ifdef LEAN_DEBUG
// lean_assert(vectors_are_equal<T>(clone_w, w, get_number_of_rows()));
// delete clone_w;
#endif
}
template <typename T, typename X>
void eta_matrix<T, X>::apply_from_right(indexed_vector<T> & w) {
if (w.m_index.size() == 0)
return;
#ifdef LEAN_DEBUG
// vector<T> wcopy(w.m_data);
// apply_from_right(wcopy);
#endif
T & t = w[m_column_index];
t /= m_diagonal_element;
bool was_in_index = (!numeric_traits<T>::is_zero(t));
for (auto & it : m_column_vector.m_data) {
t += w[it.first] * it.second;
}
if (numeric_traits<T>::precise() ) {
if (!numeric_traits<T>::is_zero(t)) {
if (!was_in_index)
w.m_index.push_back(m_column_index);
} else {
if (was_in_index)
w.erase_from_index(m_column_index);
}
} else {
if (!lp_settings::is_eps_small_general(t, 1e-14)) {
if (!was_in_index)
w.m_index.push_back(m_column_index);
} else {
if (was_in_index)
w.erase_from_index(m_column_index);
t = zero_of_type<T>();
}
}
#ifdef LEAN_DEBUG
// lean_assert(w.is_OK());
// lean_assert(vectors_are_equal<T>(wcopy, w.m_data));
#endif
}
#ifdef LEAN_DEBUG
template <typename T, typename X>
T eta_matrix<T, X>::get_elem(unsigned i, unsigned j) const {
if (j == m_column_index){
if (i == j) {
return 1 / m_diagonal_element;
}
return m_column_vector[i];
}
return i == j ? numeric_traits<T>::one() : numeric_traits<T>::zero();
}
#endif
template <typename T, typename X>
void eta_matrix<T, X>::conjugate_by_permutation(permutation_matrix<T, X> & p) {
// this = p * this * p(-1)
#ifdef LEAN_DEBUG
// auto rev = p.get_reverse();
// auto deb = ((*this) * rev);
// deb = p * deb;
#endif
m_column_index = p.get_rev(m_column_index);
for (auto & pair : m_column_vector.m_data) {
pair.first = p.get_rev(pair.first);
}
#ifdef LEAN_DEBUG
// lean_assert(deb == *this);
#endif
}
}

28
src/util/lp/eta_matrix_instances.cpp Normal file
View file

@ -0,0 +1,28 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#include <memory>
#include "util/vector.h"
#include "util/lp/numeric_pair.h"
#include "util/lp/eta_matrix.hpp"
#ifdef LEAN_DEBUG
template double lean::eta_matrix<double, double>::get_elem(unsigned int, unsigned int) const;
template lean::mpq lean::eta_matrix<lean::mpq, lean::mpq>::get_elem(unsigned int, unsigned int) const;
template lean::mpq lean::eta_matrix<lean::mpq, lean::numeric_pair<lean::mpq> >::get_elem(unsigned int, unsigned int) const;
#endif
template void lean::eta_matrix<double, double>::apply_from_left(vector<double>&, lean::lp_settings&);
template void lean::eta_matrix<double, double>::apply_from_right(vector<double>&);
template void lean::eta_matrix<double, double>::conjugate_by_permutation(lean::permutation_matrix<double, double>&);
template void lean::eta_matrix<lean::mpq, lean::mpq>::apply_from_left(vector<lean::mpq>&, lean::lp_settings&);
template void lean::eta_matrix<lean::mpq, lean::mpq>::apply_from_right(vector<lean::mpq>&);
template void lean::eta_matrix<lean::mpq, lean::mpq>::conjugate_by_permutation(lean::permutation_matrix<lean::mpq, lean::mpq>&);
template void lean::eta_matrix<lean::mpq, lean::numeric_pair<lean::mpq> >::apply_from_left(vector<lean::numeric_pair<lean::mpq> >&, lean::lp_settings&);
template void lean::eta_matrix<lean::mpq, lean::numeric_pair<lean::mpq> >::apply_from_right(vector<lean::mpq>&);
template void lean::eta_matrix<lean::mpq, lean::numeric_pair<lean::mpq> >::conjugate_by_permutation(lean::permutation_matrix<lean::mpq, lean::numeric_pair<lean::mpq> >&);
template void lean::eta_matrix<double, double>::apply_from_left_local<double>(lean::indexed_vector<double>&, lean::lp_settings&);
template void lean::eta_matrix<lean::mpq, lean::mpq>::apply_from_left_local<lean::mpq>(lean::indexed_vector<lean::mpq>&, lean::lp_settings&);
template void lean::eta_matrix<lean::mpq, lean::numeric_pair<lean::mpq> >::apply_from_left_local<lean::mpq>(lean::indexed_vector<lean::mpq>&, lean::lp_settings&);
template void lean::eta_matrix<lean::mpq, lean::numeric_pair<lean::mpq> >::apply_from_right(lean::indexed_vector<lean::mpq>&);
template void lean::eta_matrix<lean::mpq, lean::mpq>::apply_from_right(lean::indexed_vector<lean::mpq>&);
template void lean::eta_matrix<double, double>::apply_from_right(lean::indexed_vector<double>&);

39
src/util/lp/hash_helper.h Normal file
View file

@ -0,0 +1,39 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#pragma once
#include <utility>
#include <functional>
#include "util/numerics/mpq.h"
#ifdef __CLANG__
#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Wmismatched-tags"
#endif
namespace std {
template<>
struct hash<lean::mpq> {
inline size_t operator()(const lean::mpq & v) const {
return v.hash();
}
};
}
template <class T>
inline void hash_combine(std::size_t & seed, const T & v) {
seed ^= std::hash<T>()(v) + 0x9e3779b9 + (seed << 6) + (seed >> 2);
}
namespace std {
template<typename S, typename T> struct hash<pair<S, T>> {
inline size_t operator()(const pair<S, T> & v) const {
size_t seed = 0;
hash_combine(seed, v.first);
hash_combine(seed, v.second);
return seed;
}
};
}
#ifdef __CLANG__
#pragma clang diagnostic pop
#endif

42
src/util/lp/implied_bound.h Normal file
View file

@ -0,0 +1,42 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#pragma once
#include "util/lp/lp_settings.h"
#include "util/lp/lar_constraints.h"
namespace lean {
struct implied_bound {
mpq m_bound;
unsigned m_j; // the column for which the bound has been found
bool m_is_low_bound;
bool m_coeff_before_j_is_pos;
unsigned m_row_or_term_index;
bool m_strict;
lconstraint_kind kind() const {
lconstraint_kind k = m_is_low_bound? GE : LE;
if (m_strict)
k = static_cast<lconstraint_kind>(k / 2);
return k;
}
bool operator==(const implied_bound & o) const {
return m_j == o.m_j && m_is_low_bound == o.m_is_low_bound && m_bound == o.m_bound &&
m_coeff_before_j_is_pos == o.m_coeff_before_j_is_pos &&
m_row_or_term_index == o.m_row_or_term_index && m_strict == o.m_strict;
}
implied_bound(){}
implied_bound(const mpq & a,
unsigned j,
bool low_bound,
bool coeff_before_j_is_pos,
unsigned row_or_term_index,
bool strict):
m_bound(a),
m_j(j),
m_is_low_bound(low_bound),
m_coeff_before_j_is_pos(coeff_before_j_is_pos),
m_row_or_term_index(row_or_term_index),
m_strict(strict) {}
};
}

55
src/util/lp/indexed_value.h Normal file
View file

@ -0,0 +1,55 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#pragma once
namespace lean {
template <typename T>
class indexed_value {
public:
T m_value;
// the idea is that m_index for a row element gives its column, and for a column element its row
unsigned m_index;
// m_other point is the offset of the corresponding element in its vector : for a row element it point to the column element offset,
// for a column element it points to the row element offset
unsigned m_other;
indexed_value() {}
indexed_value(T v, unsigned i) : m_value(v), m_index(i) {}
indexed_value(T v, unsigned i, unsigned other) :
m_value(v), m_index(i), m_other(other) {
}
indexed_value(const indexed_value & iv) {
m_value = iv.m_value;
m_index = iv.m_index;
m_other = iv.m_other;
}
indexed_value & operator=(const indexed_value & right_side) {
m_value = right_side.m_value;
m_index = right_side.m_index;
m_other = right_side.m_other;
return *this;
}
const T & value() const {
return m_value;
}
void set_value(T val) {
m_value = val;
}
};
#ifdef LEAN_DEBUG
template <typename X>
bool check_vector_for_small_values(indexed_vector<X> & w, lp_settings & settings) {
for (unsigned i : w.m_index) {
const X & v = w[i];
if ((!is_zero(v)) && settings.abs_val_is_smaller_than_drop_tolerance(v))
return false;
}
return true;
}
#endif
}

169
src/util/lp/indexed_vector.h Normal file
View file

@ -0,0 +1,169 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#pragma once
#include "util/vector.h"
#include "util/debug.h"
#include <string>
#include <iomanip>
#include "util/lp/lp_utils.h"
#include "util/lp/lp_settings.h"
#include <unordered_set>
namespace lean {
template <typename T> void print_vector(const vector<T> & t, std::ostream & out);
template <typename T> void print_vector(const buffer<T> & t, std::ostream & out);
template <typename T> void print_sparse_vector(const vector<T> & t, std::ostream & out);
void print_vector(const vector<mpq> & t, std::ostream & out);
template <typename T>
class indexed_vector {
public:
// m_index points to non-zero elements of m_data
vector<T> m_data;
vector<unsigned> m_index;
indexed_vector(unsigned data_size) {
m_data.resize(data_size, numeric_traits<T>::zero());
}
indexed_vector& operator=(const indexed_vector<T>& y) {
for (unsigned i: m_index)
m_data[i] = zero_of_type<T>();
m_index = y.m_index;
m_data.resize(y.data_size());
for (unsigned i : m_index)
m_data[i] = y[i];
return *this;
}
bool operator==(const indexed_vector<T>& y) const {
std::unordered_set<unsigned> y_index;
for (unsigned i : y.m_index)
y_index.insert(i);
std::unordered_set<unsigned> this_index;
for (unsigned i : m_index)
this_index.insert(i);
for (unsigned i : y.m_index) {
if (this_index.find(i) == this_index.end())
return false;
}
for (unsigned i : m_index) {
if (y_index.find(i) == y_index.end())
return false;
}
return vectors_are_equal(m_data, m_data);
}
indexed_vector() {}
void resize(unsigned data_size);
unsigned data_size() const {
return m_data.size();
}
unsigned size() {
return m_index.size();
}
void set_value(const T& value, unsigned index);
void set_value_as_in_dictionary(unsigned index) {
lean_assert(index < m_data.size());
T & loc = m_data[index];
if (is_zero(loc)) {
m_index.push_back(index);
loc = one_of_type<T>(); // use as a characteristic function
}
}
void clear();
void clear_all();
const T& operator[] (unsigned i) const {
return m_data[i];
}
T& operator[] (unsigned i) {
return m_data[i];
}
void clean_up() {
#if 0==1
for (unsigned k = 0; k < m_index.size(); k++) {
unsigned i = m_index[k];
T & v = m_data[i];
if (lp_settings::is_eps_small_general(v, 1e-14)) {
v = zero_of_type<T>();
m_index.erase(m_index.begin() + k--);
}
}
#endif
vector<unsigned> index_copy;
for (unsigned i : m_index) {
T & v = m_data[i];
if (!lp_settings::is_eps_small_general(v, 1e-14)) {
index_copy.push_back(i);
} else if (!numeric_traits<T>::is_zero(v)) {
v = zero_of_type<T>();
}
}
m_index = index_copy;
}
void erase_from_index(unsigned j);
void add_value_at_index_with_drop_tolerance(unsigned j, const T& val_to_add) {
T & v = m_data[j];
bool was_zero = is_zero(v);
v += val_to_add;
if (lp_settings::is_eps_small_general(v, 1e-14)) {
v = zero_of_type<T>();
if (!was_zero) {
erase_from_index(j);
}
} else {
if (was_zero)
m_index.push_back(j);
}
}
void add_value_at_index(unsigned j, const T& val_to_add) {
T & v = m_data[j];
bool was_zero = is_zero(v);
v += val_to_add;
if (is_zero(v)) {
if (!was_zero)
erase_from_index(j);
} else {
if (was_zero)
m_index.push_back(j);
}
}
void restore_index_and_clean_from_data() {
m_index.resize(0);
for (unsigned i = 0; i < m_data.size(); i++) {
T & v = m_data[i];
if (lp_settings::is_eps_small_general(v, 1e-14)) {
v = zero_of_type<T>();
} else {
m_index.push_back(i);
}
}
}
#ifdef LEAN_DEBUG
bool is_OK() const;
void print(std::ostream & out);
#endif
};
}

110
src/util/lp/indexed_vector.hpp Normal file
View file

@ -0,0 +1,110 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#include "util/vector.h"
#include "util/lp/indexed_vector.h"
#include "util/lp/lp_settings.h"
namespace lean {
template <typename T>
void print_vector(const vector<T> & t, std::ostream & out) {
for (unsigned i = 0; i < t.size(); i++)
out << t[i] << " ";
out << std::endl;
}
template <typename T>
void print_vector(const buffer<T> & t, std::ostream & out) {
for (unsigned i = 0; i < t.size(); i++)
out << t[i] << " ";
out << std::endl;
}
template <typename T>
void print_sparse_vector(const vector<T> & t, std::ostream & out) {
for (unsigned i = 0; i < t.size(); i++) {
if (is_zero(t[i]))continue;
out << "[" << i << "] = " << t[i] << ", ";
}
out << std::endl;
}
void print_vector(const vector<mpq> & t, std::ostream & out) {
for (unsigned i = 0; i < t.size(); i++)
out << t[i].get_double() << std::setprecision(3) << " ";
out << std::endl;
}
template <typename T>
void indexed_vector<T>::resize(unsigned data_size) {
clear();
m_data.resize(data_size, numeric_traits<T>::zero());
lean_assert(is_OK());
}
template <typename T>
void indexed_vector<T>::set_value(const T& value, unsigned index) {
m_data[index] = value;
lean_assert(std::find(m_index.begin(), m_index.end(), index) == m_index.end());
m_index.push_back(index);
}
template <typename T>
void indexed_vector<T>::clear() {
for (unsigned i : m_index)
m_data[i] = numeric_traits<T>::zero();
m_index.resize(0);
}
template <typename T>
void indexed_vector<T>::clear_all() {
unsigned i = m_data.size();
while (i--) m_data[i] = numeric_traits<T>::zero();
m_index.resize(0);
}
template <typename T>
void indexed_vector<T>::erase_from_index(unsigned j) {
auto it = std::find(m_index.begin(), m_index.end(), j);
if (it != m_index.end())
m_index.erase(it);
}
#ifdef LEAN_DEBUG
template <typename T>
bool indexed_vector<T>::is_OK() const {
return true;
const double drop_eps = 1e-14;
for (unsigned i = 0; i < m_data.size(); i++) {
if (!is_zero(m_data[i]) && lp_settings::is_eps_small_general(m_data[i], drop_eps)) {
return false;
}
if (lp_settings::is_eps_small_general(m_data[i], drop_eps) != (std::find(m_index.begin(), m_index.end(), i) == m_index.end())) {
return false;
}
}
std::unordered_set<unsigned> s;
for (unsigned i : m_index) {
//no duplicates!!!
if (s.find(i) != s.end())
return false;
s.insert(i);
if (i >= m_data.size())
return false;
}
return true;
}
template <typename T>
void indexed_vector<T>::print(std::ostream & out) {
out << "m_index " << std::endl;
for (unsigned i = 0; i < m_index.size(); i++) {
out << m_index[i] << " ";
}
out << std::endl;
print_vector(m_data, out);
}
#endif
}

36
src/util/lp/indexed_vector_instances.cpp Normal file
View file

@ -0,0 +1,36 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#include "util/vector.h"
#include "util/lp/indexed_vector.hpp"
namespace lean {
template void indexed_vector<double>::clear();
template void indexed_vector<double>::clear_all();
template void indexed_vector<double>::erase_from_index(unsigned int);
template void indexed_vector<double>::set_value(const double&, unsigned int);
template void indexed_vector<mpq>::clear();
template void indexed_vector<unsigned>::clear();
template void indexed_vector<mpq>::clear_all();
template void indexed_vector<mpq>::erase_from_index(unsigned int);
template void indexed_vector<mpq>::resize(unsigned int);
template void indexed_vector<unsigned>::resize(unsigned int);
template void indexed_vector<mpq>::set_value(const mpq&, unsigned int);
template void indexed_vector<unsigned>::set_value(const unsigned&, unsigned int);
#ifdef LEAN_DEBUG
template bool indexed_vector<double>::is_OK() const;
template bool indexed_vector<mpq>::is_OK() const;
template bool indexed_vector<lean::numeric_pair<mpq> >::is_OK() const;
template void lean::indexed_vector< lean::mpq>::print(std::basic_ostream<char,struct std::char_traits<char> > &);
template void lean::indexed_vector<double>::print(std::basic_ostream<char,struct std::char_traits<char> > &);
template void lean::indexed_vector<lean::numeric_pair<lean::mpq> >::print(std::ostream&);
#endif
}
template void lean::print_vector<double>(vector<double> const&, std::ostream&);
template void lean::print_vector<unsigned int>(vector<unsigned int> const&, std::ostream&);
template void lean::print_vector<std::string>(vector<std::string> const&, std::ostream&);
template void lean::print_vector<lean::numeric_pair<lean::mpq> >(vector<lean::numeric_pair<lean::mpq>> const&, std::ostream&);
template void lean::indexed_vector<double>::resize(unsigned int);
template void lean::print_vector< lean::mpq>(vector< lean::mpq> const &, std::basic_ostream<char, std::char_traits<char> > &);
template void lean::print_vector<std::pair<lean::mpq, unsigned int> >(vector<std::pair<lean::mpq, unsigned int>> const&, std::ostream&);
template void lean::indexed_vector<lean::numeric_pair<lean::mpq> >::erase_from_index(unsigned int);

576
src/util/lp/init_lar_solver.h Normal file
View file

@ -0,0 +1,576 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
// here we are inside lean::lar_solver class
bool strategy_is_undecided() const {
return m_settings.simplex_strategy() == simplex_strategy_enum::undecided;
}
var_index add_var(unsigned ext_j) {
var_index i;
lean_assert (ext_j < m_terms_start_index);
if (ext_j >= m_terms_start_index)
throw 0; // todo : what is the right way to exit?
if (try_get_val(m_ext_vars_to_columns, ext_j, i)) {
return i;
}
lean_assert(m_vars_to_ul_pairs.size() == A_r().column_count());
i = A_r().column_count();
m_vars_to_ul_pairs.push_back (ul_pair(static_cast<unsigned>(-1)));
add_non_basic_var_to_core_fields(ext_j);
lean_assert(sizes_are_correct());
return i;
}
void register_new_ext_var_index(unsigned ext_v) {
lean_assert(!contains(m_ext_vars_to_columns, ext_v));
unsigned j = static_cast<unsigned>(m_ext_vars_to_columns.size());
m_ext_vars_to_columns[ext_v] = j;
lean_assert(m_columns_to_ext_vars_or_term_indices.size() == j);
m_columns_to_ext_vars_or_term_indices.push_back(ext_v);
}
void add_non_basic_var_to_core_fields(unsigned ext_j) {
register_new_ext_var_index(ext_j);
m_mpq_lar_core_solver.m_column_types.push_back(column_type::free_column);
m_columns_with_changed_bound.increase_size_by_one();
add_new_var_to_core_fields_for_mpq(false);
if (use_lu())
add_new_var_to_core_fields_for_doubles(false);
}
void add_new_var_to_core_fields_for_doubles(bool register_in_basis) {
unsigned j = A_d().column_count();
A_d().add_column();
lean_assert(m_mpq_lar_core_solver.m_d_x.size() == j);
// lean_assert(m_mpq_lar_core_solver.m_d_low_bounds.size() == j && m_mpq_lar_core_solver.m_d_upper_bounds.size() == j); // restore later
m_mpq_lar_core_solver.m_d_x.resize(j + 1 );
m_mpq_lar_core_solver.m_d_low_bounds.resize(j + 1);
m_mpq_lar_core_solver.m_d_upper_bounds.resize(j + 1);
lean_assert(m_mpq_lar_core_solver.m_d_heading.size() == j); // as A().column_count() on the entry to the method
if (register_in_basis) {
A_d().add_row();
m_mpq_lar_core_solver.m_d_heading.push_back(m_mpq_lar_core_solver.m_d_basis.size());
m_mpq_lar_core_solver.m_d_basis.push_back(j);
}else {
m_mpq_lar_core_solver.m_d_heading.push_back(- static_cast<int>(m_mpq_lar_core_solver.m_d_nbasis.size()) - 1);
m_mpq_lar_core_solver.m_d_nbasis.push_back(j);
}
}
void add_new_var_to_core_fields_for_mpq(bool register_in_basis) {
unsigned j = A_r().column_count();
A_r().add_column();
lean_assert(m_mpq_lar_core_solver.m_r_x.size() == j);
// lean_assert(m_mpq_lar_core_solver.m_r_low_bounds.size() == j && m_mpq_lar_core_solver.m_r_upper_bounds.size() == j); // restore later
m_mpq_lar_core_solver.m_r_x.resize(j + 1);
m_mpq_lar_core_solver.m_r_low_bounds.increase_size_by_one();
m_mpq_lar_core_solver.m_r_upper_bounds.increase_size_by_one();
m_mpq_lar_core_solver.m_r_solver.m_inf_set.increase_size_by_one();
m_mpq_lar_core_solver.m_r_solver.m_costs.resize(j + 1);
m_mpq_lar_core_solver.m_r_solver.m_d.resize(j + 1);
lean_assert(m_mpq_lar_core_solver.m_r_heading.size() == j); // as A().column_count() on the entry to the method
if (register_in_basis) {
A_r().add_row();
m_mpq_lar_core_solver.m_r_heading.push_back(m_mpq_lar_core_solver.m_r_basis.size());
m_mpq_lar_core_solver.m_r_basis.push_back(j);
if (m_settings.bound_propagation())
m_rows_with_changed_bounds.insert(A_r().row_count() - 1);
} else {
m_mpq_lar_core_solver.m_r_heading.push_back(- static_cast<int>(m_mpq_lar_core_solver.m_r_nbasis.size()) - 1);
m_mpq_lar_core_solver.m_r_nbasis.push_back(j);
}
}
var_index add_term_undecided(const vector<std::pair<mpq, var_index>> & coeffs,
const mpq &m_v) {
m_terms.push_back(new lar_term(coeffs, m_v));
m_orig_terms.push_back(new lar_term(coeffs, m_v));
return m_terms_start_index + m_terms.size() - 1;
}
// terms
var_index add_term(const vector<std::pair<mpq, var_index>> & coeffs,
const mpq &m_v) {
if (strategy_is_undecided())
return add_term_undecided(coeffs, m_v);
m_terms.push_back(new lar_term(coeffs, m_v));
m_orig_terms.push_back(new lar_term(coeffs, m_v));
unsigned adjusted_term_index = m_terms.size() - 1;
var_index ret = m_terms_start_index + adjusted_term_index;
if (use_tableau() && !coeffs.empty()) {
add_row_for_term(m_orig_terms.back(), ret);
if (m_settings.bound_propagation())
m_rows_with_changed_bounds.insert(A_r().row_count() - 1);
}
lean_assert(m_ext_vars_to_columns.size() == A_r().column_count());
return ret;
}
void add_row_for_term(const lar_term * term, unsigned term_ext_index) {
lean_assert(sizes_are_correct());
add_row_from_term_no_constraint(term, term_ext_index);
lean_assert(sizes_are_correct());
}
void add_row_from_term_no_constraint(const lar_term * term, unsigned term_ext_index) {
register_new_ext_var_index(term_ext_index);
// j will be a new variable
unsigned j = A_r().column_count();
ul_pair ul(j);
m_vars_to_ul_pairs.push_back(ul);
add_basic_var_to_core_fields();
if (use_tableau()) {
auto it = iterator_on_term_with_basis_var(*term, j);
A_r().fill_last_row_with_pivoting(it,
m_mpq_lar_core_solver.m_r_solver.m_basis_heading);
m_mpq_lar_core_solver.m_r_solver.m_b.resize(A_r().column_count(), zero_of_type<mpq>());
} else {
fill_last_row_of_A_r(A_r(), term);
}
m_mpq_lar_core_solver.m_r_x[j] = get_basic_var_value_from_row_directly(A_r().row_count() - 1);
if (use_lu())
fill_last_row_of_A_d(A_d(), term);
}
void add_basic_var_to_core_fields() {
bool use_lu = m_mpq_lar_core_solver.need_to_presolve_with_double_solver();
lean_assert(!use_lu || A_r().column_count() == A_d().column_count());
m_mpq_lar_core_solver.m_column_types.push_back(column_type::free_column);
m_columns_with_changed_bound.increase_size_by_one();
m_rows_with_changed_bounds.increase_size_by_one();
add_new_var_to_core_fields_for_mpq(true);
if (use_lu)
add_new_var_to_core_fields_for_doubles(true);
}
constraint_index add_var_bound(var_index j, lconstraint_kind kind, const mpq & right_side) {
constraint_index ci = m_constraints.size();
if (!is_term(j)) { // j is a var
auto vc = new lar_var_constraint(j, kind, right_side);
m_constraints.push_back(vc);
update_column_type_and_bound(j, kind, right_side, ci);
} else {
add_var_bound_on_constraint_for_term(j, kind, right_side, ci);
}
lean_assert(sizes_are_correct());
return ci;
}
void update_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_index) {
switch(m_mpq_lar_core_solver.m_column_types[j]) {
case column_type::free_column:
update_free_column_type_and_bound(j, kind, right_side, constr_index);
break;
case column_type::boxed:
update_boxed_column_type_and_bound(j, kind, right_side, constr_index);
break;
case column_type::low_bound:
update_low_bound_column_type_and_bound(j, kind, right_side, constr_index);
break;
case column_type::upper_bound:
update_upper_bound_column_type_and_bound(j, kind, right_side, constr_index);
break;
case column_type::fixed:
update_fixed_column_type_and_bound(j, kind, right_side, constr_index);
break;
default:
lean_assert(false); // cannot be here
}
}
void add_var_bound_on_constraint_for_term(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
lean_assert(is_term(j));
unsigned adjusted_term_index = adjust_term_index(j);
unsigned term_j;
if (try_get_val(m_ext_vars_to_columns, j, term_j)) {
mpq rs = right_side - m_orig_terms[adjusted_term_index]->m_v;
m_constraints.push_back(new lar_term_constraint(m_orig_terms[adjusted_term_index], kind, right_side));
update_column_type_and_bound(term_j, kind, rs, ci);
}
else {
add_constraint_from_term_and_create_new_column_row(j, m_orig_terms[adjusted_term_index], kind, right_side);
}
}
void add_constraint_from_term_and_create_new_column_row(unsigned term_j, const lar_term* term,
lconstraint_kind kind, const mpq & right_side) {
add_row_from_term_no_constraint(term, term_j);
unsigned j = A_r().column_count() - 1;
update_column_type_and_bound(j, kind, right_side - term->m_v, m_constraints.size());
m_constraints.push_back(new lar_term_constraint(term, kind, right_side));
lean_assert(A_r().column_count() == m_mpq_lar_core_solver.m_r_solver.m_costs.size());
}
void decide_on_strategy_and_adjust_initial_state() {
lean_assert(strategy_is_undecided());
if (m_vars_to_ul_pairs.size() > m_settings.column_number_threshold_for_using_lu_in_lar_solver) {
m_settings.simplex_strategy() = simplex_strategy_enum::lu;
} else {
m_settings.simplex_strategy() = simplex_strategy_enum::tableau_rows; // todo: when to switch to tableau_costs?
}
adjust_initial_state();
}
void adjust_initial_state() {
switch (m_settings.simplex_strategy()) {
case simplex_strategy_enum::lu:
adjust_initial_state_for_lu();
break;
case simplex_strategy_enum::tableau_rows:
adjust_initial_state_for_tableau_rows();
break;
case simplex_strategy_enum::tableau_costs:
lean_assert(false); // not implemented
case simplex_strategy_enum::undecided:
adjust_initial_state_for_tableau_rows();
break;
}
}
void adjust_initial_state_for_lu() {
copy_from_mpq_matrix(A_d());
unsigned n = A_d().column_count();
m_mpq_lar_core_solver.m_d_x.resize(n);
m_mpq_lar_core_solver.m_d_low_bounds.resize(n);
m_mpq_lar_core_solver.m_d_upper_bounds.resize(n);
m_mpq_lar_core_solver.m_d_heading = m_mpq_lar_core_solver.m_r_heading;
m_mpq_lar_core_solver.m_d_basis = m_mpq_lar_core_solver.m_r_basis;
/*
unsigned j = A_d().column_count();
A_d().add_column();
lean_assert(m_mpq_lar_core_solver.m_d_x.size() == j);
// lean_assert(m_mpq_lar_core_solver.m_d_low_bounds.size() == j && m_mpq_lar_core_solver.m_d_upper_bounds.size() == j); // restore later
m_mpq_lar_core_solver.m_d_x.resize(j + 1 );
m_mpq_lar_core_solver.m_d_low_bounds.resize(j + 1);
m_mpq_lar_core_solver.m_d_upper_bounds.resize(j + 1);
lean_assert(m_mpq_lar_core_solver.m_d_heading.size() == j); // as A().column_count() on the entry to the method
if (register_in_basis) {
A_d().add_row();
m_mpq_lar_core_solver.m_d_heading.push_back(m_mpq_lar_core_solver.m_d_basis.size());
m_mpq_lar_core_solver.m_d_basis.push_back(j);
}else {
m_mpq_lar_core_solver.m_d_heading.push_back(- static_cast<int>(m_mpq_lar_core_solver.m_d_nbasis.size()) - 1);
m_mpq_lar_core_solver.m_d_nbasis.push_back(j);
}*/
}
void adjust_initial_state_for_tableau_rows() {
for (unsigned j = 0; j < m_terms.size(); j++) {
if (contains(m_ext_vars_to_columns, j + m_terms_start_index))
continue;
add_row_from_term_no_constraint(m_terms[j], j + m_terms_start_index);
}
}
// this fills the last row of A_d and sets the basis column: -1 in the last column of the row
void fill_last_row_of_A_d(static_matrix<double, double> & A, const lar_term* ls) {
lean_assert(A.row_count() > 0);
lean_assert(A.column_count() > 0);
unsigned last_row = A.row_count() - 1;
lean_assert(A.m_rows[last_row].empty());
for (auto & t : ls->m_coeffs) {
lean_assert(!is_zero(t.second));
var_index j = t.first;
A.set(last_row, j, - t.second.get_double());
}
unsigned basis_j = A.column_count() - 1;
A.set(last_row, basis_j, - 1 );
}
void update_free_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_ind) {
mpq y_of_bound(0);
switch (kind) {
case LT:
y_of_bound = -1;
case LE:
m_mpq_lar_core_solver.m_column_types[j] = column_type::upper_bound;
lean_assert(m_mpq_lar_core_solver.m_column_types()[j] == column_type::upper_bound);
lean_assert(m_mpq_lar_core_solver.m_r_upper_bounds.size() > j);
{
auto up = numeric_pair<mpq>(right_side, y_of_bound);
m_mpq_lar_core_solver.m_r_upper_bounds[j] = up;
}
set_upper_bound_witness(j, constr_ind);
break;
case GT:
y_of_bound = 1;
case GE:
m_mpq_lar_core_solver.m_column_types[j] = column_type::low_bound;
lean_assert(m_mpq_lar_core_solver.m_r_upper_bounds.size() > j);
{
auto low = numeric_pair<mpq>(right_side, y_of_bound);
m_mpq_lar_core_solver.m_r_low_bounds[j] = low;
}
set_low_bound_witness(j, constr_ind);
break;
case EQ:
m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
m_mpq_lar_core_solver.m_r_low_bounds[j] = m_mpq_lar_core_solver.m_r_upper_bounds[j] = numeric_pair<mpq>(right_side, zero_of_type<mpq>());
set_upper_bound_witness(j, constr_ind);
set_low_bound_witness(j, constr_ind);
break;
default:
lean_unreachable();
}
m_columns_with_changed_bound.insert(j);
}
void update_upper_bound_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
lean_assert(m_mpq_lar_core_solver.m_column_types()[j] == column_type::upper_bound);
mpq y_of_bound(0);
switch (kind) {
case LT:
y_of_bound = -1;
case LE:
{
auto up = numeric_pair<mpq>(right_side, y_of_bound);
if (up < m_mpq_lar_core_solver.m_r_upper_bounds()[j]) {
m_mpq_lar_core_solver.m_r_upper_bounds[j] = up;
set_upper_bound_witness(j, ci);
m_columns_with_changed_bound.insert(j);
}
}
break;
case GT:
y_of_bound = 1;
case GE:
m_mpq_lar_core_solver.m_column_types[j] = column_type::boxed;
{
auto low = numeric_pair<mpq>(right_side, y_of_bound);
m_mpq_lar_core_solver.m_r_low_bounds[j] = low;
set_low_bound_witness(j, ci);
m_columns_with_changed_bound.insert(j);
if (low > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
m_status = INFEASIBLE;
m_infeasible_column_index = j;
} else {
m_mpq_lar_core_solver.m_column_types[j] = m_mpq_lar_core_solver.m_r_low_bounds()[j] < m_mpq_lar_core_solver.m_r_upper_bounds()[j]? column_type::boxed : column_type::fixed;
}
}
break;
case EQ:
{
auto v = numeric_pair<mpq>(right_side, zero_of_type<mpq>());
if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
m_status = INFEASIBLE;
set_low_bound_witness(j, ci);
m_infeasible_column_index = j;
} else {
m_mpq_lar_core_solver.m_r_low_bounds[j] = m_mpq_lar_core_solver.m_r_upper_bounds[j] = v;
m_columns_with_changed_bound.insert(j);
set_low_bound_witness(j, ci);
set_upper_bound_witness(j, ci);
m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
}
break;
}
break;
default:
lean_unreachable();
}
}
void update_boxed_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
lean_assert(m_status == INFEASIBLE || (m_mpq_lar_core_solver.m_column_types()[j] == column_type::boxed && m_mpq_lar_core_solver.m_r_low_bounds()[j] < m_mpq_lar_core_solver.m_r_upper_bounds()[j]));
mpq y_of_bound(0);
switch (kind) {
case LT:
y_of_bound = -1;
case LE:
{
auto up = numeric_pair<mpq>(right_side, y_of_bound);
if (up < m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
m_mpq_lar_core_solver.m_r_upper_bounds[j] = up;
set_upper_bound_witness(j, ci);
m_columns_with_changed_bound.insert(j);
}
if (up < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
m_status = INFEASIBLE;
lean_assert(false);
m_infeasible_column_index = j;
} else {
if (m_mpq_lar_core_solver.m_r_low_bounds()[j] == m_mpq_lar_core_solver.m_r_upper_bounds()[j])
m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
}
}
break;
case GT:
y_of_bound = 1;
case GE:
{
auto low = numeric_pair<mpq>(right_side, y_of_bound);
if (low > m_mpq_lar_core_solver.m_r_low_bounds[j]) {
m_mpq_lar_core_solver.m_r_low_bounds[j] = low;
m_columns_with_changed_bound.insert(j);
set_low_bound_witness(j, ci);
}
if (low > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
m_status = INFEASIBLE;
m_infeasible_column_index = j;
} else if ( low == m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
}
}
break;
case EQ:
{
auto v = numeric_pair<mpq>(right_side, zero_of_type<mpq>());
if (v < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
m_status = INFEASIBLE;
m_infeasible_column_index = j;
set_upper_bound_witness(j, ci);
} else if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
m_status = INFEASIBLE;
m_infeasible_column_index = j;
set_low_bound_witness(j, ci);
} else {
m_mpq_lar_core_solver.m_r_low_bounds[j] = m_mpq_lar_core_solver.m_r_upper_bounds[j] = v;
set_low_bound_witness(j, ci);
set_upper_bound_witness(j, ci);
m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
m_columns_with_changed_bound.insert(j);
}
break;
}
default:
lean_unreachable();
}
}
void update_low_bound_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
lean_assert(m_mpq_lar_core_solver.m_column_types()[j] == column_type::low_bound);
mpq y_of_bound(0);
switch (kind) {
case LT:
y_of_bound = -1;
case LE:
{
auto up = numeric_pair<mpq>(right_side, y_of_bound);
m_mpq_lar_core_solver.m_r_upper_bounds[j] = up;
set_upper_bound_witness(j, ci);
m_columns_with_changed_bound.insert(j);
if (up < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
m_status = INFEASIBLE;
m_infeasible_column_index = j;
} else {
m_mpq_lar_core_solver.m_column_types[j] = m_mpq_lar_core_solver.m_r_low_bounds()[j] < m_mpq_lar_core_solver.m_r_upper_bounds()[j]? column_type::boxed : column_type::fixed;
}
}
break;
case GT:
y_of_bound = 1;
case GE:
{
auto low = numeric_pair<mpq>(right_side, y_of_bound);
if (low > m_mpq_lar_core_solver.m_r_low_bounds[j]) {
m_mpq_lar_core_solver.m_r_low_bounds[j] = low;
m_columns_with_changed_bound.insert(j);
set_low_bound_witness(j, ci);
}
}
break;
case EQ:
{
auto v = numeric_pair<mpq>(right_side, zero_of_type<mpq>());
if (v < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
m_status = INFEASIBLE;
m_infeasible_column_index = j;
set_upper_bound_witness(j, ci);
} else {
m_mpq_lar_core_solver.m_r_low_bounds[j] = m_mpq_lar_core_solver.m_r_upper_bounds[j] = v;
set_low_bound_witness(j, ci);
set_upper_bound_witness(j, ci);
m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
}
m_columns_with_changed_bound.insert(j);
break;
}
default:
lean_unreachable();
}
}
void update_fixed_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
lean_assert(m_status == INFEASIBLE || (m_mpq_lar_core_solver.m_column_types()[j] == column_type::fixed && m_mpq_lar_core_solver.m_r_low_bounds()[j] == m_mpq_lar_core_solver.m_r_upper_bounds()[j]));
lean_assert(m_status == INFEASIBLE || (m_mpq_lar_core_solver.m_r_low_bounds()[j].y.is_zero() && m_mpq_lar_core_solver.m_r_upper_bounds()[j].y.is_zero()));
auto v = numeric_pair<mpq>(right_side, mpq(0));
mpq y_of_bound(0);
switch (kind) {
case LT:
if (v <= m_mpq_lar_core_solver.m_r_low_bounds[j]) {
m_status = INFEASIBLE;
m_infeasible_column_index = j;
set_upper_bound_witness(j, ci);
}
break;
case LE:
{
if (v < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
m_status = INFEASIBLE;
m_infeasible_column_index = j;
set_upper_bound_witness(j, ci);
}
}
break;
case GT:
{
if (v >= m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
m_status = INFEASIBLE;
m_infeasible_column_index =j;
set_low_bound_witness(j, ci);
}
}
break;
case GE:
{
if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
m_status = INFEASIBLE;
m_infeasible_column_index = j;
set_low_bound_witness(j, ci);
}
}
break;
case EQ:
{
if (v < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
m_status = INFEASIBLE;
m_infeasible_column_index = j;
set_upper_bound_witness(j, ci);
} else if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
m_status = INFEASIBLE;
m_infeasible_column_index = j;
set_low_bound_witness(j, ci);
}
break;
}
default:
lean_unreachable();
}
}

66
src/util/lp/int_set.h Normal file
View file

@ -0,0 +1,66 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#pragma once
#include "util/vector.h"
#include "util/lp/indexed_vector.h"
#include <ostream>
namespace lean {
// serves at a set of non-negative integers smaller than the set size
class int_set {
vector<int> m_data;
public:
vector<int> m_index;
int_set(unsigned size): m_data(size, -1) {}
int_set() {}
bool contains(unsigned j) const {
if (j >= m_data.size())
return false;
return m_data[j] >= 0;
}
void insert(unsigned j) {
lean_assert(j < m_data.size());
if (contains(j)) return;
m_data[j] = m_index.size();
m_index.push_back(j);
}
void erase(unsigned j) {
if (!contains(j)) return;
unsigned pos_j = m_data[j];
unsigned last_pos = m_index.size() - 1;
int last_j = m_index[last_pos];
if (last_pos != pos_j) {
// move last to j spot
m_data[last_j] = pos_j;
m_index[pos_j] = last_j;
}
m_index.pop_back();
m_data[j] = -1;
}
void resize(unsigned size) {
m_data.resize(size, -1);
}
void increase_size_by_one() {
resize(m_data.size() + 1);
}
unsigned data_size() const { return m_data.size(); }
unsigned size() const { return m_index.size();}
bool is_empty() const { return size() == 0; }
void clear() {
for (unsigned j : m_index)
m_data[j] = -1;
m_index.resize(0);
}
void print(std::ostream & out ) const {
for (unsigned j : m_index) {
out << j << " ";
}
out << std::endl;
}
};
}

50
src/util/lp/iterator_on_column.h Normal file
View file

@ -0,0 +1,50 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#pragma once
#include "util/lp/linear_combination_iterator.h"
#include "util/lp/static_matrix.h"
#include "util/lp/lar_term.h"
namespace lean {
template <typename T, typename X>
struct iterator_on_column:linear_combination_iterator<T> {
const vector<column_cell>& m_column; // the offset in term coeffs
const static_matrix<T, X> & m_A;
int m_i; // the initial offset in the column
unsigned size() const { return m_column.size(); }
iterator_on_column(const vector<column_cell>& column, const static_matrix<T,X> & A) // the offset in term coeffs
:
m_column(column),
m_A(A),
m_i(-1) {}
bool next(mpq & a, unsigned & i) {
if (++m_i >= static_cast<int>(m_column.size()))
return false;
const column_cell& c = m_column[m_i];
a = m_A.get_val(c);
i = c.m_i;
return true;
}
bool next(unsigned & i) {
if (++m_i >= static_cast<int>(m_column.size()))
return false;
const column_cell& c = m_column[m_i];
i = c.m_i;
return true;
}
void reset() {
m_i = -1;
}
linear_combination_iterator<mpq> * clone() {
iterator_on_column * r = new iterator_on_column(m_column, m_A);
return r;
}
};
}

38
src/util/lp/iterator_on_indexed_vector.h Normal file
View file

@ -0,0 +1,38 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#pragma once
#include "util/lp/linear_combination_iterator.h"
namespace lean {
template <typename T>
struct iterator_on_indexed_vector:linear_combination_iterator<T> {
const indexed_vector<T> & m_v;
unsigned m_offset;
iterator_on_indexed_vector(const indexed_vector<T> & v) :
m_v(v),
m_offset(0)
{}
unsigned size() const { return m_v.m_index.size(); }
bool next(T & a, unsigned & i) {
if (m_offset >= m_v.m_index.size())
return false;
i = m_v.m_index[m_offset++];
a = m_v.m_data[i];
return true;
}
bool next(unsigned & i) {
if (m_offset >= m_v.m_index.size())
return false;
i = m_v.m_index[m_offset++];
return true;
}
void reset() {
m_offset = 0;
}
linear_combination_iterator<T>* clone() {
return new iterator_on_indexed_vector(m_v);
}
};
}

44
src/util/lp/iterator_on_pivot_row.h Normal file
View file

@ -0,0 +1,44 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#pragma once
#include "util/lp/iterator_on_indexed_vector.h"
namespace lean {
template <typename T>
struct iterator_on_pivot_row:linear_combination_iterator<T> {
bool m_basis_returned;
const indexed_vector<T> & m_v;
unsigned m_basis_j;
iterator_on_indexed_vector<T> m_it;
unsigned size() const { return m_it.size(); }
iterator_on_pivot_row(const indexed_vector<T> & v, unsigned basis_j) :
m_basis_returned(false),
m_v(v), m_basis_j(basis_j), m_it(v) {}
bool next(T & a, unsigned & i) {
if (m_basis_returned == false) {
m_basis_returned = true;
a = one_of_type<T>();
i = m_basis_j;
return true;
}
return m_it.next(a, i);
}
bool next(unsigned & i) {
if (m_basis_returned == false) {
m_basis_returned = true;
i = m_basis_j;
return true;
}
return m_it.next(i);
}
void reset() {
m_basis_returned = false;
m_it.reset();
}
linear_combination_iterator<T> * clone() {
iterator_on_pivot_row * r = new iterator_on_pivot_row(m_v, m_basis_j);
return r;
}
};
}

37
src/util/lp/iterator_on_row.h Normal file
View file

@ -0,0 +1,37 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#pragma once
#include "util/lp/linear_combination_iterator.h"
namespace lean {
template <typename T>
struct iterator_on_row:linear_combination_iterator<T> {
const vector<row_cell<T>> & m_row;
unsigned m_i; // offset
iterator_on_row(const vector<row_cell<T>> & row) : m_row(row), m_i(0)
{}
unsigned size() const { return m_row.size(); }
bool next(T & a, unsigned & i) {
if (m_i == m_row.size())
return false;
auto &c = m_row[m_i++];
i = c.m_j;
a = c.get_val();
return true;
}
bool next(unsigned & i) {
if (m_i == m_row.size())
return false;
auto &c = m_row[m_i++];
i = c.m_j;
return true;
}
void reset() {
m_i = 0;
}
linear_combination_iterator<T>* clone() {
return new iterator_on_row(m_row);
}
};
}

57
src/util/lp/iterator_on_term_with_basis_var.h Normal file
View file

@ -0,0 +1,57 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#pragma once
#include "util/lp/linear_combination_iterator.h"
#include "util/lp/numeric_pair.h"
#include "util/lp/lar_term.h"
namespace lean {
struct iterator_on_term_with_basis_var:linear_combination_iterator<mpq> {
const lar_term & m_term;
std::unordered_map<unsigned, mpq>::const_iterator m_i; // the offset in term coeffs
bool m_term_j_returned;
unsigned m_term_j;
unsigned size() const {return static_cast<unsigned>(m_term.m_coeffs.size() + 1);}
iterator_on_term_with_basis_var(const lar_term & t, unsigned term_j) :
m_term(t),
m_i(t.m_coeffs.begin()),
m_term_j_returned(false),
m_term_j(term_j) {}
bool next(mpq & a, unsigned & i) {
if (m_term_j_returned == false) {
m_term_j_returned = true;
a = - one_of_type<mpq>();
i = m_term_j;
return true;
}
if (m_i == m_term.m_coeffs.end())
return false;
i = m_i->first;
a = m_i->second;
m_i++;
return true;
}
bool next(unsigned & i) {
if (m_term_j_returned == false) {
m_term_j_returned = true;
i = m_term_j;
return true;
}
if (m_i == m_term.m_coeffs.end())
return false;
i = m_i->first;
m_i++;
return true;
}
void reset() {
m_term_j_returned = false;
m_i = m_term.m_coeffs.begin();
}
linear_combination_iterator<mpq> * clone() {
iterator_on_term_with_basis_var * r = new iterator_on_term_with_basis_var(m_term, m_term_j);
return r;
}
};
}

86
src/util/lp/lar_constraints.h Normal file
View file

@ -0,0 +1,86 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#pragma once
#include "util/vector.h"
#include <utility>
#include <unordered_map>
#include <string>
#include <algorithm>
#include "util/lp/lp_utils.h"
#include "util/lp/ul_pair.h"
#include "util/lp/lar_term.h"
namespace lean {
inline lconstraint_kind flip_kind(lconstraint_kind t) {
return static_cast<lconstraint_kind>( - static_cast<int>(t));
}
inline std::string lconstraint_kind_string(lconstraint_kind t) {
switch (t) {
case LE: return std::string("<=");
case LT: return std::string("<");
case GE: return std::string(">=");
case GT: return std::string(">");
case EQ: return std::string("=");
}
lean_unreachable();
return std::string(); // it is unreachable
}
class lar_base_constraint {
public:
lconstraint_kind m_kind;
mpq m_right_side;
virtual vector<std::pair<mpq, var_index>> get_left_side_coefficients() const = 0;
lar_base_constraint() {}
lar_base_constraint(lconstraint_kind kind, const mpq& right_side) :m_kind(kind), m_right_side(right_side) {}
virtual unsigned size() const = 0;
virtual ~lar_base_constraint(){}
virtual mpq get_free_coeff_of_left_side() const { return zero_of_type<mpq>();}
};
struct lar_var_constraint: public lar_base_constraint {
unsigned m_j;
vector<std::pair<mpq, var_index>> get_left_side_coefficients() const {
vector<std::pair<mpq, var_index>> ret;
ret.push_back(std::make_pair(one_of_type<mpq>(), m_j));
return ret;
}
unsigned size() const { return 1;}
lar_var_constraint(unsigned j, lconstraint_kind kind, const mpq& right_side) : lar_base_constraint(kind, right_side), m_j(j) { }
};
struct lar_term_constraint: public lar_base_constraint {
const lar_term * m_term;
vector<std::pair<mpq, var_index>> get_left_side_coefficients() const {
return m_term->coeffs_as_vector();
}
unsigned size() const { return m_term->size();}
lar_term_constraint(const lar_term *t, lconstraint_kind kind, const mpq& right_side) : lar_base_constraint(kind, right_side), m_term(t) { }
virtual mpq get_free_coeff_of_left_side() const { return m_term->m_v;}
};
class lar_constraint : public lar_base_constraint {
public:
vector<std::pair<mpq, var_index>> m_coeffs;
lar_constraint() {}
lar_constraint(const vector<std::pair<mpq, var_index>> & left_side, lconstraint_kind kind, const mpq & right_side)
: lar_base_constraint(kind, right_side), m_coeffs(left_side) {}
lar_constraint(const lar_base_constraint & c) {
lean_assert(false); // should not be called : todo!
}
unsigned size() const {
return static_cast<unsigned>(m_coeffs.size());
}
vector<std::pair<mpq, var_index>> get_left_side_coefficients() const { return m_coeffs; }
};
}

801
src/util/lp/lar_core_solver.h Normal file
View file

@ -0,0 +1,801 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#pragma once
#include "util/vector.h"
#include <string>
#include <utility>
#include "util/lp/lp_core_solver_base.h"
#include <algorithm>
#include "util/lp/indexed_vector.h"
#include "util/lp/binary_heap_priority_queue.h"
#include "util/lp/breakpoint.h"
#include "util/lp/stacked_unordered_set.h"
#include "util/lp/lp_primal_core_solver.h"
#include "util/lp/stacked_vector.h"
#include "util/lp/lar_solution_signature.h"
#include "util/lp/iterator_on_column.h"
#include "util/lp/iterator_on_indexed_vector.h"
#include "util/lp/stacked_value.h"
namespace lean {
class lar_core_solver {
// m_sign_of_entering is set to 1 if the entering variable needs
// to grow and is set to -1 otherwise
int m_sign_of_entering_delta;
vector<std::pair<mpq, unsigned>> m_infeasible_linear_combination;
int m_infeasible_sum_sign; // todo: get rid of this field
vector<numeric_pair<mpq>> m_right_sides_dummy;
vector<mpq> m_costs_dummy;
vector<double> m_d_right_sides_dummy;
vector<double> m_d_costs_dummy;
public:
stacked_value<simplex_strategy_enum> m_stacked_simplex_strategy;
stacked_vector<column_type> m_column_types;
// r - solver fields, for rational numbers
vector<numeric_pair<mpq>> m_r_x; // the solution
stacked_vector<numeric_pair<mpq>> m_r_low_bounds;
stacked_vector<numeric_pair<mpq>> m_r_upper_bounds;
static_matrix<mpq, numeric_pair<mpq>> m_r_A;
stacked_vector<unsigned> m_r_pushed_basis;
vector<unsigned> m_r_basis;
vector<unsigned> m_r_nbasis;
vector<int> m_r_heading;
stacked_vector<unsigned> m_r_columns_nz;
stacked_vector<unsigned> m_r_rows_nz;
// d - solver fields, for doubles
vector<double> m_d_x; // the solution in doubles
vector<double> m_d_low_bounds;
vector<double> m_d_upper_bounds;
static_matrix<double, double> m_d_A;
stacked_vector<unsigned> m_d_pushed_basis;
vector<unsigned> m_d_basis;
vector<unsigned> m_d_nbasis;
vector<int> m_d_heading;
lp_primal_core_solver<mpq, numeric_pair<mpq>> m_r_solver; // solver in rational numbers
lp_primal_core_solver<double, double> m_d_solver; // solver in doubles
lar_core_solver(
lp_settings & settings,
const column_namer & column_names
);
lp_settings & settings() { return m_r_solver.m_settings;}
const lp_settings & settings() const { return m_r_solver.m_settings;}
int get_infeasible_sum_sign() const { return m_infeasible_sum_sign; }
const vector<std::pair<mpq, unsigned>> & get_infeasibility_info(int & inf_sign) const {
inf_sign = m_infeasible_sum_sign;
return m_infeasible_linear_combination;
}
void fill_not_improvable_zero_sum_from_inf_row();
column_type get_column_type(unsigned j) { return m_column_types[j];}
void init_costs(bool first_time);
void init_cost_for_column(unsigned j);
// returns m_sign_of_alpha_r
int column_is_out_of_bounds(unsigned j);
void calculate_pivot_row(unsigned i);
void print_pivot_row(std::ostream & out, unsigned row_index) const { // remove later debug !!!!
for (unsigned j : m_r_solver.m_pivot_row.m_index) {
if (numeric_traits<mpq>::is_pos(m_r_solver.m_pivot_row.m_data[j]))
out << "+";
out << m_r_solver.m_pivot_row.m_data[j] << m_r_solver.column_name(j) << " ";
}
out << " +" << m_r_solver.column_name(m_r_solver.m_basis[row_index]) << std::endl;
for (unsigned j : m_r_solver.m_pivot_row.m_index) {
m_r_solver.print_column_bound_info(j, out);
}
m_r_solver.print_column_bound_info(m_r_solver.m_basis[row_index], out);
}
void advance_on_sorted_breakpoints(unsigned entering);
void change_slope_on_breakpoint(unsigned entering, breakpoint<numeric_pair<mpq>> * b, mpq & slope_at_entering);
bool row_is_infeasible(unsigned row);
bool row_is_evidence(unsigned row);
bool find_evidence_row();
void prefix_r();
void prefix_d();
unsigned m_m() const {
return m_r_A.row_count();
}
unsigned m_n() const {
return m_r_A.column_count();
}
bool is_tiny() const { return this->m_m() < 10 && this->m_n() < 20; }
bool is_empty() const { return this->m_m() == 0 && this->m_n() == 0; }
template <typename L>
int get_sign(const L & v) {
return v > zero_of_type<L>() ? 1 : (v < zero_of_type<L>() ? -1 : 0);
}
void fill_evidence(unsigned row);
void solve();
bool low_bounds_are_set() const { return true; }
const indexed_vector<mpq> & get_pivot_row() const {
return m_r_solver.m_pivot_row;
}
void fill_not_improvable_zero_sum();
void pop_basis(unsigned k) {
if (!settings().use_tableau()) {
m_r_pushed_basis.pop(k);
m_r_basis = m_r_pushed_basis();
m_r_solver.init_basis_heading_and_non_basic_columns_vector();
m_d_pushed_basis.pop(k);
m_d_basis = m_d_pushed_basis();
m_d_solver.init_basis_heading_and_non_basic_columns_vector();
} else {
m_d_basis = m_r_basis;
m_d_nbasis = m_r_nbasis;
m_d_heading = m_r_heading;
}
}
void push() {
lean_assert(m_r_solver.basis_heading_is_correct());
lean_assert(!need_to_presolve_with_double_solver() || m_d_solver.basis_heading_is_correct());
lean_assert(m_column_types.size() == m_r_A.column_count());
m_stacked_simplex_strategy = settings().simplex_strategy();
m_stacked_simplex_strategy.push();
m_column_types.push();
// rational
if (!settings().use_tableau())
m_r_A.push();
m_r_low_bounds.push();
m_r_upper_bounds.push();
if (!settings().use_tableau()) {
push_vector(m_r_pushed_basis, m_r_basis);
push_vector(m_r_columns_nz, m_r_solver.m_columns_nz);
push_vector(m_r_rows_nz, m_r_solver.m_rows_nz);
}
m_d_A.push();
if (!settings().use_tableau())
push_vector(m_d_pushed_basis, m_d_basis);
}
template <typename K>
void push_vector(stacked_vector<K> & pushed_vector, const vector<K> & vector) {
lean_assert(pushed_vector.size() <= vector.size());
for (unsigned i = 0; i < vector.size();i++) {
if (i == pushed_vector.size()) {
pushed_vector.push_back(vector[i]);
} else {
pushed_vector[i] = vector[i];
}
}
pushed_vector.push();
}
void pop_markowitz_counts(unsigned k) {
m_r_columns_nz.pop(k);
m_r_rows_nz.pop(k);
m_r_solver.m_columns_nz.resize(m_r_columns_nz.size());
m_r_solver.m_rows_nz.resize(m_r_rows_nz.size());
for (unsigned i = 0; i < m_r_columns_nz.size(); i++)
m_r_solver.m_columns_nz[i] = m_r_columns_nz[i];
for (unsigned i = 0; i < m_r_rows_nz.size(); i++)
m_r_solver.m_rows_nz[i] = m_r_rows_nz[i];
}
void pop(unsigned k) {
// rationals
if (!settings().use_tableau())
m_r_A.pop(k);
m_r_low_bounds.pop(k);
m_r_upper_bounds.pop(k);
m_column_types.pop(k);
if (m_r_solver.m_factorization != nullptr) {
delete m_r_solver.m_factorization;
m_r_solver.m_factorization = nullptr;
}
m_r_x.resize(m_r_A.column_count());
m_r_solver.m_costs.resize(m_r_A.column_count());
m_r_solver.m_d.resize(m_r_A.column_count());
if(!settings().use_tableau())
pop_markowitz_counts(k);
m_d_A.pop(k);
if (m_d_solver.m_factorization != nullptr) {
delete m_d_solver.m_factorization;
m_d_solver.m_factorization = nullptr;
}
m_d_x.resize(m_d_A.column_count());
pop_basis(k);
m_stacked_simplex_strategy.pop(k);
settings().simplex_strategy() = m_stacked_simplex_strategy;
lean_assert(m_r_solver.basis_heading_is_correct());
lean_assert(!need_to_presolve_with_double_solver() || m_d_solver.basis_heading_is_correct());
}
bool need_to_presolve_with_double_solver() const {
return settings().simplex_strategy() == simplex_strategy_enum::lu;
}
template <typename L>
bool is_zero_vector(const vector<L> & b) {
for (const L & m: b)
if (!is_zero(m)) return false;
return true;
}
bool update_xj_and_get_delta(unsigned j, non_basic_column_value_position pos_type, numeric_pair<mpq> & delta) {
auto & x = m_r_x[j];
switch (pos_type) {
case at_low_bound:
if (x == m_r_solver.m_low_bounds[j])
return false;
delta = m_r_solver.m_low_bounds[j] - x;
m_r_solver.m_x[j] = m_r_solver.m_low_bounds[j];
break;
case at_fixed:
case at_upper_bound:
if (x == m_r_solver.m_upper_bounds[j])
return false;
delta = m_r_solver.m_upper_bounds[j] - x;
x = m_r_solver.m_upper_bounds[j];
break;
case free_of_bounds: {
return false;
}
case not_at_bound:
switch (m_column_types[j]) {
case column_type::free_column:
return false;
case column_type::upper_bound:
delta = m_r_solver.m_upper_bounds[j] - x;
x = m_r_solver.m_upper_bounds[j];
break;
case column_type::low_bound:
delta = m_r_solver.m_low_bounds[j] - x;
x = m_r_solver.m_low_bounds[j];
break;
case column_type::boxed:
if (x > m_r_solver.m_upper_bounds[j]) {
delta = m_r_solver.m_upper_bounds[j] - x;
x += m_r_solver.m_upper_bounds[j];
} else {
delta = m_r_solver.m_low_bounds[j] - x;
x = m_r_solver.m_low_bounds[j];
}
break;
case column_type::fixed:
delta = m_r_solver.m_low_bounds[j] - x;
x = m_r_solver.m_low_bounds[j];
break;
default:
lean_assert(false);
}
break;
default:
lean_unreachable();
}
m_r_solver.remove_column_from_inf_set(j);
return true;
}
void prepare_solver_x_with_signature_tableau(const lar_solution_signature & signature) {
lean_assert(m_r_solver.inf_set_is_correct());
for (auto &t : signature) {
unsigned j = t.first;
if (m_r_heading[j] >= 0)
continue;
auto pos_type = t.second;
numeric_pair<mpq> delta;
if (!update_xj_and_get_delta(j, pos_type, delta))
continue;
for (const auto & cc : m_r_solver.m_A.m_columns[j]){
unsigned i = cc.m_i;
unsigned jb = m_r_solver.m_basis[i];
m_r_solver.m_x[jb] -= delta * m_r_solver.m_A.get_val(cc);
m_r_solver.update_column_in_inf_set(jb);
}
lean_assert(m_r_solver.A_mult_x_is_off() == false);
}
lean_assert(m_r_solver.inf_set_is_correct());
}
template <typename L, typename K>
void prepare_solver_x_with_signature(const lar_solution_signature & signature, lp_primal_core_solver<L,K> & s) {
for (auto &t : signature) {
unsigned j = t.first;
lean_assert(m_r_heading[j] < 0);
auto pos_type = t.second;
switch (pos_type) {
case at_low_bound:
s.m_x[j] = s.m_low_bounds[j];
break;
case at_fixed:
case at_upper_bound:
s.m_x[j] = s.m_upper_bounds[j];
break;
case free_of_bounds: {
s.m_x[j] = zero_of_type<K>();
continue;
}
case not_at_bound:
switch (m_column_types[j]) {
case column_type::free_column:
lean_assert(false); // unreachable
case column_type::upper_bound:
s.m_x[j] = s.m_upper_bounds[j];
break;
case column_type::low_bound:
s.m_x[j] = s.m_low_bounds[j];
break;
case column_type::boxed:
if (settings().random_next() % 2) {
s.m_x[j] = s.m_low_bounds[j];
} else {
s.m_x[j] = s.m_upper_bounds[j];
}
break;
case column_type::fixed:
s.m_x[j] = s.m_low_bounds[j];
break;
default:
lean_assert(false);
}
break;
default:
lean_unreachable();
}
}
lean_assert(is_zero_vector(s.m_b));
s.solve_Ax_eq_b();
}
template <typename L, typename K>
void catch_up_in_lu_in_reverse(const vector<unsigned> & trace_of_basis_change, lp_primal_core_solver<L,K> & cs) {
// recover the previous working basis
for (unsigned i = trace_of_basis_change.size(); i > 0; i-= 2) {
unsigned entering = trace_of_basis_change[i-1];
unsigned leaving = trace_of_basis_change[i-2];
cs.change_basis_unconditionally(entering, leaving);
}
cs.init_lu();
}
//basis_heading is the basis heading of the solver owning trace_of_basis_change
// here we compact the trace as we go to avoid unnecessary column changes
template <typename L, typename K>
void catch_up_in_lu(const vector<unsigned> & trace_of_basis_change, const vector<int> & basis_heading, lp_primal_core_solver<L,K> & cs) {
if (cs.m_factorization == nullptr || cs.m_factorization->m_refactor_counter + trace_of_basis_change.size()/2 >= 200) {
for (unsigned i = 0; i < trace_of_basis_change.size(); i+= 2) {
unsigned entering = trace_of_basis_change[i];
unsigned leaving = trace_of_basis_change[i+1];
cs.change_basis_unconditionally(entering, leaving);
}
if (cs.m_factorization != nullptr)
delete cs.m_factorization;
cs.m_factorization = nullptr;
} else {
indexed_vector<L> w(cs.m_A.row_count());
// the queues of delayed indices
std::queue<unsigned> entr_q, leav_q;
auto * l = cs.m_factorization;
lean_assert(l->get_status() == LU_status::OK);
for (unsigned i = 0; i < trace_of_basis_change.size(); i+= 2) {
unsigned entering = trace_of_basis_change[i];
unsigned leaving = trace_of_basis_change[i+1];
bool good_e = basis_heading[entering] >= 0 && cs.m_basis_heading[entering] < 0;
bool good_l = basis_heading[leaving] < 0 && cs.m_basis_heading[leaving] >= 0;
if (!good_e && !good_l) continue;
if (good_e && !good_l) {
while (!leav_q.empty() && cs.m_basis_heading[leav_q.front()] < 0)
leav_q.pop();
if (!leav_q.empty()) {
leaving = leav_q.front();
leav_q.pop();
} else {
entr_q.push(entering);
continue;
}
} else if (!good_e && good_l) {
while (!entr_q.empty() && cs.m_basis_heading[entr_q.front()] >= 0)
entr_q.pop();
if (!entr_q.empty()) {
entering = entr_q.front();
entr_q.pop();
} else {
leav_q.push(leaving);
continue;
}
}
lean_assert(cs.m_basis_heading[entering] < 0);
lean_assert(cs.m_basis_heading[leaving] >= 0);
if (l->get_status() == LU_status::OK) {
l->prepare_entering(entering, w); // to init vector w
l->replace_column(zero_of_type<L>(), w, cs.m_basis_heading[leaving]);
}
cs.change_basis_unconditionally(entering, leaving);
}
if (l->get_status() != LU_status::OK) {
delete l;
cs.m_factorization = nullptr;
}
}
if (cs.m_factorization == nullptr) {
if (numeric_traits<L>::precise())
init_factorization(cs.m_factorization, cs.m_A, cs.m_basis, settings());
}
}
bool no_r_lu() const {
return m_r_solver.m_factorization == nullptr || m_r_solver.m_factorization->get_status() == LU_status::Degenerated;
}
void solve_on_signature_tableau(const lar_solution_signature & signature, const vector<unsigned> & changes_of_basis) {
r_basis_is_OK();
lean_assert(settings().use_tableau());
bool r = catch_up_in_lu_tableau(changes_of_basis, m_d_solver.m_basis_heading);
if (!r) { // it is the case where m_d_solver gives a degenerated basis
prepare_solver_x_with_signature_tableau(signature); // still are going to use the signature partially
m_r_solver.find_feasible_solution();
m_d_basis = m_r_basis;
m_d_heading = m_r_heading;
m_d_nbasis = m_r_nbasis;
delete m_d_solver.m_factorization;
m_d_solver.m_factorization = nullptr;
} else {
prepare_solver_x_with_signature_tableau(signature);
m_r_solver.start_tracing_basis_changes();
m_r_solver.find_feasible_solution();
if (settings().get_cancel_flag())
return;
m_r_solver.stop_tracing_basis_changes();
// and now catch up in the double solver
lean_assert(m_r_solver.total_iterations() >= m_r_solver.m_trace_of_basis_change_vector.size() /2);
catch_up_in_lu(m_r_solver.m_trace_of_basis_change_vector, m_r_solver.m_basis_heading, m_d_solver);
}
lean_assert(r_basis_is_OK());
}
bool adjust_x_of_column(unsigned j) {
/*
if (m_r_solver.m_basis_heading[j] >= 0) {
return false;
}
if (m_r_solver.column_is_feasible(j)) {
return false;
}
m_r_solver.snap_column_to_bound_tableau(j);
lean_assert(m_r_solver.column_is_feasible(j));
m_r_solver.m_inf_set.erase(j);
*/
lean_assert(false);
return true;
}
bool catch_up_in_lu_tableau(const vector<unsigned> & trace_of_basis_change, const vector<int> & basis_heading) {
lean_assert(r_basis_is_OK());
// the queues of delayed indices
std::queue<unsigned> entr_q, leav_q;
for (unsigned i = 0; i < trace_of_basis_change.size(); i+= 2) {
unsigned entering = trace_of_basis_change[i];
unsigned leaving = trace_of_basis_change[i+1];
bool good_e = basis_heading[entering] >= 0 && m_r_solver.m_basis_heading[entering] < 0;
bool good_l = basis_heading[leaving] < 0 && m_r_solver.m_basis_heading[leaving] >= 0;
if (!good_e && !good_l) continue;
if (good_e && !good_l) {
while (!leav_q.empty() && m_r_solver.m_basis_heading[leav_q.front()] < 0)
leav_q.pop();
if (!leav_q.empty()) {
leaving = leav_q.front();
leav_q.pop();
} else {
entr_q.push(entering);
continue;
}
} else if (!good_e && good_l) {
while (!entr_q.empty() && m_r_solver.m_basis_heading[entr_q.front()] >= 0)
entr_q.pop();
if (!entr_q.empty()) {
entering = entr_q.front();
entr_q.pop();
} else {
leav_q.push(leaving);
continue;
}
}
lean_assert(m_r_solver.m_basis_heading[entering] < 0);
lean_assert(m_r_solver.m_basis_heading[leaving] >= 0);
m_r_solver.change_basis_unconditionally(entering, leaving);
if(!m_r_solver.pivot_column_tableau(entering, m_r_solver.m_basis_heading[entering])) {
// unroll the last step
m_r_solver.change_basis_unconditionally(leaving, entering);
#ifdef LEAN_DEBUG
bool t =
#endif
m_r_solver.pivot_column_tableau(leaving, m_r_solver.m_basis_heading[leaving]);
#ifdef LEAN_DEBUG
lean_assert(t);
#endif
return false;
}
}
lean_assert(r_basis_is_OK());
return true;
}
bool r_basis_is_OK() const {
#ifdef LEAN_DEBUG
if (!m_r_solver.m_settings.use_tableau())
return true;
for (unsigned j : m_r_solver.m_basis) {
lean_assert(m_r_solver.m_A.m_columns[j].size() == 1);
lean_assert(m_r_solver.m_A.get_val(m_r_solver.m_A.m_columns[j][0]) == one_of_type<mpq>());
}
for (unsigned j =0; j < m_r_solver.m_basis_heading.size(); j++) {
if (m_r_solver.m_basis_heading[j] >= 0) continue;
if (m_r_solver.m_column_types[j] == column_type::fixed) continue;
lean_assert(static_cast<unsigned>(- m_r_solver.m_basis_heading[j] - 1) < m_r_solver.m_column_types.size());
lean_assert( m_r_solver.m_basis_heading[j] <= -1);
}
#endif
return true;
}
void solve_on_signature(const lar_solution_signature & signature, const vector<unsigned> & changes_of_basis) {
lean_assert(!settings().use_tableau());
if (m_r_solver.m_factorization == nullptr) {
for (unsigned j = 0; j < changes_of_basis.size(); j+=2) {
unsigned entering = changes_of_basis[j];
unsigned leaving = changes_of_basis[j + 1];
m_r_solver.change_basis_unconditionally(entering, leaving);
}
init_factorization(m_r_solver.m_factorization, m_r_A, m_r_basis, settings());
} else {
catch_up_in_lu(changes_of_basis, m_d_solver.m_basis_heading, m_r_solver);
}
if (no_r_lu()) { // it is the case where m_d_solver gives a degenerated basis, we need to roll back
// std::cout << "no_r_lu" << std::endl;
catch_up_in_lu_in_reverse(changes_of_basis, m_r_solver);
m_r_solver.find_feasible_solution();
m_d_basis = m_r_basis;
m_d_heading = m_r_heading;
m_d_nbasis = m_r_nbasis;
delete m_d_solver.m_factorization;
m_d_solver.m_factorization = nullptr;
} else {
prepare_solver_x_with_signature(signature, m_r_solver);
m_r_solver.start_tracing_basis_changes();
m_r_solver.find_feasible_solution();
if (settings().get_cancel_flag())
return;
m_r_solver.stop_tracing_basis_changes();
// and now catch up in the double solver
lean_assert(m_r_solver.total_iterations() >= m_r_solver.m_trace_of_basis_change_vector.size() /2);
catch_up_in_lu(m_r_solver.m_trace_of_basis_change_vector, m_r_solver.m_basis_heading, m_d_solver);
}
}
void create_double_matrix(static_matrix<double, double> & A) {
for (unsigned i = 0; i < m_r_A.row_count(); i++) {
auto & row = m_r_A.m_rows[i];
for (row_cell<mpq> & c : row) {
A.set(i, c.m_j, c.get_val().get_double());
}
}
}
void fill_basis_d(
vector<unsigned>& basis_d,
vector<int>& heading_d,
vector<unsigned>& nbasis_d){
basis_d = m_r_basis;
heading_d = m_r_heading;
nbasis_d = m_r_nbasis;
}
template <typename L, typename K>
void extract_signature_from_lp_core_solver(const lp_primal_core_solver<L, K> & solver, lar_solution_signature & signature) {
signature.clear();
lean_assert(signature.size() == 0);
for (unsigned j = 0; j < solver.m_basis_heading.size(); j++) {
if (solver.m_basis_heading[j] < 0) {
signature[j] = solver.get_non_basic_column_value_position(j);
}
}
}
void get_bounds_for_double_solver() {
unsigned n = m_n();
m_d_low_bounds.resize(n);
m_d_upper_bounds.resize(n);
double delta = find_delta_for_strict_boxed_bounds().get_double();
if (delta > 0.000001)
delta = 0.000001;
for (unsigned j = 0; j < n; j++) {
if (low_bound_is_set(j)) {
const auto & lb = m_r_solver.m_low_bounds[j];
m_d_low_bounds[j] = lb.x.get_double() + delta * lb.y.get_double();
}
if (upper_bound_is_set(j)) {
const auto & ub = m_r_solver.m_upper_bounds[j];
m_d_upper_bounds[j] = ub.x.get_double() + delta * ub.y.get_double();
lean_assert(!low_bound_is_set(j) || (m_d_upper_bounds[j] >= m_d_low_bounds[j]));
}
}
}
void scale_problem_for_doubles(
static_matrix<double, double>& A,
vector<double> & low_bounds,
vector<double> & upper_bounds) {
vector<double> column_scale_vector;
vector<double> right_side_vector(A.column_count());
settings().reps_in_scaler = 5;
scaler<double, double > scaler(right_side_vector,
A,
settings().scaling_minimum,
settings().scaling_maximum,
column_scale_vector,
settings());
if (! scaler.scale()) {
// the scale did not succeed, unscaling
A.clear();
create_double_matrix(A);
} else {
for (unsigned j = 0; j < A.column_count(); j++) {
if (m_r_solver.column_has_upper_bound(j)) {
upper_bounds[j] /= column_scale_vector[j];
}
if (m_r_solver.column_has_low_bound(j)) {
low_bounds[j] /= column_scale_vector[j];
}
}
}
}
// returns the trace of basis changes
vector<unsigned> find_solution_signature_with_doubles(lar_solution_signature & signature) {
if (m_d_solver.m_factorization == nullptr || m_d_solver.m_factorization->get_status() != LU_status::OK) {
vector<unsigned> ret;
return ret;
}
get_bounds_for_double_solver();
extract_signature_from_lp_core_solver(m_r_solver, signature);
prepare_solver_x_with_signature(signature, m_d_solver);
m_d_solver.start_tracing_basis_changes();
m_d_solver.find_feasible_solution();
if (settings().get_cancel_flag())
return vector<unsigned>();
m_d_solver.stop_tracing_basis_changes();
extract_signature_from_lp_core_solver(m_d_solver, signature);
return m_d_solver.m_trace_of_basis_change_vector;
}
bool low_bound_is_set(unsigned j) const {
switch (m_column_types[j]) {
case column_type::free_column:
case column_type::upper_bound:
return false;
case column_type::low_bound:
case column_type::boxed:
case column_type::fixed:
return true;
default:
lean_assert(false);
}
return false;
}
bool upper_bound_is_set(unsigned j) const {
switch (m_column_types[j]) {
case column_type::free_column:
case column_type::low_bound:
return false;
case column_type::upper_bound:
case column_type::boxed:
case column_type::fixed:
return true;
default:
lean_assert(false);
}
return false;
}
void update_delta(mpq& delta, numeric_pair<mpq> const& l, numeric_pair<mpq> const& u) const {
lean_assert(l <= u);
if (l.x < u.x && l.y > u.y) {
mpq delta1 = (u.x - l.x) / (l.y - u.y);
if (delta1 < delta) {
delta = delta1;
}
}
lean_assert(l.x + delta * l.y <= u.x + delta * u.y);
}
mpq find_delta_for_strict_boxed_bounds() const{
mpq delta = numeric_traits<mpq>::one();
for (unsigned j = 0; j < m_r_A.column_count(); j++ ) {
if (m_column_types()[j] != column_type::boxed)
continue;
update_delta(delta, m_r_low_bounds[j], m_r_upper_bounds[j]);
}
return delta;
}
mpq find_delta_for_strict_bounds(const mpq & initial_delta) const{
mpq delta = initial_delta;
for (unsigned j = 0; j < m_r_A.column_count(); j++ ) {
if (low_bound_is_set(j))
update_delta(delta, m_r_low_bounds[j], m_r_x[j]);
if (upper_bound_is_set(j))
update_delta(delta, m_r_x[j], m_r_upper_bounds[j]);
}
return delta;
}
void init_column_row_nz_for_r_solver() {
m_r_solver.init_column_row_non_zeroes();
}
linear_combination_iterator<mpq> * get_column_iterator(unsigned j) {
if (settings().use_tableau()) {
return new iterator_on_column<mpq, numeric_pair<mpq>>(m_r_solver.m_A.m_columns[j], m_r_solver.m_A);
} else {
m_r_solver.solve_Bd(j);
return new iterator_on_indexed_vector<mpq>(m_r_solver.m_ed);
}
}
};
}

293
src/util/lp/lar_core_solver.hpp Normal file
View file

@ -0,0 +1,293 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#include <string>
#include "util/vector.h"
#include "util/lp/lar_core_solver.h"
#include "util/lp/lar_solution_signature.h"
namespace lean {
lar_core_solver::lar_core_solver(
lp_settings & settings,
const column_namer & column_names
):
m_infeasible_sum_sign(0),
m_r_solver(m_r_A,
m_right_sides_dummy,
m_r_x,
m_r_basis,
m_r_nbasis,
m_r_heading,
m_costs_dummy,
m_column_types(),
m_r_low_bounds(),
m_r_upper_bounds(),
settings,
column_names),
m_d_solver(m_d_A,
m_d_right_sides_dummy,
m_d_x,
m_d_basis,
m_d_nbasis,
m_d_heading,
m_d_costs_dummy,
m_column_types(),
m_d_low_bounds,
m_d_upper_bounds,
settings,
column_names){}
void lar_core_solver::init_costs(bool first_time) {
lean_assert(false); // should not be called
// lean_assert(this->m_x.size() >= this->m_n());
// lean_assert(this->m_column_types.size() >= this->m_n());
// if (first_time)
// this->m_costs.resize(this->m_n());
// X inf = this->m_infeasibility;
// this->m_infeasibility = zero_of_type<X>();
// for (unsigned j = this->m_n(); j--;)
// init_cost_for_column(j);
// if (!(first_time || inf >= this->m_infeasibility)) {
// LP_OUT(this->m_settings, "iter = " << this->total_iterations() << std::endl);
// LP_OUT(this->m_settings, "inf was " << T_to_string(inf) << " and now " << T_to_string(this->m_infeasibility) << std::endl);
// lean_assert(false);
// }
// if (inf == this->m_infeasibility)
// this->m_iters_with_no_cost_growing++;
}
void lar_core_solver::init_cost_for_column(unsigned j) {
/*
// If j is a breakpoint column, then we set the cost zero.
// When anylyzing an entering column candidate we update the cost of the breakpoints columns to get the left or the right derivative if the infeasibility function
const numeric_pair<mpq> & x = this->m_x[j];
// set zero cost for each non-basis column
if (this->m_basis_heading[j] < 0) {
this->m_costs[j] = numeric_traits<T>::zero();
return;
}
// j is a basis column
switch (this->m_column_types[j]) {
case fixed:
case column_type::boxed:
if (x > this->m_upper_bounds[j]) {
this->m_costs[j] = 1;
this->m_infeasibility += x - this->m_upper_bounds[j];
} else if (x < this->m_low_bounds[j]) {
this->m_infeasibility += this->m_low_bounds[j] - x;
this->m_costs[j] = -1;
} else {
this->m_costs[j] = numeric_traits<T>::zero();
}
break;
case low_bound:
if (x < this->m_low_bounds[j]) {
this->m_costs[j] = -1;
this->m_infeasibility += this->m_low_bounds[j] - x;
} else {
this->m_costs[j] = numeric_traits<T>::zero();
}
break;
case upper_bound:
if (x > this->m_upper_bounds[j]) {
this->m_costs[j] = 1;
this->m_infeasibility += x - this->m_upper_bounds[j];
} else {
this->m_costs[j] = numeric_traits<T>::zero();
}
break;
case free_column:
this->m_costs[j] = numeric_traits<T>::zero();
break;
default:
lean_assert(false);
break;
}*/
}
// returns m_sign_of_alpha_r
int lar_core_solver::column_is_out_of_bounds(unsigned j) {
/*
switch (this->m_column_type[j]) {
case fixed:
case column_type::boxed:
if (this->x_below_low_bound(j)) {
return -1;
}
if (this->x_above_upper_bound(j)) {
return 1;
}
return 0;
case low_bound:
if (this->x_below_low_bound(j)) {
return -1;
}
return 0;
case upper_bound:
if (this->x_above_upper_bound(j)) {
return 1;
}
return 0;
default:
return 0;
break;
}*/
lean_assert(false);
return true;
}
void lar_core_solver::calculate_pivot_row(unsigned i) {
lean_assert(!m_r_solver.use_tableau());
lean_assert(m_r_solver.m_pivot_row.is_OK());
m_r_solver.m_pivot_row_of_B_1.clear();
m_r_solver.m_pivot_row_of_B_1.resize(m_r_solver.m_m());
m_r_solver.m_pivot_row.clear();
m_r_solver.m_pivot_row.resize(m_r_solver.m_n());
if (m_r_solver.m_settings.use_tableau()) {
unsigned basis_j = m_r_solver.m_basis[i];
for (auto & c : m_r_solver.m_A.m_rows[i]) {
if (c.m_j != basis_j)
m_r_solver.m_pivot_row.set_value(c.get_val(), c.m_j);
}
return;
}
m_r_solver.calculate_pivot_row_of_B_1(i);
m_r_solver.calculate_pivot_row_when_pivot_row_of_B1_is_ready(i);
}
void lar_core_solver::prefix_r() {
if (!m_r_solver.m_settings.use_tableau()) {
m_r_solver.m_copy_of_xB.resize(m_r_solver.m_n());
m_r_solver.m_ed.resize(m_r_solver.m_m());
m_r_solver.m_pivot_row.resize(m_r_solver.m_n());
m_r_solver.m_pivot_row_of_B_1.resize(m_r_solver.m_m());
m_r_solver.m_w.resize(m_r_solver.m_m());
m_r_solver.m_y.resize(m_r_solver.m_m());
m_r_solver.m_rows_nz.resize(m_r_solver.m_m(), 0);
m_r_solver.m_columns_nz.resize(m_r_solver.m_n(), 0);
init_column_row_nz_for_r_solver();
}
m_r_solver.m_b.resize(m_r_solver.m_m());
if (m_r_solver.m_settings.simplex_strategy() != simplex_strategy_enum::tableau_rows) {
if(m_r_solver.m_settings.use_breakpoints_in_feasibility_search)
m_r_solver.m_breakpoint_indices_queue.resize(m_r_solver.m_n());
m_r_solver.m_costs.resize(m_r_solver.m_n());
m_r_solver.m_d.resize(m_r_solver.m_n());
m_r_solver.m_using_infeas_costs = true;
}
}
void lar_core_solver::prefix_d() {
m_d_solver.m_b.resize(m_d_solver.m_m());
m_d_solver.m_breakpoint_indices_queue.resize(m_d_solver.m_n());
m_d_solver.m_copy_of_xB.resize(m_d_solver.m_n());
m_d_solver.m_costs.resize(m_d_solver.m_n());
m_d_solver.m_d.resize(m_d_solver.m_n());
m_d_solver.m_ed.resize(m_d_solver.m_m());
m_d_solver.m_pivot_row.resize(m_d_solver.m_n());
m_d_solver.m_pivot_row_of_B_1.resize(m_d_solver.m_m());
m_d_solver.m_w.resize(m_d_solver.m_m());
m_d_solver.m_y.resize(m_d_solver.m_m());
m_d_solver.m_steepest_edge_coefficients.resize(m_d_solver.m_n());
m_d_solver.m_column_norms.clear();
m_d_solver.m_column_norms.resize(m_d_solver.m_n(), 2);
m_d_solver.m_inf_set.clear();
m_d_solver.m_inf_set.resize(m_d_solver.m_n());
}
void lar_core_solver::fill_not_improvable_zero_sum_from_inf_row() {
lean_assert(m_r_solver.A_mult_x_is_off() == false);
unsigned bj = m_r_basis[m_r_solver.m_inf_row_index_for_tableau];
m_infeasible_sum_sign = m_r_solver.inf_sign_of_column(bj);
m_infeasible_linear_combination.clear();
for (auto & rc : m_r_solver.m_A.m_rows[m_r_solver.m_inf_row_index_for_tableau]) {
m_infeasible_linear_combination.push_back(std::make_pair( rc.get_val(), rc.m_j));
}
}
void lar_core_solver::fill_not_improvable_zero_sum() {
if (m_r_solver.m_settings.simplex_strategy() == simplex_strategy_enum::tableau_rows) {
fill_not_improvable_zero_sum_from_inf_row();
return;
}
// reusing the existing mechanism for row_feasibility_loop
m_infeasible_sum_sign = m_r_solver.m_settings.use_breakpoints_in_feasibility_search? -1 : 1;
m_infeasible_linear_combination.clear();
for (auto j : m_r_solver.m_basis) {
const mpq & cost_j = m_r_solver.m_costs[j];
if (!numeric_traits<mpq>::is_zero(cost_j)) {
m_infeasible_linear_combination.push_back(std::make_pair(cost_j, j));
}
}
// m_costs are expressed by m_d ( additional costs), substructing the latter gives 0
for (unsigned j = 0; j < m_r_solver.m_n(); j++) {
if (m_r_solver.m_basis_heading[j] >= 0) continue;
const mpq & d_j = m_r_solver.m_d[j];
if (!numeric_traits<mpq>::is_zero(d_j)) {
m_infeasible_linear_combination.push_back(std::make_pair(-d_j, j));
}
}
}
void lar_core_solver::solve() {
lean_assert(m_r_solver.non_basic_columns_are_set_correctly());
lean_assert(m_r_solver.inf_set_is_correct());
if (m_r_solver.current_x_is_feasible() && m_r_solver.m_look_for_feasible_solution_only) {
m_r_solver.set_status(OPTIMAL);
return;
}
++settings().st().m_need_to_solve_inf;
lean_assert(!m_r_solver.A_mult_x_is_off());
lean_assert((!settings().use_tableau()) || r_basis_is_OK());
if (need_to_presolve_with_double_solver()) {
prefix_d();
lar_solution_signature solution_signature;
vector<unsigned> changes_of_basis = find_solution_signature_with_doubles(solution_signature);
if (m_d_solver.get_status() == TIME_EXHAUSTED) {
m_r_solver.set_status(TIME_EXHAUSTED);
return;
}
if (settings().use_tableau())
solve_on_signature_tableau(solution_signature, changes_of_basis);
else
solve_on_signature(solution_signature, changes_of_basis);
lean_assert(!settings().use_tableau() || r_basis_is_OK());
} else {
if (!settings().use_tableau()) {
bool snapped = m_r_solver.snap_non_basic_x_to_bound();
lean_assert(m_r_solver.non_basic_columns_are_set_correctly());
if (snapped)
m_r_solver.solve_Ax_eq_b();
}
if (m_r_solver.m_look_for_feasible_solution_only)
m_r_solver.find_feasible_solution();
else
m_r_solver.solve();
lean_assert(!settings().use_tableau() || r_basis_is_OK());
}
if (m_r_solver.get_status() == INFEASIBLE) {
fill_not_improvable_zero_sum();
} else if (m_r_solver.get_status() != UNBOUNDED) {
m_r_solver.set_status(OPTIMAL);
}
lean_assert(r_basis_is_OK());
lean_assert(m_r_solver.non_basic_columns_are_set_correctly());
lean_assert(m_r_solver.inf_set_is_correct());
}
}

10
src/util/lp/lar_core_solver_instances.cpp Normal file
View file

@ -0,0 +1,10 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#include <utility>
#include <memory>
#include <string>
#include "util/vector.h"
#include <functional>
#include "util/lp/lar_core_solver.hpp"

13
src/util/lp/lar_solution_signature.h Normal file
View file

@ -0,0 +1,13 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#pragma once
#include "util/vector.h"
#include "util/debug.h"
#include "util/lp/lp_settings.h"
#include <unordered_map>
namespace lean {
typedef std::unordered_map<unsigned, non_basic_column_value_position> lar_solution_signature;
}

1549
src/util/lp/lar_solver.h Normal file
View file

File diff suppressed because it is too large Load diff

Some files were not shown because too many files have changed in this diff Show more