[adegenet-commits] r96 - in pkg: R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Apr 15 13:57:51 CEST 2008
Author: jombart
Date: 2008-04-15 13:57:51 +0200 (Tue, 15 Apr 2008)
New Revision: 96
Modified:
pkg/R/monmonier.R
pkg/src/monmonier-utils.c
Log:
Fixed the problems occuring with coords from a regular grid.
C code reviewed and tuned (debugging messages available but commented out).
Pkg not checked.
Still needs to fix optimize.monmonier.
Modified: pkg/R/monmonier.R
===================================================================
--- pkg/R/monmonier.R 2008-04-11 15:54:12 UTC (rev 95)
+++ pkg/R/monmonier.R 2008-04-15 11:57:51 UTC (rev 96)
@@ -27,9 +27,8 @@
}
## PRECISION of the xy coordinates (in digits)
-## used to retrieve edges from their coordinates
-## do not go over 10
-PRECISION=8
+## used when coordinates are inputed in C code.
+PRECISION=12
## conversion of the connection network
cn.nb <- cn
@@ -56,7 +55,7 @@
if(scanthres){
- plot(sort(unique(D)[unique(D) > 0],dec=TRUE),main="Local distances plot",type="l",xlab="rank")
+ plot(sort(unique(D)[unique(D) > 0],dec=TRUE),main="Local distances plot",type="l",xlab="rank",ylab="Sorted local distances")
abline(h=Dlim,lty=2)
mtext("Dashed line indicates present threshold")
cat("Indicate the threshold (\'d\' for default): ")
@@ -138,8 +137,8 @@
## A,B,C,D are not coordinates, but indices of vertices
## -AB is the edge whose middle is M ; given to avoid checking MN vs AB
## -CD is the edge whose middle is N ; given to avoid checking MN vs CD
-##
-checkNext <- function(M,N,A,B,C,D,segMat=allSeg){
+## - toRemove: a vector of characters naming other edges in segMat to be removed before checking
+checkNext <- function(M,N,A,B,C,D,segMat=allSeg,toRemove=NULL){
## orientation of the segment MN
xmin <- min(M[1],N[1])
xmax <- max(M[1],N[1])
@@ -167,6 +166,13 @@
subsetSeg <- subsetSeg[!(subsetSeg[,2] < ymin & subsetSeg[,4] < ymin),,drop=FALSE]
subsetSeg <- subsetSeg[!(subsetSeg[,2] > ymax & subsetSeg[,4] > ymax),,drop=FALSE]
+ ## handle toRemove here
+ if(!is.null(toRemove)){
+ idx <- match(toRemove, rownames(subsetSeg))
+ idx <- idx[!is.na(idx)]
+ if(length(idx)>0) { subsetSeg <- subsetSeg[-idx,,drop=FALSE] }
+ }
+
## temp is used for the output of CheckAllSeg
## initialized at 10, which is never returned by CheckAllSeg
temp <- as.integer(10)
@@ -184,6 +190,9 @@
## =======================
##
## restore the matrix type if there is only one segment to check for crossing in subsetSeg
+
+ ## round down coordinates in subsetSeg
+ subsetSeg <- round(subsetSeg, digits=PRECISION)
if(nrow(subsetSeg)>0) {
temp <- .C("CheckAllSeg",as.integer(nrow(subsetSeg)),as.integer(ncol(subsetSeg)),
@@ -194,11 +203,11 @@
if(DEBUG) {
if(temp==1) cat("\n can't go there (code",temp,")") else cat("\n new segment ok (code",temp,")")
}
- ## if a code 1 is returned, CheckAllSeg returns FALSE
+ ## if a code 1 or 3 is returned, CheckAllSeg returns FALSE
## else it returns TRUE
## additional control used (code 10, CheckAllSeg failed)
if(temp==10) stop("CheckAllSeg failure (returned value=10, i.e. unchanged, not computed). Please report the bug.")
- if(temp==1) return(FALSE) else return(TRUE)
+ if(temp==1 | temp==2 | temp==3) return(FALSE) else return(TRUE)
} # end of checkNext
## result object is created
@@ -288,6 +297,7 @@
s1 <- 1
## update existing segments
allSeg <- rbind(allSeg,c(result[[run]]$dir1[[i1-1]]$M,result[[run]]$dir1[[i1]]$M) )
+ rownames(allSeg) <- c(rownames(allSeg)[-nrow(allSeg)] , paste("dir1",i1-1,i1,sep="-"))
## add 1 to the boundary length
current.bd.length <- current.bd.length + 1
hasExpanded <- TRUE
@@ -322,6 +332,7 @@
s2 <- 1
## update existing segments
allSeg <- rbind(allSeg,c(result[[run]]$dir2[[i2-1]]$M,result[[run]]$dir2[[i2]]$M) )
+ rownames(allSeg) <- c(rownames(allSeg)[-nrow(allSeg)] , paste("dir2",i2-1,i2,sep="-"))
## add 1 to the boundary length
current.bd.length <- current.bd.length + 1
hasExpanded <- TRUE
@@ -341,12 +352,17 @@
## handle the looping of a boundary
if(hasExpanded && (current.bd.length>3) && allowLoop){
## check if the two ends of the boundary can be joined
+ ## remove segments ending each direction (to avoid messy code 2 in checkNext)
+ terminalEdges <- c(paste("dir1",i1-1,i1,sep="-"),
+ paste("dir2",i2-1,i2,sep="-"))
+
canLoop <- checkNext(result[[run]]$dir1[[length(result[[run]]$dir1)]]$M,
result[[run]]$dir2[[length(result[[run]]$dir2)]]$M,
result[[run]]$dir1[[length(result[[run]]$dir1)]]$A[1],
result[[run]]$dir1[[length(result[[run]]$dir1)]]$B[1],
result[[run]]$dir2[[length(result[[run]]$dir2)]]$A[1],
- result[[run]]$dir2[[length(result[[run]]$dir2)]]$B[1]
+ result[[run]]$dir2[[length(result[[run]]$dir2)]]$B[1],
+ toRemove=terminalEdges
)
if(canLoop) {
@@ -465,8 +481,8 @@
cex.bwd.1 <- cex.bwd.1/max(cex.bwd.1)
cex.bwd.2 <- cex.bwd.2/max(cex.bwd.2)
## amplify the differences
- cex.bwd.1 <- cex.bwd.1^1.5
- cex.bwd.2 <- cex.bwd.2^1.5
+ ## cex.bwd.1 <- cex.bwd.1^1.5
+ ## cex.bwd.2 <- cex.bwd.2^1.5
if(add.arrows) {
@@ -582,7 +598,7 @@
if(is.null(threshold) || (threshold<0)) {Dlim <- summary(unique(D[D>0]))[5]} else {Dlim <- threshold}
if(scanthres){
- plot(sort(unique(D)[unique(D) > 0],dec=TRUE),main="Local distances plot", type="l", xlab="rank")
+ plot(sort(unique(D)[unique(D) > 0],dec=TRUE),main="Local distances plot", type="l", xlab="rank",ylab="Sorted local distances")
abline(h=Dlim,lty=2)
mtext("Dashed line indicates present threshold")
cat("Indicate the threshold (\'d\' for default): ")
Modified: pkg/src/monmonier-utils.c
===================================================================
--- pkg/src/monmonier-utils.c 2008-04-11 15:54:12 UTC (rev 95)
+++ pkg/src/monmonier-utils.c 2008-04-15 11:57:51 UTC (rev 96)
@@ -53,9 +53,25 @@
int Collinear( tPointd a, tPointd b, tPointd c );
int AreaSign( tPointd a, tPointd b, tPointd c );
void CheckAllSeg(int *nrow, int *ncol, double *tab, tPointd a, tPointd b, int *answer);
+double dAbs(double a);
+int dEqual(double a, double b);
/*-------------------------------------------------------------------*/
+/* dAbs returns the absolute value of a double */
+double dAbs(double a)
+{
+ if(a>=0.0) return a;
+ else return -a;
+}
+/* dEqual returns 1 if two doubles are equal, and 0 otherwise */
+int dEqual(double a, double b)
+{
+ if(dAbs(a-b) < NEARZERO) return 1;
+ else return 0;
+}
+
+
/*-------------------------------------------------------------------
CheckAllSeg: performs the intersection test of a segment with
a given set of segments. Calls SegSeg to perform 2-segments tests.
@@ -91,7 +107,7 @@
We stop as soon as one segment is crossed */
temp = 0;
i = 1;
- while(temp!=1 && i<=n){
+ while(temp ==0 && i<=n){
c[X] = mat[i][1];
c[Y] = mat[i][2];
d[X] = mat[i][3];
@@ -123,7 +139,10 @@
double num, denom; /* Numerator and denoninator of equations. */
int code = 10; /* returned value, default 10 is a failure */
- /* Initialisation of the interesection point 'p' */
+ /* For debugging
+ printf("\n!!! SegSeg: code initialized at %d\n",code);*/
+
+ /* Initialization of the intersection point 'p' */
tPointd p;
p[X] = -1;
@@ -133,44 +152,63 @@
b[X] * (double)( c[Y] - d[Y] ) +
d[X] * (double)( b[Y] - a[Y] ) +
c[X] * (double)( a[Y] - b[Y] );
- /* If denom is zero, then segments are parallel: handle separately. Beware to avoid ... == 0 with doubles */
- if ((denom < NEARZERO) && (denom > -NEARZERO)) code = Parallel(a, b, c, d, p);
+ /* If denom is zero, then segments are parallel: handle separately.
+ Beware to avoid ... == 0 with doubles,
+ as well as ...==... */
+ if (dAbs(denom) < NEARZERO) {
+ code = Parallel(a, b, c, d, p);
+ /* For debugging
+ printf("\n!!! SegSeg: call to Parallel (denom=%f)\n",denom);*/
+ }
else{
- num = a[X] * (double)( d[Y] - c[Y] ) +
- c[X] * (double)( a[Y] - d[Y] ) +
- d[X] * (double)( c[Y] - a[Y] );
- /* code 2 handled here */
- if ( ((num < NEARZERO) && (num > -NEARZERO)) || (num == denom) ) code = 2;
- s = num / denom;
-
- num = -( a[X] * (double)( c[Y] - b[Y] ) +
- b[X] * (double)( a[Y] - c[Y] ) +
- c[X] * (double)( b[Y] - a[Y] ) );
+ num = a[X] * (double)( d[Y] - c[Y] ) +
+ c[X] * (double)( a[Y] - d[Y] ) +
+ d[X] * (double)( c[Y] - a[Y] );
+ /* code 2 handled here */
+ /*if ( ((num < NEARZERO) && (num > -NEARZERO)) || (num == denom) ) code = 2;*/
+ if ( (dAbs(num) < NEARZERO) || (dEqual(num,denom)) ) code = 2;
+
+ /* Debugging step 1
+ printf("\n!!! SegSeg step1: dAbs(num)=%f, dEqual(num,denom)=%d), code=%d\n",dAbs(num),dEqual(num,denom),code);
+ printf("\nNEARZERO=%f\n",NEARZERO);*/
- t = num / denom;
-
- if ( ((num < NEARZERO) && (num > -NEARZERO)) || (num == denom) ) code = 2;
-
- if ( (NEARZERO < s) && (s < 1.0) &&
- (NEARZERO < t) && (t < 1.0) )
- code = 1;
- else if ( (-NEARZERO > s) || (s > 1.0) ||
- (-NEARZERO > t) || (t > 1.0) )
- code = 0;
-
- p[X] = a[X] + s * ( b[X] - a[X] );
- p[Y] = a[Y] + s * ( b[Y] - a[Y] );
- }
+ s = num / denom;
+
+ num = -( a[X] * (double)( c[Y] - b[Y] ) +
+ b[X] * (double)( a[Y] - c[Y] ) +
+ c[X] * (double)( b[Y] - a[Y] ) );
+
+ t = num / denom;
+
+ /* if ( ((num < NEARZERO) && (num > -NEARZERO)) || (num == denom) ) code = 2;*/
+ if ( (dAbs(num) < NEARZERO) || (dEqual(num,denom)) ) code = 2;
+ /* Debugging step 2
+ printf("\n!!! SegSeg step2: dAbs(num)=%f, dEqual(num,denom)=%d), code=%d\n",dAbs(num),dEqual(num,denom),code);
+ printf("\nNEARZERO=%f\n",NEARZERO);*/
+
+ if ( (NEARZERO < s) && (s < 1.0) &&
+ (NEARZERO < t) && (t < 1.0) )
+ code = 1;
+ else if ( (-NEARZERO > s) || (s > 1.0) ||
+ (-NEARZERO > t) || (t > 1.0) )
+ code = 0;
+
+ p[X] = a[X] + s * ( b[X] - a[X] );
+ p[Y] = a[Y] + s * ( b[Y] - a[Y] );
+ }
+
+ /* Debugging step 3
+ printf("\n!!! SegSeg step3: final value of code=%d\n",code);*/
return code;
}
int Parallel( tPointd a, tPointd b, tPointd c, tPointd d, tPointd p )
{
/* Avoid to consider segments as parallel whenever two points are the same. */
- if ( (a[X]==b[X] && a[Y]==b[Y]) || (c[X]==d[X] && c[Y]==d[Y]) ) return 0;
+/*if ( (a[X]==b[X] && a[Y]==b[Y]) || (c[X]==d[X] && c[Y]==d[Y]) ) return 0;*/
+if ( (dEqual(a[X],b[X]) && dEqual(a[Y],b[Y])) || (dEqual(c[X],d[X]) && dEqual(c[Y],d[Y])) ) return 0;
+if ( Collinear( a, b, c)==0 ) return 0;
- if ( Collinear( a, b, c)==0 ) return 0;
-
if ( Between( a, b, c ) ) {
Assignpx( p, c );
return 3;
@@ -205,11 +243,11 @@
/* If ab not vertical, check betweenness on x; else on y. */
if ( a[X] != b[X] )
- return ((a[X] <= c[X]) && (c[X] <= b[X])) ||
- ((a[X] >= c[X]) && (c[X] >= b[X]));
+ return ((a[X] <= c[X]) && (c[X] <= b[X])) ||
+ ((a[X] >= c[X]) && (c[X] >= b[X]));
else
- return ((a[Y] <= c[Y]) && (c[Y] <= b[Y])) ||
- ((a[Y] >= c[Y]) && (c[Y] >= b[Y]));
+ return ((a[Y] <= c[Y]) && (c[Y] <= b[Y])) ||
+ ((a[Y] >= c[Y]) && (c[Y] >= b[Y]));
}
int Collinear( tPointd a, tPointd b, tPointd c )
More information about the adegenet-commits
mailing list