Diffing the solution of an m-by-n vectorized model with a 1-by-m-by-n vectorized model doesn't work as expected

[pgkirsch]

Not unrelated to #1465 (which I promise I will reply to soon).

As part of converting a model from m x n vectorization to l x m x n vectorization, I want to baby-step via a 1 x m x n solution and diff.

I would expect the diff to be of length m x n for each var, instead I seem to be getting m x m x n length diffs.

MWE (adapted from @1ozturkbe's example, thank you):

In [1]: from gpkit import VectorVariable, Model, Vectorize 
   ...: import numpy as np                                                                                                              

In [2]: with Vectorize(2):     
   ...:     v = VectorVariable(5, 'v') 
   ...:     u = VectorVariable(5, 'u ') 
   ...: m = Model(np.prod(v), [v >= u]) 
   ...: m.substitutions.update({u: np.random.rand(5, 2)}) 
   ...: sol = m.solve() 
   ...: print(sol.table())                                                                                                              
Using solver 'mosek_cli'
Solving for 10 variables.
Solving took 0.118 seconds.

Cost
----
 0.001461

Free Variables
--------------
v : [ 0.408     0.745     0.552     0.65      0.235
      0.467     0.799     0.613     0.364     0.686     ]

Constants
---------
u  : [ 0.408     0.745     0.552     0.65      0.235
       0.467     0.799     0.613     0.364     0.686     ]

...

In [3]: with Vectorize(1): 
   ...:     with Vectorize(2): 
   ...:         v = VectorVariable(5, 'v') 
   ...:         u = VectorVariable(5, 'u ') 
   ...: m = Model(np.prod(v), [v >= u]) 
   ...: m.substitutions.update({u: np.random.rand(5, 2, 1)}) 
   ...: sol2 = m.solve() 
   ...: print(sol2.table())                                                                                                             
Using solver 'mosek_cli'
Solving for 10 variables.
Solving took 0.133 seconds.

Cost
----
 9.472e-05

Free Variables
--------------
v : [ 0.509     0.271     0.336     0.327     0.736
      0.597     0.889     0.0734    0.308     0.71      ]

Constants
---------
u  : [ 0.509     0.271     0.336     0.327     0.736
       0.597     0.889     0.0734    0.308     0.71      ]

...

In [4]: print(sol.diff(sol2))                                                                                                           
Solution difference
(positive means the argument is bigger)
-------------------
u  : [  +24.9%    -31.7%     -7.9%    -21.7%   +116.7%
        -42.0%    -66.1%    -55.8%    -25.6%    -60.6%
        -17.5%    -54.9%    -39.1%    -48.2%    +43.3%
        -30.1%    -59.1%    -46.7%    -10.2%    -52.4%
        +80.5%     -1.3%    +33.2%    +13.3%   +213.3%
        +27.8%    -25.3%     -2.6%    +64.0%    -13.1%
       +118.0%    +19.2%    +60.9%    +36.8%   +278.4%
        -84.3%    -90.8%    -88.0%    -79.8%    -89.3%
        -24.5%    -58.7%    -44.3%    -52.6%    +31.0%
        +52.1%    -11.1%    +15.9%    +95.2%     +3.4%   ]
 v : [  +24.9%    -31.7%     -7.9%    -21.7%   +116.7%
        -42.0%    -66.1%    -55.8%    -25.6%    -60.6%
        -17.5%    -54.9%    -39.1%    -48.2%    +43.3%
        -30.1%    -59.1%    -46.7%    -10.2%    -52.4%
        +80.5%     -1.3%    +33.2%    +13.3%   +213.3%
        +27.8%    -25.3%     -2.6%    +64.0%    -13.1%
       +118.0%    +19.2%    +60.9%    +36.8%   +278.4%
        -84.3%    -90.8%    -88.0%    -79.8%    -89.3%
        -24.5%    -58.7%    -44.3%    -52.6%    +31.0%
        +52.1%    -11.1%    +15.9%    +95.2%     +3.4%   ]

Solution sensitivity delta
--------------------------
The largest sensitivity delta is +0

[bqpd]

It doesn't even diff for me; looking into it.

[bqpd]

fixed in https://github.com/convexengineering/gpkit/pull/1468, where

from gpkit import VectorVariable, Model, Vectorize
import numpy as np

with Vectorize(2):
    v = VectorVariable(5, "v", "m")
    u = VectorVariable(5, "u", "ft")
m = Model(np.prod(v), [v >= u])
m.substitutions.update({u: np.random.rand(5, 2)})
sol = m.solve(verbosity=0)

with Vectorize(3):
    with Vectorize(2):
        v = VectorVariable(5, "v", "m")
        u = VectorVariable(5, "u", "ft")
m = Model(np.prod(v), [v >= u])
m.substitutions.update({u: np.random.rand(5, 2, 3)})
sol2 = m.solve(verbosity=0)

print("First u")
print(sol.table(tables=["constants"]))
print("Second u")
print(sol2.table(tables=["constants"]))
print(sol.diff(sol2, absdiff=True))

produces

First u
Constants
---------
u : [ 0.878     0.00168   0.341     0.465     0.312       [ft]
      0.271     0.75      0.201     0.251     0.356     ]

Second u
Constants
---------
u : [ 0.178     0.531     0.131     0.115     0.363       [ft]
      0.129     0.938     0.612     0.835     0.628
      0.402     0.143     0.0631    0.854     0.0219
      0.0785    0.874     0.353     0.075     0.965
      0.237     0.981     0.361     0.162     0.0275
      0.902     0.794     0.595     0.213     0.154     ]

Solution difference
(positive means the argument is smaller)
-------------------
u : [ +393.9%    -49.0%    -98.7%   +553.0%     -5.9%
       +55.4%    -50.4%    -59.0%    -62.6%    -43.3%
      +118.3%    +89.0%    -97.3%    -12.2%   +1460.3%
      +155.5%    -46.8%    -28.9%   +315.9%    -63.1%
      +271.0%    -72.4%    -99.5%   +362.3%   +1140.7%
       -77.7%    -41.5%    -57.8%    +46.2%   +131.0%   ]
v : [ +393.9%    -49.0%    -98.7%   +553.0%     -5.9%
       +55.4%    -50.4%    -59.0%    -62.6%    -43.3%
      +118.3%    +89.0%    -97.3%    -12.2%   +1460.3%
      +155.5%    -46.8%    -28.9%   +315.9%    -63.1%
      +271.0%    -72.4%    -99.5%   +362.3%   +1140.7%
       -77.7%    -41.5%    -57.8%    +46.2%   +131.0%   ]

Absolute solution difference
----------------------------
u : [  0.7      -0.26     -0.13      0.64     -0.021      [ft]
       0.071    -0.47     -0.36     -0.52     -0.27
       0.48      0.13     -0.061    -0.1       0.32
       0.12     -0.41     -0.1       0.24     -0.61
       0.64     -0.71     -0.36      0.59      0.31
      -0.7      -0.33     -0.34      0.099     0.2      ]
v : [  0.21     -0.079    -0.039     0.19     -0.0065     [m]
       0.022    -0.14     -0.11     -0.16     -0.083
       0.14      0.039    -0.019    -0.032     0.097
       0.037    -0.12     -0.031     0.072    -0.19
       0.2      -0.22     -0.11      0.18      0.096
      -0.21     -0.1      -0.1       0.03      0.062    ]

Solution sensitivity delta
--------------------------
The largest sensitivity delta is +0

trying to compare the absolute deltas to the constants tables will show you just how these are broadcast :p

[pgkirsch]

Woop thank you! And +1 for supporting "tiled"(?) comparisons.