Merge pull request #99 from convexengineering/issue_95

v0.1 refactor

LICENSE

@@ -1,6 +1,6 @@
 The MIT License (MIT)
 
-Copyright (c) 2014 MIT Applied Optimization Group
+Copyright (c) 2021 Convex Engineering Group
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal

README.md

@@ -1,8 +1,10 @@
 <img src="http://gpfit.readthedocs.io/en/latest/_images/GPfit_logo.png" width=110 alt="GPfit" />
 
+**[Documentation](http://gpfit.readthedocs.org/)** | [Examples](http://gpfit.readthedocs.org/en/latest/examples.html) |
+
 GPfit is a Python package for fitting geometric programming models to data.
-It requires installation of [GPkit](http://gpkit.readthedocs.org/en/latest/).
-This [paper](http://hoburg.mit.edu/publications/gpfitting.pdf)
-describes the approach.
+This [paper](https://convex.mit.edu/publications/gpfitting.pdf) describes the
+approach.
+GPfit requires installation of [GPkit](http://gpkit.readthedocs.org/en/latest/).
 
 [![Build Status](https://acdl.mit.edu/csi/buildStatus/icon?job=CE_gpfit_Push)](https://acdl.mit.edu/csi/job/CE_gpfit_Push/)

docs/source/citinggpfit.rst

@@ -3,10 +3,13 @@ Citing GPfit
 
 If you use GPfit, please cite it with the following bibtex::
 
-    @Misc{gpfit,
-          author={MIT Convex Optimization Group},
-          title={GPfit},
-          howpublished={\url{https://github.com/convexopt/gpfit}},
-          year={2015},
-          note={Version 0.1}
-         }
+        @article{hoburg2016data,
+          title={Data fitting with geometric-programming-compatible softmax functions},
+          author={Hoburg, Warren and Kirschen, Philippe and Abbeel, Pieter},
+          journal={Optimization and Engineering},
+          volume={17},
+          number={4},
+          pages={897--918},
+          year={2016},
+          publisher={Springer}
+        }

docs/source/conf.py

@@ -49,7 +49,7 @@
 
 # General information about the project.
 project = 'GPfit'
-copyright = '2015, MIT Convex Optimization Group'
+copyright = '2021, Convex Engineering Group'
 
 # The version info for the project you're documenting, acts as replacement for
 # |version| and |release|, also used in various other places throughout the

docs/source/examples.rst

@@ -1,14 +1,25 @@
 Examples
 ********
 
+Both examples come from Section 6 of the paper "`Fitting geometric programming
+models to data <http://web.mit.edu/~whoburg/www/papers/gp_fitting.pdf>`_", by
+W. Hoburg et. al.
 
 Example 1
 =========
-This example comes from Section 6 of the paper "`Fitting geometric programming models to data <http://web.mit.edu/~whoburg/www/papers/gp_fitting.pdf>`_", by W. Hoburg et. al.
 Fit convex portion of :math:`w = \frac{u^2 + 3}{(u+1)^2}` on :math:`1 \leq u \leq 3`.
 
 .. literalinclude:: examples/hoburgabbeel_ex6_1.py
 
 Output:
 
 .. literalinclude:: examples/hoburgabbeel_ex6_1_output.txt
+
+Example 2
+=========
+
+.. literalinclude:: examples/hoburgabbeel_ex6_3.py
+
+Output:
+
+.. literalinclude:: examples/hoburgabbeel_ex6_3_output.txt

docs/source/examples/hoburgabbeel_ex6_1.py

@@ -1,13 +1,11 @@
 "Fits an example function"
-from __future__ import division
 from numpy import logspace, log, log10, random
 from gpfit.fit import fit
 
-# fixed initial guess for fitting
 random.seed(33404)
 
 u = logspace(0, log10(3), 101)
-w = (u**2 + 3) / (u + 1)**2
+w = (u**2 + 3)/(u + 1)**2
 x = log(u)
 y = log(w)
 K = 3

docs/source/examples/hoburgabbeel_ex6_1_output.txt

@@ -7,9 +7,9 @@ w**3.4411 = 0.422736 * (u_1)**-2.14843
     + 0.424169 * (u_1)**-2.14784
     + 0.15339 * (u_1)**0.584654
 ISMA fit from params
-1 = (0.994797/w**0.238961) * (u_1)**-0.138389
-    + (0.949519/w**0.0924504) * (u_1)**0.0165798
-    + (0.967646/w**0.115505) * (u_1)**-0.0131876
+1 = (0.992648/w**0.35353) * (u_1)**-0.204093
+    + (0.947302/w**0.0920266) * (u_1)**0.017725
+    + (0.961409/w**0.11673) * (u_1)**-0.011164
 MA RMS Error: 0.0023556
 SMA RMS Error: 2.3856e-05
-ISMA RMS Error: 1.0765e-06
+ISMA RMS Error: 8.0757e-07

docs/source/examples/hoburgabbeel_ex6_3.py

@@ -0,0 +1,22 @@
+"Example 6.3 from Hoburg/Abbeel GPfit paper"
+import numpy as np
+from numpy.random import seed, random_sample
+from gpfit.fit import fit
+
+seed(33404)
+
+Vdd = random_sample(1000) + 1
+Vth = 0.2*random_sample(1000) + 0.2
+P = Vdd**2 + 30*Vdd*np.exp(-(Vth - 0.06*Vdd)/0.039)
+u = np.vstack((Vdd, Vth))
+x = np.log(u)
+y = np.log(P)
+K = 4
+
+_, errorMA = fit(x, y, K, "MA")
+_, errorSMA = fit(x, y, K, "SMA")
+_, errorISMA = fit(x, y, K, "ISMA")
+
+print("MA RMS Error: %.5g" % errorMA)
+print("SMA RMS Error: %.5g" % errorSMA)
+print("ISMA RMS Error: %.5g" % errorISMA)

docs/source/examples/hoburgabbeel_ex6_3_output.txt

@@ -0,0 +1,18 @@
+MA fit from params
+w = 0.0877784 * (u_1)**2.3499 * (u_2)**-1.86683
+w = 0.824389 * (u_1)**2.02321 * (u_2)**-0.203195
+w = 0.0175649 * (u_1)**2.84017 * (u_2)**-2.78935
+w = 0.38541 * (u_1)**2.13119 * (u_2)**-0.826143
+SMA fit from params
+w**2.22467 = 5.42035e-05 * (u_1)**5.24451 * (u_2)**-6.68343
+    + 3.61239e-08 * (u_1)**9.613 * (u_2)**-9.4893
+    + 1.07674e-05 * (u_1)**0.538839 * (u_2)**-6.12098
+    + 1.0777 * (u_1)**4.43371 * (u_2)**0.0978425
+ISMA fit from params
+1 = (0.309986/w**0.309511) * (u_1)**0.6327 * (u_2)**-0.813931
+    + (0.012176/w**0.67353) * (u_1)**2.61311 * (u_2)**-2.59033
+    + (0.889774/w**0.242028) * (u_1)**0.485343 * (u_2)**-0.0942005
+    + (1.00588/w**0.0753052) * (u_1)**0.149934 * (u_2)**0.0116639
+MA RMS Error: 0.007382
+SMA RMS Error: 0.00016261
+ISMA RMS Error: 0.00017923

docs/source/index.rst

@@ -12,7 +12,6 @@ Table of Contents
        :maxdepth: 2
 
        gpfit101
-       gettingstarted
        examples
        citinggpfit
        glossary.rst

gpfit/classes.py

@@ -0,0 +1,135 @@
+"""
+Implements the Max Affine, Softmax Affine, and Implicit Softmax Affine function
+classes
+"""
+import numpy as np
+from .logsumexp import lse_scaled, lse_implicit
+
+
+def max_affine(x, ba):
+    """
+    Evaluates max affine function at values of x, given a set of
+    max affine fit parameters.
+
+    Arguments
+    ---------
+        x: 2D array [nPoints x nDim]
+            Independent variable data
+
+        ba: 2D array
+            max affine fit parameters
+            [[b1, a11, ... a1k]
+             [ ....,          ]
+             [bk, ak1, ... akk]]
+
+    Returns
+    -------
+        y: 1D array [nPoints]
+            Max affine output
+        dydba: 2D array [nPoints x (nDim + 1)*K]
+            dydba
+    """
+    npt, dimx = x.shape
+    K = ba.size // (dimx + 1)
+    ba = np.reshape(ba, (dimx + 1, K), order="F")  # 'F' gives Fortran indexing
+    X = np.hstack((np.ones((npt, 1)), x))  # augment data with column of ones
+    y, partition = np.dot(X, ba).max(1), np.dot(X, ba).argmax(1)
+
+    dydba = np.zeros((npt, (dimx + 1)*K))
+    for k in range(K):
+        inds = np.equal(partition, k)
+        indadd = (dimx + 1)*k
+        ixgrid = np.ix_(inds.nonzero()[0], indadd + np.arange(dimx + 1))
+        dydba[ixgrid] = X[inds, :]
+
+    return y, dydba
+
+
+# pylint: disable=too-many-locals
+def softmax_affine(x, params):
+    """
+    Evaluates softmax affine function at values of x, given a set of
+    SMA fit parameters.
+
+    Arguments:
+    ----------
+            x:      Independent variable data
+                        2D numpy array [nPoints x nDimensions]
+
+            params: Fit parameters
+                        1D numpy array [(nDim + 2)*K,]
+                        [b1, a11, .. a1d, b2, a21, .. a2d, ...
+                         bK, aK1, aK2, .. aKd, alpha]
+
+    Returns:
+    --------
+            y:      SMA approximation to log transformed data
+                        1D numpy array [nPoints]
+
+            dydp:   Jacobian matrix
+    """
+
+    npt, dimx = x.shape
+    ba = params[0:-1]
+    softness = params[-1]
+    alpha = 1/softness
+    if alpha <= 0:
+        return np.inf*np.ones((npt, 1)), np.nan
+    K = np.size(ba) // (dimx + 1)
+    ba = ba.reshape(dimx + 1, K, order="F")
+    X = np.hstack((np.ones((npt, 1)), x))  # augment data with column of ones
+    z = np.dot(X, ba)  # compute affine functions
+    y, dydz, dydsoftness = lse_scaled(z, alpha)
+
+    dydsoftness = -dydsoftness*(alpha**2)
+    nrow, ncol = dydz.shape
+    repmat = np.tile(dydz, (dimx + 1, 1)).reshape(nrow, ncol*(dimx + 1), order="F")
+    dydba = repmat*np.tile(X, (1, K))
+    dydsoftness.shape = (dydsoftness.size, 1)
+    dydp = np.hstack((dydba, dydsoftness))
+
+    return y, dydp
+
+
+# pylint: disable=too-many-locals
+def implicit_softmax_affine(x, params):
+    """
+    Evaluates implicit softmax affine function at values of x, given a set of
+    ISMA fit parameters.
+
+    Arguments:
+    ----------
+            x:      Independent variable data
+                        2D numpy array [nPoints x nDimensions]
+
+            params: Fit parameters
+                        1D numpy array [(nDim + 2)*K,]
+                        [b1, a11, .. a1d, b2, a21, .. a2d, ...
+                         bK, aK1, aK2, .. aKd, alpha1, alpha2, ... alphaK]
+
+    Returns:
+    --------
+            y:      ISMA approximation to log transformed data
+                        1D numpy array [nPoints]
+
+            dydp:   Jacobian matrix
+
+    """
+
+    npt, dimx = x.shape
+    K = params.size // (dimx + 2)
+    ba = params[0:-K]
+    alpha = params[-K:]
+    if any(alpha <= 0):
+        return np.inf*np.ones((npt, 1)), np.nan
+    ba = ba.reshape(dimx + 1, K, order="F")  # reshape ba to matrix
+    X = np.hstack((np.ones((npt, 1)), x))  # augment data with column of ones
+    z = np.dot(X, ba)  # compute affine functions
+    y, dydz, dydalpha = lse_implicit(z, alpha)
+
+    nrow, ncol = dydz.shape
+    repmat = np.tile(dydz, (dimx + 1, 1)).reshape(nrow, ncol*(dimx + 1), order="F")
+    dydba = repmat*np.tile(X, (1, K))
+    dydp = np.hstack((dydba, dydalpha))
+
+    return y, dydp