548 lines
21 KiB
C
548 lines
21 KiB
C
#include <stdio.h>//printf
|
|
#include <stdlib.h>//qsort
|
|
#include <math.h>//fabs
|
|
#include <stdbool.h>//bool data type
|
|
#include <strings.h>//strcasecmp
|
|
#include <assert.h>//assert, necessary for random integer selection
|
|
|
|
unsigned int * seq_len(const unsigned int START, const unsigned int END) {
|
|
//named after R function of same name, but simpler function
|
|
unsigned start = (unsigned)START;
|
|
unsigned end = (unsigned)END;
|
|
if (START == END) {
|
|
unsigned int *restrict sequence = malloc( (end+1) * sizeof(unsigned int));
|
|
if (sequence == NULL) {
|
|
printf("malloc failed at %s line %u\n", __FILE__, __LINE__);
|
|
perror("");
|
|
exit(EXIT_FAILURE);
|
|
}
|
|
for (unsigned i = 0; i < end; i++) {
|
|
sequence[i] = i+1;
|
|
}
|
|
return sequence;
|
|
}
|
|
if (START > END) {
|
|
end = (unsigned)START;
|
|
start = (unsigned)END;
|
|
}
|
|
const unsigned LENGTH = end - start ;
|
|
unsigned int *restrict sequence = malloc( (1+LENGTH) * sizeof(unsigned int));
|
|
if (sequence == NULL) {
|
|
printf("malloc failed at %s line %u\n", __FILE__, __LINE__);
|
|
perror("");
|
|
exit(EXIT_FAILURE);
|
|
}
|
|
if (START < END) {
|
|
for (unsigned index = 0; index <= LENGTH; index++) {
|
|
sequence[index] = start + index;
|
|
}
|
|
} else {
|
|
for (unsigned index = 0; index <= LENGTH; index++) {
|
|
sequence[index] = end - index;
|
|
}
|
|
}
|
|
return sequence;
|
|
}
|
|
|
|
//modified from https://phoxis.org/2012/07/12/get-sorted-index-orderting-of-an-array/
|
|
|
|
double *restrict base_arr = NULL;
|
|
|
|
static int compar_increase (const void *restrict a, const void *restrict b) {
|
|
int aa = *((int *restrict ) a), bb = *((int *restrict) b);
|
|
if (base_arr[aa] < base_arr[bb]) {
|
|
return 1;
|
|
} else if (base_arr[aa] == base_arr[bb]) {
|
|
return 0;
|
|
} else {
|
|
return -1;
|
|
}
|
|
}
|
|
|
|
static int compar_decrease (const void *restrict a, const void *restrict b) {
|
|
int aa = *((int *restrict ) a), bb = *((int *restrict) b);
|
|
if (base_arr[aa] < base_arr[bb]) {
|
|
return -1;
|
|
} else if (base_arr[aa] == base_arr[bb]) {
|
|
return 0;
|
|
} else {
|
|
return 1;
|
|
}
|
|
}
|
|
|
|
unsigned int * order (const double *restrict ARRAY, const unsigned int SIZE, const bool DECREASING) {
|
|
//this has the same name as the same R function
|
|
unsigned int *restrict idx = malloc(SIZE * sizeof(unsigned int));
|
|
if (idx == NULL) {
|
|
printf("failed to malloc at %s line %u.\n", __FILE__, __LINE__);
|
|
perror("");
|
|
exit(EXIT_FAILURE);
|
|
}
|
|
base_arr = malloc(sizeof(double) * SIZE);
|
|
if (base_arr == NULL) {
|
|
printf("failed to malloc at %s line %u.\n", __FILE__, __LINE__);
|
|
perror("");
|
|
exit(EXIT_FAILURE);
|
|
}
|
|
for (unsigned int i = 0; i < SIZE; i++) {
|
|
base_arr[i] = ARRAY[i];
|
|
idx[i] = i;
|
|
}
|
|
if (DECREASING == false) {
|
|
qsort(idx, SIZE, sizeof(unsigned int), compar_decrease);
|
|
} else if (DECREASING == true) {
|
|
qsort(idx, SIZE, sizeof(unsigned int), compar_increase);
|
|
}
|
|
free(base_arr); base_arr = NULL;
|
|
return idx;
|
|
}
|
|
|
|
double * cummin(const double *restrict ARRAY, const unsigned int NO_OF_ARRAY_ELEMENTS) {
|
|
//this takes the same name of the R function which it copies
|
|
//this requires a free() afterward where it is used
|
|
if (NO_OF_ARRAY_ELEMENTS < 1) {
|
|
puts("cummin function requires at least one element.\n");
|
|
printf("Failed at %s line %u\n", __FILE__, __LINE__);
|
|
exit(EXIT_FAILURE);
|
|
}
|
|
double *restrict output = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS);
|
|
if (output == NULL) {
|
|
printf("failed to malloc at %s line %u.\n", __FILE__, __LINE__);
|
|
perror("");
|
|
exit(EXIT_FAILURE);
|
|
}
|
|
double cumulative_min = ARRAY[0];
|
|
for (unsigned int i = 0; i < NO_OF_ARRAY_ELEMENTS; i++) {
|
|
if (ARRAY[i] < cumulative_min) {
|
|
cumulative_min = ARRAY[i];
|
|
}
|
|
output[i] = cumulative_min;
|
|
}
|
|
return output;
|
|
}
|
|
|
|
double * cummax(const double *restrict ARRAY, const unsigned int NO_OF_ARRAY_ELEMENTS) {
|
|
//this takes the same name of the R function which it copies
|
|
//this requires a free() afterward where it is used
|
|
if (NO_OF_ARRAY_ELEMENTS < 1) {
|
|
puts("function requires at least one element.\n");
|
|
printf("Failed at %s line %u\n", __FILE__, __LINE__);
|
|
exit(EXIT_FAILURE);
|
|
}
|
|
double *restrict output = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS);
|
|
if (output == NULL) {
|
|
printf("failed to malloc at %s line %u.\n", __FILE__, __LINE__);
|
|
perror("");
|
|
exit(EXIT_FAILURE);
|
|
}
|
|
double cumulative_max = ARRAY[0];
|
|
for (unsigned int i = 0; i < NO_OF_ARRAY_ELEMENTS; i++) {
|
|
if (ARRAY[i] > cumulative_max) {
|
|
cumulative_max = ARRAY[i];
|
|
}
|
|
output[i] = cumulative_max;
|
|
}
|
|
return output;
|
|
}
|
|
|
|
double * pminx(const double *restrict ARRAY, const unsigned int NO_OF_ARRAY_ELEMENTS, const double X) {
|
|
//named after the R function pmin
|
|
if (NO_OF_ARRAY_ELEMENTS < 1) {
|
|
puts("pmin requires at least one element.\n");
|
|
printf("Failed at %s line %u\n", __FILE__, __LINE__);
|
|
exit(EXIT_FAILURE);
|
|
}
|
|
double *restrict pmin_array = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS);
|
|
if (pmin_array == NULL) {
|
|
printf("failed to malloc at %s line %u.\n", __FILE__, __LINE__);
|
|
perror("");
|
|
exit(EXIT_FAILURE);
|
|
}
|
|
for (unsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
|
|
if (ARRAY[index] < X) {
|
|
pmin_array[index] = ARRAY[index];
|
|
} else {
|
|
pmin_array[index] = X;
|
|
}
|
|
}
|
|
return pmin_array;
|
|
}
|
|
|
|
void double_say (const double *restrict ARRAY, const size_t NO_OF_ARRAY_ELEMENTS) {
|
|
printf("[1] %e", ARRAY[0]);
|
|
for (unsigned int i = 1; i < NO_OF_ARRAY_ELEMENTS; i++) {
|
|
printf(" %.10f", ARRAY[i]);
|
|
if (((i+1) % 5) == 0) {
|
|
printf("\n[%u]", i+1);
|
|
}
|
|
}
|
|
puts("\n");
|
|
}
|
|
|
|
/*void uint_say (const unsigned int *restrict ARRAY, const size_t NO_OF_ARRAY_ELEMENTS) {
|
|
//for debugging
|
|
printf("%u", ARRAY[0]);
|
|
for (size_t i = 1; i < NO_OF_ARRAY_ELEMENTS; i++) {
|
|
printf(",%u", ARRAY[i]);
|
|
}
|
|
puts("\n");
|
|
}*/
|
|
|
|
double * uint2double (const unsigned int *restrict ARRAY, const unsigned int NO_OF_ARRAY_ELEMENTS) {
|
|
double *restrict doubleArray = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS);
|
|
if (doubleArray == NULL) {
|
|
printf("Failure to malloc at %s line %u.\n", __FILE__, __LINE__);
|
|
perror("");
|
|
exit(EXIT_FAILURE);
|
|
}
|
|
for (unsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
|
|
doubleArray[index] = (double)ARRAY[index];
|
|
}
|
|
return doubleArray;
|
|
}
|
|
|
|
double min2 (const double N1, const double N2) {
|
|
if (N1 < N2) {
|
|
return N1;
|
|
} else {
|
|
return N2;
|
|
}
|
|
}
|
|
|
|
double * p_adjust (const double *restrict PVALUES, const unsigned int NO_OF_ARRAY_ELEMENTS, const char *restrict STRING) {
|
|
//this function is a translation of R's p.adjust "BH" method
|
|
// i is always i[index] = NO_OF_ARRAY_ELEMENTS - index - 1
|
|
if (NO_OF_ARRAY_ELEMENTS < 1) {
|
|
puts("p_adjust requires at least one element.\n");
|
|
printf("Failed at %s line %u\n", __FILE__, __LINE__);
|
|
exit(EXIT_FAILURE);
|
|
}
|
|
short int TYPE = -1;
|
|
if (STRING == NULL) {
|
|
TYPE = 0;
|
|
} else if (strcasecmp(STRING, "BH") == 0) {
|
|
TYPE = 0;
|
|
} else if (strcasecmp(STRING, "fdr") == 0) {
|
|
TYPE = 0;
|
|
} else if (strcasecmp(STRING, "by") == 0) {
|
|
TYPE = 1;
|
|
} else if (strcasecmp(STRING, "Bonferroni") == 0) {
|
|
TYPE = 2;
|
|
} else if (strcasecmp(STRING, "hochberg") == 0) {
|
|
TYPE = 3;
|
|
} else if (strcasecmp(STRING, "holm") == 0) {
|
|
TYPE = 4;
|
|
} else if (strcasecmp(STRING, "hommel") == 0) {
|
|
TYPE = 5;
|
|
} else {
|
|
printf("%s doesn't match any accepted FDR methods.\n", STRING);
|
|
printf("Failed at %s line %u\n", __FILE__, __LINE__);
|
|
exit(EXIT_FAILURE);
|
|
}
|
|
//---------------------------------------------------------------------------
|
|
//---------------------------------------------------------------------------
|
|
if (TYPE == 2) {//Bonferroni method
|
|
double *restrict bonferroni = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS);
|
|
if (bonferroni == NULL) {
|
|
printf("failed to malloc at %s line %u.\n", __FILE__, __LINE__);
|
|
perror("");
|
|
exit(EXIT_FAILURE);
|
|
}
|
|
for (unsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
|
|
const double BONFERRONI = PVALUES[index] * NO_OF_ARRAY_ELEMENTS;
|
|
if (BONFERRONI >= 1.0) {
|
|
bonferroni[index] = 1.0;
|
|
} else if ((0.0 <= BONFERRONI) && (BONFERRONI < 1.0)) {
|
|
bonferroni[index] = BONFERRONI;
|
|
} else {
|
|
printf("%g is outside of the interval I planned.\n", BONFERRONI);
|
|
printf("Failure at %s line %u\n", __FILE__, __LINE__);
|
|
exit(EXIT_FAILURE);
|
|
}
|
|
}
|
|
return bonferroni;
|
|
//---------------------------------------------------------------------------
|
|
//---------------------------------------------------------------------------
|
|
} else if (TYPE == 4) {//Holm method
|
|
/*these values are computed separately from BH, BY, and Hochberg because they are
|
|
computed differently*/
|
|
unsigned int *restrict o = order(PVALUES, NO_OF_ARRAY_ELEMENTS, false);
|
|
//sorted in reverse of methods 0-3
|
|
double *restrict o2double = uint2double(o, NO_OF_ARRAY_ELEMENTS);
|
|
double *restrict cummax_input = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS);
|
|
for (unsigned index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
|
|
cummax_input[index] = (NO_OF_ARRAY_ELEMENTS - index ) * (double)PVALUES[o[index]];
|
|
// printf("cummax_input[%zu] = %e\n", index, cummax_input[index]);
|
|
}
|
|
free(o); o = NULL;
|
|
unsigned int *restrict ro = order(o2double, NO_OF_ARRAY_ELEMENTS, false);
|
|
free(o2double); o2double = NULL;
|
|
|
|
double *restrict cummax_output = cummax(cummax_input, NO_OF_ARRAY_ELEMENTS);
|
|
free(cummax_input); cummax_input = NULL;
|
|
|
|
double *restrict pmin = pminx(cummax_output, NO_OF_ARRAY_ELEMENTS, 1);
|
|
free(cummax_output); cummax_output = NULL;
|
|
double *restrict qvalues = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS);
|
|
for (unsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
|
|
qvalues[index] = pmin[ro[index]];
|
|
}
|
|
free(pmin); pmin = NULL;
|
|
free(ro); ro = NULL;
|
|
return qvalues;
|
|
//---------------------------------------------------------------------------
|
|
//---------------------------------------------------------------------------
|
|
} else if (TYPE == 5) {//Hommel method
|
|
//i <- seq_len(n)
|
|
//o <- order(p)
|
|
unsigned int *restrict o = order(PVALUES, NO_OF_ARRAY_ELEMENTS, false);//false is R's default
|
|
//p <- p[o]
|
|
double *restrict p = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS);
|
|
if (p == NULL) {
|
|
printf("failed to malloc at %s line %u.\n", __FILE__, __LINE__);
|
|
perror("");
|
|
exit(EXIT_FAILURE);
|
|
}
|
|
for (unsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
|
|
p[index] = PVALUES[o[index]];
|
|
}
|
|
//ro <- order(o)
|
|
double *restrict o2double = uint2double(o, NO_OF_ARRAY_ELEMENTS);
|
|
free(o); o = NULL;
|
|
unsigned int *restrict ro = order(o2double, NO_OF_ARRAY_ELEMENTS, false);
|
|
free(o2double); o2double = NULL;
|
|
// puts("ro");
|
|
//q <- pa <- rep.int(min(n * p/i), n)
|
|
double *restrict q = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS);
|
|
if (q == NULL) {
|
|
printf("failed to malloc at %s line %u.\n", __FILE__, __LINE__);
|
|
perror("");
|
|
exit(EXIT_FAILURE);
|
|
}
|
|
double *restrict pa = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS);
|
|
if (pa == NULL) {
|
|
printf("failed to malloc at %s line %u.\n", __FILE__, __LINE__);
|
|
perror("");
|
|
exit(EXIT_FAILURE);
|
|
}
|
|
double min = (double)NO_OF_ARRAY_ELEMENTS * p[0];
|
|
for (unsigned index = 1; index < NO_OF_ARRAY_ELEMENTS; index++) {
|
|
const double TEMP = (double)NO_OF_ARRAY_ELEMENTS * p[index] / (double)(1+index);
|
|
if (TEMP < min) {
|
|
min = TEMP;
|
|
}
|
|
}
|
|
for (unsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
|
|
pa[index] = min;
|
|
q[index] = min;
|
|
}
|
|
// puts("q & pa");
|
|
// double_say(q, NO_OF_ARRAY_ELEMENTS);
|
|
/*for (j in (n - 1):2) {
|
|
ij <- seq_len(n - j + 1)
|
|
i2 <- (n - j + 2):n
|
|
q1 <- min(j * p[i2]/(2:j))
|
|
q[ij] <- pmin(j * p[ij], q1)
|
|
q[i2] <- q[n - j + 1]
|
|
pa <- pmax(pa, q)
|
|
}
|
|
*/
|
|
for (unsigned j = (NO_OF_ARRAY_ELEMENTS-1); j >= 2; j--) {
|
|
// printf("j = %zu\n", j);
|
|
unsigned int *restrict ij = seq_len(0,NO_OF_ARRAY_ELEMENTS - j);
|
|
const size_t I2_LENGTH = j - 1;
|
|
unsigned int *restrict i2 = malloc(I2_LENGTH * sizeof(unsigned int));
|
|
for (unsigned i = 0; i < I2_LENGTH; i++) {
|
|
i2[i] = NO_OF_ARRAY_ELEMENTS-j+2+i-1;
|
|
//R's indices are 1-based, C's are 0-based, I added the -1
|
|
}
|
|
|
|
double q1 = (double)j * p[i2[0]] / 2.0;
|
|
for (unsigned int i = 1; i < I2_LENGTH; i++) {//loop through 2:j
|
|
const double TEMP_Q1 = (double)j * p[i2[i]] / (double)(2 + i);
|
|
if (TEMP_Q1 < q1) {
|
|
q1 = TEMP_Q1;
|
|
}
|
|
}
|
|
|
|
for (unsigned int i = 0; i < (NO_OF_ARRAY_ELEMENTS - j + 1); i++) {//q[ij] <- pmin(j * p[ij], q1)
|
|
q[ij[i]] = min2( (double)j*p[ij[i]], q1);
|
|
}
|
|
free(ij); ij = NULL;
|
|
|
|
for (unsigned int i = 0; i < I2_LENGTH; i++) {//q[i2] <- q[n - j + 1]
|
|
q[i2[i]] = q[NO_OF_ARRAY_ELEMENTS - j];//subtract 1 because of starting index difference
|
|
}
|
|
free(i2); i2 = NULL;
|
|
|
|
for (unsigned int i = 0; i < NO_OF_ARRAY_ELEMENTS; i++) {//pa <- pmax(pa, q)
|
|
if (pa[i] < q[i]) {
|
|
pa[i] = q[i];
|
|
}
|
|
}
|
|
// printf("j = %zu, pa = \n", j);
|
|
// double_say(pa, N);
|
|
}//end j loop
|
|
free(p); p = NULL;
|
|
for (unsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
|
|
q[index] = pa[ro[index]];//Hommel q-values
|
|
}
|
|
//now free memory
|
|
free(ro); ro = NULL;
|
|
free(pa); pa = NULL;
|
|
return q;
|
|
}
|
|
//The methods are similarly computed and thus can be combined for clarity
|
|
unsigned int *restrict o = order(PVALUES, NO_OF_ARRAY_ELEMENTS, true);
|
|
if (o == NULL) {
|
|
printf("failed to malloc at %s line %u.\n", __FILE__, __LINE__);
|
|
perror("");
|
|
exit(EXIT_FAILURE);
|
|
}
|
|
double *restrict o_double = uint2double(o, NO_OF_ARRAY_ELEMENTS);
|
|
for (unsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
|
|
if ((PVALUES[index] < 0) || (PVALUES[index] > 1)) {
|
|
printf("array[%u] = %lf, which is outside the interval [0,1]\n", index, PVALUES[index]);
|
|
printf("died at %s line %u\n", __FILE__, __LINE__);
|
|
exit(EXIT_FAILURE);
|
|
}
|
|
}
|
|
|
|
unsigned int *restrict ro = order(o_double, NO_OF_ARRAY_ELEMENTS, false);
|
|
if (ro == NULL) {
|
|
printf("failed to malloc at %s line %u.\n", __FILE__, __LINE__);
|
|
perror("");
|
|
exit(EXIT_FAILURE);
|
|
}
|
|
free(o_double); o_double = NULL;
|
|
double *restrict cummin_input = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS);
|
|
if (TYPE == 0) {//BH method
|
|
for (unsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
|
|
const double NI = (double)NO_OF_ARRAY_ELEMENTS / (double)(NO_OF_ARRAY_ELEMENTS - index);// n/i simplified
|
|
cummin_input[index] = NI * PVALUES[o[index]];//PVALUES[o[index]] is p[o]
|
|
}
|
|
} else if (TYPE == 1) {//BY method
|
|
double q = 1.0;
|
|
for (unsigned int index = 2; index < (1+NO_OF_ARRAY_ELEMENTS); index++) {
|
|
q += 1.0/(double)index;
|
|
}
|
|
for (unsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
|
|
const double NI = (double)NO_OF_ARRAY_ELEMENTS / (double)(NO_OF_ARRAY_ELEMENTS - index);// n/i simplified
|
|
cummin_input[index] = q * NI * PVALUES[o[index]];//PVALUES[o[index]] is p[o]
|
|
}
|
|
} else if (TYPE == 3) {//Hochberg method
|
|
for (unsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
|
|
// pmin(1, cummin((n - i + 1L) * p[o]))[ro]
|
|
cummin_input[index] = (double)(index + 1) * PVALUES[o[index]];
|
|
}
|
|
}
|
|
free(o); o = NULL;
|
|
double *restrict cummin_array = NULL;
|
|
cummin_array = cummin(cummin_input, NO_OF_ARRAY_ELEMENTS);
|
|
free(cummin_input); cummin_input = NULL;//I don't need this anymore
|
|
double *restrict pmin = pminx(cummin_array, NO_OF_ARRAY_ELEMENTS, 1);
|
|
free(cummin_array); cummin_array = NULL;
|
|
double *restrict q_array = malloc(NO_OF_ARRAY_ELEMENTS*sizeof(double));
|
|
for (unsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
|
|
q_array[index] = pmin[ro[index]];
|
|
}
|
|
|
|
free(ro); ro = NULL;
|
|
free(pmin); pmin = NULL;
|
|
return q_array;
|
|
}
|
|
|
|
|
|
int main(void) {
|
|
const double PVALUES[] = {4.533744e-01, 7.296024e-01, 9.936026e-02, 9.079658e-02, 1.801962e-01,
|
|
8.752257e-01, 2.922222e-01, 9.115421e-01, 4.355806e-01, 5.324867e-01,
|
|
4.926798e-01, 5.802978e-01, 3.485442e-01, 7.883130e-01, 2.729308e-01,
|
|
8.502518e-01, 4.268138e-01, 6.442008e-01, 3.030266e-01, 5.001555e-02,
|
|
3.194810e-01, 7.892933e-01, 9.991834e-01, 1.745691e-01, 9.037516e-01,
|
|
1.198578e-01, 3.966083e-01, 1.403837e-02, 7.328671e-01, 6.793476e-02,
|
|
4.040730e-03, 3.033349e-04, 1.125147e-02, 2.375072e-02, 5.818542e-04,
|
|
3.075482e-04, 8.251272e-03, 1.356534e-03, 1.360696e-02, 3.764588e-04,
|
|
1.801145e-05, 2.504456e-07, 3.310253e-02, 9.427839e-03, 8.791153e-04,
|
|
2.177831e-04, 9.693054e-04, 6.610250e-05, 2.900813e-02, 5.735490e-03};//just the pvalues
|
|
const double CORRECT_ANSWERS[6][50] = {//each first index is type
|
|
{6.126681e-01, 8.521710e-01, 1.987205e-01, 1.891595e-01, 3.217789e-01,
|
|
9.301450e-01, 4.870370e-01, 9.301450e-01, 6.049731e-01, 6.826753e-01,
|
|
6.482629e-01, 7.253722e-01, 5.280973e-01, 8.769926e-01, 4.705703e-01,
|
|
9.241867e-01, 6.049731e-01, 7.856107e-01, 4.887526e-01, 1.136717e-01,
|
|
4.991891e-01, 8.769926e-01, 9.991834e-01, 3.217789e-01, 9.301450e-01,
|
|
2.304958e-01, 5.832475e-01, 3.899547e-02, 8.521710e-01, 1.476843e-01,
|
|
1.683638e-02, 2.562902e-03, 3.516084e-02, 6.250189e-02, 3.636589e-03,
|
|
2.562902e-03, 2.946883e-02, 6.166064e-03, 3.899547e-02, 2.688991e-03,
|
|
4.502862e-04, 1.252228e-05, 7.881555e-02, 3.142613e-02, 4.846527e-03,
|
|
2.562902e-03, 4.846527e-03, 1.101708e-03, 7.252032e-02, 2.205958e-02},//Benjamini-Hochberg
|
|
{1.000000e+00, 1.000000e+00, 8.940844e-01, 8.510676e-01, 1.000000e+00,
|
|
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
|
|
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
|
|
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 5.114323e-01,
|
|
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
|
|
1.000000e+00, 1.000000e+00, 1.754486e-01, 1.000000e+00, 6.644618e-01,
|
|
7.575031e-02, 1.153102e-02, 1.581959e-01, 2.812089e-01, 1.636176e-02,
|
|
1.153102e-02, 1.325863e-01, 2.774239e-02, 1.754486e-01, 1.209832e-02,
|
|
2.025930e-03, 5.634031e-05, 3.546073e-01, 1.413926e-01, 2.180552e-02,
|
|
1.153102e-02, 2.180552e-02, 4.956812e-03, 3.262838e-01, 9.925057e-02},//Benjamini & Yekutieli
|
|
{1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
|
|
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
|
|
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
|
|
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
|
|
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
|
|
1.000000e+00, 1.000000e+00, 7.019185e-01, 1.000000e+00, 1.000000e+00,
|
|
2.020365e-01, 1.516674e-02, 5.625735e-01, 1.000000e+00, 2.909271e-02,
|
|
1.537741e-02, 4.125636e-01, 6.782670e-02, 6.803480e-01, 1.882294e-02,
|
|
9.005725e-04, 1.252228e-05, 1.000000e+00, 4.713920e-01, 4.395577e-02,
|
|
1.088915e-02, 4.846527e-02, 3.305125e-03, 1.000000e+00, 2.867745e-01},//Bonferroni
|
|
{9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01,
|
|
9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01,
|
|
9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01,
|
|
9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01,
|
|
9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01,
|
|
9.991834e-01, 9.991834e-01, 4.632662e-01, 9.991834e-01, 9.991834e-01,
|
|
1.575885e-01, 1.383967e-02, 3.938014e-01, 7.600230e-01, 2.501973e-02,
|
|
1.383967e-02, 3.052971e-01, 5.426136e-02, 4.626366e-01, 1.656419e-02,
|
|
8.825610e-04, 1.252228e-05, 9.930759e-01, 3.394022e-01, 3.692284e-02,
|
|
1.023581e-02, 3.974152e-02, 3.172920e-03, 8.992520e-01, 2.179486e-01},//Hochberg
|
|
{1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
|
|
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
|
|
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
|
|
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
|
|
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
|
|
1.000000e+00, 1.000000e+00, 4.632662e-01, 1.000000e+00, 1.000000e+00,
|
|
1.575885e-01, 1.395341e-02, 3.938014e-01, 7.600230e-01, 2.501973e-02,
|
|
1.395341e-02, 3.052971e-01, 5.426136e-02, 4.626366e-01, 1.656419e-02,
|
|
8.825610e-04, 1.252228e-05, 9.930759e-01, 3.394022e-01, 3.692284e-02,
|
|
1.023581e-02, 3.974152e-02, 3.172920e-03, 8.992520e-01, 2.179486e-01},//Holm
|
|
{ 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.987624e-01, 9.991834e-01,
|
|
9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01,
|
|
9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01,
|
|
9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.595180e-01,
|
|
9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01,
|
|
9.991834e-01, 9.991834e-01, 4.351895e-01, 9.991834e-01, 9.766522e-01,
|
|
1.414256e-01, 1.304340e-02, 3.530937e-01, 6.887709e-01, 2.385602e-02,
|
|
1.322457e-02, 2.722920e-01, 5.426136e-02, 4.218158e-01, 1.581127e-02,
|
|
8.825610e-04, 1.252228e-05, 8.743649e-01, 3.016908e-01, 3.516461e-02,
|
|
9.582456e-03, 3.877222e-02, 3.172920e-03, 8.122276e-01, 1.950067e-01}//Hommel
|
|
};
|
|
//the following loop checks each type with R's answers
|
|
const char *restrict TYPES[] = {"bh", "by", "bonferroni", "hochberg", "holm", "hommel"};
|
|
for (unsigned short int type = 0; type <= 5; type++) {
|
|
double *restrict q = p_adjust(PVALUES, sizeof(PVALUES) / sizeof(*PVALUES), TYPES[type]);
|
|
double error = fabs(q[0] - CORRECT_ANSWERS[type][0]);
|
|
// printf("%e - %e = %g\n", q[0], CORRECT_ANSWERS[type][0], error);
|
|
// puts("p q");
|
|
// printf("%g\t%g\n", pvalues[0], q[0]);
|
|
for (unsigned int i = 1; i < sizeof(PVALUES) / sizeof(*PVALUES); i++) {
|
|
const double this_error = fabs(q[i] - CORRECT_ANSWERS[type][i]);
|
|
// printf("%e - %e = %g\n", q[i], CORRECT_ANSWERS[type][i], error);
|
|
error += this_error;
|
|
}
|
|
double_say(q, sizeof(PVALUES) / sizeof(*PVALUES));
|
|
free(q); q = NULL;
|
|
printf("\ntype %u = '%s' has cumulative error of %g\n", type, TYPES[type], error);
|
|
}
|
|
|
|
return 0;
|
|
}
|