Repository: hoburg/gpkit
Pull request #571 (2016-03-10 19:36:41)

add boundedness debugging tool
  - [x] Turn into `BoundedConstraintSet` and move functions into `__init__` and `sens_from_dual/sens_from_gpconstr`, with a verbosity flag to print during solve.

bqpd: @whoburg, ready for review.
whoburg: initial review complete. I tend to think this could be even simpler than it is now -- just a dumb but reliable check on the value of the variables in the solution and whether they slid out to the bounds -- no sensitivities needed.

That would reduce the number of tolerance knobs and I think make this more intuitive.

?
whoburg: I see the checkbox above - are we aiming to merge this and do that later, or get that done before merging? I'm fine with either, once the other review items are addressed/decided on.
bqpd: I believe the current test is an example of a place where two of the
variables don't hit the bounds but are sensitive to them.
On Mar 25, 2016 12:56, "Warren Hoburg" <notifications@github.com> wrote:

> I see the checkbox above - are we aiming to merge this and do that later,
> or get that done before merging? I'm fine with either, once the other
> review items are addressed/decided on.
>
> —
> You are receiving this because you authored the thread.
> Reply to this email directly or view it on GitHub
> <https://github.com/hoburg/gpkit/pull/571#issuecomment-201362373>
>

bqpd: @whoburg, I made the output dictionary clear as to what the entries mean. Ready for review.
whoburg: @bqpd, review my latest commit to this branch, and then I think this is ready to merge
bqpd: Do we want to make this a constraint set before merging? Also there's a missing file docstring for `unbounded.py`
pgkirsch: Is this ready to merge @bqpd?
bqpd: I was thinking of turning it into a ConstraintSet first, sorry.
whoburg: also just noticed that `determine_unbounded_variables` seems to call `localsolve` -- need to make it work with gps that need `solve`?
bqpd: Yeah, that's why it should be a ConstraintSet. :p
mayork: I get this error when running this line of code bounds = determine_unbounded_variables(design, "cvxopt") where design is a model

```Traceback (most recent call last):
  File "C:\Users\mayork\Documents\GitHub\gpkit-models\gpkitmodels\aircraft\engine_component_testing.py", line 182, in <module>
    bounds = determine_unbounded_variables(design, "cvxopt")
  File "C:\Python27\lib\site-packages\gpkit\tools.py", line 139, in determine_unbounded_variables
    sol = m.solve(solver, verbosity, **kwargs)
  File "C:\Python27\lib\site-packages\gpkit\constraints\prog_factories.py", line 79, in solvefn
    self.program, solvefn = genfunction(self, verbosity-1)
  File "C:\Python27\lib\site-packages\gpkit\constraints\prog_factories.py", line 34, in programify
    prog = program(self.cost, self, substitutions, verbosity)
  File "C:\Python27\lib\site-packages\gpkit\geometric_program.py", line 72, in __init__
    super(GeometricProgram, self).init_from_nomials(self.posynomials)
  File "C:\Python27\lib\site-packages\gpkit\nomials\data.py", line 78, in init_from_nomials
    exps = functools_reduce(add, (tuple(s.exps) for s in nomials))
  File "C:\Python27\lib\site-packages\gpkit\nomials\data.py", line 78, in <genexpr>
    exps = functools_reduce(add, (tuple(s.exps) for s in nomials))
AttributeError: 'numpy.float64' object has no attribute 'exps'```
whoburg: shall we try to rebase this onto master?
mayork: if this could be a priority that would be awesome, I think both would really help debug the engine code
bqpd: Sure thing! I'll rewrite it as a constraint set
On Jun 20, 2016 13:07, "mayork" <notifications@github.com> wrote:

> if this could be a priority that would be awesome, I think both would
> really help debug the engine code
>
> —
> You are receiving this because you were mentioned.
> Reply to this email directly, view it on GitHub
> <https://github.com/hoburg/gpkit/pull/571#issuecomment-227269708>, or mute
> the thread
> <https://github.com/notifications/unsubscribe/ABagGNIGA_mo28u5hDQqmQ0YMBSZ87TVks5qNwD3gaJpZM4HuEOB>
> .
>

mayork: I know ned's writing the gpkit paper right now so this could take a little bit, but is there a fast local hack I could do to get this running? If not its not a big deal
whoburg: I'll take a quick look at this. Getting on an airplane so response will be ~3 hours from now
whoburg: I'm going to force push a rebased version of this (with a bug fix) to this PR. I still have a local copy of the previous version in case @bqpd needs it.
whoburg: conflicts when rebasing were limited to import errors.
whoburg: @mayork, give this a try. I wasn't sure exactly what to run to reproduce the error you saw, so not sure yet whether or not this will fix it.
whoburg: Another todo before this is closed -- the example `docs/source/examples/unbounded.py` has an inline comment

`# sol = m.solve("mosek")  # returns UNKNOWN`

... but this isn't what happens on my machine. I get a solver cost did not match computed cost error from gpkit's solution checking. We should 

- [ ] figure out why this changed, and
- [ ] unit test that the inline comment is accurate
whoburg: The example should probably be changed to not depend on mosek -- this fails the windows machine that only has cvxopt installed.
mayork: @whoburg I think this fixed the error. I actually had that weird importing problem again but worked around it. Might be worth discussing that problem at some point

Also, I was thinking that before this is merged in to master there should be a way to handle signomial programs. In determine_unbounded_variables solve is called, maybe an if else to check if signomials are enabled then call localsolve?
whoburg: I'm going to rebase this again once #751 is merged -- that should allow me to clear up the windows unit test failure
whoburg: I'm going to rebase force push this again to get #751 
whoburg: @galbramc curious to get your thoughts on how to deal with this. Look at my last commit to this PR and you'll see I'm just trying to run this test only if mosek is installed. We're now failing because the output changed. I think it has something to do with the fact that the CI script does `echo 'installed_solvers : cvxopt'`. Any thoughts on best way to make this work?
whoburg: I'm really confused about why all tests but windows succeeded here: https://acdl.mit.edu/csi/job/gpkit_PullRequest_Unit/1331/ -- but are now failing.
galbramc: Ahh, the ```echo 'installed_solvers : cvxopt'``` was a hack that @bqpd suggested we use so we would only test cvxopt even when both solvers are installed. I'm all for changing the hack if you have a better solution.
galbramc: In regards to 1331, I don't know why you get ```ImportError: Could not load MOSEK library: ImportError('No module named lib',)```. I don't think mosek is even installed on that VM...
whoburg: right, we got it in 1331 because `docs/source/examples/unbounded.py` explicitly called `solve("mosek")`. I tried to fix that in the last commit, but now we have this output changing issue.

I guess maybe we need to refactor the example so that it's solver-independent. that's probably the right solution.
galbramc: You know better than me how to fix that particular problem...
whoburg: yep -- sorry for the distraction
galbramc: No problem. Do you want me to change the script in regards to ```echo 'installed_solvers : cvxopt'```?
whoburg: Not yet. That should happen as part of the refactor proposed in #258. Maybe @aechevarria1 is interested.
galbramc: There is also an issue (I can't remember the number) with storing the settings file in the pip cache while you are at it.
whoburg: #679 
mayork: I was talking to Woody one day and mentioned I once saw a model that solved when it wasn't bounded but didn't solve bounded....I just ran into something similar while doing some debugging. When bounded the model at commit https://github.com/hoburg/gpkit-models/pull/73/commits/891bded301f7a05625937297c88c66941aedeaea solves, unbounded it remains unsolved after 200 iterations. I'm not sure what's going on but I thought maybe one of y'all would want to look at it to understand better why this occurs
galbramc: Sounds like you should create an issue with enough information ( including an example that reproduces the problem ) so they can sort out the possible bug.
whoburg: Noticed using the boundedness debugging tool that it's not possible to access the solution of the bounded problem (since the method returns the boundedness dict and doesn't save the solution array).

Also, not sure if this is intended, but (on the new unit test test_unbounded_debugging, for example), after running `bounds = determine_unbounded_variables(m)` on a model `m`, `m.solution` is None.

Also, for an SP model, after a determine_unbounded_variables solve, `m.program` ends up being None as well, and there is no attribute `m.gps`.

