451868044

Build Type

push

github

Committed by Jeffrey Hokanson

Commit Message

Bumped version

Run Details

307 of 377 branches covered (81.43%)

Branch coverage included in aggregate %.

1 of 1 new or added line in 1 file covered. (100.0%)

2 existing lines in 1 file now uncovered.

1311 of 1422 relevant lines covered (92.19%)

2.73 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

56.95

/polyrat/slp.py

""" Sequential linear programing approach for minimax optimization

"""

import numpy as np
from iterprinter import IterationPrinter
from scipy.optimize import OptimizeResult, approx_fprime, linprog, Bounds
from scipy.optimize._constraints import new_bounds_to_old
import cvxopt as co

def slp(fun, x0, args = (), jac = None, bounds = None, constraints = (), hess = None, hessp = None, verbose = True,
        maxiter = 100, xtol = 1e-6, c_armijo = 0.01, bt_maxiter = 30, **options):
        r""" Sequential linear programming approach for constrained optimization problems

        Solves 
        
        min_x  f(x)
        s.t.   lbc <= c(x) <= ubc 
           lb <= x <= ub

        """

        nfev = 0
        njev = 0

        if jac is None or jac is False:
                assert False

        # Normalize constraints format as a list
        if not isinstance(constraints, list):
                constraints = [constraints]

        if verbose:
                iterprinter = IterationPrinter(it = '4d', obj = '20.14e', alpha = '8.3e', delta = '8.3e')
                iterprinter.print_header(it = 'iter', obj = "objective fun.", alpha = 'step length', delta = 'trust region radius')

        x = np.copy(x0)
        delta = 1e0

        stop = False
        it = 0

        fx = fun(x)


        if verbose:
                iterprinter.print_iter(it = it, obj = fx)

        while it < maxiter and not stop:
                it += 1
                # Generate linearized constraints
                Acon = []
                bcon = []
                for con in constraints:
                        fx = con.fun(x)
                        Jx = con.jac(x)

                        # fx + Jx @ p <= ub
                        I = np.isfinite(con.ub)
                        Acon.append(Jx[I])
                        bcon.append(con.ub[I] - fx[I])
                        # lb <= fx + Jx @ p
                        I = np.isfinite(con.lb)
                        Acon.append(-Jx[I])
                        bcon.append(fx[I] - con.lb[I])        


                fx = fun(x)
                # Generate search direction
                g = jac(x)        
                
                # Convert the bounds into a format for linprog
                if isinstance(bounds, Bounds):
                        #linprog_bounds = [(lbi, ubi) for lbi, ubi in zip(bounds.lb, bounds.ub)]
                        linprog_bounds = new_bounds_to_old(np.maximum(bounds.lb - x, -delta), np.minimum(bounds.ub - x, delta), len(x0)) 
                elif bounds is None:
                        # Since Linprog assumes non-negative, we have to specify infinite lower bounds 
                        linprog_bounds = [(None, None) for i in enumerate(x0)]        
                else:
                        raise NotImplementedError
                        #linprog_bounds = bounds

                # Solve LP
                if False:
                        res = linprog(g, A_ub = np.vstack(Acon), b_ub = np.hstack(bcon), bounds = linprog_bounds)
                        if not res.success:
                                if verbose: print("linear program failed with message", res.message)
                                return OptimizeResult(
                                                x = x, 
                                                success = False,
                                                status = -1,
                                                message = 'LP solve failed with message ' + res.message,
                                                fun = fx,
                                                jac = g,         
                                                nit = it,
                                                nfev = nfev,
                                                njev = njev,)

                        p = res.x
                else:
                        co.solvers.options['show_progress'] = False
                        G = co.matrix(np.vstack(Acon + [np.eye(len(x)), -np.eye(len(x))] ))
                        h = co.matrix(np.hstack(bcon + [b[1] for b in linprog_bounds] + [-b[0] for b in linprog_bounds])) 
                        sol = co.solvers.lp(co.matrix(g),G, h)
                        #print(sol)
                        p = np.array(sol['x']).flatten()

        
                # Backtracking line search        
                pred_decrease = g.T @ p
                for alpha in [0.5**k for k in range(bt_maxiter)]:
                        xn = x + alpha*p
                        fxn = fun(xn)
                        # Check if satisfy Armijo conditions
                        improved = (fxn        <= fx + c_armijo*alpha*pred_decrease)        
                        feasible = True
                        for con in constraints:
                                cxn = con.fun(xn)
                                feasible = feasible & (np.all(cxn <= con.ub) and np.all(con.lb <= cxn))
                        
                        if feasible and improved:
                                break

                if not np.isclose(alpha,1):
                        delta = np.max(np.abs(alpha*p))
                else:
                        delta *=2
                

                if not (feasible and improved):
                        return OptimizeResult(
                                        x = x, 
                                        success = False,
                                        status = -2,
                                        message = 'unable to find a successful step',
                                        fun = fx,
                                        jac = g,         
                                        nit = it,
                                        nfev = nfev,
                                        njev = njev,)
                                
                if verbose:
                        iterprinter.print_iter(it=it, obj = fxn, alpha =alpha, delta = delta)

                if np.max(np.abs(alpha*p)) < xtol:
                        if verbose: print("small change in parameters")
                        return OptimizeResult(
                                        x = x, 
                                        success = True,
                                        status = 0,
                                        message = 'small change in objective',
                                        fun = fx,
                                        jac = g,         
                                        nit = it,
                                        nfev = nfev,
                                        njev = njev,)
                x = xn
                        
        return OptimizeResult(
                        x = x, 
                        success = True,
                        status = 0,
                        message = 'maximum number of iterations exceeded',
                        fun = fx,
                        jac = g,         
                        nit = it,
                        nfev = nfev,
                        njev = njev,)


