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