Theta no longer vector var, no substitutions, pictures from tests, ge…

…neral cleanup

FeasibilityDemo.py

@@ -1,8 +1,5 @@
-from __future__ import print_function
-from __future__ import absolute_import
 from gpkit import Variable, Model
 import numpy as np
-
 from robust.robust import RobustModel
 from plot_feasibilities import plot_feasibilities
 
@@ -62,47 +59,41 @@
 sol = m.solve()
 
 
-def plot_feasibility_simple_Wing(type_of_uncertainty_set, x, y, str1, val1, str2, val2):
-    plot_feasibilities(x, y, m)
-    x.key.descr[str1] = val1
-    y.key.descr[str2] = val2
+def plot_feasibility_simple_Wing(type_of_uncertainty_set, axisvars, str1, val1, str2, val2):
+    x_name = axisvars[0].key.descr["name"].replace("\\", "")
+    y_name = axisvars[1].key.descr["name"].replace("\\", "")
+    orig_plt = plot_feasibilities(axisvars, m)
+    orig_plt.savefig("test_figs/%s_orig_%s_%s.png" % (type_of_uncertainty_set, x_name, y_name))
+    axisvars[0].key.descr[str1] = val1
+    axisvars[1].key.descr[str2] = val2
     RM = RobustModel(m, type_of_uncertainty_set, linearizationTolerance=1e-4)
     RMsol = RM.robustsolve(verbosity=0, minNumOfLinearSections=20, maxNumOfLinearSections=40)
     print("nominal: ", {k: v for k, v in list(sol["freevariables"].items())
                         if k in m.varkeys and k.key.fix is True})
     print("robust: ", {k: v for k, v in list(RMsol["freevariables"].items())
                        if k in m.varkeys and k.key.fix is True})
     print('cost', RMsol['cost'])
-    plot_feasibilities(x, y, m, RM, numberofsweeps=120, skipfailures=True)
-    del x.key.descr[str1]
-    del y.key.descr[str2]
+    robust_plt = plot_feasibilities(axisvars, m, RM, numberofsweeps=120, skipfailures=True)
+    robust_plt.savefig("test_figs/%s_robust_%s_%s.png" % (type_of_uncertainty_set, x_name, y_name))
+    del axisvars[0].key.descr[str1]
+    del axisvars[1].key.descr[str2]
 
 
 if __name__ == "__main__":
