call_R is deprecated (has been but now CRAN is enforcing it). Investigate switching to pracma::romberg or re-writing the C-code. (spoiler alert: we re-write the C-code) #21

swihart · 2022-01-21T21:07:50Z

I've been learning about call_R.

Deprecated as per R News in June 2012:

• call_R and call_S are deprecated. They still
exist in the headers and as entry points, but are
no longer documented and should not be used
for new code

Cool. Cool, cool, cool. I have no idea what call_R does. It was mentioned in this old version of base::Foreign.
For repeated/rmutil/gnlrim it seems call_R is only called in the romberg files, specifically evalRfn. I think Jim wrote evalRfn because I can only find references to it via google searches to my github repos and the following from R-help:

R Help

Jim Lindsey's non-CRAN "rmutil" package has `int()' which looks great,
judged from the help (int.Rd), allowing to choose between Romberg and
TOMS 614 methods, allowing for indefinite integrals and singularities at
the end point(s).
On the other hand, the C code seems somewhat complicated, and (for romberg.c),
quite a bit nested (romberg -> evalRfn -> ..; interp) and completely
without comments.

So how about this? Experiment making gnlrim free of romberg in anticipation for doing the same rmutil and repeated. We can see how easy it would be to use pracma::romberg.

The text was updated successfully, but these errors were encountered:

swihart · 2022-01-21T21:10:11Z

l'd like to do speed trials. So maybe we can make a new function using pracma::romberg and verify it is not too too much slower and still gives same results. Then we can work on gutting out the .c files etc so that there are no vestigial appendages.

swihart · 2022-01-22T02:58:16Z

I've tried some things on a branch. They involve vectorizing/not vectorizing, cubature::hcubature and/or pracma::romberg and/or stats::integrate(). Nothing is working.

I then found a wonderful repo that shows how to do .C() and call_R() and .Call() and eval. I am going to write things out carefully by hand and then see if I can rewrite the guts myself.

CJ Geyer's qux repo that shows call_R and eval

swihart · 2022-01-23T00:12:33Z

It looks like I'm working up two solutions. I have a branch pracma_romberg and had to finagle an apply statement to work through the ff function which is returning a value for each id. I've come to understand that ff vectorized in the sense that it could handle each subject in parallel, but it is not vectorized for the argument x. On the pracma_romberg branch, I have an hcubature solution which is painfully slow. I have an apply-pracma::romberg solution that does well for small sample sizes and bonks on a big sample size example. To try to speed that one up, I waded into parallel:mclapply but it ironically did not speed up the small sample size example despite different number of cores and toggling prescheduling.

The writing is on the wall. I gotta rewrite some stuff in /src (Solution 2 -- as opposed to Solution 1 in the paragraph above).

God help us.

The options for this are the following:

just rewrite evalRfn (and just add an environ argument to romberg.c to pass to evalRfn)
rewrite all 3 functions.

Some things I found that might help think through SEXP-ing it up:

Section 5.09
Section 5.11.1 Finding Zeros
R-ext R Internal Structures
CJGeyer's .C / call_R examples in R file / C file
CJGeyer's .Call / eval examples in R file / C file

And my 3 files that I want to rewrite:

int1 inside gnlrim.R
src/romberg.c has 3 functions. They all might need to be rewritten. Lean on CJ Geyer's and the Finding Zero example. USE BRANCHES.

swihart · 2022-01-23T01:23:43Z

Made branch replace_dotc_with_dotcall

The goal is to simply rewrite the .C into a .Call:

int1 <- function(ff, aa, bb){
	.C("romberg",
		ff,
		as.double(aa),
		as.double(bb),
		len=as.integer(nnest),
		eps=as.double(eps),
		pts=as.integer(points),
		max=as.integer(steps),
		err=integer(1),
		res=double(nnest),
		PACKAGE="gnlrim")$res



	}

int1 <- function(ff, aa, bb){

         envir2 <- environment(fun=ff)
	.Call("romberg_sexp",
		ff,
		as.double(aa),
		as.double(bb),
		len=as.integer(nnest),
		eps=as.double(eps),
		pts=as.integer(points),
		max=as.integer(steps),
		err=integer(1),
		res=double(nnest),
               envir2,
		PACKAGE="gnlrim")



	}

All I did was

change .C to .Call
erase $res
set envir2 and make it final argument before PACKAGE
I copied contents of src/romberg.c into romberg_sexp.c and changed the romberg.c in the following way:
- change romberg to romberg_sexp
- add last argument SEXP envir

void romberg(void *fcn, double *a, double *b, int *len, double *eps,
             int *pts, int *max, int *err, double sumlen[])
{
  int i,j,j1,finish;
  double errsum,*tab1,*tab2,*x,*fx,*sum,*tmpsum,*zz,*pnt1,*pnt2;

  /* allocate vectors according to number of points and length of vector */
  tab1=(double*)R_alloc((size_t)(*pts),sizeof(double));
  tab2=(double*)R_alloc((size_t)(*pts),sizeof(double));
  x=(double*)R_alloc((size_t)(*max*(*len+1)),sizeof(double));
  fx=(double*)R_alloc((size_t)(*max*(*len+1)),sizeof(double));
  sum=(double*)R_alloc((size_t)(*len),sizeof(double));
  tmpsum=(double*)R_alloc((size_t)(*len),sizeof(double));
  zz=(double*)R_alloc((size_t)(*len),sizeof(double));
  pnt1=(double*)R_alloc((size_t)(*len),sizeof(double));
  pnt2=(double*)R_alloc((size_t)(*len),sizeof(double));
  if(!tab1||!tab2||!x||!fx||!sum||!tmpsum||!zz||!pnt1||!pnt2){
    *err=1;
    return;}
  *err=0;
  for(i=0;i<*len;i++)x[i**max]=1.0;
  /* iterate, decreasing step size, until convergence or max number of steps */
  for(j=0;j<*max;j++){
    j1=j+1;
    /* evaluate R function and calculate sum of trapezoids */
    evalRfn(fcn,a,b,j1,*len,sum,tmpsum,zz,pnt1,pnt2);
    finish=(j1>=*pts?1:0);
    /* repeatedly call polynomial interpolation routine */
    for(i=0;i<*len;i++){
      fx[j+i**max]=sum[i];
      if(j1>=*pts){
        interp(&x[j1-*pts+i**max],&fx[j1-*pts+i**max],*pts,tab1,tab2,&sumlen[i],&errsum,err);
        if(*err)return;
        /* check convergence */
        if(fabs(errsum)>*eps*fabs(sumlen[i]))finish=0;}
      /* decrease step size */
      x[j1+i**max]=x[j+i**max]/9.0;
      fx[j1+i**max]=fx[j+i**max];}
    if(finish)return;}
  *err=3;
  return;}

to

