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 )
> 
> 
> 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" )
> 
> #g.zonolist[[1]] = zonohedron( colorimetry.genlist[[2]] )
> #names( g.zonolist )[1]  =  attr( g.zonolist[[1]], "fullname" )
> 
> #g.zonolist[[4]] = spherize( g.zonolist[[1]] )
> #names( g.zonolist )[4]  = paste( names[1], "after spherized", sep=' ', collapse='' )
> 
> 
> 
> #print( str(g.zonolist) )
> 
> 
> 
> 
> 
> 
> 
> #   returns a data frame with these columns
> #       idxpair
> #       alpha
> #       direction
> 
> randomDirs <- function( matgen, base, level=0 )
+     {
+     out = randomAtLevel( matgen, level )
+ 
+     out$direction   = out$XYZ  - matrix( base, nrow=nrow(out$XYZ), ncol=ncol(out$XYZ), byrow=TRUE )
+     rownames( out$direction )   = 1:nrow(out$direction)
+ 
+     out$XYZ     = NULL
+ 
+     #rownames(out)   = 1:nrow(out)
+ 
+     return( out )
+     }
> 
> 
> 
> randomAtLevel <- function( matgen, level=0 )
+     {
+     set.seed(0)
+ 
+     n   = ncol(matgen)
+ 
+     level   = as.integer(level)
+ 
+     ok  = (0 <= level) && (level <= n-2)
+     if( ! ok )  return(NULL)
+ 
+     #white   = rowSums( matgen )
+ 
+     #   out = matrix( 0, 2*n, 3 )
+ 
+     idxpair = matrix( NA_integer_, n, 2 )
+     alpha   = matrix( NA_real_, n, 2 )
+     XYZ     = matrix( NA_real_, n, 3 )
+ 
+     for( k in 1:n )
+         {
+         shift   = k-1L
+ 
+         #if( shift == 0 )
+         #    perm    = 1:n
+         #else
+         #    perm    = c( (shift+1):n,1:shift )   # c( (n-shift+1):n, 1:(n-shift) )
+ 
+         #print( perm )
+ 
+         j1 = 1L
+         j2 = j1 + level + 1L
+ 
+         idxpair[k, ]    = ( ( c(j1,j2) + shift-1L ) %% n ) + 1L     #     perm[ c(j1,j2) ]
+         alpha[k, ]      = runif(2)
+ 
+         #pcube   = numeric(n)
+         #pcube[ c(j1,j2) ]   = alpha[k, ]
+         #if( 0 < level ) pcube[ (j1+1):(j2-1) ]  = 1
+         #pcube[perm]   = pcube   # rotate cyclic
+ 
+         pcube   = zonohedra:::pcubefromdata( idxpair[k, ], alpha[k, ], n )
+         #print( pcube )
+ 
+         # and multiply
+         XYZ[k, ] = as.double( matgen %*% pcube )
+ 
+         #   now add the antipodal data
+         #idxpair[k+n, ]  = perm[ c(j2,j1) ]
+         #alpha[k+n, ]    = 1 - alpha[k, ]
+         #XYZ[k+n, ]      = white - XYZ[k, ]
+         }
+ 
+     out = data.frame( row.names=1:(n) )
+     out$idxpair = idxpair
+     out$alpha   = alpha
+     out$XYZ     = XYZ
+ 
+     return( out )
+     }
> 
> 
> 
> randomCubeInterior <- function( n, reset=TRUE )
+     {
+     if( reset ) set.seed(0)
+ 
+     eps = 5.e-4
+ 
+     out = matrix( runif(n*3,min=eps,max=1-eps), n, 3 )
+ 
+     return( out )
+     }
> 
> 
> randomCubeBoundary <- function( n )
+     {
+     set.seed(0)
+ 
+     coord   = as.integer( runif(n,min=1,max=4) )
+ 
+     value   = ( sign( runif(n,min=-1,max=1) ) + 1 ) / 2
+ 
+     out = randomCubeInterior( n, FALSE )
+ 
+     out[ cbind(1:n,coord) ]    = value
+ 
+     return( out )
+     }
> 
> 
> 
> 
> 
> 
> testRaytrace <- function( tol=5.e-6 )
+     {
+     set.seed(0)
+ 
+     g.zonolist[[4]] = zonohedron( colorimetry.genlist[['xyz1931.1nm']] )
+     names( g.zonolist )[4]  =  attr( g.zonolist[[4]], "fullname" )
+ 
+ 
+     for( k in 1:length(g.zonolist) )
+         {
+         zono    = g.zonolist[[k]]
+ 
+         if( is.null(zono) ) next
+ 
+         cat( "\n---------  testRaytrace() ", k, ")  ", names(g.zonolist)[k], '  --------------\n', sep='' )
+         flush.console()
+ 
+         basepoints  = 3
+         #rays        = 10
+ 
+         whitept = 2 * getcenter( zono )
+ 
+         matsimp = getsimplified( getmatroid(zono) )
+ 
+         matgen  = getmatrix( matsimp )
+ 
+         m   = ncol(matgen)
+ 
+         idxfromgnd  = zonohedra:::idxfromgroundfun( getground( matsimp ) )
+ 
+         levelvec = c( 0, 1, as.integer(m/2), m-4 )
+ 
+         levelvec = c( 1, as.integer(m/4), as.integer(m/2), as.integer(2*m/3), m-4 )
+         #   levelvec = c( 1 )
+ 
+         timetrace   = numeric(0)
+ 
+         for( i in 1:basepoints )
+             {
+             s   = i / (basepoints+1)
+ 
+             base = s * whitept
+ 
+             for( level in levelvec )
+                 {
+                 dftest  = randomDirs( matgen, base, level )
+ 
+                 df  = raytrace2trans( zono, base, dftest$direction )     #; print( str(df) )
+ 
+                 #   check for failures
+                 if( any( is.na(df$tmax) ) )
+                     {
+                     cat( "testRaytrace(). base=", base, "  level=", level, '\n' )
+                     cat( "   tmax failures=", sum( is.na(df$tmax) ), '\n' )
+                     return(FALSE)
+                     }
+ 
+                 #   compare idxpair values
+                 idxpair     = idxfromgnd[ df$gndpair ]
+ 
+                 if( any( dftest$idxpair != idxpair ) )
+                     {
+                     cat( "testRaytrace(). base=", base, "  level=", level, '\n' )
+                     cat( "   idxpair failures=", sum( dftest$idxpair != idxpair ), '\n' )
+                     return(FALSE)
+                     }
+ 
+                 #   compare alpha values
+                 delta   = dftest$alpha - df$alpha
+ 
+                 good    = abs( delta ) <= tol
+                 if( any( ! good ) )
+                     {
+                     cat( "testRaytrace(). base=", base, "  level=", level, '\n' )
+                     cat( "    alpha failures=", sum( ! good ), "  tol=", tol, '\n' )
+                     cat( "    range(delta) = ", range(delta), '\n' )
+                     return(FALSE)
+                     }
+ 
+                 delta   = df$tmax - 1
+                 if( any( tol < abs(delta) ) )
+                     {
+                     count   = sum( tol < abs(delta) )
+                     mess    = sprintf( "%d of the %d tmax values are outside tol=%g.\n",
+                                             count, length(df$tmax), tol )
+                     cat( mess )
+                     return(FALSE)
+                     }
+ 
+                 timetrace   = c( timetrace, df$timetrace )
+ 
+                 if( FALSE )
+                     {
+                     cat( '\n' )
+                     mess = sprintf( "    fivenum summary of iters  --  total rays=%d  --\n", length(df$iters) )
+                     cat( mess )
+                     print( fivenum(df$iters) )
+                     cat( "    mean(iters) =", mean(df$iters), '\n' )
+                     flush.console()
+                     }
+                 }
+             }
+ 
+         cat( '\n' )
+         mess = sprintf( "fivenum summary of timetrace  --  total rays=%d  --\n", length(timetrace) )
+         cat( mess )
+         print( fivenum(timetrace) )
+ 
+         cat( "    mean(timetrace) =", mean(timetrace), '\n' )
+         }
+ 
+     return( TRUE )
+     }
> 
> 
> testCubeRays <- function( basepoints=50, rays=200,  tol=5.e-13 )
+     {
+     cat( "\n----   testCubeRays()  -------\n" ) ; flush.console()
+ 
+     cube    = zonohedron( classics.genlist[["C"]] )
+ 
+     #   make basepoints
+     basemat = randomCubeInterior( basepoints )
+ 
+     #   make boundary points
+     boundarymat     = randomCubeBoundary( rays )
+ 
+     df  = NULL
+     for( k in 1:basepoints )
+         {
+         base    = basemat[k, ]
+ 
+         dirmat  = boundarymat - matrix( base, nrow(boundarymat), length(base), byrow=TRUE )
+ 
+         dfsub   = raytrace2trans( cube, base, dirmat )
+ 
+         #   check for any errors
+         if( any( is.na(dfsub$tmax) ) )
+             {
+             count   = sum( is.na(dfsub$tmax) )
+             mess    = sprintf( "For basepoint %g,%g,%g,  there were %d failures, out of %d.\n",
+                                     base[1], base[2], base[3], count, length(dfsub$tmax) )
+             cat( mess )
+             return(FALSE)
+             }
+ 
+ 
+         if( tol < max( abs( boundarymat - dfsub$point ) ) )
+             {
+             count   = sum( tol < abs( boundarymat - dfsub$point ) )
+             mess    = sprintf( "For basepoint %g,%g,%g,  %d of the %d boundary points are outside tol=%g.\n",
+                                     base[1], base[2], base[3], count, nrow(boundarymat), tol )
+             cat( mess )
+             return(FALSE)
+             }
+ 
+         df  = rbind( df, dfsub )
+         }
+ 
+     delta   = df$tmax - 1
+     if( any( tol < abs(delta) ) )
+         {
+         count   = sum( tol < abs(delta) )
+         mess    = sprintf( "%d of the %d tmax values are outside tol=%g.\n",
+                                 count, length(df$tmax), tol )
+         cat( mess )
+         return(FALSE)
+         }
+ 
+     #print( df )
+ 
+ 
+     mess = sprintf( "    fivenum summary of timetrace  --  total rays=%d  --\n", length(df$timetrace) )
+     cat( mess )
+     print( fivenum(df$timetrace) )
+     cat( "    mean(timetrace) =", mean(df$timetrace), '\n' )
+ 
+     cat( '\n' )
+     mess = sprintf( "    fivenum summary of iters  --  total rays=%d  --\n", length(df$iters) )
+     cat( mess )
+     print( fivenum(df$iters) )
+     cat( "    mean(iters) =", mean(df$iters), '\n' )
+ 
+     return(TRUE)
+     }
> 
> 
> testPoles <- function( basepoints=50, tol=5.e-9 )
+     {
+     cat( "\n----   testPoles()  -------\n" ) ; flush.console()
+ 
+     set.seed(0)
+ 
+     zono    = zonohedron( classics.genlist[["BD"]] )
+ 
+     matgen  = getmatrix(zono)
+ 
+     n   = ncol( matgen )
+ 
+     #   make random points in n-cube
+     alpha   = matrix(  runif( n*basepoints,min=0.001, max=0.999 ), nrow=n, ncol=basepoints )
+ 
+     #   make basepoints
+     basemat = t( matgen %*% alpha )     #; print( basemat )
+ 
+     #   make rays through "white"
+     white   = 2 * getcenter(zono)
+     #direction   = matrix( white, nrow=nrow(basemat), ncol=3, byrow=TRUE ) - basemat
+ 
+     #   make rays through 0
+     #direction   = rbind( direction, -basemat )
+ 
+     for( k in 1:nrow(basemat) )
+         {
+         base    = basemat[k, ]
+ 
+         direction   = rbind( white - base, -base )
+ 
+         df  = raytrace2trans( zono, base, direction, tol=tol )  #; print( df )
+ 
+         #   gndpair must be all NA
+         if( ! all( is.na(df$gndpair) ) )
+             {
+             cat( "testPoles().  gndpair is invalid.\n" )
+             print( df$gndpair )
+             return(FALSE)
+             }
+ 
+         delta   = df$tmax - 1   #; print( range(delta) )
+ 
+         if( any( tol < abs(delta) ) )
+             {
+             count   = sum( tol < abs(delta) )
+             mess    = sprintf( "testPoles().  %d of the %d tmax values are outside tol=%g.\n",
+                                     count, length(df$tmax), tol )
+             cat( mess )
+             return(FALSE)
+             }
+         }
+ 
+     return( TRUE )
+     }
> 
> #   1600 munsell chips
> 
> testMunsell <- function( path=NULL, tol=5.e-6  )
+     {
+     cat( "\n----   testMunsell()  -------\n" ) ; flush.console()
+ 
+     if( is.null(path) ) path = system.file( "extdata/Munsell1600.txt", package='zonohedra' )
+ 
+     timestart   = zonohedra:::gettime()
+ 
+     zono = zonohedron( colorimetry.genlist[['xyz1931.1nm']] )
+ 
+     #   calibrate
+     sums    = rowSums( getmatrix(zono) )  #; print(sums)
+ 
+     #   use a diagonal scaling to match Illuminant E
+     zono = lintransform( zono, diag(100/sums) )
+ 
+     XYZ = as.matrix( read.table( path ) ) # ; print( str(XYZ) )
+ 
+     #XYZ = XYZ[ 1492:1541, ]
+     #print( XYZ / XYZ[ ,2] )
+ 
+     base    = getcenter(zono) #; print(2*base)
+ 
+     direction   = XYZ - matrix(base,nrow=nrow(XYZ),ncol=ncol(XYZ),byrow=TRUE)
+ 
+     df  = raytrace2trans( zono, base, direction, tol=tol )  #; print( df )
+ 
+     #   check for any errors
+     if( any( is.na(df$tmax) ) )
+         {
+         count   = sum( is.na(df$tmax) )
+         mess    = sprintf( "For basepoint %g,%g,%g,  there were %d failures, out of %d.\n",
+                                 base[1], base[2], base[3], count, length(df$tmax) )
+         cat( mess )
+         return(FALSE)
+         }
+ 
+     cat( "time elapsed = ", zonohedra:::gettime() - timestart, "sec  \n" )
+     cat( '\n' )
+ 
+     mess = sprintf( "    fivenum summary of timetrace  --  total rays=%d  --\n", length(df$timetrace) )
+     cat( mess )
+     print( fivenum(df$timetrace) )
+     cat( "    mean(timetrace) =", mean(df$timetrace), "   sum(timetrace) =", sum(df$timetrace),'\n' )
+ 
+     cat( '\n' )
+     mess = sprintf( "    fivenum summary of iters  --  total rays=%d  --\n", length(df$iters) )
+     cat( mess )
+     print( fivenum(df$iters) )
+     cat( "    mean(iters) =", mean(df$iters), '\n' )
+ 
+ 
+     return(TRUE)
+     }
> 
> 
> testSections <- function( genlimit=100, tol=5.e-9 )
+     {
+     set.seed(0)
+ 
+     g.zonolist[[1]] = zonohedron( classics.genlist[[3]] )
+ 
+     g.zonolist[[2]] = zonohedron( classics.genlist[[6]] )
+ 
+     g.zonolist[[3]] = zonohedron( classics.genlist[[8]] )
+ 
+     g.zonolist[[4]] = zonohedron( colorimetry.genlist[['xyz1931.1nm']] )
+ 
+     for( k in 1:length(g.zonolist) )
+         {
+         zono    = g.zonolist[[k]]
+ 
+         if( is.null(zono) ) next
+ 
+         cat( '\n\n' )
+         cat( "---------  testSections() ", k, ")   ", attr(zono,"fullname"), ' --------------\n', sep='' )
+         flush.console()
+ 
+         directions  = 20
+ 
+         if( genlimit <= ncol( getmatrix(zono) ) )   directions = 2
+ 
+         sections    = 5
+ 
+         #   make uniformly random directions
+         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 = section2trans( zono, normal, beta, invert=TRUE )
+ 
+             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) )
+                 {
+                 df      = res[[j]]
+ 
+                 point   = df$point
+ 
+                 m   = nrow(point)
+ 
+                 #   cat( "    direction=", i, "  beta=", beta[j], "  points=", m, '\n' )
+ 
+                 #   all points must be not NA
+                 mask        = is.finite( point )
+                 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)
+                     }
+ 
+ 
+                 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)
+                     }
+ 
+                 #   check inversion accuracy
+                 delta   =  pcube  %*%  t(zono$matroid$matrix) -  point   #; 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)
+     }
> 
> 
> if( ! testPoles() )     stop( "testPoles() failed !" )