-    # plot_feasibility_simple_Wing('elliptical', W_0, W_W_coeff2, 'r', 1.6, 'r', 1.66, True)
-    # plot_feasibility_simple_Wing('elliptical', W_0, W_W_coeff2, 'r', 1.6, 'r', 1.66, False)
-    # plot_feasibility_simple_Wing('box', W_0, W_W_coeff2, 'pr', 60, 'pr', 66, True)
-    # plot_feasibility_simple_Wing('box', W_0, W_W_coeff2, 'pr', 60, 'pr', 66, False)
-
-    # plot_feasibility_simple_Wing('elliptical', W_W_coeff1, W_W_coeff2, 'r', 1.6, 'r', 1.66, True)
-    # plot_feasibility_simple_Wing('elliptical', W_W_coeff1, W_W_coeff2, 'r', 1.6, 'r', 1.66, False)
-    # plot_feasibility_simple_Wing('box', W_W_coeff1, W_W_coeff2, 'pr', 60, 'pr', 66, True)
-    # plot_feasibility_simple_Wing('box', W_W_coeff1, W_W_coeff2, 'pr', 60, 'pr', 66, False)
-
-    plot_feasibility_simple_Wing('elliptical', C_Lmax, V_min, 'r', 1.25, 'r', 1.2)
-    plot_feasibility_simple_Wing('box', C_Lmax, V_min, 'pr', 25, 'pr', 20)
-
-    # plot_feasibility_simple_Wing('elliptical', dummy, e, 'r', 1.12, 'r', 1.31, True)
-    # plot_feasibility_simple_Wing('elliptical', dummy, e, 'r', 1.12, 'r', 1.31, False)
-    # plot_feasibility_simple_Wing('box', dummy, e, 'pr', 12, 'pr', 31, True)
-    # plot_feasibility_simple_Wing('box', dummy, e, 'pr', 12, 'pr', 31, False)
-
-    # plot_feasibility_simple_Wing('elliptical', rho, dummy, 'r', 1.1, 'r', 1.12, True)
-    # plot_feasibility_simple_Wing('elliptical', rho, dummy, 'r', 1.1, 'r', 1.12, False)
-    # plot_feasibility_simple_Wing('box', rho, dummy, 'pr', 10, 'pr', 12, True)
-    # plot_feasibility_simple_Wing('box', rho, dummy, 'pr', 10, 'pr', 12, False)
-
-    # plot_feasibility_simple_Wing('elliptical', rho, W_0, 'r', 1.1, 'r', 1.6, True)
-    # plot_feasibility_simple_Wing('elliptical', rho, W_0, 'r', 1.1, 'r', 1.6, False)
-    # plot_feasibility_simple_Wing('box', rho, W_0, 'r', 1.1, 'r', 1.6, True)
-    # plot_feasibility_simple_Wing('box', rho, W_0, 'r', 1.1, 'r', 1.6, False)
+    # plot_feasibility_simple_Wing('elliptical', [W_0, W_W_coeff2], 'r', 1.6, 'r', 1.66)
+    # plot_feasibility_simple_Wing('box', [W_0, W_W_coeff2], 'pr', 60, 'pr', 66)
+
+    # plot_feasibility_simple_Wing('elliptical', [W_W_coeff1, W_W_coeff2], 'r', 1.6, 'r', 1.66)
+    # plot_feasibility_simple_Wing('box', [W_W_coeff1, W_W_coeff2], 'pr', 60, 'pr', 66)
+
+    # plot_feasibility_simple_Wing('elliptical', [C_Lmax, V_min], 'r', 1.25, 'r', 1.2)
+    # plot_feasibility_simple_Wing('box', [C_Lmax, V_min], 'pr', 25, 'pr', 20)
+
+    # plot_feasibility_simple_Wing('elliptical', [dummy, e], 'r', 1.12, 'r', 1.31)
+    # plot_feasibility_simple_Wing('box', [dummy, e], 'pr', 12, 'pr', 31)
+
+    # plot_feasibility_simple_Wing('elliptical', [rho, dummy], 'r', 1.1, 'r', 1.12)
+    # plot_feasibility_simple_Wing('box', [rho, dummy], 'pr', 10, 'pr', 12)
+
+    plot_feasibility_simple_Wing('elliptical', [rho, W_0], 'r', 1.1, 'r', 1.6)
+    plot_feasibility_simple_Wing('box', [rho, W_0], 'r', 1.1, 'r', 1.6)

README.md

@@ -6,10 +6,9 @@ Given a GPkit model, this repository will visualize its solution to allow for be
 ## Usage
 
 ```python
-plot_feasibilities(x, y, m, [rm, iterate, design_feasibility, skipfailures, numberofsweeps])
+plot_feasibilities(axisvariables, m, [rm, iterate, design_feasibility, skipfailures, numberofsweeps])
 ```
-x = GPkit variable to plot on the x axis  
-y = GPkit variable to plot on the y axis  
+axisvariables = GPkit variables to plot (first x, then y)
 rm = Robust Model, if it exists (defaults to None)  
 iterate = Boolean of whether or not to use an iterated plot to a specific resolution (defaults to False)  
 skipfailures = Boolean of whether or not to skip errors during sweeps (defaults to False)  

plot_feasibilities.py

@@ -1,24 +1,21 @@
 from __future__ import division
-from builtins import map, range
 import numpy as np
 from matplotlib import pyplot as plt
 from matplotlib.patches import Polygon
-from gpkit import Model, Variable, GPCOLORS, GPBLU, VectorVariable
+from gpkit import Model, Variable, ConstraintSet, GPCOLORS, VectorVariable
 from gpkit.small_scripts import mag
 