void romberg(void *fcn, double *a, double *b, int *len, double *eps,
             int *pts, int *max, int *err, double sumlen[])
{
  int i,j,j1,finish;
  double errsum,*tab1,*tab2,*x,*fx,*sum,*tmpsum,*zz,*pnt1,*pnt2;

  /* allocate vectors according to number of points and length of vector */
  tab1=(double*)R_alloc((size_t)(*pts),sizeof(double));
  tab2=(double*)R_alloc((size_t)(*pts),sizeof(double));
  x=(double*)R_alloc((size_t)(*max*(*len+1)),sizeof(double));
  fx=(double*)R_alloc((size_t)(*max*(*len+1)),sizeof(double));
  sum=(double*)R_alloc((size_t)(*len),sizeof(double));
  tmpsum=(double*)R_alloc((size_t)(*len),sizeof(double));
  zz=(double*)R_alloc((size_t)(*len),sizeof(double));
  pnt1=(double*)R_alloc((size_t)(*len),sizeof(double));
  pnt2=(double*)R_alloc((size_t)(*len),sizeof(double));
  if(!tab1||!tab2||!x||!fx||!sum||!tmpsum||!zz||!pnt1||!pnt2){
    *err=1;
    return;}
  *err=0;
  for(i=0;i<*len;i++)x[i**max]=1.0;
  /* iterate, decreasing step size, until convergence or max number of steps */
  for(j=0;j<*max;j++){
    j1=j+1;
    /* evaluate R function and calculate sum of trapezoids */
    evalRfn(fcn,a,b,j1,*len,sum,tmpsum,zz,pnt1,pnt2);
    finish=(j1>=*pts?1:0);
    /* repeatedly call polynomial interpolation routine */
    for(i=0;i<*len;i++){
      fx[j+i**max]=sum[i];
      if(j1>=*pts){
        interp(&x[j1-*pts+i**max],&fx[j1-*pts+i**max],*pts,tab1,tab2,&sumlen[i],&errsum,err);
        if(*err)return;
        /* check convergence */
        if(fabs(errsum)>*eps*fabs(sumlen[i]))finish=0;}
      /* decrease step size */
      x[j1+i**max]=x[j+i**max]/9.0;
      fx[j1+i**max]=fx[j+i**max];}
    if(finish)return;}
  *err=3;
  return;}

And then we do CMD+SHFT+B. Does it work?
A: NO. First delete the src/*.o files.

Q: Ok, .o files gone, does CMB+SHFT+B work?
A: No, need interp and evalRfn to be different than other file, so add _sexp to them for now.

Q: Ok, we have deleted .o and made sure we have uniquely named funcitons in .c files. Does CMD+SHFT+B work?
A: Yes.

Q: Run the first example in ?gnlrim. What happens?
A: Error:

> library(gnlrim)
> ?gnlrim
> dose <- c(9,12,4,9,11,10,2,11,12,9,9,9,4,9,11,9,14,7,9,8)
> y <- c(8.674419, 11.506066, 11.386742, 27.414532, 12.135699,  4.359469,
+        1.900681, 17.425948,  4.503345,  2.691792,  5.731100, 10.534971,
+        11.220260,  6.968932,  4.094357, 16.393806, 14.656584,  8.786133,
+        20.972267, 17.178012)
> id <- rep(1:4, each=5)
> gnlrim(y, mu=~a+b*dose+rand, random="rand", nest=id, pmu=c(a=8.7,b=0.25),
+        pshape=3.44, pmix=2.3)
 Error in int1(ff, neg1, zero) : 
  cannot allocate memory block of size 4917249 Tb 
4.
int1(ff, neg1, zero) at gnlrim.R#275
3.
inta(fn) at gnlrim.R#1013
2.
like(p) at gnlrim.R#1032
1.
gnlrim(y, mu = ~a + b * dose + rand, random = "rand", nest = id, 
    pmu = c(a = 8.7, b = 0.25), pshape = 3.44, pmix = 2.3)

---> Maybe something to do with R_alloc()?

Q: Great! Before anything else, How about CMD+SHFT+E?
A: NOPE -- but everything was good until the examples. In next comment, let's think if we really need R_alloc.

✓  checking compiled code ...
E  checking examples (1.8s)
   Running examples in ‘gnlrim-Ex.R’ failed
   The error most likely occurred in:
   
   > base::assign(".ptime", proc.time(), pos = "CheckExEnv")
   > ### Name: gnlrim
   > ### Title: Generalized Nonlinear Regression with a Random Parameter
   > ### Aliases: gnlrim
   > ### Keywords: models
   > 
   > ### ** Examples
   > 
   > 
   > dose <- c(9,12,4,9,11,10,2,11,12,9,9,9,4,9,11,9,14,7,9,8)
   > y <- c(8.674419, 11.506066, 11.386742, 27.414532, 12.135699,  4.359469,
   +       1.900681, 17.425948,  4.503345,  2.691792,  5.731100, 10.534971,
   +       11.220260,  6.968932,  4.094357, 16.393806, 14.656584,  8.786133,
   +       20.972267, 17.178012)
   > id <- rep(1:4, each=5)
   > 
   > gnlrim(y, mu=~a+b*dose+rand, random="rand", nest=id, pmu=c(a=8.7,b=0.25),
   +        pshape=3.44, pmix=2.3)
   Error in int1(ff, neg1, zero) : 
     cannot allocate memory block of size 4917249 Tb
   Calls: gnlrim -> like -> inta -> int1
   Execution halted
✓  checking for non-standard things in the check directory
✓  checking for detritus in the temp directory
   
   See
     ‘/Users/swihartbj/github/gnlrim.Rcheck/00check.log’
   for details.
   
── R CMD check results ─────────────────────────────────────── gnlrim 0.1.0 ────
Duration: 39.4s

> checking examples ... ERROR
  Running examples in ‘gnlrim-Ex.R’ failed
  The error most likely occurred in:
  
  > base::assign(".ptime", proc.time(), pos = "CheckExEnv")
  > ### Name: gnlrim
  > ### Title: Generalized Nonlinear Regression with a Random Parameter
  > ### Aliases: gnlrim
  > ### Keywords: models
  > 
  > ### ** Examples
  > 
  > 
  > dose <- c(9,12,4,9,11,10,2,11,12,9,9,9,4,9,11,9,14,7,9,8)
  > y <- c(8.674419, 11.506066, 11.386742, 27.414532, 12.135699,  4.359469,
  +       1.900681, 17.425948,  4.503345,  2.691792,  5.731100, 10.534971,
  +       11.220260,  6.968932,  4.094357, 16.393806, 14.656584,  8.786133,
  +       20.972267, 17.178012)
  > id <- rep(1:4, each=5)
  > 
  > gnlrim(y, mu=~a+b*dose+rand, random="rand", nest=id, pmu=c(a=8.7,b=0.25),
  +        pshape=3.44, pmix=2.3)
  Error in int1(ff, neg1, zero) : 
    cannot allocate memory block of size 4917249 Tb
  Calls: gnlrim -> like -> inta -> int1
  Execution halted

1 error x | 0 warnings ✓ | 0 notes ✓
Error: R CMD check found ERRORs
Execution halted

Exited with status 1.

swihart · 2022-01-23T01:49:21Z

Found this from Hadley: Vectors.md and more generally R internals.

From the R Internals page:

Here we focus on best practices and modern tools. To wit, we recommend that you use R_NO_REMAP so all API functions have the prefix R_ or Rf_:

#define R_NO_REMAP
#include <R.h>
#include <Rinternals.h>

swihart · 2022-01-23T02:17:34Z

Okay I deleted branch replace_dotc_with_dotcall.

And I'm going to remake it. I'm going to rebuild this stuff line by line. My hopes is that once I figure out the allocation stuff, the guts of the code just runs.

swihart · 2022-01-23T02:18:01Z

Rebirth: Made branch replace_dotc_with_dotcall.

Let's just make it all SEXP and only have romberg_sexp().

1st hard lesson: If you include the #define R_NO_REMAP, the code snippets you were dinking with from CJ Geyer will need a bunch of Rf_ and R_ prefixes. I think I got it precariously set up. We can CMD+SHFT+B.

This returns the length which is the number of ids nnest.

#define R_NO_REMAP
#include <R.h>
#include <Rinternals.h>
#include <math.h>
#include <stddef.h>


SEXP romberg_sexp(SEXP fcn, SEXP state, SEXP envir)
{


  SEXP call, result;
  PROTECT(call   = Rf_lang2(fcn,state));
  PROTECT(result = Rf_eval(call,envir));
  SEXP foo;
  PROTECT(foo = Rf_coerceVector(result, REALSXP));
  SEXP len;
  PROTECT(len = Rf_ScalarInteger(LENGTH(foo)));
  /* for (int i = 0; i < len; i++) */
  /*   if (! R_finite(REAL(foo)[i])) */
  /*     Rf_error("function returned vector with non-finite element"); */
    UNPROTECT(4);
    return len;
  }

