-
Notifications
You must be signed in to change notification settings - Fork 525
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Pynumero.sparse Updates #1285
Pynumero.sparse Updates #1285
Conversation
- more strict dimension requirements in BlockMatrix - removed BlockSymMatrix and MPIBlockSymMatrix - adding some tests for BlockMatrix - more careful copying within BlockMatrix methods - Simplified BlockMatrix.copyfrom significantly; hopefully these changes will also improve performance, but I don't have any data to support this - Improved caching of block row and block column dimensions within BlockMatrix
@santiagoropb Can you also please take a look at this PR? |
Codecov Report
@@ Coverage Diff @@
## master #1285 +/- ##
=======================================
Coverage 69.96% 69.96%
=======================================
Files 536 536
Lines 82759 82759
=======================================
Hits 57903 57903
Misses 24856 24856 Continue to review full report at Codecov.
|
@michaelbynum sure! I'll take a look during the weekend. From a quick look looks great! Thank you! One quick thing I noticed is that test_nlp_composition.py was partially updated to call toarray instead of todense. That file is probably not run after Carls refactor but if you want to keep it there I suggest to update the set and get block notation too. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This looks overall pretty good! I had some minor comments regarding slicing and matrix products.
One thing that I always wonder is if we could have a decorator on the blockvector and blockmatrix classes to use numba just-in-time compiling. The decorator should check if numba is installed and if so it should just-in-time compile. That could potentially speed up all the python loops we have every time we go through the blocks. The additional work to support that seems very minor. We can also add numba to the conda recipe to always have it installed.
@@ -961,7 +961,7 @@ def test_jacobian_g(self): | |||
B2[i, i] = 1.0 | |||
Ji[1, 0] = coo_matrix(B1) | |||
Ji[1, 1] = coo_matrix(B2) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
You may want to change these to set_block across this whole file
Structured Vector interface | ||
Structured vector interface. This interface can be used to | ||
performe operations on vectors composed by vectors. For example | ||
bv = BlockVector([v1, v2, v3]), where vi are numpy.ndarrays or BlockVectors. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This doc string is out of date. After your changes we no longer have a constructor that takes a list of blocks
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@santiagoropb Good catch. Sorry I removed this functionality from the constructor, but it made MPIBlockVector
a lot simpler.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Let me clarify. It was simpler to not allow this functionality in MPIBlockVector
and I wanted the API to be consistent across all of the block vectors and matrices.
if x1.__class__.__name__ == 'MPIBlockVector': | ||
raise RuntimeError('Operation not supported by BlockVector') | ||
if x2.__class__.__name__ == 'MPIBlockVector': | ||
raise RuntimeError('Operation not supported by BlockVector') | ||
raise NotImplementedError() | ||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I like how you call now the assert_block_structure method. Just a side comment. The asserts can be turned of with python -0. That can speed up the code since some block asserts can be expensive. Specially the block matrix ones
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
With the stricter dimension requirements, I was able to make all of the assert_block_structure
methods quite fast. All it does is check the length of a set.
self._block_mask[key] = True | ||
self._brow_lengths[key] = value.size | ||
def __getitem__(self, item): | ||
raise NotImplementedError('BlockVector does not support __getitem__.') |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think we can still support this as it was a normal array. We should check if all blocks are defined, then given the "flat index" we identify which block contains the index and we return the value. If the block is a BlockVector we recurse
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The tricky part may be when the user passes a slice. We can raise an exception or as in the scalar case gather elements from each block and put them into a numpy array
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I agree. I suggest we add an issue and revisit this later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I created a new issue.
|
||
def __le__(self, other): | ||
def __setitem__(self, key, value): |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Similar comments as in getitem
prod = m * m2 | ||
|
||
self.assertIsInstance(prod, BlockMatrix) | ||
self.assertTrue(np.allclose(flat_prod.toarray(), prod.toarray())) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
If we decide to support product with numpy array we should add a test here
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
an easy way to support multiplication with 2d arrays would be to call the constructor of scipy.sparse on the array to get a "sparse copy" and use that sparse copy within our multiplication. This is important because some times scipy returns matrices of only one entry as 2d arrays.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
we should add a test anyway to make sure the correct error is thrown.
if isinstance(other, MPIBlockVector): | ||
raise NotImplementedError('Vector-Matrix multiply not supported yet') | ||
if isinstance(other, BlockVector): | ||
raise NotImplementedError('Vector-Matrix multiply not supported yet') |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
At some point @carldlaird mentioned we could get some performance benefits if we support xTA. Not sure if thats still in his plans?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
We should support this when possible. Currently, we are doing a transpose (which copies) of the matrix which could be expensive. Can we add an issue to revisit?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I added an issue.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Most of these can be pushed to discussion, todo, or a new issue. The API changes and comment changes should be done - hopefully the "easy" things can be addressed now.
Feel free to add a single issue with a list of things that were not yet addressed.
|
||
return result | ||
|
||
def make_local_copy(self): |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
same question as above. I think I made a comment about wanting this on MPIBlockMatrix.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I created an issue for this.
self._set_block_size(key, value.size) | ||
|
||
def __getitem__(self, item): | ||
raise NotImplementedError('MPIBlockVector does not support __getitem__.') |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Should we point the user to get_block and set_block in these error messages?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Done.
# checks only the mask of one of them since all must have the same | ||
global_mask[bid] = processor_to_mask[0][bid] | ||
|
||
disp_owner = self._rank_owner[bid] if self._rank_owner[bid] >= 0 else 'A' |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Just think it is better to be explicit rather than relying on "codes" to interpret the output.
prod = m * m2 | ||
|
||
self.assertIsInstance(prod, BlockMatrix) | ||
self.assertTrue(np.allclose(flat_prod.toarray(), prod.toarray())) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
we should add a test anyway to make sure the correct error is thrown.
@@ -0,0 +1,4 @@ | |||
|
|||
|
|||
class MPISpaceWarning(Warning): |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Was this used anywhere?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Good question. I don't think so. I'll remove it.
@carldlaird and @santiagoropb I think a addressed all of your concerns. Please have a look. |
Summary/Motivation:
This PR makes updates to
pyomo.contrib.pynumero.sparse
.Changes proposed in this PR:
BlockVector
andBlockMatrix
MPIBlockVector
andMPIBlockMatrix
Legal Acknowledgement
By contributing to this software project, I agree to the following terms and conditions for my contribution: