Allowing substituted, unitless signomials in exponents.

gpkit/constraints/gp.py

@@ -5,7 +5,7 @@
 from time import time
 from collections import defaultdict
 import numpy as np
-from ..nomials import NomialData
+from ..nomials import NomialData, NomialMap
 from ..small_classes import CootMatrix, SolverLog, Numbers, FixedScalar
 from ..keydict import KeyDict
 from ..small_scripts import mag
@@ -122,7 +122,7 @@ def gen(self):
             self._cs.extend(hmap.values())
         self.vks = self.varlocs
         self.A, self.missingbounds = genA(self.exps, self.varlocs,
-                                          self.meq_idxs)
+                                          self.meq_idxs, self.substitutions)
 
     # pylint: disable=too-many-statements, too-many-locals
     def solve(self, solver=None, verbosity=1, warn_on_check=False,
@@ -394,7 +394,7 @@ def _almost_equal(num1, num2):
                                  " cost %s" % (np.exp(dual_cost), cost))
 
 
-def genA(exps, varlocs, meq_idxs):  # pylint: disable=invalid-name
+def genA(exps, varlocs, meq_idxs, substitutions=None):  # pylint: disable=invalid-name
     """Generates A matrix
 
     Returns
@@ -407,24 +407,50 @@ def genA(exps, varlocs, meq_idxs):  # pylint: disable=invalid-name
     """
     missingbounds = {}
     row, col, data = [], [], []
+    bte = set() # variables bounded through equalities or sig exp substitutions
     for j, var in enumerate(varlocs):
+        # print(var)
         upperbound, lowerbound = False, False
         row.extend(varlocs[var])
         col.extend([j]*len(varlocs[var]))
-        data.extend(exps[i][var] for i in varlocs[var])
-        for i in varlocs[var]:
-            if i not in meq_idxs:
+        exp_arr = []
+        # Adding data to A matrix
+        for k, i in enumerate(varlocs[var]):
+            if isinstance(exps[i][var], NomialMap):
+                varkeydict = KeyDict({key:key for item in exps[i][var].keys() \
+                                      for key in item.keys()})
+                subbed_exp = exps[i][var].sub(substitutions, varkeydict)
+                if list(subbed_exp.keys())[0]:
+                    raise ValueError("Signomial exponent %s has variables that " 
+                                     "have not been fully substituted. Complete "
+                                     "substitutions and try again." % \
+                                     (exps[i])) #TODO: improve error.
+                exp_arr.extend([list(subbed_exp.values())[0]])
+                # Add substituted variables in signomial exponent to bte
+                # to make them bounded.
+                bte = bte.union(varkeydict.keys())
+            else:
+                exp_arr.extend([exps[i][var]])
+        # print(exp_arr)
+        data.extend(exp_arr)
+        # Checking boundedness
+        for k, i in enumerate(varlocs[var]):
+            # Checking variables subbed in exponent
+            if var in bte:
+                lowerbound = True
+                upperbound = True
+                break
+            elif i not in meq_idxs:
                 if upperbound and lowerbound:
                     break
-                elif exps[i][var] > 0:  # pylint:disable=simplifiable-if-statement
+                elif exp_arr[k] > 0:
                     upperbound = True
                 else:
                     lowerbound = True
         if not upperbound:
             missingbounds[(var, "upper")] = ""
         if not lowerbound:
             missingbounds[(var, "lower")] = ""
-
     check_mono_eq_bounds(missingbounds, gen_mono_eq_bounds(exps, meq_idxs))
 
     # space the matrix out for trailing constant terms
@@ -433,6 +459,9 @@ def genA(exps, varlocs, meq_idxs):  # pylint: disable=invalid-name
             row.append(i)
             col.append(0)
             data.append(0)
+    if len(row) != len(col) != len(data):
+        raise ValueError("The A matrix generated does not have the right "
+                         "dimensions.")
     A = CootMatrix(row, col, data)
 
     return A, missingbounds

gpkit/constraints/sgp.py

@@ -194,6 +194,7 @@ def penalty_ccp_solve(self, solver=None, verbosity=1, x0=None, reltol=1e-4,
                                                x0, reltol, iteration_limit,
                                                mutategp, **kwargs)
         self.gps = relaxed_model.gps
+        self.solver_outs = relaxed_model.solver_outs
         return self.result
 
     @property

gpkit/nomials/core.py

@@ -1,9 +1,8 @@
 "The shared non-mathematical backbone of all Nomials"
 from __future__ import unicode_literals, print_function
 from .data import NomialData
-from ..small_classes import Numbers, FixedScalar
-from ..small_scripts import nomial_latex_helper
-
+from ..small_classes import Numbers, FixedScalar, HashVector
+from ..small_scripts import nomial_latex_helper, try_str_without
 
 class Nomial(NomialData):
     "Shared non-mathematical properties of all nomials"
@@ -24,7 +23,9 @@ def str_without(self, excluded=()):
                 for (var, x) in exp.items():
                     if x != 0:
                         varstr = var.str_without(excluded)
-                        if x != 1:
+                        if isinstance(x, (HashVector)):
+                            varstr += "^(%s)" % try_str_without(x, []) #TODO: fix this string...
+                        elif x != 1:
                             varstr += "^%.2g" % x
                         varstrs.append(varstr)
                 varstrs.sort()

gpkit/nomials/map.py

@@ -20,8 +20,9 @@ class NomialMap(HashVector):
     x and y are VarKey objects.
     """
     units = None
-    expmap = None  # used for monomial-mapping postsubstitution; see .mmap()
-    csmap = None   # used for monomial-mapping postsubstitution; see .mmap()
+    expmap = None   # used for monomial-mapping postsubstitution; see .mmap()
+    csmap = None    # used for monomial-mapping postsubstitution; see .mmap()
+    varexps = False # used for signomial exponent operations
 
     def copy(self):
         "Return a copy of this"

gpkit/nomials/math.py

@@ -5,6 +5,7 @@
 from .core import Nomial
 from .array import NomialArray
 from .. import units
+from ..keydict import KeyDict
 from ..constraints import SingleEquationConstraint
 from ..globals import SignomialsEnabled
 from ..small_classes import Numbers
@@ -63,8 +64,22 @@ def __init__(self, hmap=None, cs=1, require_positive=True):  # pylint: disable=t
             self.__class__ = Monomial
         else:
             self.__class__ = Posynomial
+        self.check_exps()
         self.ast = ()
 
+    def check_exps(self):
+        """Checking units and adding varkeys of signomial exponents."""
+        for exp in self.exps:
+            for _, v in exp.items():
+                if isinstance(v, HashVector):
+                    if v.units:
+                        raise ValueError("Exponent %s is united. Please "
+                                         "normalize or remove units." % v)
+                    else:
+                        self.hmap.varexps = True
+                        for key, _ in v.items():
+                            self.vks.update([k for k in key.keys()])
+
     def diff(self, var):
         """Derivative of this with respect to a Variable
 
@@ -119,7 +134,11 @@ def mono_approximation(self, x0):
         Monomial (unless self(x0) < 0, in which case a Signomial is returned)
         """
         x0, _, _ = parse_subs(self.varkeys, x0)  # use only varkey keys
-        psub = self.hmap.sub(x0, self.varkeys, parsedsubs=True)
+        if self.hmap.varexps:
+            raise ValueError("Varexps only works for GPs for now... "
+                             "stick to GPs!")
+        else:
+            psub = self.hmap.sub(x0, self.varkeys, parsedsubs=True)
         if EMPTY_HV not in psub or len(psub) > 1:
             raise ValueError("Variables %s remained after substituting x0=%s"
                              " into %s" % (psub, x0, self))
@@ -329,10 +348,27 @@ def __rtruediv__(self, other, rev=True):
         return self.__rdiv__(other)
 
     def __pow__(self, expo):
-        if isinstance(expo, Numbers):
+        if isinstance(expo, Numbers+(Signomial,)):
             (exp, c), = self.hmap.items()
-            exp = exp*expo if expo else EMPTY_HV
-            hmap = NomialMap({exp: c**expo})
+            if isinstance(expo, Signomial):
+                if expo.units:
+                    raise ValueError("Exponents cannot be united. Please "
+                                     "remove units from Signomial %s." %
+                                     expo.str_without())
+                exp = HashVector({k: expo.hmap*v for k, v in exp.items()})
+            else:
+                exp = exp*expo if expo else EMPTY_HV
+            if c != 1 and isinstance(expo, Signomial):
+                raise ValueError("Float %s raised to Signomial %s is "
+                                 "not currently supported. Please replace %s "
+                                 "with a substituted "
+                                 "Variable." % (str(c), expo.str_without(),
+                                                str(c)))
+            elif c != 1:
+                newc = c**expo
+            else:
+                newc = c
+            hmap = NomialMap({exp: newc})
             if expo and self.hmap.units:
                 hmap.units = self.hmap.units**expo
             else:
@@ -438,9 +474,9 @@ def __init__(self, left, oper, right):
         if self.unsubbed:
             for exp in self.unsubbed[0].hmap:
                 for key, e in exp.items():
-                    if e > 0:
+                    if isinstance(e, Numbers) and e > 0:
                         self.bounded.add((key, "upper"))
-                    if e < 0:
+                    if isinstance(e, Numbers) and e < 0:
                         self.bounded.add((key, "lower"))
         for key in self.substitutions:
             for bound in ("upper", "lower"):
@@ -537,10 +573,19 @@ def sens_from_dual(self, la, nu, result):  # pylint: disable=unused-argument
                 for idx, percentage in self.const_mmap.items():
                     nu_[idx] += percentage * la*scale
             nu = nu_
-        return {var: sum([presub.exps[i][var]*nu[i]
-                          for i in presub.varlocs[var]])
-                for var in self.varkeys}  # Constant sensitivities
-
+        sens_dict = KeyDict({v: 0*v.units for v in self.varkeys if v.units})
+        sens_dict.update({v: 0 for v in self.varkeys if not v.units})
+        for var in self.varkeys:
+            for i in presub.varlocs[var]:
+                if not isinstance(presub.exps[i][var], HashVector):
+                    sens_dict[var] += presub.exps[i][var]*nu[i]
+                else:
+                    exps = presub.exps[i][var]
+                    varkeydict = KeyDict({key:key for item in exps.keys() \
+                                      for key in item.keys()})
+                    subbed_exp = exps.sub(result, varkeydict).values()[0]
+                    sens_dict[var] += subbed_exp*nu[i]
+        return sens_dict
     def as_gpconstr(self, x0):  # pylint: disable=unused-argument
         "The GP version of a Posynomial constraint is itself"
         return self
@@ -641,9 +686,9 @@ def __init__(self, left, oper, right):
         if self.unsubbed:
             for exp, c in self.unsubbed[0].hmap.items():
                 for key, e in exp.items():
-                    if e*c > 0:
+                    if isinstance(e, Numbers) and e*c > 0:
                         self.bounded.add((key, "upper"))
-                    if e*c < 0:
+                    if isinstance(e, Numbers) and e*c < 0:
                         self.bounded.add((key, "lower"))
         for key in self.substitutions:
             for bound in ("upper", "lower"):

gpkit/repr_conventions.py

@@ -123,10 +123,11 @@ def parse_ast(self, excluded=("units")):
             aststr = "-%s" % parenthesize(strify(values, excluded), mult=False)
         elif oper == "pow":
             left = parenthesize(strify(values[0], excluded))
-            x = values[1]
+            x = parenthesize(strify(values[1], excluded))
             if left == "1":
                 aststr = "1"
-            elif UNICODE_EXPONENTS and int(x) == x and x >= 2 and x <= 9:
+            elif isinstance(x, Numbers) and UNICODE_EXPONENTS and \
+                int(x) == x and x >= 2 and x <= 9:
                 if int(x) in (2, 3):
                     aststr = "%s%s" % (left, unichr(176+x))
                 elif int(x) in (4, 5, 6, 7, 8, 9):

gpkit/small_classes.py

@@ -192,17 +192,17 @@ def __neg__(self):
 
     def __pow__(self, other):
         "Accepts scalars. Return Hashvector with each value put to a power."
-        if isinstance(other, Numbers):
-            return self.__class__({key: val**other
-                                   for (key, val) in self.items()})
+        return self.__class__({key: val**other
+                               for (key, val) in self.items()})
         return NotImplemented
 
     def __mul__(self, other):
         """Accepts scalars and dicts. Returns with each value multiplied.
 
         If the other object inherits from dict, multiplication is element-wise
         and their key's intersection will form the new keys."""
-        if isinstance(other, Numbers):
+        from gpkit.nomials.math import Signomial
+        if isinstance(other, Numbers + (Signomial,)):
             return self.__class__({key: val*other
                                    for (key, val) in self.items()})
         elif isinstance(other, dict):

gpkit/tests/t_nomials.py

@@ -2,8 +2,9 @@
 import math
 import sys
 import unittest
+import numpy as np
 from gpkit import Variable, Monomial, Posynomial, Signomial, SignomialsEnabled
-from gpkit import VectorVariable, NomialArray
+from gpkit import VectorVariable, NomialArray, Model, units
 from gpkit.nomials import NomialMap
 from gpkit.small_classes import HashVector
 from gpkit.exceptions import InvalidPosynomial
@@ -273,6 +274,34 @@ def test_eq_ne(self):
             self.assertEqual(Signomial(-3), -3)
             self.assertNotEqual(Signomial(-3), 3)
 
+    def test_subbed_sig_exp(self):
+        a = Variable('a')
+        b = Variable('b')
+        c = Variable('c')
+        x = Variable('x')
+        y = Variable('y', 'm')
+        z = Variable('z', 'm^2')
+        with SignomialsEnabled():
+            with self.assertRaises(ValueError):
+                #pylint: disable=unused-variable
+                mony = (2*a)**c # float**signomial check
+                mony = a**y # united signomial check
+                mony = a**(2*b + z*y) # dimension check
+            # mony = a**(2*b + z/y**2)
+            # subs = {a: 3, b: 2, z: 4, y: 0.3}
+            # self.assertEqual()
+            constraints = [b*x**(y*z**(-0.5)) >= 1]
+            constraints_subbed = [3*x**(1*2**(-0.5)) >= 1]
+            m = Model(x, constraints)
+            m_subbed = Model(x, constraints_subbed)
+            m.substitutions.update({y: 1*units('m')})
+            self.assertRaises(ValueError, m.solve, verbosity=0) # substitutions check
+            m.substitutions.update({b:3, z: 2*units('m^2')})
+            sol = m.solve(verbosity=0)
+            self.assertEqual(sol(x),
+                             m_subbed.solve(verbosity=0)(x))
+            self.assertAlmostEqual(sol['sensitivities']['constants'][b],
+                                   -np.sqrt(2), places=5)
 
 class TestPosynomial(unittest.TestCase):
     """TestCase for the Posynomial class"""