swihart · 2022-01-23T02:58:48Z

And this pair works as expected:

the romberg_sexp:

SEXP romberg_sexp(SEXP fcn, SEXP state, SEXP envir)
{
  
  SEXP call, result;
  PROTECT(call   = Rf_lang2(fcn,state));
  PROTECT(result = Rf_eval(call,envir));
  SEXP foo;
  PROTECT(foo = Rf_coerceVector(result, REALSXP));
  int len = LENGTH(foo);
  for (int i = 0; i < len; i++)
    if (! R_finite(REAL(foo)[i]))
      Rf_error("function returned vector with non-finite element");
    UNPROTECT(3);
    return foo;
  }

and the int1 from within glnrim -- the print statements show I can evaluate the function as expected.

int1 <- function(ff, aa, bb){
	## .C("romberg",
	## 	ff,
	## 	as.double(aa),
	## 	as.double(bb),
	## 	len=as.integer(nnest),
	## 	eps=as.double(eps),
	## 	pts=as.integer(points),
	## 	max=as.integer(steps),
	## 	err=integer(1),
	## 	res=double(nnest),
	## 	PACKAGE="gnlrim")$res

xx <- c(1,2,3,4)
envir2 <- environment(fun=ff)
print("xx values:")
print(xx)
print("ff(xx) value:")
print(ff(xx))
print(".Call call:")
print(	.Call("romberg_sexp", ff, xx, envir2, PACKAGE="gnlrim"))


	}

swihart · 2022-01-23T03:16:02Z

Ok, I'm to the point of writing line by line and found a whole section on R_alloc! Apparently, it is recommended when we need storage for C objects.

From:

Memory Allocation - all about R_alloc

We see

6.1.1 Transient storage allocation
Here R will reclaim the memory at the end of the call to .C, .Call or .External. Use

char *R_alloc(size_t n, int size)

which allocates n units of size bytes each. A typical usage (from package stats) is

x = (int *) R_alloc(nrows(merge)+2, sizeof(int));

(size_t is defined in stddef.h which the header defining R_alloc includes.)

so if my first R_alloc() in romberg_sexp is

tab1=(double*)R_alloc((size_t)(*pts),sizeof(double));

Then this is saying tab1 allocates pts units of sizeof(double) each. I wonder if this is the thing giving me the error about allocating a vector that is a Tb big. And I wonder if it has to do with pts being *pts (a pointer...?)

...hours later...
i'm having R_alloc vector size problems and I read trhis near the bottom of the section:

These functions should only be used in code called by .C etc, never from front-ends. They are not thread-safe.

Maybe I need a different way to do vectors?

swihart · 2022-01-23T03:46:10Z

SECTION 5.13 EXTERNAL POINTERS AND WEAK REFERENCES

SO THINGS were going well until I tried to do math on *max and *len. What to do? I think I can just pass max and len and they don't need to be pointers in this .Call way of doing things?

Seems risky. I know pointers are important in C.

And that's about all I know.

Ok. so I change the SEXP max to int *max and int *len and did a cmd+shft+b and got the ol' Tb allocation error:

cannot allocate memory block of size 4917249 Tb

So, try removing pointers? Calculate these things within?

swihart · 2022-01-23T05:12:10Z

A blog post on switching from .C to .Call

Hadley's .Call() tutorial

swihart · 2022-01-23T05:53:03Z

Okay I'm going to bed. These things seem to be working well and the next line I write within romberg_sexp will be evalRfn

#define R_NO_REMAP
#include <R.h>
#include <Rinternals.h>
#include <math.h>
#include <stddef.h>


SEXP romberg_sexp(SEXP fcn, SEXP state, SEXP len,
		  SEXP pts, SEXP max, int *err, SEXP envir)
{

  int i, j, j1, finish;
  double errsum, *tab1, *tab2, *x, *fx, *sum, *tmpsum, *zz, *pnt1, *pnt2;

  int *PTS = INTEGER(pts);
  int *MAX = INTEGER(max);
  int *LEN = INTEGER(len);
  
  /*allocate vectors according to member of points and length of vector*/
    tab1 = (double*)R_alloc((size_t)(*PTS)         ,sizeof(double));
    tab2 = (double*)R_alloc((size_t)(*PTS)         ,sizeof(double));
       x = (double*)R_alloc((size_t)(*MAX*(*LEN+1)),sizeof(double));
      fx = (double*)R_alloc((size_t)(*MAX*(*LEN+1)),sizeof(double));
     sum = (double*)R_alloc((size_t)(*LEN)         ,sizeof(double));
  tmpsum = (double*)R_alloc((size_t)(*LEN)         ,sizeof(double));
      zz = (double*)R_alloc((size_t)(*LEN)         ,sizeof(double));
    pnt1 = (double*)R_alloc((size_t)(*LEN)         ,sizeof(double));
    pnt2 = (double*)R_alloc((size_t)(*LEN)         ,sizeof(double));

    if(!tab1||!tab2||!x||!fx||!sum||!tmpsum||!zz||!pnt1||!pnt2){
      *err=1;
      Rf_error("*err is now 1");}
    *err=0;
    for(i=0;i<*LEN;i++)x[i**MAX]=1.0;

  
  SEXP call, result;
  PROTECT(call   = Rf_lang2(fcn,state));
  PROTECT(result = Rf_eval(call,envir));
  SEXP foo;
  PROTECT(foo = Rf_coerceVector(result, REALSXP));
  int len2 = LENGTH(foo);
  for (int i = 0; i < len2; i++)
    if (! R_finite(REAL(foo)[i]))
      Rf_error("function returned vector with non-finite element");
    UNPROTECT(3);
    return foo;
  }

