From 5a476db6d2a907b4de9a42a141936bdbf32df3b7 Mon Sep 17 00:00:00 2001
From: alex <argriffi@ncsu.edu>
Date: Wed, 28 Nov 2012 12:45:08 -0500
Subject: [PATCH 1/8] 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 <argriffi@ncsu.edu>
Date: Wed, 28 Nov 2012 13:10:45 -0500
Subject: [PATCH 2/8] 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 <argriffi@ncsu.edu>
Date: Wed, 28 Nov 2012 18:39:16 -0500
Subject: [PATCH 3/8] 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 <argriffi@ncsu.edu>
Date: Wed, 5 Dec 2012 10:07:23 -0500
Subject: [PATCH 4/8] 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 <argriffi@ncsu.edu>
Date: Fri, 7 Dec 2012 09:02:44 -0500
Subject: [PATCH 5/8] 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 <argriffi@ncsu.edu>
Date: Fri, 7 Dec 2012 20:02:58 -0500
Subject: [PATCH 6/8] 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)
 

From 37b06d82094718aa227d9eff6e001d09fde4301a Mon Sep 17 00:00:00 2001
From: alex <argriffi@ncsu.edu>
Date: Sat, 15 Dec 2012 17:26:42 -0500
Subject: [PATCH 7/8] add a logistic regression example

---
 .../sphinx/examples/logistic_regression.py    | 307 ++++++++++++++++++
 .../sphinx/examples/logistic_regression.rst   |  39 +++
 documentation/sphinx/index.rst                |   1 +
 3 files changed, 347 insertions(+)
 create mode 100644 documentation/sphinx/examples/logistic_regression.py
 create mode 100644 documentation/sphinx/examples/logistic_regression.rst

