Vector substitutions for vectorized models

[pgkirsch]

(Forgive me, I know this has been discussed before but I can't find a ticket that answers my question)

When a scalar substitution is made in a vectorized model, it seems that substitution is "tiled" for any arbitrary model vectorization.

But when a multi-element substitution is made for a vectorized variable, it only seems to work as long as the parent model is not vectorized.

Here is a MWE of something that works:

import numpy as np                                                               
from gpkit import Variable, Model, ConstraintSet, Vectorize                      
                                                                                 
class Pie(Model):                                                                
    def setup(self):                                                             
        self.x = x = Variable("x")                                               
        z = Variable("z")                                                        
        constraints = [                                                          
            x >= z,                                                              
        ]                                                                        
        substitutions = {'z': 1}                                                 
        return constraints, substitutions                                        
                                                                                 
class Cake(Model):                                                               
    def setup(self):                                                             
        self.y = y = Variable("y")                                               
        with Vectorize(2):                                                       
            s = Pie()                                                            
        constraints = [y >= s.x]                                                 
        constraints += [s]                                                       
        subs = {'x': np.array([[2], [3]])}                                       
        return constraints, subs                                                 
                                                                                 
class Yum1(Model):                                                               
    def setup(self):                                                             
        with Vectorize(1):                                                       
            cake = Cake()                                                        
        y = cake.y                                                               
        self.cost = sum(y)                                                       
        constraints = ConstraintSet([cake])                                      
        return constraints                                                       
                                                                                 
m = Yum1()                                                                       
sol = m.solve()                                                                  
print(sol.table())

Optimal Cost
------------
 3

Free Variables
--------------
  | Yum1.Cake
y : [ 3        ]

Fixed Variables
---------------
  | Yum1.Cake.Pie
x : [ 2         3        ]
z : [ 1         1        ]

Here is a MWE of something that breaks:

# same code above
class Yum2(Model):                                                               
    def setup(self):                                                             
        with Vectorize(2):      # <--- only difference                                                 
            cake = Cake()                                                        
        y = cake.y                                                               
        self.cost = sum(y)                                                       
        constraints = ConstraintSet([cake])                                      
        return constraints                                                       
                                                                                 
m = Yum2()                                                                       
sol = m.solve()                                                                  
print(sol.table())


---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
~/Documents/Optimization/HOPS/hops/t2.py in <module>
     45 
     46 m = Yum2()
---> 47 sol = m.solve()
     48 print(sol.table())

~/Documents/Optimization/gpkit/gpkit/constraints/prog_factories.py in solvefn(self, solver, verbosity, skipsweepfailures, **solveargs)
    114          RuntimeWarning if an error occurs in solving or parsing the solution.
    115          """
--> 116         constants, sweep, linked = parse_subs(self.varkeys, self.substitutions)
    117         solution = SolutionArray()
    118         solution.modelstr = str(self)

~/Documents/Optimization/gpkit/gpkit/nomials/substitution.py in parse_subs(varkeys, substitutions, clean)
     20                     sub = dict.__getitem__(substitutions, var)
     21                     keys = varkeys.keymap[var]
---> 22                     append_sub(sub, keys, constants, sweep, linkedsweep)
     23         else:
     24             for var in substitutions:

~/Documents/Optimization/gpkit/gpkit/nomials/substitution.py in append_sub(sub, keys, constants, sweep, linkedsweep)
     65                 raise ValueError("cannot substitute array of shape %s for"
     66                                  " variable %s of shape %s." %
---> 67                                  (sub.shape, key.veckey, key.shape))
     68         if hasattr(value, "__call__") and not hasattr(value, "key"):
     69             linkedsweep[key] = value

ValueError: cannot substitute array of shape (2, 1) for variable Yum2.Cake.Pie.x[:] of shape (2, 2).

The only difference is that Cake is vectorized in Yum2. The workaround is to use np.tile to make the substitution the correct shape (using the global vectorization) but that's inelegant (in the non-toy version of my problem) and I'd like to move away from it if at all possible.

Is there code that correctly tiles scalars but not vectors or am I thinking about it wrong? Is there something ambiguous about the substitutions I'm trying to make that would prevent this from being possible?

[whoburg]

Just want to say, this is a really well-written issue. Assuming this prompts a code change, let's definitely add these MWE's as a unit test.

[pgkirsch]

Thanks @whoburg!

[pgkirsch]

For the sake of completeness (and for any users experiencing this issue in the future), it's worth including two workarounds in case this doesn't prompt a code change.

Option 1: Direct substitution

Downside: doesn't address post-model-creation substitutions

import numpy as np                                                               
from gpkit import Variable, Model, ConstraintSet, Vectorize                      
                                                                                 
class Pie(Model):                                                                
    def setup(self):                                                             
        # self.x = x = Variable("x")            <-- before                  
        self.x = x = Variable("x", [[2], [3]])      
        z = Variable("z")                                                        
        constraints = [                                                          
            x >= z,                                                              
        ]                                                                        
        substitutions = {'z': 1}                                                 
        return constraints, substitutions                                        
                                                                                 
class Cake(Model):                                                               
    def setup(self):                                                             
        self.y = y = Variable("y")                                               
        with Vectorize(2):                                                       
            s = Pie()                                                            
        constraints = [y >= s.x]                                                 
        constraints += [s]                                                       
        #subs = {'x': np.array([[2], [3]])}     <-- before                                                                                               
        #return constraints, subs               <-- before
        return constraints                                                                                                           
                                                                                 
class Yum2(Model):                                                               
    def setup(self):                                                             
        with Vectorize(2):                                                       
            cake = Cake()                                                        
        y = cake.y                                                               
        self.cost = sum(y)                                                       
        constraints = ConstraintSet([cake])                                      
        return constraints                                                       
                                                                                 
m = Yum2()                                                                       
sol = m.solve()                                                                  
print(sol.table())                                                               

Option 2: numpy.tile

Downside: requires sub-model-level knowledge of its own vectorization. Also, gross.

import numpy as np                                                               
from gpkit import Variable, Model, ConstraintSet, Vectorize                      
                                                                                 
class Pie(Model):                                                                
    def setup(self):                                                             
        self.x = x = Variable("x")                  
        z = Variable("z")                                                        
        constraints = [                                                          
            x >= z,                                                              
        ]                                                                        
        substitutions = {'z': 1}                                                 
        return constraints, substitutions                                        
                                                                                 
class Cake(Model):                                                               
    def setup(self):                                                             
        self.y = y = Variable("y")                                               
        with Vectorize(2):                                                       
            s = Pie()                                                            
        constraints = [y >= s.x]                                                 
        constraints += [s]                                                       
        #subs = {'x': np.array([[2], [3]])}     <-- before
        numberofcakes = Vectorize.vectorization[0]                               
        subs = {'x': np.tile(np.array([[2], [3]]), numberofcakes)}                                                                                                
        return constraints, subs
                                          
                                                                                 
class Yum2(Model):                                                               
    def setup(self):                                                             
        with Vectorize(2):                                                       
            cake = Cake()                                                        
        y = cake.y                                                               
        self.cost = sum(y)                                                       
        constraints = ConstraintSet([cake])                                      
        return constraints                                                       
                                                                                 
m = Yum2()                                                                       
sol = m.solve()                                                                  
print(sol.table())                                                               

[bqpd]

Right now the cleanest fix seems to me to be vectorizing that substitution upon finishing Cake's setup; this would require going through cake.substitutions to find substitutions that match the last dimensions (including any vectorizations within Cake) of their variable but are missing all (and only) those dimensions in the vectorization environment.

[bqpd]

in re: workarounds, option 1 is very fragile; I would recommend using option 2 (with a big #FIXME comment) instead of giving a vector substitution to what would look to any other code reader like a scalar variable.

[bqpd]

@pgkirsch I'm thinking of implementing the fix mentioned two comments above; does that seem like it addresses your issue?

[pgkirsch]

If that seems like the cleanest solution to you -- and it doesn't feel like too much of a hack to appease a maybe-weird use-case -- then that sounds good to me! I guess I was expecting the change to be to Vectorize rather than setup but you obviously know way better than me!

[bqpd]

@pgkirsch looking at this again, I there's an Option 3: you don't have to propagate the number of cakes (as in your Option 2 above), because that information is already encoded in x. By replacing the substitution in Cake with:

        xsub = np.broadcast_to([2, 3], reversed(s.x.shape)).T
        subs = {'x': xsub}

it works for any number of cakes.

This is obviously a bit clunky. The alternative is this being added to model.py:

        CostedConstraintSet.__init__(self, cost, constraints)
        if substitutions:  # check if we have to broadcast any vecsubs
            for key, value in substitutions.items():
                try:
                    vk, idx = self.substitutions.parse_and_index(key)
                    if not idx and vk.shape != getattr(value, "shape", None):
                        substitutions[key] = np.broadcast_to(value,
                                                             reversed(vk.shape)).T
                except KeyError:  # didn't parse, presume fine
                    pass
            self.substitutions.update(substitutions)

replacing CostedConstraintSet.__init__(self, cost, constraints, substitutions).

This code adds more surface area for bugs...but I really like being able to say subs = {'x': [2, 3]} and have it Just Work. Thoughts? Especially given that this is specific to the returned-substitutions-dictionary that you're the only heavy user of?

[bqpd]

I tried running tests with the modification and found it interacted with boundedness tests in a negative way :(

Could be fixed by moving that part of ConstraintSet into a method, but this leans me further towards not implementing this for returned substitutions dictionaries.

[pgkirsch]

Thanks for looking into this more @bqpd! I don't want to break other functionality with this (for now) very uncommon use-case. That xsub line seems like a great option. Given that I almost certainly won't remember it, how would you feel about making a method (e.g. vector_sub(substitution, varkey) or some better name) in small_scripts.py?

Use would be:

subs = {x: vector_sub([2, 3], x)}

Though I'm sure there's probably a more elegant solution that doesn't involve redundant x...

[bqpd]

yeah that's a great idea! called it broadcast_substitutions