ff(): Simpson integration weights off by n/(n-1) and wrong even-n pattern

[fabian-s]

ff() (and sff()) build their integration operator L with
compute_integration_weights(., "simpson") (R/pffr-utilities.R), which is
not composite Simpson's rule:

• the [1, 4, 2, ..., 4, 1] pattern is scaled by (b - a) / (3 * n) instead
  of (b - a) / (3 * (n - 1)), i.e. it divides by the number of grid
  points rather than the number of intervals;
• for even n the alternation ends in 2 before the closing 1
  ([1, 4, 2, 4, 2, 1] at n = 6), a pattern composite Simpson does not
  admit.

Consequently the weights do not sum to the length of the integration domain.
A constant on [0, 1] integrates to

n	30	31	60	61	93
sum(w)	0.9556	0.9677	0.9778	0.9836	0.9892

instead of 1. n = 93 is the CCA grid of the bundled DTI data.

Since the design matrix carries the factor, the estimated coefficient surface
beta(s, t) comes out rescaled by the reciprocal — inflated by 4.6% at
n = 30 and by 1.1% at n = 93. For odd n it is exactly a factor
n / (n - 1). Simulation studies that use the same weights to generate and
to fit the data are unaffected (the factor cancels); real-data fits are not.

The same expression is inlined in af(), af_old(), lf(), lf_old() and
lf.vd(), so pfr() terms are affected as well.

Fixed in branch fix-ff-weights: composite Simpson with
h = (b - a) / (n - 1), using Simpson's 3/8 rule on the last three intervals
for even n, so a constant integrates to exactly b - a and cubics are exact
for every n >= 3. The old weights stay reachable as
integration = "simpson_legacy" (documented as deprecated) for reproducing
results from earlier versions. Regression tests in
tests/testthat/test-pffr-integration-weights.R.

The "trapezoidal" weights were checked and are correct. The "riemann"
weights (used whenever limits = is given) sum to b - a plus one mean grid
spacing, which is inherent to a first-order rule giving all n points
positive weight; left unchanged but now documented.