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

Fixes bug where multiple identical indexed slices are ignored #1361

Merged
merged 2 commits into from
Sep 25, 2017

Conversation

AllenHW
Copy link
Contributor

@AllenHW AllenHW commented Sep 21, 2017

Motivation and context:
As per #947, we have a bug where advanced indexing works when the indices are different but does not work when same index repeats multiple times. Eric found the api .at(), and I use it to fix the bug.

How has this been tested?
I added a test test_set_weight_solver in tests/test_connection.py.

How long should this take to review?

  • Quick (less than 40 lines changed or changes are straightforward)

Types of changes:

  • Bug fix (non-breaking change which fixes an issue)

Checklist:

  • I have read the CONTRIBUTING.rst document.
  • I have updated the documentation accordingly.
  • I have included a changelog entry.
  • I have added tests to cover my changes.
  • All new and existing tests passed.

Still to do:

@AllenHW AllenHW self-assigned this Sep 21, 2017
@AllenHW
Copy link
Contributor Author

AllenHW commented Sep 22, 2017

Some of the tests are failing but can't reproduce them locally because I don't have the same numpy version (I should get virtualenv..). Not sure why the case with numpy 1.8 is failing because .at() supported by that version. Will look into it

@AllenHW
Copy link
Contributor Author

AllenHW commented Sep 22, 2017

It looks like the issue was that for numpy 1.8, .at() does not take Ellipsis as a slice. So for instance np.add.at(np.array([1,2]), Ellipsis, np.array([1,1])) doesn't work. Hence I had to add a silly branch statement where if the slice is Ellipsis we fall to the previous method.
For numpy >=1.9 this is not an issue. So the branch condition is just to support numpy 1.8

@@ -388,7 +388,10 @@ def make_step(self, signals, dt, rng):

if inc:
def step_copy():
dst[dst_slice] += src[src_slice]
if type(dst_slice) == type(Ellipsis):
Copy link
Contributor Author

Choose a reason for hiding this comment

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

type(dst_slice) == type(Ellipsis) if no slice is provided. So we'd think that dst = dst[dst_slice]
Yet, dst += src[src_slice] sometimes fails because "dst referenced before assignment" because of some weird numpy internal

Copy link
Collaborator

Choose a reason for hiding this comment

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

I think isinstance(dst_slice, Ellipsis) is usually preferred for type checks (except when derived types should be explicitly excluded which is very rare).

Copy link
Collaborator

Choose a reason for hiding this comment

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

Also, the slices should not change during the simulation. That should allow to move the if-branching outside of the step_copy function. That might have some performance impact, but no idea whether that's significant.

I'm also curious if there is a performance impact replacing += with at.

Copy link
Member

Choose a reason for hiding this comment

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

I agree with Jan, the slices won't change during the simulation so the if branch should be outside the step_copy function. The step functions returned in the builder are what are actually run when you do sim.run, so they should be written to be as fast as possible.

Additionally, it should only be necessary to use at when dst_slice or src_slice are lists (or numpy arrays) of integers, correct? You could check for that also before defining the step function so that we only use at when we need to.

Copy link
Collaborator

Choose a reason for hiding this comment

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

If we really want to be specific, the only time we need at is when dst_slice is a list or ndarray with repeated indices. We could check for this with something like len(np.unique(dst_slice)) < len(dst_slice).

Copy link
Member

@tbekolay tbekolay left a comment

Choose a reason for hiding this comment

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

Looks good! Made some inline comments for minor changes.

@@ -388,7 +388,10 @@ def make_step(self, signals, dt, rng):

if inc:
def step_copy():
dst[dst_slice] += src[src_slice]
if type(dst_slice) == type(Ellipsis):
Copy link
Member

Choose a reason for hiding this comment

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

I agree with Jan, the slices won't change during the simulation so the if branch should be outside the step_copy function. The step functions returned in the builder are what are actually run when you do sim.run, so they should be written to be as fast as possible.

Additionally, it should only be necessary to use at when dst_slice or src_slice are lists (or numpy arrays) of integers, correct? You could check for that also before defining the step function so that we only use at when we need to.

d_data = sim.data[d_probe]