def sl1lp(fun, x0, args = (), jac = None, bounds = None, constraints = (), verbose = True,
        maxiter = 100, **options):
        r"""

        Solve a constrained optimization problem in the l1 relaxation (see, Nocedal and Wright Sec. 18.5)

        """ 

        if jac is None or jac is False:
                assert False

        # Convert the bounds into a format for linprog
        if isinstance(bounds, Bounds):
                #linprog_bounds = [(lbi, ubi) for lbi, ubi in zip(bounds.lb, bounds.ub)]
                linprog_bounds = new_bounds_to_old(bounds.lb, bounds.ub, len(x0)) 
        elif bounds is None:
                # Since Linprog assumes non-negative, we have to specify infinite lower bounds 
                linprog_bounds = [(None, None) for i in enumerate(x0)]        
        else:
                linprog_bounds = bounds



        # Normalize constraints format as a list
        if not isinstance(constraints, list):
                constraints = [constraints]

        if verbose:
                printer = IterationPrinter(it = '4d', obj = '20.8e',  con= '20.8e', alpha = '8.3e')
                printer.print_header(it = 'iter', obj = "objective fun.", con = 'constraint violation',  alpha = 'step length')


        x = np.copy(x0)
        fx = fun(x)
        con_violation = 0
        for con in constraints:
                cx = con.fun(x)
                con_violation += np.sum(np.maximum(con.lb - cx, 0))
                con_violation += np.sum(np.maximum(cx - con.ub, 0))


        if verbose:
                printer.print_iter(it = 0, obj = fx, con = con_violation)
        mu = 1
        it = 0
        stop = False

        while it <= maxiter and not stop:
                it += 1
                
                assert False

