Repository: convexengineering/gpfit (language: python)

Issue:
v0.1 refactor
I propose we create a 0.0 release of GPfit and then perform a pretty wholesale refactor of the code. Here is the list of changes I would like to make:

- [x] General code tightening (more concise code, less unnecessary whitespace)
- [x] Py3 only (e.g. delete `__future__` imports)
- [x] Delete legacy Matlab code
- [x] Better variable naming
    - [x] Appropriate use of case (rename upper case and mixed case vars unless global)
- [x] Better function naming
    - [x] Rename `get_params`
    - [x] Rename `ba_init`
- [x] Better class naming
    - [x] Rename `FitCS`
- [x] Code restructuring
    - [x] MA, SMA, ISMA should all be in one module
    - [x] Reduce confusing overlap of naming of `fit.py` and `fit_constraintset.py`
    - [x] Separate FitCS and XfoilFit (move to `xfoil.py`)
    - [x] Group `lse_scaled` and `lse_implicit` into one module
    - ~Fit should be a class not a function~
    - ~Eliminate `plot_fit.py` and `print_fit.py` (they should be class methods)~
- [x] Remove examples/tests redundancy
    - [x] Better naming of the examples too
- #98
- ~Write the Getting Started docs page~

There is a little overlap with #73 but would be good to include those items too!

@bqpd @whoburg any objections or things you want to add?


Discussion / hints:
Those look good to me! I might just add "write the Getting Started docs page"?
Good call. The docs need a lot of love but that seems like an appropriate piece to include in this effort.
The cleaner the API the shorter that page will be, too!

@pgkirsch this sounds fantastic! What are your thoughts on using a code formatter like `black` to accomplish some of the earlier bullets?
Thanks @whoburg! `black` seems cool, I'll give it a go.

Fix patch:
diff --git a/LICENSE b/LICENSE
index 375fa23..af0974d 100644
--- a/LICENSE
+++ b/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
diff --git a/README.md b/README.md
index af6d726..b7f87f8 100644
--- a/README.md
+++ b/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/)
diff --git a/docs/source/citinggpfit.rst b/docs/source/citinggpfit.rst
index 446fde5..a982232 100644
--- a/docs/source/citinggpfit.rst
+++ b/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}
-         }
\ No newline at end of file
+        @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}
+        }
diff --git a/docs/source/conf.py b/docs/source/conf.py
index 915d7f5..b42595c 100644
--- a/docs/source/conf.py
+++ b/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
diff --git a/docs/source/examples.rst b/docs/source/examples.rst
index 23f1600..5870ac9 100644
--- a/docs/source/examples.rst
+++ b/docs/source/examples.rst
@@ -1,10 +1,12 @@
 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
