scipy: BUG: `special.nctdtr` broken in main

Describe your issue.

special.nctdtr does not relturn correct values in main. Likely cause is incorrect wrapping in #19560.

Additional small bug introduced is that the cythonized code is also returning the incorrect types, e.g., what should be a python float (scalar) is being returned as np.float64 (numpy float64 scalar).

Reproducing Code Example

froms scipy import special
special.nctdtr(58, -3.8, 1.67)


return `nan`, when is should return something on the order of `0.9999999675827712`

Error message

No error, just a wrong result bug.

SciPy/NumPy/Python version and system information

1.13.0.dev0+1181.2353b15 2.0.0.dev0+git20240113.d2f60ff sys.version_info(major=3, minor=12, micro=0, releaselevel='final', serial=0)
Build Dependencies:
  blas:
    detection method: pkgconfig
    found: true
    include directory: /usr/local/include
    lib directory: /usr/local/lib
    name: openblas
    openblas configuration: USE_64BITINT= DYNAMIC_ARCH=1 DYNAMIC_OLDER= NO_CBLAS=
      NO_LAPACK= NO_LAPACKE= NO_AFFINITY=1 USE_OPENMP= HASWELL MAX_THREADS=2
    pc file directory: /usr/local/lib/pkgconfig
    version: 0.3.21.dev
  lapack:
    detection method: pkgconfig
    found: true
    include directory: /usr/local/include
    lib directory: /usr/local/lib
    name: openblas
    openblas configuration: USE_64BITINT= DYNAMIC_ARCH=1 DYNAMIC_OLDER= NO_CBLAS=
      NO_LAPACK= NO_LAPACKE= NO_AFFINITY=1 USE_OPENMP= HASWELL MAX_THREADS=2
    pc file directory: /usr/local/lib/pkgconfig
    version: 0.3.21.dev
  pybind11:
    detection method: config-tool
    include directory: unknown
    name: pybind11
    version: 2.11.1
Compilers:
  c:
    commands: cc
    linker: ld.bfd
    name: gcc
    version: 10.2.1
  c++:
    commands: c++
    linker: ld.bfd
    name: gcc
    version: 10.2.1
  cython:
    commands: cython
    linker: cython
    name: cython
    version: 3.0.8
  fortran:
    commands: gfortran
    linker: ld.bfd
    name: gcc
    version: 10.2.1
  pythran:
    include directory: ../../opt/python/cp312-cp312/lib/python3.12/site-packages/pythran
    version: 0.15.0
Machine Information:
  build:
    cpu: x86_64
    endian: little
    family: x86_64
    system: linux
  cross-compiled: false
  host:
    cpu: x86_64
    endian: little
    family: x86_64
    system: linux
Python Information:
  path: /opt/python/cp312-cp312/bin/python
  version: '3.12'

About this issue

Original URL
State: closed
Created 5 months ago
Comments: 31 (31 by maintainers)

Most upvoted comments

That’s hardly possible to achieve. These functions call each other and you have to find the right argument to traverse everything to hit the right branch. A bit of educated brute force would clean up the obvious ones.

It isn’t that hard. I have written two projects with a total of around 80k LOC that have 99.5%+ coverage.

I suspect that the existing test probably hit most of the paths, and so only need coverage for the rest. One big advantage to using Cython over fortran is improved access to coverage, and a better path to improving coverage. Unfortunately, it looks like SciPy doesn’t do coverage in CI anymore.

bashtage on Jan 18, 2024

All of the cdftnc_whichX for x > 1 have cumf in place of cumtnc

https://github.com/scipy/scipy/blob/2a4d6e978baa015935e009cd00537d2ddeee2531/scipy/special/_cdflib.pyx#L3612

steppi on Jan 17, 2024

@ilayn, everywhere that has twoi *= xi below, should actually have twoi = 2.*xi. A consequence of the Fortran code using two for 2.0d0 in twoi = two*xi.

https://github.com/scipy/scipy/blob/2a4d6e978baa015935e009cd00537d2ddeee2531/scipy/special/_cdflib.pyx#L4614-L4638

steppi on Jan 17, 2024

special.nctdtrinc is also broken and silently returns wrong values.

from scipy import special

df, t_stat = (98, 1.5)
alpha = 1/20
alpha_half = alpha / 2
nc_median = special.nctdtrinc(df, 0.5, t_stat)
ci = special.nctdtrinc(df, [1 - alpha_half, alpha_half], t_stat)
print(nc_median)
print(ci)

On 1.11.4

1.4961639036840118
[-0.47493335  3.46737272]

On nightly (1.13.0.dev0+1181.2353b15)

1.793278926568212
[9.24724715e+01 9.80628003e-03]

bashtage on Jan 17, 2024

@ilayn

problem is here

https://github.com/scipy/scipy/blob/736216d5ea2645cb76f4468d02addc2d18acd5c0/scipy/special/_cdflib.pyx#L3571-L3572

it should be

    if not (-1.e6 <= pnonc <= 1.e6):
        return (0., 0., -3, 1.e6 if pnonc > -1e6 else -1.e6)

you should probably check other bound conditions.

steppi on Jan 17, 2024