P-value correction: Difference between revisions
Content deleted Content added
→{{header|Perl}}: removed unused variable and unnecessary checks |
→Version 1: new version does not show compiler warnings about ignored qualifiers or possible conversion errors. Eliminated reindexing step with Hommel method, so code is shorter and should be slightly faster with that method |
||
Line 58: | Line 58: | ||
#include <stdbool.h>//bool data type |
#include <stdbool.h>//bool data type |
||
#include <strings.h>//strcasecmp |
#include <strings.h>//strcasecmp |
||
#include <assert.h>//assert, necessary for random integer selection |
|||
unsigned int * |
unsigned int * seq_len(const unsigned int START, const unsigned int END) { |
||
//named after R function of same name, but simpler function |
//named after R function of same name, but simpler function |
||
unsigned start = (unsigned)START; |
|||
unsigned end = (unsigned)END; |
|||
if (START == END) { |
if (START == END) { |
||
unsigned int *restrict sequence = malloc( (end+1) * sizeof(unsigned int)); |
unsigned int *restrict sequence = malloc( (end+1) * sizeof(unsigned int)); |
||
Line 70: | Line 71: | ||
exit(EXIT_FAILURE); |
exit(EXIT_FAILURE); |
||
} |
} |
||
for ( |
for (unsigned i = 0; i < end; i++) { |
||
sequence[i] = i+1; |
sequence[i] = i+1; |
||
} |
} |
||
Line 76: | Line 77: | ||
} |
} |
||
if (START > END) { |
if (START > END) { |
||
end = START; |
end = (unsigned)START; |
||
start = END; |
start = (unsigned)END; |
||
} |
} |
||
const |
const unsigned LENGTH = end - start ; |
||
unsigned int *restrict sequence = malloc( (1+LENGTH) * sizeof(unsigned int)); |
unsigned int *restrict sequence = malloc( (1+LENGTH) * sizeof(unsigned int)); |
||
if (sequence == NULL) { |
if (sequence == NULL) { |
||
Line 87: | Line 88: | ||
} |
} |
||
if (START < END) { |
if (START < END) { |
||
for ( |
for (unsigned index = 0; index <= LENGTH; index++) { |
||
sequence[index] = start + index; |
sequence[index] = start + index; |
||
} |
} |
||
} else { |
} else { |
||
for ( |
for (unsigned index = 0; index <= LENGTH; index++) { |
||
sequence[index] = end - index; |
sequence[index] = end - index; |
||
} |
} |
||
Line 124: | Line 125: | ||
} |
} |
||
unsigned int * |
unsigned int * order (const double *restrict ARRAY, const unsigned int SIZE, const bool DECREASING) { |
||
//this has the same name as the same R function |
//this has the same name as the same R function |
||
unsigned int *restrict idx = malloc(SIZE * sizeof(unsigned int)); |
unsigned int *restrict idx = malloc(SIZE * sizeof(unsigned int)); |
||
Line 151: | Line 152: | ||
} |
} |
||
double * |
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 takes the same name of the R function which it copies |
||
//this requires a free() afterward where it is used |
//this requires a free() afterward where it is used |
||
Line 175: | Line 176: | ||
} |
} |
||
double * |
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 takes the same name of the R function which it copies |
||
//this requires a free() afterward where it is used |
//this requires a free() afterward where it is used |
||
Line 190: | Line 191: | ||
} |
} |
||
double cumulative_max = ARRAY[0]; |
double cumulative_max = ARRAY[0]; |
||
for ( |
for (unsigned int i = 0; i < NO_OF_ARRAY_ELEMENTS; i++) { |
||
if (ARRAY[i] > cumulative_max) { |
if (ARRAY[i] > cumulative_max) { |
||
cumulative_max = ARRAY[i]; |
cumulative_max = ARRAY[i]; |
||
Line 199: | Line 200: | ||
} |
} |
||
double * |
double * pminx(const double *restrict ARRAY, const unsigned int NO_OF_ARRAY_ELEMENTS, const double X) { |
||
//named after the R function pmin |
//named after the R function pmin |
||
if (NO_OF_ARRAY_ELEMENTS < 1) { |
if (NO_OF_ARRAY_ELEMENTS < 1) { |
||
Line 224: | Line 225: | ||
void double_say (const double *restrict ARRAY, const size_t NO_OF_ARRAY_ELEMENTS) { |
void double_say (const double *restrict ARRAY, const size_t NO_OF_ARRAY_ELEMENTS) { |
||
printf("[1] %e", ARRAY[0]); |
printf("[1] %e", ARRAY[0]); |
||
for ( |
for (unsigned int i = 1; i < NO_OF_ARRAY_ELEMENTS; i++) { |
||
printf(" %.10f", ARRAY[i]); |
printf(" %.10f", ARRAY[i]); |
||
if (((i+1) % 5) == 0) { |
if (((i+1) % 5) == 0) { |
||
printf("\n[% |
printf("\n[%u]", i+1); |
||
} |
} |
||
} |
} |
||
Line 242: | Line 243: | ||
}*/ |
}*/ |
||
double * |
double * uint2double (const unsigned int *restrict ARRAY, const unsigned int NO_OF_ARRAY_ELEMENTS) { |
||
double *restrict doubleArray = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS); |
double *restrict doubleArray = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS); |
||
if (doubleArray == NULL) { |
if (doubleArray == NULL) { |
||
Line 249: | Line 250: | ||
exit(EXIT_FAILURE); |
exit(EXIT_FAILURE); |
||
} |
} |
||
for ( |
for (unsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) { |
||
doubleArray[index] = (double)ARRAY[index]; |
doubleArray[index] = (double)ARRAY[index]; |
||
} |
} |
||
Line 263: | Line 264: | ||
} |
} |
||
double * |
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 |
//this function is a translation of R's p.adjust "BH" method |
||
// i is always i[index] = NO_OF_ARRAY_ELEMENTS - index - 1 |
// i is always i[index] = NO_OF_ARRAY_ELEMENTS - index - 1 |
||
Line 272: | Line 273: | ||
} |
} |
||
short int TYPE = -1; |
short int TYPE = -1; |
||
if |
if (STRING == NULL) { |
||
TYPE = 0; |
|||
} else if (strcasecmp(STRING, "BH") == 0) { |
|||
TYPE = 0; |
TYPE = 0; |
||
} else if (strcasecmp(STRING, "fdr") == 0) { |
} else if (strcasecmp(STRING, "fdr") == 0) { |
||
Line 300: | Line 303: | ||
exit(EXIT_FAILURE); |
exit(EXIT_FAILURE); |
||
} |
} |
||
for ( |
for (unsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) { |
||
const double BONFERRONI = PVALUES[index] * NO_OF_ARRAY_ELEMENTS; |
const double BONFERRONI = PVALUES[index] * NO_OF_ARRAY_ELEMENTS; |
||
if (BONFERRONI >= 1.0) { |
if (BONFERRONI >= 1.0) { |
||
Line 322: | Line 325: | ||
double *restrict o2double = uint2double(o, NO_OF_ARRAY_ELEMENTS); |
double *restrict o2double = uint2double(o, NO_OF_ARRAY_ELEMENTS); |
||
double *restrict cummax_input = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS); |
double *restrict cummax_input = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS); |
||
for ( |
for (unsigned index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) { |
||
cummax_input[index] = (NO_OF_ARRAY_ELEMENTS - index ) * PVALUES[o[index]]; |
cummax_input[index] = (NO_OF_ARRAY_ELEMENTS - index ) * (double)PVALUES[o[index]]; |
||
// printf("cummax_input[%zu] = %e\n", index, cummax_input[index]); |
// printf("cummax_input[%zu] = %e\n", index, cummax_input[index]); |
||
} |
} |
||
Line 336: | Line 339: | ||
free(cummax_output); cummax_output = NULL; |
free(cummax_output); cummax_output = NULL; |
||
double *restrict qvalues = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS); |
double *restrict qvalues = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS); |
||
for ( |
for (unsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) { |
||
qvalues[index] = pmin[ro[index]]; |
qvalues[index] = pmin[ro[index]]; |
||
} |
} |
||
Line 355: | Line 358: | ||
exit(EXIT_FAILURE); |
exit(EXIT_FAILURE); |
||
} |
} |
||
for ( |
for (unsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) { |
||
p[index] = PVALUES[o[index]]; |
p[index] = PVALUES[o[index]]; |
||
} |
} |
||
Line 378: | Line 381: | ||
} |
} |
||
double min = (double)NO_OF_ARRAY_ELEMENTS * p[0]; |
double min = (double)NO_OF_ARRAY_ELEMENTS * p[0]; |
||
for ( |
for (unsigned index = 1; index < NO_OF_ARRAY_ELEMENTS; index++) { |
||
const double TEMP = (double)NO_OF_ARRAY_ELEMENTS * p[index] / (1+index); |
const double TEMP = (double)NO_OF_ARRAY_ELEMENTS * p[index] / (double)(1+index); |
||
if (TEMP < min) { |
if (TEMP < min) { |
||
min = TEMP; |
min = TEMP; |
||
} |
} |
||
} |
} |
||
for ( |
for (unsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) { |
||
pa[index] = min; |
pa[index] = min; |
||
q[index] = min; |
q[index] = min; |
||
Line 399: | Line 402: | ||
} |
} |
||
*/ |
*/ |
||
for ( |
for (unsigned j = (NO_OF_ARRAY_ELEMENTS-1); j >= 2; j--) { |
||
// printf("j = %zu\n", j); |
// printf("j = %zu\n", j); |
||
unsigned int *restrict ij = seq_len( |
unsigned int *restrict ij = seq_len(0,NO_OF_ARRAY_ELEMENTS - j); |
||
for (size_t i = 0; i < NO_OF_ARRAY_ELEMENTS - j + 1; i++) { |
|||
ij[i]--;//R's indices are 1-based, C's are 0-based |
|||
} |
|||
const size_t I2_LENGTH = j - 1; |
const size_t I2_LENGTH = j - 1; |
||
unsigned int *restrict i2 = malloc(I2_LENGTH * sizeof(unsigned int)); |
unsigned int *restrict i2 = malloc(I2_LENGTH * sizeof(unsigned int)); |
||
for ( |
for (unsigned i = 0; i < I2_LENGTH; i++) { |
||
i2[i] = NO_OF_ARRAY_ELEMENTS-j+2+i-1; |
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 |
//R's indices are 1-based, C's are 0-based, I added the -1 |
||
} |
} |
||
double q1 = j * p[i2[0]] / 2.0; |
double q1 = (double)j * p[i2[0]] / 2.0; |
||
for ( |
for (unsigned int i = 1; i < I2_LENGTH; i++) {//loop through 2:j |
||
const double TEMP_Q1 = (double)j * p[i2[i]] / (2 + i); |
const double TEMP_Q1 = (double)j * p[i2[i]] / (double)(2 + i); |
||
if (TEMP_Q1 < q1) { |
if (TEMP_Q1 < q1) { |
||
q1 = TEMP_Q1; |
q1 = TEMP_Q1; |
||
Line 420: | Line 420: | ||
} |
} |
||
for ( |
for (unsigned int i = 0; i < (NO_OF_ARRAY_ELEMENTS - j + 1); i++) {//q[ij] <- pmin(j * p[ij], q1) |
||
q[ij[i]] = min2( j*p[ij[i]], q1); |
q[ij[i]] = min2( (double)j*p[ij[i]], q1); |
||
} |
} |
||
free(ij); ij = NULL; |
free(ij); ij = NULL; |
||
for ( |
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 |
q[i2[i]] = q[NO_OF_ARRAY_ELEMENTS - j];//subtract 1 because of starting index difference |
||
} |
} |
||
free(i2); i2 = NULL; |
free(i2); i2 = NULL; |
||
for ( |
for (unsigned int i = 0; i < NO_OF_ARRAY_ELEMENTS; i++) {//pa <- pmax(pa, q) |
||
if (pa[i] < q[i]) { |
if (pa[i] < q[i]) { |
||
pa[i] = q[i]; |
pa[i] = q[i]; |
||
Line 439: | Line 439: | ||
}//end j loop |
}//end j loop |
||
free(p); p = NULL; |
free(p); p = NULL; |
||
for ( |
for (unsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) { |
||
q[index] = pa[ro[index]];//Hommel q-values |
q[index] = pa[ro[index]];//Hommel q-values |
||
} |
} |
||
Line 472: | Line 472: | ||
double *restrict cummin_input = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS); |
double *restrict cummin_input = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS); |
||
if (TYPE == 0) {//BH method |
if (TYPE == 0) {//BH method |
||
for ( |
for (unsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) { |
||
const double NI = (double)NO_OF_ARRAY_ELEMENTS / (NO_OF_ARRAY_ELEMENTS - index);// n/i simplified |
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] |
cummin_input[index] = NI * PVALUES[o[index]];//PVALUES[o[index]] is p[o] |
||
} |
} |
||
} else if (TYPE == 1) {//BY method |
} else if (TYPE == 1) {//BY method |
||
double q = 1.0; |
double q = 1.0; |
||
for ( |
for (unsigned int index = 2; index < (1+NO_OF_ARRAY_ELEMENTS); index++) { |
||
q += |
q += 1.0/(double)index; |
||
} |
} |
||
for ( |
for (unsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) { |
||
const double NI = (double)NO_OF_ARRAY_ELEMENTS / (NO_OF_ARRAY_ELEMENTS - index);// n/i simplified |
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] |
cummin_input[index] = q * NI * PVALUES[o[index]];//PVALUES[o[index]] is p[o] |
||
} |
} |
||
} else if (TYPE == 3) {//Hochberg method |
} else if (TYPE == 3) {//Hochberg method |
||
for ( |
for (unsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) { |
||
// pmin(1, cummin((n - i + 1L) * p[o]))[ro] |
// pmin(1, cummin((n - i + 1L) * p[o]))[ro] |
||
cummin_input[index] = (index + 1) * PVALUES[o[index]]; |
cummin_input[index] = (double)(index + 1) * PVALUES[o[index]]; |
||
} |
} |
||
} |
} |
||
Line 506: | Line 506: | ||
return q_array; |
return q_array; |
||
} |
} |
||
int main(void) { |
int main(void) { |