|
|
|
|
@ -1,732 +0,0 @@
|
|
|
|
|
From 40d377b86c856f5a4510a6f5cd56be004873ad77 Mon Sep 17 00:00:00 2001
|
|
|
|
|
From: =?UTF-8?q?Marcus=20M=C3=BCller?= <mueller@kit.edu>
|
|
|
|
|
Date: Mon, 12 Oct 2020 21:30:12 +0200
|
|
|
|
|
Subject: [PATCH] Remove defunct rgb_operm
|
|
|
|
|
|
|
|
|
|
[Retrieved from:
|
|
|
|
|
https://github.com/eddelbuettel/dieharder/pull/2/commits/40d377b86c856f5a4510a6f5cd56be004873ad77]
|
|
|
|
|
Signed-off-by: Fabrice Fontaine <fontaine.fabrice@gmail.com>
|
|
|
|
|
---
|
|
|
|
|
include/Makefile.am | 1 -
|
|
|
|
|
include/dieharder/rgb_operm.h | 38 --
|
|
|
|
|
include/dieharder/tests.h | 2 -
|
|
|
|
|
libdieharder/rgb_operm.c | 633 ----------------------------------
|
|
|
|
|
4 files changed, 674 deletions(-)
|
|
|
|
|
delete mode 100644 include/dieharder/rgb_operm.h
|
|
|
|
|
delete mode 100644 libdieharder/rgb_operm.c
|
|
|
|
|
|
|
|
|
|
diff --git a/include/Makefile.am b/include/Makefile.am
|
|
|
|
|
index f80b4ff..e4659cd 100644
|
|
|
|
|
--- a/include/Makefile.am
|
|
|
|
|
+++ b/include/Makefile.am
|
|
|
|
|
@@ -33,7 +33,6 @@ nobase_include_HEADERS = dieharder/copyright.h \
|
|
|
|
|
dieharder/rgb_lagged_sums.h \
|
|
|
|
|
dieharder/rgb_lmn.h \
|
|
|
|
|
dieharder/rgb_minimum_distance.h \
|
|
|
|
|
- dieharder/rgb_operm.h \
|
|
|
|
|
dieharder/rgb_persist.h \
|
|
|
|
|
dieharder/rgb_permutations.h \
|
|
|
|
|
dieharder/rgb_timing.h \
|
|
|
|
|
diff --git a/include/dieharder/rgb_operm.h b/include/dieharder/rgb_operm.h
|
|
|
|
|
deleted file mode 100644
|
|
|
|
|
index c48fa37..0000000
|
|
|
|
|
--- a/include/dieharder/rgb_operm.h
|
|
|
|
|
+++ /dev/null
|
|
|
|
|
@@ -1,38 +0,0 @@
|
|
|
|
|
-/*
|
|
|
|
|
- * rgb_operm test header.
|
|
|
|
|
- */
|
|
|
|
|
-
|
|
|
|
|
-/*
|
|
|
|
|
- * function prototype
|
|
|
|
|
- */
|
|
|
|
|
-int rgb_operm(Test **test,int irun);
|
|
|
|
|
-
|
|
|
|
|
-static Dtest rgb_operm_dtest __attribute__((unused)) = {
|
|
|
|
|
- "RGB Overlapping Permuations Test",
|
|
|
|
|
- "rgb_operm",
|
|
|
|
|
- "\n\
|
|
|
|
|
-#========================================================================\n\
|
|
|
|
|
-# RGB Overlapping Permutations Test\n\
|
|
|
|
|
-# Forms both the exact (expected) covariance matrix for overlapping\n\
|
|
|
|
|
-# permutations of random integer and an empirical covariance matrix\n\
|
|
|
|
|
-# formed from a long string of samples. The difference is expected\n\
|
|
|
|
|
-# to have a chisq distribution and hence can be transformed into a\n\
|
|
|
|
|
-# sample p-value. Note that this is one possible functional replacement\n\
|
|
|
|
|
-# for the broken/defunct diehard operm5 test, but one that permits k (the\n\
|
|
|
|
|
-# number of numbers in the overlapping permutation window) to be varied\n\
|
|
|
|
|
-# from 2 to perhaps 8.\n\
|
|
|
|
|
-#\n",
|
|
|
|
|
- 100, /* Default psamples */
|
|
|
|
|
- 100000, /* Default tsamples */
|
|
|
|
|
- 1, /* We magically make all the bit tests return a single histogram */
|
|
|
|
|
- rgb_operm,
|
|
|
|
|
- 0
|
|
|
|
|
-};
|
|
|
|
|
-
|
|
|
|
|
-/*
|
|
|
|
|
- * Global variables.
|
|
|
|
|
- *
|
|
|
|
|
- * rgb_operm_k is the size of the overlapping window that is slid along
|
|
|
|
|
- * a data stream of rands from x_i to x_{i+k} to compute c[][].
|
|
|
|
|
- */
|
|
|
|
|
-unsigned int rgb_operm_k;
|
|
|
|
|
diff --git a/include/dieharder/tests.h b/include/dieharder/tests.h
|
|
|
|
|
index 1674aed..b50dbe3 100644
|
|
|
|
|
--- a/include/dieharder/tests.h
|
|
|
|
|
+++ b/include/dieharder/tests.h
|
|
|
|
|
@@ -11,7 +11,6 @@
|
|
|
|
|
#include <dieharder/rgb_kstest_test.h>
|
|
|
|
|
#include <dieharder/rgb_lagged_sums.h>
|
|
|
|
|
#include <dieharder/rgb_minimum_distance.h>
|
|
|
|
|
-#include <dieharder/rgb_operm.h>
|
|
|
|
|
#include <dieharder/rgb_permutations.h>
|
|
|
|
|
#include <dieharder/dab_bytedistrib.h>
|
|
|
|
|
#include <dieharder/dab_dct.h>
|
|
|
|
|
@@ -80,7 +79,6 @@
|
|
|
|
|
RGB_PERMUTATIONS,
|
|
|
|
|
RGB_LAGGED_SUMS,
|
|
|
|
|
RGB_LMN,
|
|
|
|
|
- RGB_OPERM,
|
|
|
|
|
DAB_BYTEDISTRIB,
|
|
|
|
|
DAB_DCT,
|
|
|
|
|
DAB_FILLTREE,
|
|
|
|
|
diff --git a/libdieharder/rgb_operm.c b/libdieharder/rgb_operm.c
|
|
|
|
|
deleted file mode 100644
|
|
|
|
|
index 15f8e9a..0000000
|
|
|
|
|
--- a/libdieharder/rgb_operm.c
|
|
|
|
|
+++ /dev/null
|
|
|
|
|
@@ -1,633 +0,0 @@
|
|
|
|
|
-/*
|
|
|
|
|
- * ========================================================================
|
|
|
|
|
- * $Id: rgb_operm.c 252 2006-10-10 13:17:36Z rgb $
|
|
|
|
|
- *
|
|
|
|
|
- * See copyright in copyright.h and the accompanying file COPYING
|
|
|
|
|
- * ========================================================================
|
|
|
|
|
- */
|
|
|
|
|
-
|
|
|
|
|
-/*
|
|
|
|
|
- * ========================================================================
|
|
|
|
|
- * This is the revised Overlapping Permutations test. It directly
|
|
|
|
|
- * simulates the covariance matrix of overlapping permutations. The way
|
|
|
|
|
- * this works below (tentatively) is:
|
|
|
|
|
- *
|
|
|
|
|
- * For a bit ntuple of length N, slide a window of length N to the
|
|
|
|
|
- * right one bit at a time. Compute the permutation index of the
|
|
|
|
|
- * original ntuple, the permutation index of the window ntuple, and
|
|
|
|
|
- * accumulate the covariance matrix of the two positions. This
|
|
|
|
|
- * can be directly and precisely computed as well. The simulated
|
|
|
|
|
- * result should be distributed according to the chisq distribution,
|
|
|
|
|
- * so we subtract the two and feed it into the chisq program as a
|
|
|
|
|
- * vector to compute p.
|
|
|
|
|
- *
|
|
|
|
|
- * This MAY NOT BE RIGHT. I'm working from both Marsaglia's limited
|
|
|
|
|
- * documentation (in a program that doesn't do ANYTHING like what the
|
|
|
|
|
- * documentation says it does) and from Nilpotent Markov Processes.
|
|
|
|
|
- * But I confess to not quite understand how to actually perform the
|
|
|
|
|
- * test in the latter -- it is very good at describing the construction
|
|
|
|
|
- * of the target matrix, not so good at describing how to transform
|
|
|
|
|
- * this into a chisq and p.
|
|
|
|
|
- *
|
|
|
|
|
- * FWIW, as I get something that actually works here, I'm going to
|
|
|
|
|
- * THOROUGHLY document it in the book that will accompany the test.
|
|
|
|
|
- *========================================================================
|
|
|
|
|
- */
|
|
|
|
|
-
|
|
|
|
|
-#include <dieharder/libdieharder.h>
|
|
|
|
|
-#define RGB_OPERM_KMAX 10
|
|
|
|
|
-
|
|
|
|
|
-/*
|
|
|
|
|
- * Some globals that will eventually go in the test include where they
|
|
|
|
|
- * arguably belong.
|
|
|
|
|
- */
|
|
|
|
|
-double fpipi(int pi1,int pi2,int nkp);
|
|
|
|
|
-uint piperm(size_t *data,int len);
|
|
|
|
|
-void make_cexact();
|
|
|
|
|
-void make_cexpt();
|
|
|
|
|
-int nperms,noperms;
|
|
|
|
|
-double **cexact,**ceinv,**cexpt,**idty;
|
|
|
|
|
-double *cvexact,*cvein,*cvexpt,*vidty;
|
|
|
|
|
-
|
|
|
|
|
-int rgb_operm(Test **test,int irun)
|
|
|
|
|
-{
|
|
|
|
|
-
|
|
|
|
|
- int i,j,n,nb,iv,s;
|
|
|
|
|
- uint csamples; /* rgb_operm_k^2 is vector size of cov matrix */
|
|
|
|
|
- uint *count,ctotal; /* counters */
|
|
|
|
|
- uint size;
|
|
|
|
|
- double pvalue,ntuple_prob,pbin; /* probabilities */
|
|
|
|
|
- Vtest *vtest; /* Chisq entry vector */
|
|
|
|
|
-
|
|
|
|
|
- gsl_matrix_view CEXACT,CEINV,CEXPT,IDTY;
|
|
|
|
|
-
|
|
|
|
|
- /*
|
|
|
|
|
- * For a given n = ntuple size in bits, there are n! bit orderings
|
|
|
|
|
- */
|
|
|
|
|
- MYDEBUG(D_RGB_OPERM){
|
|
|
|
|
- printf("#==================================================================\n");
|
|
|
|
|
- printf("# rgb_operm: Running rgb_operm verbosely for k = %d.\n",rgb_operm_k);
|
|
|
|
|
- printf("# rgb_operm: Use -v = %d to focus.\n",D_RGB_OPERM);
|
|
|
|
|
- printf("# rgb_operm: ======================================================\n");
|
|
|
|
|
- }
|
|
|
|
|
-
|
|
|
|
|
- /*
|
|
|
|
|
- * Sanity check first
|
|
|
|
|
- */
|
|
|
|
|
- if((rgb_operm_k < 0) || (rgb_operm_k > RGB_OPERM_KMAX)){
|
|
|
|
|
- printf("\nError: rgb_operm_k must be a positive integer <= %u. Exiting.\n",RGB_OPERM_KMAX);
|
|
|
|
|
- exit(0);
|
|
|
|
|
- }
|
|
|
|
|
-
|
|
|
|
|
- nperms = gsl_sf_fact(rgb_operm_k);
|
|
|
|
|
- noperms = gsl_sf_fact(3*rgb_operm_k-2);
|
|
|
|
|
- csamples = rgb_operm_k*rgb_operm_k;
|
|
|
|
|
- gsl_permutation * p = gsl_permutation_alloc(nperms);
|
|
|
|
|
-
|
|
|
|
|
- /*
|
|
|
|
|
- * Allocate memory for value_max vector of Vtest structs and counts,
|
|
|
|
|
- * PER TEST. Note that we must free both of these when we are done
|
|
|
|
|
- * or leak.
|
|
|
|
|
- */
|
|
|
|
|
- vtest = (Vtest *)malloc(csamples*sizeof(Vtest));
|
|
|
|
|
- count = (uint *)malloc(csamples*sizeof(uint));
|
|
|
|
|
- Vtest_create(vtest,csamples+1);
|
|
|
|
|
-
|
|
|
|
|
- /*
|
|
|
|
|
- * We have to allocate and free the cexact and cexpt matrices here
|
|
|
|
|
- * or they'll be forgotten when these routines return.
|
|
|
|
|
- */
|
|
|
|
|
- MYDEBUG(D_RGB_OPERM){
|
|
|
|
|
- printf("# rgb_operm: Creating and zeroing cexact[][] and cexpt[][].\n");
|
|
|
|
|
- }
|
|
|
|
|
- cexact = (double **)malloc(nperms*sizeof(double*));
|
|
|
|
|
- ceinv = (double **)malloc(nperms*sizeof(double*));
|
|
|
|
|
- cexpt = (double **)malloc(nperms*sizeof(double*));
|
|
|
|
|
- idty = (double **)malloc(nperms*sizeof(double*));
|
|
|
|
|
- cvexact = (double *)malloc(nperms*nperms*sizeof(double));
|
|
|
|
|
- cvein = (double *)malloc(nperms*nperms*sizeof(double));
|
|
|
|
|
- cvexpt = (double *)malloc(nperms*nperms*sizeof(double));
|
|
|
|
|
- vidty = (double *)malloc(nperms*nperms*sizeof(double));
|
|
|
|
|
- for(i=0;i<nperms;i++){
|
|
|
|
|
- /* Here we pack addresses to map the matrix addressing onto the vector */
|
|
|
|
|
- cexact[i] = &cvexact[i*nperms];
|
|
|
|
|
- ceinv[i] = &cvein[i*nperms];
|
|
|
|
|
- cexpt[i] = &cvexpt[i*nperms];
|
|
|
|
|
- idty[i] = &vidty[i*nperms];
|
|
|
|
|
- for(j = 0;j<nperms;j++){
|
|
|
|
|
- cexact[i][j] = 0.0;
|
|
|
|
|
- ceinv[i][j] = 0.0;
|
|
|
|
|
- cexpt[i][j] = 0.0;
|
|
|
|
|
- idty[i][j] = 0.0;
|
|
|
|
|
- }
|
|
|
|
|
- }
|
|
|
|
|
-
|
|
|
|
|
- make_cexact();
|
|
|
|
|
- make_cexpt();
|
|
|
|
|
-
|
|
|
|
|
- iv=0;
|
|
|
|
|
- for(i=0;i<nperms;i++){
|
|
|
|
|
- for(j=0;j<nperms;j++){
|
|
|
|
|
- cvexact[iv] = cexact[i][j];
|
|
|
|
|
- cvexpt[iv] = cexpt[i][j];
|
|
|
|
|
- vidty[iv] = 0.0;
|
|
|
|
|
- }
|
|
|
|
|
- }
|
|
|
|
|
-
|
|
|
|
|
- CEXACT = gsl_matrix_view_array(cvexact, nperms, nperms);
|
|
|
|
|
- CEINV = gsl_matrix_view_array(cvein , nperms, nperms);
|
|
|
|
|
- CEXPT = gsl_matrix_view_array(cvexpt , nperms, nperms);
|
|
|
|
|
- IDTY = gsl_matrix_view_array(vidty , nperms, nperms);
|
|
|
|
|
-
|
|
|
|
|
- /*
|
|
|
|
|
- * Hmmm, looks like cexact isn't invertible. Duh. So it has eigenvalues.
|
|
|
|
|
- * This seems to be important (how, I do not know) so let's find out.
|
|
|
|
|
- * Here is the gsl ritual for evaluating eigenvalues etc.
|
|
|
|
|
- */
|
|
|
|
|
-
|
|
|
|
|
- gsl_vector *eval = gsl_vector_alloc (nperms);
|
|
|
|
|
- gsl_matrix *evec = gsl_matrix_alloc (nperms,nperms);
|
|
|
|
|
- /*
|
|
|
|
|
- gsl_eigen_nonsymm_workspace* w = gsl_eigen_nonsymmv_alloc(nperms);
|
|
|
|
|
- gsl_eigen_nonsymm_params (1,0,w);
|
|
|
|
|
- gsl_eigen_nonsymmv(&CEXACT.matrix, eval, evec, w);
|
|
|
|
|
- gsl_eigen_nonsymmv_free (w);
|
|
|
|
|
- */
|
|
|
|
|
- gsl_eigen_symmv_workspace* w = gsl_eigen_symmv_alloc(nperms);
|
|
|
|
|
- gsl_eigen_symmv(&CEXACT.matrix, eval, evec, w);
|
|
|
|
|
- gsl_eigen_symmv_free (w);
|
|
|
|
|
- gsl_eigen_symmv_sort (eval, evec, GSL_EIGEN_SORT_ABS_ASC);
|
|
|
|
|
-
|
|
|
|
|
- {
|
|
|
|
|
- int i;
|
|
|
|
|
-
|
|
|
|
|
- printf("#==================================================================\n");
|
|
|
|
|
- for (i = 0; i < nperms; i++) {
|
|
|
|
|
- double eval_i = gsl_vector_get (eval, i);
|
|
|
|
|
- gsl_vector_view evec_i = gsl_matrix_column (evec, i);
|
|
|
|
|
- printf ("eigenvalue[%u] = %g\n", i, eval_i);
|
|
|
|
|
- printf ("eigenvector[%u] = \n",i);
|
|
|
|
|
- gsl_vector_fprintf (stdout,&evec_i.vector, "%10.5f");
|
|
|
|
|
- }
|
|
|
|
|
- printf("#==================================================================\n");
|
|
|
|
|
- }
|
|
|
|
|
-
|
|
|
|
|
- gsl_vector_free (eval);
|
|
|
|
|
- gsl_matrix_free (evec);
|
|
|
|
|
-
|
|
|
|
|
-/*
|
|
|
|
|
- gsl_linalg_LU_decomp(&CEXACT.matrix, p, &s);
|
|
|
|
|
- gsl_linalg_LU_invert(&CEXACT, p, &CEINV);
|
|
|
|
|
- gsl_permutation_free(p);
|
|
|
|
|
- gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, &CEINV.matrix, &CEXPT.matrix, 0.0, &IDTY.matrix);
|
|
|
|
|
- printf("#==================================================================\n");
|
|
|
|
|
- printf("# Should be inverse of C, assuming it is invertible:\n");
|
|
|
|
|
- for(i=0;i<nperms;i++){
|
|
|
|
|
- printf("# ");
|
|
|
|
|
- for(j = 0;j<nperms;j++){
|
|
|
|
|
- printf("%8.3f ",idty[i][j]);
|
|
|
|
|
- }
|
|
|
|
|
- printf("\n");
|
|
|
|
|
- }
|
|
|
|
|
- printf("#==================================================================\n");
|
|
|
|
|
- printf("#==================================================================\n");
|
|
|
|
|
- printf("# Should be normal on identity:\n");
|
|
|
|
|
- for(i=0;i<nperms;i++){
|
|
|
|
|
- printf("# ");
|
|
|
|
|
- for(j = 0;j<nperms;j++){
|
|
|
|
|
- printf("%8.3f ",idty[i][j]);
|
|
|
|
|
- }
|
|
|
|
|
- printf("\n");
|
|
|
|
|
- }
|
|
|
|
|
- printf("#==================================================================\n");
|
|
|
|
|
- */
|
|
|
|
|
-
|
|
|
|
|
-
|
|
|
|
|
-
|
|
|
|
|
- /*
|
|
|
|
|
- * OK, at this point we have two matrices: cexact[][] is filled with
|
|
|
|
|
- * the exact covariance matrix expected for the overlapping permutations.
|
|
|
|
|
- * cexpt[][] has been filled numerically by generating strings of random
|
|
|
|
|
- * uints or floats, generating sort index permutations, and
|
|
|
|
|
- * using them to IDENTICALLY generate an "experimental" version of c[][].
|
|
|
|
|
- * The two should correspond, in the limit of large tsamples. IF I
|
|
|
|
|
- * understand Alhakim, Kawczak and Molchanov, then the way to implement
|
|
|
|
|
- * the simplest possible chisq test is to evaluate:
|
|
|
|
|
- * cexact^-1 cexpt \approx I
|
|
|
|
|
- * where the diagonal terms should form a vector that is chisq distributed?
|
|
|
|
|
- * Let's try this...
|
|
|
|
|
- */
|
|
|
|
|
-
|
|
|
|
|
-
|
|
|
|
|
-
|
|
|
|
|
- /*
|
|
|
|
|
- * Free cexact[][] and cexpt[][]
|
|
|
|
|
- * Fix this when we're done so we don't leak; for now to much trouble.
|
|
|
|
|
- for(i=0;i<nperms;i++){
|
|
|
|
|
- free(cexact[i]);
|
|
|
|
|
- free(cexpt[i]);
|
|
|
|
|
- }
|
|
|
|
|
- free(cexact);
|
|
|
|
|
- free(cexpt);
|
|
|
|
|
- */
|
|
|
|
|
-
|
|
|
|
|
- return(0);
|
|
|
|
|
-
|
|
|
|
|
-}
|
|
|
|
|
-
|
|
|
|
|
-void make_cexact()
|
|
|
|
|
-{
|
|
|
|
|
-
|
|
|
|
|
- int i,j,k,ip,t,nop;
|
|
|
|
|
- double fi,fj;
|
|
|
|
|
- /*
|
|
|
|
|
- * This is the test vector.
|
|
|
|
|
- */
|
|
|
|
|
- double testv[RGB_OPERM_KMAX*2]; /* easier than malloc etc, but beware length */
|
|
|
|
|
- /*
|
|
|
|
|
- * pi[] is the permutation index of a sample. ps[] holds the
|
|
|
|
|
- * actual sample.
|
|
|
|
|
- */
|
|
|
|
|
- size_t pi[4096],ps[4096];
|
|
|
|
|
- /*
|
|
|
|
|
- * We seem to have made a mistake of sorts. We actually have to sum
|
|
|
|
|
- * BOTH the forward AND the backward directions. That means that the
|
|
|
|
|
- * permutation vector has to be of length 3k-1, with the pi=1 term
|
|
|
|
|
- * corresponding to the middle. So for k=2, instead of 0,1,2 we need
|
|
|
|
|
- * 0 1 2 3 4 and we'll have to do 23, 34 in the leading direction and
|
|
|
|
|
- * 21, 10 in the trailing direction.
|
|
|
|
|
- */
|
|
|
|
|
- gsl_permutation **operms;
|
|
|
|
|
-
|
|
|
|
|
- MYDEBUG(D_RGB_OPERM){
|
|
|
|
|
- printf("#==================================================================\n");
|
|
|
|
|
- printf("# rgb_operm: Running cexact()\n");
|
|
|
|
|
- }
|
|
|
|
|
-
|
|
|
|
|
- /*
|
|
|
|
|
- * Test fpipi(). This is probably cruft, actually.
|
|
|
|
|
- MYDEBUG(D_RGB_OPERM){
|
|
|
|
|
- printf("# rgb_operm: Testing fpipi()\n");
|
|
|
|
|
- for(i=0;i<nperms;i++){
|
|
|
|
|
- for(j = 0;j<nperms;j++){
|
|
|
|
|
- printf("# rgb_operm: fpipi(%u,%u,%u) = %f\n",i,j,nperms,fpipi(i,j,nperms));
|
|
|
|
|
- }
|
|
|
|
|
- }
|
|
|
|
|
- }
|
|
|
|
|
- */
|
|
|
|
|
-
|
|
|
|
|
- MYDEBUG(D_RGB_OPERM){
|
|
|
|
|
- printf("#==================================================================\n");
|
|
|
|
|
- printf("# rgb_operm: Forming set of %u overlapping permutations\n",noperms);
|
|
|
|
|
- printf("# rgb_operm: Permutations\n");
|
|
|
|
|
- printf("# rgb_operm:==============================\n");
|
|
|
|
|
- }
|
|
|
|
|
- operms = (gsl_permutation**) malloc(noperms*sizeof(gsl_permutation*));
|
|
|
|
|
- for(i=0;i<noperms;i++){
|
|
|
|
|
- operms[i] = gsl_permutation_alloc(3*rgb_operm_k - 2);
|
|
|
|
|
- /* Must quiet down
|
|
|
|
|
- MYDEBUG(D_RGB_OPERM){
|
|
|
|
|
- printf("# rgb_operm: ");
|
|
|
|
|
- }
|
|
|
|
|
- */
|
|
|
|
|
- if(i == 0){
|
|
|
|
|
- gsl_permutation_init(operms[i]);
|
|
|
|
|
- } else {
|
|
|
|
|
- gsl_permutation_memcpy(operms[i],operms[i-1]);
|
|
|
|
|
- gsl_permutation_next(operms[i]);
|
|
|
|
|
- }
|
|
|
|
|
- /*
|
|
|
|
|
- MYDEBUG(D_RGB_OPERM){
|
|
|
|
|
- gsl_permutation_fprintf(stdout,operms[i]," %u");
|
|
|
|
|
- printf("\n");
|
|
|
|
|
- }
|
|
|
|
|
- */
|
|
|
|
|
- }
|
|
|
|
|
-
|
|
|
|
|
- /*
|
|
|
|
|
- * We now form c_exact PRECISELY the same way that we do c_expt[][]
|
|
|
|
|
- * below, except that instead of pulling random samples of integers
|
|
|
|
|
- * or floats and averaging over the permutations thus represented,
|
|
|
|
|
- * we iterate over the complete set of equally weighted permutations
|
|
|
|
|
- * to get an exact answer. Note that we have to center on 2k-1 and
|
|
|
|
|
- * go both forwards and backwards.
|
|
|
|
|
- */
|
|
|
|
|
- for(t=0;t<noperms;t++){
|
|
|
|
|
- /*
|
|
|
|
|
- * To sort into a perm, test vector needs to be double.
|
|
|
|
|
- */
|
|
|
|
|
- for(k=0;k<3*rgb_operm_k - 2;k++) testv[k] = (double) operms[t]->data[k];
|
|
|
|
|
-
|
|
|
|
|
- /* Not cruft, but quiet...
|
|
|
|
|
- MYDEBUG(D_RGB_OPERM){
|
|
|
|
|
- printf("#------------------------------------------------------------------\n");
|
|
|
|
|
- printf("# Generating offset sample permutation pi's\n");
|
|
|
|
|
- }
|
|
|
|
|
- */
|
|
|
|
|
- for(k=0;k<2*rgb_operm_k - 1;k++){
|
|
|
|
|
- gsl_sort_index((size_t *) ps,&testv[k],1,rgb_operm_k);
|
|
|
|
|
- pi[k] = piperm((size_t *) ps,rgb_operm_k);
|
|
|
|
|
-
|
|
|
|
|
- /* Not cruft, but quiet...
|
|
|
|
|
- MYDEBUG(D_RGB_OPERM){
|
|
|
|
|
- printf("# %u: ",k);
|
|
|
|
|
- for(ip=k;ip<rgb_operm_k+k;ip++){
|
|
|
|
|
- printf("%.1f ",testv[ip]);
|
|
|
|
|
- }
|
|
|
|
|
- printf("\n# ");
|
|
|
|
|
- for(ip=0;ip<rgb_operm_k;ip++){
|
|
|
|
|
- printf("%u ",ps[ip]);
|
|
|
|
|
- }
|
|
|
|
|
- printf(" = %u\n",pi[k]);
|
|
|
|
|
- }
|
|
|
|
|
- */
|
|
|
|
|
-
|
|
|
|
|
- }
|
|
|
|
|
-
|
|
|
|
|
- /*
|
|
|
|
|
- * This is the business end of things. The covariance matrix is the
|
|
|
|
|
- * the sum of a central function of the permutation indices that yields
|
|
|
|
|
- * nperms-1/nperms on diagonal, -1/nperms off diagonal, for all the
|
|
|
|
|
- * possible permutations, for the FIRST permutation in a sample (fi)
|
|
|
|
|
- * times the sum of the same function over all the overlapping permutations
|
|
|
|
|
- * drawn from the same sample. Quite simple, really.
|
|
|
|
|
- */
|
|
|
|
|
- for(i=0;i<nperms;i++){
|
|
|
|
|
- fi = fpipi(i,pi[rgb_operm_k-1],nperms);
|
|
|
|
|
- for(j=0;j<nperms;j++){
|
|
|
|
|
- fj = 0.0;
|
|
|
|
|
- for(k=0;k<rgb_operm_k;k++){
|
|
|
|
|
- fj += fpipi(j,pi[rgb_operm_k - 1 + k],nperms);
|
|
|
|
|
- if(k != 0){
|
|
|
|
|
- fj += fpipi(j,pi[rgb_operm_k - 1 - k],nperms);
|
|
|
|
|
- }
|
|
|
|
|
- }
|
|
|
|
|
- cexact[i][j] += fi*fj;
|
|
|
|
|
- }
|
|
|
|
|
- }
|
|
|
|
|
-
|
|
|
|
|
- }
|
|
|
|
|
-
|
|
|
|
|
- MYDEBUG(D_RGB_OPERM){
|
|
|
|
|
- printf("# rgb_operm:==============================\n");
|
|
|
|
|
- printf("# rgb_operm: cexact[][] = \n");
|
|
|
|
|
- }
|
|
|
|
|
- for(i=0;i<nperms;i++){
|
|
|
|
|
- MYDEBUG(D_RGB_OPERM){
|
|
|
|
|
- printf("# ");
|
|
|
|
|
- }
|
|
|
|
|
- for(j=0;j<nperms;j++){
|
|
|
|
|
- cexact[i][j] /= noperms;
|
|
|
|
|
- MYDEBUG(D_RGB_OPERM){
|
|
|
|
|
- printf("%10.6f ",cexact[i][j]);
|
|
|
|
|
- }
|
|
|
|
|
- }
|
|
|
|
|
- MYDEBUG(D_RGB_OPERM){
|
|
|
|
|
- printf("\n");
|
|
|
|
|
- }
|
|
|
|
|
- }
|
|
|
|
|
-
|
|
|
|
|
- /*
|
|
|
|
|
- * Free operms[]
|
|
|
|
|
- */
|
|
|
|
|
- for(i=0;i<noperms;i++){
|
|
|
|
|
- gsl_permutation_free(operms[i]);
|
|
|
|
|
- }
|
|
|
|
|
- free(operms);
|
|
|
|
|
-
|
|
|
|
|
-}
|
|
|
|
|
-
|
|
|
|
|
-void make_cexpt()
|
|
|
|
|
-{
|
|
|
|
|
-
|
|
|
|
|
- int i,j,k,ip,t;
|
|
|
|
|
- double fi,fj;
|
|
|
|
|
- /*
|
|
|
|
|
- * This is the test vector.
|
|
|
|
|
- */
|
|
|
|
|
- double testv[RGB_OPERM_KMAX*2]; /* easier than malloc etc, but beware length */
|
|
|
|
|
- /*
|
|
|
|
|
- * pi[] is the permutation index of a sample. ps[] holds the
|
|
|
|
|
- * actual sample.
|
|
|
|
|
- */
|
|
|
|
|
- int pi[4096],ps[4096];
|
|
|
|
|
-
|
|
|
|
|
- MYDEBUG(D_RGB_OPERM){
|
|
|
|
|
- printf("#==================================================================\n");
|
|
|
|
|
- printf("# rgb_operm: Running cexpt()\n");
|
|
|
|
|
- }
|
|
|
|
|
-
|
|
|
|
|
- /*
|
|
|
|
|
- * We evaluate cexpt[][] by sampling. In a nutshell, this involves
|
|
|
|
|
- * a) Filling testv[] with 2*rgb_operm_k - 1 random uints or doubles
|
|
|
|
|
- * It clearly cannot matter which we use, as long as the probability of
|
|
|
|
|
- * exact duplicates in a sample is very low.
|
|
|
|
|
- * b) Using gsl_sort_index the exact same way it was used in make_cexact()
|
|
|
|
|
- * to generate the pi[] index, using ps[] as scratch space for the sort
|
|
|
|
|
- * indices.
|
|
|
|
|
- * c) Evaluating fi and fj from the SAMPLED result, tsamples times.
|
|
|
|
|
- * d) Normalizing.
|
|
|
|
|
- * Note that this is pretty much identical to the way we formed c_exact[][]
|
|
|
|
|
- * except that we are determining the relative frequency of each sort order
|
|
|
|
|
- * permutation 2*rgb_operm_k-1 long.
|
|
|
|
|
- *
|
|
|
|
|
- * NOTE WELL! I honestly think that it is borderline silly to view
|
|
|
|
|
- * this as a matrix and to go through all of this nonsense. The theoretical
|
|
|
|
|
- * c_exact[][] is computed from the observation that all the permutations
|
|
|
|
|
- * of n objects have equal weight = 1/n!. Consequently, they should
|
|
|
|
|
- * individually be binomially distributed, tending to normal with many
|
|
|
|
|
- * samples. Collectively they should be distributed like a vector of
|
|
|
|
|
- * equal binomial probabilities and a p-value should follow either from
|
|
|
|
|
- * chisq on n!-1 DoF or for that matter a KS test. I see no way that
|
|
|
|
|
- * making it into a matrix can increase the sensitivity of the test -- if
|
|
|
|
|
- * the p-values are well defined in the two cases they can only be equal
|
|
|
|
|
- * by their very definition.
|
|
|
|
|
- *
|
|
|
|
|
- * If you are a statistician reading these words and disagree, please
|
|
|
|
|
- * communicate with me and explain why I'm wrong. I'm still very much
|
|
|
|
|
- * learning statistics and would cherish gentle correction.
|
|
|
|
|
- */
|
|
|
|
|
- for(t=0;t<tsamples;t++){
|
|
|
|
|
- /*
|
|
|
|
|
- * To sort into a perm, test vector needs to be double.
|
|
|
|
|
- */
|
|
|
|
|
- for(k=0;k<3*rgb_operm_k - 2;k++) testv[k] = (double) gsl_rng_get(rng);
|
|
|
|
|
-
|
|
|
|
|
- /* Not cruft, but quiet...
|
|
|
|
|
- MYDEBUG(D_RGB_OPERM){
|
|
|
|
|
- printf("#------------------------------------------------------------------\n");
|
|
|
|
|
- printf("# Generating offset sample permutation pi's\n");
|
|
|
|
|
- }
|
|
|
|
|
- */
|
|
|
|
|
- for(k=0;k<2*rgb_operm_k-1;k++){
|
|
|
|
|
- gsl_sort_index(ps,&testv[k],1,rgb_operm_k);
|
|
|
|
|
- pi[k] = piperm(ps,rgb_operm_k);
|
|
|
|
|
-
|
|
|
|
|
- /* Not cruft, but quiet...
|
|
|
|
|
- MYDEBUG(D_RGB_OPERM){
|
|
|
|
|
- printf("# %u: ",k);
|
|
|
|
|
- for(ip=k;ip<rgb_operm_k+k;ip++){
|
|
|
|
|
- printf("%.1f ",testv[ip]);
|
|
|
|
|
- }
|
|
|
|
|
- printf("\n# ");
|
|
|
|
|
- for(ip=0;ip<rgb_operm_k;ip++){
|
|
|
|
|
- printf("%u ",permsample->data[ip]);
|
|
|
|
|
- }
|
|
|
|
|
- printf(" = %u\n",pi[k]);
|
|
|
|
|
- }
|
|
|
|
|
- */
|
|
|
|
|
- }
|
|
|
|
|
-
|
|
|
|
|
- /*
|
|
|
|
|
- * This is the business end of things. The covariance matrix is the
|
|
|
|
|
- * the sum of a central function of the permutation indices that yields
|
|
|
|
|
- * nperms-1/nperms on diagonal, -1/nperms off diagonal, for all the
|
|
|
|
|
- * possible permutations, for the FIRST permutation in a sample (fi)
|
|
|
|
|
- * times the sum of the same function over all the overlapping permutations
|
|
|
|
|
- * drawn from the same sample. Quite simple, really.
|
|
|
|
|
- */
|
|
|
|
|
- for(i=0;i<nperms;i++){
|
|
|
|
|
- fi = fpipi(i,pi[rgb_operm_k-1],nperms);
|
|
|
|
|
- for(j=0;j<nperms;j++){
|
|
|
|
|
- fj = 0.0;
|
|
|
|
|
- for(k=0;k<rgb_operm_k;k++){
|
|
|
|
|
- fj += fpipi(j,pi[rgb_operm_k - 1 + k],nperms);
|
|
|
|
|
- if(k != 0){
|
|
|
|
|
- fj += fpipi(j,pi[rgb_operm_k - 1 - k],nperms);
|
|
|
|
|
- }
|
|
|
|
|
- }
|
|
|
|
|
- cexpt[i][j] += fi*fj;
|
|
|
|
|
- }
|
|
|
|
|
- }
|
|
|
|
|
-
|
|
|
|
|
- }
|
|
|
|
|
-
|
|
|
|
|
- MYDEBUG(D_RGB_OPERM){
|
|
|
|
|
- printf("# rgb_operm:==============================\n");
|
|
|
|
|
- printf("# rgb_operm: cexpt[][] = \n");
|
|
|
|
|
- }
|
|
|
|
|
- for(i=0;i<nperms;i++){
|
|
|
|
|
- MYDEBUG(D_RGB_OPERM){
|
|
|
|
|
- printf("# ");
|
|
|
|
|
- }
|
|
|
|
|
- for(j=0;j<nperms;j++){
|
|
|
|
|
- cexpt[i][j] /= tsamples;
|
|
|
|
|
- MYDEBUG(D_RGB_OPERM){
|
|
|
|
|
- printf("%10.6f ",cexpt[i][j]);
|
|
|
|
|
- }
|
|
|
|
|
- }
|
|
|
|
|
- MYDEBUG(D_RGB_OPERM){
|
|
|
|
|
- printf("\n");
|
|
|
|
|
- }
|
|
|
|
|
- }
|
|
|
|
|
-
|
|
|
|
|
-}
|
|
|
|
|
-
|
|
|
|
|
-uint piperm(size_t *data,int len)
|
|
|
|
|
-{
|
|
|
|
|
-
|
|
|
|
|
- uint i,j,k,max,min;
|
|
|
|
|
- uint pindex,uret,tmp;
|
|
|
|
|
- static gsl_permutation** lookup = 0;
|
|
|
|
|
-
|
|
|
|
|
- /*
|
|
|
|
|
- * Allocate space for lookup table and fill it.
|
|
|
|
|
- */
|
|
|
|
|
- if(lookup == 0){
|
|
|
|
|
- lookup = (gsl_permutation**) malloc(nperms*sizeof(gsl_permutation*));
|
|
|
|
|
- MYDEBUG(D_RGB_OPERM){
|
|
|
|
|
- printf("# rgb_operm: Allocating piperm lookup table of perms.\n");
|
|
|
|
|
- }
|
|
|
|
|
- for(i=0;i<nperms;i++){
|
|
|
|
|
- lookup[i] = gsl_permutation_alloc(rgb_operm_k);
|
|
|
|
|
- }
|
|
|
|
|
- for(i=0;i<nperms;i++){
|
|
|
|
|
- if(i == 0){
|
|
|
|
|
- gsl_permutation_init(lookup[i]);
|
|
|
|
|
- } else {
|
|
|
|
|
- gsl_permutation_memcpy(lookup[i],lookup[i-1]);
|
|
|
|
|
- gsl_permutation_next(lookup[i]);
|
|
|
|
|
- }
|
|
|
|
|
- }
|
|
|
|
|
-
|
|
|
|
|
- /*
|
|
|
|
|
- * This method yields a mirror symmetry in the permutations top to
|
|
|
|
|
- * bottom.
|
|
|
|
|
- for(i=0;i<nperms/2;i++){
|
|
|
|
|
- if(i == 0){
|
|
|
|
|
- gsl_permutation_init(lookup[i]);
|
|
|
|
|
- for(j=0;j<rgb_operm_k;j++){
|
|
|
|
|
- lookup[nperms-i-1]->data[rgb_operm_k-j-1] = lookup[i]->data[j];
|
|
|
|
|
- }
|
|
|
|
|
- } else {
|
|
|
|
|
- gsl_permutation_memcpy(lookup[i],lookup[i-1]);
|
|
|
|
|
- gsl_permutation_next(lookup[i]);
|
|
|
|
|
- for(j=0;j<rgb_operm_k;j++){
|
|
|
|
|
- lookup[nperms-i-1]->data[rgb_operm_k-j-1] = lookup[i]->data[j];
|
|
|
|
|
- }
|
|
|
|
|
- }
|
|
|
|
|
- }
|
|
|
|
|
- */
|
|
|
|
|
- MYDEBUG(D_RGB_OPERM){
|
|
|
|
|
- for(i=0;i<nperms;i++){
|
|
|
|
|
- printf("# rgb_operm: %u => ",i);
|
|
|
|
|
- gsl_permutation_fprintf(stdout,lookup[i]," %u");
|
|
|
|
|
- printf("\n");
|
|
|
|
|
- }
|
|
|
|
|
- }
|
|
|
|
|
-
|
|
|
|
|
- }
|
|
|
|
|
-
|
|
|
|
|
- for(i=0;i<nperms;i++){
|
|
|
|
|
- if(memcmp(data,lookup[i]->data,len*sizeof(uint))==0){
|
|
|
|
|
- /* Not cruft, but off:
|
|
|
|
|
- MYDEBUG(D_RGB_OPERM){
|
|
|
|
|
- printf("# piperm(): ");
|
|
|
|
|
- gsl_permutation_fprintf(stdout,lookup[i]," %u");
|
|
|
|
|
- printf(" = %u\n",i);
|
|
|
|
|
- }
|
|
|
|
|
- */
|
|
|
|
|
- return(i);
|
|
|
|
|
- }
|
|
|
|
|
- }
|
|
|
|
|
- printf("We'd better not get here...\n");
|
|
|
|
|
-
|
|
|
|
|
- return(0);
|
|
|
|
|
-
|
|
|
|
|
-}
|
|
|
|
|
-
|
|
|
|
|
-double fpipi(int pi1,int pi2,int nkp)
|
|
|
|
|
-{
|
|
|
|
|
-
|
|
|
|
|
- int i;
|
|
|
|
|
- double fret;
|
|
|
|
|
-
|
|
|
|
|
- /*
|
|
|
|
|
- * compute the k-permutation index from iperm for the window
|
|
|
|
|
- * at data[offset] of length len. If it matches pind, return
|
|
|
|
|
- * the first quantity, otherwise return the second.
|
|
|
|
|
- */
|
|
|
|
|
- if(pi1 == pi2){
|
|
|
|
|
-
|
|
|
|
|
- fret = (double) (nkp - 1.0)/nkp;
|
|
|
|
|
- if(verbose < 0){
|
|
|
|
|
- printf(" f(%d,%d) = %10.6f\n",pi1,pi2,fret);
|
|
|
|
|
- }
|
|
|
|
|
- return(fret);
|
|
|
|
|
-
|
|
|
|
|
- } else {
|
|
|
|
|
-
|
|
|
|
|
- fret = (double) (-1.0/nkp);
|
|
|
|
|
- if(verbose < 0){
|
|
|
|
|
- printf(" f(%d,%d) = %10.6f\n",pi1,pi2,fret);
|
|
|
|
|
- }
|
|
|
|
|
- return(fret);
|
|
|
|
|
-
|
|
|
|
|
- }
|
|
|
|
|
-
|
|
|
|
|
-
|
|
|
|
|
-}
|
|
|
|
|
-
|
|
|
|
|
-
|
|
|
|
|
-
|
|
|
|
|
-
|
Fabrice Fontaine
Thomas Petazzoni