----   testPoles()  -------
> 
> if( ! testCubeRays() )  stop( "testCubeRays() failed !" )

----   testCubeRays()  -------
    fivenum summary of timetrace  --  total rays=10000  --
[1] 9.380002e-05 1.542500e-04 2.060001e-04 2.361001e-04 3.640130e-02
    mean(timetrace) = 0.0002115886 

    fivenum summary of iters  --  total rays=10000  --
[1] 1 2 3 5 6
    mean(iters) = 3.42 
> 
> if( ! testRaytrace() )  stop( "testRaytrace() failed !" )

---------  testRaytrace() 4)  xyz at 1nm step  --------------

fivenum summary of timetrace  --  total rays=5100  --
[1] 0.0008896 0.0027791 0.0035153 0.0042231 0.0839710
    mean(timetrace) = 0.004178462 
> 
> if( ! testMunsell() )  stop( "testMunsell() failed !" )

----   testMunsell()  -------
time elapsed =  8.153702 sec  

    fivenum summary of timetrace  --  total rays=1600  --
[1] 0.0001637999 0.0029206000 0.0037648501 0.0044273501 0.0532092000
    mean(timetrace) = 0.004443158    sum(timetrace) = 7.109054 

    fivenum summary of iters  --  total rays=1600  --
[1]      0.0  34549.5  48606.5  73629.0 111353.0
    mean(iters) = 52043.72 
> 
> if( ! testSections() )  stop( "testSections() failed !" )


---------  testSections() 1)   Bilinski dodecahedron --------------
    sections=100  time=0.4208038  0.004208038 per section  points=648  0.0006493886 per point.


---------  testSections() 2)   rhombic triacontahedron --------------
    sections=100  time=0.0931396  0.000931396 per section  points=1092  8.529267e-05 per point.


---------  testSections() 3)   truncated rhombic dodecahedron --------------
    sections=100  time=0.1080512  0.001080512 per section  points=1374  7.863988e-05 per point.


---------  testSections() 4)   xyz at 1nm step --------------
    sections=10  time=5.320953  0.5320953 per section  points=7420  0.0007171096 per point.
> 
> 
> cat( "\nPassed all 2-transition tests !\n", file=stderr() )

Passed all 2-transition tests !
> 
> proc.time()
   user  system elapsed 
  26.62   15.04   41.76