@@ -12,3 +14,12 @@ Fit convex portion of :math:`w = \frac{u^2 + 3}{(u+1)^2}` on :math:`1 \leq u \le
 Output:
 
 .. literalinclude:: examples/hoburgabbeel_ex6_1_output.txt
+
+Example 2
+=========
+
+.. literalinclude:: examples/hoburgabbeel_ex6_3.py
+
+Output:
+
+.. literalinclude:: examples/hoburgabbeel_ex6_3_output.txt
diff --git a/docs/source/examples/hoburgabbeel_ex6_1.py b/docs/source/examples/hoburgabbeel_ex6_1.py
index 0df7d66..749469d 100644
--- a/docs/source/examples/hoburgabbeel_ex6_1.py
+++ b/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
diff --git a/docs/source/examples/hoburgabbeel_ex6_1_output.txt b/docs/source/examples/hoburgabbeel_ex6_1_output.txt
index 5443bc8..d0bbd7c 100644
--- a/docs/source/examples/hoburgabbeel_ex6_1_output.txt
+++ b/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
diff --git a/docs/source/examples/hoburgabbeel_ex6_3.py b/docs/source/examples/hoburgabbeel_ex6_3.py
new file mode 100644
index 0000000..783cbc8
--- /dev/null
+++ b/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)
diff --git a/docs/source/examples/hoburgabbeel_ex6_3_output.txt b/docs/source/examples/hoburgabbeel_ex6_3_output.txt
new file mode 100644
index 0000000..beeeb39
--- /dev/null
+++ b/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
diff --git a/docs/source/gettingstarted.rst b/docs/source/gettingstarted.rst
deleted file mode 100644
index 671023d..0000000
--- a/docs/source/gettingstarted.rst
+++ /dev/null
@@ -1,2 +0,0 @@
-Getting Started
-***************
\ No newline at end of file
diff --git a/docs/source/index.rst b/docs/source/index.rst
index fd81956..6ec95b0 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -12,7 +12,6 @@ Table of Contents
        :maxdepth: 2
 
        gpfit101
-       gettingstarted
        examples
        citinggpfit
        glossary.rst
diff --git a/gpfit/classes.py b/gpfit/classes.py
new file mode 100644
index 0000000..8c761ea
--- /dev/null
+++ b/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
diff --git a/gpfit/constraint_set.py b/gpfit/constraint_set.py
new file mode 100644
index 0000000..932f90a
--- /dev/null
+++ b/gpfit/constraint_set.py
@@ -0,0 +1,134 @@
+"""Fit constraint set"""
+from numpy import amax, array, hstack
+from gpkit import ConstraintSet
+from gpkit import Variable, NomialArray, NamedVariables, VectorVariable
+from gpkit.small_scripts import initsolwarning, appendsolwarning
+
+
+# pylint: disable=too-many-instance-attributes, too-many-locals,
+# pylint: disable=too-many-branches, no-member, import-error
+# pylint: disable=too-many-arguments, import-outside-toplevel
+class FitConstraintSet(ConstraintSet):
+    """
+    Constraint set for fitted functions
+
+    Arguments
+    ---------
+    fitdata : dict
+        dictionary of fit parameters
+    ivar : gpkit Variable, Monomial, or NomialArray
+        independent variable
+    dvars : list of gpkit Variables, Monomials, or NomialArrays
+        dependent variables
+    err_margin : string, either "max" or "rms"
+        flag to set margin factor using RMS or max error
+
+    """
+
+    def __init__(self, fitdata, ivar=None, dvars=None, name="fit", err_margin=None):
+
+        self.fitdata = fitdata
+
+        if ivar is None:
+            with NamedVariables("fit"):
+                dvars = VectorVariable(fitdata["d"], "u")
+                ivar = Variable("w")
+
+        self.dvars = dvars
+        self.ivar = ivar
+        self.rms_err = fitdata["rms_err"]
+        self.max_err = fitdata["max_err"]
+
+        monos = [
+            fitdata["c%d" % k]
+            *NomialArray(
+                array(dvars).T
+                **array([fitdata["e%d%d" % (k, i)] for i in range(fitdata["d"])])
+            ).prod(NomialArray(dvars).ndim - 1)
+            for k in range(fitdata["K"])
+        ]
+
+        with NamedVariables(name):
+            self.mfac = Variable("mfac", 1.0, "-", "fit factor")
+        if err_margin == "max":
+            self.mfac.key.descr["value"] = 1 + self.max_err
+        elif err_margin == "rms":
+            self.mfac.key.descr["value"] = 1 + self.rms_err
+
+        if fitdata["ftype"] == "ISMA":
+            # constraint of the form 1 >= c1*u1^exp1*u2^exp2*w^(-alpha) + ....
+            alpha = array([fitdata["a%d" % k] for k in range(fitdata["K"])])
+            lhs, rhs = 1, NomialArray(monos/(ivar/self.mfac)**alpha).sum(0)
+        elif fitdata["ftype"] == "SMA":
+            # constraint of the form w^alpha >= c1*u1^exp1 + c2*u2^exp2 +....
+            alpha = fitdata["a1"]
+            lhs, rhs = (ivar/self.mfac)**alpha, NomialArray(monos).sum(0)
+        elif fitdata["ftype"] == "MA":
+            # constraint of the form w >= c1*u1^exp1, w >= c2*u2^exp2, ....
+            lhs, rhs = (ivar/self.mfac), NomialArray(monos).T
+
+        if fitdata["K"] == 1:
+            # when possible, return an equality constraint
+            if hasattr(rhs, "shape"):
+                if rhs.ndim > 1:
+                    self.constraint = [(lh == rh) for lh, rh in zip(lhs, rhs)]
+                else:
+                    self.constraint = lhs == rhs
+            else:
+                self.constraint = lhs == rhs
+        else:
+            if hasattr(rhs, "shape"):
+                if rhs.ndim > 1:
+                    self.constraint = [(lh >= rh) for lh, rh in zip(lhs, rhs)]
+                else:
+                    self.constraint = lhs >= rhs
+            else:
+                self.constraint = lhs >= rhs
+
+        self.bounds = {}
+        for i, dvar in enumerate(self.dvars):
+            self.bounds[dvar] = [fitdata["lb%d" % i], fitdata["ub%d" % i]]
+
+        ConstraintSet.__init__(self, [self.constraint])
+
+    def get_dataframe(self):
+        "return a pandas DataFrame of fit parameters"
+        import pandas as pd
+
+        df = pd.DataFrame(list(self.fitdata.values())).transpose()
+        df.columns = list(self.fitdata.keys())
+        return df
+
+    def process_result(self, result):
+        """
+        make sure fit result is within bounds of fitted data
+        """
+        super().process_result(result)
+        initsolwarning(result, "Fit Out-of-Bounds")
+
+        if self.mfac not in result["sensitivities"]["constants"]:
+            return
+        if amax([abs(result["sensitivities"]["constants"][self.mfac])]) < 1e-5:
+            return
+
+        for dvar in self.dvars:
+            if isinstance(dvar, NomialArray):
+                num = [result(x) for x in dvar]
+            else:
+                num = result(dvar)
+            direct = None
+            if any(x < self.bounds[dvar][0] for x in hstack([num])):
+                direct, state = "lower", "below"
+                bnd = self.bounds[dvar][0]
+            if any(x > self.bounds[dvar][1] for x in hstack([num])):
+                direct, state = "upper", "above"
+                bnd = self.bounds[dvar][1]
+
+            if direct:
+                msg = (
+                    "Variable %.100s could cause inaccurate result"
+                    " because it is %s" % (dvar, state)
+                    + " %s bound. Solution is %.4f but"
+                    " bound is %.4f" % (direct, amax([num]), bnd)
+                )
+                appendsolwarning(msg, self, result, "Fit Out-of-Bounds")
diff --git a/gpfit/examples/ex61setup.py b/gpfit/examples/ex61setup.py
deleted file mode 100644
index 5ba96cf..0000000
--- a/gpfit/examples/ex61setup.py
+++ /dev/null
@@ -1,13 +0,0 @@
-from __future__ import division
-from numpy import linspace, logspace, log, exp, log10
-from gpfit.fit import fit
-
-m = 501
-
-u = logspace(0,log10(3),501)
-#u = u.reshape(u.size,1)
-w = (u**2 + 3)/(u+1)**2
-#w = w.reshape(w.size,1)
-
-x = log(u)
-y = log(w)
diff --git a/gpfit/examples/ex63.py b/gpfit/examples/ex63.py
deleted file mode 100644
index 006f2b3..0000000
--- a/gpfit/examples/ex63.py
+++ /dev/null
@@ -1,14 +0,0 @@
-import unittest
-from numpy import log, exp, log10, vstack
-from gpfit.fit import fit
-from numpy.random import random_sample
-
-Vdd = random_sample(1000,) + 1
-Vth = 0.2*random_sample(1000,) + 0.2
-P = Vdd**2 + 30*Vdd*exp(-(Vth-0.06*Vdd)/0.039)
-u = vstack((Vdd,Vth))
-x = log(u)
-y = log(P)
-K = 4
-
-cstrt, rmsErr = fit(x,y,K,"ISMA")
\ No newline at end of file
diff --git a/gpfit/fit.py b/gpfit/fit.py
index f67c4a6..a5b4cae 100644
--- a/gpfit/fit.py
+++ b/gpfit/fit.py
@@ -1,46 +1,27 @@
 "Implements the all-important 'fit' function."
 from numpy import ones, exp, sqrt, mean, square, hstack
-from .implicit_softmax_affine import implicit_softmax_affine
-from .softmax_affine import softmax_affine
-from .max_affine import max_affine
-from .levenberg_marquardt import levenberg_marquardt
-from .ba_init import ba_init
-from .print_fit import print_ISMA, print_SMA, print_MA
-from .fit_constraintset import FitCS
-
-ALPHA_INIT = 10
-RFUN = {"ISMA": implicit_softmax_affine,
-        "SMA": softmax_affine,
-        "MA": max_affine}
-
-# pylint: disable=invalid-name
-def get_params(ftype, K, xdata, ydata):
-    "Perform least-squares fitting optimization."
-    def rfun(params):
-        "A specific residual function."
-        [yhat, drdp] = RFUN[ftype](xdata, params)
-        r = yhat - ydata
-        return r, drdp
+from .classes import max_affine, softmax_affine, implicit_softmax_affine
+from .least_squares import levenberg_marquardt
+from .initialize import get_initial_parameters
+from .print_fit import print_isma, print_sma, print_ma
+from .constraint_set import FitConstraintSet
 
-    ba = ba_init(xdata, ydata.reshape(ydata.size, 1), K).flatten('F')
+CLASSES = {
+    "ISMA": implicit_softmax_affine,
+    "SMA": softmax_affine,
+    "MA": max_affine,
+}
 
-    if ftype == "ISMA":
-        params, _ = levenberg_marquardt(rfun, hstack((ba, ALPHA_INIT*ones(K))))
-    elif ftype == "SMA":
-        params, _ = levenberg_marquardt(rfun, hstack((ba, ALPHA_INIT)))
-    else:
-        params, _ = levenberg_marquardt(rfun, ba)
-
-    return params
 
 # pylint: disable=too-many-locals
 # pylint: disable=too-many-branches
 # pylint: disable=import-error
-def fit(xdata, ydata, K, ftype="ISMA"):
+def fit(xdata, ydata, K, ftype="ISMA", alpha0=10):
     """
     Fit a log-convex function to multivariate data, returning a FitConstraintSet
 
-    INPUTS
+    Arguments
+    ---------
         xdata:      Independent variable data
                         2D numpy array [nDim, nPoints]
                         [[<--------- x1 ------------->]
@@ -50,25 +31,27 @@ def fit(xdata, ydata, K, ftype="ISMA"):
                         1D numpy array [nPoints,]
                         [<---------- y ------------->]
 
+        K:          Number of terms
+
+        ftype:      Function class ["MA", "SMA", "ISMA"]
+
+        alpha0:     Initial guess for smoothing parameter alpha
+
+    Returns
+    -------
+        FitConstraintSet
+
         rms_error:  float
             Root mean square error between original (not log transformed)
             data and function fit.
     """
 
-    # Check data is in correct form
     if ydata.ndim > 1:
-        raise ValueError('Dependent data should be a 1D numpy array')
+        raise ValueError("Dependent data should be a 1D numpy array")
 
-    # Transform to column vector
     xdata = xdata.reshape(xdata.size, 1) if xdata.ndim == 1 else xdata.T
-
-    # Dimension of function (number of independent variables)
-    d = int(xdata.shape[1])
-
-    # begin fit data dict to generate fit constraint set
+    d = int(xdata.shape[1])  # Number of dimensions
     fitdata = {"ftype": ftype, "K": K, "d": d}
-
-    # find bounds of x
     if d == 1:
         fitdata["lb0"] = exp(min(xdata.reshape(xdata.size)))
         fitdata["ub0"] = exp(max(xdata.reshape(xdata.size)))
@@ -77,36 +60,49 @@ def fit(xdata, ydata, K, ftype="ISMA"):
             fitdata["lb%d" % i] = exp(min(xdata.T[i]))
             fitdata["ub%d" % i] = exp(max(xdata.T[i]))
 
-    params = get_params(ftype, K, xdata, ydata)
+    def residual(params):
+        "A specific residual function."
+        [yhat, drdp] = CLASSES[ftype](xdata, params)
+        r = yhat - ydata
+        return r, drdp
+
+    ba = get_initial_parameters(xdata, ydata.reshape(ydata.size, 1),
+                                K).flatten("F")
+    if ftype == "ISMA":
+        params, _ = levenberg_marquardt(residual, hstack((ba, alpha0*ones(K))))
+    elif ftype == "SMA":
+        params, _ = levenberg_marquardt(residual, hstack((ba, alpha0)))
+    else:
+        params, _ = levenberg_marquardt(residual, ba)
 
     # A: exponent parameters, B: coefficient parameters
-    A = params[[i for i in range(K*(d+1)) if i % (d + 1) != 0]]
-    B = params[[i for i in range(K*(d+1)) if i % (d + 1) == 0]]
+    A = params[[i for i in range(K*(d + 1)) if i % (d + 1) != 0]]
+    B = params[[i for i in range(K*(d + 1)) if i % (d + 1) == 0]]
 
     if ftype == "ISMA":
-        alpha = 1./params[list(range(-K, 0))]
+        alpha = 1/params[list(range(-K, 0))]
         for k in range(K):
             fitdata["c%d" % k] = exp(alpha[k]*B[k])
             fitdata["a%d" % k] = alpha[k]
             for i in range(d):
-                fitdata["e%d%d" % (k, i)] = alpha[k]*A[d*k+i]
-        print_ISMA(A, B, alpha, d, K)
+                fitdata["e%d%d" % (k, i)] = alpha[k]*A[d*k + i]
+        print_isma(A, B, alpha, d, K)
     elif ftype == "SMA":
-        alpha = 1./params[-1]
+        alpha = 1/params[-1]
         fitdata["a1"] = alpha
         for k in range(K):
             fitdata["c%d" % k] = exp(alpha*B[k])
             for i in range(d):
-                fitdata["e%d%d" % (k, i)] = alpha*A[d*k+i]
-        print_SMA(A, B, alpha, d, K)
+                fitdata["e%d%d" % (k, i)] = alpha*A[d*k + i]
+        print_sma(A, B, alpha, d, K)
     elif ftype == "MA":
         alpha = 1
         fitdata["a1"] = 1
         for k in range(K):
             fitdata["c%d" % k] = exp(B[k])
             for i in range(d):
-                fitdata["e%d%d" % (k, i)] = A[d*k+i]
-        print_MA(A, B, d, K)
+                fitdata["e%d%d" % (k, i)] = A[d*k + i]
+        print_ma(A, B, d, K)
 
     if min(exp(B*alpha)) < 1e-100:
         raise ValueError("Fitted constraint contains too small a value...")
@@ -114,26 +110,14 @@ def fit(xdata, ydata, K, ftype="ISMA"):
         raise ValueError("Fitted constraint contains too large a value...")
 
     def evaluate(xdata):
-        """
-        Evaluate the y of this fit over a range of xdata.
-
-        INPUTS
-            xdata:      Independent variable data
-                            2D numpy array [nDim, nPoints]
-                            [[<--------- x1 ------------->]
-                             [<--------- x2 ------------->]]
-
-        OUTPUT
-            ydata:      Dependent variable data in 1D numpy array
-        """
+        """Evaluate the fit over a range of xdata."""
         xdata = xdata.reshape(xdata.size, 1) if xdata.ndim == 1 else xdata.T
-        return RFUN[ftype](xdata, params)[0]
+        return CLASSES[ftype](xdata, params)[0]
 
-    # cstrt.evaluate = evaluate
-    fitdata["rms_err"] = sqrt(mean(square(evaluate(xdata.T)-ydata)))
-    fitdata["max_err"] = sqrt(max(square(evaluate(xdata.T)-ydata)))
+    fitdata["rms_err"] = sqrt(mean(square(evaluate(xdata.T) - ydata)))
+    fitdata["max_err"] = sqrt(max(square(evaluate(xdata.T) - ydata)))
 
-    cs = FitCS(fitdata)
+    cs = FitConstraintSet(fitdata)
     cs.evaluate = evaluate
 
     return cs, cs.rms_err
diff --git a/gpfit/fit_constraintset.py b/gpfit/fit_constraintset.py
deleted file mode 100644
index 17835cf..0000000
--- a/gpfit/fit_constraintset.py
+++ /dev/null
@@ -1,190 +0,0 @@
-" fit constraint set "
-from __future__ import print_function
-from __future__ import division
-from gpkit import ConstraintSet
-from gpkit import Variable, NomialArray, NamedVariables, VectorVariable
-from numpy import amax, array, hstack, where
-
-# pylint: disable=too-many-instance-attributes, too-many-locals,
-# pylint: disable=too-many-branches, no-member, import-error
-# pylint: disable=too-many-arguments
-
-
-class FitCS(ConstraintSet):
-    """Constraint set for fitted functions
-
-    Arguments
-    ---------
-    fitdata : dict
-        dictionary of fit parameters
-    ivar : gpkit Variable, Monomial, or NomialArray
-        independent variable
-    dvars : list of gpkit Variables, Monomials, or NomialArrays
-        dependent variables
-    err_margin : string, either "Max" or "RMS"
-        flag to set margin factor using RMS or max error
-
-    """
-    def __init__(self, fitdata, ivar=None, dvars=None, name="",
-                 err_margin=None):
-
-        self.fitdata = fitdata
-
-        if ivar is None:
-            with NamedVariables("fit"):
-                dvars = VectorVariable(fitdata["d"], "u")
-                ivar = Variable("w")
-
-        self.dvars = dvars
-        self.ivar = ivar
-        self.rms_err = fitdata["rms_err"]
-        self.max_err = fitdata["max_err"]
-
-        monos = [fitdata["c%d" % k]*NomialArray(array(dvars).T**array(
-            [fitdata["e%d%d" % (k, i)] for i in
-             range(fitdata["d"])])).prod(NomialArray(dvars).ndim - 1)
-                 for k in range(fitdata["K"])]
-
-        if err_margin == "Max":
-            self.mfac = Variable("m_{fac-" + name + "-fit}",
-                                 1 + self.max_err, "-", "max error of fit")
-        elif err_margin == "RMS":
-            self.mfac = Variable("m_{fac-" + name + "-fit}",
-                                 1 + self.rms_err, "-", "RMS error of fit")
-        elif err_margin is None:
-            self.mfac = Variable("m_{fac-" + name + "-fit}", 1.0, "-",
-                                 "fit factor")
-        else:
-            raise ValueError("Invalid name for err_margin: valid inputs Max, "
-                             "RMS")
-
-        if fitdata["ftype"] == "ISMA":
-            # constraint of the form 1 >= c1*u1^exp1*u2^exp2*w^(-alpha) + ....
-            alpha = array([fitdata["a%d" % k] for k in
-                           range(fitdata["K"])])
-            lhs, rhs = 1, NomialArray(monos/(ivar/self.mfac)**alpha).sum(0)
-        elif fitdata["ftype"] == "SMA":
-            # constraint of the form w^alpha >= c1*u1^exp1 + c2*u2^exp2 +....
-            alpha = fitdata["a1"]
-            lhs, rhs = (ivar/self.mfac)**alpha, NomialArray(monos).sum(0)
-        elif fitdata["ftype"] == "MA":
-            # constraint of the form w >= c1*u1^exp1, w >= c2*u2^exp2, ....
-            lhs, rhs = (ivar/self.mfac), NomialArray(monos).T
-
-        if fitdata["K"] == 1:
-            # when possible, return an equality constraint
-            if hasattr(rhs, "shape"):
-                if rhs.ndim > 1:
-                    self.constraint = [(lh == rh) for lh, rh in zip(lhs, rhs)]
-                else:
-                    self.constraint = (lhs == rhs)
-            else:
-                self.constraint = (lhs == rhs)
-        else:
-            if hasattr(rhs, "shape"):
-                if rhs.ndim > 1:
-                    self.constraint = [(lh >= rh) for lh, rh in zip(lhs, rhs)]
-                else:
-                    self.constraint = (lhs >= rhs)
-            else:
-                self.constraint = (lhs >= rhs)
-
-        self.bounds = {}
-        for i, dvar in enumerate(self.dvars):
-            self.bounds[dvar] = [fitdata["lb%d" % i],
-                                 fitdata["ub%d" % i]]
-
-        ConstraintSet.__init__(self, [self.constraint])
-
-    def get_fitdata(self):
-        " return fit data "
-        return self.fitdata
-
-    def get_dataframe(self):
-        " return a pandas DataFrame of fit parameters "
-        import pandas as pd
-        df = pd.DataFrame(list(self.fitdata.values())).transpose()
-        df.columns = list(self.fitdata.keys())
-        return df
-
-    def process_result(self, result):
-        """
-        make sure fit result is within bounds of fitted data
-        """
-        super(FitCS, self).process_result(result)
-
-        if self.mfac not in result["sensitivities"]["constants"]:
-            return
-        if amax([abs(result["sensitivities"]["constants"][self.mfac])]) < 1e-5:
-            return
-
-        for dvar in self.dvars:
-            if isinstance(dvar, NomialArray):
-                num = [result(x) for x in dvar]
-            else:
-                num = result(dvar)
-            direct = None
-            if any(x < self.bounds[dvar][0] for x in hstack([num])):
-                direct, state = "lower", "below"
-                bnd = self.bounds[dvar][0]
-            if any(x > self.bounds[dvar][1] for x in hstack([num])):
-                direct, state = "upper", "above"
-                bnd = self.bounds[dvar][1]
-
-            if direct:
-                msg = ("Variable %.100s could cause inaccurate result"
-                       " because it is %s" % (dvar, state)
-                       + " %s bound. Solution is %.4f but"
-                       " bound is %.4f" %
-                       (direct, amax([num]), bnd))
-                print("Warning: " + msg)
-
-
-class XfoilFit(FitCS):
-    """Special FitCS that can post-solve compare result to XFOIL
-
-    Arguments (in addition to the arguments to FitCS)
-    ---------
-    airfoil:                            airfoil of fitted data
-                                            str (e.g. "xxx.dat", "naca xxxx")
-
-    """
-    def __init__(self, fitdata, ivar=None, dvars=None, name="",
-                 err_margin=None, airfoil=False):
-
-        super(XfoilFit, self).__init__(fitdata, ivar=ivar, dvars=dvars,
-                                       name=name, err_margin=err_margin)
-
-        self.airfoil = airfoil
-
-    def process_result(self, result):
-        """
-        if data comes from Xfoil and airfoil is provided check against xfoil
-        """
-        super(XfoilFit, self).process_result(result)
-
-        if self.mfac not in result["sensitivities"]["constants"]:
-            return
-        if amax([abs(result["sensitivities"]["constants"][self.mfac])]) < 1e-5:
-            return
-        if not self.airfoil:
-            return
-
-        from .xfoilWrapper import xfoil_comparison
-        cl, re = 0.0, 0.0
-        for dvar in self.dvars:
-            if "Re" in str(dvar):
-                re = result(dvar)
-            if "C_L" in str(dvar):
-                cl = result(dvar)
-        cd = result(self.ivar)
-        if not hasattr(cl, "__len__") and hasattr(re, "__len__"):
-            cl = [cl]*len(re)
-        err, cdx = xfoil_comparison(self.airfoil, cl, re, cd)
-        ind = where(err > 0.05)[0]
-        for i in ind:
-            msg = ("Drag error for %s is %.2f. Re=%.1f; CL=%.4f;"
-                   " Xfoil cd=%.6f, GP sol cd=%.6f" %
-                   (", ".join(self.ivar.descr["models"]), err[i], re[i],
-                    cl[i], cd[i], cdx[i]))
-            print("Warning: %s" % msg)
diff --git a/gpfit/implicit_softmax_affine.py b/gpfit/implicit_softmax_affine.py
deleted file mode 100644
index 57671b6..0000000
--- a/gpfit/implicit_softmax_affine.py
+++ /dev/null
@@ -1,48 +0,0 @@
-"Implements ISMA residual function"
-from __future__ import division
-from numpy import ones, nan, inf, hstack, dot, tile
-from .lse_implicit import lse_implicit
-
-
-# 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.
-
-    INPUTS:
-            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]
-
-    OUTPUTS:
-            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 inf*ones((npt, 1)), nan
-    ba = ba.reshape(dimx+1, K, order='F')  # reshape ba to matrix
-
-    X = hstack((ones((npt, 1)), x))  # augment data with column of ones
-    z = dot(X, ba)  # compute affine functions
-
-    y, dydz, dydalpha = lse_implicit(z, alpha)
-
-    nrow, ncol = dydz.shape
-    repmat = tile(dydz, (dimx+1, 1)).reshape(nrow, ncol*(dimx+1), order='F')
-    dydba = repmat * tile(X, (1, K))
-    dydp = hstack((dydba, dydalpha))
-
-    return y, dydp
diff --git a/gpfit/ba_init.py b/gpfit/initialize.py
similarity index 70%
rename from gpfit/ba_init.py
rename to gpfit/initialize.py
index 0493368..82d02ef 100644
--- a/gpfit/ba_init.py
+++ b/gpfit/initialize.py
@@ -1,48 +1,43 @@
-"Implements ba_init"
-from __future__ import print_function
+"Implements get_initial_parameters"
 from numpy import ones, hstack, zeros, tile, argmin
 from numpy.linalg import lstsq, matrix_rank
 from numpy.random import permutation as randperm
 
 
 # pylint: disable=too-many-locals
-def ba_init(x, y, K):
+def get_initial_parameters(x, y, K):
     """
     Initializes max-affine fit to data (y, x)
     ensures that initialization has at least K+1 points per partition (i.e.
     per affine function)
 
-    INPUTS:
+    Arguments:
+    ----------
         x:      Independent variable data
                     2D column vector [nPoints x nDims]
 
         y:      Dependent variable data
                     2D column vector [nPoints x 1]
 
-    OUTPUTS:
+    Returns:
+    --------
         ba:     Initial b and a parameters
                     2D array [(dimx+1) x k]
 
     """
-    defaults = {}
-    defaults['bverbose'] = False
-    options = defaults
 
     npt, dimx = x.shape
+    if K*(dimx + 1) > npt:
+        raise ValueError("Not enough data points")
 
     X = hstack((ones((npt, 1)), x))
-    b = zeros((dimx+1, K))
-
-    if K*(dimx+1) > npt:
-        raise ValueError('Not enough data points')
-
-    # Choose K unique indices
-    randinds = randperm(npt)[0:K]
+    b = zeros((dimx + 1, K))
+    randinds = randperm(npt)[0:K]  # Choose K unique indices
 
     # partition based on distances
     sqdists = zeros((npt, K))
     for k in range(K):
-        sqdists[:, k] = ((x - tile(x[randinds[k], :], (npt, 1))) ** 2).sum(1)
+        sqdists[:, k] = ((x - tile(x[randinds[k], :], (npt, 1)))**2).sum(1)
 
     # index to closest k for each data pt
     mindistind = argmin(sqdists, axis=1)
@@ -50,7 +45,8 @@ def ba_init(x, y, K):
     # loop through each partition, making local fits
     # note we expand partitions that result in singular least squares problems
     # why this way? some points will be shared by multiple partitions, but
-    # resulting max-affine fit will tend to be good. (as opposed to solving least-norm version)
+    # resulting max-affine fit will tend to be good. (as opposed to solving
+    # least-norm version)
     for k in range(K):
         inds = mindistind == k
 
@@ -60,25 +56,20 @@ def ba_init(x, y, K):
             sortdistind = sqdists[:, k].argsort()
 
             i = sum(inds)  # number of points in partition
-            iinit = i
 
-            if i < dimx+1:
+            if i < dimx + 1:
                 # obviously, at least need dimx+1 points. fill these in before
                 # checking any ranks
-                inds[sortdistind[i+1:dimx+1]] = 1  # TODO: check index
-                i = dimx+1  # TODO: check index
+                inds[sortdistind[i + 1 : dimx + 1]] = 1  # TODO: check index
+                i = dimx + 1  # TODO: check index
 
             # now add points until rank condition satisfied
-            while matrix_rank(X[inds, :]) < dimx+1:
-                i = i+1
+            while matrix_rank(X[inds, :]) < dimx + 1:
+                i = i + 1
                 inds[sortdistind[i]] = 1
 
-            if options['bverbose']:
-                print("ba_init: Added %s points to partition %s to maintain"
-                      "full rank for local fitting." % (i-iinit, k))
         # now create the local fit
         b[:, k] = lstsq(X[inds.nonzero()], y[inds.nonzero()], rcond=-1)[0][:, 0]
         # Rank condition specified to default for python upgrades
 
-
     return b
diff --git a/gpfit/levenberg_marquardt.py b/gpfit/least_squares.py
similarity index 73%
rename from gpfit/levenberg_marquardt.py
rename to gpfit/least_squares.py
index b750cb2..3cc4a18 100644
--- a/gpfit/levenberg_marquardt.py
+++ b/gpfit/least_squares.py
@@ -1,6 +1,4 @@
-"Implements LM"
-from __future__ import print_function
-from __future__ import division
+"Implements the Levenberg-Marquardt algorithm"
 from time import time
 from sys import float_info
 import numpy as np
@@ -8,20 +6,23 @@
 from scipy.sparse import spdiags, issparse
 
 
-# pylint: disable=too-many-locals,too-many-arguments,too-many-branches,too-many-statements
-def levenberg_marquardt(residfun, initparams,
-                        verbose=False,
-                        lambdainit=0.02,
-                        maxiter=5000,
-                        maxtime=5.,
-                        tolgrad=np.sqrt(float_info.epsilon),
-                        tolrms=1e-7):
+# pylint: disable=too-many-locals,too-many-arguments,too-many-branches,too-many-statements,no-else-break
+def levenberg_marquardt(
+    residfun,
+    initparams,
+    verbose=False,
+    lambdainit=0.02,
+    maxiter=5000,
+    maxtime=5.0,
+    tolgrad=np.sqrt(float_info.epsilon),
+    tolrms=1e-7,
+):
     """
     Levenberg-Marquardt alogrithm
     Minimizes sum of squared error of residual function
 
-    INPUTS
-    ------
+    Arguments
+    ---------
     residfun: function
         Mapping from parameters to residuals, of the form
         (r, drdp) = residfun(params)
@@ -43,11 +44,11 @@ def levenberg_marquardt(residfun, initparams,
     tolrms: float
         Tolerance on change in rms error per iteration
 
-    OUTPUTS
+    Returns
     -------
     params: np.array (1D)
         Parameter vector that locally minimizes norm(residfun, 2)
-    RMStraj: np.array
+    rmstraj: np.array
         History of RMS errors after each step (first point is initialization)
     """
 
@@ -57,12 +58,13 @@ def levenberg_marquardt(residfun, initparams,
     nparam = initparams.size
 
     if initparams.ndim > 1:
-        raise ValueError('params should be a vector')
+        raise ValueError("params should be a vector")
 
     # Define display formatting if required
-    formatstr1 = '%5.0f        %9.6g        %9.3g\n'
-    formatstr = ('%5.0f        %9.6g        %9.3g        '
-                 '%12.4g        %12.4g        %8.4g\n')
+    formatstr1 = "%5.0f        %9.6g        %9.3g\n"
+    formatstr = (
+        "%5.0f        %9.6g        %9.3g        %12.4g        %12.4g        %8.4g\n"
+    )
 
     # Get residual values and jacobian at initial point; extract size info
     params = initparams
@@ -72,8 +74,9 @@ def levenberg_marquardt(residfun, initparams,
     r.shape = (npt, 1)  # Make r into column vector
 
     if J.shape != (npt, nparam):
-        raise ValueError('Jacobian size %s inconsistent with (%s, %s)' %
-                         (J.shape, npt, nparam))
+        raise ValueError(
+            "Jacobian size %s inconsistent with (%s, %s)" % (J.shape, npt, nparam)
+        )
 
     # "Accept" initial point
     rms = norm(r)/np.sqrt(npt)  # 2-norm
@@ -86,35 +89,32 @@ def levenberg_marquardt(residfun, initparams,
     diagJJ = sum(J*J, 0).T
     zeropad = np.zeros((nparam, 1))
     lamb = lambdainit
-    RMStraj = [rms]
+    rmstraj = [rms]
 
     # Display info for 1st iter
     if verbose:
-        print('\n                    First-Order                        '
-              'Norm of \n')
-        print('Iter        Residual        optimality            Lambda'
-              '            step        Jwarp \n')
+        print("\n                    First-Order                        " "Norm of \n")
+        print(
+            "Iter        Residual        optimality            Lambda"
+            "            step        Jwarp \n"
+        )
         print(formatstr1 % (itr, rms, maxgrad))
 
     # Main Loop
     while True:
-
         if itr == maxiter:
             if verbose:
-                print('Reached maximum number of iterations')
+                print("Reached maximum number of iterations")
             break
-
         elif time() - t > maxtime:
             if verbose:
-                print('Reached maxtime (%s seconds)' % maxtime)
+                print("Reached maxtime (%s seconds)" % maxtime)
             break
-        elif (itr >= 2 and
-              abs(RMStraj[itr] - RMStraj[itr-2]) <
-              RMStraj[itr]*tolrms):
+        elif itr >= 2 and abs(rmstraj[itr] - rmstraj[itr - 2]) < rmstraj[itr]*tolrms:
             # Should really only allow this exit case
             # if trust region constraint is slack
             if verbose:
-                print('RMS changed less than tolrms')
+                print("RMS changed less than tolrms")
             break
 
         itr += 1
@@ -145,7 +145,7 @@ def levenberg_marquardt(residfun, initparams,
         # Check function value at trialp
         trialr, trialJ = residfun(trialp)
         trialrms = norm(trialr)/np.sqrt(npt)
-        RMStraj.append(trialrms)
+        rmstraj.append(trialrms)
 
         # Accept or reject trial params
         if trialrms < rms:
@@ -162,11 +162,11 @@ def levenberg_marquardt(residfun, initparams,
 
             if maxgrad < tolgrad:
                 if verbose:
-                    print('1st order optimality attained')
+                    print("1st order optimality attained")
                 break
 
             if prev_trial_accepted and itr > 1:
-                lamb = lamb/10.
+                lamb = lamb/10
 
             prev_trial_accepted = True
             params_updated = True
@@ -178,9 +178,9 @@ def levenberg_marquardt(residfun, initparams,
             prev_trial_accepted = False
             params_updated = False
 
-    assert len(RMStraj) == itr + 1
-    RMStraj = np.array(RMStraj)
+    assert len(rmstraj) == itr + 1
+    rmstraj = np.array(rmstraj)
     if verbose:
-        print('Final RMS: ' + repr(rms))
+        print("Final RMS: " + repr(rms))
 
-    return params, RMStraj
+    return params, rmstraj
diff --git a/gpfit/lse_implicit.py b/gpfit/logsumexp.py
similarity index 54%
rename from gpfit/lse_implicit.py
rename to gpfit/logsumexp.py
index 4907ea6..afd3f92 100644
--- a/gpfit/lse_implicit.py
+++ b/gpfit/logsumexp.py
@@ -1,6 +1,4 @@
-"Implements lse_implicit"
-from __future__ import print_function
-from __future__ import division
+"Module for log-sum-exponential functions"
 from numpy import zeros, spacing, exp, log, tile
 
 
@@ -13,14 +11,16 @@ def lse_implicit(x, alpha):
     - lse_implicit is a mapping R^n --> R, where n number of dimensions
     - implementation: newton raphson steps to find f(x,y) = 0
 
-    INPUTS:
+    Arguments:
+    ----------
             x:      independent variable data
                         2D numpy array [nPoints x nDim]
 
             alpha:  local softness parameter
                         1D array [K] (K=number of terms)
 
-    OUTPUTS:
+    Returns:
+    --------
             y:      ISMA approximation to log transformed data
                         1D numpy array [nPoints]
 
@@ -31,29 +31,22 @@ def lse_implicit(x, alpha):
                         2D array [nPoints x nDim]
     """
 
-    bverbose = False
-
     tol = 10*spacing(1)
     npt, nx = x.shape
 
     if nx != alpha.size:
-        raise ValueError('alpha size mismatch')
+        raise ValueError("alpha size mismatch")
 
     alphamat = tile(alpha, (npt, 1))
 
     m = x.max(1)  # maximal x values
     # distance from m; note h <= 0 for all entries
     h = x - (tile(m, (nx, 1))).T
-    # (unless x has infs; in this case h has nans;
-    # can prob deal w/ this gracefully)
-    # should also deal with alpha==inf case gracefully
-
     L = zeros((npt,))  # initial guess. note y = m + L
-
     Lmat = (tile(L, (nx, 1))).T
 
     # initial eval
-    expo = exp(alphamat*(h-Lmat))
+    expo = exp(alphamat*(h - Lmat))
     alphaexpo = alphamat*expo
     sumexpo = expo.sum(axis=1)
     sumalphaexpo = alphaexpo.sum(axis=1)
@@ -63,27 +56,60 @@ def lse_implicit(x, alpha):
     i = abs(f) > tol  # inds to update
 
     while any(i):
-        L[i] = L[i] - f[i]/dfdL[i]    # newton step
+        L[i] = L[i] - f[i]/dfdL[i]  # newton step
         # re-evaluate
         Lmat[i, :] = (tile(L[i], (nx, 1))).T
-        expo[i, :] = exp(alphamat[i, :] * (h[i, :] - Lmat[i, :]))
-        alphaexpo[i, :] = alphamat[i, :] * expo[i, :]
+        expo[i, :] = exp(alphamat[i, :]*(h[i, :] - Lmat[i, :]))
+        alphaexpo[i, :] = alphamat[i, :]*expo[i, :]
         sumexpo[i] = expo[i, :].sum(axis=1)
         sumalphaexpo[i, ] = alphaexpo[i, :].sum(axis=1)
         f[i] = log(sumexpo[i])
         dfdL[i] = -sumalphaexpo[i, ]/sumexpo[i]
         neval = neval + 1
-
-        # update inds that need to be evaluated
-        i[i] = abs(f[i]) > tol
-
-    if bverbose:
-        print('lse_implicit converged in ' +
-              repr(neval) + ' newton-raphson steps')
+        i[i] = abs(f[i]) > tol  # update inds that need to be evaluated
 
     y = m + L
-
     dydx = alphaexpo/(tile(sumalphaexpo, (nx, 1))).T
     dydalpha = (h - Lmat)*expo/(tile(sumalphaexpo, (nx, 1))).T
 
     return y, dydx, dydalpha
+
+
+def lse_scaled(x, alpha):
+    """
+    Log-sum-exponential function with derivatives
+    - sums across the second dimension of x
+    - note that lse_scaled is a mapping R^n --> R
+
+    Arguments:
+    ----------
+            x:      independent variable data
+                        2D numpy array [nPoints x nDim]
+
+            alpha:  local softness parameter
+                        1D array [K] (K=number of terms)
+
+    Returns:
+    --------
+            y:      ISMA approximation to log transformed data
+                        1D numpy array [nPoints]
+
+            dydx:   Deriv of y wrt each x
+                        2D numpy array[nPoints x nDim]
+
+            dydalpha:
+                    [n-element 1D array], n is number of data points
+    """
+
+    _, n = x.shape
+    m = x.max(axis=1)  # maximal x values
+    h = x - (tile(m, (n, 1))).T  # distance from m; note h <= 0 for all entries
+    expo = exp(alpha*h)
+    sumexpo = expo.sum(axis=1)
+    L = log(sumexpo)/alpha
+    y = L + m
+    dydx = expo/(tile(sumexpo, (n, 1))).T
+    # note that sum(dydx,2)==1, i.e. dydx is a probability distribution
+    dydalpha = ((h*expo).sum(axis=1)/sumexpo - L)/alpha
+
+    return y, dydx, dydalpha
diff --git a/gpfit/lse_scaled.py b/gpfit/lse_scaled.py
deleted file mode 100644
index 3220e5f..0000000
--- a/gpfit/lse_scaled.py
+++ /dev/null
@@ -1,50 +0,0 @@
-"Implements lse_scaled"
-from __future__ import division
-from numpy import tile, exp, log
-
-
-def lse_scaled(x, alpha):
-    """
-    Log-sum-exponential function with derivatives
-    - sums across the second dimension of x
-    - note that lse_scaled is a mapping R^n --> R
-
-    INPUTS:
-            x:      independent variable data
-                        2D numpy array [nPoints x nDim]
-
-            alpha:  local softness parameter
-                        1D array [K] (K=number of terms)
-
-    OUTPUTS:
-            y:      ISMA approximation to log transformed data
-                        1D numpy array [nPoints]
-
-            dydx:   Deriv of y wrt each x
-                        2D numpy array[nPoints x nDim]
-
-            dydalpha:
-                    [n-element 1D array], n is number of data points
-    """
-
-    _, n = x.shape
-
-    m = x.max(axis=1)  # maximal x values
-
-    h = x - (tile(m, (n, 1))).T  # distance from m; note h <= 0 for all entries
-    # (unless x has infs; in this case h has nans;
-    # can prob deal w/ this gracefully)
-    # should also deal with alpha==inf case gracefully
-
-    expo = exp(alpha*h)
-    sumexpo = expo.sum(axis=1)
-
-    L = log(sumexpo)/alpha
-
-    y = L + m
-
-    dydx = expo/(tile(sumexpo, (n, 1))).T
-    # note that sum(dydx,2)==1, i.e. dydx is a probability distribution
-    dydalpha = ((h*expo).sum(axis=1)/sumexpo - L)/alpha
-
-    return y, dydx, dydalpha
diff --git a/gpfit/max_affine.py b/gpfit/max_affine.py
deleted file mode 100644
index fef9433..0000000
--- a/gpfit/max_affine.py
+++ /dev/null
@@ -1,46 +0,0 @@
-"Implements MA residual function"
-from __future__ import division
-import numpy as np
-
-
-def max_affine(x, ba):
-    """
-    Evaluates max affine function at values of x, given a set of
-    max affine fit parameters.
-
-    INPUTS
-    ------
-        x: 2D array [nPoints x nDim]
-            Independent variable data
-
-        ba: 2D array
-            max affine fit parameters
-            [[b1, a11, ... a1k]
-             [ ....,          ]
-             [bk, ak1, ... akk]]
-
-    OUTPUTS
-    -------
-        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
-
-    # augment data with column of ones
-    X = np.hstack((np.ones((npt, 1)), x))
-
-    y, partition = np.dot(X, ba).max(1), np.dot(X, ba).argmax(1)
-
-    # The not-sparse sparse version
-    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
diff --git a/gpfit/plot_fit.py b/gpfit/plot_fit.py
index fe28bc2..9c050c7 100644
--- a/gpfit/plot_fit.py
+++ b/gpfit/plot_fit.py
@@ -1,59 +1,59 @@
 "Fit plotting"
-from __future__ import division
 import numpy as np
 import matplotlib.pyplot as plt
 from gpfit.fit import fit
-from gpfit.print_fit import print_MA, print_SMA
+from gpfit.print_fit import print_ma, print_sma
+
 
 # pylint: disable=invalid-name
 # pylint: disable=too-many-locals
-def plot_fit_1d(udata, wdata, K=1, fitclass='SMA', plotspace='log'):
+def plot_fit_1d(udata, wdata, K=1, fitclass="SMA", plotspace="log"):
     "Finds and plots a fit (MA or SMA) for 1D data"
 
     cstrt, _ = fit(np.log(udata), np.log(wdata), K, fitclass)
     uu = np.linspace(min(udata), max(udata), 1000)
 
-    if fitclass == 'MA':
-        uvarkey, = cstrt[0].left.varkeys
+    if fitclass == "MA":
+        (uvarkey,) = cstrt[0].left.varkeys
         A = [c.left.exps[0][uvarkey] for c in cstrt]
         B = np.log([c.left.cs[0] for c in cstrt])
         WW = []
         for k in range(K):
             WW += [np.exp(B[k])*uu**A[k]]
-        stringlist = print_MA(A, B, 1, K)
-
+        stringlist = print_ma(A, B, 1, K)
 
-    if fitclass == 'SMA':
-        wexps, = cstrt[0].left.exps
+    if fitclass == "SMA":
+        (wexps,) = cstrt[0].left.exps
         alpha = list(wexps.values())[0]
-        uvarkey, = cstrt[0].right.varkeys
+        (uvarkey,) = cstrt[0].right.varkeys
         A = [d[uvarkey]/alpha for d in cstrt[0].right.exps]
-        B = np.log(cstrt[0].right.cs) / alpha
+        B = np.log(cstrt[0].right.cs)/alpha
 
         ww = 0
         for k in range(K):
-            ww += (np.exp(alpha*B[k])*uu**(alpha*A[k]))
-        WW = [ww**(1./alpha)]
+            ww += np.exp(alpha*B[k])*uu**(alpha*A[k])
+        WW = [ww**(1/alpha)]
 
-        print_str = print_SMA(A, B, alpha, 1, K)
-        stringlist = [''.join(print_str)]
+        print_str = print_sma(A, B, alpha, 1, K)
+        stringlist = ["".join(print_str)]
 
     f, ax = plt.subplots()
-    if plotspace == 'log':
-        ax.loglog(udata, wdata, '+r')
+    if plotspace == "log":
+        ax.loglog(udata, wdata, "+r")
         for ww in WW:
             ax.loglog(uu, ww)
-    elif plotspace == 'linear':
-        ax.plot(udata, wdata, '+r')
+    elif plotspace == "linear":
+        ax.plot(udata, wdata, "+r")
         for ww in WW:
             ax.plot(uu, ww)
-    ax.set_xlabel('u')
-    ax.legend(['Data'] + stringlist, loc='best')
+    ax.set_xlabel("u")
+    ax.legend(["Data"] + stringlist, loc="best")
 
     return f, ax
 
+
 if __name__ == "__main__":
     N = 51
     U = np.logspace(0, np.log10(3), N)
-    W = (U**2+3) / (U+1)**2
-    plot_fit_1d(U, W, K=2, fitclass='SMA', plotspace="linear")
+    W = (U**2 + 3)/(U + 1)**2
+    plot_fit_1d(U, W, K=2, fitclass="SMA", plotspace="linear")
diff --git a/gpfit/print_fit.py b/gpfit/print_fit.py
index 1dcf696..a55d656 100644
--- a/gpfit/print_fit.py
+++ b/gpfit/print_fit.py
@@ -1,67 +1,60 @@
 "Implements functions for raw fit printing from params"
-from __future__ import print_function
 from numpy import exp
 
 
 # pylint: disable=invalid-name
-def print_ISMA(A, B, alpha, d, K):
-    "prints ISMA fit from params"
-    print("ISMA fit from params")
-    stringList = [None]*K
+def print_isma(A, B, alpha, d, K):
+    """prints ISMA fit from params"""
 
-    printString = '1 = '
+    print("ISMA fit from params")
+    string_list = [None]*K
+    print_string = "1 = "
     for k in range(K):
         if k > 0:
-            print(printString)
-            printString = '    + '
-
-        printString += '({0:.6g}/w**{1:.6g})'.format(exp(alpha[k] * B[k]), alpha[k])
-
+            print(print_string)
+            print_string = "    + "
+        print_string += "({0:.6g}/w**{1:.6g})".format(exp(alpha[k]*B[k]), alpha[k])
         for i in range(d):
-            printString += ' * (u_{0:d})**{1:.6g}'.format(i+1, alpha[k] * A[d*k + i])
+            print_string += " * (u_{0:d})**{1:.6g}".format(
+                i + 1, alpha[k]*A[d*k + i]
+            )
+        string_list[k] = print_string
+    print(print_string)
 
-        stringList[k] = printString
-
-    print(printString)
-    return stringList
+    return string_list
 
 
 # pylint: disable=invalid-name
-def print_SMA(A, B, alpha, d, K):
-    "prints SMA fit from params"
-    print("SMA fit from params")
-    stringList = [None]*K
+def print_sma(A, B, alpha, d, K):
+    """prints SMA fit from params"""
 
-    printString = 'w**{0:.6g} = '.format(alpha)
+    print("SMA fit from params")
+    string_list = [None]*K
+    print_string = "w**{0:.6g} = ".format(alpha)
     for k in range(K):
         if k > 0:
-            print(printString)
-            printString = '    + '
-
-        printString += '{0:.6g}'.format(exp(alpha * B[k]))
-
+            print(print_string)
+            print_string = "    + "
+        print_string += "{0:.6g}".format(exp(alpha*B[k]))
         for i in range(d):
-            printString += ' * (u_{0:d})**{1:.6g}'.format(i+1, alpha * A[d*k + i])
+            print_string += " * (u_{0:d})**{1:.6g}".format(i + 1, alpha*A[d*k + i])
+        string_list[k] = print_string
+    print(print_string)
 
-        stringList[k] = printString
-
-    print(printString)
-    return stringList
+    return string_list
 
 
 # pylint: disable=invalid-name
-def print_MA(A, B, d, K):
-    "prints MA fit from params"
-    print("MA fit from params")
-    stringList = [None]*K
+def print_ma(A, B, d, K):
+    """prints MA fit from params"""
 
+    print("MA fit from params")
+    string_list = [None]*K
     for k in range(K):
-        printString = 'w = {0:.6g}'.format(exp(B[k]))
-
+        print_string = "w = {0:.6g}".format(exp(B[k]))
         for i in range(d):
-            printString += ' * (u_{0:d})**{1:.6g}'.format(i+1, A[d*k + i])
-
-        stringList[k] = printString
-        print(printString)
+            print_string += " * (u_{0:d})**{1:.6g}".format(i + 1, A[d*k + i])
+        string_list[k] = print_string
+        print(print_string)
 
-    return stringList
+    return string_list
diff --git a/gpfit/softmax_affine.py b/gpfit/softmax_affine.py
deleted file mode 100644
index 0d5e8f5..0000000
--- a/gpfit/softmax_affine.py
+++ /dev/null
@@ -1,51 +0,0 @@
-"Implements SMA residual function"
-from __future__ import division
-from numpy import size, inf, nan, ones, hstack, dot, tile
-from .lse_scaled import lse_scaled
-
-
-# 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.
-
-    INPUTS:
-            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]
-
-    OUTPUTS:
-            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 inf*ones((npt, 1)), nan
-    K = size(ba)//(dimx+1)
-    ba = ba.reshape(dimx+1, K, order='F')
-
-    X = hstack((ones((npt, 1)), x))  # augment data with column of ones
-    z = dot(X, ba)  # compute affine functions
-
-    y, dydz, dydsoftness = lse_scaled(z, alpha)
-
-    dydsoftness = -dydsoftness*(alpha**2)
-
-    nrow, ncol = dydz.shape
-    repmat = tile(dydz, (dimx+1, 1)).reshape(nrow, ncol*(dimx+1), order='F')
-    dydba = repmat * tile(X, (1, K))
-    dydsoftness.shape = (dydsoftness.size, 1)
-    dydp = hstack((dydba, dydsoftness))
-
-    return y, dydp
diff --git a/gpfit/xfoil/constraint_set.py b/gpfit/xfoil/constraint_set.py
new file mode 100644
index 0000000..51aa706
--- /dev/null
+++ b/gpfit/xfoil/constraint_set.py
@@ -0,0 +1,63 @@
+from gpfit.constraint_set import FitConstraintSet
+
+
+class XfoilFit(FitConstraintSet):
+    """Special FitConstraintSet that can post-solve compare result to XFOIL
+
+    Arguments (in addition to the arguments to FitConstraintSet)
+    ---------
+    airfoil:                            airfoil of fitted data
+                                            str (e.g. "xxx.dat", "naca xxxx")
+
+    """
+
+    def __init__(
+        self, fitdata, ivar=None, dvars=None, name="", err_margin=None, airfoil=False
+    ):
+
+        super(XfoilFit, self).__init__(
+            fitdata, ivar=ivar, dvars=dvars, name=name, err_margin=err_margin
+        )
+
+        self.airfoil = airfoil
+
+    def process_result(self, result):
+        """
+        if data comes from Xfoil and airfoil is provided check against xfoil
+        """
+        super(XfoilFit, self).process_result(result)
+
+        if self.mfac not in result["sensitivities"]["constants"]:
+            return
+        if amax([abs(result["sensitivities"]["constants"][self.mfac])]) < 1e-5:
+            return
+        if not self.airfoil:
+            return
+
+        from .xfoilWrapper import xfoil_comparison
+
+        cl, re = 0.0, 0.0
+        for dvar in self.dvars:
+            if "Re" in str(dvar):
+                re = result(dvar)
+            if "C_L" in str(dvar):
+                cl = result(dvar)
+        cd = result(self.ivar)
+        if not hasattr(cl, "__len__") and hasattr(re, "__len__"):
+            cl = [cl]*len(re)
+        err, cdx = xfoil_comparison(self.airfoil, cl, re, cd)
+        ind = where(err > 0.05)[0]
+        for i in ind:
+            msg = (
+                "Drag error for %s is %.2f. Re=%.1f; CL=%.4f;"
+                " Xfoil cd=%.6f, GP sol cd=%.6f"
+                % (
+                    ", ".join(self.ivar.descr["models"]),
+                    err[i],
+                    re[i],
+                    cl[i],
+                    cd[i],
+                    cdx[i],
+                )
+            )
+            print("Warning: %s" % msg)
diff --git a/gpfit/xfoil_wrapper.py b/gpfit/xfoil/wrapper.py
similarity index 73%
rename from gpfit/xfoil_wrapper.py
rename to gpfit/xfoil/wrapper.py
index 59a28ac..7fd661f 100644
--- a/gpfit/xfoil_wrapper.py
+++ b/gpfit/xfoil/wrapper.py
@@ -1,17 +1,18 @@
 " python wrapper to call xfoil "
-from __future__ import print_function
-from __future__ import division
 import subprocess
 import numpy as np
 
-#pylint: disable=invalid-name, bare-except, too-many-locals
-#pylint: disable=too-many-arguments
+
+# pylint: disable=invalid-name, bare-except, too-many-locals
+# pylint: disable=too-many-arguments
+
 
 def xfoil_comparison(airfoil, Cl, Re, Cd):
     """
-    comparison of XFOIL cd to input cd for given cl and re
-    INPUTS
-    ------
+    Comparison of XFOIL Cd to input Cd for given Cl and Re.
+
+    Arguments
+    ---------
     airfoil:                    airfoil name
                                     str - ("xxx.dat", "naca xxxx")
     Cl:                         lift coefficient
@@ -21,7 +22,7 @@ def xfoil_comparison(airfoil, Cl, Re, Cd):
     Cd:                         drag coefficient
                                     float or 1d-list or array
 
-    OUTPUTS
+    Returns
     -------
     err:                        error between Cd and XFOIL computed Cd
                                     1-d numpy array
@@ -58,12 +59,13 @@ def xfoil_comparison(airfoil, Cl, Re, Cd):
 
     return np.array(err), np.array(cdxs)
 
-def single_call(topline, cl, Re, M, max_iter=100,
-                pathname="/usr/local/bin/xfoil"):
+
+def single_call(topline, cl, Re, M, max_iter=100, pathname="/usr/local/bin/xfoil"):
     """
-    single XFOIL call for given cl, re and mach number
-    INPUTS
-    ------
+    Single XFOIL call for given Cl, Re and Mach number.
+
+    Arguments
+    ---------
     topline:                        load airfoil call in XFOIL
                                         str - e.g. "load xxx.dat"
     cl:                             lift coefficient
@@ -76,7 +78,7 @@ def single_call(topline, cl, Re, M, max_iter=100,
                                         int: default is 100
     pathname:                       system path to XFOIL
                                         str
-    OUTPUTS
+    Returns
     -------
     cd:                         XFOIL drag coefficient
                                     float
@@ -88,38 +90,39 @@ def single_call(topline, cl, Re, M, max_iter=100,
                                     str
     """
 
-    proc = subprocess.Popen([pathname], stdout=subprocess.PIPE,
-                            stdin=subprocess.PIPE)
-    proc.stdin.write(topline +
-                     'oper \n' +
-                     "iter %d\n" %(max_iter)+
-                     'visc \n' +
-                     "%.2e \n" %(Re) +
-                     "M \n" +
-                     "%.2f \n" %(M) +
-                     "a 2.0 \n" +
-                     "cl %.4f \n" %(cl) +
-                     '\n' +
-                     'quit \n')
+    proc = subprocess.Popen([pathname], stdout=subprocess.PIPE, stdin=subprocess.PIPE)
+    proc.stdin.write(
+        topline
+        + "oper \n"
+        + "iter %d\n" % (max_iter)
+        + "visc \n"
+        + "%.2e \n" % (Re)
+        + "M \n"
+        + "%.2f \n" % (M)
+        + "a 2.0 \n"
+        + "cl %.4f \n" % (cl)
+        + "\n"
+        + "quit \n"
+    )
     stdout_val = proc.communicate()[0]
     proc.stdin.close()
 
-    if ("VISCAL:  Convergence failed\n" in stdout_val):
+    if "VISCAL:  Convergence failed\n" in stdout_val:
         return stdout_val
 
     res = {}
     ostr = stdout_val.split()
     ctr = 0
     for i in range(0, len(ostr)):
-        ix = len(ostr)-(i+1)
+        ix = len(ostr) - (i + 1)
         vl = ostr[ix]
-        if vl in ['a', 'CL', 'CD', 'Cm']:
+        if vl in ["a", "CL", "CD", "Cm"]:
             res[vl] = ostr[ix + 2]
             ctr += 1
         if ctr >= 4:
             break
-    cd = res['CD']
-    cl = res['CL']
+    cd = res["CD"]
+    cl = res["CL"]
     # alpha_ret = res['a']
-    cm = res['Cm']
+    cm = res["Cm"]
     return float(cd), float(cl), float(cm), stdout_val
diff --git a/legacy_m_code/LM.m b/legacy_m_code/LM.m
deleted file mode 100755
index 61157bb..0000000
--- a/legacy_m_code/LM.m
+++ /dev/null
@@ -1,166 +0,0 @@
-function [params, RMStraj] = LM(residfun, initparams, varargin)
-
-%Levenberg-Marquardt algorithm
-%minimizes sum of squared error of residfun(params)
-%INPUTS
-%residfun should have the form [r, drdp] = residfun(params)
-%for curve fitting,
-% if residfun is (ydata - y(params)), drdp is -dydp
-% if residfun is (y(params) - ydata), drdp is dydp
-%params is a column vector of initial param guesses
-%varargin is a list of option, value pairs
-%OUTPUTS
-%params: best params found
-%RMStraj: history of RMS errors after each step (first point is initialization)
-
-tic;
-%check incoming params
-[nparam, sb1] = size(initparams);
-if sb1 > 1
-    error('params should be a column vector');
-end
-
-%set defaults; incorporate incoming options
-defaults.bverbose = true;
-%defaults.bplot = true;
-defaults.lambdainit = 0.02;
-defaults.maxiter = 5000;
-defaults.maxtime = 5;
-defaults.tolgrad = sqrt(eps);
-defaults.tolrms = 1e-7;
-options = process_options(defaults, varargin{:});
-
-%define display formatting if required
-formatstr1 = ' %5.0f   %9.6g    %9.3g\n';
-formatstr =  ' %5.0f   %9.6g    %9.3g %12.4g  %12.4g   %8.4g\n';
-
-%get residual values and jacobian at initial point; extract size info
-params = initparams;
-params_updated = true;
-[r, J] = residfun(params);
-[npt, sb1] = size(r);
-if sb1 > 1
-    error('residfun should return a column vector');
-elseif any(size(J) ~= [npt, nparam])
-    error('Jacobian size inconsistent');
-end
-
-%"accept" initial point
-rms = norm(r, 2)/sqrt(npt);
-maxgrad = norm(r'*J, Inf);
-prev_trial_accepted = false;
-
-%initializations
-iter = 1;
-Jissparse = issparse(J);
-diagJJ = sum(J.*J, 1)';
-zeropad = zeros(nparam,1);
-lambda = options.lambdainit;
-RMStraj = zeros(options.maxiter,1);
-RMStraj(1) = rms;
-gradcutoff = options.tolgrad;  %*max(maxgrad, 1); 
-
-%display info for 1st iter
-if(options.bverbose)
-fprintf( ...
-    ['\n                      First-Order                   Norm of \n', ...
-        '   Iter    Residual    optimality      Lambda        step     Jwarp\n']);
-        fprintf(formatstr1,iter,rms,maxgrad);
-end
-
-%main loop
-while true
-    
-    %test for exit cases
-    if iter == options.maxiter
-        if options.bverbose, disp('Reached maxiter'); end
-        break;
-    elseif toc > options.maxtime
-        if options.bverbose, disp(['Reached maxtime (', num2str(options.maxtime), ' seconds)']); end
-        break;
-    elseif iter > 2 && abs(RMStraj(iter)-RMStraj(iter-2)) < RMStraj(iter)*options.tolrms
-        %should really only allow this exit case if trust region
-        %constraint is slack
-        if options.bverbose, disp(['RMS changed less than tolrms']); end
-        break;
-    end
-    
-    iter = iter + 1;
-    
-    %compute diagonal scaling matrix based on current lambda and J
-    %note this matrix changes every iteration, since either lambda or J
-    %(or both) change on every iteration
-    if Jissparse
-        D = spdiags(sqrt(lambda*diagJJ), 0, nparam, nparam);
-    else
-        D = diag(sqrt(lambda*diagJJ));
-    end
-    
-    %update augmented least squares system
-    if params_updated
-        %J = trialJ;
-        diagJJ = sum(J.*J, 1)';
-        augJ = [J; D];
-        augr = [-r; zeropad];
-    else
-        augJ(npt+1:end, :) = D;
-    end
-    
-    %compute step for this lambda
-    step = augJ\augr;
-    trialp = params + step;
-    
-    %check function value at trialp
-    [trialr, trialJ] = residfun(trialp);
-    trialrms = norm(trialr)/sqrt(npt);
-    RMStraj(iter) = trialrms;
-    
-    %accept or reject trial params
-    if trialrms < rms
-        params = trialp;
-        J = trialJ;
-        r = trialr;
-        rms = trialrms;
-        maxgrad = norm(r'*J, Inf);
-        dsp();  %dsp here so that all grad info is for updated point, but lambda not yet updated
-        if maxgrad < gradcutoff
-            if options.bverbose, disp('1st order optimality attained'); end
-            break;
-        end
-        if prev_trial_accepted && iter > 1
-            lambda = lambda/10;
-        end
-        prev_trial_accepted = true;
-        params_updated = true;
-    else
-        dsp();
-        lambda = lambda*10;
-        prev_trial_accepted = false;
-        params_updated = false;
-    end
-    
-end
-
-RMStraj(iter+1:end) = [];
-disp(['Final RMS: ', num2str(rms)]);
-% if options.bplot
-%     fs = 14;
-%     figure(1); clf(1);
-%     semilogy(RMStraj);
-%     xlabel('iteration', 'fontsize', fs);
-%     ylabel('RMS error', 'fontsize', fs);
-% end
-
-    function dsp()
-        %display current progress
-        %note that dsp should be called immediately after evaluating a
-        %trial point (and before updating lambda) in order for lambda to
-        %correspond to the trust region used in the current iterate
-        if(options.bverbose)
-            fprintf(formatstr,iter,trialrms,maxgrad,lambda,norm(step),full(max(diagJJ)/min(diagJJ)));
-            %str = ['i=', num2str(iter), ' | bestrms=', num2str(rms), ' | lambda=', num2str(lambda), ' | maxg=', num2str(norm(r'*J, 'inf'))];
-            %disp(str);
-        end
-    end
-
-end
diff --git a/legacy_m_code/LMTR.m b/legacy_m_code/LMTR.m
deleted file mode 100755
index 9d173de..0000000
--- a/legacy_m_code/LMTR.m
+++ /dev/null
@@ -1,173 +0,0 @@
-function [params, RMStraj] = LMTR(residfun, initparams, varargin)
-
-%locally solve nonlinear least squares problem
-%using trust region algorithm motivated by Levenberg-Marquardt
-
-%INPUTS
-%residfun should have the form [r, drdp] = residfun(params)
-%for curve fitting,
-% if residfun is (ydata - y(params)), drdp is -dydp
-% if residfun is (y(params) - ydata), drdp is dydp
-%params is a column vector of initial param guesses
-%varargin is a list of option, value pairs
-
-%OUTPUTS
-%params: best params found
-%RMStraj: history of RMS errors after each step (first point is initialization)
-
-tic;
-%check incoming params
-[nparam, sb1] = size(initparams);
-if sb1 > 1
-    error('params should be a column vector');
-end
-
-%set defaults; incorporate incoming options
-defaults.bverbose = true;
-defaults.maxiter = 80;
-defaults.maxtime = inf;
-defaults.initdelta = 1;
-defaults.rhomin = 0.1;
-defaults.deltamax = Inf;
-options = process_options(defaults, varargin{:});
-
-if options.bverbose
-    %define display formatting if required
-    formatstr1 = ' %5.0f   %9.6g    %9.3g\n';
-    formatstr =  ' %5.0f   %9.6g    %9.3g    %8.4g %12.4g  %12.4g   %8.4g   %8.4g\n';
-end
-
-%get residual values and size at initial point; extract size info
-params = initparams;
-[r, J] = residfun(params);
-[npt, sb1] = size(r);
-if sb1 > 1
-    error('residfun should return a column vector');
-elseif any(size(J) ~= [npt, nparam])
-    error('Jacobian size inconsistent');
-end
-
-%"accept" initial point
-rr = r'*r;
-rms = sqrt(rr/npt);
-point_accepted = true;  %accept initial point so model gets updated
-maxgrad = 2*norm(J'*r, inf);
-
-%initializations
-iter = 1;   %number of function eval trials
-naccept = 0;  %number of accepted steps
-RMStraj = zeros(options.maxiter, 1);
-RMStraj(1) = rms;
-delta = options.initdelta;
-
-%trust region subproblem: initialize terms that don't change iter to iter
-prob.qcsubk = ones(nparam,1);
-prob.qcsubi = (1:nparam)';
-prob.qcsubj = (1:nparam)';
-
-if options.bverbose
-    %display info for 1st iter
-    fprintf( ...
-        ['\n                      First-Order                            Norm of \n', ...
-        '   Iter    Residual    optimality     rho        Delta        step     Jwarp        alpha\n']);
-    fprintf(formatstr1,iter,rms,maxgrad);
-end
-
-while true
-    %test for exit cases
-    if iter == options.maxiter
-        if options.bverbose, disp('Reached maxiter'); end
-        break;
-    elseif toc > options.maxtime
-        if options.bverbose, disp(['Reached maxtime (', num2str(options.maxtime), ' seconds)']); end
-        break;
-    end
-    
-    %if point accepted, update model, check gradient for exit case
-    if(point_accepted)
-        %model is that r = r0+J*step, so rr ~ rr0 + step'*JJ*step + 2J'*step
-        %we actually optimize 1/2*step'JJstep + J'*step 
-        %mosek includes 1/2 in quadratic terms automatically
-        % optimal model objective *2 is predicted decrease in r'r
-        prob.c = J'*r;  % == -g, the gradient of objective wrt delta
-        %check gradient for exit case here
-        JJ = J'*J;
-        D = diag(JJ);
-        
-        %set up trust region subproblem
-        %form
-        [prob.qosubi, prob.qosubj, prob.qoval] = find(sparse(tril(JJ)));
-        prob.a = sparse(1,size(JJ,2));  %no linear terms in constraints
-        %can add bound constraints on step *and* params using blx, bux
-        
-        %elliptical step bound (2-norm)
-        prob.qcval = D; 
-        
-    end
-    
-    iter = iter + 1;
-    
-    %solve trust region subproblem as a SOCP
-    prob.buc = 1/2*delta;   %mosek assumes 1/2 in objective, so add here too
-    [return_code, res] = mosekopt('minimize echo(0)', prob);
-    if return_code ~= 0
-        error(['mosek returned return code ', res.rcodestr]);
-    end
-    
-    
-    step = res.sol.itr.xx;
-    trialp = params + step;
-    [trialr, trialJ] = residfun(trialp);
-    trialrr = trialr'*trialr;
-    rho = (rr - trialrr)/(-2*res.sol.itr.pobjval);
-    RMStraj(iter) = sqrt(trialrr/npt);
-    
-    %next bit of logic follows NumericalOptimization algorithm 4.1
-    if rho > options.rhomin     %also represented by 'eta'
-        %accept step
-        params = trialp;
-        J = trialJ;
-        r = trialr;
-        rr = trialrr;
-        rms = sqrt(rr/npt);
-        maxgrad = 2*norm(J'*r, inf);
-        point_accepted = true;
-    else
-        point_accepted = false;
-    end
-    
-    %display info
-    dsp();
-        
-    if rho < .25    %actual decrease too small; shrink trust region
-        delta = delta/4;
-    else %step decreased residual appreciably
-        if rho > .75 && abs(sum(D.*step.^2)/delta - 1) < 1e-6
-            %if model in good agreement and trust region constrained step
-            delta = min(2*delta, options.deltamax);
-            %else leave delta alone
-        end
-    end
-        
-    
-    %be sure to update point_accepted based on trial
-    
-end
-
-RMStraj(iter+1:end) = [];
-
-
-    function dsp()
-        %display current progress
-        %note that dsp should be called immediately after evaluating a
-        %trial point (and before updating delta) in order for delta to
-        %correspond to the trust region used in the current iterate
-        if(options.bverbose)
-            fprintf(formatstr,iter,RMStraj(iter),maxgrad,rho,delta,norm(step),cond(full(JJ)), 1/params(end));
-            %str = ['i=', num2str(iter), ' | bestrms=', num2str(rms), ' | lambda=', num2str(lambda), ' | maxg=', num2str(norm(r'*J, 'inf'))];
-            %disp(str);
-        end
-    end
-
-
-end
\ No newline at end of file
diff --git a/legacy_m_code/albatross.m b/legacy_m_code/albatross.m
deleted file mode 100644
index 1ef0f0e..0000000
--- a/legacy_m_code/albatross.m
+++ /dev/null
@@ -1,47 +0,0 @@
-[X, Y] = meshgrid(1:5,1:5);
-Z = X.^2+Y.^2;
-% surf(X,Y,Z);
-% hold all;
-w = Z(:);
-X = X(:);
-Y = Y(:);
-u = [X,Y];
-
-x = log(u);
-y = log(w);
-
-s = compare_fits(x,y, 2,1);
-
-[u1, u2] = meshgrid(1:5,1:5);
-
-% Max Affine Fitting
-PAR_MA = s.maxaffine.params{1};
-A = PAR_MA([2,3,5,6]);
-B = PAR_MA([1,4]);
-w_MA = exp(B(1)) .* u1.^A(1) .* u2.^A(2);
-% surf(u1,u2,w_MA);
-% hold all;
-w_MA = exp(B(2)) .* u1.^A(3) .* u2.^A(4);
-% surf(u1,u2,w_MA);
-% hold all;
-
-% Softmax Affine Fitting
-PAR_SMA = s.softmax_optMAinit.params{1};
-A = PAR_SMA([2,3,5,6]);
-B = PAR_SMA([1,4]);
-alpha = PAR_SMA(end);
-w_SMA = (exp(alpha.*B(1)) .* u1.^(alpha*A(1)) .* u2.^(alpha*A(2)) +...
-        exp(alpha.*B(2)) .* u1.^(alpha*A(3)) .* u2.^(alpha*A(4))...
-        ).^(1/alpha);
-% surf(u1,u2,w_SMA);
-
-% Implicit Softmax Affine Fitting
-PAR_ISMA = s.implicit_originit.params{1};
-% alpha = PAR_SMA(end-1:end);
-% A = PAR_SMA([2,3,5,6]);
-% B = PAR_SMA([1,4]);
-% f_ISMA = @(w_ISMA) exp(alpha(1).*B(1))/w_ISMA.^alpha(1) .* uu1.^(alpha(1)*A(1)) .* uu2.^(alpha(1)*A(2)) +...
-%     exp(alpha(2).*B(2))/w_ISMA.^alpha(2) .* uu1.^(alpha(2)*A(3)) .* uu2.^(alpha(2)*A(4)) - 1;
-w_ISMA = implicit_softmax_affine(x,PAR_ISMA);
-w_ISMA = exp(reshape(w_ISMA,5,5))
-% surf(u1,u2,w_ISMA);
\ No newline at end of file
diff --git a/legacy_m_code/alligator.m b/legacy_m_code/alligator.m
deleted file mode 100644
index c904e29..0000000
--- a/legacy_m_code/alligator.m
+++ /dev/null
@@ -1,37 +0,0 @@
-[X, Y] = meshgrid(1:5,1:5);
-Z = (X.^2+Y.^2).^(1/2);
-surf(X,Y,Z);
-hold all;
-w = Z(:);
-X = X(:);
-Y = Y(:);
-
-u = [X,Y];
-
-x = log(u);
-y = log(w);
-
-s = compare_fits(x,y, 2,1);
-
-[u1, u2] = meshgrid(1:5,1:5);
-
-% Max Affine Fitting
-PAR_MA = s.maxaffine.params{1};
-A = PAR_MA([2,3,5,6]);
-B = PAR_MA([1,4]);
-w_MA = exp(B(1)) .* u1.^A(1) .* u2.^A(2);
-% surf(u1,u2,w_MA);
-% hold all;
-w_MA = exp(B(2)) .* u1.^A(3) .* u2.^A(4);
-% surf(u1,u2,w_MA);
-% hold all;
-
-% Softmax Affine Fitting
-PAR_SMA = s.softmax_optMAinit.params{1};
-A = PAR_SMA([2,3,5,6]);
-B = PAR_SMA([1,4]);
-alpha = 1/PAR_SMA(end);
-w_SMA = (exp(alpha.*B(1)) .* u1.^(alpha*A(1)) .* u2.^(alpha*A(2)) +...
-        exp(alpha.*B(2)) .* u1.^(alpha*A(3)) .* u2.^(alpha*A(4))...
-        ).^(1/alpha);
-% surf(u1,u2,w_SMA);
\ No newline at end of file
diff --git a/legacy_m_code/angelfish.m b/legacy_m_code/angelfish.m
deleted file mode 100644
index 98434a7..0000000
--- a/legacy_m_code/angelfish.m
+++ /dev/null
@@ -1,38 +0,0 @@
-[X, Y] = meshgrid(1:0.1:5,1:0.1:5);
-Z = X./Y + Y./X;
-surf(X,Y,Z);
-
-hold all;
-w = Z(:);
-X = X(:);
-Y = Y(:);
-
-u = [X,Y];
-
-x = log(u);
-y = log(w);
-
-s = compare_fits(x,y, 2,1);
-
-[u1, u2] = meshgrid(1:5,1:5);
-
-% Max Affine Fitting
-PAR_MA = s.maxaffine.params{1};
-A = PAR_MA([2,3,5,6]);
-B = PAR_MA([1,4]);
-w_MA = exp(B(1)) .* u1.^A(1) .* u2.^A(2);
-% surf(u1,u2,w_MA);
-% hold all;
-w_MA = exp(B(2)) .* u1.^A(3) .* u2.^A(4);
-% surf(u1,u2,w_MA);
-% hold all;
-
-% Softmax Affine Fitting
-PAR_SMA = s.softmax_optMAinit.params{1};
-A = PAR_SMA([2,3,5,6]);
-B = PAR_SMA([1,4]);
-alpha = 1/PAR_SMA(end);
-w_SMA = (exp(alpha.*B(1)) .* u1.^(alpha*A(1)) .* u2.^(alpha*A(2)) +...
-        exp(alpha.*B(2)) .* u1.^(alpha*A(3)) .* u2.^(alpha*A(4))...
-        ).^(1/alpha);
-surf(u1,u2,w_SMA);
\ No newline at end of file
diff --git a/legacy_m_code/anteater.m b/legacy_m_code/anteater.m
deleted file mode 100644
index 8e07da6..0000000
--- a/legacy_m_code/anteater.m
+++ /dev/null
@@ -1,38 +0,0 @@
-[X, Y] = meshgrid(1:0.1:5,1:0.1:5);
-Z = (X./Y + Y./X).^0.5;
-surf(X,Y,Z);
-
-hold all;
-w = Z(:);
-X = X(:);
-Y = Y(:);
-
-u = [X,Y];
-
-x = log(u);
-y = log(w);
-
-s = compare_fits(x,y, 2,1);
-
-[u1, u2] = meshgrid(1:5,1:5);
-
-% Max Affine Fitting
-PAR_MA = s.maxaffine.params{1};
-A = PAR_MA([2,3,5,6]);
-B = PAR_MA([1,4]);
-w_MA = exp(B(1)) .* u1.^A(1) .* u2.^A(2);
-surf(u1,u2,w_MA);
-hold all;
-w_MA = exp(B(2)) .* u1.^A(3) .* u2.^A(4);
-surf(u1,u2,w_MA);
-hold all;
-
-% Softmax Affine Fitting
-% % PAR_SMA = s.softmax_optMAinit.params{1};
-% % A = PAR_SMA([2,3,5,6]);
-% % B = PAR_SMA([1,4]);
-% % alpha = 1/PAR_SMA(end);
-% % w_SMA = (exp(alpha.*B(1)) .* u1.^(alpha*A(1)) .* u2.^(alpha*A(2)) +...
-% %         exp(alpha.*B(2)) .* u1.^(alpha*A(3)) .* u2.^(alpha*A(4))...
-% %         ).^(1/alpha);
-% % surf(u1,u2,w_SMA);
\ No newline at end of file
diff --git a/legacy_m_code/antelope.m b/legacy_m_code/antelope.m
deleted file mode 100644
index 8db8317..0000000
--- a/legacy_m_code/antelope.m
+++ /dev/null
@@ -1,39 +0,0 @@
-[X, Y] = meshgrid(1:0.01:2,0.2:0.002:0.4);
-Z = (1.09*X.^4.27.*Y.^0.112 + (7.79e-5)*X.^4.75./Y.^6.44 + (1.36e-7)*X.^8.94./Y.^8.89).^(1/2.14);
-surf(X,Y,Z);
-
-hold all;
-w = Z(:);
-X = X(:);
-Y = Y(:);
-
-u = [X,Y];
-
-x = log(u);
-y = log(w);
-
-s = compare_fits(x,y, 2,1);
-
-[u1, u2] = meshgrid(1:0.1:2,0.2:0.02:0.4);
-
-% Max Affine Fitting
-% PAR_MA = s.maxaffine.params{1};
-% A = PAR_MA([2,3,5,6]);
-% B = PAR_MA([1,4]);
-% w_MA = exp(B(1)) .* u1.^A(1) .* u2.^A(2);
-% surf(u1,u2,w_MA);
-% hold all;
-% w_MA = exp(B(2)) .* u1.^A(3) .* u2.^A(4);
-% surf(u1,u2,w_MA);
-% hold all;
-
-% Softmax Affine Fitting
-PAR_SMA = s.softmax_optMAinit.params{1};
-A = PAR_SMA([2,3,5,6]);
-B = PAR_SMA([1,4]);
-alpha = 1/PAR_SMA(end);
-% alpha = PAR_SMA(end);
-w_SMA = (exp(alpha.*B(1)) .* u1.^(alpha*A(1)) .* u2.^(alpha*A(2)) +...
-        exp(alpha.*B(2)) .* u1.^(alpha*A(3)) .* u2.^(alpha*A(4))...
-        ).^(1/alpha);
-surf(u1,u2,w_SMA);
\ No newline at end of file
diff --git a/legacy_m_code/appenzeller.m b/legacy_m_code/appenzeller.m
deleted file mode 100644
index 4d67c06..0000000
--- a/legacy_m_code/appenzeller.m
+++ /dev/null
@@ -1,39 +0,0 @@
-[VDD, VTH] = meshgrid(1:0.01:2,0.2:0.002:0.4);
-Z = VDD.^2 + 30.*VDD.*exp(-(VTH-0.06.*VDD)/0.039);
-surf(X,Y,Z);
-
-hold all;
-w = Z(:);
-X = X(:);
-Y = Y(:);
-
-u = [X,Y];
-
-x = log(u);
-y = log(w);
-
-s = compare_fits(x,y, 2,1);
-
-[u1, u2] = meshgrid(1:0.1:2,0.2:0.02:0.4);
-
-% Max Affine Fitting
-% PAR_MA = s.maxaffine.params{1};
-% A = PAR_MA([2,3,5,6]);
-% B = PAR_MA([1,4]);
-% w_MA = exp(B(1)) .* u1.^A(1) .* u2.^A(2);
-% surf(u1,u2,w_MA);
-% hold all;
-% w_MA = exp(B(2)) .* u1.^A(3) .* u2.^A(4);
-% surf(u1,u2,w_MA);
-% hold all;
-
-% Softmax Affine Fitting
-PAR_SMA = s.softmax_optMAinit.params{1};
-A = PAR_SMA([2,3,5,6]);
-B = PAR_SMA([1,4]);
-alpha = 1/PAR_SMA(end);
-% alpha = PAR_SMA(end);
-w_SMA = (exp(alpha.*B(1)) .* u1.^(alpha*A(1)) .* u2.^(alpha*A(2)) +...
-        exp(alpha.*B(2)) .* u1.^(alpha*A(3)) .* u2.^(alpha*A(4))...
-        ).^(1/alpha);
-surf(u1,u2,w_SMA);
\ No newline at end of file
diff --git a/legacy_m_code/armadillo.m b/legacy_m_code/armadillo.m
deleted file mode 100644
index e9b2b4a..0000000
--- a/legacy_m_code/armadillo.m
+++ /dev/null
@@ -1,40 +0,0 @@
-[X, Y] = meshgrid(1:0.01:2,0.2:0.002:0.4);
-Z = X.^2 + 30.*X.*exp(-(Y-0.06.*X)/0.039);
-surf(X,Y,Z);
-
-hold all;
-w = Z(:);
-X = X(:);
-Y = Y(:);
-
-u = [X,Y];
-
-x = log(u);
-y = log(w);
-
-s = compare_fits(x,y, 3,1);
-
-[u1, u2] = meshgrid(1:0.1:2,0.2:0.02:0.4);
-
-% Max Affine Fitting
-PAR_MA = s.maxaffine.params{1}
-A = PAR_MA([2,3,5,6])
-B = PAR_MA([1,4])
-w_MA_1 = exp(B(1)) .* u1.^A(1) .* u2.^A(2);
-% surf(u1,u2,w_MA);
-% hold all;
-w_MA_2 = exp(B(2)) .* u1.^A(3) .* u2.^A(4);
-% surf(u1,u2,w_MA);
-% hold all;
-
-% Softmax Affine Fitting
-PAR_SMA = s.softmax_optMAinit.params{1};
-A = PAR_SMA([2,3,5,6,8,9]);
-B = PAR_SMA([1,4,7]);
-alpha = 1/PAR_SMA(end);
-% alpha = PAR_SMA(end);
-w_SMA = (exp(alpha.*B(1)) .* u1.^(alpha*A(1)) .* u2.^(alpha*A(2)) +...
-        exp(alpha.*B(2)) .* u1.^(alpha*A(3)) .* u2.^(alpha*A(4)) +...
-        exp(alpha.*B(3)) .* u1.^(alpha*A(5)) .* u2.^(alpha*A(6))...
-        ).^(1/alpha);
-surf(u1,u2,w_SMA);
\ No newline at end of file
diff --git a/legacy_m_code/compare_fits.m b/legacy_m_code/compare_fits.m
deleted file mode 100755
index 5f690ee..0000000
--- a/legacy_m_code/compare_fits.m
+++ /dev/null
@@ -1,46 +0,0 @@
-function s = compare_fits(xdata, ydata, Ks, ntry)
-
-[npt, dimx] = size(xdata);
-s = struct();
-s.Ks = Ks;
-alphainit = 10;
-
-for ik = 1:length(Ks)
-    for t = 1:ntry
-        k = Ks(ik);
-        
-        bainit = max_affine_init(xdata, ydata, k);
-        
-        rfun = @(p) generic_resid_fun(@max_affine, xdata, ydata, p);
-        tic;
-        [params, RMStraj] = LM(rfun, [bainit(:)]);
-        store_results('maxaffine');
-        
-        rfun = @(p) generic_resid_fun(@softmax_affine, xdata, ydata, p);
-        tic;
-        [params, RMStraj] = LM(rfun, [params(:); alphainit]);
-        store_results('softmax_optMAinit');
-        
-        rfun = @(p) generic_resid_fun(@softmax_affine, xdata, ydata, p);
-        tic;
-        [params, RMStraj] = LM(rfun, [bainit(:); alphainit]);
-        store_results('softmax_originit');
-
-        rfun = @(p) generic_resid_fun(@implicit_softmax_affine, xdata, ydata, p);
-        tic;
-        [params, RMStraj] = LM(rfun, [bainit(:); alphainit*ones(k,1)]);
-        store_results('implicit_originit');
-        
-    end
-end
-
-
-    function store_results(fieldname)
-        s.(fieldname).resid(ik,t) = min(RMStraj);
-        s.(fieldname).iter(ik,t) = length(RMStraj);
-        s.(fieldname).params{ik,t} = params;
-        s.(fieldname).time(ik,t) = toc;
-    end
-
-
-end
\ No newline at end of file
diff --git a/legacy_m_code/ex6_1.m b/legacy_m_code/ex6_1.m
deleted file mode 100644
index 4835199..0000000
--- a/legacy_m_code/ex6_1.m
+++ /dev/null
@@ -1,11 +0,0 @@
-m = 501;
-
-u = logspace(0,log10(3),501);
-u=u';
-w = (u.^2 + 3)./(u+1).^2;
-
-x = log(u);
-y = log(w);
-
-s = compare_fits(x, y, 3, 1);
-
diff --git a/legacy_m_code/generic_resid_fun.m b/legacy_m_code/generic_resid_fun.m
deleted file mode 100755
index 0072180..0000000
--- a/legacy_m_code/generic_resid_fun.m
+++ /dev/null
@@ -1,14 +0,0 @@
-function [r, drdp] = generic_resid_fun(yfun, xdata, ydata, params)
-    %generic residual function -- converts yfun(xdata,params)
-    %to a residfun [r, drdp] = generic_resid_fun(yfun, xdata, ydata, params)
-    %used by nonlinear least squares fitting algorithms
-    %to get a residual function [r,drdp] = residfun(params),
-    %use rfun = @(p) generic_resid_fun(@yfun, xdata, ydata, p)
-    
-    %note this function defines resids as + when yhat > ydata 
-    %(opposite of typical conventions, but eliminates need for sign change)
-
-    [yhat, drdp] = yfun(xdata, params);
-    r = yhat - ydata;
-    
-end
\ No newline at end of file
diff --git a/legacy_m_code/implicit_softmax_affine.m b/legacy_m_code/implicit_softmax_affine.m
deleted file mode 100755
index 50632bc..0000000
--- a/legacy_m_code/implicit_softmax_affine.m
+++ /dev/null
@@ -1,47 +0,0 @@
-function [y, dydp] = implicit_softmax_affine(x, params)
-%%params may come in as
-%  1) a cell {ba, alpha}, or 
-%  2) a vector (with alpha last)
-%after reshaping (column-major), ba is dimx+1 by K
-%first row is b
-%rest is a
-
-[npt, dimx] = size(x);
-
-if iscell(params)
-    ba = params{1};
-    alpha = params{2};
-    K = numel(ba)/(dimx+1);
-else
-    K = numel(params)/(dimx+2);
-    ba = params(1:end-K);
-    alpha = params(end-K+1:end);
-end
-
-%reshape ba to matrix
-ba = reshape(ba, dimx+1, K);
-
-if any(alpha <= 0)
-    y = Inf*ones(npt,1);
-    dydp = nan;
-    return;
-end
-
-%augment data with column of ones
-X = [ones(npt,1), x];
-
-%compute affine functions
-z = X*ba;
-
-if nargout == 1
-    y = lse_implicit(z, alpha);
-else
-    [y, dydz, dydalpha] = lse_implicit(z, alpha);
-    dydba = repcols(dydz, dimx+1) .* repmat(X, 1, K);
-    dydp = [dydba, dydalpha];
-end
-
-end
-
-
-
diff --git a/legacy_m_code/lse_implicit.m b/legacy_m_code/lse_implicit.m
deleted file mode 100755
index 7beb4ff..0000000
--- a/legacy_m_code/lse_implicit.m
+++ /dev/null
@@ -1,71 +0,0 @@
-function [y, dydx, dydalpha] = lse_implicit(x, alpha)
-%implicit log sum exp function with derivatives
-%sums across the second dimension of x;
-%alpha should be a rox vector with length(alpha)==size(x,2)
-%returns one y for every row of x
-%dydx gives deriv of each y wrt each x
-%alpha is local softness parameter for each column of x
-
-%lse_implicit is a mapping R^n --> R, where n==size(x,2)
-%implementation: newton raphson steps to find f(x,y) = 0
-
-tol = 10*eps;
-[npt, nx] = size(x);
-
-if all(size(alpha) == [nx, 1])
-    alpha = alpha';
-end
-if ~all(size(alpha) == [1, nx])
-    error('alpha size mismatch');
-end
-alphamat = repmat(alpha, npt, 1);
-
-m = max(x, [], 2);  %maximal x values
-h = x - repmat(m, 1, nx);   %distance from m; note h <= 0 for all entries
-%(unless x has infs; in this case h has nans; can prob deal w/ this gracefully)
-%should also deal with alpha==inf case gracefully
-
-L = zeros(npt, 1);  %initial guess. note y = m + L
-Lmat = repmat(L, 1, nx);
-
-%initial eval
-expo = exp(alphamat.*(h-Lmat));
-alphaexpo = alphamat.*expo;
-sumexpo = sum(expo, 2);
-sumalphaexpo = sum(alphaexpo, 2);
-f = log(sumexpo);
-dfdL = -sum(alphaexpo, 2)./sumexpo;
-neval = 1;
-i = abs(f) > tol;   %inds to update
-%disp(['max newton-raphson residual: ', num2str(max(abs(f)))]);
-
-while any(i)
-    L(i) = L(i) - f(i)./dfdL(i);    %newton step
-    
-    %re-evaluate
-    Lmat(i,:) = repmat(L(i), 1, nx);
-    expo(i,:) = exp(alphamat(i,:).*(h(i,:)-Lmat(i,:)));
-    alphaexpo(i,:) = alphamat(i,:).*expo(i,:);
-    sumexpo(i) = sum(expo(i,:), 2);
-    sumalphaexpo(i,:) = sum(alphaexpo(i,:), 2);
-    f(i) = log(sumexpo(i));
-    dfdL(i) = -sumalphaexpo(i,:)./sumexpo(i);
-    neval = neval + 1;
-    if neval > 40
-        disp('');
-    end
-    
-    %update inds that need to be evaluated
-    i(i) = abs(f(i)) > tol;
-    %disp(['max newton-raphson residual: ', num2str(max(abs(f)))]);
-end
-
-disp(['lse_implicit converged in ', num2str(neval), ' newton-raphson steps']);
-
-y = m + L;
-
-if nargout > 1
-    dydx = alphaexpo./repmat(sumalphaexpo, 1, nx);
-    dydalpha = (h - Lmat).*expo./repmat(sumalphaexpo, 1, nx);
-end
-
diff --git a/legacy_m_code/lse_scaled.m b/legacy_m_code/lse_scaled.m
deleted file mode 100755
index 399d0a1..0000000
--- a/legacy_m_code/lse_scaled.m
+++ /dev/null
@@ -1,29 +0,0 @@
-function [y, dydx, dydalpha] = lse_scaled(x, alpha)
-%log sum exp function with derivatives
-%sums across the second dimension of x; 
-%returns one y for every row of x
-%dydx gives the deriv of each y wrt each x 
-%alpha is softness parameter \in R
-
-%note that lse_scaled is a mapping R^n --> R, where n==size(x,2)
-
-[~, n] = size(x);
-
-m = max(x, [], 2);  %maximal x values
-h = x - repmat(m, 1, n);   %distance from m; note h <= 0 for all entries
-%(unless x has infs; in this case h has nans; can prob deal w/ this gracefully)
-%should also deal with alpha==inf case gracefully
-
-expo = exp(alpha*h);
-sumexpo = sum(expo, 2);
-
-L = log(sumexpo)/alpha;
-y = L + m;
-
-if nargout > 1
-    dydx = expo./repmat(sumexpo, 1, n);
-    %note that sum(dydx,2)==1, i.e. dydx is a probability distribution
-    dydalpha = (sum(h.*expo, 2)./sumexpo - L)/alpha;
-end
-
-end
\ No newline at end of file
diff --git a/legacy_m_code/magnani_boyd.m b/legacy_m_code/magnani_boyd.m
deleted file mode 100755
index ed7f942..0000000
--- a/legacy_m_code/magnani_boyd.m
+++ /dev/null
@@ -1,58 +0,0 @@
-function [ba, residtraj] = magnani_boyd(x, y, K, varargin)
-
-%carries out the Magnani-Boyd Least Squares Partition Algorithm to fit a
-%convex max-affine model to data (x, y)
-
-defaults.maxiter = 100;
-defaults.relcutoff = 1e-6; %stop when abs(resid-oldresid)/oldresid < relcutoff
-defaults.bverbose = true;
-options = process_options(defaults, varargin{:});
-
-npt = size(x, 1);
-dimx = size(x,2);
-X = [ones(npt,1), x];
-residfun = @(yhat) norm(yhat - y, 2)/sqrt(npt);
-
-iter = 0;
-%oldresid = Inf; %used to keep track of changes in resid
-residtraj = zeros(options.maxiter, 1);
-while(iter < options.maxiter)
-    iter = iter + 1;
-    
-    %start by updating ba (or initializing it if first iter)
-    if(iter == 1)
-        ba = max_affine_init(x, y, K);
-    else
-        %fit each local partition with an affine plane
-        for k = 1:K
-            inds = partition==k;
-            if rank( X(inds, :) ) >= dimx+1
-                %only do the update for partitions whose rank is large
-                %enough
-                ba(:, k) = X(inds, :)\y(inds);  %MATLAB warning says size changes, but this is false.
-            end
-        end
-    end
-
-    %calculate yhat and new partitioning
-    [yhat, partition] = max(X*ba, [], 2);
-    residtraj(iter) = residfun(yhat);
-    
-    if options.bverbose
-        str = ['i=', num2str(iter), ': resid=', num2str(residtraj(iter))];
-        disp(str);
-    end
-    
-    if iter > 1 && abs(residtraj(iter)-residtraj(iter-1))/residtraj(iter) < options.relcutoff
-        if options.bverbose
-            disp('magnani_boyd: exiting on relative change in resid');
-        end
-        break;
-    end
-    
-end
-
-%update residtraj to cut off unused entries at end
-residtraj(iter+1:end) = [];
-
-end
\ No newline at end of file
diff --git a/legacy_m_code/max_affine.m b/legacy_m_code/max_affine.m
deleted file mode 100755
index 1f3b63a..0000000
--- a/legacy_m_code/max_affine.m
+++ /dev/null
@@ -1,33 +0,0 @@
-function [y, dydba] = max_affine(x,ba)
-%ba may come in as a matrix or as a vector
-%after reshaping (column-major), ba is dimx+1 by K
-%first row is b
-%rest is a
-
-[npt, dimx] = size(x);
-K = numel(ba)/(dimx+1);
-ba = reshape(ba,dimx+1,K);
-
-%augment data with column of ones
-X = [ones(npt,1), x];
-
-if(nargout == 1)
-    y = max(X*ba, [], 2);
-else
-    [y, partition] = max(X*ba, [], 2);
-    
-%   (non-sparse version)
-%     dydba = zeros(size(x,2)+1, K, npt);
-%     %todo - make this sparse
-%     for k = 1:K
-%         inds = partition == k;
-%         dydba(:,k,inds) = [ones(1,nnz(inds)); x(inds,:)'];
-%     end
-
-    dydba = spalloc(npt, (dimx+1)*K, npt*(dimx+1));
-    for k = 1:K
-        inds = partition==k;
-        indadd = (dimx+1)*(k-1);
-        dydba(inds, indadd+(1:dimx+1)) = X(inds,:);
-    end
-end
\ No newline at end of file
diff --git a/legacy_m_code/max_affine_init.m b/legacy_m_code/max_affine_init.m
deleted file mode 100755
index ac11827..0000000
--- a/legacy_m_code/max_affine_init.m
+++ /dev/null
@@ -1,63 +0,0 @@
-function ba = max_affine_init(x, y, K, varargin)
-%initializes max-affine fit to data (y, x)
-%ensures that initialization has at least K+1 points per partition (i.e.
-%per affine function)
-
-defaults.bverbose = true;
-options = process_options(defaults, varargin{:});
-
-npt = size(x,1);
-dimx = size(x,2);
-X = [ones(npt, 1), x];
-ba = zeros(dimx+1,K);
-if K*(dimx+1) > npt
-    error('Not enough data points');
-end
-
-%choose K unique indices
-randinds = randperm(npt, K);
-
-%partition based on distances
-sqdists = zeros(npt, K);
-for k = 1:K
-    sqdists(:,k) = sum((x - repmat(x(randinds(k),:), npt, 1)).^2, 2);
-end
-[~, mindistind] = min(sqdists, [], 2);  %index to closest k for each data pt
-
-%loop through each partition, making local fits
-%note we expand partitions that result in singular least squares problems
-%why this way? some points will be shared by multiple partitions, but
-%resulting max-affine fit will tend to be good. (as opposed to solving least-norm version)
-for k = 1:k
-    inds = mindistind == k;
-    
-    %before fitting, check rank and increase partition size if necessary
-    %(this does create overlaps)
-    if rank(X(inds, :)) < dimx + 1
-        [~, sortdistind] = sort(sqdists(:,k));
-        i = sum(inds);  %i is number of points in partition
-        iinit = i;
-        if i < dimx+1
-            %obviously, at least need dimx+1 points. fill these in before
-            %checking any ranks
-            inds(sortdistind(i+1:dimx+1)) = 1;
-            i = dimx+1;
-        end
-        %now add points until rank condition satisfied
-        while rank(X(inds, :)) < dimx+1
-            i = i+1;
-            inds(sortdistind(i)) = 1;
-        end
-        
-        if options.bverbose
-            disp(['max_affine_init: Added ', num2str(i-iinit), ' points to partition ', num2str(k), ' to maintain full rank for local fitting.']);
-        end
-        
-    end
-    
-    %now create the local fit
-    ba(:,k) = X(inds, :)\y(inds);
-    
-end
-
-end
\ No newline at end of file
diff --git a/legacy_m_code/process_options.m b/legacy_m_code/process_options.m
deleted file mode 100755
index 35935d8..0000000
--- a/legacy_m_code/process_options.m
+++ /dev/null
@@ -1,11 +0,0 @@
-function options = process_options(defaults, varargin)
-
-options = defaults;
-%process inputs
-i = 1;
-while(i < length(varargin))
-    options.(varargin{i}) = varargin{i+1};
-    i = i+2;
-end
-
-end
\ No newline at end of file
diff --git a/legacy_m_code/profile_drag_data.mat b/legacy_m_code/profile_drag_data.mat
deleted file mode 100755
index d63e9d1..0000000
Binary files a/legacy_m_code/profile_drag_data.mat and /dev/null differ
diff --git a/legacy_m_code/qcqp_scratch.m b/legacy_m_code/qcqp_scratch.m
deleted file mode 100755
index 7bd7c1f..0000000
--- a/legacy_m_code/qcqp_scratch.m
+++ /dev/null
@@ -1,35 +0,0 @@
-% cvx_begin sdp quiet
-% variable t
-% variable x(4,1)
-% minimize t
-% subject to
-% [speye(10), A*x-b; (A*x-b)' t] > 0
-% [speye(4), x; x', 1] > 0
-% cvx_end
-
-% cvx_begin quiet
-% variable x(4,1)
-% minimize norm(A*x-b)
-% subject to 
-% norm(x, 2) < 1
-% cvx_end
-
-clear prob;
-
-AA = A'*A;
-prob.c = -A'*b;
-[prob.qosubi, prob.qosubj, prob.qoval] = find(sparse(tril(AA)));
-prob.a = sparse(1,size(AA,2));
-prob.buc = 1/2;
-
-% prob.cones = cell(1,1);
-% prob.cones{1}.type = 'MSK_CT_QUAD';
-% prob.cones{1}.sub = [size(AA,2)+1, 1:size(AA,2)];
-
-%quadratic constraints on x
-prob.qcsubk = ones(size(AA,2),1);
-prob.qcsubi = (1:size(AA,2))';
-prob.qcsubj = (1:size(AA,2))';
-prob.qcval = diag(AA);
-
-[r, res] = mosekopt('minimize echo(0)', prob);
\ No newline at end of file
diff --git a/legacy_m_code/repcols.m b/legacy_m_code/repcols.m
deleted file mode 100755
index 89e36f4..0000000
--- a/legacy_m_code/repcols.m
+++ /dev/null
@@ -1,9 +0,0 @@
-function matout = repcols(matin, n)
-%replicate columns of a matrix
-%returns a matrix with n times as many columns as matin
-%example: if matin is [a b c], repcols(matin, 2) returns [a a b b c c].
-[nrow, ncol] = size(matin);
-
-matout = reshape(repmat(matin, n, 1), nrow, ncol*n);
-
-end
\ No newline at end of file
diff --git a/legacy_m_code/softmax_affine.m b/legacy_m_code/softmax_affine.m
deleted file mode 100755
index 7c5d1e2..0000000
--- a/legacy_m_code/softmax_affine.m
+++ /dev/null
@@ -1,56 +0,0 @@
-function [y, dydp] = softmax_affine(x,params,softness_param_is_alpha)
-%params may come in as
-%  1) a cell {ba, alpha}, or 
-%  2) a vector (with alpha last)
-%after reshaping (column-major), ba is dimx+1 by K
-%first row is b
-%rest is a
-
-if nargin == 2
-    softness_param_is_alpha = false;    %then it's gamma (gamma = 1/alpha)
-end
-
-if iscell(params)
-    ba = params{1};
-    softness = params{2};
-else
-    ba = params(1:end-1);
-    softness = params(end);
-end
-
-if softness_param_is_alpha
-    alpha = softness;
-else
-    alpha = 1/softness;
-end
-
-%check sizes
-[npt, dimx] = size(x);
-K = numel(ba)/(dimx+1);
-ba = reshape(ba,dimx+1,K);
-
-if alpha <= 0
-    y = Inf*ones(npt,1);
-    dydp = nan;
-    return;
-end
-
-%augment data with column of ones
-X = [ones(npt,1), x];
-
-%compute affine functions
-z = X*ba;
-
-if nargout == 1
-    y = lse_scaled(z, alpha);
-else
-    [y, dydz, dydsoftness] = lse_scaled(z, alpha);
-    if ~softness_param_is_alpha
-        dydsoftness = -dydsoftness*(alpha^2);
-    end
-    dydba = repcols(dydz, dimx+1) .* repmat(X, 1, K);
-    dydp = [dydba, dydsoftness];
-end
-
-end
-
diff --git a/legacy_m_code/softmax_affine_standalone.m b/legacy_m_code/softmax_affine_standalone.m
deleted file mode 100755
index 86c97e6..0000000
--- a/legacy_m_code/softmax_affine_standalone.m
+++ /dev/null
@@ -1,46 +0,0 @@
-function [y, dydp] = softmax_affine(x,params)
-%params may come in as
-%  1) a cell {ba, alpha}, or 
-%  2) a vector (with alpha last)
-%after reshaping (column-major), ba is dimx+1 by K
-%first row is b
-%rest is a
-
-if iscell(params)
-    ba = params{1};
-    alpha = params{2};
-else
-    ba = params(1:end-1);
-    alpha = params(end);
-end
-
-%check sizes
-[npt, dimx] = size(x);
-K = numel(ba)/(dimx+1);
-ba = reshape(ba,dimx+1,K);
-
-%augment data with column of ones
-X = [ones(npt,1), x];
-
-%compute affine functions
-z = X*ba;
-m = max(z,[],2);
-delta = z - repmat(m, 1, K);
-
-%function value
-expo = exp(alpha*delta);
-sumexpo = sum(expo,2);
-lse = log(sum(expo,2))/alpha;
-y = m + lse;
-
-if(nargout > 1)
-	%todo deal with sparse case
-	%todo deal with iscell(params) case
-	%for each K, ba is dimx+1
-	dydba = repcols(expo./repmat(sumexpo, 1, K), dimx+1) .* repmat(X, 1, K);
-    %still need to deal with alpha
-    dydalpha = (sum(expo.*delta,2)./sumexpo - lse)./alpha;
-    dydp = [dydba, dydalpha];
-end
-
-
diff --git a/setup.py b/setup.py
index c1f3429..0f2ee5c 100644
--- a/setup.py
+++ b/setup.py
@@ -3,7 +3,7 @@
 
 LICENSE = """The MIT License (MIT)
 
-Copyright (c) 2016 MIT Hoburg Research 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
@@ -26,11 +26,12 @@
 setup(
     name="gpfit",
     description="Package for fitting geometric programming models to data",
-    author="MIT Department of Aeronautics and Astronautics",
+    author="Convex Engineering Group",
     author_email="gpkit@mit.edu",
-    url="https://github.com/hoburg/gpfit",
-    install_requires=["numpy", "scipy", "gpkit"],
-    version="0.1",
+    url="https://github.com/convexengineering/gpfit",
+    python_requires=">=3.4",
+    install_requires=["numpy", "scipy", "gpkit", "matplotlib"],
+    version="0.1.0",
     packages=["gpfit", "gpfit.tests"],
     license=LICENSE,
 )


Test patch:
diff --git a/gpfit/tests/run_tests.py b/gpfit/tests/run_tests.py
index 1a423a8..42b220a 100644
--- a/gpfit/tests/run_tests.py
+++ b/gpfit/tests/run_tests.py
@@ -1,27 +1,19 @@
 "Runs all tests"
 from gpkit.tests.helpers import run_tests
-TESTS = []
-
-from gpfit.tests import t_LM
-TESTS += t_LM.TESTS
-
-from gpfit.tests import t_lse_implicit
-TESTS += t_lse_implicit.TESTS
 
-from gpfit.tests import t_lse_scaled
-TESTS += t_lse_scaled.TESTS
+TESTS = []
 
-from gpfit.tests import t_ba_init
-TESTS += t_ba_init.TESTS
+from gpfit.tests import t_least_squares
+TESTS += t_least_squares.TESTS
 
-from gpfit.tests import t_max_affine
-TESTS += t_max_affine.TESTS
+from gpfit.tests import t_logsumexp
+TESTS += t_logsumexp.TESTS
 
-from gpfit.tests import t_softmax_affine
-TESTS += t_softmax_affine.TESTS
+from gpfit.tests import t_initialize
+TESTS += t_initialize.TESTS
 
-from gpfit.tests import t_implicit_softmax_affine
-TESTS += t_implicit_softmax_affine.TESTS
+from gpfit.tests import t_classes
+TESTS += t_classes.TESTS
 
 from gpfit.tests import t_print_fit
 TESTS += t_print_fit.TESTS
@@ -29,15 +21,15 @@
 from gpfit.tests import t_plot_fit
 TESTS += t_plot_fit.TESTS
 
-from gpfit.tests import t_ex6_3
-TESTS += t_ex6_3.TESTS
-
 from gpfit.tests import t_examples
 TESTS += t_examples.TESTS
 
 from gpfit.tests import t_fit
 TESTS += t_fit.TESTS
 
+from gpfit.tests import t_constraint_set
+TESTS += t_constraint_set.TESTS
+
 
 def run(xmloutput=False):
     """Run all gpfit unit tests.
@@ -48,10 +40,10 @@ def run(xmloutput=False):
         If true, generate xml output files for continuous integration
     """
     if xmloutput:
-        run_tests(TESTS, xmloutput='test_reports')
+        run_tests(TESTS, xmloutput="test_reports")
     else:
         run_tests(TESTS)
 
 
-if __name__ == '__main__':
+if __name__ == "__main__":
     run()
diff --git a/gpfit/tests/t_classes.py b/gpfit/tests/t_classes.py
new file mode 100644
index 0000000..a0e686e
--- /dev/null
+++ b/gpfit/tests/t_classes.py
@@ -0,0 +1,83 @@
+"Test max_affine"
+import unittest
+from numpy import arange, newaxis
+from gpfit.classes import max_affine, softmax_affine, implicit_softmax_affine
+
+
+class t_max_affine(unittest.TestCase):
+    "Test max_affine"
+
+    x = arange(0.0, 16.0)[:, newaxis]
+    ba = arange(1.0, 7.0).reshape(2, 3)
+
+    y, dydba = max_affine(x, ba)
+
+    def test_y_size(self):
+        self.assertEqual(self.y.size, self.x.size)
+
+    def test_y_ndim(self):
+        self.assertEqual(self.y.ndim, 1)
+
+    def test_dydba_shape(self):
+        self.assertEqual(self.dydba.shape, (self.x.size, self.ba.size))
+
+    def test_dydba_ndim(self):
+        self.assertEqual(self.dydba.ndim, 2)
+
+
+class t_softmax_affine(unittest.TestCase):
+    "Tests softmax_affine"
+
+    x = arange(0.0, 16.0)[:, newaxis]
+    params = arange(1.0, 6.0)
+
+    y, dydp = softmax_affine(x, params)
+
+    def test_y_size(self):
+        self.assertEqual(self.y.size, self.x.size)
+
+    def test_y_ndim(self):
+        self.assertEqual(self.y.ndim, 1)
+
+    def test_dydp_shape(self):
+        self.assertEqual(self.dydp.shape, (self.x.size, self.params.size))
+
+    def test_dydp_ndim(self):
+        self.assertEqual(self.dydp.ndim, 2)
+
+
+class t_implicit_softmax_affine(unittest.TestCase):
+    "Tests implicit_softmax_affine"
+
+    x = arange(0.0, 16.0)[:, newaxis]
+    params = arange(1.0, 7.0)
+
+    y, dydp = implicit_softmax_affine(x, params)
+
+    def test_y_size(self):
+        self.assertEqual(self.y.size, self.x.size)
+
+    def test_y_ndim(self):
+        self.assertEqual(self.y.ndim, 1)
+
+    def test_dydp_shape(self):
+        self.assertEqual(self.dydp.shape, (self.x.size, self.params.size))
+
+    def test_dydp_ndim(self):
+        self.assertEqual(self.dydp.ndim, 2)
+
+
+TESTS = [
+    t_max_affine,
+    t_softmax_affine,
+    t_implicit_softmax_affine,
+]
+
+if __name__ == "__main__":
+    SUITE = unittest.TestSuite()
+    LOADER = unittest.TestLoader()
+
+    for t in TESTS:
+        SUITE.addTests(LOADER.loadTestsFromTestCase(t))
+
+    unittest.TextTestRunner(verbosity=2).run(SUITE)
diff --git a/gpfit/tests/t_constraint_set.py b/gpfit/tests/t_constraint_set.py
new file mode 100644
index 0000000..9b753cc
--- /dev/null
+++ b/gpfit/tests/t_constraint_set.py
@@ -0,0 +1,43 @@
+"unit tests for gpfit.print_fit module"
+import unittest
+import numpy as np
+from gpkit import Variable, Model
+from gpfit.fit import fit
+from gpfit.constraint_set import FitConstraintSet
+
+
+class TestFitConstraintSet(unittest.TestCase):
+    "Unit tests for FitConstraintSet"
+
+    u = np.logspace(0, np.log10(3), 501)
+    w = (u**2 + 3)/(u + 1)**2
+    x = np.log(u)
+    y = np.log(w)
+    K = 3
+    cs, _ = fit(x, y, K, "SMA")
+    uvar = Variable("u")
+    wvar = Variable("w")
+    fcs = FitConstraintSet(cs.fitdata, wvar, [uvar])
+
+    def test_fit_constraint_set(self):
+        m = Model(self.wvar, [self.fcs, self.uvar <= 2])
+        sol = m.solve(verbosity=0)
+        self.assertAlmostEqual(sol["cost"], 0.7777593)
+
+    def test_outside_bounds(self):
+        m = Model(self.wvar, [self.fcs, self.uvar <= 4])
+        sol = m.solve(verbosity=0)
+        self.assertTrue(len(sol['warnings']["Fit Out-of-Bounds"]) == 1)
+        self.assertAlmostEqual(sol["cost"], 0.7499633)
+
+
+TESTS = [TestFitConstraintSet]
+
+if __name__ == '__main__':
+    SUITE = unittest.TestSuite()
+    LOADER = unittest.TestLoader()
+
+    for t in TESTS:
+        SUITE.addTests(LOADER.loadTestsFromTestCase(t))
+
+    unittest.TextTestRunner(verbosity=2).run(SUITE)
diff --git a/gpfit/tests/t_ex6_3.py b/gpfit/tests/t_ex6_3.py
deleted file mode 100644
index 896cfce..0000000
--- a/gpfit/tests/t_ex6_3.py
+++ /dev/null
@@ -1,76 +0,0 @@
-"unit tests based on example 6.3 from Hoburg/Abbeel GPfit paper"
-import unittest
-import numpy as np
-from numpy import log, exp, vstack
-from numpy.random import random_sample
-from gpfit.fit import fit
-
-SEED = 33404
-
-
-class TestEx63ISMA(unittest.TestCase):
-    "ISMA unit tests based on example 6.3 from GPfit paper"
-
-    def test_rms_error(self):
-        np.random.seed(SEED)
-        Vdd = random_sample(1000,) + 1
-        Vth = 0.2*random_sample(1000,) + 0.2
-        P = Vdd**2 + 30*Vdd*exp(-(Vth-0.06*Vdd)/0.039)
-        u = vstack((Vdd, Vth))
-        x = log(u)
-        y = log(P)
-        K = 4
-
-        _, rms_error = fit(x, y, K, "ISMA")
-
-        self.assertTrue(rms_error < 1e-3)
-
-
-class TestEx63SMA(unittest.TestCase):
-    "SMA unit tests based on example 6.3 from GPfit paper"
-
-    def test_rms_error(self):
-        np.random.seed(SEED)
-        Vdd = random_sample(1000,) + 1
-        Vth = 0.2*random_sample(1000,) + 0.2
-        P = Vdd**2 + 30*Vdd*exp(-(Vth-0.06*Vdd)/0.039)
-        u = vstack((Vdd, Vth))
-        x = log(u)
-        y = log(P)
-        K = 4
-
-        _, rms_error = fit(x, y, K, "SMA")
-
-        self.assertTrue(rms_error < 1e-2)
-
-
-class TestEx63MA(unittest.TestCase):
-    "MA unit tests based on example 6.3 from GPfit paper"
-
-    def test_rms_error(self):
-        np.random.seed(SEED)
-        Vdd = random_sample(1000,) + 1
-        Vth = 0.2*random_sample(1000,) + 0.2
-        P = Vdd**2 + 30*Vdd*exp(-(Vth-0.06*Vdd)/0.039)
-        u = vstack((Vdd, Vth))
-        x = log(u)
-        y = log(P)
-        K = 4
-
-        _, rms_error = fit(x, y, K, "MA")
-
-        self.assertTrue(rms_error < 1e-2)
-
-
-TESTS = [TestEx63ISMA,
-         TestEx63SMA,
-         TestEx63MA]
-
-if __name__ == '__main__':
-    SUITE = unittest.TestSuite()
-    LOADER = unittest.TestLoader()
-
-    for t in TESTS:
-        SUITE.addTests(LOADER.loadTestsFromTestCase(t))
-
-    unittest.TextTestRunner(verbosity=2).run(SUITE)
diff --git a/gpfit/tests/t_examples.py b/gpfit/tests/t_examples.py
index 97612fe..cdb269a 100644
--- a/gpfit/tests/t_examples.py
+++ b/gpfit/tests/t_examples.py
@@ -34,14 +34,22 @@ def test_hoburgabbeel_ex6_1(self, example):
         self.assertTrue(example.errorSMA < 1e-4)
         self.assertTrue(example.errorISMA < 1e-5)
 
+    def test_hoburgabbeel_ex6_3(self, example):
+        """test_hoburgabbeel_ex6_3"""
+        self.assertTrue(example.errorMA < 1e-2)
+        self.assertTrue(example.errorSMA < 1e-3)
+        self.assertTrue(example.errorISMA < 1e-3)
+
 
 FILE_DIR = os.path.dirname(os.path.realpath(__file__))
-EXAMPLE_DIR = os.path.abspath(FILE_DIR + '../../../docs/source/examples')
+EXAMPLE_DIR = os.path.abspath(FILE_DIR + "../../../docs/source/examples")
 # use gpkit.tests.helpers.generate_example_tests default: only default solver
-TESTS = generate_example_tests(EXAMPLE_DIR, [TestExamples],
-                               settings["installed_solvers"])
+TESTS = generate_example_tests(
+    EXAMPLE_DIR, [TestExamples], settings["installed_solvers"]
+)
 
 
 if __name__ == "__main__":
     from gpkit.tests.helpers import run_tests
+
     run_tests(TESTS)
diff --git a/gpfit/tests/t_fit.py b/gpfit/tests/t_fit.py
index 469b5ca..1ab240a 100644
--- a/gpfit/tests/t_fit.py
+++ b/gpfit/tests/t_fit.py
@@ -1,22 +1,24 @@
-#unit tests for gpfit.fit module"
+"""unit tests for fit module"""
 import unittest
 from numpy import logspace, log10, log, vstack
 from gpfit.fit import fit
 
-class t_fit(unittest.TestCase):
+
+class TestFit(unittest.TestCase):
+    """Test fit function"""
 
     u = logspace(0, log10(3), 501)
-    w = (u**2 + 3)/(u+1)**2
+    w = (u**2 + 3)/(u + 1)**2
     x = log(u)
     y = log(w)
     K = 3
 
     def test_rms_error(self):
-        cstrt, rms_error = fit(self.x, self.y, self.K, "SMA")
+        _, rms_error = fit(self.x, self.y, self.K, "SMA")
         self.assertTrue(rms_error < 1e-4)
-        cstrt, rms_error = fit(self.x, self.y, self.K, "ISMA")
+        _, rms_error = fit(self.x, self.y, self.K, "ISMA")
         self.assertTrue(rms_error < 1e-5)
-        cstrt, rms_error = fit(self.x, self.y, self.K, "MA")
+        _, rms_error = fit(self.x, self.y, self.K, "MA")
         self.assertTrue(rms_error < 1e-2)
 
     def test_incorrect_inputs(self):
@@ -24,7 +26,7 @@ def test_incorrect_inputs(self):
             fit(self.x, vstack((self.y, self.y)), self.K, "MA")
 
 
-TESTS = [t_fit]
+TESTS = [TestFit]
 
 if __name__ == '__main__':
     SUITE = unittest.TestSuite()
diff --git a/gpfit/tests/t_implicit_softmax_affine.py b/gpfit/tests/t_implicit_softmax_affine.py
deleted file mode 100644
index 90b1fe7..0000000
--- a/gpfit/tests/t_implicit_softmax_affine.py
+++ /dev/null
@@ -1,36 +0,0 @@
-"Tests implicit_softmax_affine"
-import unittest
-from numpy import arange, newaxis
-from gpfit.implicit_softmax_affine import implicit_softmax_affine
-
-
-class t_implicit_softmax_affine(unittest.TestCase):
-    "Tests implicit_softmax_affine"
-
-    x = arange(0., 16.)[:, newaxis]
-    params = arange(1., 7.)
-
-    y, dydp = implicit_softmax_affine(x, params)
-
-    def test_y_size(self):
-        self.assertEqual(self.y.size, self.x.size)
-
-    def test_y_ndim(self):
-        self.assertEqual(self.y.ndim, 1)
-
-    def test_dydp_shape(self):
-        self.assertEqual(self.dydp.shape, (self.x.size, self.params.size))
-
-    def test_dydp_ndim(self):
-        self.assertEqual(self.dydp.ndim, 2)
-
-TESTS = [t_implicit_softmax_affine]
-
-if __name__ == '__main__':
-    SUITE = unittest.TestSuite()
-    LOADER = unittest.TestLoader()
-
-    for t in TESTS:
-        SUITE.addTests(LOADER.loadTestsFromTestCase(t))
-
-    unittest.TextTestRunner(verbosity=2).run(SUITE)
diff --git a/gpfit/tests/t_ba_init.py b/gpfit/tests/t_initialize.py
similarity index 66%
rename from gpfit/tests/t_ba_init.py
rename to gpfit/tests/t_initialize.py
index 4803384..751ae3c 100644
--- a/gpfit/tests/t_ba_init.py
+++ b/gpfit/tests/t_initialize.py
@@ -1,9 +1,9 @@
-"unit tests for ba_init function"
+"unit tests for get_initial_parameters function"
 import unittest
 import numpy as np
 from numpy import arange, newaxis, vstack, log, exp
 from numpy.random import random_sample
-from gpfit.ba_init import ba_init
+from gpfit.initialize import get_initial_parameters
 
 SEED = 33404
 
@@ -13,10 +13,11 @@ class TestMaxAffineInitK2(unittest.TestCase):
     This unit test ensures that max affine init produces an array
     of the expected shape and size
     """
-    x = arange(0., 16.)[:, newaxis]
-    y = arange(0., 16.)[:, newaxis]
+
+    x = arange(0.0, 16.0)[:, newaxis]
+    y = arange(0.0, 16.0)[:, newaxis]
     K = 2
-    ba = ba_init(x, y, K)
+    ba = get_initial_parameters(x, y, K)
 
     def test_ba_ndim_k2(self):
         self.assertEqual(self.ba.ndim, 2)
@@ -24,7 +25,7 @@ def test_ba_ndim_k2(self):
     def test_ba_shape_k2(self):
 
         _, dimx = self.x.shape
-        self.assertEqual(self.ba.shape, (dimx+1, self.K))
+        self.assertEqual(self.ba.shape, (dimx + 1, self.K))
 
 
 class TestMaxAffineInitK4(unittest.TestCase):
@@ -32,10 +33,11 @@ class TestMaxAffineInitK4(unittest.TestCase):
     This unit test ensures that max affine init produces an array
     of the expected shape and size
     """
+
     np.random.seed(SEED)
-    Vdd = random_sample(1000,) + 1
-    Vth = 0.2*random_sample(1000,) + 0.2
-    P = Vdd**2 + 30*Vdd*exp(-(Vth-0.06*Vdd)/0.039)
+    Vdd = random_sample(1000) + 1
+    Vth = 0.2*random_sample(1000) + 0.2
+    P = Vdd**2 + 30*Vdd*exp(-(Vth - 0.06*Vdd)/0.039)
     u = vstack((Vdd, Vth))
     x = log(u)
     y = log(P)
@@ -43,15 +45,15 @@ class TestMaxAffineInitK4(unittest.TestCase):
     y = y.reshape(y.size, 1)
     K = 4
 
-    ba = ba_init(x, y, K)
+    ba = get_initial_parameters(x, y, K)
 
     def test_ba_shape_k4(self):
         self.assertEqual(self.ba.shape, (3, 4))
 
-TESTS = [TestMaxAffineInitK2,
-         TestMaxAffineInitK4]
 
-if __name__ == '__main__':
+TESTS = [TestMaxAffineInitK2, TestMaxAffineInitK4]
+
+if __name__ == "__main__":
     SUITE = unittest.TestSuite()
     LOADER = unittest.TestLoader()
 
diff --git a/gpfit/tests/t_LM.py b/gpfit/tests/t_least_squares.py
similarity index 76%
rename from gpfit/tests/t_LM.py
rename to gpfit/tests/t_least_squares.py
index e0afe91..26c9848 100644
--- a/gpfit/tests/t_LM.py
+++ b/gpfit/tests/t_least_squares.py
@@ -1,20 +1,20 @@
 "Tests levenberg_marquardt"
 import unittest
 from numpy import arange, newaxis
-from gpfit.levenberg_marquardt import levenberg_marquardt
-from gpfit.max_affine import max_affine
+from gpfit.least_squares import levenberg_marquardt
+from gpfit.classes import max_affine
 
 
 def rfun(params):
     "A specific residual function."
-    [yhat, drdp] = max_affine(arange(0., 16.)[:, newaxis], params)
-    r = yhat - arange(0., 16.)[:, newaxis].T[0]
+    [yhat, drdp] = max_affine(arange(0.0, 16.0)[:, newaxis], params)
+    r = yhat - arange(0.0, 16.0)[:, newaxis].T[0]
     return r, drdp
 
 
 class t_levenberg_marquardt(unittest.TestCase):
     "Tests levenberg_marquardt"
-    initparams = arange(1., 5.)
+    initparams = arange(1.0, 5.0)
     params, RMStraj = levenberg_marquardt(rfun, initparams)
 
     def test_params_size(self):
@@ -31,9 +31,10 @@ def test_rmstraj_ndim(self):
 
         self.assertEqual(self.RMStraj.ndim, 1)
 
+
 TESTS = [t_levenberg_marquardt]
 
-if __name__ == '__main__':
+if __name__ == "__main__":
     SUITE = unittest.TestSuite()
     LOADER = unittest.TestLoader()
 
diff --git a/gpfit/tests/t_lse_implicit.py b/gpfit/tests/t_logsumexp.py
similarity index 57%
rename from gpfit/tests/t_lse_implicit.py
rename to gpfit/tests/t_logsumexp.py
index 2f778cd..d1d1d51 100644
--- a/gpfit/tests/t_lse_implicit.py
+++ b/gpfit/tests/t_logsumexp.py
@@ -1,8 +1,8 @@
-"unit tests for implicit log-sum-exp implementation"
+"unit tests for log-sum-exp functions"
 import unittest
 import numpy as np
 from numpy import array, arange
-from gpfit.lse_implicit import lse_implicit
+from gpfit.logsumexp import lse_implicit, lse_scaled
 
 SEED = 33404
 
@@ -10,8 +10,8 @@
 class TestLSEimplicit1D(unittest.TestCase):
     "tests with one-dimensional input"
 
-    x = arange(1., 31.).reshape(15, 2)
-    alpha = arange(1., 3.)
+    x = arange(1.0, 31.0).reshape(15, 2)
+    alpha = arange(1.0, 3.0)
     y, dydx, dydalpha = lse_implicit(x, alpha)
 
     def test_y_ndim(self):
@@ -32,6 +32,7 @@ def test_dydalpha_ndim(self):
     def test_dydalpha_size(self):
         self.assertEqual(self.dydalpha.shape[0], self.x.shape[0])
 
+
 class TestLSEimplicit2D(unittest.TestCase):
     "tests with 2D input"
 
@@ -49,10 +50,42 @@ def test_dydx_shape(self):
     def test_dydalpha_shape(self):
         self.assertEqual(self.dydalpha.shape, self.x.shape)
 
-TESTS = [TestLSEimplicit1D,
-         TestLSEimplicit2D]
 
-if __name__ == '__main__':
+class TestLSEScaled(unittest.TestCase):
+    "Test lse_implicit"
+
+    x = arange(1.0, 31.0).reshape(15, 2)
+    alpha = 1.0
+    y, dydx, dydalpha = lse_scaled(x, alpha)
+
+    def test_y_ndim(self):
+        self.assertEqual(self.y.ndim, 1)
+
+    def test_y_size(self):
+        self.assertEqual(self.y.size, self.x.shape[0])
+
+    def test_dydx_ndim(self):
+        self.assertEqual(self.dydx.ndim, 2)
+
+    def test_dydx_shape_0(self):
+        self.assertEqual(self.dydx.shape[0], self.x.shape[0])
+
+    def test_dydalpha_ndim(self):
+        self.assertEqual(self.dydalpha.ndim, 1)
+
+    def test_dydalpha_size(self):
+        self.assertEqual(self.dydalpha.size, self.x.shape[0])
+
+    # test alpha is integer? negative? 0? array?
+
+
+TESTS = [
+    TestLSEimplicit1D,
+    TestLSEimplicit2D,
+    TestLSEScaled,
+]
+
+if __name__ == "__main__":
     SUITE = unittest.TestSuite()
     LOADER = unittest.TestLoader()
 
diff --git a/gpfit/tests/t_lse_scaled.py b/gpfit/tests/t_lse_scaled.py
deleted file mode 100644
index abd3333..0000000
--- a/gpfit/tests/t_lse_scaled.py
+++ /dev/null
@@ -1,48 +0,0 @@
-"Test lse_scaled"
-import unittest
-from numpy import arange
-from gpfit.lse_scaled import lse_scaled
-
-
-class t_lse_scaled(unittest.TestCase):
-    "Test lse_implicit"
-
-    x = arange(1., 31.).reshape(15, 2)
-    alpha = 1.
-    y, dydx, dydalpha = lse_scaled(x, alpha)
-
-    def test_y_ndim(self):
-        self.assertEqual(self.y.ndim, 1)
-
-    def test_y_size(self):
-        self.assertEqual(self.y.size, self.x.shape[0])
-
-    def test_dydx_ndim(self):
-        self.assertEqual(self.dydx.ndim, 2)
-
-    def test_dydx_shape_0(self):
-        self.assertEqual(self.dydx.shape[0], self.x.shape[0])
-
-    def test_dydx_shape_1(self):
-        pass
-        # self.assertEqual(self.dydx.shape[0], ???????)
-
-    def test_dydalpha_ndim(self):
-        self.assertEqual(self.dydalpha.ndim, 1)
-
-    def test_dydalpha_size(self):
-        self.assertEqual(self.dydalpha.size, self.x.shape[0])
-
-    # test alpha is integer? negative? 0? array?
-
-
-TESTS = [t_lse_scaled]
-
-if __name__ == '__main__':
-    SUITE = unittest.TestSuite()
-    LOADER = unittest.TestLoader()
-
-    for t in TESTS:
-        SUITE.addTests(LOADER.loadTestsFromTestCase(t))
-
-    unittest.TextTestRunner(verbosity=2).run(SUITE)
diff --git a/gpfit/tests/t_max_affine.py b/gpfit/tests/t_max_affine.py
deleted file mode 100644
index 2ba0457..0000000
--- a/gpfit/tests/t_max_affine.py
+++ /dev/null
@@ -1,36 +0,0 @@
-"Test max_affine"
-import unittest
-from numpy import arange, newaxis
-from gpfit.max_affine import max_affine
-
-
-class t_max_affine(unittest.TestCase):
-    "Test max_affine"
-
-    x = arange(0., 16.)[:, newaxis]
-    ba = arange(1., 7.).reshape(2, 3)
-
-    y, dydba = max_affine(x, ba)
-
-    def test_y_size(self):
-        self.assertEqual(self.y.size, self.x.size)
-
-    def test_y_ndim(self):
-        self.assertEqual(self.y.ndim, 1)
-
-    def test_dydba_shape(self):
-        self.assertEqual(self.dydba.shape, (self.x.size, self.ba.size))
-
-    def test_dydba_ndim(self):
-        self.assertEqual(self.dydba.ndim, 2)
-
-TESTS = [t_max_affine]
-
-if __name__ == '__main__':
-    SUITE = unittest.TestSuite()
-    LOADER = unittest.TestLoader()
-
-    for t in TESTS:
-        SUITE.addTests(LOADER.loadTestsFromTestCase(t))
-
-    unittest.TextTestRunner(verbosity=2).run(SUITE)
diff --git a/gpfit/tests/t_plot_fit.py b/gpfit/tests/t_plot_fit.py
index 12e6e38..097a36d 100644
--- a/gpfit/tests/t_plot_fit.py
+++ b/gpfit/tests/t_plot_fit.py
@@ -6,12 +6,12 @@
 
 class TestPlotFit(unittest.TestCase):
     "Unit tests for plot_fit_1d"
+    N = 51
+    u = np.logspace(0, np.log10(3), N)
+    w = (u**2 + 3)/(u + 1)**2
 
     def test_plot_fit_1d(self):
-        N = 51
-        U = np.logspace(0, np.log10(3), N)
-        W = (U**2+3) / (U+1)**2
-        plot_fit_1d(U, W, K=2, fitclass='SMA', plotspace="linear")
+        plot_fit_1d(self.u, self.w, K=2, fitclass='SMA', plotspace="linear")
 
 
 TESTS = [TestPlotFit]
diff --git a/gpfit/tests/t_print_fit.py b/gpfit/tests/t_print_fit.py
index 0f9d1fb..e8fc36c 100644
--- a/gpfit/tests/t_print_fit.py
+++ b/gpfit/tests/t_print_fit.py
@@ -1,49 +1,61 @@
 "unit tests for gpfit.print_fit module"
 import unittest
 from numpy import array
-from gpfit.print_fit import print_MA, print_SMA, print_ISMA
+from gpfit.print_fit import print_ma, print_sma, print_isma
 
 A = array([2, 3, 4, 6, 7, 8, 10, 11, 12])
 B = array([1, 5, 9])
-ALPHA = array([1./13., 1./14., 1./15.])
+ALPHA = array([1/13, 1/14, 1/15])
 DIM = 3
 K = 3
 
 
 class TestPrintFit(unittest.TestCase):
-    "Unit tests for print_MA, print_SMA, print_ISMA"
+    "Unit tests for print_ma, print_sma, print_isma"
 
     def test_ma(self):
-        strings = print_MA(A, B, DIM, K)
+        strings = print_ma(A, B, DIM, K)
         self.assertEqual(
             strings,
-            ['w = 2.71828 * (u_1)**2 * (u_2)**3 * (u_3)**4',
-             'w = 148.413 * (u_1)**6 * (u_2)**7 * (u_3)**8',
-             'w = 8103.08 * (u_1)**10 * (u_2)**11 * (u_3)**12'])
+            [
+                "w = 2.71828 * (u_1)**2 * (u_2)**3 * (u_3)**4",
+                "w = 148.413 * (u_1)**6 * (u_2)**7 * (u_3)**8",
+                "w = 8103.08 * (u_1)**10 * (u_2)**11 * (u_3)**12",
+            ],
+        )
 
     def test_sma(self):
-        strings = print_SMA(A, B, ALPHA[0], DIM, K)
-        self.assertEqual(strings, ['w**0.0769231 = 1.07996 * (u_1)**0.153846 '
-                                   '* (u_2)**0.230769 * (u_3)**0.307692',
-                                   '    + 1.46905 * (u_1)**0.461538 * '
-                                   '(u_2)**0.538462 * (u_3)**0.615385',
-                                   '    + 1.99832 * (u_1)**0.769231 * '
-                                   '(u_2)**0.846154 * (u_3)**0.923077'])
+        strings = print_sma(A, B, ALPHA[0], DIM, K)
+        self.assertEqual(
+            strings,
+            [
+                "w**0.0769231 = 1.07996 * (u_1)**0.153846 "
+                "* (u_2)**0.230769 * (u_3)**0.307692",
+                "    + 1.46905 * (u_1)**0.461538 * "
+                "(u_2)**0.538462 * (u_3)**0.615385",
+                "    + 1.99832 * (u_1)**0.769231 * "
+                "(u_2)**0.846154 * (u_3)**0.923077",
+            ],
+        )
 
     def test_isma(self):
-        strings = print_ISMA(A, B, ALPHA, DIM, K)
+        strings = print_isma(A, B, ALPHA, DIM, K)
         self.assertEqual(
             strings,
-            ['1 = (1.07996/w**0.0769231) * (u_1)**0.153846 * '
-             '(u_2)**0.230769 * (u_3)**0.307692',
-             '    + (1.42924/w**0.0714286) * (u_1)**0.428571 * '
-             '(u_2)**0.5 * (u_3)**0.571429',
-             '    + (1.82212/w**0.0666667) * (u_1)**0.666667 * '
-             '(u_2)**0.733333 * (u_3)**0.8'])
+            [
+                "1 = (1.07996/w**0.0769231) * (u_1)**0.153846 * "
+                "(u_2)**0.230769 * (u_3)**0.307692",
+                "    + (1.42924/w**0.0714286) * (u_1)**0.428571 * "
+                "(u_2)**0.5 * (u_3)**0.571429",
+                "    + (1.82212/w**0.0666667) * (u_1)**0.666667 * "
+                "(u_2)**0.733333 * (u_3)**0.8",
+            ],
+        )
+
 
 TESTS = [TestPrintFit]
 
-if __name__ == '__main__':
+if __name__ == "__main__":
     SUITE = unittest.TestSuite()
     LOADER = unittest.TestLoader()
 
diff --git a/gpfit/tests/t_softmax_affine.py b/gpfit/tests/t_softmax_affine.py
deleted file mode 100644
index adfd117..0000000
--- a/gpfit/tests/t_softmax_affine.py
+++ /dev/null
@@ -1,36 +0,0 @@
-"Tests softmax_affine"
-import unittest
-from numpy import arange, newaxis
-from gpfit.softmax_affine import softmax_affine
-
-
-class t_softmax_affine(unittest.TestCase):
-    "Tests softmax_affine"
-
-    x = arange(0., 16.)[:, newaxis]
-    params = arange(1., 6.)
-
-    y, dydp = softmax_affine(x, params)
-
-    def test_y_size(self):
-        self.assertEqual(self.y.size, self.x.size)
-
-    def test_y_ndim(self):
-        self.assertEqual(self.y.ndim, 1)
-
-    def test_dydp_shape(self):
-        self.assertEqual(self.dydp.shape, (self.x.size, self.params.size))
-
-    def test_dydp_ndim(self):
-        self.assertEqual(self.dydp.ndim, 2)
-
-TESTS = [t_softmax_affine]
-
-if __name__ == '__main__':
-    SUITE = unittest.TestSuite()
-    LOADER = unittest.TestLoader()
-
-    for t in TESTS:
-        SUITE.addTests(LOADER.loadTestsFromTestCase(t))
-
-    unittest.TextTestRunner(verbosity=2).run(SUITE)
diff --git a/legacy_m_code/LM_test.m b/legacy_m_code/LM_test.m
deleted file mode 100755
index 70426b8..0000000
--- a/legacy_m_code/LM_test.m
+++ /dev/null
@@ -1,31 +0,0 @@
-function [params, RMStraj, X,RESNORM,RESIDUAL,EXITFLAG,OUTPUT] = LM_test()
-
-npt = 1000;
-dimx = 2;
-x = randn(npt,dimx);
-y = sum(x.^2, 2);
-%y = max(sum(x.^2, 2), 5*max(x,[],2));
-
-K = 4;
-ba = max_affine_init(x, y, K);
-alphainit = 10;
-
-rfun = @(p) generic_resid_fun(@softmax_affine, x, y, p);
-
-%tic;
-[params, RMStraj] = LM(rfun, [ba(:); alphainit]);
-myLMtime = toc;
-tic;
-[X,RESNORM,RESIDUAL,EXITFLAG,OUTPUT] = lsqnonlin(rfun,[ba(:); alphainit],[],[],optimset('Jacobian', 'on'));
-trusttime = toc;
-
-disp('        finalRMS    iterations     time');
-disp(['my LM   ', num2str(norm(rfun(params))/sqrt(npt)), '       ', num2str(length(RMStraj)), '        ', num2str(myLMtime)]);
-disp(['trust   ', num2str(sqrt(RESNORM/npt)), '       ', num2str(OUTPUT.iterations), '        ', num2str(trusttime)]);
-
-figure(1); clf(1);
-plot(x(:,1), x(:,1).^2, '.'); hold on;
-xplo = [linspace(min(x(:,1)), max(x(:,1)), 1000)', zeros(1000, dimx-1)];
-plot(xplo(:,1), softmax_affine(xplo, params), 'r');
-
-end
\ No newline at end of file
diff --git a/legacy_m_code/circuit_test.m b/legacy_m_code/circuit_test.m
deleted file mode 100644
index 3c7afb6..0000000
--- a/legacy_m_code/circuit_test.m
+++ /dev/null
@@ -1,42 +0,0 @@
-close all
-vdd = linspace(1.0,2.0,30);
-vth = linspace(0.2,0.4,30);
-[VDD, VTH] = meshgrid(vdd, vth);
-
-
-P = VDD.^2 + 30.*VDD.*exp(-(VTH-0.06.*VDD)/0.039);
-surf(VDD,VTH,P)
-hold on;
-% surf(log(VDD),log(VTH),log(P))
-w = P(:);
-% [w, i] = sort(w);
-
-
-VDD = VDD(:);
-% VDD = VDD(i);
-VTH = VTH(:);
-% VTH = VTH(i);
-u = [VDD,VTH];
-
-x = log(u);
-y = log(w);
-
-% plot(x(:,2),y)
-s = compare_fits(x,y, 3,1);
-PAR = s.softmax_optMAinit.params{1};
-% B = PAR(1:3:end-1);
-% A = PAR([2,3,5,6,8,9]);
-B = PAR(1:3:end-2);
-A = PAR([2,3,5,6,8,9]);
-
-alpha = 1/PAR(end)
-coeff = exp(alpha*B)
-powers = alpha.*A
-
-[VDD, VTH] = meshgrid(vdd, vth);
-P_approx = (...
-            coeff(1).*VDD.^(powers(1)).*VTH.^(powers(2))+ ...
-            coeff(2).*VDD.^(powers(3)).*VTH.^(powers(4)) + ...
-            coeff(3).*VDD.^(powers(5)).*VTH.^(powers(6))...
-            ).^(1/alpha);
-surf(VDD,VTH,P_approx)
\ No newline at end of file
diff --git a/legacy_m_code/testScript.m b/legacy_m_code/testScript.m
deleted file mode 100644
index c978945..0000000
--- a/legacy_m_code/testScript.m
+++ /dev/null
@@ -1,23 +0,0 @@
-% x = linspace(0,10,11);
-% x = x';
-x1 = 0:15;
-x2 = x1.^2;
-x3 = x1.^3;
-% x = x1';
-x = [x1', x2'];
-% x = [x1', x2', x3'];
-
-% y1 = x(1:round(end/2));
-% y2 = 2*x(round(end/2) + 1:end);
-% y = [y1; y2];
-y = [0,1,2,3,4,5,6,7,8,9,10,12,14,16,18,20]';
-
-Ks = 2;
-ntry = 1;
-
-s = compare_fits(x,y,Ks,ntry);
-
-disp(s.maxaffine.params{1})
-disp(s.softmax_optMAinit.params{1})
-disp(s.softmax_originit.params{1})
-disp(s.implicit_originit.params{1})
\ No newline at end of file