R Under development (unstable) (2024-07-05 r86875 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > require( zonohedra, quiet=TRUE ) Package: zonohedra. Author: Glenn Davis. Version: 0.3-0. Built: R 4.5.0; x86_64-w64-mingw32; 2024-07-07 02:06:25 UTC; windows Attaching package: 'zonohedra' The following object is masked from 'package:base': rank > > options( width=160 ) > > set.seed(0) > > > g.zonolist = list() > > g.zonolist[[1]] = zonohedron( classics.genlist[[3]] ) > names( g.zonolist )[1] = attr( g.zonolist[[1]], "fullname" ) > > g.zonolist[[2]] = zonohedron( classics.genlist[[6]] ) > names( g.zonolist )[2] = attr( g.zonolist[[2]], "fullname" ) > > g.zonolist[[3]] = zonohedron( classics.genlist[[13]] ) > names( g.zonolist )[3] = attr( g.zonolist[[3]], "fullname" ) > > #g.zonolist[[4]] = zonohedron( colorimetry.genlist[[1]] ) > #names( g.zonolist )[3] = attr( g.zonolist[[4]], "fullname" ) > > # for the 1nm data, inside.zonohedron() takes too long ! > g.zonolist[[4]] = zonohedron( colorimetry.genlist[[2]] ) > names( g.zonolist )[4] = attr( g.zonolist[[4]], "fullname" ) > > #g.zonolist[[4]] = spherize( g.zonolist[[1]] ) > #names( g.zonolist )[4] = paste( names[1], "after spherized", sep=' ', collapse='' ) > > > > #print( str(g.zonolist) ) > > > testRaytrace <- function( tol=5.e-13 ) + { + set.seed(0) + + for( k in 1:length(g.zonolist) ) # ) + { + cat( '\n' ) + cat( "--------- testRaytrace() ", k, ") ", names(g.zonolist)[k], '--------------\n' ) + + zono = g.zonolist[[k]] + + basepoints = 8 + rays = 10 + + whitept = 2 * zono$center + + for( i in 1:basepoints ) + { + s = i / (basepoints+1) + + base = s * whitept + + direction = matrix( rnorm(3*rays), nrow=rays ) + + df = raytrace( zono, base, direction, invert=TRUE ) #; print( str(df) ) + + point = df$point + + # check that distance of these boundary points is very small + distance = inside( zono, point )$distance #; print(distance) + + mask = abs(distance) <= tol + + ok = all( mask, na.rm=TRUE ) + if( ! ok ) + { + cat( "testRaytrace(). k=", k, " i=", i, " tol=", tol, "\n" ) + cat( " boundary failures=", sum( ! mask, na.rm=TRUE ), '\n' ) + return(FALSE) + } + + + # test inversion of these points on the boundary + pcube = df$pcube # invertboundary( zono, point, tol=tol )$pcube + if( is.null(pcube) ) + { + cat( "testRaytrace(). k=", k, " i=", i, " tol=", tol, "\n" ) + cat( " pcube is NULL, for boundary points.\n" ) + return(FALSE) + } + + delta = pcube %*% t(zono$matroid$matrix) - point #; print(delta) + ok = all( abs(delta) < tol ) #; print(ok) + if( ! ok ) + { + cat( "testRaytrace(). k=", k, " i=", i, "\n" ) + cat( " inversion failures for boundary points=", + sum( ! (abs(delta) < tol), na.rm=TRUE ), " max=", max(abs(delta)), '\n' ) + print( pcube ) + return(FALSE) + } + + + + # move a little from the boundary to the base and make sure the new points are all inside + test = 0.1 * matrix( base, nrow=nrow(point), ncol=3, byrow=TRUE ) + 0.9 * point + inside = inside( zono, test )$inside #; print(inside) + + ok = all( inside, na.rm=TRUE ) + if( ! ok ) + { + cat( "testRaytrace(). k=", k, " i=", i, "\n" ) + cat( " inside failures=", sum( ! inside, na.rm=TRUE ), '\n' ) + return(FALSE) + } + + # move a little past the boundary and make sure the new points are all outside + test = -0.1 * matrix( base, nrow=nrow(point), ncol=3, byrow=TRUE ) + 1.1 * point + outside = ! inside( zono, test )$inside #; print(outside) + + ok = all( outside, na.rm=TRUE ) + if( ! ok ) + { + cat( "testRaytrace(). k=", k, " i=", i, "\n" ) + cat( " outside failures=", sum( ! outside, na.rm=TRUE ), '\n' ) + return(FALSE) + } + } + } + + return(TRUE) + } > > > testSections <- function( genlimit=100, tol=5.e-9 ) + { + set.seed(0) + + + for( k in 1:length(g.zonolist) ) + { + cat( '\n\n' ) + cat( "--------- testSections() ", k, ") ", names(g.zonolist)[k], ' --------------\n', sep='' ) + flush.console() + + zono = g.zonolist[[k]] + + directions = 20 + + if( genlimit <= ncol( getmatrix(zono) ) ) directions = 1 + + sections = 5 + + direction = matrix( rnorm(3*directions), nrow=directions ) + + timeelapsed = 0 + timestart = zonohedra:::gettime() + + points = 0 + + for( i in 1:directions ) + { + normal = direction[i, ] + + # find max and min of the functional over the zonohedron + df = support( zono, rbind(normal,-normal) ) + valmax = df$value[1] + valmin = -df$value[2] + + # cat( "normal=", normal, " valmin=", valmin, "valmax=", valmax, '\n' ) + + # generate sections between min and max + s = (1:sections)/(sections+1) + + beta = (1-s)*valmin + s*valmax + + res = section( zono, normal, beta ) + + count = length(res) + + # count the number of points + for( j in seq_len(count) ) + points = points + nrow(res[[j]]$point) + + # if( genlimit <= ncol( getmatrix(zono) ) ) count=0 # too big for verification and inversion + + # do some verification and inversion + for( j in seq_len(count) ) + { + section = res[[j]]$point + + m = nrow(section) + + # cat( " direction=", i, " beta=", beta[j], " points=", m, '\n' ) + + # all points must be not NA + mask = is.finite( section ) + ok = all( mask ) + if( ! ok ) + { + cat( "testSections(). k=", k, " i=", i, " j=", j, " normal=", normal, "\n" ) + cat( " section finite failures=", sum( ! mask ), '\n' ) + return(FALSE) + } + + # check that distance of these boundary points is very small + distance = inside( zono, section )$distance #; print(distance) + + mask = abs(distance) <= tol * max( abs(valmax), abs(valmin) ) + + ok = all( mask ) + if( ! ok ) + { + cat( "testSections(). k=", k, " i=", i, " j=", j, " normal=", normal, " tol=", tol, "\n" ) + cat( " boundary failures=", sum( ! mask ), '\n' ) + cat( " distance=", distance, '\n' ) + return(FALSE) + } + + + # test inversion of these boundary points + df = invertboundary( zono, section, tol=tol ) + + if( is.null(df) ) + { + cat( "testSections(). k=", k, " i=", i, " tol=", tol, "\n" ) + cat( " pcube is NULL, for boundary points.\n" ) + return(FALSE) + } + + pcube = df$pcube + + if( any( is.na(pcube) ) ) + { + cat( "testSections(). k=", k, " i=", i, "\n" ) + cat( " normal =", normal, " beta =", beta[j], '\n' ) + mess = sprintf( " # of NAs in pcube is %d, out of %d values.\n", + sum(is.na(pcube)), length(pcube) ) + cat( mess ) + return(FALSE) + } + + if( df$transitions[1] %% 2 == 1 ) + { + cat( " transitions[1]=", df$transitions[1], '\n' ) + cat( " pcube=", pcube[1, ], '\n' ) + } + + delta = pcube %*% t(zono$matroid$matrix) - section #; print(delta) + ok = all( abs(delta) < tol ) #; print(ok) + if( ! ok ) + { + cat( "testSections(). k=", k, " i=", i, "\n" ) + cat( " inversion failures for boundary points=", + sum( ! (abs(delta) < tol), na.rm=TRUE ), " max=", max(abs(delta)), '\n' ) + # print( pcube ) + return(FALSE) + } + + } + } + + timeelapsed = zonohedra:::gettime() - timestart + + sections = directions * length(beta) + + cat( " sections=", sections, " time=", timeelapsed, " ", timeelapsed/sections, " per section", + " points=", points, " ", timeelapsed/points, " per point.\n", sep='' ) + + flush.console() + } + + return(TRUE) + } > > > > library( arrangements ) > library( zonohedra ) > > > ### test raytracing of points in non-trivial facets of a zonohedron > > # thematroid rank 3 matroid, does not have be simple > # subgrnd a subset of the ground set, must be contiguous (not verified) > # pairs the number of pairs to generate > # > # returns a data.frame with pairs rows, and these columns: > > buildmulti2trans <- function( thematroid, subgrnd, pairs ) + { + set.seed(0) + + ground = getground(thematroid) + + idxfromgnd = zonohedra:::idxfromgroundfun( ground ) + + idxpair = arrangements::combinations( idxfromgnd[subgrnd], 2, nsample=pairs ) #; print( idxpair ) + + dim(idxpair) = c(pairs,2) + + if( FALSE ) idxpair = rbind( idxpair, idxpair[ ,2:1,drop=F] ) + + n = length(ground) #; print( n ) + + idxgnd = ground[ idxpair ] + dim(idxgnd) = dim(idxpair) + + pcube = matrix( NA_real_, nrow=nrow(idxpair), ncol=n ) + colnames(pcube) = ground + + rank = rep( NA_integer_, nrow(idxpair) ) + + for( k in 1:nrow(idxpair) ) + { + alpha = runif(2) + + pcube[k, ] = zonohedra:::pcubefromdata( idxpair[k, ], alpha, n ) + + rank[k] = rank( thematroid, idxgnd[k,1] : idxgnd[k,2] ) + } + + out = data.frame( row.names=1:nrow(idxpair) ) + out$idxpair = idxpair + out$idxgnd = idxgnd + out$pcube = pcube + out$rank = rank + + return(out) + } > > > testRaytraceNT <- function( zono, pairs=200, tol1=1.e-4, tol2=1.e-2 ) + { + cat( '\n' ) + cat( "--------- testRaytraceNT() ", attr(zono,"fullname"), ' --------------\n' ) + flush.console() + + + ntcount = length(zono$zonogon) + if( ntcount == 0 ) return(FALSE) + + #matsimple = getsimplified(zono$matroid) + + thematroid = getmatroid(zono) + + df = buildmulti2trans( thematroid, thematroid$hyperplane[[1]], pairs ) + + matgen = getmatrix(thematroid) + + df$XYZ = df$pcube %*% t(matgen) + + # print( str(df) ) + + base = getcenter(zono) + direction = df$XYZ - matrix(base,nrow=nrow(df$XYZ),ncol=3,byrow=T) + + dftest = raytrace( zono, base, direction, invert=TRUE ) + + dftest$idxgnd = df$idxgnd + + dftest$rank = df$rank + + matgenorg = getmatrix(zono) + + XYZ = dftest$pcube %*% t(matgenorg) + + dftest$deltaXYZ = apply( abs(df$XYZ - XYZ), 1L, max ) # rowMeans( abs(df$XYZ - XYZ) ) + + dftest$deltacube = apply( abs(df$pcube - dftest$pcube), 1L, max ) # rowMeans( abs(df$pcube - dftest$pcube) ) + + # change deltacube entries with rank=1 to NA. These are degenerate. + mask = df$rank<2 + dftest$deltacube[ mask ] = NA_real_ + + #print(dftest) + + mask = dftest$transitions == 2 + ok = all( mask ) + if( ! ok ) + { + cat( "testraytraceNT(). FAIL. transitions==2 failures=", sum( ! mask, na.rm=TRUE ), + " max(transitions) =", max( dftest$transitions ), '\n' ) + return(FALSE) + } + + deltamax = max( dftest$deltaXYZ, na.rm=T ) + if( tol1 < deltamax ) + { + mess = sprintf( "testraytraceNT(). FAIL. max(deltaXYZ)=%g > %g.\n", deltamax, tol1 ) + cat( mess ) + return(FALSE) + } + + deltamax = max( dftest$deltacube, na.rm=T ) + if( tol2 < deltamax ) + { + kmax = which.max(dftest$deltacube) + mess = sprintf( "testraytraceNT(). FAIL. max(deltacube)=%g > %g. argmax=%d,%d\n", + deltamax, tol2, dftest$idxgnd[kmax,1], dftest$idxgnd[kmax,2] ) + cat( mess ) + return(FALSE) + } + + return( TRUE ) + } > > > if( ! testRaytrace() ) stop( "testRaytrace() failed !" ) --------- testRaytrace() 1 ) Bilinski dodecahedron -------------- --------- testRaytrace() 2 ) rhombic triacontahedron -------------- --------- testRaytrace() 3 ) truncated small rhombicosidodecahedron -------------- --------- testRaytrace() 4 ) xyz at 1nm step -------------- > > # for testing non-trivial facets, use the one facet in zonocol.1nm > if( ! testRaytraceNT( g.zonolist[[4]] ) ) stop( "testRaytraceNT() failed !" ) --------- testRaytraceNT() xyz at 1nm step -------------- > > if( ! testSections( 400 ) ) stop( "testSections() failed !" ) --------- testSections() 1) Bilinski dodecahedron -------------- sections=100 time=0.4591105 0.004591105 per section points=648 0.0007085039 per point. --------- testSections() 2) rhombic triacontahedron -------------- sections=100 time=0.4145967 0.004145967 per section points=1082 0.0003831762 per point. --------- testSections() 3) truncated small rhombicosidodecahedron -------------- sections=100 time=9.759187 0.09759187 per section points=1898 0.005141827 per point. --------- testSections() 4) xyz at 1nm step -------------- sections=5 time=21.85773 4.371546 per section points=3544 0.006167531 per point. > > cat( "\nPassed all zonohedron tests !\n", file=stderr() ) Passed all zonohedron tests ! > > proc.time() user system elapsed 33.76 12.10 45.92