int1 <- function(ff, aa, bb){
	## .C("romberg",
	## 	ff,
	## 	as.double(aa),
	## 	as.double(bb),
	## 	len=as.integer(nnest),
	## 	eps=as.double(eps),
	## 	pts=as.integer(points),
	## 	max=as.integer(steps),
	## 	err=integer(1),
	## 	res=double(nnest),
	## 	PACKAGE="gnlrim")$res

xx <- c(1,2,3,4)
envir2 <- environment(fun=ff)
print("xx values:")
print(xx)
print("ff(xx) value:")
print(ff(xx))
print(".Call call:")
    print(	.Call("romberg_sexp",
                      ff,
                      xx,
                      len=as.integer(nnest),
                      as.integer(points),
                      max=as.integer(steps),
                      err=integer(1),
                      envir2,
                      PACKAGE="gnlrim"))


	}

swihart · 2022-01-24T23:58:36Z

Ergh. I've been struggling on this for a while. From I can tell, the call_R() in romberg.c was nice because the void zz was made into the argument for the R-land ff and could ff could be given zz values that were generated in C. The example I found from CJ Geyer's page was where state was essentially the x-values from R-land being passed into C-land. So now I am tinkering the with Section 5.11.1 zero finding code that I linked above. Since I'm using the R_NO_Map I had to put in some Rf_ prefixes. Essentially the example then is:

R code (which I put right inside gnlrim and am printing values in teh same dev workflow I had going when I was working directly on romberg_sexp.c):

zero <- function(f, guesses, tol = 1e-7) {
    f.check <- function(x) {
        x <- f(x)
        if(!is.numeric(x)) stop("Need a numeric result")
        as.double(x)
    }
    .Call("zero", body(f.check), as.double(guesses), as.double(tol),
          new.env(), PACKAGE="gnlrim")
}

cube1 <- function(x) (x^2 + 1) * (x - 1.5)
print("value of `zero(cube1, c(-4, 50))`: ")
print(zero(cube1,c(-4,50)))

And saved in src/zero.c:

#define R_NO_REMAP
#include <R.h>
#include <Rinternals.h>
#include <math.h>
#include <stddef.h>


/* zero <- function(f, guesses, tol = 1e-7) { */
/*     f.check <- function(x) { */
/*         x <- f(x) */
/*         if(!is.numeric(x)) stop("Need a numeric result") */
/*         as.double(x) */
/*     } */
/*     .Call("zero", body(f.check), as.double(guesses), as.double(tol), */
/*           new.env()) */
/* } */

SEXP mkans(double x)
{
    // no need for PROTECT() here, as REAL(.) does not allocate:
    SEXP ans = Rf_allocVector(REALSXP, 1);
    REAL(ans)[0] = x;
    return ans;
}

double feval(double x, SEXP f, SEXP rho)
{
    // a version with (too) much PROTECT()ion .. "better safe than sorry"
    SEXP symbol, value;
     Rf_protect(symbol =  Rf_install("x"));
     Rf_protect(value = mkans(x));
     Rf_defineVar(symbol, value, rho);
     Rf_unprotect(2);
    return(REAL(Rf_eval(f, rho))[0]);
}

SEXP zero(SEXP f, SEXP guesses, SEXP stol, SEXP rho)
{
    double x0 = REAL(guesses)[0], x1 = REAL(guesses)[1],
           tol = REAL(stol)[0];
    double f0, f1, fc, xc;

    if(tol <= 0.0)  Rf_error("non-positive tol value");
    f0 = feval(x0, f, rho); f1 = feval(x1, f, rho);
    if(f0 == 0.0) return mkans(x0);
    if(f1 == 0.0) return mkans(x1);
    if(f0*f1 > 0.0)  Rf_error("x[0] and x[1] have the same sign");

    for(;;) {
        xc = 0.5*(x0+x1);
        if(fabs(x0-x1) < tol) return  mkans(xc);
        fc = feval(xc, f, rho);
        if(fc == 0) return  mkans(xc);
        if(f0*fc > 0.0) {
            x0 = xc; f0 = fc;
        } else {
            x1 = xc; f1 = fc;
        }
    }
}

We can see there is a process to make the passed function into being a C-land object that can take values generated in C.

I have made an analogous zero_ff, and my goal is to modify it to where it will return values on my ff as per gnlrim. See next comment.

swihart · 2022-01-25T15:55:35Z

BREAKTHROUGH CITY.

I realized that doing the finding-zeroes with the feval and whatever was rough-going because I want it to evaluate as it would in R. It's a little brainy that it returns a vector and each element is the id (nnest). If I just pass in one value at a time in a for loop in C it is just going to evaluate the first slot of ff at that value. So the epiphany came when I corrected a ill-posed assumption of mine: state doesn't have to originate from R. I can generate it from within C and then use the CJ Geyer method. To test this out though ,I needed away to come up with a vector of 4 in C that wouldn't crash everything. (aside: an annoying thing has been that when I edit .c code and try to do another cycle in R of testing I get the R Abort screen of death and have to restart, close Rstudio, delete all my src/*.o files, and then double click .Rproj to climb the mountain again with my fingers crossed it doesn't bonk. It started getting bad when I put in zero.c and zero_ff.c so I might try ousting those in the near future if this test works).

Narrator: the tests worked.

To see how I can make an R vector from within C without crashing everything I started at the top and read down through 5.09. And that was a treat:

Section 5.09.1 show this example:

#include <R.h>
#include <Rinternals.h>

    SEXP ab;
      ....
    ab = PROTECT(allocVector(REALSXP, 2));
    REAL(ab)[0] = 123.45;
    REAL(ab)[1] = 67.89;
    UNPROTECT(1);

So I did this in my burgeoning R-code:

int1 <- function(ff, aa, bb){
	## .C("romberg",
	## 	ff,
	## 	as.double(aa),
	## 	as.double(bb),
	## 	len=as.integer(nnest),
	## 	eps=as.double(eps),
	## 	pts=as.integer(points),
	## 	max=as.integer(steps),
	## 	err=integer(1),
	## 	res=double(nnest),
	## 	PACKAGE="gnlrim")$res

xx <- c(1,2,2.5,3.01)
envir2 <- environment(fun=ff)
print("xx values:")
print(xx)
print("ff(xx) value:")
print(ff(xx))
print(".Call call:")
print(
        .Call("romberg_sexp",
                      ff,
                      as.double(aa),
 	              as.double(bb),
                      xx,
                      len=as.integer(nnest),
                      eps=as.double(eps),
                      as.integer(points),
                      max=as.integer(steps),
                      err=integer(1),
                      envir2,
                      PACKAGE="gnlrim")
     )


	}

