Breaking: Internal function log1p now carries a
cubic Taylor term. The prior quadratic form left a truncation error of
order \frac{\beta^3}{3}. For large \alpha with
very small \beta this reached the sixth significant figure of
ordinary probabilities. Results for \beta > 10^{-4} are
bitwise-unchanged.
Breaking: Non-finite parameters now return NaN with
a warning across d/p/q/rdelap instead of clamp-laundered 0s and 1s.
Breaking: ddelap, pdelap, qdelap, and
rdelap with singleton parameters now first attempt to build their
PMF and CDF tables using a three-term recurrence
relationship derived from the Delaporte distribution's probability
generating function. This takes O(K) arithmetic operations instead
of O(K^2) log_gamma/exp evaluations. Masses whose
leading term underflows—roughly \lambda + \alpha\log(1+\beta) >
745—are handled by internal rescaling, allowing accuracy to be retained
for very large \lambda. This is done natively in linear- or
log-space as per the call. This may shift results relative to the old
“always sum” algorithm at the level of accumulated rounding: one or
two ulps for short tables, growing with table length to roughly
the 13th significant digit for very long ones. Both parameters extreme
enough to risk overflow in the recurrence fall back to the summation
algorithm and vector-valued parameters continue to use the per-element
path.
Breaking: ddelap and pdelap now treat
arguments below the support the same way base R does. ddelap
returns 0 (-Inf with log = TRUE) and pdelap
returns a lower-tail CDF of 0 and a survival of 1.
Previously these inputs returned NaN with a “NaNs produced”
warning.
Breaking: ddelap's internal summand formula
subtracted \lambda from its running sum before a large canceling
pair, \log\Gamma(\alpha + i) - \log\Gamma(\alpha) (exactly 0 at i = 0), had resolved. When
\alpha is very small, \log\Gamma(\alpha) is
itself large, and \lambda, when smaller than that intermediate
value's rounding unit, was silently discarded before the cancellation
completed. This was invisible with log = FALSE, since exp()
of the resulting error and exp(0) = 1 are both indistinguishable in
linear-space double precision at the affected magnitudes. Parenthesizing
the canceling pair so it resolves first fixes both the log-space and
linear-space summand formulas. Ordinary (non-degenerate) parameters are
bitwise-unchanged.
Breaking: ddelap with log = TRUE now
accumulates the PMF in log space using a streaming log-sum-exp instead of
taking the log of the linear-space result. This is technically breaking
as deep-tail log-probabilities which previously returned -Inf
because the linear PMF underflowed below the smallest representable double
are now accurate to near full precision. For example,
ddelap(2000, 1, 1, 1, log = TRUE) returns approximately -1385.99
instead of -Inf. Hoisting the variables for both log and normal
space routines makes them slightly faster than before.
Breaking: pdelap with lower.tail = FALSE now
computes the upper tail by direct summation of the PMF when doing so is
needed to preserve accuracy, instead of 1 - \textrm{CDF}. Survival
probabilities near or below machine epsilon previously returned only
floating-point rounding noise. For example,
pdelap(60, 1, 1, 1, lower.tail = FALSE) returned roughly 3.3e-16
when the true value is roughly 1.18e-18. The algorithms are now accurate
to near full double precision. This is technically breaking, but users
should see no change outside of the deep tail which was likely not of
interest to most. This summation is used only once the survival
probability is small enough that the ordinary complement
1 - \textrm{CDF} would lose meaningful precision: below
\sqrt{\epsilon} \approx 1.49 \times 10^{-8}.
Since the summation's remainder-bound safeguard gives it an
O(\beta^2) worst-case cost, it is reserved for where it
is genuinely needed.
Breaking: pdelap with log.p = TRUE previously
computed the linear-space CDF/survival first and took log() of the
result, so any value that underflowed to exactly 0 in the deep tail
silently returned -Inf even though the true log-probability is
finite. log.p = TRUE is now computed directly in log space
throughout, matching the log-space guarantee
ddelap(..., log = TRUE) already provides, and is accurate to near
full double precision even where the linear-space value underflows.
pdelap derives log-probabilities in the saturated bands—log-CDF
values within about sqrt(.Machine$double.eps) of 0, and their
log-survival mirror—from the opposite, relatively accurate side via
log1p(-exp(.)) rather than reporting the accumulation's own
rounding noise. On the log scale those values now agree with an
independent dnbinom\ast
dpois convolution oracle to about
1e-13 relative. In linear space, the original algorithm was already
accurate to an ulp, so linear-scale results are unaffected, as
is everything outside the bands. All three computational
routes—recurrence table, degenerate-ratio elemental fill, and
vector-parameter recycling—apply the same repair. A consequence is
that qdelap/pdelap round trips with log.p = TRUE are
now exact in both tails wherever the input is representable.
Breaking: qdelap for probabilities within
accumulated rounding of 1 combined with highly overdispersed
parameters previously ground its lookup-table build indefinitely,
as the cumulative sum plateaus below such targets no matter how far it
extends. These targets are now resolved from the accurate upper-tail
survival table, returning the true finite quantile instantly.
For example, qdelap(1 - 1e-12, 5, 5, 7000) now returns 7630.
Ordinary probabilities are unaffected.
Breaking: qdelap with lower.tail = FALSE
and/or log.p = TRUE now searches in the native space of the
supplied probabilities instead of transforming them to linear
lower-tail space first. The old 1 - p flip collapsed survival
probabilities below machine epsilon to exactly 1, returning
Inf where qpois resolves a finite
quantile, and the old exp() of log-probabilities collapsed deep
log-space targets the same way; near-1 targets were inflated by an
ulp, stepping round-tripped quantiles up by one. Upper-tail
targets are now searched against a survival table built as pdelap,
and log-space targets against log-probability tables accumulated by
streaming log-sum-exp. The default (lower.tail = TRUE,
log.p = FALSE) path is bitwise-unchanged. Interior search targets
carry one-sided tolerances sized to their measured noise sources, in the
spirit of base R's 64 * EPS fuzz in ‘qpois.c’. The
vector-parameter path retains base R's fuzz verbatim.
Breaking: qdelap no longer hangs when passed extreme
parameter sets whose quantiles are computationally unreachable. The
original explicit-summation algorithm, now used when the recurrence table
would overflow, costs O(\textrm{answer}^2). Past
2^{15} summation steps qdelap now returns
NaN for the unresolved targets with the warning
“quantile too large to compute exactly for these parameter
values; NaN returned”, as a truncated cumulative sum carries no
information about a quantile beyond it. The parameter sets involved
generally put the true quantile beyond 2^{53}, where a count
is no longer exactly representable as a double in any case. Targets the
partial sum did resolve, including the p = 0 and p = 1
boundaries, still return normally.
Breaking: The approximate form of qdelap, which has
much less utility now (see documentation), will always generate at least
10,000 variates for pooling. It remains capped at no more than 10M
variates. Prior behavior for pathologically tiny means (~1e-12) could have
resulted in no variants being generated and an NA return.
Breaking: The approximate form of qdelap now
pre-allocates the result vector at full length first and then assigns the
NaN, Inf, 0, or proper value for each category. This ensures
that element i of the output always corresponds to element i
of the input. Previously, the return may have scrambled the "valid" inputs
if there were non-valid values in the input.
For x \geq 2^{63}, including +Inf, the functions no
longer convert x to an integer kind which it then overflows.
Instead, they return 0 or -Inf directly.
Zero-length parameters now return zero-length results in
p/d/qdelap per base R. This is not considered breaking as it used
to fail with a segfault from Fortran.
Zero-length parameters in rdelap now return NaNs per
base R. This is not considered breaking as it used to fail with a
segfault from Fortran.
qdelap and rdelap with a NaN scalar parameter
previously looped indefinitely while consuming memory. Now they return
NaN with a warning, matching pdelap.
Prevented qdelap from overwriting passed p vector
within Fortran, preserving the .Call contract.
Breaking: Both the exact and approximate versions of
qdelap now return NaN with a warning for probabilities
strictly greater than one, following base R convention. The function
returns Inf only for probabilities exactly equal to one.
Prior behavior was that all probabilities \ge 1 returned Inf.
Breaking: setDelapThreads no longer calls
omp_set_num_threads, and getDelapThreads no longer calls
omp_get_max_threads. The thread count is package-local and is
passed to each parallelized routine, whose num_threads clauses
already scope it per region. Process-global OpenMP state is never
modified, so Delaporte neither affects nor is affected by the
threading of other OpenMP-using packages. On load, the count is seeded
from the OpenMP runtime default, which respects OMP_NUM_THREADS and
OMP_THREAD_LIMIT.
As the singleton-parameter CDF build is no longer
quadratic, the maximum q for pdelap's fast path was
raised from 2^{14} to 2^{24}, and the interactive
confirmation prompt now guards only the genuinely slow routes:
vector-valued parameters with q at or above 2^{15}, or any
q at or above 2^{24}. When not interactive, pdelap
checks the new (and undocumented) Delaporte.maxq option in the
local environment and stops if the largest value exceeds the option or
warns if it does not. The default value is 2^{15}.
Changing qdelap to use the lookup table most of the time
exact = TRUE makes it orders of magnitude faster than
exact = FALSE for essentially all practical parameter ranges.
The only time where exact = FALSE may still be worth considering is
when the expected table length, K = \alpha\beta + \lambda +
10\sqrt{\alpha\beta(1+\beta) + \lambda}, exceeds roughly 1.5\times
10^{7}. Beyond that point, the exact path's memory use of ~17
bytes per unit of K begins to exceed the bounded ~330MB/2.5s cost of
the Monte Carlo approximation. The table size K depends on the
distribution's dispersion, not just its mean. A low mean with a high
overdispersion can require as large a table as a much larger,
near-Poisson mean. The exact = FALSE remains unavailable for
vector-valued parameters regardless of this threshold, as approximate
quantiles cannot be pooled across distinct parameter triplets.
qdelap and rdelap with singleton parameters now size
their quantile lookup table from a moment-based estimate and grow it
geometrically, replacing element-by-element growth whose copying was
O(K^2). The table is built with the identical routine and
accumulation order as pdelap, so the qdelap(pdelap(k)) == k
round trip is preserved bitwise. They now also locate quantiles in the
sorted CDF vector by a binary search instead of a linear scan.
MoMdelap now checks that at least 3 values are passed to it,
none of which are NA, otherwise it will error. This is not
considered breaking as it used to fail with a different error which
propagated from Fortran (missing value where TRUE/FALSE needed). Now it
never triggers a Fortran call.
Cleaned ‘Makevars’ from unneeded calls to BLAS and LAPACK libraries.
The guidance added to Writing R Extensions in r88406 is better adhered to by storing the package-specific state (the thread count and detected CPU count) in an environment inside the package namespace. The environment is no longer created in the user's workspace, and is no longer exported.
Refactored C and Fortran sources to be more manageable using header files, splitting initialization from active function calls, hardening with explicit "only" lists on the Fortran side, removing unnecessary includes, and using the modern bind(C) interface to call back-and-forth between C and Fortran instead of the older F77 idiom.
The manual now documents the package's accuracy limits at the
representation edges: the context dependence of lower.tail =
FALSE values at the level of their accumulated rounding, the
\pm 1 round-trip shoulder for survival targets within about
sqrt(.Machine$double.eps) of 1, the linear-space
convention (and consequent deep-tail saturation) of vector-parameter
qdelap, and the cost profile of deep upper-tail log-scale
quantiles at high dispersion.
Updated unit tests and documentation for the changes and fixed a number of spelling and grammar mistakes.
Breaking: The sOMPT_C function and entry point was
removed, including its R_RegisterCCallable registration, as were
the .onDetach thread reset and its Fortran backend sOMPT_f.
Breaking: DelaporteEnv is no longer exported.
Corrected improper link in ‘SECURITY.md’.
Other spelling errors corrected.
Changes related to Github workflows and requirements.
Minor tweaks to R and Fortran code. No change to calculations.
Corrected package metadata unit tests and switched sapply to
vapply for type safety.
Updated Github action YAMLs.
Updated ‘NEWS’ markup.
Switched “max cpu” holding variable location from
options to an environment variable per R commit 88406 by Brian
Ripley.
Converted some architecture-specific tests to use looser tolerances based on Tomas Kalibera's blog post.
Made checks of the lt and lg variables passed to
Fortran from C to be against c_int variables, which they should be.
Testing for NaN within Fortran now more efficient.
Other changes to C and Fortran code for efficiency and readability.
Added package metadata tests.
Set up new version of rhub checking.
Use native CRAN DOI.
qdelap now automatically calls the exact = TRUE
version when passed inappropriate parameters instead of just stopping with
an error.
Refactor ‘README.md’ into component MD files for Github purposes.
Capture and restore original thread count after running tests.
Efficiency tweaks to underlying Fortran code.
Converted log1p implementation to one based on its Taylor
sequence, instead of catastrophic substitution elimination, as some
compilers may optimize out that subtraction. Since any difference in
result should be at machine precision or lower, this is considered a
minor update.
Removed back keyword from minloc call in the Fortran
implementation of qdelap. Technically unneeded since its value is
the default FALSE. This allows for systems with an only partial
implementation of Fortran 2008 to compile (see
Issue 3 and thanks to
David Pascall).
Correct error message of setDelapThreads to reflect that
setting thread count to 0 is also an error. Functionality has not changed.
Removed unnecessary explicit check for NaN in qdelap
as it will be covered by check using anyNA.
Tweaked documentation for clarity.
Added unit test for new log1p function.
Added LTO_OPT and LTO_FC_OPT flags to ‘Makevars’.
Added getDelapThreads and setDelapThreads to
explicitly control OpenMP thread usage based on
WRE 1.2.1.1.
Added parallel as an import for detectCores.
Added SPDX header to more files.
Added more unit tests.
Changed Fortran suffix to .f90 for better compatibility with
the Intel Fortran compiler, per CRAN.
Tweaked minor linting and documentation issues.
Explicitly enumerated exports instead of using wide-ranging regex.
Removed -std=f2008 on the advice of Professor Ripley.
Fortran source code now reliant on Fortran 2008.
Fortran module ieee_arithmetic used.
Updated help documentation.
Updated Github actions.
Updated internal code to be more compliant with personal style guide.
Updated ‘CITATION’ to use bibentry per CRAN
request.
gamln replaced with Fortran 2008 log_gamma.
position replaced with Fortran minloc.
The set_nan and set_inf procedures were replaced with
ieee_arithmetic calls.
Checks for NaN now use ieee_arithmetic calls.
Updated package with information necessary to achieve OpenSSF silver status.
Converted unit testing framework to tinytest to reduce dependencies.
Augmented citation with DOI and CFF file.
Corrected help file.
Removed SIMD calls added in 8.0.0 due to Solaris.
Removed LTO mentions from Makevars added in 8.0.0 per Brian Ripley.
Pass more R-consistent error messages.
Made more robust to improper inputs such as negative, NA, or
NaN parameters or passed values.
R versions of qdelap and rdelap behave more like the
Fortran versions.
Exact code more in line with R defaults for d/p/q/r functions in
that NaN inputs return NaN outputs.
Handling of Inf and -Inf now more in line with R
defaults for d/p/q/r functions.
Parameter errors now check for < 0 instead of < EPS.
ddelap and pdelap now use 8-byte integers
(INT64) for their counters which allows calculation for values
up to 2^{63}.
ddelap will return 0 for all values > 2^{63}
under the assumption that the PMF is spread too thinly.
pdelap will not run for values > 2^{63}.
ddelap of Inf returns 0 and pdelap of
Inf returns 1.
pdelap gives a warning for finite values
> 2^{15} as, depending on the parameters, this may take
hours to finish (it's a double summation!). The user has to respond with
y to continue.
Updated help and tests.
Allow compiler to use OpenMP SIMD constructs.
Allow using LTO.
Switched continuous integration from Travis/AppVeyor to Github actions (thanks to Jim Hester for his help in fixing the coverage check issue).
Corrected CRAN requests.
Cleaned up ‘NEWS’.
Converted tests to be compliant with testthat: 3rd edition.
Cleaned up prose.
Added a citation to ‘DESCRIPTION’.
Added a test to increase code coverage.
‘CITATION’ file.
Explicit citation information to ‘README’.
Travis, AppVeyor, and CodeCov integration.
Added R_forceSymbols(dll, TRUE) to source per WRE
5.4.
Tests for package versioning.
Moved canonical repository to Github.
Use FFLAGS instead of FCFLAGS in Makevars.
Removed depsy badge from ‘README’ as site is now defunct.
Updated ‘NEWS’ to reflect recent changes.
Updated various file versions to reflect recent changes.
Makevars adjusted by CRAN to prevent issues with parallel make.
Changed Makevars flags and targets and forced dependence on R 3.6+ to comply with new CRAN policies on Fortran and OpenMP.
Preventing spurious floating point errors by implementing hard floor of 0 and hard ceiling of 1.
Updated comments in ‘delaporte.f95’ source code.
Changed rdelap testing to always use a fixed seed.
Documentation cleaned up.
Removed the cleanzeros function.
Removed unneeded Makevars target cleanup.
Removed GNU-specific Makevars target .NOTPARALLEL.
Added cleanzeros function that zaps to 0 any value less than
or equal to EPS (thanks Kipper Fletez-Brant).
Test for zeros (thanks Kipper Fletez-Brant).
Added items to ‘README’ to comply with CII best practices and added badge.
Added various badges to ‘README’.
Added copyright and license information inside R source code file.
MoMdelap now allows a skew type selection, similar to
that of the skewness function in the e1071 package
(thanks to William Cipolli for the suggestion). Skew defaults to former
selection so there should be no change to any values.
Tests added for all three skewness variants.
ORCID added to author description.
Copyright date corrected to reflect actual copyright and not date of licensure change.
There was an error in the calculation of beta, and thus all
three parameters, in MoMdelap (thanks to William Cipolli).
The approximate version of qdelap was not correct. It pooled
together ALL generated random variables. As it is impractical to generate
millions of random variates for each unique parameter triplet, the
functionality was changed so that only the “exact” version of
qdelap may be used if any of the parameters are a vector.
In the continued struggle to have the package pass all its tests on Solaris SPARC, all logicals needing to be passed to Fortran are now passed as integers and not booleans. While all other platforms had no problem using Fortran 2003's ISO_C_bindings, Solaris SPARC did.
Added note to ‘README’ about Solaris SPARC.
More tests including for singleton approximate qdelap.
Older implementation of qdelap and rdelap
approximations are now completely removed from code and documentation.
Solaris, being big-endian, appears unable to use the same bit-mask
transfer for NaN and Inf as any other R platform.
Therefore, the package has to revert to calling C to call the R magic
words which should recognize the platform.
Changed handling of non-integers to be in line with R convention
that they have 0 probability instead of being considered the next integer.
ddelap will return 0 for each non-integral value and will pass one
warning, unlike base R which passes one warning per non-integral value.
qdelap will just return the value for the last integer and no
warnings, similar to base R.
Use bit-representation of NaN and Inf to remove need
to call C for setting these values. This allows the log-gamma and
singleton d/p/qdelap functions to be elemental.
Changed registration of the C SEXP functions seen by the
package to use R_CallMethodDef as now preferred by CRAN.
Explicitly listed the five C SEXP functions seen by the
package as R_RegisterCCallable for potential downstream package
use.
Corrected some prose in documentation page.
Older ‘NEWS’ entries refactored for accuracy.
Explicitly cast n as an integer in loggamma routine. This
should suppress warning on OSX mavericks.
Older implementation of qdelap and rdelap
approximations is now defunct. Only the new faster approximation based on
section 3.9 of Karlis & Xekalaki (2005) remains. The next minor version
may remove the .Defunct calls and the references to the old version
in the documentation.
Removed colClasses = 'double' from tests. This should solve
Sun SPARC issues (It didn't).
Specified, clarified, and added many more unit tests.
Back end code ported from C++ & Rcpp to Fortran and C. In most cases this alone showed an increase in speed. There should be no noticeable difference in results. Random numbers generated with the same seeds should be the same as in previous packages.
OpenMP calls activated for systems which have it which allows parallel processing of certain loops resulting in further speed increases over and above those supplied by Fortran. OpenMP compatibility is approximately version 3.1; no SIMD instructions are used.
License changed from GPL2+/LGPL3+ to BSD-2 clause.
Includes various emergency fixes to address parallel make and other compilation issues in versions 4.0.0–4.0.2.
Did not build properly on some systems using a parallel make.
Did not build properly on Fedora using Clang.
Did not build properly on various platforms such as Fedora and Solaris.
The distributional functions now return NaN where the passed
alpha, beta, or lambda is close to 0 (less than
machine epsilon). Within vectors, valid triplets should still return
proper values. The approximate versions of qdelap and rdelap
do not handle 0 properly, and an error is thrown directing the user to the
exact version (thanks are due to Jonathan Cairns and Jin-Rui Xu for their
help in testing).
Faster versions of the approximate qdelap and rdelap
have been implemented.
Tests for near 0 handling on exact version.
The old approximation versions of qdelap and rdelap
can still be used by passing old = TRUE. This will be removed
eventually.
Reverted back to 2.2-3 as change to trapping zeros caused serious downstream issues. Will work on more elegant correction shortly (thanks to Mikhail Spivakov, Jonathan Cairns, and Jin-Rui Xu).
In the next release (3.0.0), passing FALSE to
exact in rdelap and qdelap will use the method
described in section 3.9 of Karlis & Xekalaki (2005) to generate
approximately Delaporte random variates much more quickly for large
n and large values of the parameters. However, when implemented,
this will break reproducible even for specifically set seeds. The new code
is commented out in the function.
Change non-exact qdelap mode to numeric instead of
integer to handle returned Inf.
Refactor ‘NEWS’ to more closely comply with keep a changelog suggestions.
Check for any parameters not strictly greater than 0 (thanks to Mikhail Spivakov and Jonathan Cairns).
Explicitly import code from package stats.
Updated code in unit tests to reflect change in
testthat package (move to expect_equal).
Minor tweaks to prose and grammar and expanded description.
Updated ‘DESCRIPTION’ file to comply with recent changes.
Trapped machine precision issues when calling right-tailed version of cumulative distribution function (thanks to Mikhail Spivakov and Jonathan Cairns).
Update documentation to reflect rounding procedure when non-integer passed as value.
Vectorized d/p/q/rdelap; alpha, beta, and
lambda can now be vectors.
Added more robust tests.
Minor optimizations for speed.
Corrected lower.tail and log.p which were not
implemented properly.
Adjusted unit test so it can be run more easily on 32 bit machines with limited memory.
Fixed spelling and other prose.
Added function to calculate method of moments estimates for the parameters. These estimates should also serve as good starting points for maximum likelihood estimation.
Added unit testing.
Minor cleanup after Rcpp update (remove ‘Makevar’ files and tweak documentation).
Updating for changes made to Rcpp 0.11.0. Previous versions may no longer compile properly and the current version may no longer work with older versions of Rcpp.
Change “Depends” to “LinkingTo” in line with CRAN philosophy on package dependencies.
Update documentation to give rough estimate as to when
exact = FALSE should be used (rarely).
Have R code be written in a consistent style.
Update ‘README’ now that CRAN uses it as well.
Minor grammar and typography tweaks to documentation.
Poisson-negative binomial approximation to quantile and random
variate functions. This greatly speeds up calculation for large
\alpha, \beta, or \lambda.
Converted ‘NEWS’ to ‘NEWS.Rd’.
Added repository and bug report URL to DESCRIPTION.
Help documentation corrections, comments, and tweaks.
Initial version.