diff --git a/documentation/sphinx/examples/logistic_regression.py b/documentation/sphinx/examples/logistic_regression.py
new file mode 100644
index 0000000..5aed109
--- /dev/null
+++ b/documentation/sphinx/examples/logistic_regression.py
@@ -0,0 +1,307 @@
+"""
+Logistic Regression
+
+This is an algopy implementation of a scicomp stackexchange question
+http://scicomp.stackexchange.com/questions/4826/logistic-regression-with-python
+"""
+
+import functools
+
+import numpy
+import scipy.optimize
+import algopy
+
+
+g_data = """\
+low,age,lwt,race,smoke,ptl,ht,ui
+0,19,182,2,0,0,0,1
+0,33,155,3,0,0,0,0
+0,20,105,1,1,0,0,0
+0,21,108,1,1,0,0,1
+0,18,107,1,1,0,0,1
+0,21,124,3,0,0,0,0
+0,22,118,1,0,0,0,0
+0,17,103,3,0,0,0,0
+0,29,123,1,1,0,0,0
+0,26,113,1,1,0,0,0
+0,19,95,3,0,0,0,0
+0,19,150,3,0,0,0,0
+0,22,95,3,0,0,1,0
+0,30,107,3,0,1,0,1
+0,18,100,1,1,0,0,0
+0,18,100,1,1,0,0,0
+0,15,98,2,0,0,0,0
+0,25,118,1,1,0,0,0
+0,20,120,3,0,0,0,1
+0,28,120,1,1,0,0,0
+0,32,121,3,0,0,0,0
+0,31,100,1,0,0,0,1
+0,36,202,1,0,0,0,0
+0,28,120,3,0,0,0,0
+0,25,120,3,0,0,0,1
+0,28,167,1,0,0,0,0
+0,17,122,1,1,0,0,0
+0,29,150,1,0,0,0,0
+0,26,168,2,1,0,0,0
+0,17,113,2,0,0,0,0
+0,17,113,2,0,0,0,0
+0,24,90,1,1,1,0,0
+0,35,121,2,1,1,0,0
+0,25,155,1,0,0,0,0
+0,25,125,2,0,0,0,0
+0,29,140,1,1,0,0,0
+0,19,138,1,1,0,0,0
+0,27,124,1,1,0,0,0
+0,31,215,1,1,0,0,0
+0,33,109,1,1,0,0,0
+0,21,185,2,1,0,0,0
+0,19,189,1,0,0,0,0
+0,23,130,2,0,0,0,0
+0,21,160,1,0,0,0,0
+0,18,90,1,1,0,0,1
+0,18,90,1,1,0,0,1
+0,32,132,1,0,0,0,0
+0,19,132,3,0,0,0,0
+0,24,115,1,0,0,0,0
+0,22,85,3,1,0,0,0
+0,22,120,1,0,0,1,0
+0,23,128,3,0,0,0,0
+0,22,130,1,1,0,0,0
+0,30,95,1,1,0,0,0
+0,19,115,3,0,0,0,0
+0,16,110,3,0,0,0,0
+0,21,110,3,1,0,0,1
+0,30,153,3,0,0,0,0
+0,20,103,3,0,0,0,0
+0,17,119,3,0,0,0,0
+0,17,119,3,0,0,0,0
+0,23,119,3,0,0,0,0
+0,24,110,3,0,0,0,0
+0,28,140,1,0,0,0,0
+0,26,133,3,1,2,0,0
+0,20,169,3,0,1,0,1
+0,24,115,3,0,0,0,0
+0,28,250,3,1,0,0,0
+0,20,141,1,0,2,0,1
+0,22,158,2,0,1,0,0
+0,22,112,1,1,2,0,0
+0,31,150,3,1,0,0,0
+0,23,115,3,1,0,0,0
+0,16,112,2,0,0,0,0
+0,16,135,1,1,0,0,0
+0,18,229,2,0,0,0,0
+0,25,140,1,0,0,0,0
+0,32,134,1,1,1,0,0
+0,20,121,2,1,0,0,0
+0,23,190,1,0,0,0,0
+0,22,131,1,0,0,0,0
+0,32,170,1,0,0,0,0
+0,30,110,3,0,0,0,0
+0,20,127,3,0,0,0,0
+0,23,123,3,0,0,0,0
+0,17,120,3,1,0,0,0
+0,19,105,3,0,0,0,0
+0,23,130,1,0,0,0,0
+0,36,175,1,0,0,0,0
+0,22,125,1,0,0,0,0
+0,24,133,1,0,0,0,0
+0,21,134,3,0,0,0,0
+0,19,235,1,1,0,1,0
+0,25,95,1,1,3,0,1
+0,16,135,1,1,0,0,0
+0,29,135,1,0,0,0,0
+0,29,154,1,0,0,0,0
+0,19,147,1,1,0,0,0
+0,19,147,1,1,0,0,0
+0,30,137,1,0,0,0,0
+0,24,110,1,0,0,0,0
+0,19,184,1,1,0,1,0
+0,24,110,3,0,1,0,0
+0,23,110,1,0,0,0,0
+0,20,120,3,0,0,0,0
+0,25,241,2,0,0,1,0
+0,30,112,1,0,0,0,0
+0,22,169,1,0,0,0,0
+0,18,120,1,1,0,0,0
+0,16,170,2,0,0,0,0
+0,32,186,1,0,0,0,0
+0,18,120,3,0,0,0,0
+0,29,130,1,1,0,0,0
+0,33,117,1,0,0,0,1
+0,20,170,1,1,0,0,0
+0,28,134,3,0,0,0,0
+0,14,135,1,0,0,0,0
+0,28,130,3,0,0,0,0
+0,25,120,1,0,0,0,0
+0,16,95,3,0,0,0,0
+0,20,158,1,0,0,0,0
+0,26,160,3,0,0,0,0
+0,21,115,1,0,0,0,0
+0,22,129,1,0,0,0,0
+0,25,130,1,0,0,0,0
+0,31,120,1,0,0,0,0
+0,35,170,1,0,1,0,0
+0,19,120,1,1,0,0,0
+0,24,116,1,0,0,0,0
+0,45,123,1,0,0,0,0
+1,28,120,3,1,1,0,1
+1,29,130,1,0,0,0,1
+1,34,187,2,1,0,1,0
+1,25,105,3,0,1,1,0
+1,25,85,3,0,0,0,1
+1,27,150,3,0,0,0,0
+1,23,97,3,0,0,0,1
+1,24,128,2,0,1,0,0
+1,24,132,3,0,0,1,0
+1,21,165,1,1,0,1,0
+1,32,105,1,1,0,0,0
+1,19,91,1,1,2,0,1
+1,25,115,3,0,0,0,0
+1,16,130,3,0,0,0,0
+1,25,92,1,1,0,0,0
+1,20,150,1,1,0,0,0
+1,21,200,2,0,0,0,1
+1,24,155,1,1,1,0,0
+1,21,103,3,0,0,0,0
+1,20,125,3,0,0,0,1
+1,25,89,3,0,2,0,0
+1,19,102,1,0,0,0,0
+1,19,112,1,1,0,0,1
+1,26,117,1,1,1,0,0
+1,24,138,1,0,0,0,0
+1,17,130,3,1,1,0,1
+1,20,120,2,1,0,0,0
+1,22,130,1,1,1,0,1
+1,27,130,2,0,0,0,1
+1,20,80,3,1,0,0,1
+1,17,110,1,1,0,0,0
+1,25,105,3,0,1,0,0
+1,20,109,3,0,0,0,0
+1,18,148,3,0,0,0,0
+1,18,110,2,1,1,0,0
+1,20,121,1,1,1,0,1
+1,21,100,3,0,1,0,0
+1,26,96,3,0,0,0,0
+1,31,102,1,1,1,0,0
+1,15,110,1,0,0,0,0
+1,23,187,2,1,0,0,0
+1,20,122,2,1,0,0,0
+1,24,105,2,1,0,0,0
+1,15,115,3,0,0,0,1
+1,23,120,3,0,0,0,0
+1,30,142,1,1,1,0,0
+1,22,130,1,1,0,0,0
+1,17,120,1,1,0,0,0
+1,23,110,1,1,1,0,0
+1,17,120,2,0,0,0,0
+1,26,154,3,0,1,1,0
+1,20,106,3,0,0,0,0
+1,26,190,1,1,0,0,0
+1,14,101,3,1,1,0,0
+1,28,95,1,1,0,0,0
+1,14,100,3,0,0,0,0
+1,23,94,3,1,0,0,0
+1,17,142,2,0,0,1,0
+1,21,130,1,1,0,1,0
+"""
+
+
+########################################################################
+# boilerplate functions for algopy
+
+def eval_grad(f, theta):
+    theta = algopy.UTPM.init_jacobian(theta)
+    return algopy.UTPM.extract_jacobian(f(theta))
+
+def eval_hess(f, theta):
+    theta = algopy.UTPM.init_hessian(theta)
+    return algopy.UTPM.extract_hessian(len(theta), f(theta))
+
+
+########################################################################
+# application specific functions
+
+def get_neg_ll(vY, mX, vBeta):
+    """
+    @param vY: predefined numpy array
+    @param mX: predefined numpy array
+    @param vBeta: parameters of the likelihood function
+    """
+    #FIXME: algopy could benefit from the addition of a logsumexp function...
+    alpha = algopy.dot(mX, vBeta)
+    return -algopy.sum(
+            vY*algopy.log(algopy.special.expit(alpha)) +
+            (1-vY)*(algopy.log(algopy.special.expit(-alpha))))
+
+def preprocess_data():
+    """
+    Convert the data from a hardcoded string into something nicer.
+    """
+    data = numpy.loadtxt(g_data.splitlines(), delimiter=',', skiprows=1)
+    vY = data[:, 0]
+    mX = data[:, 1:]
+    intercept = numpy.ones(mX.shape[0]).reshape(mX.shape[0], 1)
+    mX = numpy.concatenate((intercept, mX), axis=1)
+    iK = mX.shape[1]
+    iN = mX.shape[0]
+    return vY, mX
+
+def main():
+
+    # extract the data from the hardcoded string
+    vY, mX = preprocess_data()
+
+    # this is a hardcoded point which is supposed to have a good likelihood
+    vBeta_star = numpy.array([
+        -.10296645, -.0332327, -.01209484, .44626211,
+        .92554137, .53973828, 1.7993371, .7148045,
+        ])
+
+    # these are arbitrary parameter values somwhat near the good values
+    vBeta_0 = numpy.array([-.1, -.03, -.01, .44, .92, .53, 1.8, .71])
+
+    # define the objective function and the autodiff gradient and hessian
+    f = functools.partial(get_neg_ll, vY, mX)
+    g = functools.partial(eval_grad, f)
+    h = functools.partial(eval_hess, f)
+
+    # show the neg log likelihood for the good parameters
+    print 'hardcoded good values:'
+    print vBeta_star
+    print
+    print 'neg log likelihood for good values:'
+    print f(vBeta_star)
+    print
+    print
+    print 'hardcoded okay values:'
+    print vBeta_0
+    print
+    print 'neg log likelihood for okay values:'
+    print f(vBeta_0)
+    print
+    print
+
+    # do the max likelihood search
+    results = scipy.optimize.fmin_ncg(
+            f,
+            vBeta_0,
+            fprime=g,
+            fhess=h,
+            avextol=1e-6,
+            disp=0,
+            )
+
+    # extract the max likelihood values
+    vBeta_mle = results
+
+    print 'maximum likelihood estimates:'
+    print vBeta_mle
+    print
+    print 'neg log likelihood for maximum likelihood estimates:'
+    print f(vBeta_mle)
+    print
+
+
+if __name__ == '__main__':
+    main()
+
diff --git a/documentation/sphinx/examples/logistic_regression.rst b/documentation/sphinx/examples/logistic_regression.rst
new file mode 100644
index 0000000..83ac49a
--- /dev/null
+++ b/documentation/sphinx/examples/logistic_regression.rst
@@ -0,0 +1,39 @@
+Logistic Regression
+----------------------------
+
+In this example we want to use AlgoPy to help compute the
+maximum likelihood estimates for a nonlinear model.
+It is based on this
+`question <http://scicomp.stackexchange.com/questions/4826>`_
+on the scicomp stackexchange.
+
+
+Here is the python code:
+
+.. literalinclude:: logistic_regression.py
+
+
+Here is its output::
+
+    hardcoded good values:
+    [-0.10296645 -0.0332327  -0.01209484  0.44626211  0.92554137  0.53973828
+      1.7993371   0.7148045 ]
+
+    neg log likelihood for good values:
+    102.173732637
+
+
+    hardcoded okay values:
+    [-0.1  -0.03 -0.01  0.44  0.92  0.53  1.8   0.71]
+
+    neg log likelihood for okay values:
+    104.084160515
+
+
+    maximum likelihood estimates:
+    [-0.10296656 -0.0332327  -0.01209484  0.44626209  0.92554132  0.53973824
+      1.79933692  0.71480445]
+
+    neg log likelihood for maximum likelihood estimates:
+    102.173732637
+
diff --git a/documentation/sphinx/index.rst b/documentation/sphinx/index.rst
index 3d39bb1..4d0cfe9 100644
--- a/documentation/sphinx/index.rst
+++ b/documentation/sphinx/index.rst
@@ -337,6 +337,7 @@ Application Examples:
    examples/hessian_of_potential_function.rst
    examples/minimal_surface.rst
    examples/neg_binom_regression.rst
