From d4dedb8bf51fc75777375fd8fe5192de20fd90c9 Mon Sep 17 00:00:00 2001
From: stepan <stepan.sindelar@oracle.com>
Date: Thu, 22 Dec 2016 20:08:11 +0100
Subject: [PATCH] Stat externals: ppois, qpois
---
.../oracle/truffle/r/library/stats/PPois.java | 37 +++++
.../oracle/truffle/r/library/stats/QPois.java | 129 ++++++++++++++++++
.../foreign/CallAndExternalFunctions.java | 6 +
.../test/library/stats/TestStatFunctions.java | 2 +-
4 files changed, 173 insertions(+), 1 deletion(-)
create mode 100644 com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/PPois.java
create mode 100644 com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/QPois.java
diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/PPois.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/PPois.java
new file mode 100644
index 0000000000..cb67510db3
--- /dev/null
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/PPois.java
@@ -0,0 +1,37 @@
+/*
+ * This material is distributed under the GNU General Public License
+ * Version 2. You may review the terms of this license at
+ * http://www.gnu.org/licenses/gpl-2.0.html
+ *
+ * Copyright (C) 1998 Ross Ihaka
+ * Copyright (c) 2000, The R Core Team
+ * Copyright (c) 2016, Oracle and/or its affiliates
+ *
+ * All rights reserved.
+ */
+package com.oracle.truffle.r.library.stats;
+
+import static com.oracle.truffle.r.library.stats.GammaFunctions.pgamma;
+
+import com.oracle.truffle.r.library.stats.StatsFunctions.Function2_2;
+
+public class PPois implements Function2_2 {
+ @Override
+ public double evaluate(double x, double lambda, boolean lowerTail, boolean logP) {
+ if (Double.isNaN(x) || Double.isNaN(lambda) || lambda < 0.) {
+ return RMath.mlError();
+ }
+ if (x < 0) {
+ return DPQ.rdt0(lowerTail, logP);
+ }
+ if (lambda == 0.) {
+ return DPQ.rdt1(lowerTail, logP);
+ }
+ if (!Double.isFinite(x)) {
+ return DPQ.rdt1(lowerTail, logP);
+ }
+ x = Math.floor(x + 1e-7);
+
+ return pgamma(lambda, x + 1, 1., !lowerTail, logP);
+ }
+}
diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/QPois.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/QPois.java
new file mode 100644
index 0000000000..9803bad38d
--- /dev/null
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/QPois.java
@@ -0,0 +1,129 @@
+/*
+ * This material is distributed under the GNU General Public License
+ * Version 2. You may review the terms of this license at
+ * http://www.gnu.org/licenses/gpl-2.0.html
+ *
+ * Copyright (C) 1998 Ross Ihaka
+ * Copyright (c) 2000-2014, The R Core Team
+ * Copyright (c) 2016, Oracle and/or its affiliates
+ *
+ * All rights reserved.
+ */
+package com.oracle.truffle.r.library.stats;
+
+import static com.oracle.truffle.r.library.stats.MathConstants.DBL_EPSILON;
+
+import com.oracle.truffle.r.library.stats.DPQ.EarlyReturn;
+import com.oracle.truffle.r.library.stats.StatsFunctions.Function2_2;
+
+public class QPois implements Function2_2 {
+ private final Qnorm qnorm = new Qnorm();
+ private final PPois ppois = new PPois();
+
+ @Override
+ public double evaluate(double p, double lambda, boolean lowerTail, boolean logP) {
+ if (Double.isNaN(p) || Double.isNaN(lambda)) {
+ return p + lambda;
+ }
+ if (!Double.isFinite(lambda)) {
+ return RMath.mlError();
+ }
+ if (lambda < 0) {
+ return RMath.mlError();
+ }
+ if (lambda == 0) {
+ return 0;
+ }
+
+ try {
+ DPQ.rqp01boundaries(p, 0, Double.POSITIVE_INFINITY, lowerTail, logP);
+ } catch (EarlyReturn e) {
+ return e.result;
+ }
+
+ double mu = lambda;
+ double sigma = Math.sqrt(lambda);
+ /* gamma = sigma; PR#8058 should be kurtosis which is mu^-0.5 */
+ double gamma = 1.0 / sigma;
+
+ /*
+ * Note : "same" code in qpois.c, qbinom.c, qnbinom.c -- FIXME: This is far from optimal
+ * [cancellation for p ~= 1, etc]:
+ */
+ if (!lowerTail || logP) {
+ p = DPQ.rdtqiv(p, lowerTail, logP); /* need check again (cancellation!): */
+ if (p == 0.) {
+ return 0;
+ }
+ if (p == 1.) {
+ return Double.POSITIVE_INFINITY;
+ }
+ }
+ /* temporary hack --- FIXME --- */
+ if (p + 1.01 * DBL_EPSILON >= 1.) {
+ return Double.POSITIVE_INFINITY;
+ }
+
+ /* y := approx.value (Cornish-Fisher expansion) : */
+ double z = qnorm.evaluate(p, 0., 1., /* lower_tail */true, /* log_p */false);
+ // #ifdef HAVE_NEARBYINT
+ // y = nearbyint(mu + sigma * (z + gamma * (z*z - 1) / 6));
+ // #else
+ double y = Math.round(mu + sigma * (z + gamma * (z * z - 1) / 6));
+
+ z = ppois.evaluate(y, lambda, /* lower_tail */true, /* log_p */false);
+
+ /* fuzz to ensure left continuity; 1 - 1e-7 may lose too much : */
+ p *= 1 - 64 * DBL_EPSILON;
+
+ /* If the mean is not too large a simple search is OK */
+ if (lambda < 1e5) {
+ return search(y, z, p, lambda, 1).y;
+ } else {
+ /* Otherwise be a bit cleverer in the search */
+ double incr = Math.floor(y * 0.001);
+ double oldincr;
+ do {
+ oldincr = incr;
+ SearchResult searchResult = search(y, z, p, lambda, incr);
+ y = searchResult.y;
+ z = searchResult.z;
+ incr = RMath.fmax2(1, Math.floor(incr / 100));
+ } while (oldincr > 1 && incr > lambda * 1e-15);
+ return y;
+ }
+ }
+
+ private SearchResult search(double yIn, double zIn, double p, double lambda, double incr) {
+ if (zIn >= p) {
+ /* search to the left */
+ double y = yIn;
+ for (;;) {
+ double z = zIn;
+ if (y == 0 || (z = ppois.evaluate(y - incr, lambda, /* l._t. */true, /* log_p */false)) < p) {
+ return new SearchResult(y, z);
+ }
+ y = RMath.fmax2(0, y - incr);
+ }
+ } else { /* search to the right */
+ double y = yIn;
+ for (;;) {
+ y = y + incr;
+ double z = zIn;
+ if ((z = ppois.evaluate(y, lambda, /* l._t. */true, /* log_p */false)) >= p) {
+ return new SearchResult(y, z);
+ }
+ }
+ }
+ }
+
+ private static final class SearchResult {
+ final double y;
+ final double z;
+
+ SearchResult(double y, double z) {
+ this.y = y;
+ this.z = z;
+ }
+ }
+}
diff --git a/com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/foreign/CallAndExternalFunctions.java b/com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/foreign/CallAndExternalFunctions.java
index cb63d6058c..0c2e5d5d50 100644
--- a/com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/foreign/CallAndExternalFunctions.java
+++ b/com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/foreign/CallAndExternalFunctions.java
@@ -71,11 +71,13 @@ import com.oracle.truffle.r.library.stats.LogNormal.QLNorm;
import com.oracle.truffle.r.library.stats.Logis;
import com.oracle.truffle.r.library.stats.Logis.DLogis;
import com.oracle.truffle.r.library.stats.Logis.RLogis;
+import com.oracle.truffle.r.library.stats.PPois;
import com.oracle.truffle.r.library.stats.Pbeta;
import com.oracle.truffle.r.library.stats.Pbinom;
import com.oracle.truffle.r.library.stats.Pf;
import com.oracle.truffle.r.library.stats.Pnorm;
import com.oracle.truffle.r.library.stats.Pt;
+import com.oracle.truffle.r.library.stats.QPois;
import com.oracle.truffle.r.library.stats.Qbinom;
import com.oracle.truffle.r.library.stats.Qnorm;
import com.oracle.truffle.r.library.stats.Qt;
@@ -309,6 +311,10 @@ public class CallAndExternalFunctions {
return StatsFunctionsFactory.Function3_1NodeGen.create(new DUnif());
case "qunif":
return StatsFunctionsFactory.Function3_2NodeGen.create(new QUnif());
+ case "ppois":
+ return StatsFunctionsFactory.Function2_2NodeGen.create(new PPois());
+ case "qpois":
+ return StatsFunctionsFactory.Function2_2NodeGen.create(new QPois());
case "rbinom":
return RandFunction2Node.createInt(new Rbinom());
case "pbinom":
diff --git a/com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/library/stats/TestStatFunctions.java b/com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/library/stats/TestStatFunctions.java
index 29626c0df3..1e355ca8a4 100644
--- a/com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/library/stats/TestStatFunctions.java
+++ b/com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/library/stats/TestStatFunctions.java
@@ -61,7 +61,7 @@ public class TestStatFunctions extends TestBase {
assertEval(Output.IgnoreWhitespace, template("set.seed(1); %0(%1)", FUNCTION2_1_NAMES, FUNCTION2_1_PARAMS));
}
- private static final String[] FUNCTION2_2_NAMES = {"pchisq", "pexp", "qexp", "qgeom", "pgeom", "qt", "pt"};
+ private static final String[] FUNCTION2_2_NAMES = {"pchisq", "pexp", "qexp", "qgeom", "pgeom", "qt", "pt", "qpois", "ppois"};
private static final String[] FUNCTION2_2_PARAMS = {
"0, 10",
"c(-1, 0, 0.2, 2), rep(c(-1, 0, 0.1, 0.9, 3), 4)",
--
GitLab