1	""" Sequential linear programing approach for minimax optimization
2
3	"""
4
5	import numpy as np	3✔
6	from iterprinter import IterationPrinter	3✔
7	from scipy.optimize import OptimizeResult, approx_fprime, linprog, Bounds	3✔
8	from scipy.optimize._constraints import new_bounds_to_old	3✔
9	import cvxopt as co	3✔
10
11	def slp(fun, x0, args = (), jac = None, bounds = None, constraints = (), hess = None, hessp = None, verbose = True,	3✔
12	maxiter = 100, xtol = 1e-6, c_armijo = 0.01, bt_maxiter = 30, **options):
13	r""" Sequential linear programming approach for constrained optimization problems
14
15	Solves
16
17	min_x f(x)
18	s.t. lbc <= c(x) <= ubc
19	lb <= x <= ub
20
21	"""
22
23	nfev = 0	3✔
24	njev = 0	3✔
25
26	if jac is None or jac is False:	3!
27	assert False	×
28
29	# Normalize constraints format as a list
30	if not isinstance(constraints, list):	3!
31	constraints = [constraints]	3✔
32
33	if verbose:	3!
34	iterprinter = IterationPrinter(it = '4d', obj = '20.14e', alpha = '8.3e', delta = '8.3e')	3✔
35	iterprinter.print_header(it = 'iter', obj = "objective fun.", alpha = 'step length', delta = 'trust region radius')	3✔
36
37	x = np.copy(x0)	3✔
38	delta = 1e0	3✔
39
40	stop = False	3✔
41	it = 0	3✔
42
43	fx = fun(x)	3✔
44
45
46	if verbose:	3!
47	iterprinter.print_iter(it = it, obj = fx)	3✔
48
49	while it < maxiter and not stop:	3✔
50	it += 1	3✔
51	# Generate linearized constraints
52	Acon = []	3✔
53	bcon = []	3✔
54	for con in constraints:	3✔
55	fx = con.fun(x)	3✔
56	Jx = con.jac(x)	3✔
57
58	# fx + Jx @ p <= ub
59	I = np.isfinite(con.ub)	3✔
60	Acon.append(Jx[I])	3✔
61	bcon.append(con.ub[I] - fx[I])	3✔
62	# lb <= fx + Jx @ p
63	I = np.isfinite(con.lb)	3✔
64	Acon.append(-Jx[I])	3✔
65	bcon.append(fx[I] - con.lb[I])	3✔
66
67
68	fx = fun(x)	3✔
69	# Generate search direction
70	g = jac(x)	3✔
71
72	# Convert the bounds into a format for linprog
73	if isinstance(bounds, Bounds):	3!
74	#linprog_bounds = [(lbi, ubi) for lbi, ubi in zip(bounds.lb, bounds.ub)]
75	linprog_bounds = new_bounds_to_old(np.maximum(bounds.lb - x, -delta), np.minimum(bounds.ub - x, delta), len(x0))	3✔
76	elif bounds is None:	×
77	# Since Linprog assumes non-negative, we have to specify infinite lower bounds
78	linprog_bounds = [(None, None) for i in enumerate(x0)]	×
79	else:
80	raise NotImplementedError	×
81	#linprog_bounds = bounds
82
83	# Solve LP
84	if False:
85	res = linprog(g, A_ub = np.vstack(Acon), b_ub = np.hstack(bcon), bounds = linprog_bounds)
86	if not res.success:
87	if verbose: print("linear program failed with message", res.message)
88	return OptimizeResult(
89	x = x,
90	success = False,
91	status = -1,
92	message = 'LP solve failed with message ' + res.message,
93	fun = fx,
94	jac = g,
95	nit = it,
96	nfev = nfev,
97	njev = njev,)
98
99	p = res.x
100	else:
101	co.solvers.options['show_progress'] = False	3✔
102	G = co.matrix(np.vstack(Acon + [np.eye(len(x)), -np.eye(len(x))] ))	3✔
103	h = co.matrix(np.hstack(bcon + [b[1] for b in linprog_bounds] + [-b[0] for b in linprog_bounds]))	3✔
104	sol = co.solvers.lp(co.matrix(g),G, h)	3✔
105	#print(sol)
106	p = np.array(sol['x']).flatten()	3✔
107
108
109	# Backtracking line search
110	pred_decrease = g.T @ p	3✔
111	for alpha in [0.5**k for k in range(bt_maxiter)]:	3!
112	xn = x + alpha*p	3✔
113	fxn = fun(xn)	3✔
114	# Check if satisfy Armijo conditions
115	improved = (fxn <= fx + c_armijoalphapred_decrease)	3✔
116	feasible = True	3✔
117	for con in constraints:	3✔
118	cxn = con.fun(xn)	3✔
119	feasible = feasible & (np.all(cxn <= con.ub) and np.all(con.lb <= cxn))	3✔
120
121	if feasible and improved:	3✔
122	break	3✔
123
124	if not np.isclose(alpha,1):	3✔
125	delta = np.max(np.abs(alpha*p))	3✔
126	else:
127	delta *=2	3✔
128
129
130	if not (feasible and improved):	3!
131	return OptimizeResult(	×
132	x = x,
133	success = False,
134	status = -2,
135	message = 'unable to find a successful step',
136	fun = fx,
137	jac = g,
138	nit = it,
139	nfev = nfev,
140	njev = njev,)
141
142	if verbose:	3!
143	iterprinter.print_iter(it=it, obj = fxn, alpha =alpha, delta = delta)	3✔
144
145	if np.max(np.abs(alpha*p)) < xtol:	3!
UNCOV 146	if verbose: print("small change in parameters")	×
UNCOV 147	return OptimizeResult(	×
148	x = x,
149	success = True,
150	status = 0,
151	message = 'small change in objective',
152	fun = fx,
153	jac = g,
154	nit = it,
155	nfev = nfev,
156	njev = njev,)
157	x = xn	3✔
158
159	return OptimizeResult(	3✔
160	x = x,
161	success = True,
162	status = 0,
163	message = 'maximum number of iterations exceeded',
164	fun = fx,
165	jac = g,
166	nit = it,
167	nfev = nfev,
168	njev = njev,)
169
170
171	def sl1lp(fun, x0, args = (), jac = None, bounds = None, constraints = (), verbose = True,	3✔
172	maxiter = 100, **options):
173	r"""
174
175	Solve a constrained optimization problem in the l1 relaxation (see, Nocedal and Wright Sec. 18.5)
176
177	"""
178
179	if jac is None or jac is False:	×
180	assert False	×
181
182	# Convert the bounds into a format for linprog
183	if isinstance(bounds, Bounds):	×
184	#linprog_bounds = [(lbi, ubi) for lbi, ubi in zip(bounds.lb, bounds.ub)]
185	linprog_bounds = new_bounds_to_old(bounds.lb, bounds.ub, len(x0))	×
186	elif bounds is None:	×
187	# Since Linprog assumes non-negative, we have to specify infinite lower bounds
188	linprog_bounds = [(None, None) for i in enumerate(x0)]	×
189	else:
190	linprog_bounds = bounds	×
191
192
193
194	# Normalize constraints format as a list
195	if not isinstance(constraints, list):	×
196	constraints = [constraints]	×
197
198	if verbose:	×
199	printer = IterationPrinter(it = '4d', obj = '20.8e', con= '20.8e', alpha = '8.3e')	×
200	printer.print_header(it = 'iter', obj = "objective fun.", con = 'constraint violation', alpha = 'step length')	×
201
202
203	x = np.copy(x0)	×
204	fx = fun(x)	×
205	con_violation = 0	×
206	for con in constraints:	×
207	cx = con.fun(x)	×
208	con_violation += np.sum(np.maximum(con.lb - cx, 0))	×
209	con_violation += np.sum(np.maximum(cx - con.ub, 0))	×
210
211
212	if verbose:	×
213	printer.print_iter(it = 0, obj = fx, con = con_violation)	×
214	mu = 1	×
215	it = 0	×
216	stop = False	×
217
218	while it <= maxiter and not stop:	×
219	it += 1	×
220
221	assert False	×

jeffrey-hokanson / polyrat / 451868044

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous