updating substitutions with a linked vectorvar

[pgkirsch]

Was a bit of a chore making a MWE for this one but think I got it 😅:

from gpkit import Variable, Model, ConstraintSet, Vectorize                      

class Simple(Model):                                                             
    def setup(self):                                                             
        self.x = x = Variable("x")                                               
        y = Variable("y", 1)                                                     
        z = Variable("z", 2)                                                     
        constraints = [                                                          
            x >= y + z,                                                          
        ]                                                                        
        return constraints                                                       

class Cake(Model):                                                               
    def setup(self):                                                             
        with Vectorize(1):                                                       
            s = Simple()                                                         
        c = ConstraintSet([s])                                                   
        self.cost = sum(s.x)                                                     
        return c                                                                 

m = Cake()                                                                       
m.substitutions.update({                                                         
    "y": ("sweep", [1, 2, 3]),                                                   
    "z": ("sweep", lambda v: v("y")**2),                                         
})                                                                               
sol = m.solve()                                                                  
print(sol.table()) 

Yields:

---> 27 print(sol.table())

~/Documents/Optimization/gpkit/gpkit/solution_array.py in table(self, showvars, tables, **kwargs)
    659                     showvars = self._parse_showvars(showvars)
    660                     data = {k: data[k] for k in showvars if k in data}