Not sure if some of these are intentional design decisions, but currently it's difficult to query the values of variables after a boundedness solve.
bqpd: @whoburg, ready for review
mayork: I haven't switched over to this yet so I've been modifying a local copy of this...one thing I just did that seems potentially useful is to only lower bound or only upper bound variables. I think this is helpful for debugging...if y'all think this could be useful I can take a stab at adding a keyword to allow only upper or lower bounding once this is merged into master
bqpd: sounds great! Make an issue for it now so we remember?

On Thu, Aug 4, 2016 at 2:23 PM, mayork <notifications@github.com> wrote:

> I haven't switched over to this yet so I've been modifying a local copy of
> this...one thing I just did that seems potentially useful is to only lower
> bound or only upper bound variables. I think this is helpful for
> debugging...if y'all think this could be useful I can take a stab at adding
> a keyword to allow only upper or lower bounding once this is merged into
> master
>
> —
> You are receiving this because you were mentioned.
> Reply to this email directly, view it on GitHub
> <https://github.com/hoburg/gpkit/pull/571#issuecomment-237687819>, or mute
> the thread
> <https://github.com/notifications/unsubscribe-auth/ABagGA7MQwrgFINhci7r16ugddgwrg_iks5qclhNgaJpZM4HuEOB>
> .
>

whoburg: code review complete.

Patch:
diff --git a/.pylintrc b/.pylintrc
index 436477b..3e2f2e2 100644
--- a/.pylintrc
+++ b/.pylintrc
@@ -96,7 +96,7 @@ msg-template={path}:{line}: [{msg_id}({symbol}), {obj}] {msg}
 bad-functions=filter,input
 
 # Good variable names which should always be accepted, separated by a comma
-good-names=i,j,k,_,gp,sp,to,A,c,cs
+good-names=i,j,k,_,gp,sp,to,A,c,cs,x0
 
 # Bad variable names which should always be refused, separated by a comma
 bad-names=foo,bar,baz,toto,tutu,tata
@@ -133,7 +133,7 @@ attr-rgx=[a-z_][a-z0-9_]{2,30}$
 attr-name-hint=[a-z_][a-z0-9_]{2,30}$
 
 # Regular expression matching correct argument names
-argument-rgx=[a-z_][a-z0-9_]{2,30}|x0|nu|la|ax$
+argument-rgx=[a-z_][a-z0-9_]{2,30}|nu|la|ax$
 
 # Naming hint for argument names
 argument-name-hint=[a-z_][a-z0-9_]{2,30}$
@@ -320,7 +320,7 @@ exclude-protected=_asdict,_fields,_replace,_source,_make
 [DESIGN]
 
 # Maximum number of arguments for function / method
-max-args=6
+max-args=7
 
 # Argument names that match this expression will be ignored. Default to name
 # with leading underscore
