read_cats_series.c

最大似然估计算法

#include "timeseries.h"

time_series read_cats_series(char *filename, options op, double scale_factor, int columns) {

	FILE *fpin;
	time_series ts;

	int j, k, n, ns, bf;

	int n_offset;
	int n_data;

	int n_alloc   = 512;
	int off_alloc = 16;

	int *flag;

	double tt, no, ea, up, dno, dea, dup;

	char line[512];

	if ((fpin = fopen(filename, "r")) == NULL) {
		fprintf(stderr, " read_cats_series : Cannot open file %s\n", filename);
		exit(EXIT_FAILURE);
	}
	if (columns == -1) columns = 7;

	if (columns == 0) {
		fprintf(stderr, " read_cats_series : Columns = 0 means there are no series to process!\n");
		exit(EXIT_FAILURE);
	}

	flag = (int *) calloc((size_t) 3, sizeof(int));
	ts.n_series = 0;
	for (k = 0; k < 3; k++) {
        	flag[k] = 1;
        	bf = (int) pow(2.0, (double) 2-k);
        	if ((columns & bf) == 0) flag[k] = 0;
		else ts.n_series++;
	}

	ns = ts.n_series;
	ts.got_errors = 1;

	ts.t            = (double *) calloc((size_t)n_alloc      , sizeof(double));
	ts.data	        = (double *) calloc((size_t)(n_alloc*ns) , sizeof(double));
	ts.formal_error = (double *) calloc((size_t)(n_alloc*ns) , sizeof(double));
	ts.offsets      = (double *) calloc((size_t)(off_alloc)  , sizeof(double));
	ts.off_code     = (int *)    calloc((size_t)(off_alloc)  , sizeof(int));

	n_data   = 0;
	n_offset = 0;

	while (fgets (line, 512, fpin)) {
		if (strncmp(line,"# offs",6) == 0) {
			n = sscanf(line, "# offset %lf %d", &ts.offsets[n_offset], &ts.off_code[n_offset]);
	       		if (n == 1) ts.off_code[n_offset] = 7;
			ts.off_code[n_offset] = ts.off_code[n_offset] > 7 ? 7 : ts.off_code[n_offset];
			ts.off_code[n_offset] = ts.off_code[n_offset] < 0 ? 0 : ts.off_code[n_offset];
			n_offset++;
			if (n_offset == off_alloc) {
				off_alloc += 10;
				ts.offsets  = (double *) realloc( (void *)ts.offsets,  off_alloc * sizeof(double) );
				ts.off_code = (int *)    realloc( (void *)ts.off_code, off_alloc * sizeof(int) );
			}
		} else if (strncmp(line,"#",1) == 0) continue;
		else {
			sscanf(line, "%lf %lf %lf %lf %lf %lf %lf", &tt, &no, &ea, &up, &dno, &dea, &dup);
			ts.t[n_data] = tt;

			k = 0;
			if (flag[0]) {
				ts.data[(n_data*ns)+k]         =  no * scale_factor;
				ts.formal_error[(n_data*ns)+k] = dno * scale_factor;
				k++;
			} 
			if (flag[1]) {
				ts.data[(n_data*ns)+k]         =  ea * scale_factor;
				ts.formal_error[(n_data*ns)+k] = dea * scale_factor;
				k++;
			}
			if (flag[2]) {
				ts.data[(n_data*ns)+k]         =  up * scale_factor;
				ts.formal_error[(n_data*ns)+k] = dup * scale_factor;
				k++;
			}
			n_data++;
		}
		if (n_data == n_alloc) {
			n_alloc += 512;
			ts.t            = (double *) realloc( (void *)ts.t,     n_alloc * sizeof(double) );
			ts.data         = (double *) realloc( (void *)ts.data,  n_alloc * ts.n_series * sizeof(double) );
			ts.formal_error = (double *) realloc( (void *)ts.formal_error,    n_alloc * ts.n_series * sizeof(double) );
		}
	}

	ts.n_data    = n_data;
	ts.n_offsets = n_offset;

	ts.index       = (int *)    calloc((size_t)(n_data),           sizeof(int));
	ts.c_id        = (int *)    calloc((size_t)(ts.n_series),      sizeof(int));

	ns = ts.n_series;
	ts.t            = (double *) realloc( (void *)ts.t,            n_data *     sizeof(double));
	ts.data         = (double *) realloc( (void *)ts.data,         n_data * ns * sizeof(double));
	ts.formal_error = (double *) realloc( (void *)ts.formal_error, n_data * ns * sizeof(double));
	ts.offsets      = (double *) realloc( (void *)ts.offsets,      n_offset *   sizeof(double));
	ts.off_code     = (int *)    realloc( (void *)ts.off_code,     n_offset *   sizeof(int));

	for (k = 0, j = 0; j < 3; j++) { 
		if (flag[j]) ts.c_id[k++] = j;
	}

	ts.fs = calculate_fs2(ts.t, ts.n_data, ts.index, &ts.count);

	ts.n_full = ts.index[ts.n_data-1]+1;
	ts.time_span = ts.t[n_data-1] - ts.t[0];

	if (op.verbose) fprintf(op.fpout, " Sampling frequency %.6g (Hz), %.2f days\n", ts.fs, 1.0 / ts.fs / 24.0 / 3600.0);
	if (op.verbose) fprintf(op.fpout, " Number of samples 1 period apart = %d of %d\n", ts.count, ts.n_data-1);
	if (op.verbose) fprintf(op.fpout, " Number of points in full series  = %d\n", ts.n_full);

	free(flag);

	return(ts);
}