/*--------------------------------------------------------------------------------
ランダム・シャッフル・サロゲート法
Random Shuffle (RS) Surrogate method
RS_surrogate.c
*--------------------------------------------------------------------------------*/
#include "RS_surrogate.h"
int main( void )
{
FILE *fo;
char filename[20];
int i, k;
int num;
unsigned long int seed = MT_SEED;
double *data;
double *RS;
/*--------------------------------------------------------------------------*/
// 乱数の初期設定
// initialisation of random numbers
init_genrand( seed );
/*--------------------------------------------------------------------------*/
if(( data = readdata( &num )) == NULL ) {
fputs( "readdata error\n", stderr );
return -1;
}
if (( RS = (double *)malloc((num) * sizeof(double))) == NULL ) {
fputs( "cannot allocate memory for RS\n", stderr );
return -1;
}
for ( k = 0; k < SURROGATE_DATA_NUM; k++ ) {
printf("# Now, generating the %dth random shuffle surrogate data", k+1 );
Sleep(10); /* 0.01秒 */
printf("\r");
if ( k < 9 ) {
sprintf(filename,"RS0%1d.dat", k+1);
}
if ( k >= 9 ) {
sprintf(filename,"RS%1d.dat", k+1);
}
if ((fo = fopen(filename,"w")) == NULL) {
printf("file open error\n");
exit(-1);
}
if ( random_shuffle_surrogate( num, data, RS ) != 0 ) {
fputs( "cannot generate random shuffle surrogate data\n", stderr );
return -1;
}
fprintf(fo, "%d\n", num );
for( i = 0; i < num; i++ ){
fprintf(fo, "%.10lf\n", RS[i]);
}
fclose(fo);
}
free( data );
free( RS );
return 0;
}
/*-------------------------------------------------------------------*/
int random_shuffle_surrogate( int num, double *data, double *RS )
{
int i, kk;
int *order;
unsigned long int j;
unsigned long int r;
unsigned long int org_order_sum;
unsigned long int RS_order_sum;
if(( order = (int*)malloc((num) * sizeof(int))) == NULL ) {
fputs("cannot allocate memory for order\n", stderr );
return -1;
}
org_order_sum = 0;
for( i = 0; i < num; i++ ) {
order[i] = i;
org_order_sum += i;
}
j = num;
RS_order_sum = 0;
for( i = 0; i < num; i++ ){
r = genrand_int32(); // 正の整数の一様分布, Uniform distribution of positive integers
r %= j;
kk = order[r];
order[r] = order[j-1];
j--;
RS[i] = data[kk];
RS_order_sum += kk;
// printf("%.10lf\n", RS[i]);
}
if ( org_order_sum != RS_order_sum ) {
// printf("%lu\n", org_order_sum);
// printf("%lu\n", RS_order_sum);
fputs( "there are problems with the generated order\n", stderr );
free( order );
return -1;
}
free( order );
return 0;
}