[spcopula-commits] r102 - / pkg pkg/R pkg/man pkg/src

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Sun Sep 1 13:10:27 CEST 2013


Author: mappe
Date: 2013-09-01 13:10:27 +0200 (Sun, 01 Sep 2013)
New Revision: 102

Added:
   pkg/man/getNeighbours.experimental.Rd
   pkg/src/
   pkg/src/Rinterface.h
   pkg/src/spatialPreparation.c
   thoughts-Marius/
Modified:
   pkg/NAMESPACE
   pkg/R/spatialPreparation.R
Log:
- added thoughts-Marius folder
- added C implementation of getNeighbours() in getNeighbours.experimental()

Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE	2013-08-01 10:31:15 UTC (rev 101)
+++ pkg/NAMESPACE	2013-09-01 11:10:27 UTC (rev 102)
@@ -27,7 +27,7 @@
 # export(setSizeLim)
 
 # spatial
-export(getNeighbours,getStNeighbours)
+export(getNeighbours,getNeighbours.experimental, getStNeighbours)
 export(calcBins)
 export(calcSpTreeDists, dropSpTree)
 

Modified: pkg/R/spatialPreparation.R
===================================================================
--- pkg/R/spatialPreparation.R	2013-08-01 10:31:15 UTC (rev 101)
+++ pkg/R/spatialPreparation.R	2013-09-01 11:10:27 UTC (rev 102)
@@ -98,6 +98,60 @@
                        predLocs, prediction, var))
 }
 
+
+
+getNeighbours.experimental = function (dataLocs, predLocs, var = names(dataLocs)[1], size = 5, 
+									   prediction = FALSE, min.dist = 0.01) 
+{
+  stopifnot((!prediction && missing(predLocs)) || (prediction && 
+                                                     !missing(predLocs)))
+  stopifnot(min.dist > 0 || prediction)
+  if (missing(predLocs) && !prediction) 
+    predLocs = dataLocs
+  stopifnot(is(predLocs, "Spatial"))
+  if (any(is.na(match(var, names(dataLocs))))) 
+    stop("At least one of the variables is unkown or is not part of the data.")
+  nLocs <- length(predLocs)
+  size <- min(size, length(dataLocs) + prediction)
+  
+  allLocs <- matrix(0, nLocs, size)
+  allDists <- matrix(0, nLocs, size - 1)
+  result = .C("getNeighbours_2d", 
+     allLocs = as.double(allLocs),
+     allDists = as.double(allDists),
+     as.double(coordinates(dataLocs)),
+     as.double(coordinates(predLocs)),
+     as.integer(length(dataLocs)),
+     as.integer(length(predLocs)),
+     as.double(min.dist),
+     as.integer(size), 
+     as.integer(prediction),
+     PACKAGE="spcopula")
+  
+  if (!prediction) allData <- matrix(dataLocs[result$allLocs, var, drop = F]@data[[1]], nLocs, size)
+  else allData <- matrix(c(rep(NA,nLocs),dataLocs[result$allLocs[(nLocs+1):(nLocs*size)], var, drop = F]@data[[1]]), nLocs, size)
+  
+  if (!prediction) 
+    predLocs <- NULL
+  
+  colnames(allData) <- paste(paste("N", rep(0:(size - 1), each = length(var)), 
+                                   sep = ""), rep(var, size), sep = ".")
+  
+  
+  return(neighbourhood(allData, matrix(result$allDists,nLocs, size - 1), matrix(result$allLocs, nLocs, size), dataLocs, 
+                        predLocs, prediction, var))
+}
+
+
+
+
+
+
+
+
+
+
+
 #############
 ## BINNING ##
 #############