And this in the romberg_sexp.c -- pay attention to the very bottom where I set the last element of the vector to a different value to make sure I had control of the evaluation:

SEXP romberg_sexp(SEXP fcn, SEXP a, SEXP b, SEXP state, SEXP len, SEXP eps,
		  SEXP pts, SEXP max, int *err, SEXP envir)
{

  int i, j, j1, finish;
  double errsum, *tab1, *tab2, *x, *fx, *sum, *tmpsum, *zz, *pnt1, *pnt2;

  int *PTS = INTEGER(pts);
  int *MAX = INTEGER(max);
  int *LEN = INTEGER(len);
  double *EPS = REAL(eps);
  double *A = REAL(a);
  double *B = REAL(b);

  /*allocate vectors according to member of points and length of vector*/
    tab1 = (double*)R_alloc((size_t)(*PTS)         ,sizeof(double));
    tab2 = (double*)R_alloc((size_t)(*PTS)         ,sizeof(double));
       x = (double*)R_alloc((size_t)(*MAX*(*LEN+1)),sizeof(double));
      fx = (double*)R_alloc((size_t)(*MAX*(*LEN+1)),sizeof(double));
     sum = (double*)R_alloc((size_t)(*LEN)         ,sizeof(double));
  tmpsum = (double*)R_alloc((size_t)(*LEN)         ,sizeof(double));
      zz = (double*)R_alloc((size_t)(*LEN)         ,sizeof(double));
    pnt1 = (double*)R_alloc((size_t)(*LEN)         ,sizeof(double));
    pnt2 = (double*)R_alloc((size_t)(*LEN)         ,sizeof(double));

    if(!tab1||!tab2||!x||!fx||!sum||!tmpsum||!zz||!pnt1||!pnt2){
      *err=1;
      Rf_error("*err is now 1");}
    *err=0;
    for(i=0;i<*LEN;i++)x[i**MAX]=1.0;

  /* iterate, decreasing step size, until convergence or max number of steps */
  for(j=0;j<*MAX;j++){
    j1=j+1;
    /* evaluate R function and calculate sum of trapezoids */
    evalRfn_sexp(fcn,A,B,j1,*LEN,sum,tmpsum,zz,pnt1,pnt2,envir);
     finish=(j1>=*PTS?1:0);
    /*repeatedly call polynomial interpolation routine*/
     for(i=0;i<*LEN;i++){
       fx[j+i**MAX]=sum[i];
       if(j1>=*PTS){
    /* 	interp(&x[j1-*PTS+i**MAX],&fx[j1-*PTS+i**MAX],*PTS,tab1,tab2,&sumlen[i],&errsum,err); */
    if(*err)Rf_error("*err is now 1");
    /*  check convergence  */
	 /* 	if(fabs(errsum)>*EPS*fabs(sumlen[i]))finish=0;*/}
       /* decrease step size */
       x[j1+i**MAX]=x[j+i**MAX]/9.0;
       fx[j1+i**MAX]=fx[j+i**MAX];}
                     }

  SEXP call, result;
  /* try a state2 def and put in c(1,2,3,9) and see if you can get it*/
  SEXP abcd;
  abcd = PROTECT(Rf_allocVector(REALSXP, 4));
  REAL(abcd)[0] = 1;
  REAL(abcd)[1] = 2;
  REAL(abcd)[2] = 2.5;
  REAL(abcd)[3] = 4.2; /*this one is different than value in R*/
  /*...*/

  PROTECT(call   = Rf_lang2(fcn,abcd));
  PROTECT(result = Rf_eval(call,envir));
  SEXP foo;
  PROTECT(foo = Rf_coerceVector(result, REALSXP));
  int len2 = LENGTH(foo);
  for (int i = 0; i < len2; i++)
    if (! R_finite(REAL(foo)[i]))
      Rf_error("function returned vector with non-finite element");
    UNPROTECT(4); /*release protections*/
    return foo;
  }

swihart · 2022-01-25T16:06:24Z