plt.plot(t, a_data)
assert np.allclose(a_data[100:], [0], atol=0.3)
Copy link
Member

Choose a reason for hiding this comment

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

These tolerances feel pretty big to me; I'll look at the plots when I rereview, but try a larger synapse value on the probes to see if the decoded value is more stable?

def test_advanced_indexing(Simulator, plt, seed):
N = 100

model = nengo.Network()
Copy link
Member

Choose a reason for hiding this comment

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

You should pass the seed into the Network (Network(seed=seed)) so that the test becomes deterministic; otherwise, we'll get random parameters each time and potentially random test failures! 😄

plt.plot(t, d_data)
assert np.allclose(c_data[100:], [-1,1], atol=0.3)
plt.plot(t, d_data)
assert np.allclose(d_data[100:], [1,1], atol=0.3)
Copy link
Member

Choose a reason for hiding this comment

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

This test is great, thanks!

We mentioned in the dev meeting that it would also be nice to support boolean arrays for indices. Have you tried testing with boolean indices? If not, it might be worth trying it out to see if it works... if it doesn't, it makes sense to implement it in another PR, but if it just works then it'd be good to add another test for that here!

Copy link
Collaborator

Choose a reason for hiding this comment

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

I tried this and it seems to work. I'll add a test for it.

@hunse
Copy link
Collaborator

hunse commented Sep 22, 2017

I did a test, and at seems to be considerably slower than += (at least three times slower, and much more if dst_slice is a slice or Ellipsis or something). So I think we should only use at if dst_slice has repeated indices.

Here's a script if anyone wants to play around.

import timeit
import numpy as np

m = int(1e6)

if 0:
    n = m
    inds = Ellipsis
elif 0:
    n = int(0.9*m)
    inds = slice(5, 5+n)
else:
    n = int(0.9*m)
    inds = np.random.permutation(m)[:n]
    # inds = np.sort(inds)

a = np.random.uniform(-1, 1, size=n)
b = np.zeros(m)

t_inc = min(timeit.repeat(
    'b[inds] += a', 'from __main__ import *; b[:] = 0', repeat=3, number=1))
t_at = min(timeit.repeat(
    'np.add.at(b, inds, a)', 'from __main__ import *; b[:] = 0', repeat=3, number=1))

print(t_inc, t_at)
print(t_at/t_inc)

if type(dst_slice) == type(Ellipsis):
dst[dst_slice] += src[src_slice]
else:
np.add.at(dst, dst_slice, src[src_slice])
else:
def step_copy():
dst[dst_slice] = src[src_slice]
Copy link
Collaborator

Choose a reason for hiding this comment

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

In this case, if there are repeated indices in dst_slice, then the behaviour is kind of undefined. I think we should have an assertion here to make sure that's not the case. I don't think a user should ever be able to trip that assertion, but it could save a developer some trouble.

@hunse
Copy link
Collaborator

hunse commented Sep 22, 2017

I'm going to go through and make some of the changes I suggested.

@AllenHW
Copy link
Contributor Author

AllenHW commented Sep 23, 2017

The below suggests that use of at is faster for repeated indices than creating say multiple connections

import timeit
import numpy as np

m = int(100)

b = np.zeros(m)
c = np.zeros(m)
inds = np.zeros(m, dtype=np.int)

t_inc = min(timeit.repeat(
    'b[[0]] += [1]', 'from __main__ import *;', repeat=m, number=m))
t_at = min(timeit.repeat(
    'np.add.at(c, inds, [1])', 'from __main__ import *;', repeat=1, number=1))

print(t_inc, t_at)
print(t_at/t_inc)

LGTM :shipit:

Copy link
Collaborator

@hunse hunse left a comment

Choose a reason for hiding this comment

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

This all LGTM, assuming that @tbekolay is able to tighten up those tolerances when he goes through it.

@tbekolay tbekolay force-pushed the advanced-index branch 2 times, most recently from b36bdb7 to 9bc1e22 Compare September 25, 2017 16:20
AllenHW and others added 2 commits September 25, 2017 12:43
Also test that boolean indexing works. With this, advanced indexing
is fully supported for ObjViews in connections.

Fixes #947.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Development

Successfully merging this pull request may close these issues.

None yet

4 participants