What did you do this week?

1. Look at the issue Is the r-value outputted by scipy.stats.linregress always the Pearson correlation coefficient?
2. WIP: ENH: improve sort_vertices_of_regions via Pythran and made it more readable
   • Tyler said test_spherical_voronoi may test inplace sort, and it is not recommended to remove a test. In this way, we’ll never pass the test.
   • For the type error, I can’t reproduce it on my computer. Is it similar to the issue BUG: RBFInterpolator fails when calling it with a slice of a (1, n) array? I encountered similar `reshaped` issues before, and found that often the type is the problem while `reshaped` is not. Once I support that type, I’ll not get the error. But in the case there they do support that type.

   
                       TypeError: Invalid call to pythranized function `sort_vertices_of_regions(int32[:, :], int32 list list)'
                       Candidates are:
                           - sort_vertices_of_regions(int64[:,:], int64 list list)
                           - sort_vertices_of_regions(int32[:,:], int32 list list)
                           - sort_vertices_of_regions(int32[:,:], int64 list list)
                           - sort_vertices_of_regions(int[:,:], int list list)                    
                   

3. Last week we concluded _spectral.pyx and _sosfilt.pyx are easy to be improved via Pythran, but later I found that _spectral.pyx already has a version in Pythran. For_sosfilt.pyx, I improved _sosfilt_float and leave _sosfilt_object in Cython. The performance for _sosfilt_float looks similar comparing Cython and Pythran. So I'm not sure whether I need to make a PR for it
4. ENH: improve siegelslopes via pythran , 10x faster. If needed, I can also improve scipy/stats/_stats_mstats_common.py ’s linregress, theilslopes and put them with siegelslopes in the same file. But other two functions do not have obvious loops so here I only improve siegelslopes.
5. ENH: improve cspline1d, qspline1d, and relative funcs via Pythran ,10x faster.
   • Segment fault on Azure pipelines. Because of calling itself in the function?
   • A lot of signatures. Any more concise way?
   • Actually, for those functions which have lots of signatures and also cause current segment faults - cspline1d_eval and qspline1d_eval , they don’t have many loops. I improved them because they are used to evaluate cspline1d and qspline1d , putting them in one file may look better. We can also leave them in the original file so that we won’t get above a.& b. problems

What is coming up next?

1. Keep working on ENH: improve cspline1d, qspline1d, and relative funcs via Pythran
2. Find more potential algorithms and improve them
3. Make a PR for _sosfilt_float and comment on it
4. keepdimsfeature support in Pythran

Did you get stuck anywhere?

I once said that np.expand_dims() does not support dim as keyword, I was wrong because the key is axis, but I still got the following error. However, np.expand_dims(x, 1) will work.

    (scipy-dev) charlotte@CHARLOTLIU-MB0 stats % pythran siegelslopes_pythran.py
    CRITICAL: I am in trouble. Your input file does not seem to match Pythran's constraints...
    siegelslopes_pythran.py:19:13 error: function uses an unknown (or unsupported) keyword argument `axis`
    ----
        deltax = np.expand_dims(x, axis=1) - x
                 ^~~~ (o_0)
    ----
     
    

What did you do this week?

• Added unit test for BUG: fix stats.binned_statistic_dd issue with values close to bin edge
• Added benchmarks for somersd BENCH: add benchmark for somersd
• Added tests in Pythran: Import test cases from scipy
• Wrote the first evaluations, will submit it later.
• WIP: ENH: improve sort_vertices_of_regions via Pythran and made it more readable However, I got some weird type error, and failed test_spherical_voronoi and test_region_types. Tyler suggested it may be not a good case for Pythran, and I'm still trying to find out why there are such errors.

What is coming up next?

• Submit the first evaluations
• Continue working on sort_vertices_of_regions(), try to fix the failures
• Look into and maybe improve some of the following algorithms: _spectral.pyx and _sosfilt.pyx

Did you get stuck anywhere?

The WIP PR mentioned above: WIP: ENH: improve sort_vertices_of_regions via Pythran and made it more readable . It fails two tests: test_spherical_voronoi and test_region_types.

What did you do this week?

First came to the old problem, bus error. It turns out that it is specific to Mac. We still don't know the cause of the problem yet.( bus error on Mac but works fine on Linux for _count_paths_outside_method pythran version)

Last week I said that the benchmark result is different from my timeit result. It is actually my mistake: I forgot to modify setup.py. After setting up correctly, the problem was fixed.

Also, for the algorithm binned_statistic_dd I was improving since last week, I have made a PR for it. At first, I improved the whole if-elif block and the benchmark shows it can make count, sum,mean 1.1x times faster, and make std, median, min, max 3x-30x faster . However, I found that Pythran can't support object type input so I failed some tests.To support object type, we need to keep the whole pure Python codes, and it will make the if-elif block duplicate and ugly. Since from the benchmark, there is not much improvement for count, sum,mean, I also tried to only improve std, median, min, max to make it look better and understandable So in the end, I only improved an small inner function but still get std, median, min, max 3x-30x faster, with no changes for count, sum,mean.(ENH: improved binned_statistic_dd via Pythran)

When I was improving binned_statistic_dd, there happened to be an open issue about float point comparision. I looked into that and fixed it.( BUG: fix stats.binned_statistic_dd issue with values close to bin edge )

Last but not least, I tried to speedup _voronoi discussed last week, and the Pythran version is 3x faster than the Cython one!

What is coming up next?

• Refer to the original Python version rather than the CPython one, make the Pythran version _voronoi more readable. After that, make a PR.
• Add test for the binned_statistic_dd bug
• Add benchmarks for somersd and _tau_b
• Prepare for the first evaluations
• In Pythran, import some scipy tests

Did you get stuck anywhere?

The bus error mentioned above, and build_docs failed on my PR recently.

What did you do this week?

As this project progressed, I came to realize that the diffculty of this project is to find good potential algorithms. This week I was searching potential algorithms in scipy.stats. At first, like the past, I used pytest --durations=50 to find the slowest 50 tests. I looked into all the functions but didn't find any suitable algorithms: some do not have obvious loops; some have loops but call another scipy method in the loop... Therefore, I began checking all the functions under scipy/stats one by one and finally found somersd and _tau_b are good candidates and I submitted a PR( gh-14308) to speedup them 4x~20x.

Besides, For the works that mentioned last week:

1. stats._moment: submitted an issue(pythran #1820) for keep_dims is not supported in np.mean()
2. stats._calc_binned_statistic: successfully improved this function and made the public function stats.binned_statistic_dd 3x-10x faster on min,max,std,median. I tried to improved the whole if-elif block but encountered some errors that I can't fix (see pythran #1819 )
3. stats._sum_abs_axis0: Thanks to Serge, the compliation errror due to variant type is fixed. I compiled and it is ~2x faster on _sum_abs_axis0 but do not have much gain on the public function onenormest. Moreover, actually there is no loop in _sum_abs_axis0 for input size smaller than 2**20(my bad!)
4. sparse.linalg.expm(_fragment_2_1): Last week I said this one is slower than the pure python version but it is not. My bad, actually with my input at that time, it will not get into _fragment_2_1 so I actually didn't test it. Only when the input size is larger than 9000 it will get into the function. Moreover, the input is csc matrix so it is not a suitable one for pythran.
5. the SciPy build error(see pythran #1815) mentioned last week: It is a really struggling problem. Serge and Ralf tried to help me fix that but it is still not working for now.

What is coming up next?

1. add benchmarks for somersd and _tau_b
2. consider merging the old benchmark PR?(gh-14228)
3. keep searching good potential algorithms to be improved.

Did you get stuck anywhere?

The problem I encountered when improving the whole if-elif block in stats.binned_statistic_dd (see pythran #1819 )

What did you do this week?

This week was quite struggling. My mentor Serge implemented supporting for `scipy.special.binom` quickly in Pythran and it shows a great improvement on the public function stats.ks_2samp( 2.62ms vs 88.2ms). However, when I was building scipy with the improved algroithm, we found that it would cause a loop problem. Serge made a PR to break the loop but in my computer it is still not working.

Then I turned to try other algorithms that I mentioned last week, however I encountered more problems:

1. stats._moment: keep_dims is not supported in np.mean()
2. stats._calc_binned_statistic: invalid pythran spec but I don't find anything wrong
3. stats.rankdata: invalid pythran spec
4. stats._sum_abs_axis0: compliation error
5. sparse.linalg.expm(_fragment_2_1): much slower than the pure python one, will keep investigating it.

What is coming up next?

• submit issue for keep_dims
• submit issue for error in stats._sum_abs_axis0
• continue improving stats._calc_binned_statistic
• Find out why sparse.linalg.expm pythran version is slower.

Did you get stuck anywhere?

Got stuck in many problems, as is written in What did you do this week section.