📄 main.cpp
字号:
free(Wj);
fclose(fp);
//#ifdef GEN_TIMING_INFO
// finish = clock();
// fprintf(fp1, "%d\n", finish - start);
// fclose(fp1);
//#endif
// return;
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
O V E R L A P P I N G T E M P L A T E T E S T
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
extern double lgam ( double x );
double Pr(int u, double eta)
{
int l;
double sum, p;
if (u == 0)
p = exp(-eta);
else {
sum = 0.0;
for(l = 1; l <= u; l++)
sum += exp(-eta-u*log(2)+l*log(eta)-lgam(l+1)+lgam(u)-lgam(l)-
lgam(u-l+1));
/*sum += exp(-eta-u*log(2)+l*log(eta)-gammln(l+1)+gammln(u)-gammln(l)-
gammln(u-l+1));*/
p = sum;
}
return p;
}
bool OverlappingTemplateMatchings(int m, int n)
{
int i, k, match, state;
double W_obs, eta, sum, chi2;
double p_value, lambda;
//FILE *output;
BitField *templates;
//char assignment[15];
int M, N, j, K = 5;
unsigned int nu[6] = {0,0,0,0,0,0};
double pi[6] = {0.143783,0.139430,0.137319,0.124314,0.106209,0.348945};
/* double U[6] = {0,1,2,3,4}; */
/*#ifdef GEN_TIMING_INFO
clock_t start, finish;
FILE *fp;
fp = fopen("Overlapping.txt", "a");
start = clock();
#endif*/
M = 1032;
N = (int)floor(n/M);
if ( (templates = (BitField*)calloc(m,sizeof(BitField))) == NULL ) {
printf("\t\t OVERLAPPING TEMPLATE OF ALL ONES TEST\n");
printf("\t\t---------------------------------------------\n");
printf("\t\tTEMPLATE DEFINITION: Insufficient memory,");
printf(" Overlapping Template Matchings test aborted!\n");
state = 0;
}
else
for( i=0; i<m; i++ )
templates[i].b = 1;
lambda = (double)(M-m+1)/pow(2,m);
eta = lambda/2.0;
sum = 0.0;
for( i=0; i<K; i++ ) { /* Compute Probabilities */
pi[i] = Pr(i,eta);
sum += pi[i];
}
pi[K] = 1 - sum;
for( i=0; i<N; i++ ) {
W_obs = 0;
for( j=0; j<M-m+1; j++ ) {
match = 1;
for( k=0; k<m; k++ )
if ( (int)templates[k].b != (int)bit[i*M+j+k] )
match = 0;
if ( match == 1 )
W_obs++;
}
if ( W_obs <= 4 )
nu[(int)W_obs]++;
else
nu[K]++;
}
sum = 0;
chi2 = 0.0; /* Compute Chi Square */
for( i=0; i<K+1; i++ ) {
chi2 += pow((double)nu[i] - (double)N*pi[i],2)/((double)N*pi[i]);
sum += nu[i];
}
p_value = igamc(K/2.,chi2/2.);
/*p_value = gammq(K/2.,chi2/2.);*/
/*if (PERIODIC_TEMPLATES) {
fprintf(stats[TESTS_OVERLAPPING_TEMPL], "\t\t OVERLAPPING TEMPLATE OF ALL ONES TEST\n");
fprintf(stats[TESTS_OVERLAPPING_TEMPL], "\t\t-----------------------------------------------\n");
fprintf(stats[TESTS_OVERLAPPING_TEMPL], "\t\tCOMPUTATIONAL INFORMATION:\n");
fprintf(stats[TESTS_OVERLAPPING_TEMPL], "\t\t-----------------------------------------------\n");
fprintf(stats[TESTS_OVERLAPPING_TEMPL], "\t\t(a) n (sequence_length) = %d\n", n);
fprintf(stats[TESTS_OVERLAPPING_TEMPL], "\t\t(b) m (block length of 1s) = %d\n", m);
fprintf(stats[TESTS_OVERLAPPING_TEMPL], "\t\t(c) M (length of substring) = %d\n", M);
fprintf(stats[TESTS_OVERLAPPING_TEMPL], "\t\t(d) N (number of substrings) = %d\n", N);
fprintf(stats[TESTS_OVERLAPPING_TEMPL], "\t\t(e) lambda [(M-m+1)/2^m] = %f\n", lambda);
fprintf(stats[TESTS_OVERLAPPING_TEMPL], "\t\t(f) eta = %f\n", eta);
fprintf(stats[TESTS_OVERLAPPING_TEMPL], "\t\t-----------------------------------------------\n");
fprintf(stats[TESTS_OVERLAPPING_TEMPL], "\t\t F R E Q U E N C Y\n");
fprintf(stats[TESTS_OVERLAPPING_TEMPL], "\t\t 0 1 2 3 4 >=5 Chi^2 P-value Assignment\n");
fprintf(stats[TESTS_OVERLAPPING_TEMPL], "\t\t-----------------------------------------------\n");
fprintf(stats[TESTS_OVERLAPPING_TEMPL], "\t\t%3d %3d %3d %3d %3d %3d %f ",nu[0],nu[1],nu[2], nu[3], nu[4],nu[5],chi2);
}*/
if ( isNegative(p_value) || isGreaterThanOne(p_value) )
printf("WARNING: P_VALUE IS OUT OF RANGE.\n");
free(templates);
if ( p_value < ALPHA ) {
return false;
}
else {
return true;
}
/*free(templates);
fprintf(stats[TESTS_OVERLAPPING_TEMPL], "%f %s\n\n", p_value, assignment); fflush(stats[TESTS_OVERLAPPING_TEMPL]);
fprintf(results[TESTS_OVERLAPPING_TEMPL], "%f\n", p_value); fflush(results[TESTS_OVERLAPPING_TEMPL]);
fprintf(grid, "%d", state); fflush(grid);
if ( p_value < tp.minimumP )
tp.minimumP = p_value;
if ( !_isnan(p_value) )
tp.lnSum += log(p_value);
tp.df++;
#ifdef GEN_TIMING_INFO
finish = clock();
fprintf(fp, "%d\n", finish - start);
fclose(fp);
#endif
return;*/
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
R A N D O M E X C U R S I O N S T E S T
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
bool RandomExcursions(int n)
{
//char assignment[7];
int b, i, j, k, J, x;
int cycleStart, cycleStop, *cycle, *S_k;
int stateX[8] = {-4,-3,-2,-1,1,2,3,4};
int counter[8] = {0, 0, 0, 0, 0, 0, 0, 0};
double p_value, sum, constraint, nu[6][8];
double pi[10][6] = {
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
{0.5000000000, 0.2500000000, 0.1250000000, 0.06250000000, 0.03125000000, 0.0312500000},
{0.7500000000, 0.06250000000, 0.04687500000, 0.03515625000, 0.02636718750, 0.0791015625},
{0.8333333333, 0.02777777778, 0.02314814815, 0.01929012346, 0.01607510288, 0.0803755143},
{0.8750000000, 0.01562500000, 0.01367187500, 0.01196289063, 0.01046752930, 0.0732727051}
};
/*#ifdef GEN_TIMING_INFO
clock_t start, finish;
FILE *fp;
fp = fopen("RndExcursion.txt", "a");
start = clock();
#endif
*/
/*#if SAVE_RANDOM_EXCURSION_PARAMETERS == 1
if ( ((cycleInfo = fopen("cycleInfo", "a")) == NULL) ) {
fprintf(stats[TESTS_RANDOM_EXCURSIONS], "\t\tRandom Excursions: File cycleInfo could not ");
fprintf(stats[TESTS_RANDOM_EXCURSIONS], "be opened.\n");
exit(-1);
}
#endif*/
if ( ((S_k = (int *) calloc(n,sizeof(int))) == NULL) ||
((cycle = (int*) calloc(MAX(1000,n/200),sizeof(int))) == NULL) ) {
printf("Random Excursions Test: Insufficent Work Space Allocated.\n");
if ( S_k != NULL )
free(S_k);
if ( cycle != NULL )
free(cycle);
return false;
}
else {
J = 0; /* DETERMINE CYCLES */
S_k[0] = 2*(int)bit[0] - 1;
for( i=1; i<n; i++ ) {
S_k[i] = S_k[i-1] + 2*bit[i] - 1;
if ( S_k[i] == 0 ) {
J++;
if ( J > MAX(1000,n/128) ) {
printf("ERROR IN FUNCTION randomExcursions: EXCEEDING THE MAX");
printf(" NUMBER OF CYCLES EXPECTED\n.");
fflush(stdout);
free(S_k);
free(cycle);
//if ( cycleInfo != NULL )
// fclose(cycleInfo);
return false;
}
cycle[J] = i;
}
}
if ( S_k[n-1] != 0 )
J++;
/*#if SAVE_RANDOM_EXCURSION_PARAMETERS == 1
fprintf(cycleInfo,"%d\n", J);
#endif*/
/*if ( RANDOM_EXCURSIONS ) {
fprintf(stats[TESTS_RANDOM_EXCURSIONS], "\t\t\t RANDOM EXCURSIONS TEST\n");
fprintf(stats[TESTS_RANDOM_EXCURSIONS], "\t\t--------------------------------------------\n");
fprintf(stats[TESTS_RANDOM_EXCURSIONS], "\t\tCOMPUTATIONAL INFORMATION:\n");
fprintf(stats[TESTS_RANDOM_EXCURSIONS], "\t\t--------------------------------------------\n");
fprintf(stats[TESTS_RANDOM_EXCURSIONS], "\t\t(a) Number Of Cycles (J) = %04d\n", J);
fprintf(stats[TESTS_RANDOM_EXCURSIONS], "\t\t(b) Sequence Length (n) = %d\n", n);
}*/
constraint = MAX(0.005*pow(n,0.5),500);
if ( J < constraint ) {
printf("\t\t---------------------------------------------");
printf("\n\t\tWARNING: TEST NOT APPLICABLE. THERE ARE ");
printf("AN\n\t\t\t INSUFFICIENT NUMBER OF CYCLES.\n");
printf("\t\t---------------------------------------------");
printf("\n");
/* for( i=0; i<8; i++ ) {
fprintf(results[TESTS_RANDOM_EXCURSIONS], "%f\n", 0.0);
fprintf(grid, "%d", 0); fflush(grid);
fprintf(pvalsRndExc, "0.000000 "); fflush(pvalsRndExc);
}*/
}
else {
if ( RANDOM_EXCURSIONS ) {
printf("\t\t(c) Rejection Constraint = %f\n\t", constraint);
printf("\t-------------------------------------------\n");
}
cycleStart = 0;
cycleStop = cycle[1];
for( k=0; k<6; k++ )
for( i=0; i<8; i++ )
nu[k][i] = 0.0;
for( j=1; j<=J; j++ ) { /* FOR EACH CYCLE */
for( i=0; i<8; i++ )
counter[i] = 0;
for( i=cycleStart; i<cycleStop; i++ ) {
if ( ((S_k[i] >= 1) && (S_k[i] <= 4)) ||
((S_k[i] >= -4) && (S_k[i] <= -1)) ) {
if ( S_k[i] < 0 )
b = 4;
else
b = 3;
counter[S_k[i]+b]++;
}
}
cycleStart = cycle[j]+1;
if ( j < J )
cycleStop = cycle[j+1];
else
cycleStop = n;
for( i=0; i<8; i++ ) {
if ( counter[i] >= 0 && counter[i] <= 4 )
nu[counter[i]][i]++;
else if (counter[i] >= 5)
nu[5][i]++;
}
}
for( i=0; i<8; i++ ) {
x = stateX[i];
sum = 0.;
for( k=0; k<6; k++ )
sum += pow(nu[k][i] - J*pi[(int)fabs(x)][k],2) / (J*pi[(int)fabs(x)][k]);
p_value = igamc(2.5, sum/2.);
if ( isNegative(p_value) || isGreaterThanOne(p_value) )
printf("WARNING: P_VALUE IS OUT OF RANGE.\n");
free(S_k);
free(cycle);
if ( p_value < ALPHA ) {
return false;
}
else {
return true;
}
/* fprintf(stats[TESTS_RANDOM_EXCURSIONS], "%s\t\tx = %2d chi^2 = %9.6f p_value = %f\n",
assignment, x, sum, p_value); fflush(stats[TESTS_RANDOM_EXCURSIONS]);
fprintf(results[TESTS_RANDOM_EXCURSIONS], "%f\n", p_value); fflush(results[TESTS_RANDOM_EXCURSIONS]);
fprintf(grid, "%d", state); fflush(grid);
fprintf(pvalsRndExc, "%f ", p_value); fflush(pvalsRndExc);
if ( p_value < tp.minimumP )
tp.minimumP = p_value;
if ( !_isnan(p_value) )
tp.lnSum += log(p_value);
tp.df++;*/
}
}
// fprintf(stats[TESTS_RANDOM_EXCURSIONS], "\n");
}
/*#ifdef GEN_TIMING_INFO
finish = clock();
fprintf(fp, "%d\n", finish - start);
fclose(fp);
#endif
return;*/
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
R A N D O M E X C U R S I O N S V A R I A N T T E S T
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
bool RandomExcursionsVariant(int n)
{
int i, p, J, x, constraint;
double p_value;
int stateX[18] = {-9, -8, -7, -6, -5, -4, -3, -2, -1,
1, 2, 3, 4, 5, 6, 7, 8, 9};
// char assignment[7];
int count;
int* S_k;
/*
#ifdef GEN_TIMING_INFO
clock_t start, finish;
FILE *fp;
fp = fopen("RndExcVar.txt", "a");
start = clock();
#endif
*/
if ( (S_k = (int *) calloc(n,sizeof(int))) == NULL ) {
printf("\t\tRANDOM EXCURSIONS VARIANT: ");
printf("Insufficent memory allocated.\n");
return false;
}
else {
J = 0;
S_k[0] = 2*(int)bit[0] - 1;
for( i=1; i<n; i++ ) {
S_k[i] = S_k[i-1] + 2*bit[i] - 1;
if ( S_k[i] == 0 )
J++;
}
if ( S_k[n-1] != 0 )
J++;
/* if ( RANDOM_EXCURSIONS_VARIANT ) {
fprintf(stats[TESTS_RANDOM_EXCUR_VAR], "\t\t\tRANDOM EXCURSIONS VARIANT TEST\n");
fprintf(stats[TESTS_RANDOM_EXCUR_VAR], "\t\t--------------------------------------------\n");
fprintf(stats[TESTS_RANDOM_EXCUR_VAR], "\t\tCOMPUTATIONAL INFORMATION:\n");
fprintf(stats[TESTS_RANDOM_EXCUR_VAR], "\t\t--------------------------------------------\n");
fprintf(stats[TESTS_RANDOM_EXCUR_VAR], "\t\t(a) Number Of Cycles (J) = %d\n", J);
fprintf(stats[TESTS_RANDOM_EXCUR_VAR], "\t\t(b) Sequence Length (n) = %d\n", n);
fprintf(stats[TESTS_RANDOM_EXCUR_VAR], "\t\t--------------------------------------------\n");
}*/
constraint = (int)MAX(0.005*pow(n,0.5),500);
if ( J < constraint ) {
printf("\n\t\tWARNING: TEST NOT APPLICABLE. THERE ARE ");
printf("AN\n\t\t\t INSUFFICIENT NUMBER OF CYCLES.\n");
printf("\t\t---------------------------------------------");
printf("\n");
/*for( i=0; i<18; i++ ) {
fprintf(results[TESTS_RANDOM_EXCUR_VAR], "%f\n", 0.0);
fprintf(grid, "%d", 0); fflush(grid);
fprintf(pvalsRndExcVar, "0.000000 "); fflush(pvalsRndExcVar);*/
//}
}
else {
for( p=0; p<=17; p++ ) {
x = stateX[p];
count = 0;
for( i=0; i<n; i++ )
if ( S_k[i] == x )
count++;
p_value = erfc(fabs(count-J)/(sqrt(2.*J*(4.*fabs(x)-2))));
free(S_k);
if ( isNegative(p_value) || isGreaterThanOne(p_value) )
printf("\t\t(b) WARNING: P_VALUE IS OUT OF RANGE.\n");
if ( p_value < ALPHA ) {
return false;
}
else {
return true;
}
//fprintf(stats[TESTS_RANDOM_EXCUR_VAR], "%s\t\t(x = %2d) Total visits = %4d; p-value = %f\n",
// assignment, x, count, p_value); fflush(stats[TESTS_RANDOM_EXCUR_VAR]);
// fprintf(results[TESTS_RANDOM_EXCUR_VAR], "%1e\n", p_value); fflush(results[TESTS_RANDOM_EXCUR_VAR]);
// fprintf(grid, "%d", state); fflush(grid);
// fprintf(pvalsRndExcVar, "%f ", p_value); fflush(pvalsRndExcVar);
/* if ( p_value < tp.minimumP )
tp.minimumP = p_value;
if ( !_isnan(p_value) )
tp.lnSum += log(p_value);
tp.df++;*/
}
}
}
// fprintf(stats[TESTS_RANDOM_EXCUR_VAR], "\n");
/*
#ifdef GEN_TIMING_INFO
finish = clock();
fprintf(fp, "%d\n", finish - start);
fclose(fp);
#endif
return;*/
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
R A N K T E S T
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
bool Rank(int n)
{
int N = (int) floor(n/(32*32)); /* NUMBER OF MATRICES */
int r;
double p_value, product;
int i, k;
double chi_squared, arg1;
double p_32, p_31, p_30; /* PROBABILITIES */
double R; /* RANK */
double F_32, F_31, F_30; /* FREQUENCIES */
BitField** matrix = create_matrix(32, 32);
//char assignment[7];
/*#ifdef GEN_TIMING_INFO
clock_t start, finish;
FILE *fp;
fp = fopen("Rank.txt", "a");
start = clock();
#endif*/
if ( isZero(N) ) {
printf("\t\t\t\tRANK TEST\n");
printf("\t\tError: Insuffucient # Of Bits To Define An 32x32 ");
printf("(%dx%d) Matrix\n", 32, 32);
/*strcpy(assignment, "SUCCESS");
p_value = 0.00;
grid = 0;
state = 0;*/
return true;
}
else {
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -