Revision e994a4415e67fa60cbfd3f208aaab20872521c0b authored by Martin Schlather on 14 February 2019, 21:02:19 UTC, committed by cran-robot on 14 February 2019, 21:02:19 UTC
1 parent e5a7a2f
bigdata.cc
/*
Authors
Martin Schlather, schlather@math.uni-mannheim.de
library for simulation of random fields
Copyright (C) 2001 -- 2017 Martin Schlather,
This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 3
of the License, or (at your option) any later version.
RO
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
*/
#include <R.h>
//#include <Rdefines.h>
//#include <R_ext/Linpack.h>
//#include <stdio.h>
//#include <stdlib.h>
//#include <unistd.h>
//#include <string.h>
#include "RF.h"
SEXP countelements(SEXP Idx, SEXP N, SEXP Totparts) {
int i, *boxes,
*idx = INTEGER(Idx),
nbox = INTEGER(Totparts)[0],
n = INTEGER(N)[0];
SEXP Boxes;
PROTECT(Boxes = allocVector(INTSXP, nbox));
boxes = INTEGER(Boxes);
for (i=0; i<nbox; i++) boxes[i]= 0;
for (i=0; i<n; i++) {
boxes[idx[i]]++;
}
UNPROTECT(1);
return Boxes;
}
SEXP countneighbours(SEXP Xdim, SEXP Parts, SEXP Squarelength,
SEXP Cumgridlen, SEXP Boxes, SEXP MAXN) {
int d, totcumlen, relstart, y, e,
dim = INTEGER(Xdim)[0],
*nb = (int*) MALLOC(sizeof(int) * dim),
*loc = (int*) MALLOC(sizeof(int) * dim),
x = 0,
sum = 0,
sl = INTEGER(Squarelength)[0],
boundary = (sl - 1) / 2,
// total = cumgridlen[dim],
maxn = INTEGER(MAXN)[0],
*parts = INTEGER(Parts),
*cumgridlen = INTEGER(Cumgridlen),
nboxes = length(Boxes),
*boxes = INTEGER(Boxes);
SEXP Neighbours;
PROTECT(Neighbours = allocVector(INTSXP, nboxes));
int *neighbours = INTEGER(Neighbours);
totcumlen = 0;
for (d=0; d<dim; d++) {
loc[d] = -boundary;
nb[d] = 0;
totcumlen += cumgridlen[d];
//print("%d %d %d\n", d, totcumlen, cumgridlen[d]);
}
relstart = totcumlen * boundary;
d = 0;
while(d < dim) {
y = x - relstart;
// print("%d\n", x);
neighbours[x] = 0;
sum = e = 0;
while(e < dim) {
bool inside = true;
int j;
for (j=0; j<dim; j++) {
double abs = loc[j] + nb[j];
if (abs < 0 || abs>=parts[j]) {inside=false; break;}
}
if (inside) {
// print("inside %d (%d, %d)\n", y, loc[0], loc[1]);
sum += boxes[y];
neighbours[x]++;
}
e = 0;
loc[e]++;
y++;
while (loc[e] > boundary) {
loc[e] = -boundary;
y -= cumgridlen[e] * sl;
if (++e >= dim) break;
loc[e]++;
y += cumgridlen[e];
}
//print("e=%d\n", e);
}
// print("sum=%d maxn=%d (%d %d) %d %d parts=%d %d; %d %dl nei=%d\n",
// sum, maxn, nb[0], nb[1], totcumlen, relstart, parts[0], parts[1],
// cumgridlen[0], cumgridlen[1], neighbours[x]);
// assert(false);
if (sum > maxn) {
Neighbours = NILSXP;
goto ErrorHandling;
}
d = 0;
nb[d]++;
x++;
while (nb[d] >= parts[d]) {
nb[d] = 0;
if (++d >= dim) break;
nb[d]++;
}
// print("d=%d (%d %d %d) [%d %d %d]\n", d, nb[0], nb[1], nb[2],
// parts[0], parts[1], parts[2]);
// assert(false);
}
ErrorHandling:
UNPROTECT(1);
FREE(nb);
FREE(loc);
return Neighbours;
}
SEXP getelements(SEXP Idx, SEXP Xdim, SEXP N, SEXP Cumgridlen, SEXP Boxes) {
int i, err = NOERROR,
*idx = INTEGER(Idx),
dim = INTEGER(Xdim)[0],
n = INTEGER(N)[0],
*cumgridlen = INTEGER(Cumgridlen),
*boxes = INTEGER(Boxes),
total = cumgridlen[dim],
*count = NULL,
**elm = NULL;
SEXP subel = R_NilValue;
if ((elm = (int **) MALLOC(sizeof(int*) * total)) == NULL ||
(count = (int*) MALLOC(sizeof(int) * total)) == NULL) {
err = ERRORMEMORYALLOCATION;
goto ErrorHandling;
}
for (i=0; i<total; i++) elm[i] = NULL;
for (i=0; i<total; i++) {
// print("%d %d l=%d \n", i, total, boxes[i]);
if ((elm[i] = (int *) MALLOC(sizeof(int) * boxes[i])) == NULL) {
err = ERRORMEMORYALLOCATION;
goto ErrorHandling;
}
count[i] = 0;
}
for (i=0; i<n; i++) {
int k=idx[i];
elm[k][(count[k])++] = i + 1;
}
PROTECT(subel = allocVector(VECSXP, total));
for (i=0; i<total; i++) {
// print("%d %d %d\n", i, total, boxes[i]);
SEXP el;
int k,
end = boxes[i];
PROTECT(el = allocVector(INTSXP, boxes[i]));
for (k=0; k<end; k++) INTEGER(el)[k] = elm[i][k];
SET_VECTOR_ELT(subel, i, el);
UNPROTECT(1);
}
UNPROTECT(1);
ErrorHandling :
if (elm != NULL) {
for (i=0; i<total; i++) FREE(elm[i]);
UNCONDFREE(elm);
}
FREE(count);
if (err!=NOERROR) XERR(err);
return subel;
}
SEXP getneighbours(SEXP Xdim, SEXP Parts, SEXP Squarelength,
SEXP Cumgridlen, SEXP Neighbours) {
int i, d, sum, totcumlen, relstart, x, y, e,
err = NOERROR,
dim = INTEGER(Xdim)[0],
*nb = (int*) MALLOC(sizeof(int) * dim),
*loc = (int*) MALLOC(sizeof(int) * dim),
*parts = INTEGER(Parts),
sl = INTEGER(Squarelength)[0],
*cumgridlen = INTEGER(Cumgridlen),
*neighbours = INTEGER(Neighbours),
boundary = (sl - 1) / 2,
total = cumgridlen[dim],
** neighb = NULL;
SEXP subnei = R_NilValue;
if ( (neighb = (int **) MALLOC(sizeof(int*) * total)) == NULL) {
err = ERRORMEMORYALLOCATION;
goto ErrorHandling;
}
// print("OK %d\n", total);
for (i=0; i<total; i++) neighb[i] = NULL;
for (i=0; i<total; i++) {
// print("%d %d l=%d %d\n", i, total, neighbours[i]);
// R_CheckUserInterrupt();
if ((neighb[i] = (int *) MALLOC(sizeof(int) * neighbours[i])) == NULL) {
err = ERRORMEMORYALLOCATION;
goto ErrorHandling;
}
}
x = 0;
totcumlen = 0;
for (d=0; d<dim; d++) {
loc[d] = -boundary;
nb[d] = 0;
totcumlen += cumgridlen[d];
}
relstart = totcumlen * boundary;
d = 0;
while(d < dim) {
y = x - relstart;
sum = e = 0;
while(e < dim) {
bool inside = true;
int j;
for (j=0; j<dim; j++) {
double abs = loc[j] + nb[j];
if (abs < 0 || abs>=parts[j]) {inside=false; break;}
}
if (inside) {
neighb[x][sum] = y + 1;
sum++;
}
e = 0;
loc[e]++;
y++;
while (loc[e] > boundary) {
loc[e] = -boundary;
y -= cumgridlen[e] * sl;
if (++e >= dim) break;
loc[e]++;
y += cumgridlen[e];
}
}
d = 0;
nb[d]++;
x++;
while (nb[d] >= parts[d]) {
nb[d] = 0;
if (++d >= dim) break;
nb[d]++;
}
}
// print("heres\n");
PROTECT(subnei = allocVector(VECSXP, total));
for (i=0; i<total; i++) {
// print("%d %d %d %d\n", i, total, neighbours[i]);
R_CheckUserInterrupt();
SEXP nei;
PROTECT(nei = allocVector(INTSXP, neighbours[i]));
int k,
end = neighbours[i];
for (k=0; k<end; k++) INTEGER(nei)[k] = neighb[i][k];
SET_VECTOR_ELT(subnei, i, nei);
UNPROTECT(1);
}
UNPROTECT(1);
ErrorHandling :
FREE(loc);
FREE(nb);
if (neighb != NULL) {
for (i=0; i<total; i++) FREE(neighb[i]);
UNCONDFREE(neighb);
}
if (err!=NOERROR) XERR(err);
return subnei;
}
Computing file changes ...