University of Iowa

College of Liberal Arts and Sciences

• STAT:7400
• Syllabus
• Outline
• Notes
• Homework
• Resources

HW7 Notes

AddOne OMP Timings

if (! dir.exists("AddOne.Rcheck/")) {
    year <- 2022
    tarball <- "AddOne_4.0-2.tar.gz"
    if (! file.exists(tarball)) {
        path <- sprintf("classes/STAT7400-%s/examples/%s", year, tarball)
        url <- sprintf("http://www.stat.uiowa.edu/~luke/%s", path)
        download.file(url, tarball)
    }
    system(sprintf("R CMD check %s", tarball))
}
library(AddOne, lib = "AddOne.Rcheck")
    
R <- 100000
K <- 3
threads <- c(1, 2, 4)
if (isTRUE(parallel::detectCores() > 4))
    threads <- c(1, 2, 4, 8)
lens <- c(8, 32, 64, 128, 256, 512, 1024, 2048)

n <- rep(rep(lens, length(threads)), K)
P <- rep(rep(threads, each = length(lens)), K)

tm <- function(n, P) {
    x <- rnorm(n)
    system.time(for (i in 1 : R) p.AddOne(x, P = P))[3]
}

system.time(time <- sapply(seq_along(n), function(i) tm(n[i], P[i])))

d <- data.frame(time, n, P = as.factor(P))

library(ggplot2)
ggplot(d, aes(n, time, color = P)) + geom_point() + geom_smooth()

Signaling NaN Warnings

The convention is to signal the "NaNs produced" warning if new NaN values are produced but not if NaN arguments are being propagated.

Pure R Code

Keep in mind:

• R logical values can be TRUE, FALSE, or NA.
• Using an NA as the test condition with if signals an error.

So you need to screen out incoming NA/NaN values before doing a comparison:

ppareto <- function(x, alpha, beta) {
    bad_alpha <- alpha[! is.na(alpha)] <= 0
    bad_beta <- beta[! is.na(beta)] <= 0
    if (any(bad_alpha, bad_beta)) {
        warning("NaNs produced")
        alpha[bad_alpha] <- NaN
        beta[bad_beta] <- NaN
    }
    x <- ifelse(x <= alpha, alpha, x)
    lp <- beta * (log(alpha) - log(x))
    1 - exp(lp)
}

Using C Code

In C code comparisons will return 1 for False and 0 for TRUE.

The IEEE floating point standard specifies that comparisons with NaN values behave like this (x can be anything, including NaN and \(\pm\infty\))

NaN >= x	NaN <= x	NaN > x	NaN < x	NaN == x	NaN != x
FALSE	FALSE	FALSE	FALSE	FALSE	TRUE

#include <Rinternals.h>

void ppareto(double *x, double *alpha, double *beta, double *p, double *dn)
{
    R_xlen_t n = (R_xlen_t) *dn;
    int na_gen = FALSE;

    for (R_xlen_t i = 0; i < n; i++) {
    if (alpha[i] <= 0 || beta[i] <= 0) {
        p[i] = R_NaN;
        na_gen = TRUE;
    }
    else if (x[i] <= alpha[i])
        p[i] = 0;
    else {
        double lp = beta[i] * (log(alpha[i]) - log(x[i]));
        p[i] = 1 - exp(lp);
    }
    }
    if (na_gen)
    warning("NaNs produced");
}

dyn.load("ppareto.so")

ppareto <- function(x, alpha, beta) {
    n <- max(length(x), length(alpha), length(beta))
    x <- rep_len(as.double(x), n)
    alpha <- rep_len(as.double(alpha), n)
    beta <- rep_len(as.double(beta), n)
    .C("ppareto", x, alpha, beta, result = double(n),
       as.double(n), NAOK = TRUE)$result
}

This work is licensed under the CC BY-NC 3.0 Creative Commons License.