+   examples/logistic_regression.rst
 
 
 Additional Information:

From 8e78e919466551001f2ecac3564e1f7a36e8d4aa Mon Sep 17 00:00:00 2001
From: alex <argriffi@ncsu.edu>
Date: Sat, 15 Dec 2012 17:44:26 -0500
Subject: [PATCH 8/8] reformulate the log likelihood for the logistic
 regression example

---
 documentation/sphinx/examples/logistic_regression.py  | 6 +++---
 documentation/sphinx/examples/logistic_regression.rst | 4 ++--
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/documentation/sphinx/examples/logistic_regression.py b/documentation/sphinx/examples/logistic_regression.py
index 5aed109..40dbe4e 100644
--- a/documentation/sphinx/examples/logistic_regression.py
+++ b/documentation/sphinx/examples/logistic_regression.py
@@ -229,9 +229,9 @@ def get_neg_ll(vY, mX, vBeta):
     """
     #FIXME: algopy could benefit from the addition of a logsumexp function...
     alpha = algopy.dot(mX, vBeta)
-    return -algopy.sum(
-            vY*algopy.log(algopy.special.expit(alpha)) +
-            (1-vY)*(algopy.log(algopy.special.expit(-alpha))))
+    return algopy.sum(
+            vY*algopy.log1p(algopy.exp(-alpha)) +
+            (1-vY)*algopy.log1p(algopy.exp(alpha)))
 
 def preprocess_data():
     """
diff --git a/documentation/sphinx/examples/logistic_regression.rst b/documentation/sphinx/examples/logistic_regression.rst
index 83ac49a..e56c64b 100644
--- a/documentation/sphinx/examples/logistic_regression.rst
+++ b/documentation/sphinx/examples/logistic_regression.rst
@@ -31,8 +31,8 @@ Here is its output::
 
 
     maximum likelihood estimates:
-    [-0.10296656 -0.0332327  -0.01209484  0.44626209  0.92554132  0.53973824
-      1.79933692  0.71480445]
+    [-0.10296655 -0.0332327  -0.01209484  0.44626209  0.92554133  0.53973824
+      1.79933696  0.71480445]
 
     neg log likelihood for maximum likelihood estimates:
     102.173732637