Skip to content

scipy compat#2295

Draft
ricardoV94 wants to merge 5 commits into
pymc-devs:mainfrom
ricardoV94:scipy_compat
Draft

scipy compat#2295
ricardoV94 wants to merge 5 commits into
pymc-devs:mainfrom
ricardoV94:scipy_compat

Conversation

@ricardoV94

Copy link
Copy Markdown
Member
  • Remove vestigial epydoc support
  • Fix SymbolicOp hash and eq before the inner graph is built
  • Add ndtr, ndtri and log_ndtr for scipy compatibility
  • Fold log(erfc(x)) into a stable LogErfc Op

The __epydoc_asRoutine attribute was set so epydoc would render scalar_elemwise
ops as functions. Epydoc has not been used for many years and nothing reads the
attribute; the accompanying comment in basic.py is dated 2008.
A SymbolicOp constructed without input_types defers building its inner graph until
the first call, so it has no fgraph yet. __hash__ and __eq__ both read self.fgraph
unconditionally and raised AttributeError, which meant a deferred op could not be
put in a set, used as a dict key, or compared. Identify it by type and props until
it is built, keeping deferred and built ops distinct.
ndtr and ndtri are exact reparametrisations of erfc and erfcinv. log_ndtr needs its
own form: log(ndtr(x)) underflows to -inf below x=-37 and saturates to 0 above x=8,
and no existing rewrite rescued it, since local_log_erfc matches log(erfc(x)) while
ndtr(x) builds 0.5 * erfc(-x / sqrt(2)).

Ndtr and LogNdtr are SymbolicOps rather than plain helpers so log(Ndtr(x)) can be
folded into LogNdtr(x), which is bit-identical to log_ndtr and holds both tails.
Only the forward pass is affected: the gradient graph is built before rewrites run,
so log_ndtr remains the differentiable path.
local_log_erfc substituted a truncated asymptotic series for log(erfc(x)) past a
hand-scanned per-dtype threshold. The series is divergent, so it needed those
thresholds, and it was less accurate than what it replaced (3.5e-14 relative at the
seam) and slower (108ms vs 65ms per 1e6 elements under numba).

erfcx(x) = exp(x**2) * erfc(x) makes log(erfc(x)) == log(erfcx(x)) - x**2 an exact
identity, so no series or threshold is needed. erfcx is unusable as x -> 0, where it
tends to 1 and the log cancels, so the branches split at 1, where both are accurate
to ~1e-16.

Emitting an Op rather than an expression containing the matched log(erfc(x)) also
removes the self-reference that the local_log_erfc_applied tag existed to guard, so
both rewrites are now plain PatternNodeRewriters.

log(y * erfc(x)) is folded too, as log(y) + LogErfc(x). This is the form a log of a
scaled erfc arrives in, e.g. ndtr(x) == 0.5 * erfc(-x / sqrt(2)).
Each of these is either a rename of an existing Op or an exact reparametrisation of
one, so none costs precision against the scipy function it claims to match:

  log_expit(x)   -> -softplus(-x)
  i0e/i1e(x)     -> ive(0/1, x)
  k0/k1(x)       -> kv(0/1, x)
  k0e/k1e(x)     -> kve(0/1, x)
  jn             -> jv, which is what it aliases in scipy too
  chdtr(v,x)     -> gammainc(v/2, x/2)
  chdtrc(v,x)    -> gammaincc(v/2, x/2)
  chdtri(v,p)    -> 2 * gammainccinv(v/2, p)
  gdtr(a,b,x)    -> gammainc(b, a*x)
  gdtrc(a,b,x)   -> gammaincc(b, a*x)
  pdtr(k,m)      -> gammaincc(floor(k)+1, m)
  pdtrc(k,m)     -> gammainc(floor(k)+1, m)
  rgamma(x)      -> 1 / gamma(x)

The gammainc reparametrisations are bit-exact, as scipy computes them the same way.
pdtr truncates k, as scipy does. rgamma needs a switch because gamma returns nan
rather than inf at the non-positive integers, where rgamma is 0.

Aliases whose only spelling loses precision are deliberately left out: ndtri_exp
(-inf below y=-745), binom (gamma overflows past n=170), sindg and friends (exact at
multiples of 90 in scipy), and betaincc (the betainc(b,a,1-x) reflection is worse
than scipy's dedicated implementation).
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant