dwt_max_level seems wrong for haar wavelet

[baogorek]

I'm trying out the haar wavelet decomposition example in the "Wavelet for Kids" tutorial. The coefficients are matching when I specify the level directly:

import numpy as np
import pywt

print pywt.__version__
0.5.0.dev0+13b8fd7

x = np.array([1, 0, -3, 2, 1, 0, 1, 2])
haar = pywt.Wavelet('haar')
coeffs = pywt.wavedec(x, haar, mode = 'sym', level = 3)
coeffs
[array([ 1.41421356]),
 array([-1.41421356]),
 array([ 1., -1.]),
 array([ 0.70710678, -3.53553391,  0.70710678, -0.70710678])]

However, there's a sign of trouble ahead when I look at the calculated max level:

print pywt.dwt_max_level(x.shape[0], haar.dec_len)
35

And if I don't specify the level argument in wavedec, it does indeed use 35 levels, which seems wrong and doesn't match the formulas in the documentation:

coeffs = pywt.wavedec(x, haar, mode = 'sym')
coeffs
[array([ 92681.90002368]),
 array([ 0.]),
 array([ 0.]),
 array([ 0.]),
...

 array([ 0.]),
 array([ 0.]),
 array([-1.41421356]),
 array([ 1., -1.]),
 array([ 0.70710678, -3.53553391,  0.70710678, -0.70710678])]

[kwohlfahrt]

I cannot reproduce, or I've misunderstood the issue. The following looks correct to me?

>>> import numpy as np
>>> import pywt
>>> x = np.array([1, 0, -3, 2, 1, 0, 1, 2])
>>> haar = pywt.Wavelet('haar')
>>> pywt.dwt_max_level(x.shape[0], haar.dec_len)
3
>>> pywt.wavedec(x, haar, mode='sym')
[array([ 1.41421356]), array([-1.41421356]), array([ 1., -1.]), array([ 0.70710678, -3.53553391,  0.70710678, -0.70710678])]
>>> print pywt.__version__
0.5.0.dev0+13b8fd7

[baogorek]

Hi @kwohlfahrt, it looks fine on your end. I'm not sure why pywt.dwt_max_level(8, 2) is evaluating to 35 on my system. I've got a work-around, obviously, but for reference I'm on Windows 10 with sys.version

2.7.12 |Anaconda 4.1.1 (64-bit)| (default, Jun 29 2016, 11:07:13) [MSC v.1500 64 bit (AMD64)]

I used conda (with pip, if I remember correctly) to install the package from github.

[grlee77]

strange, we do have some windows testing via Appveyor.

@baogorek If you run the following test, does it also fail?:

from pywt.tests.test__pywt import test_dwt_max_level
test_dwt_max_level()

[baogorek]

The test did fail, but I had to run it by copying and pasting code from github (the test module wasn't included with my installation). Attaching a zipped copy of my pywt folder from Anaconda2/Lib/site-packages in case anything is weird there.

Also, I noticed that dwt_max_level depends on the c function size_log2 in _extentions/c/common.c, and that there are bit level operations there. I've recently been trying different c/c++ compilers (to get OpenBlas to work), and the one that my path is pointing to now is mingw64. I wonder if my choice of compiler is related to the error.

[kwohlfahrt]

Yes, that seems like a good guess. Could you output the pre-processed version of size_log2 on your system? With GCC that is the -E flag, so that should work for mingw.

[grlee77]

Appveyor uses MSVC, so I don't believe we have tested with mingw64. mingw64 would not have defined _MSVC_VER, so it looks like size_log2 would be calling __builtin_clzl . I suspect that size_t and unsigned long may not be equivalent in your case?

Two possible fixes for size_log2: 1.) using an #ifdef with __MINGW64___ to call a function matching the size_t type for that setup. Maybe it is __builtin_clz instead of __builtin_clzl? 2.) explicitly cast the argument to __builtin_clzl as unsigned long

[kwohlfahrt]

I think this should be testable by checking SIZE_MAX == ULONG/UINT/ULONGLONG_MAX and selecting the appropriate function that way. Will give that a go.

[baogorek]

@kwohlfahrt here's the preprocessor output. Briefly, dtw_max_level shows up as:

unsigned char dwt_max_level(size_t input_len, size_t filter_len){
    if(filter_len <= 1 || input_len < (filter_len-1))
        return 0;

    return size_log2(input_len/(filter_len-1));
}

[kwohlfahrt]

@baogorek - I meant the contents of size_log2, not it's call site. Could I have a look at that please? Also I've made a guess at a fix in #198, could you see if that works?

[baogorek]

@kwohlfahrt , I checked out your branch and ran gcc -E common.c. The part that contains size_log2 is:

unsigned char size_log2(size_t x){
# 40 "common.c"
    return sizeof(size_t) * 8 - leading_zeros - 1;
}

In #198, I added some install errors based on this function when using pip install.