calculate_fs2.c

来自「最大似然估计算法」· C语言 代码 · 共 131 行

#include "timeseries.h"

double calculate_fs2(double *time, int n_data, int *index, int *f1_count) {

	int i, j, k, n;

	int error, count;
	int year1, year2, leap1, leap2;

	int n_attempt = 9;
	int n_sample  = 7;

	int *bad, *idt;

	double s, dp, min_dt, max_dt, scale, t, ds;
	double median, fs, rdt;

	double *dt, *tmp;

	double sample[] = {0.041666666, 1.0, 2.0, 7.0, 14.0, 30.4375, 365.25};
	double daily[]  = {365.24219, 0.0, 0.0, 0.0, 365.0, 365.25, 365.24219, 365.2564, 365.259635};

	dt     = (double *) calloc((size_t) (n_data-1),  sizeof(double));
	tmp    = (double *) calloc((size_t) (n_data-1),  sizeof(double));
	bad    = (int *)    calloc((size_t) (n_data-1),  sizeof(int));
	idt    = (int *)    calloc((size_t) (n_data-1),  sizeof(int));

	error = 0;
	for (j = 0; j < n_data-1; j++) {

		dt[j] = time[j+1] - time[j];
		tmp[j] = dt[j];

		if (dt[j] <= 0.0) bad[error++] = j;

		if (j == 0) min_dt = dt[j];
		else min_dt = dt[j] < min_dt ? dt[j] : min_dt; 
	}
	if (error) {
		fprintf(stderr, " calculate_fs : file is not in ascending");
		fprintf(stderr, " time order or has duplicated points\n");
		for (j = 0; j < error; j++) {
		       	fprintf(stderr, "              : t[%d] - t[%d] = %f\n", bad[j]+1, bad[j], dt[bad[j]]);
		}
		exit(EXIT_FAILURE);
	}
	free(bad);

	/******************************************************************/
	/* First find out the number of decimal places in the time stamps */
	/******************************************************************/

	dp = 0.0;
	do {
		dp += 1.0;
		s = 0.0;
		scale = pow(10.0,dp);
		for (j = 0; j < n_data; j++) {
			t = rint(time[j] * scale) / scale;
			s += ((t-time[j])*(t-time[j]));
		}
		s = sqrt(s / (double) n_data);	
	} while (s > 1e-13);

	/************************************************************************/
	/* The 1e-13 was arrived at partly from trial and error and it is       */
	/* related to the numerical precision of the computer, meaning it       */
	/* can only really distinguish up to 12 decimal places                  */
	/* For a specified decimal place, the differences between the full      */
	/* value and the "clipped" value should form a uniform distribution     */
	/* between -0.5 / 10^dp and 0.5 / 10^dp. Therefore the square root of   */
	/* the sum of the differences squared should have a mean equal to       */
	/* (N / 12 / 10^(2*dp))^0.5. The difficultly is for small time series   */
	/* say less the 250 the error on the mean is quite large and could      */
	/* possibly overlap with the dp level below. A better method will be    */
	/* developed later, its not really critical at this stage, for well     */
	/* defined time series.                                                 */
	/************************************************************************/

	count = n_data - 1;
	qsort((char *)tmp, (size_t)count, sizeof(double), comp_double_asc);
	if ((int)fmod(count,2) == 1) {
		median = tmp[(count+1)/2 - 1];
	} else {
                median = ( tmp[count/2-1] + tmp[count/2] ) / 2.0;
	}
	free(tmp);

	n = lrint(median / min_dt);
	rdt = fmod(median,min_dt);

	if (n > 1) {
 		if (rdt < 0.1) {
			median /= (double) n;
		}
		else {
			fprintf(stderr, " Error : calculate_fs2 : min_dt and median disagree");
			fprintf(stderr, " by non-integer value (remainder %f) greater than 0.1\n", rdt);
			exit(EXIT_FAILURE);
		}
	}
		
	count = 0;
	error = 0;
	index[0] = 0;
	for (j = 0; j < n_data-1; j++) {
	       	idt[j] = lrint( dt[j] / median);
		if (idt[j] == 0) {
			fprintf(stderr, " Error : calculate_fs2 : dt (%f) is smaller", dt[j]); 
			fprintf(stderr, " than median (%f) at %4d,%4d (%f %f) \n",median, j, j+1, time[j], time[j+1]);
			 error++;
		}
		if (idt[j] == 1) count++;
		index[j+1] = index[j] + idt[j];
	}

	if (error) {
		fprintf(stderr, "%d\n", lrint(median / min_dt));
		exit(EXIT_FAILURE);
	}
	
	fs = 1.0 / median / 365.24219 / 24.0 / 3600; 
	
	(*f1_count) = count;

	free(dt);
	free(idt);

	return(fs);
}