Rosetta Code

P-value correction: Difference between revisions

Page Discussion

P-value correction (view source)

Revision as of 15:25, 1 November 2019

177 bytes added ,  4 years ago
→‎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
(→‎{{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:
#include <stdbool.h>//bool data type
#include <strings.h>//strcasecmp
#include <assert.h>//assert, necessary for random integer selection
 
unsigned int *restrict seq_len(const size_tunsigned int START, const size_tunsigned int END) {
//named after R function of same name, but simpler function
size_tunsigned start = (unsigned)START;
size_tunsigned end = (unsigned)END;
if (START == END) {
unsigned int *restrict sequence = malloc( (end+1) * sizeof(unsigned int));
Line 70 ⟶ 71:
exit(EXIT_FAILURE);
}
for (size_tunsigned i = 0; i < end; i++) {
sequence[i] = i+1;
}
Line 76 ⟶ 77:
}
if (START > END) {
end = (unsigned)START;
start = (unsigned)END;
}
const size_tunsigned LENGTH = end - start ;
unsigned int *restrict sequence = malloc( (1+LENGTH) * sizeof(unsigned int));
if (sequence == NULL) {
Line 87 ⟶ 88:
}
if (START < END) {
for (size_tunsigned index = 0; index <= LENGTH; index++) {
sequence[index] = start + index;
}
} else {
for (size_tunsigned index = 0; index <= LENGTH; index++) {
sequence[index] = end - index;
}
Line 124 ⟶ 125:
}
 
unsigned int *restrict 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));
Line 151 ⟶ 152:
}
 
double *restrict 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
Line 175 ⟶ 176:
}
 
double *restrict 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
Line 190 ⟶ 191:
}
double cumulative_max = ARRAY[0];
for (size_tunsigned int i = 0; i < NO_OF_ARRAY_ELEMENTS; i++) {
if (ARRAY[i] > cumulative_max) {
cumulative_max = ARRAY[i];
Line 199 ⟶ 200:
}
 
double *restrict pminx(const double *restrict ARRAY, const size_tunsigned int NO_OF_ARRAY_ELEMENTS, const double X) {
//named after the R function pmin
if (NO_OF_ARRAY_ELEMENTS < 1) {
Line 224 ⟶ 225:
void double_say (const double *restrict ARRAY, const size_t NO_OF_ARRAY_ELEMENTS) {
printf("[1] %e", ARRAY[0]);
for (size_tunsigned int i = 1; i < NO_OF_ARRAY_ELEMENTS; i++) {
printf(" %.10f", ARRAY[i]);
if (((i+1) % 5) == 0) {
printf("\n[%zuu]", i+1);
}
}
Line 242 ⟶ 243:
}*/
 
double *restrict uint2double (const unsigned int *restrict ARRAY, const size_tunsigned int NO_OF_ARRAY_ELEMENTS) {
double *restrict doubleArray = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS);
if (doubleArray == NULL) {
Line 249 ⟶ 250:
exit(EXIT_FAILURE);
}
for (size_tunsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
doubleArray[index] = (double)ARRAY[index];
}
Line 263 ⟶ 264:
}
 
double *restrict p_adjust (const double *restrict PVALUES, const size_tunsigned 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
Line 272 ⟶ 273:
}
short int TYPE = -1;
if (strcasecmp(STRING, "BH") == 0NULL) {
TYPE = 0;
} else if (strcasecmp(STRING, "BH") == 0) {
TYPE = 0;
} else if (strcasecmp(STRING, "fdr") == 0) {
Line 300 ⟶ 303:
exit(EXIT_FAILURE);
}
for (size_tunsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
const double BONFERRONI = PVALUES[index] * NO_OF_ARRAY_ELEMENTS;
if (BONFERRONI >= 1.0) {
Line 322 ⟶ 325:
double *restrict o2double = uint2double(o, NO_OF_ARRAY_ELEMENTS);
double *restrict cummax_input = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS);
for (size_tunsigned 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]);
}
Line 336 ⟶ 339:
free(cummax_output); cummax_output = NULL;
double *restrict qvalues = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS);
for (size_tunsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
qvalues[index] = pmin[ro[index]];
}
Line 355 ⟶ 358:
exit(EXIT_FAILURE);
}
for (size_tunsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
p[index] = PVALUES[o[index]];
}
Line 378 ⟶ 381:
}
double min = (double)NO_OF_ARRAY_ELEMENTS * p[0];
for (size_tunsigned 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 (size_tunsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
pa[index] = min;
q[index] = min;
Line 399 ⟶ 402:
}
*/
for (size_tunsigned j = (NO_OF_ARRAY_ELEMENTS-1); j >= 2; j--) {
// printf("j = %zu\n", j);
unsigned int *restrict ij = seq_len(10,NO_OF_ARRAY_ELEMENTS - j + 1);
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;
unsigned int *restrict i2 = malloc(I2_LENGTH * sizeof(unsigned int));
for (size_tunsigned 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 (size_tunsigned 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;
Line 420:
}
 
for (size_tunsigned 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 (size_tunsigned 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 (size_tunsigned int i = 0; i < NO_OF_ARRAY_ELEMENTS; i++) {//pa <- pmax(pa, q)
if (pa[i] < q[i]) {
pa[i] = q[i];
Line 439:
}//end j loop
free(p); p = NULL;
for (size_tunsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
q[index] = pa[ro[index]];//Hommel q-values
}
Line 472:
double *restrict cummin_input = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS);
if (TYPE == 0) {//BH method
for (size_tunsigned 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 (size_tunsigned int index = 2; index < (1+NO_OF_ARRAY_ELEMENTS); index++) {
q += (double) 1.0/(double)index;
}
for (size_tunsigned 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 (size_tunsigned 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]];
}
}
Line 506:
return q_array;
}
 
 
int main(void) {
Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/278888"