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=144 ) > > set.seed(0) > > > g.zonolist = list() > > g.zonolist[[1]] = zonogon( matrix( 1:4, nrow=2, ncol=2 ) ) > names( g.zonolist )[1] = "parallelogram" > > g.zonolist[[2]] = symmetrize( g.zonolist[[1]] ) > names( g.zonolist )[2] = "symmparallelogram" > > g.zonolist[[3]] = polarzonogon( 17 ) > names( g.zonolist )[3] = "polar" > > g.zonolist[[4]] = polarzonogon( 17, 9 ) > names( g.zonolist )[4] = "halfpolar" > > mat = matrix( rnorm(2*50), nrow=2 ) > g.zonolist[[5]] = zonogon( mat ) > names( g.zonolist )[5] = "random" > > mat = rbind( rnorm(25), runif(25) ) > g.zonolist[[6]] = zonogon( mat ) > names( g.zonolist )[6] = "halfrandom" > > > #print( str(g.zonolist) ) > > > testInsideAndInversion <- function( tol=5.e-12 ) # before May 2022 it was smaller 1.e-14, but this failed for the "random" case + { + set.seed(0) + + for( k in 1:length(g.zonolist) ) + { + cat( '\n' ) + cat( "--------- testInsideAndInversion() ", k, ") ", names(g.zonolist)[k], '--------------\n' ) + + zono = g.zonolist[[k]] + + basepoints = 9 + rays = 10 + + whitept = 2 * zono$center + + for( i in 1:(basepoints-1) ) + { + s = i / (basepoints) + + base = s * whitept + + direction = matrix( rnorm(2*rays), nrow=rays ) + + df = raytrace( zono, base, direction ) #; print( str(df) ) + point = df$point + + #print( df ) + + # 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( "testInsideAndInversion(). k=", k, " i=", i, " tol=", tol, "\n" ) + cat( " raytrace boundary failures=", sum( ! mask, na.rm=TRUE ), " max(abs(distance))=", max(abs(distance)), '\n' ) + return(FALSE) + } + + # test inversion of these boundary points + #print( point ) + pcube = invert( zono, point, tol=tol )$pcube + if( is.null(pcube) ) + { + cat( "testInsideAndInversion(). k=", k, " i=", i, " tol=", tol, "\n" ) + cat( "pcube is NULL, for boundary points.\n" ) + return(FALSE) + } + + #print( pcube ) + + delta = pcube %*% t(zono$matroid$matrix) - point #; print(delta) + ok = all( abs(delta) < tol ) #; print(ok) + if( ! ok ) + { + cat( "testInsideAndInversion(). 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) + } + + + # on each segment from base to boundary, + # pick a point at random and make sure the new points are all inside + lambda = runif( nrow(point) ) + pmat = lambda * matrix( base, nrow=nrow(point), ncol=2, byrow=TRUE ) + (1-lambda) * point + inside = inside( zono, pmat )$inside #; print(inside) + + ok = all( inside, na.rm=TRUE ) + if( ! ok ) + { + cat( "testInsideAndInversion(). k=", k, " i=", i, "\n" ) + cat( " inside failures=", sum( ! inside, na.rm=TRUE ), '\n' ) + return(FALSE) + } + + # test inversion of these inside points + pcube = invert( zono, pmat, tol=tol )$pcube + if( is.null(pcube) ) + { + cat( "testInsideAndInversion(). k=", k, " i=", i, " tol=", tol, "\n" ) + cat( "pcube is NULL, for inside points.\n" ) + return(FALSE) + } + + delta = pcube %*% t(zono$matroid$matrix) - pmat #; print(delta) + ok = all( abs(delta) < tol ) #; print(ok) + if( ! ok ) + { + cat( "testInsideAndInversion(). k=", k, " i=", i, "\n" ) + cat( " inversion failures on inside =", sum( ! (abs(delta) < tol), na.rm=TRUE ), " max=", max(abs(delta)), '\n' ) + print( pcube ) + 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=2, byrow=TRUE ) + 1.1 * point + outside = ! inside( zono, test )$inside #; print(outside) + + ok = all( outside, na.rm=TRUE ) + if( ! ok ) + { + cat( "testInsideAndInversion(). k=", k, " i=", i, "\n" ) + cat( " outside failures=", sum( ! outside, na.rm=TRUE ), '\n' ) + return(FALSE) + } + } + } + + return(TRUE) + } > > > testSections <- function( tol = 1.e-14) + { + set.seed(0) + + for( k in 1:length(g.zonolist) ) + { + cat( '\n' ) + cat( "--------- testSections() ", k, ") ", names(g.zonolist)[k], '--------------\n' ) + + zono = g.zonolist[[k]] + + directions = 100 + + direction = matrix( rnorm(2*directions), nrow=directions ) + + for( i in 1:directions ) + { + normal = direction[i, ] + + # find max and min of the functional over the zonogon + df = support( zono, rbind(normal,-normal) ) + valmax = df$value[1] + valmin = -df$value[2] + + # cat( "normal=", normal, " valmin=", valmin, "valmax=", valmax, '\n' ) + + # generate 9 sections between min and max + s = (1:9)/10 + beta = (1-s)*valmin + s*valmax + df = section( zono, normal, beta ) + boundary1 = df$boundary1 + boundary2 = df$boundary2 + + # all points must be not NA + point = rbind(boundary1,boundary2) + mask = is.finite( point ) + ok = all( mask ) + if( ! ok ) + { + cat( "testSections(). k=", k, " i=", i, " normal=", normal, "\n" ) + cat( " section failures=", sum( ! mask ), '\n' ) + return(FALSE) + } + + # check that distance of these boundary points is very small + distance = inside( zono, point )$distance #; print(distance) + + mask = abs(distance) <= tol + + ok = all( mask ) + if( ! ok ) + { + cat( "testSections(). k=", k, " i=", i, " normal=", normal, " tol=", tol, "\n" ) + cat( " boundary failures=", sum( ! mask ), '\n' ) + return(FALSE) + } + + # make a few interpolation points, they must all be inside + for( s in c(0.01,0.5,0.99) ) + { + test = (1-s)*boundary1 + s*boundary2 + inside = inside( zono, test )$inside + + ok = all( inside ) + if( ! ok ) + { + cat( "testSections(). k=", k, " i=", i, " normal=", normal, "\n" ) + cat( " inside failures=", sum( ! inside ), '\n' ) + return(FALSE) + } + } + + # make a few extrapolation points, they must all be outside + for( s in c(-1,-0.01,1.01,2) ) + { + test = (1-s)*boundary1 + s*boundary2 + outside = ! inside( zono, test )$inside + + ok = all( outside ) + if( ! ok ) + { + cat( "testSections(). k=", k, " i=", i, " normal=", normal, "\n" ) + cat( " outside failures=", sum( ! outside ), '\n' ) + return(FALSE) + } + } + } + } + + return(TRUE) + } > > testSums <- function() + { + thesum = g.zonolist[[1]] %+% g.zonolist[[2]] %+% g.zonolist[[3]] %+% g.zonolist[[4]] %+% g.zonolist[[5]] %+% g.zonolist[[6]] + + if( is.null(thesum) ) return(FALSE) + + print(thesum) + + return(TRUE) + } > > if( ! testSums() ) stop( "testSums() failed !" ) zonogon: number of facets: 188 [ 94 antipodal facet-pairs] generators with mixed-directions: 2 { 1 2 } center: 2.353629 10.15435 facets that contain 0: ( 0 facets) { } pointed: FALSE salient: FALSE perimeter: 238.4721 area: 4449.789 matroid: ground set: 107 points {1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107} hyperplanes: 94 [1-point: 83] [2-point: 9] [3-point: 2] rank: 2 loops: 0 {} multiple groups: 11 [2-point: 9] [3-point: 2] uniform: FALSE paving: TRUE simple: FALSE contiguous: FALSE This matroid is constructed from a 2x107 real matrix. The summary of the simplified matroid is: ground set: 94 points {1 2 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107} Point 1 corresponds to the multiple group {1 3 5} in the original matroid. Point 2 corresponds to the multiple group {2 4 6} in the original matroid. Point 7 corresponds to the multiple group {7 24} in the original matroid. Point 8 corresponds to the multiple group {8 25} in the original matroid. Point 9 corresponds to the multiple group {9 26} in the original matroid. Point 10 corresponds to the multiple group {10 27} in the original matroid. Point 11 corresponds to the multiple group {11 28} in the original matroid. Point 12 corresponds to the multiple group {12 29} in the original matroid. Point 13 corresponds to the multiple group {13 30} in the original matroid. Point 14 corresponds to the multiple group {14 31} in the original matroid. Point 15 corresponds to the multiple group {15 32} in the original matroid. hyperplanes: 94 [1-point: 94] rank: 2 loops: 0 {} multiple groups: 0 {} uniform: TRUE paving: TRUE simple: TRUE contiguous: TRUE This matroid is constructed from a 2x94 real matrix. > > if( ! testInsideAndInversion() ) stop( "testInsideAndInversion() failed !" ) --------- testInsideAndInversion() 1 ) parallelogram -------------- --------- testInsideAndInversion() 2 ) symmparallelogram -------------- --------- testInsideAndInversion() 3 ) polar -------------- --------- testInsideAndInversion() 4 ) halfpolar -------------- --------- testInsideAndInversion() 5 ) random -------------- --------- testInsideAndInversion() 6 ) halfrandom -------------- > > if( ! testSections() ) stop( "testSections() failed !" ) --------- testSections() 1 ) parallelogram -------------- --------- testSections() 2 ) symmparallelogram -------------- --------- testSections() 3 ) polar -------------- --------- testSections() 4 ) halfpolar -------------- --------- testSections() 5 ) random -------------- --------- testSections() 6 ) halfrandom -------------- > > cat( "\nPassed all zonogon tests !\n", file=stderr() ) Passed all zonogon tests ! > > warnings() > > proc.time() user system elapsed 7.31 0.20 7.51