From 5a476db6d2a907b4de9a42a141936bdbf32df3b7 Mon Sep 17 00:00:00 2001 From: alex Date: Wed, 28 Nov 2012 12:45:08 -0500 Subject: [PATCH 1/6] add ipopt to the three minimization examples, but the .rst files are not updated with the new output because the ipopt output is too verbose with the current settings --- .../sphinx/examples/minimization/minhelper.py | 150 ++++++++++++++++++ 1 file changed, 150 insertions(+) diff --git a/documentation/sphinx/examples/minimization/minhelper.py b/documentation/sphinx/examples/minimization/minhelper.py index ee409ba..0c4a504 100644 --- a/documentation/sphinx/examples/minimization/minhelper.py +++ b/documentation/sphinx/examples/minimization/minhelper.py @@ -8,6 +8,129 @@ import scipy.optimize import algopy import numdifftools +import pyipopt + + +####################################################################### +# FIXME: +# The variables and functions in this section are ipopt-specific +# and should be moved somewhere else. + + +# http://www.coin-or.org/Ipopt/documentation/node35.html +g_ipopt_pos_inf = 2e19 +g_ipopt_neg_inf = -g_ipopt_pos_inf + +def my_ipopt_dummy_eval_g(X, user_data=None): + return numpy.array([], dtype=float) + +def my_ipopt_dummy_eval_jac_g(X, flag, user_data=None): + """ + We assume that the optimization is unconstrained. + Therefore the constraint Jacobian is not so interesting, + but we still need it to exist so that ipopt does not complain. + """ + rows = numpy.array([], dtype=int) + cols = numpy.array([], dtype=int) + if flag: + return (rows, cols) + else: + raise Exception('this should not be called') + #return numpy.array([], dtype=float) + +def my_ipopt_eval_h( + h, nvar, + X, lagrange, obj_factor, flag, user_data=None): + """ + Evaluate the sparse hessian of the Lagrangian. + The first group of parameters should be applied using functools.partial. + The second group of parameters are passed from ipopt. + @param h: a function to compute the hessian. + @param nvar: the number of parameters + @param X: parameter values + @param lagrange: something about the constraints + @param obj_factor: no clue what this is + @param flag: this asks for the sparsity structure + """ + + # Get the nonzero (row, column) entries of a lower triangular matrix. + # This is related to the fact that the Hessian is symmetric, + # and that ipopt is designed to work with sparse matrices. + row_list = [] + col_list = [] + for row in range(nvar): + for col in range(row+1): + row_list.append(row) + col_list.append(col) + rows = numpy.array(row_list, dtype=int) + cols = numpy.array(col_list, dtype=int) + + if flag: + return (rows, cols) + else: + if nvar != len(X): + raise Exception('parameter count mismatch') + if lagrange: + raise Exception('only unconstrained is implemented for now...') + values = numpy.zeros(len(rows), dtype=float) + H = h(X) + for i, (r, c) in enumerate(zip(rows, cols)): + #FIXME: am I using obj_factor correctly? + # I don't really know what it is... + values[i] = H[r, c] * obj_factor + return values + +def my_ipopt_apply_new(X): + #FIXME: I don't really know what this does, but ipopt wants it. + return True + +def my_fmin_ipopt(f, x0, fprime, fhess=None): + nvar = len(x0) + x_L = numpy.array([g_ipopt_neg_inf]*nvar, dtype=float) + x_U = numpy.array([g_ipopt_pos_inf]*nvar, dtype=float) + ncon = 0 + g_L = numpy.array([], dtype=float) + g_U = numpy.array([], dtype=float) + nnzj = 0 + nnzh = (nvar * (nvar + 1)) // 2 + + # define the callbacks + eval_f = f + eval_grad_f = fprime + eval_g = my_ipopt_dummy_eval_g + eval_jac_g = my_ipopt_dummy_eval_jac_g + if fhess: + eval_h = functools.partial(my_ipopt_eval_h, fhess, nvar) + apply_new = my_ipopt_apply_new + + # compute the results using ipopt + nlp_args = [ + nvar, + x_L, + x_U, + ncon, + g_L, + g_U, + nnzj, + nnzh, + eval_f, + eval_grad_f, + eval_g, + eval_jac_g, + ] + if fhess: + nlp_args.extend([ + eval_h, + apply_new, + ]) + nlp = pyipopt.create(*nlp_args) + results = nlp.solve(x0) + nlp.close() + return results + + +####################################################################### +# The following two functions are algopy-specific. def eval_grad(f, theta): @@ -18,6 +141,8 @@ def eval_hess(f, theta): theta = algopy.UTPM.init_hessian(theta) return algopy.UTPM.extract_hessian(len(theta), f(theta)) + + def show_local_curvature(f, g, h, x0): print 'point:' print x0 @@ -169,6 +294,31 @@ def do_searches(f, g, h, x0): print results print + print 'strategy:', 'ipopt' + print 'options:', 'default' + print 'gradient:', 'autodiff' + print 'hessian:', 'autodiff' + results = my_fmin_ipopt( + f, + x0, + fprime=g, + fhess=h, + ) + print results + print + + print 'strategy:', 'ipopt' + print 'options:', 'default' + print 'gradient:', 'autodiff' + print 'hessian:', 'finite differences' + results = my_fmin_ipopt( + f, + x0, + fprime=g, + ) + print results + print + def show_minimization_results(f, target_in, easy_init_in, hard_init_in): """ From 299596e7561681e9f6f1b01440306ea60d8f2d6c Mon Sep 17 00:00:00 2001 From: alex Date: Wed, 28 Nov 2012 13:10:45 -0500 Subject: [PATCH 2/6] add overlooked gammaln pullback test --- algopy/utpm/tests/test_utpm.py | 13 +++++++++++++ algopy/utpm/utpm.py | 2 +- 2 files changed, 14 insertions(+), 1 deletion(-) diff --git a/algopy/utpm/tests/test_utpm.py b/algopy/utpm/tests/test_utpm.py index 810a8f7..f263825 100644 --- a/algopy/utpm/tests/test_utpm.py +++ b/algopy/utpm/tests/test_utpm.py @@ -653,6 +653,19 @@ def test_gammaln(self): b = UTPM.gammaln(x + 1) - UTPM.log(x) assert_allclose(a.data, b.data) + def test_gammaln_pullback(self): + D,P = 2,1 + + # forward + x = UTPM(numpy.random.random((D,P))) + y = UTPM.gammaln(x) + + # reverse + ybar = UTPM(numpy.random.random((D,P))) + xbar = UTPM.pb_gammaln(ybar, x, y) + + assert_array_almost_equal(ybar.data[0]*y.data[1], xbar.data[0]*x.data[1]) + def test_hyperu(self): D,P,N,M = 5,1,3,3 diff --git a/algopy/utpm/utpm.py b/algopy/utpm/utpm.py index 1a6a967..b6de34f 100644 --- a/algopy/utpm/utpm.py +++ b/algopy/utpm/utpm.py @@ -831,7 +831,7 @@ def pb_gammaln(cls, ybar, x, y, out=None): xbar, = out cls._pb_gammaln(ybar.data, x.data, y.data, out = xbar.data) - return out + return xbar @classmethod def erf(cls, x): From cb7d7bfdcd418d314a1910d2ad91c4c2ae27002d Mon Sep 17 00:00:00 2001 From: alex Date: Wed, 28 Nov 2012 18:39:16 -0500 Subject: [PATCH 3/6] push a bunch of ipopt interface code back into my branch of pyipopt --- .../sphinx/examples/minimization/minhelper.py | 128 +----------------- 1 file changed, 5 insertions(+), 123 deletions(-) diff --git a/documentation/sphinx/examples/minimization/minhelper.py b/documentation/sphinx/examples/minimization/minhelper.py index 0c4a504..3339e40 100644 --- a/documentation/sphinx/examples/minimization/minhelper.py +++ b/documentation/sphinx/examples/minimization/minhelper.py @@ -10,127 +10,9 @@ import numdifftools import pyipopt - -####################################################################### -# FIXME: -# The variables and functions in this section are ipopt-specific -# and should be moved somewhere else. - - -# http://www.coin-or.org/Ipopt/documentation/node35.html -g_ipopt_pos_inf = 2e19 -g_ipopt_neg_inf = -g_ipopt_pos_inf - -def my_ipopt_dummy_eval_g(X, user_data=None): - return numpy.array([], dtype=float) - -def my_ipopt_dummy_eval_jac_g(X, flag, user_data=None): - """ - We assume that the optimization is unconstrained. - Therefore the constraint Jacobian is not so interesting, - but we still need it to exist so that ipopt does not complain. - """ - rows = numpy.array([], dtype=int) - cols = numpy.array([], dtype=int) - if flag: - return (rows, cols) - else: - raise Exception('this should not be called') - #return numpy.array([], dtype=float) - -def my_ipopt_eval_h( - h, nvar, - X, lagrange, obj_factor, flag, user_data=None): - """ - Evaluate the sparse hessian of the Lagrangian. - The first group of parameters should be applied using functools.partial. - The second group of parameters are passed from ipopt. - @param h: a function to compute the hessian. - @param nvar: the number of parameters - @param X: parameter values - @param lagrange: something about the constraints - @param obj_factor: no clue what this is - @param flag: this asks for the sparsity structure - """ - - # Get the nonzero (row, column) entries of a lower triangular matrix. - # This is related to the fact that the Hessian is symmetric, - # and that ipopt is designed to work with sparse matrices. - row_list = [] - col_list = [] - for row in range(nvar): - for col in range(row+1): - row_list.append(row) - col_list.append(col) - rows = numpy.array(row_list, dtype=int) - cols = numpy.array(col_list, dtype=int) - - if flag: - return (rows, cols) - else: - if nvar != len(X): - raise Exception('parameter count mismatch') - if lagrange: - raise Exception('only unconstrained is implemented for now...') - values = numpy.zeros(len(rows), dtype=float) - H = h(X) - for i, (r, c) in enumerate(zip(rows, cols)): - #FIXME: am I using obj_factor correctly? - # I don't really know what it is... - values[i] = H[r, c] * obj_factor - return values - -def my_ipopt_apply_new(X): - #FIXME: I don't really know what this does, but ipopt wants it. - return True - -def my_fmin_ipopt(f, x0, fprime, fhess=None): - nvar = len(x0) - x_L = numpy.array([g_ipopt_neg_inf]*nvar, dtype=float) - x_U = numpy.array([g_ipopt_pos_inf]*nvar, dtype=float) - ncon = 0 - g_L = numpy.array([], dtype=float) - g_U = numpy.array([], dtype=float) - nnzj = 0 - nnzh = (nvar * (nvar + 1)) // 2 - - # define the callbacks - eval_f = f - eval_grad_f = fprime - eval_g = my_ipopt_dummy_eval_g - eval_jac_g = my_ipopt_dummy_eval_jac_g - if fhess: - eval_h = functools.partial(my_ipopt_eval_h, fhess, nvar) - apply_new = my_ipopt_apply_new - - # compute the results using ipopt - nlp_args = [ - nvar, - x_L, - x_U, - ncon, - g_L, - g_U, - nnzj, - nnzh, - eval_f, - eval_grad_f, - eval_g, - eval_jac_g, - ] - if fhess: - nlp_args.extend([ - eval_h, - apply_new, - ]) - nlp = pyipopt.create(*nlp_args) - results = nlp.solve(x0) - nlp.close() - return results - - -####################################################################### -# The following two functions are algopy-specific. +# Suppress log spam from pyipopt. +# But ipopt itself will stil spam... +pyipopt.set_loglevel(0) def eval_grad(f, theta): @@ -298,7 +180,7 @@ def do_searches(f, g, h, x0): print 'options:', 'default' print 'gradient:', 'autodiff' print 'hessian:', 'autodiff' - results = my_fmin_ipopt( + results = pyipopt.fmin_unconstrained( f, x0, fprime=g, @@ -311,7 +193,7 @@ def do_searches(f, g, h, x0): print 'options:', 'default' print 'gradient:', 'autodiff' print 'hessian:', 'finite differences' - results = my_fmin_ipopt( + results = pyipopt.fmin_unconstrained( f, x0, fprime=g, From c5c2e17733b555f9be9d6e3a994bb1b609fbdcfd Mon Sep 17 00:00:00 2001 From: alex Date: Wed, 5 Dec 2012 10:07:23 -0500 Subject: [PATCH 4/6] comment out the pyipopt stuff --- .../sphinx/examples/minimization/minhelper.py | 59 +++++++++++-------- 1 file changed, 33 insertions(+), 26 deletions(-) diff --git a/documentation/sphinx/examples/minimization/minhelper.py b/documentation/sphinx/examples/minimization/minhelper.py index 3339e40..aa03389 100644 --- a/documentation/sphinx/examples/minimization/minhelper.py +++ b/documentation/sphinx/examples/minimization/minhelper.py @@ -2,17 +2,24 @@ This is a helper module for the minimization examples. """ +#FIXME: the pyipopt stuff has been commented out for now; +#FIXME: update this code to use a better python ipopt interface when available +# https://github.com/xuy/pyipopt +# http://gitview.danfis.cz/pipopt +# https://bitbucket.org/amitibo +# https://github.com/casadi/casadi + import functools import numpy import scipy.optimize import algopy import numdifftools -import pyipopt +#import pyipopt # Suppress log spam from pyipopt. # But ipopt itself will stil spam... -pyipopt.set_loglevel(0) +#pyipopt.set_loglevel(0) def eval_grad(f, theta): @@ -176,30 +183,30 @@ def do_searches(f, g, h, x0): print results print - print 'strategy:', 'ipopt' - print 'options:', 'default' - print 'gradient:', 'autodiff' - print 'hessian:', 'autodiff' - results = pyipopt.fmin_unconstrained( - f, - x0, - fprime=g, - fhess=h, - ) - print results - print - - print 'strategy:', 'ipopt' - print 'options:', 'default' - print 'gradient:', 'autodiff' - print 'hessian:', 'finite differences' - results = pyipopt.fmin_unconstrained( - f, - x0, - fprime=g, - ) - print results - print + #print 'strategy:', 'ipopt' + #print 'options:', 'default' + #print 'gradient:', 'autodiff' + #print 'hessian:', 'autodiff' + #results = pyipopt.fmin_unconstrained( + #f, + #x0, + #fprime=g, + #fhess=h, + #) + #print results + #print + + #print 'strategy:', 'ipopt' + #print 'options:', 'default' + #print 'gradient:', 'autodiff' + #print 'hessian:', 'finite differences' + #results = pyipopt.fmin_unconstrained( + #f, + #x0, + #fprime=g, + #) + #print results + #print def show_minimization_results(f, target_in, easy_init_in, hard_init_in): From 9a1ca740099ddf5975e8fea39a4626b89c5c983d Mon Sep 17 00:00:00 2001 From: alex Date: Fri, 7 Dec 2012 09:02:44 -0500 Subject: [PATCH 5/6] Function.__truediv__ --- algopy/tracer/tracer.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/algopy/tracer/tracer.py b/algopy/tracer/tracer.py index b84651b..1880981 100644 --- a/algopy/tracer/tracer.py +++ b/algopy/tracer/tracer.py @@ -1017,6 +1017,9 @@ def __rdiv__(self, lhs): lhs = self.__class__.totype(lhs) return lhs/self + __truediv__ = __div__ + __rtruediv__ = __rdiv__ + # ######################################################### # numpy functions From 2f305637d39019bde337a4088ade95ef2810adab Mon Sep 17 00:00:00 2001 From: alex Date: Fri, 7 Dec 2012 20:02:58 -0500 Subject: [PATCH 6/6] rephrase a tracer.py loop so that it uses more than one line of code so that I can trace my error more accurately --- algopy/tracer/tracer.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/algopy/tracer/tracer.py b/algopy/tracer/tracer.py index 1880981..8128de5 100644 --- a/algopy/tracer/tracer.py +++ b/algopy/tracer/tracer.py @@ -246,7 +246,10 @@ def f(x1, x2): else: x_list = [x] - utpm_x_list = [algopy.UTPM(numpy.asarray(xi).reshape((1,1) + numpy.shape(xi))) for xi in x_list] + utpm_x_list = [] + for xi in x_list: + element = numpy.asarray(xi).reshape((1,1) + numpy.shape(xi)) + utpm_x_list.append(algopy.UTPM(element)) self.pushforward(utpm_x_list)