Revision 28f4ee9e406fce786bb126198f84ce6ec43cb9c1 authored by Adrian Baddeley on 10 February 2011, 00:00:00 UTC, committed by Gabor Csardi on 10 February 2011, 00:00:00 UTC
1 parent c4ff53c
call3d.c
/*
$Revision: 1.5 $ $Date: 2010/10/24 10:57:02 $
R interface
Pass data between R and internally-defined data structures
# /////////////////////////////////////////////
# AUTHOR: Adrian Baddeley, CWI, Amsterdam, 1991.
#
# MODIFIED BY: Adrian Baddeley, Perth 2009
#
# This software is distributed free
# under the conditions that
# (1) it shall not be incorporated
# in software that is subsequently sold
# (2) the authorship of the software shall
# be acknowledged in any publication that
# uses results generated by the software
# (3) this notice shall remain in place
# in each file.
# //////////////////////////////////////////////
*/
#include <R.h>
#include "geom3.h"
#include "functable.h"
#undef DEBUG
#ifdef DEBUG
#define DEBUGMESSAGE(S) Rprintf(S);
#else
#define DEBUGMESSAGE(S)
#endif
void g3one(Point *p, int n, Box *b, Ftable *g);
void g3three(Point *p, int n, Box *b, Ftable *g);
void g3cen(Point *p, int n, Box *b, H4table *count);
void k3trans(Point *p, int n, Box *b, Ftable *k);
void k3isot(Point *p, int n, Box *b, Ftable *k);
void pcf3trans(Point *p, int n, Box *b, Ftable *pcf, double delta);
void pcf3isot(Point *p, int n, Box *b, Ftable *pcf, double delta);
void phatminus(Point *p, int n, Box *b, double vside, Itable *count);
void phatnaive(Point *p, int n, Box *b, double vside, Itable *count);
void p3hat4(Point *p, int n, Box *b, double vside, H4table *count);
/*
ALLOCATION OF SPACE FOR STRUCTURES/ARRAYS
We have defined an alloc() and free() function for each type.
However, the free() functions currently do nothing,
because we use R_alloc to allocate transient space,
which is freed automatically by R.
*/
Ftable *
allocFtable(n) /* allocate function table of size n */
int n;
{
Ftable *x;
x = (Ftable *) R_alloc(1, sizeof(Ftable));
x->n = n;
x->f = (double *) R_alloc(n, sizeof(double));
x->num = (double *) R_alloc(n, sizeof(double));
x->denom = (double *) R_alloc(n, sizeof(double));
return(x);
}
void freeFtable(x) Ftable *x; { }
Itable *
allocItable(n)
int n;
{
Itable *x;
x = (Itable *) R_alloc(1, sizeof(Itable));
x->n = n;
x->num = (int *) R_alloc(n, sizeof(int));
x->denom = (int *) R_alloc(n, sizeof(int));
return(x);
}
void freeItable(x) Itable *x; { }
H4table *
allocH4table(n)
int n;
{
H4table *x;
x = (H4table *) R_alloc(1, sizeof(H4table));
x->n = n;
x->obs = (int *) R_alloc(n, sizeof(int));
x->nco = (int *) R_alloc(n, sizeof(int));
x->cen = (int *) R_alloc(n, sizeof(int));
x->ncc = (int *) R_alloc(n, sizeof(int));
return(x);
}
void freeH4table(x) H4table *x; { }
Box *
allocBox() /* I know this is ridiculous but it's consistent. */
{
Box *b;
b = (Box *) R_alloc(1, sizeof(Box));
return(b);
}
void freeBox(x) Box *x; { }
Point *
allocParray(n) /* allocate array of n Points */
int n;
{
Point *p;
p = (Point *) R_alloc(n, sizeof(Point));
return(p);
}
void freeParray(x) Point *x; { }
/*
CREATE AND INITIALISE DATA STORAGE
*/
Ftable *
MakeFtable(t0, t1, n)
double *t0, *t1;
int *n;
{
Ftable *tab;
int i, nn;
nn = *n;
tab = allocFtable(nn);
tab->t0 = *t0;
tab->t1 = *t1;
for(i = 0; i < nn; i++) {
tab->f[i] = 0.0;
tab->num[i] = 0;
tab->denom[i] = 0;
}
return(tab);
}
Itable *
MakeItable(t0, t1, n)
double *t0, *t1;
int *n;
{
Itable *tab;
int i, nn;
nn = *n;
tab = allocItable(nn);
tab->t0 = *t0;
tab->t1 = *t1;
for(i = 0; i < nn; i++) {
tab->num[i] = 0;
tab->denom[i] = 0;
}
return(tab);
}
H4table *
MakeH4table(t0, t1, n)
double *t0, *t1;
int *n;
{
H4table *tab;
int i, nn;
nn = *n;
tab = allocH4table(nn);
tab->t0 = *t0;
tab->t1 = *t1;
for(i = 0; i < nn; i++) {
tab->obs[i] = 0;
tab->nco[i] = 0;
tab->cen[i] = 0;
tab->ncc[i] = 0;
}
tab->upperobs = 0;
tab->uppercen = 0;
return(tab);
}
/*
CONVERSION OF DATA TYPES
R -> internal
including allocation of internal data types as needed
*/
Point *
RtoPointarray(x,y,z,n)
double *x, *y, *z;
int *n;
{
int i, nn;
Point *p;
nn = *n;
p = allocParray(nn);
for(i = 0; i < nn; i++) {
p[i].x = x[i];
p[i].y = y[i];
p[i].z = z[i];
}
return(p);
}
Box *
RtoBox(x0, x1, y0, y1, z0, z1)
double *x0, *x1, *y0, *y1, *z0, *z1;
{
Box *b;
b = allocBox();
b->x0 = *x0;
b->x1 = *x1;
b->y0 = *y0;
b->y1 = *y1;
b->z0 = *z0;
b->z1 = *z1;
return(b);
}
/*
CONVERSION OF DATA TYPES
internal -> R
Note: it can generally be assumed that the R arguments
are already allocated vectors of correct length,
so we do not allocate them.
*/
void
FtabletoR(tab, t0, t1, n, f, num, denom)
/* internal */
Ftable *tab;
/* R representation */
double *t0, *t1;
int *n;
double *f, *num, *denom;
{
int i;
*t0 = tab->t0;
*t1 = tab->t1;
*n = tab->n;
for(i = 0; i < tab->n; i++) {
f[i] = tab->f[i];
num[i] = tab->num[i];
denom[i] = tab->denom[i];
}
freeFtable(tab);
}
void
ItabletoR(tab, t0, t1, m, num, denom)
/* internal */
Itable *tab;
/* R representation */
double *t0, *t1;
int *m;
int *num, *denom;
{
int i;
*t0 = tab->t0;
*t1 = tab->t1;
*m = tab->n;
for(i = 0; i < tab->n; i++) {
num[i] = tab->num[i];
denom[i] = tab->denom[i];
}
freeItable(tab);
}
void
H4tabletoR(tab, t0, t1, m, obs, nco, cen, ncc, upperobs, uppercen)
/* internal */
H4table *tab;
/* R representation */
double *t0, *t1;
int *m;
int *obs, *nco, *cen, *ncc;
int *upperobs, *uppercen;
{
int i;
*t0 = tab->t0;
*t1 = tab->t1;
*m = tab->n;
*upperobs = tab->upperobs;
*uppercen = tab->uppercen;
for(i = 0; i < tab->n; i++) {
obs[i] = tab->obs[i];
nco[i] = tab->nco[i];
cen[i] = tab->cen[i];
ncc[i] = tab->ncc[i];
}
freeH4table(tab);
}
/*
R CALLING INTERFACE
These routines are called from R by
> .C("routine-name", ....)
*/
void
RcallK3(x,y,z, n, x0, x1, y0, y1, z0, z1, t0, t1, m, f, num, denom, method)
double *x, *y, *z; /* points */
int *n;
double *x0, *x1, /* box */
*y0, *y1,
*z0, *z1;
double *t0, *t1; /* Ftable */
int *m;
double *f, *num, *denom;
int *method;
{
Point *p;
Box *b;
Ftable *tab;
p = RtoPointarray(x, y, z, n);
b = RtoBox(x0, x1, y0, y1, z0, z1);
tab = MakeFtable(t0, t1, m);
switch((int) *method) {
case 0:
k3trans(p, (int) *n, b, tab); break;
case 1:
k3isot(p, (int) *n, b, tab); break;
default:
Rprintf("Method %d not implemented: defaults to 0\n", *method);
k3trans(p, (int) *n, b, tab); break;
}
FtabletoR(tab, t0, t1, m, f, num, denom);
}
void
RcallG3(x,y,z, n, x0, x1, y0, y1, z0, z1, t0, t1, m, f, num, denom, method)
double *x, *y, *z; /* points */
int *n;
double *x0, *x1, /* box */
*y0, *y1,
*z0, *z1;
double *t0, *t1; /* Ftable */
int *m;
double *f, *num, *denom;
int *method;
{
Point *p;
Box *b;
Ftable *tab;
p = RtoPointarray(x, y, z, n);
b = RtoBox(x0, x1, y0, y1, z0, z1);
tab = MakeFtable(t0, t1, m);
switch(*method) {
case 1:
g3one(p, (int) *n, b, tab);
break;
case 3:
g3three(p, (int) *n, b, tab);
break;
default:
Rprintf("Method %d not implemented: defaults to 3\n", *method);
g3three(p, (int) *n, b, tab);
}
FtabletoR(tab, t0, t1, m, f, num, denom);
}
void
RcallG3cen(x,y,z, n, x0, x1, y0, y1, z0, z1,
t0, t1, m, obs, nco, cen, ncc, upperobs, uppercen)
double *x, *y, *z; /* points */
int *n;
double *x0, *x1, /* box */
*y0, *y1,
*z0, *z1;
double *t0, *t1;
int *m; /* H4table */
int *obs, *nco, *cen, *ncc;
int *upperobs, *uppercen;
{
Point *p;
Box *b;
H4table *count;
DEBUGMESSAGE("Inside RcallG3cen\n")
p = RtoPointarray(x, y, z, n);
b = RtoBox(x0, x1, y0, y1, z0, z1);
count = MakeH4table(t0, t1, m);
g3cen(p, (int) *n, b, count);
H4tabletoR(count, t0, t1, m, obs, nco, cen, ncc, upperobs, uppercen);
DEBUGMESSAGE("Leaving RcallG3cen\n")
}
void
RcallF3(x,y,z, n, x0, x1, y0, y1, z0, z1,
vside,
t0, t1, m, num, denom, method)
double *x, *y, *z; /* points */
int *n;
double *x0, *x1, /* box */
*y0, *y1,
*z0, *z1;
double *vside;
double *t0, *t1;
int *m; /* Itable */
int *num, *denom;
int *method;
{
Point *p;
Box *b;
Itable *count;
DEBUGMESSAGE("Inside Rcall_f3\n")
p = RtoPointarray(x, y, z, n);
b = RtoBox(x0, x1, y0, y1, z0, z1);
count = MakeItable(t0, t1, m);
switch((int) *method) {
case 0:
phatnaive(p, (int) *n, b, *vside, count);
break;
case 1:
phatminus(p, (int) *n, b, *vside, count);
break;
default:
Rprintf("Method %d not recognised: defaults to 1\n", *method);
phatminus(p, (int) *n, b, *vside, count);
}
ItabletoR(count, t0, t1, m, num, denom);
DEBUGMESSAGE("Leaving Rcall_f3\n")
}
void
RcallF3cen(x,y,z, n, x0, x1, y0, y1, z0, z1,
vside,
t0, t1, m, obs, nco, cen, ncc, upperobs, uppercen)
double *x, *y, *z; /* points */
int *n;
double *x0, *x1, /* box */
*y0, *y1,
*z0, *z1;
double *vside;
double *t0, *t1;
int *m; /* H4table */
int *obs, *nco, *cen, *ncc;
int *upperobs, *uppercen;
{
Point *p;
Box *b;
H4table *count;
DEBUGMESSAGE("Inside Rcallf3cen\n")
p = RtoPointarray(x, y, z, n);
b = RtoBox(x0, x1, y0, y1, z0, z1);
count = MakeH4table(t0, t1, m);
p3hat4(p, (int) *n, b, *vside, count);
H4tabletoR(count, t0, t1, m, obs, nco, cen, ncc, upperobs, uppercen);
DEBUGMESSAGE("Leaving Rcallf3cen\n")
}
void
Rcallpcf3(x,y,z, n, x0, x1, y0, y1, z0, z1, t0, t1, m, f, num, denom, method,
delta)
double *x, *y, *z; /* points */
int *n;
double *x0, *x1, /* box */
*y0, *y1,
*z0, *z1;
double *t0, *t1; /* Ftable */
int *m;
double *f, *num, *denom;
int *method;
double *delta; /* Epanechnikov kernel halfwidth */
{
Point *p;
Box *b;
Ftable *tab;
p = RtoPointarray(x, y, z, n);
b = RtoBox(x0, x1, y0, y1, z0, z1);
tab = MakeFtable(t0, t1, m);
switch((int) *method) {
case 0:
pcf3trans(p, (int) *n, b, tab, (double) *delta); break;
case 1:
pcf3isot(p, (int) *n, b, tab, (double) *delta); break;
default:
Rprintf("Method %d not implemented: defaults to 0\n", *method);
pcf3trans(p, (int) *n, b, tab, (double) *delta); break;
}
FtabletoR(tab, t0, t1, m, f, num, denom);
}
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...