Skip to content
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

Tensor spaces, second try #861

Closed
wants to merge 90 commits into from
Closed

Tensor spaces, second try #861

wants to merge 90 commits into from

Conversation

kohr-h
Copy link
Member

@kohr-h kohr-h commented Jan 25, 2017

Copied from #753

TODOs:

  • Implement basic structures for tensors and tensor spaces
  • Replicate functionality of NtuplesBase, FnBase, NumpyNtuples and NumpyFn in the corresponding tensor classes
  • Remove Ntuples and related stuff
  • Add former Ntuples and the likes as special case to tensor classes
  • Change names with "tensor" to something more lightweight
  • Find place for MatVecOperatorMatrixOperator, see Reload modules in Spyder causes bugs #584. Solved in Issue 910 matvec op #911
  • Improve pretty-printing of elements
  • Fix broken downstream code (lp_discr, ...)
  • Change code that was restricted to 1D data containers
  • Port unit tests
  • Write some doctests
  • Write more doctests
  • Remove TensorSet and others (see Scrap FunctionSet, DiscretizedSet and similar #860) duck-type space properties
  • Check documentation for outdated information, e.g. regarding 1D data layout (Update: can't find anything explicit)
  • Use notion "rank" instead of "number of dimensions" to better distinguish between space.ndim and space.size Not doing this
  • Implement indexing of TensorSpace. Done, but kind of wrong, see here and here Done properly now
  • Write extra documentation on tensors Should be part of Representing vector- and tensor-valued functions #908
  • Port show methods (1D from Ntuples, 2D from DiscreteLp)
  • Add unit tests for unusual data types and the likes
  • Simplify doctests of tensor spaces
  • Fix FunctionSpaceElement in 1D (again)
  • Implement tensor-valued functions, see Representing vector- and tensor-valued functions #908
  • Fix interpolation from array of points, related to Representing vector- and tensor-valued functions #908
  • Cleanup (old stuff, commit history, ...)
  • Fix large-scale tests
  • Make test with power function more general
  • Add more tests for product space multi-indexing
  • Implement indexing of FunctionSpace as
    • indexing the components of the output (i.e. raise an error for scalar-valued functions) in __getitem__
    • by axis of the domain using a byaxis attribute
  • Implement indexing of DiscreteLp as in FunctionSpace
  • Implement __array_ufunc__ interface on Tensor and NumpyTensor (the latter only differs in weight propagation)
  • Implement __array_ufunc__ interface on DiscreteLpElement
  • Write more extensive tests for __array_ufunc__ intefaces (bad input etc.)

@kohr-h
Copy link
Member Author

kohr-h commented Jan 26, 2017

Implement indexing of TensorSpace

Regarding this topic: I have implemented a simple indexing capability for spaces, but it is nowhere near as powerful and flexible as the lazy approach of simply indexing the array and appropriately wrap whatever comes out.
So the earlier suggestion

arr = self.data[indices]
return self.space[indices].element(arr)

is much much less powerful than lazy indexing, if we don't put great amounts of work into it.

So I'd say the space indexing can be used for simple things like

>>> space = odl.rn((2, 3, 4))
>>> space[[0, 2]]
rn((2, 4))
>>> space[0]
rn(2)

but not for indexing of elements.

@kohr-h
Copy link
Member Author

kohr-h commented Jan 26, 2017

Regarding this topic: I have implemented a simple indexing capability for spaces, but it is nowhere near as powerful and flexible as the lazy approach of simply indexing the array and appropriately wrap whatever comes out.
So the earlier suggestion

arr = self.data[indices]
return self.space[indices].element(arr)

is much much less powerful than lazy indexing, if we don't put great amounts of work into it.

So I'd say the space indexing can be used for simple things like

space = odl.rn((2, 3, 4))
space[[0, 2]]
rn((2, 4))
space[0]
rn(2)

but not for indexing of elements.

I realized that indexing like this is not really helpful since it's not compatible with power spaces. I'll think about it a bit more.

@adler-j
Copy link
Member

adler-j commented Feb 8, 2017

Is this in a state where it needs a review? If so I guess I need to get at it.

@kohr-h
Copy link
Member Author

kohr-h commented Feb 8, 2017

Is this in a state where it needs a review? If so I guess I need to get at it.

Yes it is. I deliberately didn't mention it so nobody burns time at it that they don't have :-)

@adler-j
Copy link
Member

adler-j commented Feb 9, 2017

Could you do the rebase from master first? It seems rather big.

@kohr-h
Copy link
Member Author

kohr-h commented Feb 9, 2017

Could you do the rebase from master first? It seems rather big.

Yes, I've also spotted some glitches that I want to fix along the way.

@kohr-h
Copy link
Member Author

kohr-h commented Feb 10, 2017

Rebased

Copy link
Member

@adler-j adler-j left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Huge review done, but mostly small things. I have not reviewed the tests yet, will do that later on.

On the larger scale, I have a few comments:

  • The support for DiscreteLp with vector/tensor valued functions needs to be added and properly supported. This is the largest out-standing problem with this.
  • That all examples default to shape=(2,3) makes the doc much harder to read and adds little. I understand that during development this is important, but in release we should trim this in general, just like we dont use dtype='float32' everywhere.
  • You have added __hash__ in a few places where it shouldn't, i.e. in mutable objects. Also I think you may have missed adding __ne__ in some places, check that.
  • The matrix_representation/MatrixOperator thing needs to be fixed, but could use a new issue.


x0 = np.arange(fn.size)
x = fn.element(x0)
x0 = np.arange(tspace.size)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

tspace.shape?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Haven't looked much into this file, probably needs more love

Adding Fn spaces
----------------
The abstract spaces `FnBase` and `NtuplesBase` are the workhorses of the ODL space machinery. They are used in both the discrete :math:`R^n` case, as well as data representation for discretized function spaces such as :math:`L^2([0, 1])` in the `DiscretizedSpace` class. These are in general created through the `rn` and `uniform_discr` functions who take an ``impl`` parameter, allowing users to select the backend to use.
Adding Tensor spaces
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Skipping the doc for now.

@@ -63,26 +63,26 @@


def import_submodules(package, name=None, recursive=True):
""" Import all submodules of a module, recursively, including subpackages
"""
"""Recursively import all submodules of ``package``."""
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Improve. I guess it returns the names etc?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

will check

@@ -11,7 +11,7 @@ Glossary

array-like
Any data structure which can be converted into a `numpy.ndarray` by the
`numpy.array` constructor. Includes all `NtuplesBaseVector` based classes.
`numpy.array` constructor. Includes all `Tensor` based classes.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The NtuplesBaseVector -> Tensor may be our best rename in a long time!

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah scrap those shitty old names!

@@ -22,7 +22,7 @@ Glossary

discretization
Structure to handle the mapping between abstract objects (e.g. functions) and concrete, finite realizations.
It encompasses an abstract `Set`, a finite data container (`NtuplesBaseVector` in general) and the mappings between them, :term:`sampling` and :term:`interpolation`.
It encompasses an abstract `Set`, a finite data container (`GeneralizedTensor` in general) and the mappings between them, :term:`sampling` and :term:`interpolation`.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Didn't we scrap GeneralizedTensor?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, needs to be replaced

return ProductSpace(*self.spaces[indices], field=self.field)
elif isinstance(indices, tuple):
if len(indices) > 1:
raise ValueError('too many indices: {}'.format(indices))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Another option is to pass these through further down, i.e.

return ProductSpace(self.spaces[indices[0]][indices[1:]], field=self.field)

@@ -848,7 +854,7 @@ def show(self, title=None, indices=None, **kwargs):
--------
odl.discr.lp_discr.DiscreteLpElement.show :
Display of a discretized function
odl.space.base_ntuples.NtuplesBaseVector.show :
odl.space.base_tensors.GeneralizedTensor.show :
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

find replace GeneralizedTensor

impl : string, optional
The backend to use. See `odl.space.entry_points.NTUPLES_IMPLS` and
`odl.space.entry_points.FN_IMPLS` for available options.
order : {'C', 'F'}, optional
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

'A'?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Clearly missing

`odl.space.entry_points.FN_IMPLS` for available options.
order : {'C', 'F'}, optional
Axis ordering of the data storage.
impl : str
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

optional

if self.const <= 0:
raise ValueError('expected positive constant, got {}'
''.format(const))
if not np.isfinite(self.const):
raise ValueError('`const` {} is invalid'.format(const))

def __hash__(self):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Prefer if __hash__ is always right after __eq__ since they are so strongly connected.

@kohr-h
Copy link
Member Author

kohr-h commented Feb 13, 2017

Thanks a lot for the detailed review @adler-j, will try to get this done asap.

"""Weighting of a `NumpyTensorSpace` with constant 1."""

# Implement singleton pattern for efficiency in the default case
__instanc
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You use a 3 dimensional array later.

Anyway, if we allow 2+, why not allow 1d as well, in which case this is an inner product.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Wrong place again, you need to hack it apparently.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The comment is w.r.t MatrixOperator where you assume the matrix is 2 or higher dimensional. I can easily see this working with 1 or even 0 dimensional arrays, so that condition is arbitrary.

The line "You use a 3 dimensional array later." refers to a comment that MatrixOperator input should be 2d, but you later use 3d. Also you later say that axis needs to be 1 or more, but the default is 0.

@adler-j
Copy link
Member

adler-j commented Mar 28, 2017

What is the review status on this?

@kohr-h
Copy link
Member Author

kohr-h commented Apr 1, 2017

What is the review status on this?

I figured that sampling would need tensor-valued functions, so I implemented them. Now I'm in the middle of writing/fixing the tests for those.

@kohr-h kohr-h mentioned this pull request Jul 10, 2017
@adler-j
Copy link
Member

adler-j commented Jul 14, 2017

Can you re-open this for a new review, if you want one, else we really want to go ahead here.

@kohr-h kohr-h mentioned this pull request Aug 7, 2017
@kohr-h
Copy link
Member Author

kohr-h commented Aug 10, 2017

After thinking about it some more, I think the following way of indexing spaces is most reasonable:

  • NumpyTensorSpace:
    • __getitem__: Numpy indexing, should wrap an array indexed in the same way
    • byaxis: index axes of the tensors
  • FunctionSpace:
    • __getitem__: Numpy indexing for the out_shape part, i.e., for tensor-valued functions. Raise an error for scalar-valued functions.
    • byaxis: index the domain per axis
  • DiscreteLp: same as FunctionSpace

Opinions?

@kohr-h
Copy link
Member Author

kohr-h commented Aug 11, 2017

The page for this PR has become unworkable once again, opening a new one.

@kohr-h kohr-h closed this Aug 11, 2017
@kohr-h kohr-h mentioned this pull request Aug 11, 2017
29 tasks
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

2 participants