[flopscope 0.10.0] ALL of `flopscope.stats` silently promotes float32 to float64, doubling the bill downstream

Category: cost model / API consistency
Version: flopscope 0.10.0 (0.10.0+np2.2.6), whestbench 0.14.0
Relevant to: the dtype-aware billing introduced in v0.10.0 (t/18125)

Summary

Every flopscope.stats.norm function returns float64 regardless of input dtype,
while the ordinary elementwise ufuncs correctly preserve float32. Under the new
dtype-aware cost model this is not a cosmetic inconsistency: a single stats call
silently moves the entire downstream computation into the 2x lane.

This directly undercuts the guidance in the v0.10.0 announcement - “Casting your hot
path to float32 can now produce a real saving”
- because any hot path that touches
stats.norm.* is quietly cast back to float64 at the first call.

Reproduction

import flopscope
from flopscope import numpy as fnp

with flopscope.budget(int(1e12)):
    x = fnp.ones((8, 8), dtype=fnp.float32) * 0.5
    print(x.dtype)                              # float32

    print(flopscope.stats.norm.ppf(x).dtype)    # float64   <-- promotes
    print(flopscope.stats.norm.cdf(x).dtype)    # float64   <-- promotes
    print(flopscope.stats.norm.pdf(x).dtype)    # float64   <-- promotes

    print(fnp.sqrt(x).dtype)                    # float32   correct
    print(fnp.exp(x).dtype)                     # float32   correct
    print(fnp.log(x).dtype)                     # float32   correct
    print(fnp.maximum(x, 0.0).dtype)            # float32   correct

Observed on flopscope 0.10.0+np2.2.6, numpy 2.2.6 backend.

Why it matters under the v0.10.0 cost model

The promotion is invisible at the call site but propagates: the promoted array becomes
an operand of the next matmul, whose output is float64, and so on for the rest of the
network. There is no warning and nothing in the budget breakdown attributes the extra
cost to the stats call.

Measured on our own estimator, a 256-wide 32-deep MLP, using
flopscope.budget(...).op_log on the real shapes the estimator uses:

FLOPs/MAC share of charged compute
hot matmul, float64 3.991 98.0%
hot matmul, float32 1.996 97.2%

A single stats.norm.ppf call generating the sample matrix was enough to put all 32
hot matmuls
in float64. Fixing only that one call site - .astype(float32) on the
result - changed the charged cost of the whole estimator by a factor of ~2 with no
change in the numerical result to within 0.0015% MSE.

Note the fix is not a one-liner in general: we also had to cast the weight matrices and
an analytic per-neuron term, because any single float64 operand re-promotes the
product. That is the shape of the problem - one promoted array anywhere in a chain
silently reprices everything after it.

Suggested fix

Have flopscope.stats.* follow the same dtype rules as the ufuncs, i.e. preserve the
input dtype (or at minimum honour NumPy’s standard promotion rules rather than
unconditionally producing float64).

If preserving the dtype is not desirable for accuracy reasons in the underlying
implementation, then a) documenting it prominently next to the v0.10.0 dtype guidance,
and b) emitting a FlopscopeWarning on promotion inside a budget context, would both
remove the silent-repricing failure mode.

Scope: the whole flopscope.stats surface, not just norm

Swept all 8 distributions x {ppf, cdf, pdf} on 0.10.0 with a float32 input:

distribution ppf cdf pdf
cauchy float64 float64 float64
expon float64 float64 float64
laplace float64 float64 float64
logistic float64 float64 float64
norm float64 float64 float64
uniform float64 float64 float64
lognorm (needs shape args - not exercised)
truncnorm (needs shape args - not exercised)

18 of 18 callable functions promote. Not one preserves the input dtype. The
remaining 6 require distribution shape parameters and were not exercised, but there is
no reason to expect them to differ.

So this is a property of the stats layer as a whole rather than a norm oversight,
which is why we think it is worth fixing centrally rather than per-function.