From f83c42a247e24b60004231aabc4abbe821b999e9 Mon Sep 17 00:00:00 2001 From: stepan <stepan.sindelar@oracle.com> Date: Fri, 9 Dec 2016 11:30:20 +0100 Subject: [PATCH] Implement pbeta and all {x}cauchy --- .../truffle/r/library/stats/Cauchy.java | 157 ++++++++++++++++++ .../oracle/truffle/r/library/stats/DPQ.java | 5 + .../r/library/stats/GammaFunctions.java | 4 +- .../oracle/truffle/r/library/stats/Pbeta.java | 11 +- .../truffle/r/library/stats/RCauchy.java | 31 ---- .../oracle/truffle/r/library/stats/RMath.java | 33 ++++ .../foreign/CallAndExternalFunctions.java | 14 +- .../r/runtime/ops/BinaryArithmetic.java | 3 +- .../test/library/stats/TestStatFunctions.java | 28 +++- mx.fastr/copyrights/overrides | 2 +- 10 files changed, 249 insertions(+), 39 deletions(-) create mode 100644 com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Cauchy.java delete mode 100644 com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RCauchy.java diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Cauchy.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Cauchy.java new file mode 100644 index 0000000000..3ab4aaf454 --- /dev/null +++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Cauchy.java @@ -0,0 +1,157 @@ +/* + * 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) 1998--2008, The R Core Team + * Copyright (c) 2016, 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.M_PI; +import static com.oracle.truffle.r.library.stats.TOMS708.fabs; + +import com.oracle.truffle.r.library.stats.DPQ.EarlyReturn; +import com.oracle.truffle.r.library.stats.RandGenerationFunctions.RandFunction2_Double; +import com.oracle.truffle.r.library.stats.RandGenerationFunctions.RandomNumberProvider; +import com.oracle.truffle.r.library.stats.StatsFunctions.Function3_1; +import com.oracle.truffle.r.library.stats.StatsFunctions.Function3_2; + +public final class Cauchy { + private Cauchy() { + // contains only static classes + } + + public static final class RCauchy implements RandFunction2_Double { + @Override + public double evaluate(double location, double scale, RandomNumberProvider rand) { + if (Double.isNaN(location) || !Double.isFinite(scale) || scale < 0) { + return RMath.mlError(); + } + if (scale == 0. || !Double.isFinite(location)) { + return location; + } else { + return location + scale * Math.tan(M_PI * rand.unifRand()); + } + } + } + + public static final class DCauchy implements Function3_1 { + @Override + public double evaluate(double x, double location, double scale, boolean giveLog) { + double y; + /* NaNs propagated correctly */ + if (Double.isNaN(x) || Double.isNaN(location) || Double.isNaN(scale)) { + return x + location + scale; + } + if (scale <= 0) { + return RMath.mlError(); + } + + y = (x - location) / scale; + return giveLog ? -Math.log(M_PI * scale * (1. + y * y)) : 1. / (M_PI * scale * (1. + y * y)); + } + } + + public static final class PCauchy implements Function3_2 { + @Override + public double evaluate(double x, double location, double scale, boolean lowerTail, boolean logP) { + if (Double.isNaN(x) || Double.isNaN(location) || Double.isNaN(scale)) { + return x + location + scale; + } + + if (scale <= 0) { + return RMath.mlError(); + } + + x = (x - location) / scale; + if (Double.isNaN(x)) { + return RMath.mlError(); + } + + if (!Double.isFinite(x)) { + if (x < 0) { + return DPQ.rdt0(lowerTail, logP); + } else { + return DPQ.rdt1(lowerTail, logP); + } + } + + if (!lowerTail) { + x = -x; + } + + /* + * for large x, the standard formula suffers from cancellation. This is from Morten + * Welinder thanks to Ian Smith's atan(1/x) : + */ + + // GnuR has #ifdef HAVE_ATANPI where it uses atanpi function, here we only implement the + // case when atanpi is not available for the moment + if (fabs(x) > 1) { + double y = Math.atan(1 / x) / M_PI; + return (x > 0) ? DPQ.rdclog(y, logP) : DPQ.rdval(-y, logP); + } else { + return DPQ.rdval(0.5 + Math.atan(x) / M_PI, logP); + } + } + } + + public static final class QCauchy implements Function3_2 { + @Override + public double evaluate(double p, double location, double scale, boolean lowerTail, boolean logP) { + if (Double.isNaN(p) || Double.isNaN(location) || Double.isNaN(scale)) { + return p + location + scale; + } + try { + DPQ.rqp01check(p, logP); + } catch (EarlyReturn e) { + return e.result; + } + if (scale <= 0 || !Double.isFinite(scale)) { + if (scale == 0) { + return location; + } + return RMath.mlError(); + } + + if (logP) { + if (p > -1) { + /* + * when ep := Math.exp(p), tan(pi*ep)= -tan(pi*(-ep))= -tan(pi*(-ep)+pi) = + * -tan(pi*(1-ep)) = = -tan(pi*(-Math.expm1(p)) for p ~ 0, Math.exp(p) ~ 1, + * tan(~0) may be better than tan(~pi). + */ + if (p == 0.) { + /* needed, since 1/tan(-0) = -Inf for some arch. */ + return location + (lowerTail ? scale : -scale) * Double.POSITIVE_INFINITY; + } + lowerTail = !lowerTail; + p = -Math.expm1(p); + } else { + p = Math.exp(p); + } + } else { + if (p > 0.5) { + if (p == 1.) { + return location + (lowerTail ? scale : -scale) * Double.POSITIVE_INFINITY; + } + p = 1 - p; + lowerTail = !lowerTail; + } + } + + if (p == 0.5) { + return location; + } // avoid 1/Inf below + if (p == 0.) { + return location + (lowerTail ? scale : -scale) * Double.NEGATIVE_INFINITY; + } // p = 1. is handled above + return location + (lowerTail ? -scale : scale) / RMath.tanpi(p); + /* -1/tan(pi * p) = -cot(pi * p) = tan(pi * (p - 1/2)) */ + } + } +} diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/DPQ.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/DPQ.java index 1ac88316bc..4925be2e76 100644 --- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/DPQ.java +++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/DPQ.java @@ -107,6 +107,11 @@ public final class DPQ { return lowerTail ? rdlexp(p, logP) : rdlog(p, logP); } + // R_D_Clog(p) (log_p ? log1p(-(p)) : (0.5 - (p) + 0.5)) /* [log](1-p) */ + public static double rdclog(double p, boolean logP) { + return logP ? RMath.log1p(-(p)) : (0.5 - (p) + 0.5); + } + // R_DT_qIv public static double rdtqiv(double p, boolean lowerTail, boolean logP) { return logP ? lowerTail ? Math.exp(p) : -Math.expm1(p) : rdlval(p, lowerTail); diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/GammaFunctions.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/GammaFunctions.java index a00d9241c0..1638c38e2c 100644 --- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/GammaFunctions.java +++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/GammaFunctions.java @@ -371,7 +371,7 @@ public abstract class GammaFunctions { } if (x <= 0 && x == (long) x) { /* Negative integer argument */ - RError.warning(RError.SHOW_CALLER2, RError.Message.VALUE_OUT_OF_RANGE, "lgamma"); + RError.warning(RError.SHOW_CALLER, RError.Message.VALUE_OUT_OF_RANGE, "lgamma"); return Double.POSITIVE_INFINITY; /* +Inf, since lgamma(x) = log|gamma(x)| */ } @@ -388,7 +388,7 @@ public abstract class GammaFunctions { */ if (y > gfn_sign_xmax) { - RError.warning(RError.SHOW_CALLER2, RError.Message.VALUE_OUT_OF_RANGE, "lgamma"); + RError.warning(RError.SHOW_CALLER, RError.Message.VALUE_OUT_OF_RANGE, "lgamma"); return Double.POSITIVE_INFINITY; } diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Pbeta.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Pbeta.java index d7bfc2f1a7..02c713c68a 100644 --- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Pbeta.java +++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Pbeta.java @@ -12,14 +12,21 @@ package com.oracle.truffle.r.library.stats; import com.oracle.truffle.api.CompilerDirectives.TruffleBoundary; import com.oracle.truffle.api.profiles.BranchProfile; +import com.oracle.truffle.r.library.stats.StatsFunctions.Function3_2; import com.oracle.truffle.r.library.stats.TOMS708.Bratio; -import com.oracle.truffle.r.nodes.builtin.RExternalBuiltinNode; import com.oracle.truffle.r.runtime.RError; import com.oracle.truffle.r.runtime.RError.Message; // transcribed from pbeta.c -public abstract class Pbeta extends RExternalBuiltinNode.Arg5 { +public final class Pbeta implements Function3_2 { + + private final BranchProfile naProfile = BranchProfile.create(); + + @Override + public double evaluate(double x, double a, double b, boolean lowerTail, boolean logP) { + return pbeta(x, a, b, lowerTail, logP, naProfile); + } @TruffleBoundary private static double pbetaRaw(double x, double a, double b, boolean lowerTail, boolean logProb) { diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RCauchy.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RCauchy.java deleted file mode 100644 index 6792699dcc..0000000000 --- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RCauchy.java +++ /dev/null @@ -1,31 +0,0 @@ -/* - * 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) 1998--2008, The R Core Team - * Copyright (c) 2016, 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.M_PI; - -import com.oracle.truffle.r.library.stats.RandGenerationFunctions.RandFunction2_Double; -import com.oracle.truffle.r.library.stats.RandGenerationFunctions.RandomNumberProvider; - -public final class RCauchy implements RandFunction2_Double { - @Override - public double evaluate(double location, double scale, RandomNumberProvider rand) { - if (Double.isNaN(location) || !Double.isFinite(scale) || scale < 0) { - return RMath.mlError(); - } - if (scale == 0. || !Double.isFinite(location)) { - return location; - } else { - return location + scale * Math.tan(M_PI * rand.unifRand()); - } - } -} diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RMath.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RMath.java index 04e41edbf8..7b0d2a4c0d 100644 --- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RMath.java +++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RMath.java @@ -16,6 +16,7 @@ import static com.oracle.truffle.r.library.stats.LBeta.lbeta; import com.oracle.truffle.api.CompilerDirectives.CompilationFinal; import com.oracle.truffle.api.CompilerDirectives.TruffleBoundary; +import com.oracle.truffle.r.runtime.RRuntime; /** * Encapsulates functions to be found in Rmath.h or in nmath directory in GnuR except for random @@ -43,6 +44,38 @@ public class RMath { return ((y >= 0) ? TOMS708.fabs(x) : -TOMS708.fabs(x)); } + public static double fmod(double a, double b) { + double q = a / b; + if (b != 0) { + double tmp = a - Math.floor(q) * b; + if (RRuntime.isFinite(q) && Math.abs(q) > 1 / RRuntime.EPSILON) { + // TODO support warning here + throw new UnsupportedOperationException(); + } + return tmp - Math.floor(tmp / b) * b; + } else { + return Double.NaN; + } + } + + public static double tanpi(double x) { + if (Double.isNaN(x)) { + return x; + } + if (!Double.isFinite(x)) { + return mlError(); + } + + x = fmod(x, 1.); // tan(pi(x + k)) == tan(pi x) for all integer k + // map (-1,1) --> (-1/2, 1/2] : + if (x <= -0.5) { + x++; + } else if (x > 0.5) { + x--; + } + return (x == 0.) ? 0. : ((x == 0.5) ? Double.NaN : Math.tan(MathConstants.M_PI * x)); + } + // // GNUR from fmin2.c and fmax2 // 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 3b4d61c879..946e1234e5 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 @@ -40,6 +40,10 @@ import com.oracle.truffle.r.library.methods.SlotFactory.R_getSlotNodeGen; import com.oracle.truffle.r.library.methods.SlotFactory.R_setSlotNodeGen; import com.oracle.truffle.r.library.methods.SubstituteDirectNodeGen; import com.oracle.truffle.r.library.parallel.ParallelFunctionsFactory.MCIsChildNodeGen; +import com.oracle.truffle.r.library.stats.Cauchy; +import com.oracle.truffle.r.library.stats.Cauchy.DCauchy; +import com.oracle.truffle.r.library.stats.Cauchy.PCauchy; +import com.oracle.truffle.r.library.stats.Cauchy.RCauchy; import com.oracle.truffle.r.library.stats.CdistNodeGen; import com.oracle.truffle.r.library.stats.Chisq; import com.oracle.truffle.r.library.stats.CompleteCases; @@ -59,13 +63,13 @@ import com.oracle.truffle.r.library.stats.GammaFunctions.QgammaFunc; import com.oracle.truffle.r.library.stats.Geom; import com.oracle.truffle.r.library.stats.Geom.DGeom; import com.oracle.truffle.r.library.stats.Geom.RGeom; +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.Qbinom; import com.oracle.truffle.r.library.stats.Qnorm; import com.oracle.truffle.r.library.stats.RBeta; -import com.oracle.truffle.r.library.stats.RCauchy; import com.oracle.truffle.r.library.stats.RChisq; import com.oracle.truffle.r.library.stats.RGamma; import com.oracle.truffle.r.library.stats.RHyper; @@ -285,6 +289,14 @@ public class CallAndExternalFunctions { return RandGenerationFunctionsFactory.Function2_IntNodeGen.create(new Rbinom()); case "pbinom": return StatsFunctionsFactory.Function3_2NodeGen.create(new Pbinom()); + case "pbeta": + return StatsFunctionsFactory.Function3_2NodeGen.create(new Pbeta()); + case "dcauchy": + return StatsFunctionsFactory.Function3_1NodeGen.create(new DCauchy()); + case "pcauchy": + return StatsFunctionsFactory.Function3_2NodeGen.create(new PCauchy()); + case "qcauchy": + return StatsFunctionsFactory.Function3_2NodeGen.create(new Cauchy.QCauchy()); case "pf": return StatsFunctionsFactory.Function3_2NodeGen.create(new Pf()); case "dgamma": diff --git a/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/ops/BinaryArithmetic.java b/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/ops/BinaryArithmetic.java index c3538fa9a7..fa6e33017b 100644 --- a/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/ops/BinaryArithmetic.java +++ b/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/ops/BinaryArithmetic.java @@ -110,7 +110,8 @@ public abstract class BinaryArithmetic extends Operation { public abstract String op(String left, String right); public static double fmod(double a, double b) { - // LICENSE: transcribed code from GNU R, which is licensed under GPL + // TODO: this is duplicated in RMath, once RMath is moved to runtime, this should be + // replaced double q = a / b; if (b != 0) { double tmp = a - Math.floor(q) * b; 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 30665246ce..49e4017735 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 @@ -30,7 +30,7 @@ import com.oracle.truffle.r.test.TestBase; * Common tests for functions implemented using {@code StatsFunctions} infrastructure. */ public class TestStatFunctions extends TestBase { - private static final String[] FUNCTION3_1_NAMES = {"dgamma", "dbeta"}; + private static final String[] FUNCTION3_1_NAMES = {"dgamma", "dbeta", "dcauchy"}; private static final String[] FUNCTION3_1_PARAMS = { "10, 10, 10, log=TRUE", "3, 3, 3, log=FALSE", @@ -74,6 +74,32 @@ public class TestStatFunctions extends TestBase { assertEval(Output.IgnoreWhitespace, template("set.seed(1); %0(%1, %2, %3)", FUNCTION2_2_NAMES, FUNCTION2_2_PARAMS, new String[]{"lower.tail=TRUE", "lower.tail=FALSE"}, new String[]{"log.p=TRUE", "log.p=FALSE"})); // the error cases (where log.p nor lower.tail should make no difference) + // first parameter wrong assertEval(Output.IgnoreWhitespace, template("set.seed(1); %0(%1)", FUNCTION2_2_NAMES, new String[]{"c(NA, 0, NaN, 1/0, -1/0), rep(c(1, 0, 0.1), 5)"})); + // second parameter wrong + assertEval(Output.IgnoreWhitespace, template("set.seed(1); %0(%1)", FUNCTION2_2_NAMES, new String[]{"rep(c(1, 0, 0.1), 5), c(NA, 0, NaN, 1/0, -1/0)"})); + } + + private static final String[] FUNCTION3_2_NAMES = {"pbeta", "pcauchy", "qcauchy"}; + private static final String[] FUNCTION3_2_PARAMS = { + "0, 10, 10", + "c(-1, 0, 0.2, 2), c(-1, 0, 0.1, 0.9, 3), rep(c(-1, 0, 1, 0.1, -0.1, 0.0001), 20)", + "0, c(0.0653, 0.000123, 32e-80, 8833, 79e70, 0, -1), rep(c(0.0653, 0.000123, 32e-80, 8833, 79e70), 7)" + }; + + @Test + public void testFunctions32() { + // first: the "normal params" with all the combinations of log.p and lower.tail + assertEval(Output.IgnoreWhitespace, template("set.seed(1); %0(%1, %2, %3)", + FUNCTION3_2_NAMES, FUNCTION3_2_PARAMS, new String[]{"lower.tail=TRUE", "lower.tail=FALSE"}, new String[]{"log.p=TRUE", "log.p=FALSE"})); + // the error cases (where log.p nor lower.tail should make no difference) + // first parameter wrong + assertEval(Output.IgnoreWarningContext, + template("set.seed(1); %0(%1)", FUNCTION3_2_NAMES, new String[]{"c(NA, 0, NaN, 1/0, -1/0), rep(c(1, 0, 0.1), 5), rep(c(1, 0, 0.1), 5)"})); + // second parameter wrong + assertEval(Output.IgnoreWarningContext, + template("set.seed(1); %0(%1)", FUNCTION3_2_NAMES, new String[]{"rep(c(1, 0, 0.1), 5), c(NA, 0, NaN, 1/0, -1/0), rep(c(1, 0, 0.1), 5)"})); + // third parameter wrong + assertEval(Output.IgnoreWhitespace, template("set.seed(1); %0(%1)", FUNCTION3_2_NAMES, new String[]{"rep(c(1, 0, 0.1), 5), rep(c(1, 0, 0.1), 5), c(NA, 0, NaN, 1/0, -1/0)"})); } } diff --git a/mx.fastr/copyrights/overrides b/mx.fastr/copyrights/overrides index 5fa2ed2ee9..88bf8f67c6 100644 --- a/mx.fastr/copyrights/overrides +++ b/mx.fastr/copyrights/overrides @@ -31,7 +31,7 @@ com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/methods/Slot.java, com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Arithmetic.java,gnu_r_gentleman_ihaka.copyright com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RBeta.java,gnu_r_gentleman_ihaka.copyright com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RMultinom.java,gnu_r_gentleman_ihaka.copyright -com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RCauchy.java,gnu_r_ihaka_core.copyright +com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Cauchy.java,gnu_r_ihaka_core.copyright com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Cdist.java,gnu_r_gentleman_ihaka.copyright com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/CompleteCases.java,gnu_r_gentleman_ihaka2.copyright com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Covcor.java,gnu_r.copyright -- GitLab