odlgroup / odl Goto Github PK
View Code? Open in Web Editor NEWOperator Discretization Library https://odlgroup.github.io/odl/
License: Mozilla Public License 2.0
Operator Discretization Library https://odlgroup.github.io/odl/
License: Mozilla Public License 2.0
I know that for large problems this is a BAD idea. However, for small test problems this might be a good things, especially to check things like conditioning of the problem (in my case it would be a 510x510 matrix, so it should be possible to handle).
What lead to the thought was a discussion with @ozanoktem. I have a nonlinear operator and want to solve f(x) = 0
. I'm trying to do this with Newton's method, however the CG used there does not seem to find a solution to gradf(x)p = -f(x)
, although gradf(x)
is a positive definite operator for all x
. It would therefore be useful to be able to investigate the matrix, to see if it is ill-conditioned.
As discussed in #21, we need the following two operators in default_ops
def ConstantOperator(odl.Operator):
def __init__(self, x):
self.x = x
def _call(self, _):
return self.x
def ResidualOperator(odl.Operator):
def __init__(self, op, x):
self.op = op
self.x = x
def _call(self, y):
return self.op(y) - self.x
we should add a way for users to determine the version of odl
>>> import odl
>>> odl.ODL_MAJOR_VERSION
0
>>> odl.ODL_MINOR_VERSION
9
Just a thought, currently many numpy users do things like:
a = x.dot(A.dot(x))
In odl this would be
a = x.inner(A(x))
or using the Vector.T
syntax
a = x.T(A(x))
do we want to somehow support the numpy
syntax? This could also allow odl
to be used with duck-typing in cases where a code expects a numpy array.
After the discussions we had about functionals and the properties they have in the light of (convex) optimization, the question is how to solve this in code. I'll make a quick list of what is needed, feel free to edit:
Suppose phi:X --> R
is a functional, x in X
grad phi(x)
= element in X
. This is special to Hilbert spaces - in Banach spaces it can be identified with, but is not equal to, a dual space element. Implementation-wise, gradient
will probably be an Operator
from the functional domain to itself.H_phi(x)
= linear functional X --> R
. Here it becomes more interesting: should it be something bi-linear (x, y) --> H_phi(x, y) in R
or rather an Operator
mapping a domain element to a linear functional?phi*
= functional X* --> R
. We need to know that phi
maps to R
. This can also be defined on general Banach spaces, probably harder to find closed forms though. Unless we introduce the notion of a dual space, we can only support Hilbert spaces (or maybe a few other cases where the dual can be identified with a simpler space, e.g. (L^p)* ~ L^q).phi^
= functional in X
?The fact that we need the real numbers as range would cause slight headaches in the current Operator
framework, and most definitions only seem to be relevant for functionals. So I don't really think it makes sense to add this stuff to the Operator
class. Better to make a subclass which gets these additional attributes.
I added intersphinx so that we can link documentation against numpy, it does however seem like Read The Docs has a old version of intersphinx that crashes, we should somehow fix this.
It seems having such a file is not pythonic, instead users should simply do
$ py.test
which should run the "usual tests". We need to update the readme aswell if we do this.
From #36 it seems like the cg method (found here) can get a division-by-zero error.
I think (but I'm not sure) that the problem is that if to many iterations are used, the exact solution is found. In this case, in the next iteration the search direction p
found is the zero-direction, and thus p.inner(Ap)
is zero, which is used in the code as alpha = sqnorm_r_old / p.inner(Ap)
.
Somehow 98 % of our issues end up having the proposal
label. Slightly pointless. My suggestion is that we follow the commit message label distinction made in NumPy. I'll post it here so we don't have to follow a link:
API: an (incompatible) API change
BLD: change related to building numpy
BUG: bug fix
DEP: deprecate something, or remove a deprecated object
DEV: development tool or utility
DOC: documentation
ENH: enhancement
MAINT: maintenance commit (refactoring, typos, etc.)
REV: revert an earlier commit
STY: style fix (whitespace, PEP8)
TST: addition or modification of tests
REL: related to releasing numpy
Out of these, I'd say we take
api change
bug
(already there)deprecation
documentation
maintenance
style
testing
release
new feature
or feature request
, (more specific)And I would also vote for labels for different parts of the library,
core
solvers
trafos
cuda
I'll label this as a proposal despite the irony.
Since operators like the ResidualOp
in the solver tests probably occur rather often, I suggest that we support adding an operator and an element in its range to create a new AffineOperator
(may suggest linearity) or ShiftedOperator
whose derivative is the one of the original operator, of course.
In the setup.py
there is a version number given. I think it would be useful, especially later on, to have this version number also explicit in the README
file (or somewhere else, easily accessible).
The off-the-shelf gitwash documentation is a good start, but some aspects the work flow are a bit unspecific or missing. We need to adapt those sections, e.g. coding style, commit messages.
Solvers are currently absent from docs.
We need a buildserver that is configured to run cu_ntuples
According to this travis, the test test_matrix_inner occationally fails.
This is because the result in almost_equal
is a float
, and our automated precision selector will select a higher precision (since float is 64 bit), but the calculation was actually in 32 bit!
I suggest that we do the following:
We make odlpp
a proper sub-package. I.e. we move all odl/space/cu_ntuples
and test/space/cu_ntuples
out of this repo and into the current odlpp
repo.
We then change odlpp
to odl-cuda
Finally we make odlpp
a Namespace package. What this means is that the users can install it separately, and it will appear just as now under odl.space.cuda
.
What do we think?
This has started to fail:
================================== FAILURES ===================================
_________ [doctest] odl.discr.discr_mappings.RawGridCollocation._call _________
245
246 >>> coll_op = RawGridCollocation(funcset, grid, rn, order='F')
247 >>> coll_op(func_elem)
248 Rn(6).element([-2.0, -1.0, -3.0, -2.0, -4.0, -3.0])
249
250 Use all free variables in functions you supply, otherwise
251 the automatic broadcasting will yield a wrong shape:
252
253 >>> func_elem = funcset.element(lambda x, y: 2 * x)
254 >>> coll_op(func_elem)
Expected:
Traceback (most recent call last):
...
ValueError: input shape (2,) not broadcastable to shape (6,).
Got:
Rn(6).element([2.0, 4.0, 2.0, 4.0, 2.0, 4.0])
E:\Github\ODL\odl\discr\discr_mappings.py:254: DocTestFailure
==================== 1 failed, 295 passed in 4.28 seconds =====================
My suggestion is that we should not have exception tests in a doctest.
Keeping make html
warning free enables us to notice when we get new warnings. We should thus try to make it warning free. I am currently on my way to do this.
This file is autogenerated and is way to long, we should be able to clear up a lot of these things.
All the cool kids seem to have a logo of some sort, we need one aswell. Any input on what would look cool would be cool.
We need one big logo (to use in presentations etc). And one small logo to use in webpage thumbnails and such.
There are a couple of commit messages like "merge" or "bugfix" in the commit history which are actually worse than no message at all. We need to become better at this if we want our Git history to be searchable and useful.
As a very good reference, one can once again at how NumPy does it. I suggest that we follow the same nomenclature, at least in the final versions of pull requests. Intermediate commits on branches may have bad messages, as long as they are cleaned up with git rebase -i
before the pull request is sent (or at the very latest before it is merged into master
).
Short summary:
BUG
for a bugfix (see the table in the link).add feature xyz
instead of added feature xyz
). A very simple reason for this style is that you save on average 2 characters for something else.
It's become relatively common that we check in things that break existing tests and examples to the main branch. This is OK right now, but as we proceed this will not be feasible.
I therefore propose that we (in some way) create a build server. This could preferably be run on some decently good computer at KTH, and should preferably run the tests with different configurations (Windows, Linux, python 2, python 3).
The merging of solvers has caused several warnings in the doc, these should be fixed.
Copied from odl-solvers issue 5
I suggest that we write a common class for step length calculation which tells the algorithm what it needs (which part of the history), such that an algorithm can determine which iterates and gradients (what else?) to keep.
A simple approach would be to give the StepLength
class attributes
iter_idcs
(iterates)func_idcs
(function evaluations)grad_idcs
(gradient evaluations / approximations)hess_idcs
(Hessian evaluations / approximations)step_idcs
(step lengths)which all return lists of indices of items to keep. An algorithm could manage lists of all these objects and pass its lists (possibly kept in a dictionary) to the step length call, which in turn does its thing to calculate the step length.
We need a public documentation page. This would include both the odl programming doc, but also issues #2, #4, #5. First we need to decide on where to put this documentation, some options are
Since we have decided to use Sphinx for the documentation, all of these are ofc compatible, and we could for example use all of them. For maintainance, I propose that we try to keep this at a minimum, at least until we get stable releases.
As noted in implementing #53 ProductSpace.__getitem__
currently returns a list if indices
is a slice
:
return self.spaces[indices]
In the case where indices
is a slice
, ProductSpace.__getitem__
should do:
return ProductSpace(self.spaces[indices])
to be sure, this should fail if the space has a custom norm (since this does not transfer trivially).
if indices
is a int
, we should simply
return self.spaces[indices]
Moved from https://github.com/odlgroup/odl-solvers/issues/3
I initialized a structure, but it is obviously far from optimal. We need to find sensible categories for solvers, e.g.
We need to reflect these categories in classes in a sensible manner.
With :automodule:
it seems class magic methods, most notably __init__
does not get documented. This causes some obvious issues.
As seen in this CI log the optimizer test occasionally fails.
When reading the discussion in #49, I realized that we could use a "project on subspace" operator
class Projection(odl.Operator):
def __init__(self, space, indices):
assert isinstance(space, ProductSpace)
self.indices = indices
super().__init__(space, space[indices], linear=True)
def _call(self, x):
return x[self.indices]
For some reason, we needed the requirements in text files instead of a string in setup.py
. What about extra dependencies? For example, the Fourier transform operator will depend on pyfftw
, and the wavelet transform will use PyWavelets
. I'd want to be able to pull the dependencies by writing
pip install --user -e .[fft,wavelets]
This actually works and asks for the extras_require
dictionary for the dependencies, so the dependency info goes there.
Question is: also into a file? Or in a string? Or is it perhaps possible to create extra sections in the requirements file? I could not find anything on the web treating the last topic.
Currently the following gives you 50 figures:
for i in range(50):
#do stuff
x.show()
a syntax like
fig=None
for i in range(50):
#do stuff
fig=x.show(figure=fig)
or possibly something better, would be very nice.
See odlgroup/odl-solvers#7 and odlgroup/odl-solvers#9.
I think that in the long run, the solvers need more structure than can be handled in a single file. They will probably also add dependencies on third-party software which we don't want in the core ODL package. So unless there are strong reasons against it now, I'll move the code to the odl-solvers repository.
Copied from odl-solvers issue 4
Likewise step length rules, it will be equally useful to define a common interface for termination rules of optimization schemes, provided that there is some common structure.
I have somewhat lost track about what expression in which combinations of operators, scalars and vectors are currently supported and what they mean mathematically.
I suggest that we write up a documentation chapter about operator expressions, possibly summarized as a table containing (example: operator * scalar
)
op * scal
x --> op(scal * x)
scal in op.domain.field
OperatorRightScalarMult
We should gather common questions and put together a FAQ. I started one in 5fe6d3a.
The docstrings in the __init__.py
files do not end up in the Sphinx documentation. I think they should, since they are a good place to put an introduction to the subpackage contents.
In contrast, I'd argue that one only writes a minimal docstring at the top of a module so that one does not have to scroll all the way down to see the code.
We need a way for users to install odl
without needing to point them towards git
.
Since we're building a new package, we have the luxury of choosing whether to have unicode literals enabled by default or not. Since they are standard in Python 3, I strongly suggest that we make an attempt at it.
A known issue is the ABCMeta
class complaining when unicode_literals
are imported. Back then, I tried to solve it by explicitly passing a byte string, but that didn't help. There is, however, good documentation available which may point to an elegant solution of this problem.
If some external module really screws up and cannot be convinced at all to work, it is still possible to remove unicode_literals
from imports selectively.
We need something that users can look at to see where to get started, similar to numpys quickstart
Lets make a short list of what this should contain (obviously inspired by numpy)
LinearSpace
, quite obviously Rn
, show that you can do what you'd expect with it.
element
, zero
.+
, -
, /
, *
), indexinglincomb
, multiply
.numpy
. asarray
, __array_wrap__
to use ufunc
'sSet
, I'd prefer Interval(0, 1)
, do some simple things with it.
__contains__
L2(Interval(0, 1))
l2_uniform_discretization(L2(Interval(0, 1)))
Rn
/CudaRn
possibilityodl.diagnostics
possibilities.Operator
on the space
Operator
Let's try to find out if there is a good and predictable way to move documentation from _apply()
and _call()
of Operator
implementations to the respective __call__()
method so it appears in the right place in the final doc. In principle, a script would for each subclass of Operator
cls._apply.__doc__
cls.__call__.__doc__
out : something
to out : something, optional
out
is provided as argument, it is returnedcls._apply.__doc__
and cls._call.__doc__
with a hint not to use the functions directly but rather call operator(x)
or operator(x, out)
, respectivelyThe question is where to put this functionality. Maybe one could implement it in Operator.__new__()
?
I'll do a simple test with Travis CI.
I tired to push to the repo and got the error message:
remote: Permission to odlgroup/odl.git denied to aringh.
fatal: unable to access 'https://github.com/odlgroup/odl.git/': The requested URL returned error: 403
Do I miss some rights, or is there somethings else wrong?
We need a build server that runs windows.
The build server should do pep8 checking and report fail if it does not work.
Something screwed up the build process of ReadTheDocs. The API part is empty.
We need some advice for people who want to contribute but do not exactly know how. Intended audience are users and developers. We should describe
There are a few Operator
's on product spaces that we could standardize and that would be usefull, see #239 for example.
The ones i can think of are the following (naming obviously TBD).
The distributive operator, which calls a series of "sub" operators given the same argument. Think of the matrix:
[A,
B]
DistributiveOperator(A,B)(x) = [A(x), B(x)]
The "diagonal" operator, which calls a series of operators "per space", think of the matrix
[A, 0,
0, B]
DiagonalOperator(A,B)((x,y)) = [A(x), B(y)]
Finally we can think of calling a sequence of operators with given (different) arguments and somehow reducing the result. Reduction by summation is the obvious choice, but other options could work, think of the matrix:
[A, B]
SumOperator(A,B)((x, y)) = A(x) + B(y)
All of these operators could obviously be generalized to tuples of n
operators, for example:
DistributiveOperator(A, B, C, D)(x) = [A(x), B(x), C(x), D(x)]
I think we should add these to odl.operator.default_ops
, but first I need some input on the names.
The first one ("Distributive") is called vstack in numpy
The second one seems to be scipy.linalg.block_diag
and the last one is hstack
The more I think of it, it would maybe be a good idea to not implement three separate classes for this, but instead to create one base class and then (maybe) have these as special cases. Something like a sparse matrix class (but for operators). That way the user could do complicated shit like
op = GeneralProductOperator([[A , None, None, B ],
[A , None, C , D ],
[A , C , None, None]])
which would have the obvious meaning:
op([x, y, z, w]) == [A(x) + B(w), A(x) + C(z) + D(w), A(x) + C(y)]
This is similar to the scip.sparse.bmat method.
Beside the pure implementation, we need to review one design decision, namely where the discrete space functions are specified. Currently, this is completely done in the data space Fn
by adding a weight to the inner product, norm and dist.
If we decide to keep it that way, we must add an exponent
parameter to the weighting classes, and for p != 2
, things become slightly more complicated:
inner
must be disablednorm
must be re-implemented, e.g. in the matrix-weighted case, the matrix power M^(1/p)
has to be (pre-)computed since the weighted norm is |||x|||_p = ||M^(1/p)x||_p
. Probably this is not feasible for larger problems, we'll see. For vector- and constant-weighted spaces, it's easy though.A declarative, efficient, and flexible JavaScript library for building user interfaces.
๐ Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. ๐๐๐
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google โค๏ธ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.