R Under development (unstable) (2023-11-26 r85638 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2023 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.

> library(diptest)
> 
> stopifnot(dip(c(1,1,2,2)) == 1/4)# the maximal value possible: two point dist
> 
> ## very first small "unimodal" example --- the  1/(2*n) result:
> n <- length(u <- cumsum(0:3))
> d <- dip(u, debug=TRUE)# shows the final if() {added by MM} is really needed
dip() in C: n = 4; starting with  2N*dip = 1.
'dip': LOOP-BEGIN: 2n*D= 1         [low,high] = [  1,  4]; l_lcm/gcm = ( 4, 2)
  while(gcm[ix] != lcm[iv]) : ..
  calculating dip .. (dip_l, dip_u) = (0, 1)
'dip': LOOP-BEGIN: 2n*D= 1         [low,high] = [  1,  3]; l_lcm/gcm = ( 3, 2)
  while(gcm[ix] != lcm[iv]) : .
  calculating dip .. (dip_l, dip_u) = (0, 1)
'dip': LOOP-BEGIN: 2n*D= 1         [low,high] = [  1,  2]; l_lcm/gcm = ( 2, 2)
  ** (l_lcm,l_gcm) = (2,2) ==> d := 1
  calculating dip .. (dip_l, dip_u) = (0, 0)
No improvement in  low = 1  nor  high = 2 --> END
> stopifnot(d == dip(-u), d == 1/(2*n))# exact "=" for n = 4 !
> ## Note that I believe this should *not* give 0 (as fmechler@.. did),
> ## but rather 1/(2n) because that's  (1/n) / 2  and
> ## (1/n) is the correct distance between LCM and GCM
> 
> ## Small example -- but MM sees difference (32-bit / 64-bit):
> x <- c(0,2:3,5:6)
> d1 <- dip(x,   full=TRUE, debug=2)
dip() in C: n = 5; starting with  2N*dip = 1.
'dip': LOOP-BEGIN: 2n*D= 1         [low,high] = [  1,  5]; l_lcm/gcm = ( 2, 4)
  while(gcm[ix] != lcm[iv]) :
 G(3,2) --> ix = 2, iv = 2
 --> ix = 1, iv = 2
  calculating dip .. (dip_l, dip_u) = (1, 0)
'dip': LOOP-BEGIN: 2n*D= 1         [low,high] = [  2,  5]; l_lcm/gcm = ( 3, 3)
  while(gcm[ix] != lcm[iv]) :
 L(3,2) --> ix = 2, iv = 3
 G(2,3) --> ix = 1, iv = 3
  calculating dip .. (dip_l, dip_u) = (1.33333, 0) -> new larger dip 1.33333 (j_best = 3)
'dip': LOOP-BEGIN: 2n*D= 1.3333    [low,high] = [  4,  5]; l_lcm/gcm = ( 2, 2)
  ** (l_lcm,l_gcm) = (2,2) ==> d := 1
> d2 <- dip(6-x, full=TRUE, debug=2)
dip() in C: n = 5; starting with  2N*dip = 1.
'dip': LOOP-BEGIN: 2n*D= 1         [low,high] = [  1,  5]; l_lcm/gcm = ( 4, 2)
  while(gcm[ix] != lcm[iv]) :
 L(2,2) --> ix = 1, iv = 3
 L(2,3) --> ix = 1, iv = 4
  calculating dip .. (dip_l, dip_u) = (0, 1)
'dip': LOOP-BEGIN: 2n*D= 1         [low,high] = [  1,  4]; l_lcm/gcm = ( 3, 3)
  while(gcm[ix] != lcm[iv]) :
 L(3,2) --> ix = 2, iv = 3
 G(2,3) --> ix = 1, iv = 3
  calculating dip .. (dip_l, dip_u) = (1.33333, 0) -> new larger dip 1.33333 (j_best = 2)
'dip': LOOP-BEGIN: 2n*D= 1.3333    [low,high] = [  3,  4]; l_lcm/gcm = ( 2, 2)
  ** (l_lcm,l_gcm) = (2,2) ==> d := 1
> str(d1)
List of 15
 $ call       : language dip(x = x, full.result = TRUE, debug = 2)
 $ x          : num [1:5] 0 2 3 5 6
 $ n          : int 5
 $ dip        : num 0.133
 $ lo.hi      : int [1:2] 4 5
 $ ifault     : int 0
 $ gcm        : int [1:2] 5 4
 $ lcm        : int [1:2] 4 5
 $ mn         : int [1:5] 1 1 2 2 4
 $ mj         : int [1:5] 5 3 5 5 5
 $ min.is.0   : logi FALSE
 $ debug      : int 2
 $ xl         : num 5
 $ xu         : num 6
 $ full.result: logi TRUE
 - attr(*, "class")= chr "dip"
> str(d2)
List of 15
 $ call       : language dip(x = 6 - x, full.result = TRUE, debug = 2)
 $ x          : num [1:5] 0 1 3 4 6
 $ n          : int 5
 $ dip        : num 0.133
 $ lo.hi      : int [1:2] 3 4
 $ ifault     : int 0
 $ gcm        : int [1:2] 4 3
 $ lcm        : int [1:2] 3 4
 $ mn         : int [1:5] 1 1 1 3 1
 $ mj         : int [1:5] 2 4 4 5 5
 $ min.is.0   : logi FALSE
 $ debug      : int 2
 $ xl         : num 3
 $ xu         : num 4
 $ full.result: logi TRUE
 - attr(*, "class")= chr "dip"
> 
> if(!dev.interactive(orNone=TRUE)) pdf("ex1.pdf")
> par(mfrow = 2:1, mar = .1+c(3,4,2,1), mgp=c(1.5,.6,0), oma = c(0,0,2.1,0))
> #
> plot(d1)
> abline(v=-1:7, h = seq(0,1,by=0.2), lty="83", col = "gray")
> #
> plot(d2)
> abline(v=-1:7, h = seq(0,1,by=0.2), lty="83", col = "gray")
> #
> ## "title" only now
> mtext("dip() problem with 'mirror x'", side=3, line = 0.8,
+       outer=TRUE, cex = 1.5, font = 2)
> 
> 
> ##  Yong Lu <lyongu+@cs.cmu.edu> example -- a bit smaller
> x2 <- c(1, rep(2, 9))
> stopifnot(dip(x2) == dip(3 - x2))
> str(dip(x2, full=TRUE))
List of 15
 $ call       : language dip(x = x2, full.result = TRUE)
 $ x          : num [1:10] 1 2 2 2 2 2 2 2 2 2
 $ n          : int 10
 $ dip        : num 0.05
 $ lo.hi      : int [1:2] 2 10
 $ ifault     : int 0
 $ gcm        : int [1:2] 10 2
 $ lcm        : int [1:2] 2 10
 $ mn         : int [1:10] 1 1 2 2 2 2 2 2 2 2
 $ mj         : int [1:10] 10 10 10 10 10 10 10 10 10 10
 $ min.is.0   : logi FALSE
 $ debug      : int 0
 $ xl         : num 2
 $ xu         : num 2
 $ full.result: logi TRUE
 - attr(*, "class")= chr "dip"
> cat('Time elapsed: ', (.pt <- proc.time()),'\n') # "stats"
Time elapsed:  0.35 0.04 0.39 NA NA 
> 
> ## Real data examples :
> 
> data(statfaculty)
> 
> str(dip(statfaculty, full = "all", debug = 3), vec.len = 8)
dip() in C: n = 63; starting with  2N*dip = 1.
'dip': LOOP-BEGIN: 2n*D= 1         [low,high] = [  1, 63] :
 gcm[1:6] = 63, 62, 7, 3, 2, 1
 lcm[1:5] = 1, 44, 58, 59, 63
  while(gcm[ix] != lcm[iv]) :
 G(5,2) --> ix = 4, iv = 2
 G(4,2) --> ix = 3, iv = 2
 G(3,2) --> ix = 2, iv = 2
 L(3,2) --> ix = 2, iv = 3
 L(3,3) --> ix = 2, iv = 4
 --> ix = 2, iv = 5
 --> ix = 1, iv = 5
  calculating dip .. (dip_l, dip_u) = (2, 2.11111) -> new larger dip 2.11111 (j_best = 61)
'dip': LOOP-BEGIN: 2n*D= 2.1111    [low,high] = [  7, 58] :
 gcm[1:5] = 58, 55, 51, 48, 7
 lcm[1:6] = 7, 11, 15, 42, 44, 58
  while(gcm[ix] != lcm[iv]) :
 L(5,2) --> ix = 4, iv = 3
 L(5,3) --> ix = 4, iv = 4
 L(5,4) --> ix = 4, iv = 5
 L(5,5) --> ix = 4, iv = 6
 --> ix = 3, iv = 6
 --> ix = 2, iv = 6
 --> ix = 1, iv = 6
  calculating dip .. (dip_l, dip_u) = (0, 7.5) -> new larger dip 7.5 (j_best = 48)
'dip': LOOP-BEGIN: 2n*D= 7.5       [low,high] = [  7, 44] :
 gcm[1:4] = 44, 43, 38, 7
 lcm[1:5] = 7, 11, 15, 42, 44
  while(gcm[ix] != lcm[iv]) :
 L(4,2) --> ix = 3, iv = 3
 L(4,3) --> ix = 3, iv = 4
 --> ix = 2, iv = 4
 --> ix = 2, iv = 5
 --> ix = 1, iv = 5
List of 17
 $ call       : language dip(x = statfaculty, full.result = "all", debug = 3)
 $ x          : num [1:63] 30 33 35 36 37 37 39 39 39 39 39 40 40 40 40 41 42 43 43 43 ...
 $ n          : int 63
 $ dip        : num 0.0595
 $ lo.hi      : int [1:2] 7 44
 $ ifault     : int 0
 $ gcm        : int [1:4] 44 43 38 7
 $ lcm        : int [1:5] 7 11 15 42 44
 $ mn         : int [1:63] 1 1 2 3 3 5 3 7 7 7 7 7 12 12 12 7 7 7 18 18 ...
 $ mj         : int [1:63] 44 44 15 15 6 15 11 11 11 11 15 15 15 15 42 42 20 20 20 42 ...
 $ min.is.0   : logi FALSE
 $ debug      : int 3
 $ xl         : num 39
 $ xu         : num 54
 $ full.result: chr "all"
 $ GCM        : int [1:6] 63 62 7 3 2 1
 $ LCM        : int [1:5] 1 44 58 59 63
 - attr(*, "class")= chr "dip"
> 
> data(faithful)
> fE <- faithful$eruptions
> str(dip(fE, full = "all", debug = 3),
+     vec.len= 8)
dip() in C: n = 272; starting with  2N*dip = 1.
'dip': LOOP-BEGIN: 2n*D= 1         [low,high] = [  1,272] :
 gcm[1:7] = 272, 135, 120, 119, 4, 2, 1
 lcm[1:10] = 1, 40, 58, 60, 66, 79, 91, 261, 268, 272
  while(gcm[ix] != lcm[iv]) :
 G(6,2) --> ix = 5, iv = 2
 G(5,2) --> ix = 4, iv = 2
 L(5,2) --> ix = 4, iv = 3
 L(5,3) --> ix = 4, iv = 4
 L(5,4) --> ix = 4, iv = 5
 L(5,5) --> ix = 4, iv = 6
 L(5,6) --> ix = 4, iv = 7
 L(5,7) --> ix = 4, iv = 8
 G(4,8) --> ix = 3, iv = 8
 G(3,8) --> ix = 2, iv = 8
 --> ix = 1, iv = 8
 --> ix = 1, iv = 9
 --> ix = 1, iv = 10
  calculating dip .. (dip_l, dip_u) = (50.2553, 3) -> new larger dip 50.2553 (j_best = 91)
'dip': LOOP-BEGIN: 2n*D= 50.255    [low,high] = [120,261] :
 gcm[1:7] = 261, 260, 252, 181, 146, 135, 120
 lcm[1:5] = 120, 124, 233, 246, 261
  while(gcm[ix] != lcm[iv]) :
 L(7,2) --> ix = 6, iv = 3
 G(6,3) --> ix = 5, iv = 3
 G(5,3) --> ix = 4, iv = 3
 G(4,3) --> ix = 3, iv = 3
 --> ix = 3, iv = 4
 --> ix = 3, iv = 5
 --> ix = 2, iv = 5
 --> ix = 1, iv = 5
List of 17
 $ call       : language dip(x = fE, full.result = "all", debug = 3)
 $ x          : num [1:272] 1.6 1.67 1.7 1.73 1.75 1.75 1.75 1.75 1.75 1.75 ...
 $ n          : int 272
 $ dip        : num 0.0924
 $ lo.hi      : int [1:2] 120 261
 $ ifault     : int 0
 $ gcm        : int [1:7] 261 260 252 181 146 135 120
 $ lcm        : int [1:5] 120 124 233 246 261
 $ mn         : int [1:272] 1 1 2 2 4 5 5 5 5 5 5 11 5 13 13 13 13 17 17 13 ...
 $ mj         : int [1:272] 40 40 40 10 10 10 10 10 10 40 12 36 16 16 16 26 19 19 26 26 ...
 $ min.is.0   : logi FALSE
 $ debug      : int 3
 $ xl         : num 3.83
 $ xu         : num 4.83
 $ full.result: chr "all"
 $ GCM        : int [1:7] 272 135 120 119 4 2 1
 $ LCM        : int [1:10] 1 40 58 60 66 79 91 261 268 272
 - attr(*, "class")= chr "dip"
> 
> data(precip)
> str(dip(precip, full = TRUE, debug = TRUE))
dip() in C: n = 70; starting with  2N*dip = 1.
'dip': LOOP-BEGIN: 2n*D= 1         [low,high] = [  1, 70]; l_lcm/gcm = ( 6, 4)
  while(gcm[ix] != lcm[iv]) : ......
  calculating dip .. (dip_l, dip_u) = (5, 2.5) -> new larger dip 5 (j_best = 13)
'dip': LOOP-BEGIN: 2n*D= 5         [low,high] = [ 19, 64]; l_lcm/gcm = ( 6, 6)
  while(gcm[ix] != lcm[iv]) : ........
  calculating dip .. (dip_l, dip_u) = (3.875, 3.44828)
'dip': LOOP-BEGIN: 2n*D= 5         [low,high] = [ 31, 55]; l_lcm/gcm = ( 4, 3)
  while(gcm[ix] != lcm[iv]) : ...
List of 15
 $ call       : language dip(x = precip, full.result = TRUE, debug = TRUE)
 $ x          : num [1:70] 7 7.2 7.8 7.8 11.5 13 14 14.6 15 15.2 ...
 $ n          : int 70
 $ dip        : num 0.0357
 $ lo.hi      : int [1:2] 31 55
 $ ifault     : int 0
 $ gcm        : int [1:3] 55 49 31
 $ lcm        : int [1:4] 31 32 35 55
 $ mn         : int [1:70] 1 1 1 3 1 1 6 7 8 9 ...
 $ mj         : int [1:70] 2 4 4 64 55 10 10 10 10 55 ...
 $ min.is.0   : logi FALSE
 $ debug      : int 1
 $ xl         : Named num 35.9
  ..- attr(*, "names")= chr "Dallas"
 $ xu         : Named num 43.4
  ..- attr(*, "names")= chr "Hartford"
 $ full.result: logi TRUE
 - attr(*, "class")= chr "dip"
> 
> cat('Time elapsed: ', proc.time() - .pt,'\n') # "stats"
Time elapsed:  0.07 0.03 0.09 NA NA 
> 
> if(!interactive()) warnings()
> 
> proc.time()
   user  system elapsed 
   0.42    0.07    0.48