Added: pkg/man/getNeighbours.experimental.Rd
===================================================================
--- pkg/man/getNeighbours.experimental.Rd	                        (rev 0)
+++ pkg/man/getNeighbours.experimental.Rd	2013-09-01 11:10:27 UTC (rev 102)
@@ -0,0 +1,47 @@
+\name{getNeighbours.experimental}
+\alias{getNeighbours.experimental}
+
+\title{
+Creating Local Neighbourhoods
+}
+\description{
+This function calculates a local neighbourhood to be used for fitting of spatial/spatio-temporal vine copulas and for prediction using spatial/spatio-temporal vine copulas.
+}
+\usage{
+getNeighbours.experimental(dataLocs, predLocs, var = names(dataLocs)[1], size = 5, prediction=FALSE, min.dist = 0.01)
+}
+\arguments{
+  \item{dataLocs}{
+some spatial data frame holding the data used for estimation/prediction
+}
+  \item{predLocs}{
+A spatial object defining the prediction locations, might be missing if the neighbourhood is used for fitting.
+}
+  \item{var}{
+the variable name of interest, by default the first variable is used.
+}
+  \item{size}{
+The size of the neighbourhood including the location of interest (for fitting as well for prediction).
+}
+  \item{prediction}{whether the neighbourhood should be used for prediction (TRUE) or spatial/Spatio-temporal vine copula fitting.}
+  \item{min.dist}{
+the minimal distance for a location to be included. Must be larger than 0 for fitting purposes and might be 0 for prediction.
+}
+}
+\value{
+An object of \code{\linkS4class{neighbourhood}}.
+}
+\author{
+Benedikt Graeler
+}
+
+\seealso{
+See \code{\link{neighbourhood}} for the native constructor of a \code{\linkS4class{neighbourhood}} class.
+}
+\examples{
+spdf <- data.frame(x=c(112,154,212,289,345),y=c(124,198,85,168,346),measure=rlnorm(5))
+coordinates(spdf) <- ~x+y
+
+getNeighbours.experimental(spdf,size=4)
+}
+\keyword{ spatial }
\ No newline at end of file

Added: pkg/src/Rinterface.h
===================================================================
--- pkg/src/Rinterface.h	                        (rev 0)
+++ pkg/src/Rinterface.h	2013-09-01 11:10:27 UTC (rev 102)
@@ -0,0 +1,33 @@
+
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+
+
+#ifdef _WIN32
+#define R_EXPORT __declspec(dllexport)
+#else
+#define R_EXPORT
+#endif
+
+
+
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/* FORWARD DECLARATIONS */
+
+void R_EXPORT distMatrix(double *out, double *xa, double *xb, int *na, int*nb, int *d);
+void R_EXPORT distMatrix_2d(double *out, double *xa, double *xb, int *na, int*nb);
+void R_EXPORT distMatrixInner(double *out, double *xa, int *na, int *d);
+void R_EXPORT distMatrixInner_2d(double *out, double *xa, int *na);
+
+void R_EXPORT getNeighbours_2d(double *out_locs, double *out_dists, double *data_locs, double *pred_locs, int *n_data, int *n_pred, double *min_dist, int *k, int *prediction);
+
+
+#ifdef __cplusplus
+}
+#endif
\ No newline at end of file

