P-value correction: Difference between revisions

← Older editNewer edit →
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 *restrict seq_len(const size_t START, const size_t END) {
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
size_t start = START;
unsigned start = (unsigned)START;
size_t end = END;
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 (size_t i = 0; i < end; i++) {
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 size_t LENGTH = end - start ;
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 (size_t index = 0; index <= LENGTH; index++) {
for (unsigned index = 0; index <= LENGTH; index++) {
sequence[index] = start + index;
sequence[index] = start + index;
}
}
} else {
} else {
for (size_t index = 0; index <= LENGTH; index++) {
for (unsigned index = 0; index <= LENGTH; index++) {
sequence[index] = end - index;
sequence[index] = end - index;
}
}
Line 124: Line 125:
}
}


unsigned int *restrict order (const double *restrict ARRAY, const unsigned int SIZE, const bool DECREASING) {
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 *restrict cummin(const double *restrict ARRAY, const unsigned int NO_OF_ARRAY_ELEMENTS) {
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 *restrict cummax(const double *restrict ARRAY, const unsigned int NO_OF_ARRAY_ELEMENTS) {
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 (size_t i = 0; i < NO_OF_ARRAY_ELEMENTS; i++) {
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 *restrict pminx(const double *restrict ARRAY, const size_t NO_OF_ARRAY_ELEMENTS, const double X) {
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 (size_t i = 1; i < NO_OF_ARRAY_ELEMENTS; i++) {
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[%zu]", i+1);
printf("\n[%u]", i+1);
}
}
}
}
Line 242: Line 243:
}*/
}*/


double *restrict uint2double (const unsigned int *restrict ARRAY, const size_t NO_OF_ARRAY_ELEMENTS) {
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 (size_t index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
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 *restrict p_adjust (const double *restrict PVALUES, const size_t NO_OF_ARRAY_ELEMENTS, const char *restrict STRING) {
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 (strcasecmp(STRING, "BH") == 0) {
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 (size_t index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
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 (size_t index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
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 (size_t index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
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 (size_t index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
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 (size_t index = 1; index < NO_OF_ARRAY_ELEMENTS; index++) {
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 (size_t index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
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 (size_t j = (NO_OF_ARRAY_ELEMENTS-1); j >= 2; j--) {
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(1,NO_OF_ARRAY_ELEMENTS - j + 1);
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 (size_t i = 0; i < I2_LENGTH; i++) {
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 (size_t i = 1; i < I2_LENGTH; i++) {//loop through 2:j
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 (size_t i = 0; i < (NO_OF_ARRAY_ELEMENTS - j + 1); i++) {//q[ij] <- pmin(j * p[ij], 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( j*p[ij[i]], q1);
q[ij[i]] = min2( (double)j*p[ij[i]], q1);
}
}
free(ij); ij = NULL;
free(ij); ij = NULL;


for (size_t i = 0; i < I2_LENGTH; i++) {//q[i2] <- q[n - j + 1]
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 (size_t i = 0; i < NO_OF_ARRAY_ELEMENTS; i++) {//pa <- pmax(pa, q)
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 (size_t index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
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 (size_t index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
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 (size_t index = 2; index < (1+NO_OF_ARRAY_ELEMENTS); index++) {
for (unsigned int index = 2; index < (1+NO_OF_ARRAY_ELEMENTS); index++) {
q += (double) 1.0/index;
q += 1.0/(double)index;
}
}
for (size_t index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
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 (size_t index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
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) {
Retrieved from "https://rosettacode.org/wiki/P-value_correction"