--> 661                 strs += var_table(data, self.table_titles[table], **kwargs)
    662         if kwargs.get("latex", None):
    663             preamble = "\n".join(("% \\documentclass[12pt]{article}",

~/Documents/Optimization/gpkit/gpkit/solution_array.py in var_table(data, title, printunits, latex, rawlines, varfmt, valfmt, vecfmt, minval, sortbyvals, hidebelowminval, included_models, excluded_models, sortbymodel, maxcolumns, skipifempty, **_)
    729         print(k)
    730         print(v)
--> 731         if np.isnan(v).all() or np.nanmax(np.abs(v)) <= minval:
    732             continue  # no values below minval
    733         if minval and hidebelowminval and getattr(v, "shape", None):

TypeError: ufunc 'isnan' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe'' 

Digging in, the derived sweep variable is an array of arrays, which I believe is the issue:

#print(k)
Cake.Simple.z[:]
# print(v)
[array([1.])]

I believe this is a bug, but maybe there is a more appropriate way of declaring the linked sweep I should be using?

[bqpd]

ohhh this is interesting. Several things are happening here that aren't bugs per se:

1) you don't need to specify "sweep" for a linked variable: linked take precedence. Should this raise an error? (it currently does not) 2) "z" specifies the whole vector, but because there's only one element GPkit sets just that. Should it set it as an array? (If not, the proper syntax would be lambda v: v("y")[0]**2)

This also uncovered bugs in linked, keydict, and sweeps, addressed in #1514

[bqpd]

The code as written above but with lambda v: v("y")[0]**2 now runs so I'll be closing this, but @pgkirsch curious if you have any thoughts on the questions above!

[pgkirsch]

@bqpd sorry for the slow response on this.

1. Before updating, when I replace the line with

   "z": lambda v: v("y")[0]**2,

   I get the following error:

   
   TypeError                                 Traceback (most recent call last)
   ~/Documents/Optimization/HOPS/hops/t.py in <module>
    22 m.substitutions.update({                                                         
    23     "y": ("sweep", [1, 2, 3]),
   ---> 24     "z": lambda v: v("y")[0]**2,
    25 })                                                                               
    26 sol = m.solve()

~/Documents/Optimization/gpkit/gpkit/keydict.py in update(self, *args, *kwargs) 155 else: 156 for k, v in dict(args, **kwargs).items(): --> 157 self[k] = v 158 159 def call(self, key): # if uniting is ever a speed hit, cache it

~/Documents/Optimization/gpkit/gpkit/keydict.py in setitem(self, key, value) 209 self.owned.add(key) 210 self._copyonwrite(key) --> 211 super().getitem(key)[idx] = value 212 return # succefully set a single index! 213 if key.shape: # now if we're setting an array...

TypeError: float() argument must be a string or a number, not 'function'

Seems to work after updating to latest master though!

2. My MWE was just for a trivial (one-element) vectorization, but for my real case that number is larger than one, i.e. I want every element of `z` to be derived from its counterpart element in `y`. When I take away the `[0]` index (with latest master) I get the following error (still with `Vectorization(1)`):
```python
Sweeping over 3 solves.
Warning: skipped auto-differentiation of linked variable Cake.Simple.z[0] because AttributeError("Neither Quantity object nor its magnitude ([ad(1)]) has attribute 'x'") was raised. Set `gpkit.settings["ad_errors_raise"] = True` to raise such Exceptions directly.

---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
~/Documents/Optimization/HOPS/hops/t.py in <module>
     24     "z": lambda v: v("y")**2,
     25 })                                                                               
---> 26 sol = m.solve()
     27 print(sol.table())

~/Documents/Optimization/gpkit/gpkit/constraints/prog_factories.py in solvefn(self, solver, verbosity, skipsweepfailures, **solveargs)
    101         if sweep:
    102             run_sweep(genfunction, self, solution, skipsweepfailures,
--> 103                       constants, sweep, linked, solver, verbosity, **solveargs)
    104         else:
    105             self.program, progsolve = genfunction(self)

~/Documents/Optimization/gpkit/gpkit/constraints/prog_factories.py in run_sweep(genfunction, self, solution, skipsweepfailures, constants, sweep, linked, solver, verbosity, **solveargs)
    140                           for (var, sweep_vect) in sweep_vects.items()})
    141         if linked:
--> 142             evaluate_linked(constants, linked)
    143         program, solvefn = genfunction(self, constants)
    144         self.program.append(program)  # NOTE: SIDE EFFECTS

~/Documents/Optimization/gpkit/gpkit/constraints/prog_factories.py in evaluate_linked(constants, linked)
     13     "Evaluates the values and gradients of linked variables."
     14     kdc = KeyDict({k: adnumber(maybe_flatten(v), k)
---> 15                    for k, v in constants.items()})
     16     kdc_plain = None
     17     array_calulated = {}

~/Documents/Optimization/gpkit/gpkit/constraints/prog_factories.py in <dictcomp>(.0)
     13     "Evaluates the values and gradients of linked variables."
     14     kdc = KeyDict({k: adnumber(maybe_flatten(v), k)
---> 15                    for k, v in constants.items()})
     16     kdc_plain = None
     17     array_calulated = {}

~/python/anaconda/envs/py37/lib/python3.7/site-packages/ad/__init__.py in adnumber(x, tag)
   1029     raise NotImplementedError(
   1030         'Automatic differentiation not yet supported for {0:} objects'.format(
-> 1031         type(x))
   1032         )
   1033 

NotImplementedError: Automatic differentiation not yet supported for <class 'pint.quantity.build_quantity_class.<locals>.Quantity'> objects

When using Vectorization(2) it looks like the other error above:

~/python/anaconda/envs/py37/lib/python3.7/site-packages/numpy/core/numeric.py in full(shape, fill_value, dtype, order)
    334         dtype = array(fill_value).dtype
    335     a = empty(shape, dtype, order)
--> 336     multiarray.copyto(a, fill_value, casting='unsafe')
    337     return a
    338 

TypeError: float() argument must be a string or a number, not 'function'

[bqpd]

Yeah, unfortunately linked subs for a vector variable are presently a bit of a pain.

Right now there are three ways to do them: 1) provide a vector of functions, each of which returns a single element 2) do that while also setting m["z"].veckey.original_fn to your function which returns an array - this triggers some optimizations 3) the easiest, just because it's what was happening in other codes and so it does both of the above automatically, is to specify that vector function as a value during variable creation

I should rework this system to allow just directly substituting a function for linked vectors, but for now you might find one of the options above useful

[pgkirsch]

Sorry..I'm confused. My interpretation of provide a vector of functions, each of which returns a single element is:

m.substitutions.update({                                                         
    "y": ("sweep", [1, 2, 3]),                                                   
    "z": ("sweep", [lambda v: v("y")[i]**2 for i in range(3)]),                  
}) 

but that yields:

ValueError: substitution {Cake.Simple.z[0]: <function <listcomp>.<lambda> at 0x7f7fb28a9b00>} has invalid value type <class 'function'>.

I unfortunately can't do 3 (easily) in this use case because we want to compare solutions of the same model, one with these linked subs, one without, so I want to avoid changing the variable declaration, if possible.

[bqpd]

this works!

from gpkit import Variable, Model, ConstraintSet, Vectorize

class Simple(Model):
    def setup(self):
        self.x = x = Variable("x")
        y = Variable("y", 1)
        z = Variable("z", 2)
        constraints = [
            x >= y + z,
        ]
        return constraints

class Cake(Model):
    def setup(self):
        with Vectorize(2):
            s = Simple()
        c = ConstraintSet([s])
        self.cost = sum(s.x)
        return c

m = Cake()
m.substitutions.update({
    "y": ("sweep", [[1, 1], [2, 2]]),
    "z": [lambda v: v("y")[0]**2, lambda v: v("y")[1]**2],
})
sol = m.solve()
print(sol.table())

But yeah with Vectorize(1) it doesn't, due to the special treatment mentioned above where one-element vectors are treated as scalars when setting substitutions. And it doesn't work with that i closure because lambdas don't work with closures like that.

[bqpd]

Closures work if you don't use lambdas:

from gpkit import Variable, Model, ConstraintSet, Vectorize

class Simple(Model):
    def setup(self):
        self.x = x = Variable("x")
        y = Variable("y", 1)
        z = Variable("z", 2)
        constraints = [
            x >= y + z,
        ]
        return constraints

class Cake(Model):
    def setup(self):
        with Vectorize(2):
            s = Simple()
        c = ConstraintSet([s])
        self.cost = sum(s.x)
        return c

m = Cake()

def zlink_factory(i):
    def zlinki(v):
        return v("y")[i]**2
    return zlinki

m.substitutions.update({
    "y": ("sweep", [[1, 1], [2, 2]]),
    "z": [zlink_factory(i) for i in range(len(m["z"]))],
})
sol = m.solve()
print(sol.table())

[pgkirsch]

Ah nice! Although it seems to interpret that as a 4-point sweep, rather than a 2-point sweep of a 2-element variable. Is that intentional?

Sweeping over 4 solves.
Sweeping took 0.668 seconds.

Optimal Cost
------------
 [ 4         8         8         12        ]

Swept Variables
---------------
y : [ 1         2         1         2
      1         1         2         2         ]

Sorry if this feels like a pointless thread. If I would be best waiting for a re-work of Vectorized subs, I can do that.

[bqpd]

Ah, yeah, that's just how vector sweeps work.

You should be able to do it before that re-work by following the above!

[bqpd]

Going back over this, I note that there's a sweep of a linked_fn - I don't think we want to support that. The rest, however, will be supported in a coming PR.

[pgkirsch]

Ok fair enough.