Added: pkg/src/spatialPreparation.c
===================================================================
--- pkg/src/spatialPreparation.c	                        (rev 0)
+++ pkg/src/spatialPreparation.c	2013-09-01 11:10:27 UTC (rev 102)
@@ -0,0 +1,244 @@
+
+#include "Rinterface.h"
+
+
+
+
+
+/* Helper Functions */
+inline double dist_2d(double ax, double ay, double bx, double by) {
+	return(sqrt( (ax-bx)*(ax-bx) + (ay-by)*(ay-by) ));
+}
+
+
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+
+
+
+
+
+/**
+* Computes the distance matrix of two sets of n-dimensional points
+* @param out result (na,nb) matrix (mast be allocated in R or in calling function before)
+* @param xa coordinates of the 1st set of points, each row of this matrix is a point
+* @param xb coordinates of the 2nd set of points, each row of this matrix is a point
+* @param na number of points in xa
+* @param nb number of points in xb
+* @param d  number of dimensions, equals the number of columns in xa and xb
+**/ 
+void R_EXPORT distMatrix(double *out, double *xa, double *xb, int *na, int*nb, int *d) {
+	if (*d == 2) {
+		distMatrix_2d(out, xa, xb, na, nb);
+	}
+	else
+	{
+		for (int i=0;i<*na; ++i) {
+			for (int j=0;j<*nb; ++j) {
+				out[j * (*na) + i] = 0;
+				
+				for (int k=0; k<*d; ++k) {
+					out[j * (*na) + i] += (xa[k * (*na) + i] -  xb[k * (*nb) + j]) * (xa[k * (*na) + i] -  xb[k * (*nb) + j]); // In R: col major
+					
+				}
+				out[j * (*na) + i] = sqrt(out[j * (*na) + i]);		
+			}
+		}
+	}
+}
+
+
+
+/**
+* Computes the distance matrix of two sets of 2-dimensional points
+* @param out result (na,nb) matrix (mast be allocated in R or in calling function before)
+* @param xa coordinates of the 1st set of points, each row of this matrix is a point
+* @param xb coordinates of the 2nd set of points, each row of this matrix is a point
+* @param na number of points in xa
+* @param nb number of points in xb
+**/ 
+void R_EXPORT distMatrix_2d(double *out, double *xa, double *xb, int *na, int*nb) {
+	for (int i=0;i<*na; ++i) {
+		for (int j=0;j<*nb; ++j) {
+			out[j * (*na) + i] =  dist_2d(xa[i],xa[*na + i], xb[j], xb[*nb + j]);	
+		}
+	}
+}
+
+
+
+/**
+* Computes the "inner" distance matrix of one set of n-dimensional points
+* @param out result (na,na) matrix (mast be allocated in R or in calling function before)
+* @param xa coordinates of the 1st set of points, each row of this matrix is a point
+* @param na number of points in xa
+* @param d  number of dimensions, equals the number of columns in xa and xb
+**/ 
+void R_EXPORT distMatrixInner(double *out, double *xa, int *na, int *d) {
+	if (*d == 2) {
+		distMatrixInner_2d(out, xa, na);
+	}
+	else
+	{
+		for (int i=0;i<*na; ++i) {
+			for (int j=i+1;j<*na; ++j) {
+				out[j * (*na) + i] = 0;
+				for (int k=0; k<*d; ++k) {
+					out[j * (*na) + i] += (xa[k * (*na) + i] -  xa[k * (*na) + j]) * (xa[k * (*na) + i] -  xa[k * (*na) + j]); // In R: col major
+					
+				}
+				out[j * (*na) + i] = sqrt(out[j * (*na) + i]);	
+				out[i * (*na) + j] = out[j * (*na) + i]; 
+			}
+		}
+		for (int i=0;i<*na; ++i) {
+			out[i * (*na) + i] = 0.0;
+		}
+	}
+}
+
+
+
+/**
+* Computes the "inner" distance matrix of one set of n-dimensional points
+* @param out result (na,na) matrix (mast be allocated in R or in calling function before)
+* @param xa coordinates of the 1st set of points, each row of this matrix is a point
+* @param na number of points in xa
+**/ 
+void R_EXPORT distMatrixInner_2d(double *out, double *xa, int *na) {
+	for (int i=0;i<*na; ++i) {
+		for (int j=i+1;j<*na; ++j) {
+			out[j * (*na) + i] =  dist_2d(xa[i],xa[*na + i], xa[j], xa[*na + j]);	
+			out[i * (*na) + j] = out[j * (*na) + i]; 
+		}
+	}
+	// set diagonal to zero
+	for (int i=0;i<*na; ++i) {
+		out[i * (*na) + i] = 0.0;
+	}
+}
+
+
+
+
+
+
+/*for (i in 1:nLocs) {
+    tempDists <- spDists(dataLocs, predLocs[i, ])
+    tempDists[tempDists < min.dist] <- Inf
+    spLocs <- order(tempDists)[1:(size - 1)]
+    
+    allLocs[i,] <- c(i, spLocs)
+    allDists[i,] <- tempDists[spLocs]
+    allData[i,(prediction+1):size] <- dataLocs[c(i[!prediction],spLocs),
+                                               var, drop = F]@data[[1]]
+  }
+  */
+
+/**
+*
+*
+*
+*
+* @param k size of the neighbourhood including the location of interest
+*/
+void R_EXPORT getNeighbours_2d(double *out_locs, double *out_dists, double *data_locs, double *pred_locs, int *n_data, int *n_pred, double *min_dist, int *k, int *prediction) {
+	double *C = (double*)malloc((*n_data) *  (*n_pred) * sizeof(double));
+	double *knn_dists = (double*)malloc((*k-1) * sizeof(double));
+	int *knn_idxs  = (int*)malloc((*k-1) * sizeof(int));
+	
+	// Depending on whether pred_locs are the same as data_locs, use optimized "inner" distance matrix calculation
+	if (!prediction) { 
+		// Compute inner distance matrix and store in C
+		distMatrixInner_2d(C, data_locs, n_data);
+	}
+	else {
+	
+		// Compute distance matrix and store in C
+		distMatrix_2d(C, data_locs, pred_locs, n_data, n_pred);
+	}	
+	
+	// for each pred location
+	for (int i=0; i<*n_pred; ++i) {
+		knn_dists[0] = INFINITY;
+		knn_idxs[0] = -1;
+		
+		for (int j=0; j<*n_data; ++j) {
+			if (C[i* (*n_data) + j] < *min_dist) continue;
+			int nbrank = 0;
+			while (C[i* (*n_data) + j] > knn_dists[nbrank] &&  nbrank < (*k-1)) ++nbrank; 
+			if (nbrank < (*k-1)) {		
+				for (int s=(*k-3); s >= nbrank; --s) {
+					knn_dists[s + 1] = knn_dists[s]; // shift array if new element in knn array
+					knn_idxs[s + 1] = knn_idxs[s];
+				}
+				knn_dists[nbrank] = C[i* (*n_data) + j];
+				knn_idxs[nbrank] = j+1; // Indices in R start with 1
+			}			
+		}
+		
+		out_locs[i] = i+1; // Indices in R start with 1
+		for (int j=1; j<(*k); ++j) {			
+			// Assign c(i, knn_idxs[1:size-1]) as i-th row in out_locs
+			out_locs[j*(*n_pred) + i] = knn_idxs[j-1];			
+			// Assign knn_dists[1:size-1] as i-th row in out_dists
+			out_dists[(j-1)*(*n_pred) + i] = knn_dists[j-1];
+		}
+		
+	}
+	
+	free(knn_idxs);
+	free(knn_dists);
+	free(C);
+}
+
+
+
+
+
+
+/*
+calcSpLagInd <- function(data, boundaries) {
+  lags <- vector("list",length(boundaries))
+  
+  dists <- spDists(data)
+  nlocs <- length(data)
+  
+  for (i in 1:(nlocs-1)) {
+    for (j in (i+1):nlocs) {
+      d <- dists[i,j]
+      for ( k in 1:length(boundaries)) {
+        if (d < boundaries[k]) {
+          lags[[k]] <- rbind(lags[[k]],c(i,j,d))
+          break()
+        }
+      }
+    }
+  }
+  return(lags)
+}
+*/
+
+/**
+*
+* @
+*
+*
+*
+* OpenMP?
+*/
+// SEXP calcSpLagInd(double *o double *idata, double *iboundaries, int *in_boundaries) {
+
+
+	// for
+
+// }
+
+
+
+#ifdef __cplusplus
+}
+#endif
\ No newline at end of file



More information about the spcopula-commits mailing list