So now the bold frontier. I gotta put this stuff (CJ Geyer's statements) where the call_R is.

swihart · 2022-01-25T16:09:13Z

I'm taking out all the zero stuff to see if my compiling issues resolve. I'm dumping them below and then I'll start on the fleshing out the code by putting CJ Geyer statements where call_R() is. First the dump:

R code (commented out):

#
#
# zero <- function(f, guesses, tol = 1e-7) {
#     f.check <- function(x) {
#         x <- f(x)
#         if(!is.numeric(x)) stop("Need a numeric result")
#         as.double(x)
#     }
#     .Call("zero", body(f.check), as.double(guesses), as.double(tol),
#           new.env(), PACKAGE="gnlrim")
# }
#
# cube1 <- function(x) (x^2 + 1) * (x - 1.5)
# print("value of `zero(cube1, c(-4, 50))`: ")
# print(zero(cube1,c(-4,50)))
#
#
# zero_ff <- function(f, guesses, tol = 1e-7) {
#    f.check <- function(x) {
#         x <- f(x)
#         if(!is.numeric(x)) stop("Need a numeric result")
#         as.double(x)
#     }
#     .Call("zero_ff", body(f.check), as.double(guesses), as.double(tol),
#           new.env(), PACKAGE="gnlrim")
# }
#
# cube2 <- function(x) as.numeric(lapply(x, function(W){W^2}))
# cube2(c(sqrt(11),sqrt(33)))
# print("value of `zero_ff(cube2, whatwhat)` from C-land: ")
# print(zero_ff(cube2,c(sqrt(11),sqrt(33))))
#
#
#

Now the src/zero.c:

#define R_NO_REMAP
#include <R.h>
#include <Rinternals.h>
#include <math.h>
#include <stddef.h>


/* zero <- function(f, guesses, tol = 1e-7) { */
/*     f.check <- function(x) { */
/*         x <- f(x) */
/*         if(!is.numeric(x)) stop("Need a numeric result") */
/*         as.double(x) */
/*     } */
/*     .Call("zero", body(f.check), as.double(guesses), as.double(tol), */
/*           new.env()) */
/* } */

SEXP mkans(double x)
{
    // no need for PROTECT() here, as REAL(.) does not allocate:
    SEXP ans = Rf_allocVector(REALSXP, 1);
    REAL(ans)[0] = x;
    return ans;
}

double feval(double x, SEXP f, SEXP rho)
{
    // a version with (too) much PROTECT()ion .. "better safe than sorry"
    SEXP symbol, value;
     Rf_protect(symbol =  Rf_install("x"));
     Rf_protect(value = mkans(x));
     Rf_defineVar(symbol, value, rho);
     Rf_unprotect(2);
    return(REAL(Rf_eval(f, rho))[0]);
}

SEXP zero(SEXP f, SEXP guesses, SEXP stol, SEXP rho)
{
    double x0 = REAL(guesses)[0], x1 = REAL(guesses)[1],
           tol = REAL(stol)[0];
    double f0, f1, fc, xc;

    if(tol <= 0.0)  Rf_error("non-positive tol value");
    f0 = feval(x0, f, rho); f1 = feval(x1, f, rho);
    if(f0 == 0.0) return mkans(x0);
    if(f1 == 0.0) return mkans(x1);
    if(f0*f1 > 0.0)  Rf_error("x[0] and x[1] have the same sign");

    for(;;) {
        xc = 0.5*(x0+x1);
        if(fabs(x0-x1) < tol) return  mkans(xc);
        fc = feval(xc, f, rho);
        if(fc == 0) return  mkans(xc);
        if(f0*fc > 0.0) {
            x0 = xc; f0 = fc;
        } else {
            x1 = xc; f1 = fc;
        }
    }
}

Now the src/zero_ff.c:

#define R_NO_REMAP
#include <R.h>
#include <Rinternals.h>
#include <math.h>
#include <stddef.h>



/* zero <- function(f, guesses, tol = 1e-7) { */
/*     f.check <- function(x) { */
/*         x <- f(x) */
/*         if(!is.numeric(x)) stop("Need a numeric result") */
/*         as.double(x) */
/*     } */
/*     .Call("zero", body(f.check), as.double(guesses), as.double(tol), */
/*           new.env()) */
/* } */


SEXP mkans_ff_vec(double x)
{
  //size_t n = sizeof(x) / sizeof(x[0]);
    // no need for PROTECT() here, as REAL(.) does not allocate:
    SEXP ans = Rf_allocVector(REALSXP, 1);
    REAL(ans)[0] = x;
    return ans;
}

SEXP mkans_ff(double x)
{
  //size_t n = sizeof(x) / sizeof(x[0]);
    // no need for PROTECT() here, as REAL(.) does not allocate:
    SEXP ans = Rf_allocVector(REALSXP, 1);
    REAL(ans)[0] = x;
    return ans;
}

double feval_ff(double x, SEXP f, SEXP rho)
{
    // a version with (too) much PROTECT()ion .. "better safe than sorry"
    SEXP symbol, value;
     Rf_protect(symbol =  Rf_install("x"));
     Rf_protect(value = mkans_ff(x));
     Rf_defineVar(symbol, value, rho);
     Rf_unprotect(2);
    return(REAL(Rf_eval(f, rho))[0]);
}

SEXP zero_ff(SEXP f, SEXP guesses, SEXP stol, SEXP rho)
{
  //R_xlen_t n = Rf_xlength(guesses);
  // return mkans_ff(n);

  int n = LENGTH(guesses);
  //SEXP returnableresults = Rf_allocVector(REALSXP,n);
  //double returnableresults[n];
  //  return mkans_ff(n);


  
    double x0 = REAL(guesses)[0], x1 = REAL(guesses)[1],
           tol = REAL(stol)[0];
    double f0, f1, fc, xc;

    //return(mkans_ff(feval_ff(x1,f,rho)));

    //for(int i; i<n; i++) returnableresults[i] = mkans_ff(feval_ff(x1,f,rho));

    /* for(int i; i<n; i++){ */
    /*   REAL(returnableresults)[i] = feval_ff(REAL(guesses)[i],f,rho); */
    /* } */
    /* return mkans_ff(returnableresults); */



    return(mkans_ff(1*feval_ff(x1,f,rho)));
    
    /* if(tol <= 0.0)  Rf_error("non-positive tol value"); */
    /* f0 = feval_ff(x0, f, rho); f1 = feval_ff(x1, f, rho); */
    /* if(f0 == 0.0) return mkans_ff(x0); */
    /* if(f1 == 0.0) return mkans_ff(x1); */
    /* if(f0*f1 > 0.0)  Rf_error("x[0] and x[1] have the same sign"); */

    /* for(;;) { */
    /*     xc = 0.5*(x0+x1); */
    /*     if(fabs(x0-x1) < tol) return  mkans_ff(xc); */
    /*     fc = feval_ff(xc, f, rho); */
    /*     if(fc == 0) return  mkans_ff(xc); */
    /*     if(f0*fc > 0.0) { */
    /*         x0 = xc; f0 = fc; */
    /*     } else { */
    /*         x1 = xc; f1 = fc; */
    /*     } */
    /* } */
}

swihart · 2022-01-25T16:29:50Z

Welp, taking the zero and zero_ff didn't fix the bonking. Might have to look more closely at the registering of functions. Let's do that .... here are the current contents of src/RegisterDynamicSymbol.c:

// RegisteringDynamic Symbols

#include <R.h>
#include <Rinternals.h>
#include <R_ext/Rdynload.h>

void R_init_markovchain(DllInfo* info) {
  R_registerRoutines(info, NULL, NULL, NULL, NULL);
  R_useDynamicSymbols(info, TRUE);
}

Oh yeah. Look at section 5.4: why do I have mylib suffix as markovchain ??? lol.

Section 5.04

Which lists this:

void
R_init_myLib(DllInfo *info)
{
   R_registerRoutines(info, cMethods, callMethods, NULL, NULL);
}

Yikes! I need to change markovchain to gnlrim and put in cMethods and callMethods...

Like this? Try it and see if things improve.

// RegisteringDynamic Symbols

#include <R.h>
#include <Rinternals.h>
#include <R_ext/Rdynload.h>

void R_init_gnlrim(DllInfo* info) {
  R_registerRoutines(info, cMethods, callMethods, NULL, NULL);
  R_useDynamicSymbols(info, TRUE);
}

Results: Still bonks. And it didn't know cMethods or call Methods. This inspired me to look at the _init files for rmutil and repeated and they look more like what the section 5.04 are talking about . So new plan! I will copy over repeated/src/repeated_init.c and make is gnlrim/src/gnlrim_init.c and register things as .C instead of .CAll. Next block -->

swihart · 2022-01-25T17:42:47Z

Prima facie gnlrim/src/gnlrim_init.c:

#include <R_ext/RS.h>
#include <stdlib.h> // for NULL
#include <R_ext/Rdynload.h>

/* .C calls */
extern void countfb_c(void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *);
extern void ddb_c(void *, void *, void *, void *, void *, void *, void *);
extern void ddp_c(void *, void *, void *, void *, void *, void *, void *);
extern void dmb_c(void *, void *, void *, void *, void *, void *, void *);
extern void dmp_c(void *, void *, void *, void *, void *, void *, void *);
extern void dpvfp_c(void *, void *, void *, void *, void *, void *, void *);
extern void gar_c(void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *);
extern void kcountb_c(void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *);
extern void krand_c(void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *);
extern void kserie_c(void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *);
extern void pginvgauss_c(void *, void *, void *, void *, void *, void *, void *, void *, void *, void *);
extern void ppowexp_c(void *, void *, void *, void *, void *, void *, void *, void *, void *, void *);
extern void psimplex_c(void *, void *, void *, void *, void *, void *, void *, void *, void *, void *);
extern void romberg_c(void *, void *, void *, void *, void *, void *, void *, void *, void *);

static const R_CMethodDef CEntries[] = {
    {"countfb_c",    (DL_FUNC) &countfb_c,    21},
    {"ddb_c",        (DL_FUNC) &ddb_c,         7},
    {"ddp_c",        (DL_FUNC) &ddp_c,         7},
    {"dmb_c",        (DL_FUNC) &dmb_c,         7},
    {"dmp_c",        (DL_FUNC) &dmp_c,         7},
    {"dpvfp_c",      (DL_FUNC) &dpvfp_c,       7},
    {"gar_c",        (DL_FUNC) &gar_c,        25},
    {"kcountb_c",    (DL_FUNC) &kcountb_c,    24},
    {"krand_c",      (DL_FUNC) &krand_c,      25},
    {"kserie_c",     (DL_FUNC) &kserie_c,     26},
    {"pginvgauss_c", (DL_FUNC) &pginvgauss_c, 10},
    {"ppowexp_c",    (DL_FUNC) &ppowexp_c,    10},
    {"psimplex_c",   (DL_FUNC) &psimplex_c,   10},
    {"romberg_c",    (DL_FUNC) &romberg_c,     9},
    {NULL, NULL, 0}
};


/* .Fortran calls */
extern void F77_NAME(binnest_f)(void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *);
extern void F77_NAME(chidden_f)(void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *);
extern void F77_NAME(cphidden_f)(void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *);
extern void F77_NAME(gcopula_f)(void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *);
extern void F77_NAME(hidden_f)(void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *);
extern void F77_NAME(logitord_f)(void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *);

static const R_FortranMethodDef FortranEntries[] = {
    {"binnest_f",  (DL_FUNC) &F77_NAME(binnest_f),  64},
    {"chidden_f",  (DL_FUNC) &F77_NAME(chidden_f),  40},
    {"cphidden_f", (DL_FUNC) &F77_NAME(cphidden_f), 40},
    {"gcopula_f",  (DL_FUNC) &F77_NAME(gcopula_f),  17},
    {"hidden_f",   (DL_FUNC) &F77_NAME(hidden_f),   35},
    {"logitord_f", (DL_FUNC) &F77_NAME(logitord_f), 19},
    {NULL, NULL, 0}
};

void R_init_repeated(DllInfo *dll)
{
    R_registerRoutines(dll, CEntries, NULL, FortranEntries, NULL);
    R_useDynamicSymbols(dll, FALSE);
}

First things first: change all the repeated to gnlrim and delete any functions that aren't in gnlrim.

The result -- much more condensed! -- I also changed argument in regesterRoutines to NULL since we eliminated all Fortran. So here's the result and this would work (?) if I still used .C:

#include <R_ext/RS.h>
#include <stdlib.h> // for NULL
#include <R_ext/Rdynload.h>

/* .C calls */
extern void romberg_c(void *, void *, void *, void *, void *, void *, void *, void *, void *);

static const R_CMethodDef CEntries[] = {
    {"romberg_c",    (DL_FUNC) &romberg_c,     9},
    {NULL, NULL, 0}
};


void R_init_repeated(DllInfo *dll)
{
    R_registerRoutines(dll, CEntries, NULL, NULL, NULL);
    R_useDynamicSymbols(dll, FALSE);
}

But we need .Call. So leaning on section 5.04, let's convert this stuff. Here it is:

#include <R_ext/RS.h>
#include <stdlib.h> // for NULL
#include <R_ext/Rdynload.h>

// /* .C() */
// extern void romberg_c(void *, void *, void *, void *, void *, void *, void *, void *, void *);
//
// static const R_CMethodDef CEntries[] = {
//     {"romberg_c",    (DL_FUNC) &romberg_c,     9},
//     {NULL, NULL, 0}
// };



/* .Call()  */
extern void romberg_sexp(void *, void *, void *, void *, void *, void *, void *, void *, void *, void *);

static const R_CallMethodDef callMethods[]  = {
    {"romberg_sexp", (DL_FUNC) &romberg_sexp, 10},
    {NULL, NULL, 0}
};

void R_init_gnlrim(DllInfo *dll)
{
    R_registerRoutines(dll, NULL, callMethods, NULL, NULL);
    R_useDynamicSymbols(dll, FALSE);
}

And I have removed registerDynamicSymbol.c from the src file. Looks like things are working better. Also, don't do the cmd+shft+D so often.

swihart · 2022-01-26T00:38:20Z

Oh wow.

I have it running! Giving same answers as gnlmix!

BUT!

[I'm getting some stack warning/errors]. As mentioned here in an issue for another R-package, it involves

Hi @ryankennedyio, thanks for reporting and working with the dev version of fst.

Warning: stack imbalance in '.Call', 30 then 29

This warning is due to an unmatched pair of PROTECT and UNPROTECT macro's (R API) in C++. I'll have to sort through the new code using that API to find that unmatched pair.

Okay. Took some care but I think I got it done. At first I was mad -- despite running error free (a dream come true that I couldn't even fathom a few days ago) -- I was timing gnlrim as 4x slower than gnlmix. Then I realized that gnlrim(with .C()) vs gnlrim(with .Call) were about the same. For some reason gnlrim is a bit slower than gnlmix. Do you think has to do with the log(dispersion)? THat was the key difference. That and truck load of hard random intercept distributions...
Oh well. This is great!

swihart · 2022-01-26T02:27:34Z

cry emoji

It runs on the first example of gnlrim; but when I tried some cloglog bridge stuff it segfaults a bunch. or gives weird, warnings/errors from C-land

swihart · 2022-01-26T04:09:07Z

I have installed R. Installed Rstudio. Still having problems. Starting to debug. Finding more good material on .Call:

(https://www.math.kit.edu/stoch/~lindner/media/.c.call%20extensions.pdf)

swihart · 2022-01-26T19:33:56Z

TECMO DANCE!

This one weird trick makes R-developers hate them.

I solved it. Things are running smoothly, not anymore crashy or slower than the .C implementation. The thing was I had an unprotected objected. The reason it was unprotected was because in the finding-zeros makeans function the comment said "no reason to protect" and maybe that is right for that particular invocation but for my usage I have to protect it and consequently unprotect it. Now things seems a lot smoother. I had this breakthough at midnight last night, and just to be safe committed it to the branch I was working on locally and just because it seemed like a good time for a Murphy's Law computer crash, I pushed the branch to github.

I think the standing issues are now to take out the "state=xx" argument (which means editing gnlrim_init registration number of arguments etc) but that should be straight forward. Then I guess I need to summarize what I learned and see if I can just plug and play in repeated and rmutil and get those on CRAN before the 4th of Feb 2022.

swihart · 2022-01-26T20:09:50Z

I did another push to the branch after updating roxygen and the 10 -> 9 arguments.

swihart · 2022-01-26T20:52:11Z

Stray observations:

I've learned a lot and I'm glad I had gnlrim to hack on. Now I need to transfer and execute what I've learned to repeated and rmutil. I browsed a little bit and here's what I'm gathering:

repeated

Should be similar to what I did in gnlrim -- the key difference is it is call "romberg_c" in certain places. Otherwise, it looks as thought there's only on .C("romberg_c" in the whole package and it is only in the int1 of gnlmix (perfectly analogous to my situation in gnlrim)
uses roxygen

rmutil

a little more involved than repeated. That's why I think I should do repeated first and get it all on CRAN safe and sound. So the edits to romberg.c should be the same as what i've completed for gnlrim. However there are call_R calls in Toms614 algorithm at the entry "inthp".
Did you know int1 is vectorized? Sweet! Did you know int2 is for two dimensions? It is a little interesting -- it looks like internally int2 does a romberg step locally and integrates that over the 2nd dimensiont first ... any way, it made my mind race and think that maybe we can get the "two random effects" case going iwth int2 and maybe it is generalizable to int3, int4, etc. But regardless, I think this will be more work becasue you haven't worked with .C("inthp" yet although I think it will be easier than the romberg stuff because romberg was calling on a subroutine evalRfn to get things done.
does NOT use roxygen (unlike repeated and gnlrim)

swihart · 2022-01-26T21:16:02Z

Check out these nice slides for a potential algo improvement:

romberg slides

Nice deep-dive explanation of Romberg and how (very) accurate it can be compared to other algos:

romberg blog

swihart · 2022-02-09T14:36:32Z

Just tie this up: I'm leaving branch replace_dotc_with_dotcall up -- but it grew a little unruly when I had some undisciplined attempts to chase down a valgrind error and compare .C with .Call functionality. For cran purposes, I can't have any call_R stuff in the c code, even if I'm .Rbuildignoring it. So, let's close this and tidy up the master branch given that repeated v1.1.4 is on CRAN.

swihart changed the title ~~call_R is deprecated (has been but now CRAN is enforcing it). switch to pracma::romberg~~ call_R is deprecated (has been but now CRAN is enforcing it). Investigate switching to pracma::romberg or re-writing the C-code. (spoiler alert: we re-write the C code) Jan 26, 2022

This was referenced Jan 27, 2022

Re-code to omit the use of call_R. swihart/repeated#16

Closed

Re-code to omit the use of call_R. swihart/rmutil#11

Closed

swihart closed this as completed Feb 9, 2022

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

call_R is deprecated (has been but now CRAN is enforcing it). Investigate switching to pracma::romberg or re-writing the C-code. (spoiler alert: we re-write the C-code) #21

call_R is deprecated (has been but now CRAN is enforcing it). Investigate switching to pracma::romberg or re-writing the C-code. (spoiler alert: we re-write the C-code) #21

swihart commented Jan 21, 2022

swihart commented Jan 21, 2022

swihart commented Jan 22, 2022

swihart commented Jan 23, 2022 •

edited

swihart commented Jan 23, 2022 •

edited

swihart commented Jan 23, 2022 •

edited

swihart commented Jan 23, 2022

swihart commented Jan 23, 2022 •

edited

swihart commented Jan 23, 2022

swihart commented Jan 23, 2022 •

edited

swihart commented Jan 23, 2022 •

edited

swihart commented Jan 23, 2022 •

edited

swihart commented Jan 23, 2022

swihart commented Jan 24, 2022 •

edited

swihart commented Jan 25, 2022 •

edited

swihart commented Jan 25, 2022

swihart commented Jan 25, 2022 •

edited

swihart commented Jan 25, 2022 •

edited

swihart commented Jan 25, 2022 •

edited

swihart commented Jan 26, 2022 •

edited

swihart commented Jan 26, 2022

swihart commented Jan 26, 2022 •

edited

swihart commented Jan 26, 2022

swihart commented Jan 26, 2022

swihart commented Jan 26, 2022 •

edited

swihart commented Jan 26, 2022 •

edited

swihart commented Feb 9, 2022

call_R is deprecated (has been but now CRAN is enforcing it). Investigate switching to pracma::romberg or re-writing the C-code. (spoiler alert: we re-write the C-code) #21

call_R is deprecated (has been but now CRAN is enforcing it). Investigate switching to pracma::romberg or re-writing the C-code. (spoiler alert: we re-write the C-code) #21

Comments

swihart commented Jan 21, 2022

swihart commented Jan 21, 2022

swihart commented Jan 22, 2022

swihart commented Jan 23, 2022 • edited

swihart commented Jan 23, 2022 • edited

swihart commented Jan 23, 2022 • edited

swihart commented Jan 23, 2022

swihart commented Jan 23, 2022 • edited

swihart commented Jan 23, 2022

swihart commented Jan 23, 2022 • edited

swihart commented Jan 23, 2022 • edited

swihart commented Jan 23, 2022 • edited

swihart commented Jan 23, 2022

swihart commented Jan 24, 2022 • edited

swihart commented Jan 25, 2022 • edited

BREAKTHROUGH CITY.

swihart commented Jan 25, 2022

swihart commented Jan 25, 2022 • edited

swihart commented Jan 25, 2022 • edited

swihart commented Jan 25, 2022 • edited

swihart commented Jan 26, 2022 • edited

Oh wow.

I have it running! Giving same answers as gnlmix!

BUT!

swihart commented Jan 26, 2022

swihart commented Jan 26, 2022 • edited

swihart commented Jan 26, 2022

TECMO DANCE!

This one weird trick makes R-developers hate them.

swihart commented Jan 26, 2022

swihart commented Jan 26, 2022 • edited

repeated

rmutil

swihart commented Jan 26, 2022 • edited

swihart commented Feb 9, 2022

swihart commented Jan 23, 2022 •

edited

swihart commented Jan 23, 2022 •

edited

swihart commented Jan 23, 2022 •

edited

swihart commented Jan 23, 2022 •

edited

swihart commented Jan 23, 2022 •

edited

swihart commented Jan 23, 2022 •

edited

swihart commented Jan 23, 2022 •

edited

swihart commented Jan 24, 2022 •

edited

swihart commented Jan 25, 2022 •

edited

swihart commented Jan 25, 2022 •

edited

swihart commented Jan 25, 2022 •

edited

swihart commented Jan 25, 2022 •

edited

swihart commented Jan 26, 2022 •

edited

swihart commented Jan 26, 2022 •

edited

swihart commented Jan 26, 2022 •

edited

swihart commented Jan 26, 2022 •

edited