Nothing
########### argument processing ##############
#
# A a non-empty numeric NxM matrix, or something that can be converted to be one
#
# Nmin the minimum allowed number of rows
#
# returns such a matrix, or NULL in case of error
prepareNxM <- function( A, M, Nmin=1 )
{
ok = is.numeric(A) && M*Nmin<=length(A) && (length(dim(A))<=2) # && (0<M)
ok = ok && ifelse( is.matrix(A), ncol(A)==M, ((length(A) %% M)==0) )
if( ! ok )
{
#print( "prepareNx3" )
#print( sys.frames() )
mess = substr( paste0(as.character(A),collapse=','), 1, 10 )
#arglist = list( ERROR, "A must be a non-empty numeric Nx3 matrix (with N>=%d). A='%s...'", mess )
#do.call( log.string, arglist, envir=parent.frame(n=3) )
#myfun = log.string
#environment(myfun) = parent.frame(3)
Aname = deparse(substitute(A))
# notice hack with 2L to make log.string() print name of parent function
#log.string( c(ERROR,2L), "Argument '%s' must be a non-empty numeric Nx%d matrix (with N>=%d). %s='%s...'",
# Aname, M, Nmin, Aname, mess )
log_level( ERROR, "Argument '%s' must be a non-empty numeric Nx%d matrix (with N>=%d). %s='%s...'",
Aname, M, Nmin, Aname, mess, .topcall=sys.call(-2L) )
return(NULL)
}
if( ! is.matrix(A) )
A = matrix( A, ncol=M, byrow=TRUE )
return( A )
}
# A an integer vector, representing a set of integers
# B an integer vector, representing a set of integers
#
# returns TRUE or FALSE, depending on whether A is a subset of B
#
# subset1() is fastest, because there is no error checking
subset1 <- function( A, B )
{
all( is.finite( match(A,B) ) )
}
subset2 <- function( A, B )
{
length( setdiff(A,B) ) == 0
}
normalizeColumns <- function( mat )
{
normvec = sqrt( .colSums( mat^2, nrow(mat), ncol(mat) ) )
return( t( t(mat) / normvec ) )
}
# A a numeric matrix
# MARGIN 1 (vectors are the rows) or 2 (vectors are the columns)
# for each vector, divide by L2 norm to get a unit vector
#
normalizeMatrix <- function( A, MARGIN )
{
ok = is.double(A) && is.matrix(A)
if( ! ok )
{
return(NULL)
}
MARGIN = as.integer(MARGIN)
ok = length(MARGIN)==1 && MARGIN %in% 1L:2L
if( ! ok )
{
return(NULL)
}
# make a deep (non-shallow) copy of A, because C_normalizeMatrix() modifies in-place
out = duplicate(A)
# change matrix out[] "in place"
ok = .Call( C_normalizeMatrix, out, MARGIN )
if( ! ok ) return(NULL)
return( out )
}
# ground positive integer vector in strictly increasing order. Not checked.
#
# returns lookup table from these integers to the order in the sequence
#
# out[idx] is equivalent to match( idx, ground ) but faster
idxfromgroundfun <- function( ground )
{
# make lookup table from ground to raw column index
numgen = length(ground)
out = integer( ground[numgen] )
out[ ground ] = 1:numgen
return( out )
}
shortlabel <- function( ivec )
{
if( length(ivec) <= 3 )
out = paste(ivec,collapse='+')
else
{
ran = range(ivec)
out = sprintf( "%g+...+%g", ran[1], ran[2] )
}
return(out)
}
# x a numeric matrix
# e0 used when nrow(x) >= 1
# e1 used when nrow(x) >= 2
# ground integer vector labeling the columns of x
#
# the columns of x are taken as generators of a matroid
#
# returns a matrix with no "loops" and "multiple groups"
simplify.matrix <- function( x, e0=0, e1=1.e-6, ground=NULL, ... )
{
ok = is.numeric(x) && is.matrix(x)
if( ! ok )
{
log_level( ERROR, "matrix x is invalid." )
return(NULL)
}
if( is.integer(x) ) storage.mode(x) = 'double'
elist = list(e0,e1)
ok = all( sapply(elist,is.numeric) ) && all( sapply(elist,length) == 1L )
if( ! ok )
{
log_level( ERROR, "One of e0,e1 is invalid." )
return(NULL)
}
if( is.null(ground) )
ground = 1L:ncol(x)
if( length(ground) != ncol(x) )
{
log_level( ERROR, "ground is invalid, because the length is incorrect." )
return(NULL)
}
if( ! all( 0 < diff(ground) ) )
{
log_level( ERROR, "ground is invalid, because it is not in strictly increasing order." )
return(NULL)
}
# look for small columns in x; these are the loops
loopmask = apply( x, MARGIN=2, function(vec) { max(abs(vec)) } ) <= e0
loopraw = which( loopmask )
colnames(x) = as.character( ground )
# convert from matrix column indexes to ground indexes
gnd.noloops = ground[ ! loopmask ]
# extract the matrix of nonloops
x.noloops = x[ , ! loopmask, drop=FALSE ]
if( ncol(x.noloops) == 0 ) return(x.noloops) # special case
xunit = normalizeColumns( x.noloops ) #; print( xunit )
grp = findColumnGroups( xunit, e1, oriented=FALSE ) #; print(grp)
if( all(grp==0) ) return(x.noloops) # special case
# nonloop is a list of vectors, and not a vector
nonloop = setlistfromvec( grp, gnd.noloops )
condata = condenseMatrix( x, ground, nonloop )
if( is.null(condata) ) return(NULL)
return( condata$matrix )
}
# A a numeric matrix
# ground integer vector labeling the columns of A
# conspec a list of integer vectors, all subsets of ground.
# the corresponding columns of A are co-directional, with tolerance e1
#
# for each vector in conspec take the dominant direction of the corresponding columns,
# and copy that direction to output matrix.
# returns a list with items:
# matrix with a column for each vector in conspec. This is the simplified matrix. It has length(conspec) columns.
# multiplesupp a data.frame with a row for each group of multiples in conspec, and these columns
# colidx index of corresponding column in output matrix - this is the index of the simplified generators
# cmax coordinate with largest absolute value, used to compute the following
# mixed logical, which is TRUE iff the group has "mixed directions"
# major the longer vector in the zonoseg spanned by the generators in the group, always non-zero
# minor the shorter vector in the zonoseg; this is non-zero iff the group has "mixed directions"
condenseMatrix <- function( A, ground, conspec )
{
ok = is.double(A) && is.matrix(A)
if( ! ok )
{
return(NULL)
}
ok = length(ground) == ncol(A)
if( ! ok )
{
return(NULL)
}
# make inverse lookup vector
colfromgnd = integer( max(ground) )
colfromgnd[ ground ] = 1L:ncol(A)
if( ! is.list(conspec) )
{
log_level( ERROR, "Argument conspec is invalid." )
return(NULL)
}
if( ! subset1(unlist(conspec),ground))
{
return(NULL)
}
lenvec = lengths( conspec ) #sapply( conspec, length )
conspec1 = conspec
# change all the non-singleton lists to any valid single value
# these columns will be overwritten later
multiplemask = 2 <= lenvec
conspec1[ multiplemask ] = ground[1]
out = list()
out$matrix = A[ , colfromgnd[ unlist(conspec1) ], drop=F ]
# replace all the multiple columns with dominant direction
colidx = which( multiplemask )
m = length(colidx)
cmax = integer(m)
major = matrix( 0, nrow=m, ncol=nrow(A) )
minor = matrix( 0, nrow=m, ncol=nrow(A) )
mixed = logical(m)
for( i in seq_len(m) )
{
k = colidx[i]
idx = colfromgnd[ conspec[[k]] ]
res = dominantDirection( A[ ,idx, drop=FALSE] )
cmax[i] = res$cmax
out$matrix[ ,k] = res$dominant # .rowSums( A[ ,idx], nrow(A), length(idx) )
major[i, ] = res$major
minor[i, ] = res$minor
mixed[i] = any( res$minor != 0 )
}
colnames(out$matrix) = sapply( conspec, shortlabel )
multiplesupp = data.frame( row.names=colidx )
multiplesupp$colidx = colidx
multiplesupp$cmax = cmax
multiplesupp$mixed = mixed
multiplesupp$major = major
multiplesupp$minor = minor
out$multiplesupp = multiplesupp
return( out )
}
emptymultiplesupp <-function( m )
{
out = data.frame( row.names=character(0) )
out$colidx = integer(0)
out$major = matrix( 0, nrow=0, ncol=m )
out$minor = matrix( 0, nrow=0, ncol=m )
out$mixed = integer(0)
return(out)
}
# A a numeric matrix with rank 1
#
# thus the columns of A are collinear and generate a zonoseg Z in space
# there may be vectors in a single direction (0 is an endpoint of Z),
# or both directions (0 is in the interior of Z)
# In the latter case, we say the generators are "mixed".
#
# returns a list with these items:
# cmax the 1-based index of the coordinate with maximum absolute value
# dominant the difference between the 2 endpoints of the zonoseg
# with the direction chosen so agree with the largest of sums in either direction
# NB: out$dominant = out$major - out$minor
# major the largest of the sums
# minor the smaller of the sums, most often it is 0, which means not "mixed"
dominantDirection <- function( A )
{
ok = is.double(A) && is.matrix(A)
if( ! ok )
{
return(NULL)
}
# find which row has the largest norm,
# since the rows are all collinear too, any norm will do so use L^inf
cmax = arrayInd( which.max(abs(A)), dim(A) )[1]
rowmax = A[ cmax, ]
sumneg = rowSums( A[ , rowmax < 0, drop=FALSE ] )
sumpos = rowSums( A[ , 0 < rowmax, drop=FALSE ] )
out = list()
out$cmax = cmax
delta = abs(sumpos[cmax]) - abs(sumneg[cmax])
if( delta == 0 )
# this means that sumpos and sumneg are exact opposites
# in this case, choose the one with positive value at imax
delta = sumpos[cmax] - sumneg[cmax]
if( 0 < delta )
{
out$major = sumpos ; out$minor = sumneg
}
else
{
out$major = sumneg ; out$minor = sumpos
}
out$dominant = out$major - out$minor
return(out)
}
# base a non-zero n-vector
# direction a mxn matrix, whose rows are taken as n-vectors, AND are all multiples of base
#
# returns an m-vector of values, all are +1, -1, or 0, depending on the multiple
#
# a halfspace test would work, but that would be a little slower
signvector <- function( base, direction )
{
if( length(base) != ncol(direction) )
{
log_level( FATAL, "Internal error. %d != %d.", length(base), ncol(direction) )
return(NULL)
}
# find which component of base has the largest norm,
# since the rows of direction are all collinear, any norm will do so use L^inf
imax = which.max( abs(base) )
out = sign( base[imax] ) * sign( direction[ , imax ] )
return(out)
}
# grp integer vector, as returned from grpDuplicated(). 0s are the singletons.
# ground integer vector of the point indexes
#
# returns a list of integer vectors, as indexed by ground[],
# each of which is a group of multiples with more than 1 point; compare with setlistfromvec().
grplistfromvec <- function( grp, ground=NULL )
{
if( is.null(ground) )
ground = 1L:length(grp)
else if( length(ground) != length(grp) )
{
log_level( ERROR, "length(ground)= %d is invalid!", length(ground) )
return(NULL)
}
n = max(grp)
out = vector( n, mode='list' )
for( i in seq_len(n) )
out[[i]] = ground[ which( grp == i ) ]
return(out)
}
# grp integer vector, as returned from grpDuplicated(). 0s are the singletons.
# ground integer vector of the point indexes, with the same length as grp
#
# returns a list of integer vectors, as indexed by ground[],
# This list of sets includes both singletons and multiples; compare with grplistfromvec().
# We always have length(out) <= length(grp)
setlistfromvec <- function( grp, ground=NULL )
{
if( is.null(ground) )
ground = 1L:length(grp)
else if( length(ground) != length(grp) )
{
log_level( ERROR, "length(ground)= %d is invalid!", length(ground) )
return(NULL)
}
out = as.list( ground )
keep = ! logical( length(out) )
n = max(grp)
for( i in seq_len(n) )
{
idx = which( grp == i )
out[[ idx[1] ]] = ground[ idx ]
# remove all except the first point in the group
keep[ idx[-1] ] = FALSE
}
out = out[ keep ]
return(out)
}
# xN one of:
# *) a list of integer vectors
# *) an integer vector
# *) NULL (ignored)
#
# returns the union of all inputs, in strictly ascending order
fastunion <- function( x1, x2=NULL, x3=NULL )
{
return( .Call( C_fastunion, x1, x2, x3 ) )
}
# hyper a list of integer vectors, each one of them defining a nontrivial hyperplane subset, and in increasing order. Not checked
# ground an integer vector - the union of all the sets in hyper - and in increasing order. Not checked
# returns a list of 2-point hyperplanes that, combined with hyper, form a 2-partition of the ground set
# if this is not possible, returns a character message with info on the problem
trivialhypers2 <- function( hyper, ground )
{
out = .Call( C_trivialhypers2, hyper, ground )
if( ! is.null(out$cmax) )
{
# ERROR. make a good error message
mess = "The hyperplanes do not satisfy the paving matroid properties for rank=3."
mess = c( mess, " the point pair %d,%d appears in %d hyperplanes." )
mess = paste0( mess, sep='\n' )
mess = sprintf( mess, out$pmax[1], out$pmax[2], out$cmax )
# mess = c( mess, " Try reducing argument e2." )
return( mess )
}
return( out )
}
matrix2list <- function( x, MARGIN )
{
return( .Call( C_matrix2list, x, as.integer(MARGIN) ) )
}
incidencedata <- function( hyper, ground=NULL )
{
if( is.null(ground) ) ground = fastunion(hyper)
return( .Call( C_incidencedata, hyper, ground ) )
}
incidencematrix <- function( hyper, ground, subset )
{
out = .Call( C_incidencematrix, hyper, ground, subset )
colnames(out) = as.character(subset)
return( out )
}
findRunsTRUE <- function( mask, periodic=FALSE )
{
# put sentinels on either end, to make things far simpler
dif = diff( c(FALSE,mask,FALSE) )
start = which( dif == 1 )
stop = which( dif == -1 )
if( length(start) != length(stop) )
{
log_level( FATAL, "Internal error. %d != %d", length(start), length(stop) )
return(NULL)
}
stop = stop - 1L
if( periodic && 2<=length(start) )
{
m = length(start)
if( start[1]==1 && stop[m]==length(mask) )
{
# merge first and last
start[1] = start[m]
start = start[ 1:(m-1) ]
stop = stop[ 1:(m-1) ]
}
}
return( cbind( start=start, stop=stop ) )
}
# A 2xn matrix, with n vectors in the columns generating a convex cone
# there are no 0s or multiples, so do not have to worry about generators in both directions
#
# if the cone is salient, returns the indexes of the 2 columns, in CCW order
# if the cone is not salient, returns integer(0)
cxconegenerators <- function( A )
{
n = ncol(A)
if( nrow(A)!=2 || n<=1 ) return(NULL)
theta = atan2( A[2, ], A[1, ] )
perm = order( theta )
theta_sorted = theta[perm]
gapvec = c( diff(theta_sorted), 2*pi - (theta_sorted[n] - theta_sorted[1]) )
kmax = which.max( gapvec )
if( gapvec[kmax] <= pi )
# the vectors generate the entire plane, NOT salient
return( integer(0) )
# there is a gap with angle > pi
if( kmax < n ) k2 = kmax+1
else k2 = 1
out = perm[ c(k2,kmax) ] # ;match( c(k2,kmax), perm )
return( out )
}
# u numeric N-vector
#
# returns Nx2 matrix with each row on the circle
tocircle <- function( u, tol=5.e-16 )
{
z = exp( u * 1i )
out = cbind( Re(z), Im(z) )
idx = which( abs(out) < tol )
if( 0 < length(idx) ) out[idx] = 0
return(out)
}
# u, v unit vectors of dimension n
#
# returns nxn rotation matrix that takes u to v, and fixes all vectors ortho to u and v
#
# Characteristic Classes, Milnor and Stasheff, p. 77
rotationshortest <- function( u, v )
{
n = length(u)
if( length(v) != n ) return(NULL)
uv = u + v
if( all(uv == 0) ) return(NULL)
out = diag(n) - (uv %o% uv)/(1 + sum(u*v)) + 2*(v %o% u)
return(out)
}
# dir non-zero 3-vector
#
# returns a 3x3 orthogonal matrix where the 3rd column is parallel to dir
# and the 1st and 2nd are orthogonal to dir
# Thus frame3x3 rotates dir to the z-axis
goodframe3x3 <- function( dir )
{
out = base::svd( dir, nu=3 )$u
# the 1st column is a multiple of dir, but it might be a negative multiple
if( sum( dir * out[ ,1] ) < 0 )
# reverse sign of 1st column
out[ ,1] = -out[ ,1]
# move 1st column to the 3rd
# and ensure that determinant is positive, so it's a rotation
if( 0 < det(out) )
perm = c(2,3,1) # even
else
perm = c(3,2,1) # odd
out = out[ , perm ]
return( out )
}
# point nx3 matrix of points on a great circle of S^2
# normal unit normal to the plane spanned by the great circle, this serves to orient the circle
#
# returns a permutation of 1:n that puts the points in counter-clockwise order,
# using the right-hand-rule with normal
orderoncircle <- function( point, normal )
{
# find a 3x3 rotation matrix that rotates normal to the north pole
pole = c(0,0,1)
if( all( abs(normal+pole) < 1.e-2 ) )
{
# tweak normal away from -pole,
# this will change the projected circle to an ellipse but the order is still the same
normal = c(1,12,-12)/17 # from a Pythagorean quadruple
}
Q = rotationshortest( normal, pole ) #; print( Q )
# compute xy, the z coordinate is very near 0 and ignored
xy = point %*% t( Q[1:2, ] ) #; print(xy) ; print( sqrt( rowSums(xy^2) ) )
theta = atan2( xy[ ,2], xy[ ,1] )
perm = order(theta)
return( perm )
}
# normal a non-zero 3-vector
#
# returns a 3x2 matrix where the columns complete normal to an orthogonal basis
# the 2 columns are unit vectors
# if normal is added as a 3rd column, the 3x3 matrix preserves orientation
# if normal is replaced by -normal, the columns of output matrix are swapped.
frame3x2fun <- function( normal, axischeck=FALSE )
{
ok = is.numeric(normal) && length(normal)==3 && 0<sum(abs(normal))
if( ! ok )
{
log_level( ERROR, "argument normal is invalid." )
return(NULL)
}
out = base::svd( normal, nu=3 )$u[ , 2:3 ] #; print( out )
test = crossproduct( out[ ,1], out[ ,2] )
if( sum(test*normal) < 0 )
{
# swap columns
out = out[ , 2L:1L ] #; cat( 'frame3x2fun(). columns swapped !\n' )
}
if( axischeck && isaxis( -out[ ,1] ) )
{
out[ ,1] = -out[ ,1]
out = out[ , 2L:1L ] # swap
}
else if( axischeck && isaxis( -out[ ,2] ) )
{
out[ ,2] = -out[ ,2]
out = out[ , 2L:1L ] # swap
}
if( FALSE )
{
# multiplication check
test = normal %*% out
if( 1.e-14 < max(abs(test)) )
{
log_level( ERROR, "frame3x2fun() failed orthogonal test = %g!", max(abs(test)) )
}
test = t(out) %*% out - diag(2)
if( 1.e-14 < max(abs(test)) )
{
log_level( ERROR, "frame3x2fun() failed product test = %g!", max(abs(test)) )
}
}
if( FALSE )
{
# orientation check
if( determinant( cbind(out,normal) )$sign <= 0 )
{
log_level( ERROR, "frame3x2fun() failed sign test. det = %g!", det( cbind(out,normal) ) )
}
}
return(out)
}
isaxis <- function( vec )
{
n = length(vec)
return( sum(vec==0)==n-1 && sum(vec==1)==1 )
}
# vec numeric vector
# skip integer coordinate to skip (ignore)
allequalskip <- function( vec, skip )
{
vec = vec[-skip]
return( all( vec==vec[1] ) )
}
gettime <- function()
{
return( microbenchmark::get_nanotime() * 1.e-9 )
}
createtimer <- function()
{
now = microbenchmark::get_nanotime() * 1.e-9
# using a list
return( list( created=now, now=now, elapsed=0, total=0 ) )
# using an array of doubles is actually slower !
#out = c(now,now,0,0)
#names(out) = c("created","now","elapsed","total")
#return( out )
}
updatetimer <- function( x, reset=TRUE )
{
now = microbenchmark::get_nanotime() * 1.e-9
out = x
out$elapsed = now - out$now
out$total = now - out$created
if( reset ) out$now = now
# using integers
#out[3L] = now - out[2L]
#out[4L] = now - out[1L]
#if( reset ) out[2L] = now
# using names
#out['elapsed'] = now - out['now']
#out['total'] = now - out['created']
#if( reset ) out['now'] = now
return(out)
}
# vertex 2n x m matrix with vertices of a convex polygon in the rows
#
# returns 4(n-1) x m matrix with tiling of the polygon into convex quadrangles
# If m=3, it is ready to pass to rgl::quad3d()
makequads <- function( vertex )
{
n = nrow(vertex) / 2L
ok = 2<=n && as.integer(n)==n
if( ! ok )
{
log_level( ERROR, "nrow(vertex)=%g is invalid.", nrow(vertex) )
return(NULL)
}
mat = matrix( c( 1:(n-1), 2:n, (2*n-1):(n+1), (2*n):(n+2) ), ncol=4 )
idx = as.integer( t(mat) ) #; print(idx)
out = vertex[ idx, ]
return(out)
}
# p a point in the n-cube, which we can think of as a transmittance spectrum
#
# returns the min number of transitions between 0 and 1 necessary to achieve such a p, including interpolation
# it always returns an even integer
#
# this is the fast C version
transitioncount <- function( p )
{
.Call( C_transitioncount, p )
}
# p a point in the n-cube, which we can think of as a transmittance spectrum
#
# returns the min number of transitions between 0 and 1 necessary to achieve such a p, including interpolation
# it always returns an even integer
#
# this is the slow R version
transitioncount_old <- function( p )
{
n = length(p)
if( n == 0 ) return(0L)
# find all runs of points in the interior of [0,1], in a periodic way
interior = 0<p & p<1
if( all(interior) )
# special case
return( as.integer( 2 * floor( (n+1)/2 ) ) )
mat = findRunsTRUE( interior, periodic=TRUE )
transitions = 0
if( 0 < nrow(mat) )
{
inext = c(2:n,1)
iprev = c(n,1:(n-1))
for( i in 1:nrow(mat) )
{
start = mat[i,1]
stop = mat[i,2]
m = stop - start + 1
if( m < 0 ) m = m + n
same = p[ iprev[start] ] == p[ inext[stop] ]
inc = ifelse( same, 2*floor( (m+1)/2 ), 2*floor( m/2 ) )
transitions = transitions + inc
}
}
# now remove all the interior coordinates
ppure = p[ ! interior ]
# add usual 0-1 transitions
transitions = transitions + sum( diff(ppure) != 0 )
# add wrap-around transition, if present
if( ppure[1] != ppure[ length(ppure) ] ) transitions = transitions + 1
return( as.integer(transitions) )
}
# this is only valid if 1 <= i < j <= n
PAIRINDEX <- function( i, j, n )
{
(i-1)*n - ((i)*(i+1))/2 + j
}
# subground increasing integer M-vector, and subvector of ground, NOT verified
# ground increasing integer N-vector, with M <= N
#
# each of these inputs has a matrix of all pairs, with standard ordering
# for subground there are M*(M-1)/2 pairs
# for ground there are N*(N-1)/2 pairs
#
# each pair for subground also appears in as a pair for ground
# The function returns an integer vector of length M*(M-1)/2
# giving the row in ground of each pair from subground
translateallpairs <- function( subground, ground )
{
# m = length(subground)
#subidxraw = .Call( C_allpairs, m )
#subidx = subground[ subidxraw ]
#dim(subidx) = dim(subidxraw)
subidx = allpairs( subground )
# subidx has ALL pairs from subground, in standard order
# each pair is also a pair from ground
# translate these to raw indexes in ground
idxfromground = idxfromgroundfun( ground )
idxraw = idxfromground[ subidx ]
dim(idxraw) = dim(subidx)
n = length(ground)
# idxraw has some pairs taken from 1:n
# find their position in the standard order of ALL pairs - N(N-1)/2 of them
out = .Call( C_pairindex, idxraw, n )
return( out )
}
# .vec1 and .vec2 non-zero vectors of the same dimension
#
angleBetween <- function( .vec1, .vec2, unitized=FALSE, eps=5.e-14 )
{
q = sum( .vec1*.vec2 )
if( ! unitized )
{
len1 = sqrt( sum(.vec1^2) )
len2 = sqrt( sum(.vec2^2) ) #; print( denom )
denom = len1 * len2
if( abs(denom) < eps ) return( NA_real_ )
q = q / denom #; print(q)
}
if( abs(q) < 0.99 )
{
# the usual case uses acos
out = acos(q)
}
else
{
# use asin instead
if( ! unitized )
{
.vec1 = .vec1 / len1
.vec2 = .vec2 / len2
}
if( q < 0 ) .vec2 = -.vec2
d = .vec1 - .vec2
d = sqrt( sum(d*d) )
out = 2 * asin( d/2 )
if( q < 0 ) out = pi - out
}
return(out)
}
# returns NULL unless ALL the values are valid integers
# if some are invalid, prints a warning message
intfromchar <- function( charvec )
{
if( is.null(charvec) || ! is.character(charvec) ) return(NULL)
bad = ! grepl( "[ ]*-?[0-9.]+[ ]*", charvec )
if( any(bad) )
{
log_level( WARN, "%d of %d values are invalid integers.", sum(bad), length(bad) )
return(NULL)
}
out = as.integer( charvec )
bad = (out != as.double(charvec))
bad[ is.na(bad) ] = TRUE
if( any(bad) )
{
log_level( WARN, "%d of %d values are invalid integers.", sum(bad), length(bad) )
return(NULL)
}
return( out )
}
################################# deadwood below ##########################
# A a numeric matrix
# ground integer vector labeling the columns of A
# conspec one of these:
# *) integer vector defining a subset of ground
# *) a list of integer vectors, all subsets of ground
# For the integer vector, just copy those columns to output.
# For a list, for each vector take the dominant direction of the corresponding columns,
# and copy that direction to output.
# returns a matrix with length(conspec) columns
condenseMatrix_old <- function( A, ground, conspec )
{
ok = is.double(A) && is.matrix(A)
if( ! ok )
{
return(NULL)
}
ok = length(ground) == ncol(A)
if( ! ok )
{
return(NULL)
}
# make inverse lookup vector
colfromgnd = integer( max(ground) )
colfromgnd[ ground ] = 1L:ncol(A)
if( is.integer(conspec) )
{
# this is the easy case
if( ! subset1(conspec,ground) )
{
return(NULL)
}
out = A[ , colfromgnd[conspec], drop=F ]
}
else if( is.list(conspec) )
{
if( ! subset1(unlist(conspec),ground))
{
return(NULL)
}
lenvec = lengths( conspec ) #sapply( conspec, length )
conspec1 = conspec
# change all the non-singleton lists to any valid single value
multiplemask = 2 <= lenvec
conspec1[ multiplemask ] = ground[1]
out = A[ , colfromgnd[ unlist(conspec1) ], drop=F ]
# replace all the multiple columns with sums
multiple = which( multiplemask )
for( k in multiple )
{
idx = colfromgnd[ conspec[[k]] ]
out[ ,k] = dominantDirection( A[ ,idx, drop=FALSE] )$dominant # .rowSums( A[ ,idx], nrow(A), length(idx) )
}
}
else
{
log_level( ERROR, "conspec is invalid." )
return(NULL)
}
colnames(out) = sapply( conspec, shortlabel )
return( out )
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.