weave.blitz produces different code for equal numpy ranges

[DavidSichau]

When blitz is used to inline numpy expressions the outcome depends on the formating of the ranges, even if they are the same. For example a different result is produced for u[0:x] and u[:x] however the result should be the same.

See also:
http://stackoverflow.com/questions/16230848/blitz-code-produces-different-output

import numpy as np
from scipy.weave import blitz

def test_blitz_bug(N=4):
    ReI = 1.2
    ux_blitz_buggy, ux_blitz_not_buggy, ux_np = np.zeros((N, N)), np.zeros((N, N)), np.zeros((N, N))
    uxold = np.random.randn(N, N)
    ux_np[0:,1:-1] = uxold[0:,1:-1] + ReI* (uxold[0:,2:] - 2*uxold[0:,1:-1] + uxold[0:,0:-2])
    expr_buggy = 'ux_blitz_buggy[0:,1:-1] = uxold[0:,1:-1] + ReI* (uxold[0:,2:] - 2*uxold[0:,1:-1] + uxold[0:,0:-2])'
    expr_not_buggy = 'ux_blitz_not_buggy[:,1:-1] = uxold[:,1:-1] + ReI* (uxold[:,2:] - 2*uxold[:,1:-1] + uxold[:,0:-2])'
    blitz(expr_buggy)
    blitz(expr_not_buggy)
    assert not np.allclose(ux_blitz_buggy, ux_np)
    assert np.allclose(ux_blitz_not_buggy, ux_np)

if __name__ == '__main__':
    test_blitz_bug()

[Jorge-C]

I've been testing this (I'm the one who answered that SO question) and I think the actual issue is that assignment inside a blitz expression to a slice of the form i: fails. A simplified test case is available here. Notice that the issue lies in the left hand side (the right hand side is the same for both expressions).

Assignment to :i works fine:

import numpy as np
from scipy.weave import blitz
import numpy.testing as npt


def test_no_bug():
    """Assignment to arr[:i] does not fail inside blitz expression."""
    N = 4
    expr_buggy = 'arr_blitz_buggy[:{0}] = arr[:{0}]'
    expr_not_buggy = 'arr_blitz_not_buggy[0:{0}] = arr[:{0}]'
    np.random.seed(7)
    arr = np.random.randn(N)
    for lim in [1, 2]:
        arr_blitz_buggy, arr_blitz_not_buggy, arr_np = np.zeros(N), np.zeros(N), np.zeros(N)
        blitz(expr_buggy.format(lim))
        blitz(expr_not_buggy.format(lim))
        arr_np[:lim] = arr[:lim]
        npt.assert_allclose(arr_blitz_buggy, arr_np)
        npt.assert_allclose(arr_blitz_not_buggy, arr_np)

Running it in scipy '0.13.0.dev-639ef30' gives a bunch of deprecation warnings, though.

[Jorge-C]

I believe I've tracked it down. The explanation is a bit wordy, but this
could be my first bug fix in Numpy!

At first, I was able to write a blitz-only example showing the problem:

#include <iostream>
#include <blitz/array.h>

int _beg = blitz::fromStart;
int _end = blitz::toEnd;
blitz::Range _all = blitz::Range::all();

int main() {
  blitz::Array<float, 1> A(4), B(4), C(4), D(4);
  A = 1, 2, 3, 4;
  B = 0, 0, 0, 0;
  C = 0, 0, 0, 0;
  D = 0, 0, 0, 0;
  // Buggy, it's similar to the code generated by the buggy blitz expr
  B(blitz::Range(0, _end)) = A(blitz::Range(0, _end));
  std::cout << B << std::endl;
  // Works OK, similar to the code generated by the non buggy expr
  C(blitz::Range(0, 3)) = A(blitz::Range(0, _end));
  std::cout << C << std::endl;
  // Works fine!? Seems equivalent to the buggy version...
  D(blitz::Range(0, blitz::toEnd)) = A(blitz::Range(0, _end));
  std::cout << D << std::endl;
}

From that, it was clear that blitz::toEnd had a different value
outside the main function, so the example could be simplified:

#include <iostream>
#include <blitz/array.h>
#include <climits>

/*
  In the blitz version shipped by scipy, `blitz::toEnd` is defined in
  https://github.com/scipy/scipy/blob/master/scipy/weave/blitz/blitz/range.h#L44
  as `INT_MIN` from `climits` or `limits.h`.

  In the latest version from sourceforge,
  it is defined as `const int toEnd = INT_MAX;`
*/
int _end = blitz::toEnd;
int mx = INT_MAX;

int main() {
  std::cout << "INT_MAX " << mx << std::endl;
  std::cout << "_end    " << _end << std::endl;
  std::cout << "toEnd   " << blitz::toEnd << std::endl;
  /*
    Output (blitz 0.10)
    INT_MAX 2147483647
    _end    0
    toEnd   2147483647

    Output (shipped by scipy)
    INT_MAX 2147483647
    _end    0
    toEnd   -2147483648
  */
}

After that, the blitz dependency could be removed in the test case
mimicking the way blitz::toEnd is defined (in the latest blitz
version, because I find it much easier to understand). See it here. I
don't understand why that doesn't work, but declaring toEnd as
const solves the problem. (Why?)

I also added a test case where it seemed more appropriate
(scipy/weave/tests/test_blitz_tools.py), but it's probably not the
perfect place. The full test suite passes (except an unrelated test
in scipy.io ImportError: cannot import name BytesIO).

I'd appreciate any comments or anything I have overlooked. Thanks!

[pv]

Good sleuthing!

This is a minimal example:

#include <iostream>
int _end = 1;
int main()
{   
    std::cout << _end << " == 1" << std::endl;
    return 0;
}

The problem seems to be that _end is apparently some (undocumented) global symbol, perhaps in the C++ standard library, and happens to have a value 0 :)

The fix is indeed to declare _end as const, static, or static const. Your fix seems fine to me, just submit it as a pull request, maybe explaining the reason of the bug in the commit message. (const implies static in C++)

[Jorge-C]

Thank you Pauli! I only kept the static qualifier for _end, and decorated the test as slow. I hope the pull request is fine.

[pv]

I'd prefer const qualifier as then the compiler can optimize the constant away, rather than having to always look for the global variable.

[Jorge-C]

OK, I changed that.

[pv]

Fixed in 7164bb6