-GPRED = GPCOLORS[1]
+GPBLU, GPRED = GPCOLORS[:2]
 
 
-def plot_feasibilities(x, y, m, rm=None, iterate=False, skipfailures=False,
+def plot_feasibilities(axisvariables, m, rm=None, iterate=False, skipfailures=False,
                        numberofsweeps=150):
     """ Plots feasibility space over two variables given a model.
 
     Arguments
     ---------
-    x : GPkit Variable (not a vector variable)
-        plotted on the x axis
-    y : GPkit Variable (not a vector variable)
-        plotted on the y axis
+    axisvariables : list of GPkit Variables (not a vector variable) 
+        currently only 2 elements, first plotted on x, second plotted on y
     m : GPkit Model
     rm : GPkit Model (defaults to None)
         robust model, if it exists
@@ -32,55 +29,53 @@ def plot_feasibilities(x, y, m, rm=None, iterate=False, skipfailures=False,
 
     Raises
     ------
-    ValueError if either x or y is a vector variable
-
+    ValueError if any of the axisvariables is a vector variable
+    
     """
-    axisvars = np.array([x, y])
+    axisvars = np.array(axisvariables)
     if rm:  # we have a robust model as well as a conventional one
         rmtype = rm.type_of_uncertainty_set
         from robust.robust_gp_tools import RobustGPTools
 
     class FeasCircle(Model):
         "Named model that will be swept to get feasibility boundary points."
 
-        def setup(self, m, sol, angles=None, rob=False):
-            r = 4  # radius of the circle in logspace
+        def setup(self, m, sol, angles=None):
+            r = 4 # radius of the circle in logspace
             n_axes = len(axisvars)
-            th = VectorVariable(n_axes-1, "theta", "-", "polar coordinates")
-            angles = angles or np.linspace(0, 2*np.pi, numberofsweeps)
-            substitutions = {var.key: ("sweep", angles) for var in th}
+            angles = angles or np.linspace(0, 2 * np.pi, numberofsweeps)
+            thetas = np.array([Variable("\\theta_%s" % i, angles, "-") for i in range(n_axes - 1)])
             slacks = VectorVariable(n_axes, "s", "-", "slack variables")
             x = np.array([None]*n_axes, "object")  # filled in the loop below
             for i, var in enumerate(axisvars):
                 if var.key.shape:
                     raise ValueError("Uncertain VectorVariables not supported")
+                x_center = np.log(mag(m.solution(axisvars[i])))
 
-                eta_max_x = 0  # TODO: what are these?
-                if rm:
-                    _, eta_max_x = RobustGPTools.generate_etas(var)
-                x_center = np.log(mag(m.solution(var)))
-
-                def f(c):
+                def f(c, index=i, x_val=x_center):
                     product = 1
-                    for j in range(i):
-                        product *= np.cos(c[th[j]])
-                    if i != n_axes - 1:  # TODO: should this be indented??
-                        product *= np.sin(c[th[i]])
-                    return np.exp(x_center) * np.exp(r * product)
-                x[i] = Variable("x_%s" % i, f, var.unitstr())
-
-            self.cost = (slacks**2).sum() * m.cost**0.1
-            fixedsubs = {k: v for k, v in list(sol["freevariables"].items())
-                         if k in m.varkeys and k.key.fix is True}
-            substitutions.update(fixedsubs)
-            return {"original model": m,
-                    "slack constraints": [slacks >= 1,
-                                          axisvars/slacks <= x,
-                                          axisvars*slacks >= x]}, substitutions
+                    for j in range(index):
+                        product *= np.cos(c[thetas[j]])
+                    if index != n_axes - 1: # TODO: should this be indented??
+                        product *= np.sin(c[thetas[index]])
+                    return np.exp(x_val) * np.exp(r * product)
+
+                x[i] = Variable('x_%s' % i, f, axisvars[i].unitstr())
+
+            constraints = [(slacks >= 1),
+                          (axisvars/slacks <= x),
+                          (axisvars*slacks >= x)]
+
+            cost_ref = Variable('cost_ref', 1, m.cost.unitstr(), "reference cost")
+            self.cost = (slacks**2).sum() * m.cost / cost_ref
+            slack_constr = ConstraintSet(constraints)
+            return ({"original model": m, "slack constraints": slack_constr},
+                {k: v for k, v in list(sol["freevariables"].items())
+                if k in m.varkeys and k.key.fix is True})
 
     def slope(a, b):
         if a[0] == b[0]:
-            return 1e10  # arbitrarily large number. TODO: should be np.inf?
+            return 1e10  # arbitrarily large number. TODO: if np.inf, modify arithmetic in angle to account for possibility
         return (b[1]-a[1])/(b[0]-a[0])
 
     def distance(a, b):
@@ -152,34 +147,33 @@ def iterate_angles(rob=False, spacing_start=150, slope_error=np.pi/4, dist_error
         b = [solved_angles[angle][1] for angle in angles]
         return a, b
 
-    # plot original feasibility set
-    # plot boundary of uncertainty set
+    # plot original feasibility set and boundary of uncertainty set
     orig_a, orig_b, a, b = [None] * 4
     if iterate:
         if rm:
             a, b = iterate_angles(rob=True)
         orig_a, orig_b = iterate_angles()
     else:
         if rm:
-            fc = FeasCircle(m, rm.get_robust_model().solution, rob=True)
+            fc = FeasCircle(m, rm.get_robust_model().solution)
             for axisvar in axisvars:
                 del fc.substitutions[axisvar]
             rmfeas = fc.solve(skipsweepfailures=skipfailures)
-            a, b = list(map(mag, list(map(rmfeas, [x, y]))))
+            a, b = list(map(mag, list(map(rmfeas, axisvars))))
         ofc = FeasCircle(m, m.solution)
         for axisvar in axisvars:
             del ofc.substitutions[axisvar]
         origfeas = ofc.solve(skipsweepfailures=skipfailures)
-        orig_a, orig_b = list(map(mag, list(map(origfeas, [x, y]))))
+        orig_a, orig_b = list(map(mag, list(map(origfeas, axisvars))))
 
     fig, axes = plt.subplots(2)
 
     def plot_uncertainty_set(ax):
-        xo, yo = list(map(mag, list(map(m.solution, [x, y]))))
+        xo, yo = list(map(mag, list(map(m.solution, axisvars))))
         ax.plot(xo, yo, "k.")
         if rm:
-            eta_min_x, eta_max_x = RobustGPTools.generate_etas(x)
-            eta_min_y, eta_max_y = RobustGPTools.generate_etas(y)
+            eta_min_x, eta_max_x = RobustGPTools.generate_etas(axisvars[0])
+            eta_min_y, eta_max_y = RobustGPTools.generate_etas(axisvars[1])
             x_center = np.log(xo)
             y_center = np.log(yo)
             ax.plot(np.exp(x_center), np.exp(y_center), "kx")
@@ -215,10 +209,10 @@ def plot_uncertainty_set(ax):
 
     plot_uncertainty_set(axes[0])
     axes[0].axis("equal")
-    axes[0].set_ylabel(y)
+    axes[0].set_ylabel(axisvars[1])
     plot_uncertainty_set(axes[1])
-    axes[1].set_xlabel(x)
-    axes[1].set_ylabel(y)
+    axes[1].set_xlabel(axisvars[0])
+    axes[1].set_ylabel(axisvars[1])
 
-    fig.suptitle("%s vs %s Feasibility Space" % (x, y))
-    plt.show()
+    fig.suptitle("%s vs %s Feasibility Space" % (axisvars[0], axisvars[1]))
+    return plt

tests.py

@@ -13,13 +13,13 @@ class TestFeasibilityDemo(unittest.TestCase):
     def test_elliptical(self):
         if settings["default_solver"] == "cvxopt":
             return
-        plot_feasibility_simple_Wing('elliptical', C_Lmax, V_min,
+        plot_feasibility_simple_Wing('elliptical', [C_Lmax, V_min],
                                      'r', 1.25, 'r', 1.2)
 
     def test_box(self):
         if settings["default_solver"] == "cvxopt":
             return
-        plot_feasibility_simple_Wing('box', C_Lmax, V_min,
+        plot_feasibility_simple_Wing('box', [C_Lmax, V_min],
                                      'pr', 25, 'pr', 20)