diff --git a/docs/apidoc.sh b/docs/apidoc.sh
index f16b03b..1c10c08 100755
--- a/docs/apidoc.sh
+++ b/docs/apidoc.sh
@@ -1,3 +1,9 @@
 rm source/autodoc/*
 sphinx-apidoc ../gpkit -o source/autodoc
-sed -i '1i Glossary\n********\n\n*For an alphabetical listing of all commands, check out the* :ref:`genindex`\n' source/autodoc/gpkit.rst
+
+# Delete fist 3 lines
+sed -i'' -e '1,3d' source/autodoc/gpkit.rst
+
+# Add header
+header="Glossary\n********\n\n*For an alphabetical listing of all commands, check out the* :ref:\`genindex\`\n"
+echo $header | cat - source/autodoc/gpkit.rst > gpkit.rst && mv gpkit.rst source/autodoc/gpkit.rst
diff --git a/docs/source/advancedcommands.rst b/docs/source/advancedcommands.rst
index ce7f460..227e653 100644
--- a/docs/source/advancedcommands.rst
+++ b/docs/source/advancedcommands.rst
@@ -120,7 +120,7 @@ Note that units are preserved, and that the value can be either a string (in whi
 Substituting with replacement
 ------------------------------
 
-Any of the substitutions above can be run with ``p.sub(*args, replace=True)`` to clobber any previously-substitued values.
+Any of the substitutions above can be run with ``p.sub(*args, replace=True)`` to clobber any previously-substituted values.
 
 Fixed Variables
 ---------------
@@ -159,7 +159,7 @@ Sweeps
 Declaring Sweeps
 ----------------
 
-Sweeps are useful for analyzing tradeoff surfaces. A sweep “value” is an Iterable of numbers, e.g. ``[1, 2, 3]``. Variables are swept when their substitution value takes the form ``('sweep', Iterable), (e.g. 'sweep', np.linspace(1e6, 1e7, 100))``. During variable declaration, giving an Iterable value for a Variable is assumed to be giving it a sweeep value: for example, ``x = Variable("x", [1, 2, 3]``. Sweeps can also be declared during later substitution (``gp.sub("x", ('sweep', [1, 2, 3]))``, or if the variable was already substituted for a constant, ``gp.sub("x", ('sweep', [1, 2, 3]), replace=True))``.
+Sweeps are useful for analyzing tradeoff surfaces. A sweep “value” is an Iterable of numbers, e.g. ``[1, 2, 3]``. Variables are swept when their substitution value takes the form ``('sweep', Iterable), (e.g. 'sweep', np.linspace(1e6, 1e7, 100))``. During variable declaration, giving an Iterable value for a Variable is assumed to be giving it a sweep value: for example, ``x = Variable("x", [1, 2, 3]``. Sweeps can also be declared during later substitution (``gp.sub("x", ('sweep', [1, 2, 3]))``, or if the variable was already substituted for a constant, ``gp.sub("x", ('sweep', [1, 2, 3]), replace=True))``.
 
 Solving Sweeps
 --------------
@@ -233,3 +233,26 @@ Example Usage
 The ``normsub`` argument specifies an expected value for your solution to normalize the different :math:`g_i` (you can also do this by hand). The feasibility of the problem should not depend on the normalization, but the spacing of the sweep will.
 
 The ``sweep`` argument specifies what points between 0 and 1 you wish to sample the weights at. If you want different resolutions or spacings for different weights, the ``sweeps`` argument accepts a list of sweep arrays.
+
+
+Debugging
+=========
+
+Unbounded variables
+-------------------
+In some cases a model will not solve because its variables are pushing to 0 or infinity. If the solver catches such behaviour it will return ``dual infeasible`` (or equivalent), but sometimes solvers do not catch it and return ``unknown``.
+
+``gpkit.constraints.bounded.BoundedConstraintSet`` is a
+simple tool that attempts to detect unbounded variables and get unbounded models to solve by adding extremely large upper bounds and extremely small lower bounds to all variables in a ConstraintSet.
+
+When a model with an BoundedConstraintSet is solved, it checks whether any variables slid off to the bounds, notes this in the solution dictionary and prints a warning (if verbosity is greater than 0).
+
+For example, Mosek returns ``DUAL_INFEAS_CER`` when attempting to solve the following model:
+
+.. literalinclude:: examples/unbounded.py
+
+Upon viewing the printed output,
+
+.. literalinclude:: examples/unbounded_output.txt
+
+it becomes clear that the problem is, unsurprisingly, that ``x`` is unbounded below in the original model.
diff --git a/docs/source/autodoc/gpkit.rst b/docs/source/autodoc/gpkit.rst
index fbf77f4..3ee63c6 100644
--- a/docs/source/autodoc/gpkit.rst
+++ b/docs/source/autodoc/gpkit.rst
@@ -3,9 +3,6 @@ Glossary
 
 *For an alphabetical listing of all commands, check out the* :ref:`genindex`
 
-gpkit package
-=============
-
 Subpackages
 -----------
 
diff --git a/docs/source/conf.py b/docs/source/conf.py
index 22b932c..9bb761c 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -18,7 +18,7 @@ import os
 # If extensions (or modules to document with autodoc) are in another directory,
 # add these directories to sys.path here. If the directory is relative to the
 # documentation root, use os.path.abspath to make it absolute, like shown here.
-sys.path.insert(0, os.path.abspath('.'))
+sys.path.insert(0, os.path.abspath('../../../gpkit'))
 
 # -- General configuration ------------------------------------------------
 
diff --git a/docs/source/examples/external_constraint.py b/docs/source/examples/external_constraint.py
index a8a98a9..dd55262 100644
--- a/docs/source/examples/external_constraint.py
+++ b/docs/source/examples/external_constraint.py
@@ -1,5 +1,6 @@
 "Can be found in gpkit/docs/source/examples/external_constraint.py"
 from gpkit import ConstraintSet
+from gpkit.exceptions import InvalidGPConstraint
 from external_function import external_code
 
 
@@ -21,7 +22,7 @@ class ExternalConstraint(ConstraintSet):
     # Prevents the ExternalConstraint class from solving in a GP, thus forcing
     # iteration
     def as_posyslt1(self):
-        raise TypeError("ExternalConstraint Model cannot solve as a GP.")
+        raise InvalidGPConstraint("ExternalConstraint cannot solve as a GP.")
 
     # Returns the ExternalConstraint class as a GP compatible constraint when
     # requested by the GPkit solver
diff --git a/docs/source/examples/unbounded.py b/docs/source/examples/unbounded.py
new file mode 100644
index 0000000..04d7159
--- /dev/null
+++ b/docs/source/examples/unbounded.py
@@ -0,0 +1,14 @@
+"Demonstrate a trivial unbounded variable"
+from gpkit import Variable, Model
+from gpkit.constraints.bounded import BoundedConstraintSet
+
+x = Variable("x")
+
+constraints = [x >= 1]
+
+# Model(1/x, constraints).solve()  # does not solve
+m = Model(1/x, BoundedConstraintSet(constraints))
+# by default, prints bounds warning during solve
+sol = m.solve(verbosity=0)
+print sol.table()
+print "sol['boundedness'] is:", sol["boundedness"]
diff --git a/docs/source/examples/unbounded_output.txt b/docs/source/examples/unbounded_output.txt
new file mode 100644
index 0000000..1daa5f8
--- /dev/null
+++ b/docs/source/examples/unbounded_output.txt
@@ -0,0 +1,15 @@
+
+UNBOUNDED VARIABLES
+   value near upper bound: [x]
+ sensitive to upper bound: [x]
+
+
+Cost
+----
+ 1e-30
+
+Free Variables
+--------------
+x : 1e+30 
+
+sol['boundedness'] is: {'value near upper bound': array([x], dtype=object), 'sensitive to upper bound': array([x], dtype=object)}
diff --git a/docs/source/gettingstarted.rst b/docs/source/gettingstarted.rst
index e2b477b..2574a32 100644
--- a/docs/source/gettingstarted.rst
+++ b/docs/source/gettingstarted.rst
@@ -158,12 +158,16 @@ We can also manually print the solution table, with the same result as if the ve
     -------------
     S : -1.5
 
+We can also print the optimal value and solved variables individually.
+
 .. code-block:: python
 
+    print "The optimal value is %s." % (sol["cost"])
     print "The x dimension is %s." % (sol(x))
 
 ::
 
+    The optimal value is 15.5884619886.
     The x dimension is 0.577351209028 meter.
 
 .. refactor this section; explain what can be done with a SolutionArray
diff --git a/gpkit/_mosek/expopt.py b/gpkit/_mosek/expopt.py
index 920552f..99d8e84 100644
--- a/gpkit/_mosek/expopt.py
+++ b/gpkit/_mosek/expopt.py
@@ -17,7 +17,7 @@ from ctypes import POINTER as ptr_factory
 from ctypes import c_double, c_int, c_void_p
 
 try:
-    from .lib import expopt_h
+    from .lib import expopt_h # pylint: disable=import-error
 except Exception as e:
     raise ImportError("Could not load MOSEK library: "+repr(e))
 
@@ -59,6 +59,10 @@ MSK = ModuleShortener("MSK", expopt_h)
 
 # lookup table from mosek.h, "MSKsolsta_enum"
 # switched where indicated because MOSEK solves the dual GP problem
+
+# Pylint disabled here because the _ naming convention is a
+# side effect of the ModuleShortener
+# pylint: disable=invalid-name,attribute-defined-outside-init,protected-access
 MSK._SOL_STA_LOOKUPTABLE = ["UNKNOWN",
                             "OPTIMAL",
                             "DUAL_FEAS",  # originally position 3
@@ -119,7 +123,7 @@ def printcb(void, msg):  # pylint: disable=unused-argument
     return 0
 
 
-# pylint: disable=unused-argument,too-many-locals
+# pylint: disable=unused-argument,too-many-locals,protected-access
 def imize(c, A, p_idxs, *args, **kwargs):
     """Interface to the MOSEK EXPOPT solver via C
 
diff --git a/gpkit/build.py b/gpkit/build.py
index ede2c0d..5b18074 100644
--- a/gpkit/build.py
+++ b/gpkit/build.py
@@ -14,7 +14,7 @@ settings = {}
 
 def log(*args):
     "Print a line and append it to the log string."
-    global LOGSTR
+    global LOGSTR  # pylint: disable=global-statement
     print(*args)
     LOGSTR += " ".join(args) + "\n"
 
@@ -111,7 +111,8 @@ class MosekCLI(SolverBackend):
             log("#   Trying to run mskexpopt...")
             if call("mskexpopt") in (1052, 28):  # 28 for MacOSX
                 return "in system path"
-        except Exception:
+        except:  # pylint: disable=bare-except
+            # exception type varies by operating system
             return
 
 
@@ -123,7 +124,8 @@ class CVXopt(SolverBackend):
         "Attempts to import mskexpopt."
         try:
             log("#   Trying to import cvxopt...")
-            import cvxopt
+            # Testing the import, so the variable is intentionally not used
+            import cvxopt # pylint: disable=unused-variable
             return "in Python path"
         except ImportError:
             return
@@ -152,6 +154,13 @@ class Mosek(SolverBackend):
         }
     }
 
+    def __init__(self):
+        super(Mosek, self).__init__()
+        self.expopt_files = None
+        self.bin_dir = None
+        self.flags = None
+        self.lib_path = None
+
     def look(self):
         "Looks in default install locations for latest mosek version."
         if sys.platform == "win32":
@@ -202,9 +211,7 @@ class Mosek(SolverBackend):
         self.bin_dir = pathjoin(lib_dir, "bin")
         self.lib_path = glob.glob(self.bin_dir+os.sep+libpattern)[0]
 
-        if not isfile(h_path):
-            return
-        if not isfile(self.lib_path):
+        if not isfile(h_path) or not isfile(self.lib_path):
             return
 
         expopt_dir = pathjoin(tools_dir, "examples", "c")
@@ -217,7 +224,7 @@ class Mosek(SolverBackend):
             if not isfile(expopt_file):
                 return
 
-        global settings
+        global settings  # pylint: disable=global-variable-not-assigned
         settings["mosek_bin_dir"] = self.bin_dir
         os.environ['PATH'] = os.environ['PATH'] + os.pathsep + self.bin_dir
 
@@ -226,7 +233,8 @@ class Mosek(SolverBackend):
     def build(self):
         "Builds a dynamic library to GPKITBUILD or $HOME/.gpkit"
         try:
-            import ctypesgencore
+            # Testing the import, so the variable is intentionally not used
+            import ctypesgencore # pylint: disable=unused-variable
         except ImportError:
             log("## SKIPPING MOSEK INSTALL: CTYPESGENCORE WAS NOT FOUND")
             return
@@ -286,7 +294,7 @@ class Mosek(SolverBackend):
 
 def build_gpkit():
     "Builds GPkit"
-    global settings
+    global settings  # pylint: disable=global-variable-not-assigned
 
     if isfile("__init__.py"):
         #call("ls")
diff --git a/gpkit/constraints/bounded.py b/gpkit/constraints/bounded.py
new file mode 100644
index 0000000..e804bc9
--- /dev/null
+++ b/gpkit/constraints/bounded.py
@@ -0,0 +1,103 @@
+"Implements BoundedConstraintSet"
+from collections import defaultdict
+import numpy as np
+
+from .. import Variable
+from .set import ConstraintSet
+from ..small_scripts import mag
+
+
+def varkey_bounds(varkeys, lower, upper):
+    """Returns constraints list bounding all varkeys.
+
+    Arguments
+    ---------
+    varkeys : iterable
+        list of varkeys to create bounds for
+
+    lower : float
+        lower bound for all varkeys
+
+    upper : float
+        upper bound for all varkeys
+    """
+    constraints = []
+    for varkey in varkeys:
+        variable = Variable(**varkey.descr)
+        units = varkey.descr.get("units", 1)
+        constraints.append([upper*units >= variable,
+                            variable >= lower*units])
+    return constraints
+
+
+class BoundedConstraintSet(ConstraintSet):
+    """Bounds contained variables so as to ensure dual feasibility.
+
+    Arguments
+    ---------
+    constraints : iterable
+        constraints whose varkeys will be bounded
+
+    substitutions : dict
+        as in ConstraintSet.__init__
+
+    verbosity : int
+        how detailed of a warning to print
+            0: nothing
+            1: print warnings
+
+    eps : float
+        default lower bound is eps, upper bound is 1/eps
+
+    lower : float
+        lower bound for all varkeys, replaces eps
+
+    upper : float
+        upper bound for all varkeys, replaces 1/eps
+    """
+
+    def __init__(self, constraints, substitutions=None, verbosity=1,
+                 eps=1e-30, lower=None, upper=None):
+        self.bound_las = None
+        self.verbosity = verbosity
+        self.lowerbound = lower if lower else eps
+        self.upperbound = upper if upper else 1/eps
+        constraintset = ConstraintSet([constraints])
+        self.bounded_varkeys = tuple(vk for vk in constraintset.varkeys
+                                     if "value" not in vk.descr)
+        bounding_constraints = varkey_bounds(self.bounded_varkeys,
+                                             self.lowerbound, self.upperbound)
+        constraints = [constraintset, bounding_constraints]
+        super(BoundedConstraintSet, self).__init__(constraints, substitutions)
+
+    def sens_from_dual(self, las, nus):
+        "Return sensitivities while capturing the relevant lambdas"
+        self.bound_las = las[-2*len(self.bounded_varkeys):]
+        return super(BoundedConstraintSet, self).sens_from_dual(las, nus)
+
+    def process_result(self, result):
+        "Creates (and potentially prints) a dictionary of unbounded variables."
+        lam = self.bound_las
+        out = defaultdict(list)
+        for i, varkey in enumerate(self.bounded_varkeys):
+            lam_gt, lam_lt = lam[2*i], lam[2*i+1]
+            if abs(lam_gt) >= 1e-7:  # arbitrary threshold
+                out["sensitive to upper bound"].append(varkey)
+            if abs(lam_lt) >= 1e-7:  # arbitrary threshold
+                out["sensitive to lower bound"].append(varkey)
+            value = mag(result["variables"][varkey])
+            distance_below = np.log(value/self.lowerbound)
+            distance_above = np.log(self.upperbound/value)
+            if distance_below <= 3:  # arbitrary threshold
+                out["value near lower bound"].append(varkey)
+            elif distance_above <= 3:  # arbitrary threshold
+                out["value near upper bound"].append(varkey)
+        if self.verbosity > 0 and out:
+            print
+            print "UNBOUNDED VARIABLES"
+            for key, value in out.items():
+                print "% 25s: %s" % (key, value)
+            print
+        if not "boundedness" in result:
+            result["boundedness"] = {}
+        result["boundedness"].update(out)
diff --git a/gpkit/constraints/costed.py b/gpkit/constraints/costed.py
index 70901a9..011c279 100644
--- a/gpkit/constraints/costed.py
+++ b/gpkit/constraints/costed.py
@@ -14,7 +14,7 @@ class CostedConstraintSet(ConstraintSet):
     """
     def __init__(self, cost, constraints, substitutions=None):
         self.cost = cost
-        subs = self.cost.values
+        subs = dict(self.cost.values)
         if substitutions:
             subs.update(substitutions)
         ConstraintSet.__init__(self, constraints, subs)
diff --git a/gpkit/constraints/set.py b/gpkit/constraints/set.py
index 07ed370..b9b4c84 100644
--- a/gpkit/constraints/set.py
+++ b/gpkit/constraints/set.py
@@ -1,16 +1,12 @@
 "Implements ConstraintSet"
+import numpy as np
+
 from ..small_classes import HashVector
 from ..keydict import KeySet, KeyDict
 from ..small_scripts import try_str_without
 from ..repr_conventions import _str, _repr, _repr_latex_
 
 
-def _sort_by_num_models(var):
-    "return integer for Variable sorting"
-    mods = var.key.models
-    return len(mods) if mods else 0
-
-
 def _sort_by_name_and_idx(var):
     "return tuplef for Variable sorting"
     return (var.key.str_without(["units", "idx"]), var.key.idx)
@@ -22,7 +18,7 @@ class ConstraintSet(list):
         if isinstance(constraints, ConstraintSet):
             constraints = [constraints]
         list.__init__(self, constraints)
-        subs = substitutions if substitutions else {}
+        subs = dict(substitutions) if substitutions else {}
         self.unused_variables = None
         if not isinstance(constraints, ConstraintSet):
             # constraintsetify everything
@@ -60,14 +56,27 @@ class ConstraintSet(list):
             return list.__getitem__(self, key)
         else:
             from ..nomials import Variable
-            variables = [Variable(**key.descr) for key in self.varkeys[key]]
+            variables = [Variable(**k.descr) for k in self.varkeys[key]]
             if len(variables) == 1:
                 return variables[0]
-            else:
-                variables.sort(key=_sort_by_num_models)
-                variable = variables[0]
-                # note: doesn't work for vector variables
-                return variable
+            elif variables[0].key.idx and variables[0].key.shape:
+                # maybe it's all one vector variable!
+                vector = np.full(variables[0].key.shape, np.nan, dtype="object")
+                goaldict = dict(variables[0].key.descr)
+                del goaldict["idx"]
+                for variable in variables:
+                    testdict = dict(variable.key.descr)
+                    del testdict["idx"]
+                    if testdict == goaldict:
+                        vector[variable.key.idx] = variable
+                    else:
+                        vector = None
+                        break
+                if vector is not None:
+                    return vector
+            raise ValueError("multiple variables are called '%s'; use"
+                             " variable_byname('%s') to see all of them"
+                             % (key, key))
 
     def variables_byname(self, key):
         "Get all variables with a given name"
diff --git a/gpkit/constraints/signomial_program.py b/gpkit/constraints/signomial_program.py
index 54dadb9..fdcbc42 100644
--- a/gpkit/constraints/signomial_program.py
+++ b/gpkit/constraints/signomial_program.py
@@ -3,6 +3,7 @@ from time import time
 from ..geometric_program import GeometricProgram
 from ..feasibility import feasibility_model
 from .costed import CostedConstraintSet
+from ..exceptions import InvalidGPConstraint
 
 
 class SignomialProgram(CostedConstraintSet):
@@ -48,10 +49,10 @@ class SignomialProgram(CostedConstraintSet):
         try:
             _ = self.as_posyslt1()  # should raise an error
             # TODO: is there a faster way to check?
-        except TypeError:
+        except InvalidGPConstraint:
             pass
         else:  # this is a GP
-            raise ValueError("""No Signomials remained after substitution.
+            raise ValueError("""Model valid as a Geometric Program.
 
     SignomialPrograms should only be created with Models containing Signomial
     Constraints, since Models without Signomials have global solutions and can
@@ -103,6 +104,7 @@ class SignomialProgram(CostedConstraintSet):
     appear to be converging, you may wish to increase the iteration limit by
     calling .localsolve(..., iteration_limit=NEWLIMIT).""" % len(self.gps))
             gp = self.gp(x0, verbosity=verbosity-1)
+            gp.x0 = x0           # NOTE: SIDE EFFECTS
             self.gps.append(gp)  # NOTE: SIDE EFFECTS
             try:
                 result = gp.solve(solver, verbosity-1, **kwargs)
diff --git a/gpkit/exceptions.py b/gpkit/exceptions.py
new file mode 100644
index 0000000..fa6e56c
--- /dev/null
+++ b/gpkit/exceptions.py
@@ -0,0 +1,4 @@
+"GPkit-specific Exception classes"
+class InvalidGPConstraint(Exception):
+    "Raised when a non-GP-compatible constraint is used in a GP"
+    pass
diff --git a/gpkit/geometric_program.py b/gpkit/geometric_program.py
index 5f4b592..9df2b29 100644
--- a/gpkit/geometric_program.py
+++ b/gpkit/geometric_program.py
@@ -10,6 +10,7 @@ from .small_classes import SolverLog
 
 DEFAULT_SOLVER_KWARGS = {"cvxopt": {"kktsolver": "ldl"}}
 
+
 class GeometricProgram(NomialData):
     # pylint: disable=too-many-instance-attributes
     """Standard mathematical representation of a GP.
@@ -66,16 +67,7 @@ class GeometricProgram(NomialData):
 
         for constraint in constraints:
             constraint.substitutions = self.substitutions
-            try:
-                constr_posys = constraint.as_posyslt1()
-            except TypeError as err:
-                if err.message == ("SignomialInequality could not simplify to"
-                                   " a PosynomialInequality"):
-                    raise ValueError("GeometricPrograms cannot contain"
-                                     " SignomialInequalities: try forming your"
-                                     " program as SignomialProgram or calling"
-                                     " .localsolve().")
-                raise
+            constr_posys = constraint.as_posyslt1()
             if not all(constr_posys):
                 raise ValueError("%s is an invalid constraint for a"
                                  " GeometricProgram" % constraint)
diff --git a/gpkit/interactive/plotting.py b/gpkit/interactive/plotting.py
index 78c425b..a37c41f 100644
--- a/gpkit/interactive/plotting.py
+++ b/gpkit/interactive/plotting.py
@@ -1,7 +1,7 @@
 """Plotting methods"""
-import numpy as np
 import matplotlib.pyplot as plt
 from matplotlib.mlab import griddata
+import numpy as np
 from .. import VarKey
 from ..small_scripts import unitstr, mag
 
@@ -37,7 +37,7 @@ def contour_array(model, xname, yname, znames, cellsize=(5, 5),
     # pylint: disable=too-many-arguments
     # pylint: disable=too-many-locals
     # pylint: disable=too-many-statements
-    """Arrar of contour plots for variables in a solved model
+    """Array of contour plots for variables in a solved model
 
     Arguments
     ---------
@@ -84,7 +84,8 @@ def contour_array(model, xname, yname, znames, cellsize=(5, 5),
         colors = LIGHT_COLORS[0], DARK_COLORS[0]
 
     fig, axes = plt.subplots(nrows, ncols, figsize=figsize,
-                             sharex=True, sharey=True)
+                             sharex=True, sharey=True, squeeze=False)
+    axes = axes[0]
     if nrows > 1 and ncols > 1:
         xlabeledaxes = axes[-1, :]
         ylabeledaxes = axes[:, 0]
@@ -92,7 +93,7 @@ def contour_array(model, xname, yname, znames, cellsize=(5, 5),
         xlabeledaxes = axes
         ylabeledaxes = [axes[0]]
     else:
-        xlabeledaxes = [axes[-1]]
+        xlabeledaxes = axes
         ylabeledaxes = axes
 
     for ax in xlabeledaxes:
@@ -138,19 +139,16 @@ def contour_array(model, xname, yname, znames, cellsize=(5, 5),
 #     pass
 
 
-def plot_frontiers(gp, znames, nrow=1, ncol=3, figsize=(15, 5)):
+def plot_frontiers(model, znames, nrow=1, ncol=3, figsize=(15, 5)):
     "Helper function to plot 2d contour plots."
-    sol = gp.solution
-    data = dict(sol["variables"])
-    data.update({"S{%s}" % k: v
-                 for (k, v) in sol["sensitivities"]["constants"].items()})
-    if len(gp.sweep) == 2:
-        contour_array(data,
-                      gp.sweep.keys()[0],
-                      gp.sweep.keys()[1],
-                      znames, nrow, ncol, figsize,
-                      xticks=gp.sweep.values()[0],
-                      yticks=gp.sweep.values()[1])
+    sweeps = model.solution['sweepvariables']
+    if len(sweeps) == 2:
+        contour_array(model,
+                      sweeps.keys()[0],
+                      sweeps.keys()[1],
+                      znames, figsize, nrow, ncol,
+                      xticks=sweeps.values()[0],
+                      yticks=sweeps.values()[1])
 
 
 def _combine_nearby_ticks(ticks, lim, ntick):
diff --git a/gpkit/interactive/ractor.py b/gpkit/interactive/ractor.py
index 034e84f..c717f30 100644
--- a/gpkit/interactive/ractor.py
+++ b/gpkit/interactive/ractor.py
@@ -39,14 +39,25 @@ def ractorjs(title, model, update_py, ranges, constraint_js=""):
     widget_id = "jswidget_"+str(JSWIDGET_ID())
     display(HTML("<script id='%s-after' type='text/throwaway'>%s</script>" %
                  (widget_id, constraint_js)))
-    display(HTML("<script>var %s = {storage: [], n:%i, ranges: {}, after: document.getElementById('%s-after').innerHTML, bases: [1] }</script>" % (widget_id, len(ranges), widget_id)))
+    display(HTML(
+        ("<script>var %s = {storage: [], n:%i, ranges: {},"
+         " after: document.getElementById('%s-after').innerHTML, bases: [1] }"
+         "</script>") %
+        (widget_id, len(ranges), widget_id)))
 
     container_id = widget_id + "_container"
-    display(HTML("<div id='%s'></div><style>#%s td {text-align: right; border: none !important;}\n#%s tr {border: none !important;}\n#%s table {border: none !important;}\n</style>" % (container_id, container_id, container_id, container_id)))
+    display(HTML(
+        ("<div id='%s'></div><style>"
+         "#%s td {text-align: right; border: none !important;}\n"
+         "#%s tr {border: none !important;}\n"
+         "#%s table {border: none !important;}\n</style>") %
+        (container_id, container_id, container_id, container_id)))
 
     template_id = widget_id + "_template"
     template = '<script id="%s" type="text/ractive"><table>' % template_id
-    ctrl_template = Template('<tr><td>$var</td><input value="{{$varname}}" type="range" min="0" max="$len" step="1"><td><span id="$w-$varname"></span></td></tr>\n')
+    ctrl_template = Template(('<tr><td>$var</td><input value="{{$varname}}" '
+                              'type="range" min="0" max="$len" step="1"><td>'
+                              '<span id="$w-$varname"></span></td></tr>\n'))
     data_init = ""
 
     subs = {}
@@ -94,7 +105,9 @@ def ractorjs(title, model, update_py, ranges, constraint_js=""):
 
     display(HTML(template + "</table></script>"))
 
-    loader = Template("""getScript('http://cdn.ractivejs.org/0.6.1/ractive.min.js', function() {
+    loader = Template(
+        """getScript('http://cdn.ractivejs.org/0.6.1/ractive.min.js',
+          function() {
           $w.ractive = new Ractive({
           el: '$container_id',
           template: '#$template_id',
@@ -106,7 +119,8 @@ def ractorjs(title, model, update_py, ranges, constraint_js=""):
                   varname = 'var'+i
                   idx = $w.ractive.data[varname]
                   console.log("$w-"+varname);
-                  document.getElementById("$w-"+varname).innerText = Math.round(100*$w.ranges[varname][idx])/100
+                  percent = Math.round(100*$w.ranges[varname][idx])/100
+                  document.getElementById("$w-"+varname).innerText = percent
                   idxsum += idx*$w.bases[i]
               }
               if ($w.storage[idxsum] === "") {
diff --git a/gpkit/interactive/sensitivity_map.py b/gpkit/interactive/sensitivity_map.py
index 6afbd59..c16b85a 100644
--- a/gpkit/interactive/sensitivity_map.py
+++ b/gpkit/interactive/sensitivity_map.py
@@ -1,7 +1,7 @@
 "Implements heatmapped equations to highlight sensitivities."
 import numpy as np
-from gpkit.small_scripts import mag, latex_num
 from gpkit.nomials import MonomialEquality
+from gpkit.small_scripts import nomial_latex_helper
 
 
 BLUE = np.array([16, 131, 246])/255.0
@@ -24,7 +24,6 @@ def colorfn_gen(scale, power=0.66):
         return "[rgb]{%.2f,%.2f,%.2f}" % tuple(blended_color)
     return colorfn
 
-
 # pylint: disable=too-many-locals
 def signomial_print(sig, sol, colorfn, paintby="constants", idx=None):
     "For pretty printing with Sympy"
@@ -44,32 +43,7 @@ def signomial_print(sig, sol, colorfn, paintby="constants", idx=None):
             elif x < 0:
                 neg_vars.append((varlatex, x))
 
-        pvarstrs = ['%s^{%.2g}' % (varl, x) if "%.2g" % x != "1" else varl
-                    for (varl, x) in pos_vars]
-        nvarstrs = ['%s^{%.2g}' % (varl, -x)
-                    if "%.2g" % -x != "1" else varl
-                    for (varl, x) in neg_vars]
-        pvarstrs.sort()
-        nvarstrs.sort()
-        pvarstr = ' '.join(pvarstrs)
-        nvarstr = ' '.join(nvarstrs)
-        c = mag(c)
-        cstr = "%.2g" % c
-        if pos_vars and (cstr == "1" or cstr == "-1"):
-            cstr = cstr[:-1]
-        else:
-            cstr = latex_num(c)
-
-        if not pos_vars and not neg_vars:
-            mstr = "%s" % cstr
-        elif pos_vars and not neg_vars:
-            mstr = "%s%s" % (cstr, pvarstr)
-        elif neg_vars and not pos_vars:
-            mstr = "\\frac{%s}{%s}" % (cstr, nvarstr)
-        elif pos_vars and neg_vars:
-            mstr = "%s\\frac{%s}{%s}" % (cstr, pvarstr, nvarstr)
-
-        mstrs.append(mstr)
+        mstrs.append(nomial_latex_helper(c, pos_vars, neg_vars))
 
     if paintby == "monomials":
         mstrs_ = []
diff --git a/gpkit/interactive/widgets.py b/gpkit/interactive/widgets.py
index fbec2d3..588cd70 100644
--- a/gpkit/interactive/widgets.py
+++ b/gpkit/interactive/widgets.py
@@ -108,7 +108,7 @@ def modelcontrolpanel(model, *args, **kwargs):
         cb = widgets.Checkbox(value=True)
         unit_latex = sl.varkey.unitstr()
         if unit_latex:
-            unit_latex = "$\scriptsize"+unit_latex+"$"
+            unit_latex = r"$\scriptsize"+unit_latex+"$"
         units = widgets.Latex(value=unit_latex)
         units.font_size = "1.16em"
         box = widgets.HBox(children=[cb, sl, units])
@@ -146,8 +146,8 @@ def modelcontrolpanel(model, *args, **kwargs):
 
 
 def create_settings(box):
-    "Creates a widget Container for settings and info of  a particular slider."
-    sl_enable, slider, sl_units = box.children
+    "Creates a widget Container for settings and info of a particular slider."
+    _, slider, sl_units = box.children
 
     enable = widgets.Checkbox(value=box.visible)
     link((box, 'visible'), (enable, 'value'))
diff --git a/gpkit/modified_ctypesgen.py b/gpkit/modified_ctypesgen.py
index f280c8d..d341b48 100644
--- a/gpkit/modified_ctypesgen.py
+++ b/gpkit/modified_ctypesgen.py
@@ -1,3 +1,4 @@
+# pylint: skip-file
 # from https://github.com/kanzure/ctypesgen
 #
 # Copyright (c) 2007-2008, Ctypesgen Developers
diff --git a/gpkit/nomials/nomial_core.py b/gpkit/nomials/nomial_core.py
index 776fe57..f9bf680 100644
--- a/gpkit/nomials/nomial_core.py
+++ b/gpkit/nomials/nomial_core.py
@@ -1,7 +1,7 @@
 "The shared non-mathematical backbone of all Nomials"
 from .data import NomialData
 from ..small_classes import Numbers
-from ..small_scripts import latex_num
+from ..small_scripts import nomial_latex_helper
 from ..small_scripts import mag, unitstr
 from ..repr_conventions import _str, _repr, _repr_latex_
 
@@ -62,7 +62,6 @@ class Nomial(NomialData):
         return " + ".join(sorted(mstrs)) + units
 
     def latex(self, excluded=None):
-        # pylint: disable=too-many-locals
         "For pretty printing with Sympy"
         if excluded is None:
             excluded = []
@@ -75,30 +74,7 @@ class Nomial(NomialData):
                 elif x < 0:
                     neg_vars.append((var.latex(excluded), x))
 
-            pvarstrs = ['%s^{%.2g}' % (varl, x) if "%.2g" % x != "1" else varl
-                        for (varl, x) in pos_vars]
-            nvarstrs = ['%s^{%.2g}' % (varl, -x)
-                        if "%.2g" % -x != "1" else varl
-                        for (varl, x) in neg_vars]
-            pvarstrs.sort()
-            nvarstrs.sort()
-            pvarstr = ' '.join(pvarstrs)
-            nvarstr = ' '.join(nvarstrs)
-            c = mag(c)
-            cstr = "%.2g" % c
-            if pos_vars and (cstr == "1" or cstr == "-1"):
-                cstr = cstr[:-1]
-            else:
-                cstr = latex_num(c)
-
-            if not pos_vars and not neg_vars:
-                mstrs.append("%s" % cstr)
-            elif pos_vars and not neg_vars:
-                mstrs.append("%s%s" % (cstr, pvarstr))
-            elif neg_vars and not pos_vars:
-                mstrs.append("\\frac{%s}{%s}" % (cstr, nvarstr))
-            elif pos_vars and neg_vars:
-                mstrs.append("%s\\frac{%s}{%s}" % (cstr, pvarstr, nvarstr))
+            mstrs.append(nomial_latex_helper(c, pos_vars, neg_vars))
 
         if "units" in excluded:
             return " + ".join(sorted(mstrs))
diff --git a/gpkit/nomials/nomial_math.py b/gpkit/nomials/nomial_math.py
index f6df7ad..3da2458 100644
--- a/gpkit/nomials/nomial_math.py
+++ b/gpkit/nomials/nomial_math.py
@@ -12,6 +12,7 @@ from ..varkey import VarKey
 from ..small_scripts import mag
 from .. import units as ureg
 from .. import DimensionalityError
+from ..exceptions import InvalidGPConstraint
 
 
 class Signomial(Nomial):
@@ -688,18 +689,34 @@ class SignomialInequality(ScalarSingleEquationConstraint):
             return self.unsubbed   # pylint: disable=no-member
 
         else:
-            raise TypeError("SignomialInequality could not simplify to"
-                            " a PosynomialInequality")
+            raise InvalidGPConstraint("SignomialInequality could not simplify"
+                                      " to a PosynomialInequality; try forming"
+                                      " your program as SignomialProgram or"
+                                      " calling .localsolve().")
+
+    def _fill_default_x0(self, x0, varkeys):
+        """For all keys in varkeys, updates x0 with default values and
+        substitutions. Returns x0.
+
+        Order of precedence for x0 default:
+            - substitution value
+            - x0 value
+            - sp_init value
+            - 1.0
+        """
+        if x0 is None:
+            x0 = {}
+        x0.update({
+            vk: self.substitutions.get(vk, x0.get(vk, vk.descr.get("sp_init",
+                                                                   1.0)))
+            for vk in varkeys})
+        return x0
 
     def as_gpconstr(self, x0):
         "Returns GP approximation of an SP constraint at x0"
         siglt0, = self.unsubbed
         posy, negy = siglt0.posy_negy()
-        if x0 is None:
-            x0 = {vk: vk.descr["sp_init"] for vk in negy.varlocs
-                  if "sp_init" in vk.descr}
-        x0.update({var: 1 for var in negy.varlocs if var not in x0})
-        x0.update(self.substitutions)
+        x0 = self._fill_default_x0(x0, negy.varlocs)
         pc = PosynomialInequality(posy, "<=", negy.mono_lower_bound(x0))
         pc.substitutions = self.substitutions
         return pc
@@ -715,17 +732,15 @@ class SignomialEquality(SignomialInequality):
     def as_posyslt1(self):
         "Returns the posys <= 1 representation of this constraint."
         # todo deal with substitutions
-        raise TypeError("SignomialEquality could not simplify to"
-                        " a PosynomialInequality")
+        raise InvalidGPConstraint("SignomialEquality could not simplify"
+                                  " to a PosynomialInequality; try forming"
+                                  " your program as SignomialProgram or"
+                                  " calling .localsolve().")
 
     def as_gpconstr(self, x0):
         "Returns GP approximation of an SP constraint at x0"
         siglt0, = self.unsubbed
-        if x0 is None:
-            x0 = {vk: vk.descr["sp_init"] for vk in siglt0.varlocs
-                  if "sp_init" in vk.descr}
-        x0.update({var: 1 for var in siglt0.varlocs if var not in x0})
-        x0.update(self.substitutions)
+        x0 = self._fill_default_x0(x0, siglt0.varlocs)
         posy, negy = siglt0.posy_negy()
         mec = (posy.mono_lower_bound(x0) == negy.mono_lower_bound(x0))
         mec.substitutions = self.substitutions
diff --git a/gpkit/nomials/substitution.py b/gpkit/nomials/substitution.py
index f949586..50b5ef4 100644
--- a/gpkit/nomials/substitution.py
+++ b/gpkit/nomials/substitution.py
@@ -116,6 +116,7 @@ def substitution(nomial, substitutions):
         cs_ = Quantity(cs_, nomial.cs.units)
     varlocs_ = defaultdict(list)
     varlocs_.update({vk: list(idxs) for (vk, idxs) in nomial.varlocs.items()})
+    # pylint: disable=too-many-nested-blocks
     for var, sub in subs.items():
         for i in nomial.varlocs[var]:
             x = exps_[i].pop(var)
diff --git a/gpkit/nomials/variables.py b/gpkit/nomials/variables.py
index 46c443e..bbb992e 100644
--- a/gpkit/nomials/variables.py
+++ b/gpkit/nomials/variables.py
@@ -102,10 +102,7 @@ class ArrayVariable(NomialArray):
         if "idx" in descr:
             raise KeyError("the description field 'idx' is reserved")
 
-        if isinstance(shape, Numbers):
-            shape = (shape,)
-        else:
-            shape = tuple(shape)
+        shape = (shape,) if isinstance(shape, Numbers) else tuple(shape)
 
         descr["shape"] = shape
 
@@ -120,7 +117,7 @@ class ArrayVariable(NomialArray):
             elif isinstance(arg, Strings) and "label" not in descr:
                 descr["label"] = arg
 
-        values = []
+        values, value_option = [], None
         for value_option in ["value", "sp_init"]:
             if value_option in descr:
                 values = descr.pop(value_option)
@@ -142,6 +139,7 @@ class ArrayVariable(NomialArray):
         if "name" not in descr:
             descr["name"] = "\\fbox{%s}" % VarKey.new_unnamed_id()
 
+        arraykey = VarKey(**descr)
         vl = np.empty(shape, dtype="object")
         it = np.nditer(vl, flags=['multi_index', 'refs_ok'])
         while not it.finished:
@@ -153,6 +151,7 @@ class ArrayVariable(NomialArray):
             elif valuetype == "list":
                 descr.update({value_option: values[i[0]]})
             vl[i] = Variable(**descr)
+            vl[i].key.arraykey = arraykey
 
         obj = np.asarray(vl).view(cls)
         obj.descr = descr
diff --git a/gpkit/small_classes.py b/gpkit/small_classes.py
index d853d81..7c3fc10 100644
--- a/gpkit/small_classes.py
+++ b/gpkit/small_classes.py
@@ -80,6 +80,7 @@ class DictOfLists(dict):
         "Appends a dict (of dicts) of lists to all held lists."
         if not hasattr(self, 'initialized'):
             _enlist_dict(sol, self)
+            # pylint: disable=attribute-defined-outside-init
             self.initialized = True
         else:
             _append_dict(sol, self)
diff --git a/gpkit/small_scripts.py b/gpkit/small_scripts.py
index e11087b..d9c46f4 100644
--- a/gpkit/small_scripts.py
+++ b/gpkit/small_scripts.py
@@ -27,19 +27,6 @@ def listify(item):
         return [item]
 
 
-def isequal(a, b):
-    "Determine if two elements are equal, recursing through Iterables"
-    # pylint: disable=invalid-name
-    if (isinstance(a, Iterable) and
-            not isinstance(a, Strings+(tuple, list, dict))):
-        for i, a_i in enumerate(a):
-            if not isequal(a_i, b[i]):
-                return False
-    elif a != b:
-        return False
-    return True
-
-
 def mag(c):
     "Return magnitude of a Number or Quantity"
     if isinstance(c, Quantity):
@@ -67,24 +54,49 @@ def unitstr(units, into="%s", options="~", dimless='-'):
         return ""
 
 
+def nomial_latex_helper(c, pos_vars, neg_vars):
+    """Combines (varlatex, exponent) tuples,
+    separated by positive vs negative exponent,
+    into a single latex string"""
+    # TODO this is awkward due to sensitivity_map, which needs a refactor
+    pvarstrs = ['%s^{%.2g}' % (varl, x) if "%.2g" % x != "1" else varl
+                for (varl, x) in pos_vars]
+    nvarstrs = ['%s^{%.2g}' % (varl, -x)
+                if "%.2g" % -x != "1" else varl
+                for (varl, x) in neg_vars]
+    pvarstrs.sort()
+    nvarstrs.sort()
+    pvarstr = ' '.join(pvarstrs)
+    nvarstr = ' '.join(nvarstrs)
+    c = mag(c)
+    cstr = "%.2g" % c
+    if pos_vars and (cstr == "1" or cstr == "-1"):
+        cstr = cstr[:-1]
+    else:
+        cstr = latex_num(c)
+
+    if not pos_vars and not neg_vars:
+        mstr = "%s" % cstr
+    elif pos_vars and not neg_vars:
+        mstr = "%s%s" % (cstr, pvarstr)
+    elif neg_vars and not pos_vars:
+        mstr = "\\frac{%s}{%s}" % (cstr, nvarstr)
+    elif pos_vars and neg_vars:
+        mstr = "%s\\frac{%s}{%s}" % (cstr, pvarstr, nvarstr)
+
+    return mstr
+
+
 def is_sweepvar(sub):
     "Determines if a given substitution indicates a sweep."
     try:
         if sub[0] == "sweep":
             if isinstance(sub[1], Iterable) or hasattr(sub[1], "__call__"):
                 return True
-    except:
+    except (TypeError, IndexError):
         return False
 
 
-def invalid_types_for_oper(oper, a, b):
-    "Raises TypeError for unsupported operations."
-    typea = a.__class__.__name__
-    typeb = b.__class__.__name__
-    raise TypeError("unsupported operand types"
-                    " for %s: '%s' and '%s'" % (oper, typea, typeb))
-
-
 def latex_num(c):
     "Returns latex string of numbers, potentially using exponential notation."
     cstr = "%.4g" % c
diff --git a/gpkit/tools/tools.py b/gpkit/tools/tools.py
index 3eb1308..1fa4b70 100644
--- a/gpkit/tools/tools.py
+++ b/gpkit/tools/tools.py
@@ -2,7 +2,6 @@
 import numpy as np
 from ..nomials import Variable, VectorVariable
 from ..nomials import NomialArray
-# pylint: disable=duplicate-code
 
 
 def te_exp_minus1(posy, nterm):


Test patch:
diff --git a/gpkit/tests/helpers.py b/gpkit/tests/helpers.py
index f253a71..1e01011 100644
--- a/gpkit/tests/helpers.py
+++ b/gpkit/tests/helpers.py
@@ -90,7 +90,7 @@ def run_tests(tests, xmloutput=None, verbosity=2):
     for t in tests:
         suite.addTests(loader.loadTestsFromTestCase(t))
     if xmloutput:
-        import xmlrunner
+        import xmlrunner  # pylint: disable=import-error
         xmlrunner.XMLTestRunner(output=xmloutput).run(suite)
     else:
         unittest.TextTestRunner(verbosity=verbosity).run(suite)
diff --git a/gpkit/tests/t_constraints.py b/gpkit/tests/t_constraints.py
index 67e523e..6e7d3d2 100644
--- a/gpkit/tests/t_constraints.py
+++ b/gpkit/tests/t_constraints.py
@@ -6,12 +6,29 @@ from gpkit.nomials import MonomialEquality
 from gpkit import LinkedConstraintSet, Model
 from gpkit.constraints.tight import TightConstraintSet
 from gpkit.tests.helpers import run_tests
+from gpkit.exceptions import InvalidGPConstraint
 import gpkit
 
 
 class TestConstraint(unittest.TestCase):
     """Tests for Constraint class"""
 
+    def test_constraintget(self):
+        x = Variable("x")
+        x_ = Variable("x", model="_")
+        xv = VectorVariable(2, "x")
+        xv_ = VectorVariable(2, "x", model="_")
+        self.assertEqual(Model(x, [x >= 1])["x"], x)
+        with self.assertRaises(ValueError):
+            _ = Model(x, [x >= 1, x_ >= 1])["x"]
+        with self.assertRaises(ValueError):
+            _ = Model(x, [x >= 1, xv >= 1])["x"]
+        self.assertTrue(all(Model(xv.prod(), [xv >= 1])["x"] == xv))
+        with self.assertRaises(ValueError):
+            _ = Model(xv.prod(), [xv >= 1, xv_ >= 1])["x"]
+        with self.assertRaises(ValueError):
+            _ = Model(xv.prod(), [xv >= 1, x_ >= 1])["x"]
+
     def test_link_conflict(self):
         "Check that substitution conflicts are flagged during linking."
         x_fx1 = Variable("x", 1, models=["fixed1"])
@@ -143,7 +160,7 @@ class TestSignomialInequality(unittest.TestCase):
         with SignomialsEnabled():
             sc = (x + y >= x*y)
         # make sure that the error type doesn't change on our users
-        with self.assertRaises(TypeError):
+        with self.assertRaises(InvalidGPConstraint):
             _ = sc.as_posyslt1()
 
 
diff --git a/gpkit/tests/t_examples.py b/gpkit/tests/t_examples.py
index 55d4f37..278d309 100644
--- a/gpkit/tests/t_examples.py
+++ b/gpkit/tests/t_examples.py
@@ -82,6 +82,9 @@ class TestExamples(unittest.TestCase):
             sol_rat = sol["sensitivities"]["constants"][key]/consenscheck[key]
             self.assertTrue(abs(1-sol_rat) < 1e-2)
 
+    def test_unbounded(self, example):
+        pass
+
 
 FILE_DIR = os.path.dirname(os.path.realpath(__file__))
 EXAMPLE_DIR = os.path.abspath(FILE_DIR + '../../../docs/source/examples')
diff --git a/gpkit/tests/t_model.py b/gpkit/tests/t_model.py
index 3c5b51f..5e11bd4 100644
--- a/gpkit/tests/t_model.py
+++ b/gpkit/tests/t_model.py
@@ -5,6 +5,7 @@ from gpkit import (Model, Monomial, settings, VectorVariable, Variable,
                    SignomialsEnabled, ArrayVariable, SignomialEquality)
 from gpkit.small_classes import CootMatrix
 from gpkit.feasibility import feasibility_model
+from gpkit.exceptions import InvalidGPConstraint
 
 NDIGS = {"cvxopt": 5, "mosek": 7, "mosek_cli": 5}
 # name: decimal places of accuracy
@@ -34,6 +35,7 @@ class TestGP(unittest.TestCase):
                      constraints=[x*y >= 1])
         sol = prob.solve(solver=self.solver, verbosity=0)
         self.assertEqual(type(prob.latex()), str)
+        # pylint: disable=protected-access
         self.assertEqual(type(prob._repr_latex_()), str)
         self.assertAlmostEqual(sol("x"), math.sqrt(2.), self.ndig)
         self.assertAlmostEqual(sol("y"), 1/math.sqrt(2.), self.ndig)
@@ -259,9 +261,7 @@ class TestSP(unittest.TestCase):
         y = Variable('y')
         with SignomialsEnabled():
             m = Model(x, [x >= 1-y, y <= 0.1])
-        with self.assertRaises(ValueError):
-            # solve should catch the TypeError raised by an SP constraints
-            # and raise its own ValueError instead
+        with self.assertRaises(InvalidGPConstraint):
             m.solve(verbosity=0)
         sol = m.localsolve(self.solver, verbosity=0)
         self.assertAlmostEqual(sol["variables"]["x"], 0.9, self.ndig)
@@ -396,7 +396,7 @@ class TestSP(unittest.TestCase):
                 _ = m.localsolve(verbosity=0)
 
     def test_partial_sub_signomial(self):
-        """Test SP partial x0 initialization"""
+        "Test SP partial x0 initialization"
         x = Variable('x')
         y = Variable('y')
         with SignomialsEnabled():
@@ -405,6 +405,20 @@ class TestSP(unittest.TestCase):
         first_gp_constr_posy = m.program.gps[0].constraints[0].as_posyslt1()[0]
         self.assertEqual(first_gp_constr_posy.exp[x.key], -1./3)
 
+    def test_unbounded_debugging(self):
+        "Test nearly-dual-feasible problems"
+        from gpkit.constraints.bounded import BoundedConstraintSet
+        x = Variable("x")
+        y = Variable("y")
+        m = Model(x*y, [x*y**1.000001 >= 100])
+        with self.assertRaises((RuntimeWarning, ValueError)):
+            m.solve(self.solver, verbosity=0)
+        m = Model(x*y, BoundedConstraintSet(m, verbosity=0))
+        sol = m.solve(self.solver, verbosity=0)
+        bounds = sol["boundedness"]
+        self.assertEqual(bounds["sensitive to upper bound"], [y.key])
+        self.assertEqual(bounds["sensitive to lower bound"], [x.key])
+
 
 class TestModelSolverSpecific(unittest.TestCase):
     """test cases run only for specific solvers"""
diff --git a/gpkit/tests/t_solution_array.py b/gpkit/tests/t_solution_array.py
index 7994110..635a809 100644
--- a/gpkit/tests/t_solution_array.py
+++ b/gpkit/tests/t_solution_array.py
@@ -79,6 +79,21 @@ class TestSolutionArray(unittest.TestCase):
         tab = sol.table()
         self.assertTrue(isinstance(tab, str))
 
+    def test_units_sub(self):
+        # issue 809
+        if not gpkit.units:
+            return
+        T = Variable("T", "N", "thrust")
+        Tmin = Variable("T_{min}", "N", "minimum thrust")
+        m = Model(T, [T >= Tmin])
+        tminsub = 1000 * gpkit.units.lbf
+        m.substitutions.update({Tmin: tminsub})
+        sol = m.solve(verbosity=0)
+        self.assertEqual(sol(Tmin), tminsub)
+        self.assertFalse(
+            "1000N" in
+            sol.table().replace(" ", "").replace("[", "").replace("]", ""))
+
 
 class TestResultsTable(unittest.TestCase):
     """TestCase for results_table()"""
diff --git a/gpkit/tests/t_sub.py b/gpkit/tests/t_sub.py
index 74cfeef..3683a5c 100644
--- a/gpkit/tests/t_sub.py
+++ b/gpkit/tests/t_sub.py
@@ -138,6 +138,18 @@ class TestNomialSubs(unittest.TestCase):
 class TestGPSubs(unittest.TestCase):
     """Test substitution for Model and GP objects"""
 
+    def test_persistence(self):
+        x = gpkit.Variable("x")
+        y = gpkit.Variable("y")
+        ymax = gpkit.Variable("y_{max}", 0.1)
+
+        with gpkit.SignomialsEnabled():
+            m = gpkit.Model(x, [x >= 1-y, y <= ymax])
+            m.substitutions[ymax] = 0.2
+            self.assertAlmostEqual(m.localsolve(verbosity=0)["cost"], 0.8, 3)
+            m = gpkit.Model(x, [x >= 1-y, y <= ymax])
+            self.assertAlmostEqual(m.localsolve(verbosity=0)["cost"], 0.9, 3)
+
     def test_united_sub_sweep(self):
         A = Variable("A", "USD")
         h = Variable("h", "USD/count")