1import pyomo.common.unittest as unittest 2import pyomo.environ as pe 3from pyomo.core.base import ConcreteModel, Var, Constraint, Objective 4from pyomo.common.dependencies import attempt_import 5 6np, numpy_available = attempt_import('numpy', 'Interior point requires numpy', 7 minimum_version='1.13.0') 8scipy, scipy_available = attempt_import('scipy', 'Interior point requires scipy') 9mumps, mumps_available = attempt_import('mumps', 'Interior point requires mumps') 10 11if not (numpy_available and scipy_available): 12 raise unittest.SkipTest('Interior point tests require numpy and scipy') 13 14from pyomo.contrib.pynumero.asl import AmplInterface 15asl_available = AmplInterface.available() 16if not asl_available: 17 raise unittest.SkipTest('Regularization tests require ASL') 18from pyomo.contrib.interior_point.interface import InteriorPointInterface 19from pyomo.contrib.interior_point.interior_point import InteriorPointSolver 20from pyomo.contrib.interior_point.linalg.mumps_interface import MumpsInterface 21 22 23def make_model_tri(n, small_val=1e-7, big_val=1e2): 24 m = ConcreteModel() 25 m.x = Var(range(n), initialize=0.5) 26 27 def c_rule(m, i): 28 return big_val*m.x[i-1] + small_val*m.x[i] + big_val*m.x[i+1] == 1 29 30 m.c = Constraint(range(1,n-1), rule=c_rule) 31 32 m.obj = Objective(expr=small_val*sum((m.x[i]-1)**2 for i in range(n))) 33 34 return m 35 36class TestReallocation(unittest.TestCase): 37 def _test_ip_with_reallocation(self, linear_solver, interface): 38 ip_solver = InteriorPointSolver(linear_solver, 39 max_reallocation_iterations=3, 40 reallocation_factor=1.1, 41 # The small factor is to ensure that multiple iterations of 42 # reallocation are performed. The bug in the previous 43 # implementation only occurred if 2+ reallocation iterations 44 # were needed (max_reallocation_iterations >= 3). 45 max_iter=1) 46 ip_solver.set_interface(interface) 47 48 ip_solver.solve(interface) 49 50 return ip_solver 51 52 @unittest.skipIf(not mumps_available, 'Mumps is not available') 53 def test_mumps(self): 54 n = 20000 55 m = make_model_tri(n, small_val=1e-7) 56 interface = InteriorPointInterface(m) 57 linear_solver = MumpsInterface() 58 # Default memory "buffer" factor: 20 59 linear_solver.set_icntl(14, 20) 60 61 kkt = interface.evaluate_primal_dual_kkt_matrix() 62 res = linear_solver.do_symbolic_factorization(kkt) 63 predicted = linear_solver.get_infog(16) 64 65 self._test_ip_with_reallocation(linear_solver, interface) 66 actual = linear_solver.get_icntl(23) 67 68 self.assertTrue(predicted == 12 or predicted == 11) 69 self.assertTrue(actual > predicted) 70 #self.assertEqual(actual, 14) 71 # NOTE: This test will break if Mumps (or your Mumps version) 72 # gets more conservative at estimating memory requirement, 73 # or if the numeric factorization gets more efficient. 74 75 76if __name__ == '__main__': 77 # 78 unittest.main() 79