cats.c

最大似然估计算法

#include "timeseries.h"
#include <getopt.h>

main(int argc, char **argv) {

int i, c, j, k, l, n, bf;
int p_alloc, n_periods;
int columns, error, n_harm;
int got_outfile, got_psdfile, got_kernelfile;
int n_models;
int filetype, no_elem;
int n_used_series;
int sm, n_psd;
int got_sf;
int estimate_trend;

char period_type;

char *outfile     = NULL;
char *psdfile     = NULL;
char *kernel_file = NULL;

char *model_string[100];
char word[512];

double a[3], p, scale_factor;

double *period, *P, *f;

time_series ts;

noise_model *nm;

options op;

data_kernel dk;

time_t start_time, end_time;
clock_t time1, time2;

struct tm *loctime;

struct utsname uts;
struct passwd *pw;

/*****************************/
/* Initialize some variables */
/*****************************/


got_sf         =  0;
filetype       =  1;
n_models       =  0;
error          =  0;
got_outfile    =  0;
got_psdfile    =  0;
got_kernelfile =  0;
estimate_trend =  1;
columns        = -1;

p_alloc   = 100;
n_periods = 0;
period = (double *) calloc((size_t)p_alloc, sizeof(double));

op.est_method   = 1;
op.cov_method   = 1;
op.verbose      = 0;
op.speed        = 0;
op.input_method = 5;

op.tol = (double *) calloc((size_t)4, sizeof(double)); 
op.tol[0]   = 400.0;
op.tol[1]   = 0.000001;
op.tol[2]   = 0.000020;
/* op.tol[3]   = 0.000001; */
op.tol[3]   = 0.0001; 

op.delta    = 1.5; 

op.fpout    = NULL;
op.fpsd     = NULL;
op.fpkernel = NULL;

scale_factor = 1.0;

/*******************/
/* Get the options */
/*******************/


while (1) {
	int option_index = 0;
	static struct option long_options[] = {
		{"model",     1, 0, 'M'},
		{"method",    1, 0, 'B'},
		{"columns",   1, 0, 'C'},
		{"scale",     1, 0, 'S'},
		{"sinusoid",  1, 0, 'A'},
		{"verbose",   0, 0, 'v'},
		{"Verbose",   0, 0, 'V'},
		{"speed",     1, 0, 'Z'},
		{"help",      0, 0, 'H'},
		{"quick",     1, 0, 'Q'},
		{"delta",     1, 0, 'D'},
		{"tolerance", 1, 0, 'T'},
		{"output",    1, 0, 'O'},
		{"filetype",  1, 0, 'F'},
		{"kernel",    1, 0, 'K'},
		{"psdfile",   1, 0, 'P'},
		{"notrend",   0, 0, 'E'},
		{"cov_form",  1, 0, 'X'}
	};

	c = getopt_long( argc, argv, "m:M:b:B:c:C:s:S:a:A:vVz:Z:hHq:Q:d:D:t:T:o:O:k:K:f:F:p:P:eEx:X:", long_options, &option_index);

	if (c == -1) break;

	switch (c) {
	
	case 0:
		break;
	case 'm':
	case 'M':
		model_string[n_models++] = optarg;
		break;
	case 'b':
	case 'B':
		switch (optarg[0]) {
			case 'm':
			case 'M':
				op.est_method = 1;
				break;
			case 's':
			case 'S':
				op.est_method = 2;
				break;
			case 'e':
			case 'E':
				op.est_method = 3;
				break;
			case 'w':
			case 'W':
				op.est_method = 4;
				break;
			default:
				fprintf(stderr, " Unknown estimation method : choose from mle, spectral, empirical, wls.\n"); 
				error++;
				break;
		}	
		break;
 	case 'c':
	case 'C':
		n = sscanf(optarg, "%d", &columns);
		if (n != 1) error++;
		break;
	case 's':
	case 'S':
		n = sscanf(optarg, "%lf", &scale_factor);
		if (n != 1) error++;
		got_sf = 1;
		break;
	case 'a':
	case 'A':
		n = sscanf(optarg, "%lf%c%d", &p, &period_type,&n_harm);
		if (n == 2) n_harm = 0;
		switch (period_type) {
			case 'd':
				for (j = 0; j <= n_harm; j++) {
					period[n_periods++] = p / (double) (j+1) / 365.24219;
					if (n_periods == p_alloc) {
						p_alloc += 16;
						period = (double *) realloc( (void *)period, p_alloc * sizeof(double) );
					}
				}
				break;
			case 'y':
				for (j = 0; j <= n_harm; j++) {
					period[n_periods++] = p / (double) (j+1);
					if (n_periods == p_alloc) {
						p_alloc += 16;
						period = (double *) realloc( (void *)period, p_alloc * sizeof(double) );
					}
				}
				break;
			case 'h':
				for (j = 0; j <= n_harm; j++) {
					period[n_periods++] = p / (double) (j+1) / 24.0 / 365.24219;
					if (n_periods == p_alloc) {
						p_alloc += 16;
						period = (double *) realloc( (void *)period, p_alloc * sizeof(double) );
					}
				}
				break;
			case 'm':
				for (j = 0; j <= n_harm; j++) {
					period[n_periods++] = p / (double) (j+1) / 24.0 / 60.0 / 365.24219;
					if (n_periods == p_alloc) {
						p_alloc += 16;
						period = (double *) realloc( (void *)period, p_alloc * sizeof(double) );
					}
				}
				break;
			case 's':
				for (j = 0; j <= n_harm; j++) {
					period[n_periods++] = p / (double) (j+1) / 24.0 / 3600.0 / 365.24219;
					if (n_periods == p_alloc) {
						p_alloc += 16;
						period = (double *) realloc( (void *)period, p_alloc * sizeof(double) );
					}
				}
				break;
			default:
				error++;
				fprintf(stderr, "Error defining a sinusoid period\n");
		}
		break;
	case 'v':
		op.verbose = 1;
		break;
	case 'V':
		op.verbose = 2;
		break;
	case 'x':
	case 'X':
		n = sscanf(optarg, "%d", &op.cov_method);
		if (n != 1) error++;
		if (op.cov_method != 1 || op.cov_method != 2) error++;
		break;
	case 'z':
	case 'Z':
		n = sscanf(optarg, "%d", &op.speed);
		if (n != 1) error++;
		if (op.speed < 0 || op.speed > 3) error++;
		break;
	
	case 'h':
	case 'H':
		error++;
		break;
	case 'd':
	case 'D':
		n = sscanf(optarg, "%lf", &op.delta);
		if (n != 1) error++;
		if (op.delta <= 1.0) error++;
		break;
	case 't':
	case 'T':
		n = sscanf (optarg, "%lf/%lf/%lf", &op.tol[1], &op.tol[2], &op.tol[3]);
		if (n != 3) {
			error++;
			fprintf(stderr, "Must specify 3 parameters -T tol1/tol2/tol3\n");
		}
		break;
	case 'o':
	case 'O':
		outfile = optarg;
		got_outfile = 1;
		break;
	case 'k':
	case 'K':
		kernel_file = optarg;
		got_kernelfile = 1;
		break;
	case 'f':
	case 'F':
		if (strncmp(optarg,"cats",4) == 0) {
			filetype = 1;
		} else if (strncmp(optarg,"psmsl",4) == 0) {
			filetype = 2;
		} else {
			fprintf(stderr, " filetype %s must be cats or psmsl.\n", optarg);
			error++;
		}
		break;
	case 'p':
	case 'P':
		psdfile = optarg;
		got_psdfile = 1;
		break;
	case 'e':
	case 'E':
		estimate_trend = 0;
		break;
	}
}

if (got_outfile) {
	if ((op.fpout = fopen(outfile,"a+")) == NULL) {
		fprintf(stderr, "%s : Cannot open file %s\n", argv[0], outfile);
		exit(EXIT_FAILURE);
	}
} else op.fpout = stdout;

if (got_psdfile) {
	if ((op.fpsd = fopen(psdfile,"w")) == NULL) {
		fprintf(stderr, "%s : Cannot open file %s\n", argv[0], psdfile);
		exit(EXIT_FAILURE);
	}
}

if (op.verbose) {
	fprintf(op.fpout, " Cats Version : %s\n", version);
	fprintf(op.fpout, " Cats command :");
	for (i = 0; i < argc; i++) fprintf(op.fpout, " %s", argv[i]);
	fprintf(op.fpout, "\n");
}

if (error) {
	fprintf(stderr, " Cats Version : %s\n", version);
	fprintf(stderr, " Cats command :");
	for (i = 0; i < argc; i++) fprintf(stderr, " %s", argv[i]);
	fprintf(stderr, "\n Please look in Manual for command syntax.\n");
	exit(EXIT_FAILURE);
}

/************************************************************************/
/* Now output some information about file, machine and operating system */
/* Useful for comparisons                                               */
/* Also who ran the program                                             */
/************************************************************************/

uname(&uts);
fprintf(op.fpout, " Data from file : %s\n", argv[optind]);
fprintf(op.fpout, " %s : running on %s\n", argv[0], uts.nodename);
fprintf(op.fpout, " %s release %s (version %s) on %s\n", uts.sysname, uts.release, uts.version, uts.machine);

if ((pw = getpwuid (getuid ())) != NULL) {
        fprintf(op.fpout, " userid : %s\n\n", pw->pw_name);
} else {
        fprintf(op.fpout, " userid : unknown\n\n");
}

no_elem = argc-optind;
if (no_elem > 0) {
	if (no_elem > 1) fprintf(stderr, " There are %d non-option elements : using the first as the input file\n", no_elem);

	if (filetype == 1) {
		if (!got_sf) scale_factor = 1000.0;
		ts = read_cats_series(argv[optind],op,scale_factor,columns);
	} else if (filetype == 2) {
		if (!got_sf) scale_factor = 1.0;
		ts = read_rlrdata(argv[optind],op,scale_factor); 
	} else {
		fprintf(stderr, "Unrecognized input file type\n");
		error++;
		/* exit(EXIT_FAILURE); */
	}
} else {
	fprintf(stderr, " An input file hasnt been specified\n");
	error++;
	/* exit(EXIT_FAILURE); */
}

if (error) {
	fprintf(stderr, " Cats Version : %s\n", version);
	fprintf(stderr, " Cats command :");
	for (i = 0; i < argc; i++) fprintf(stderr, " %s", argv[i]);
	fprintf(stderr, "\n Please look in Manual for command syntax.\n");
	exit(EXIT_FAILURE);
}

/*****************************************/
/* count the number of series being used */
/*****************************************/

fprintf(op.fpout, " Number of series to process : %d\n", ts.n_series);

/****************************************/
/* Now presumably we are ready to go... */
/****************************************/

start_time = time(NULL);
time1 = clock();
loctime = localtime(&start_time);
fprintf(op.fpout, " Start Time : %s\n", asctime(loctime));

/*******************************************/
/* The main loop, loops through the series */
/*******************************************/

for (c = 0; c < ts.n_series; c++) {

	/********************************************/
	/* Now process the stochastic models to use */
	/********************************************/

	nm = (noise_model *) calloc((size_t) (n_models), sizeof(noise_model));
	for (j = 0; j < n_models; j++) {
		nm[j] = decode_model_string(model_string[j],c);
		nm[j].c_id = ts.c_id[c];
		print_model(nm[j], op, 1.0);
	}
	
	ts.current_series = c;

	/******************************************************************/
	/* At some point add in the ability to read data kernel from file */
	/******************************************************************/

	if (got_kernelfile) {
		/* dk = create_A_and_d_from_file(ts, c, kernel_file); */
		dk = create_A_and_d(ts, period, n_periods, c, estimate_trend);
	} else {
		dk = create_A_and_d(ts, period, n_periods, c, estimate_trend);
	}

	if (got_psdfile || op.est_method == 2) {
		P = (double *) calloc((size_t) dk.n_data, sizeof(double));
		f = (double *) calloc((size_t) dk.n_data, sizeof(double));
		n_psd = create_spectrum(P, f, ts, dk, op);
		P = (double *) realloc( (void *)P, n_psd * sizeof(double));
		f = (double *) realloc( (void *)f, n_psd * sizeof(double));
		if (got_psdfile) {
			for (j = 0; j < n_psd; j++) fprintf(op.fpsd,"%g %g\n", f[j], P[j] );
			fprintf(op.fpsd, "0 0\n");
		}
	}

	
	/*********************************************************************/
	/* Now check the model options against the estimation method options */
	/*********************************************************************/

	switch(op.est_method) {
		case 1:
			mle_wrapper(nm,n_models,ts,dk,op);	
			break;
		case 2:
			psd_wrapper(nm,n_models,ts,dk,op,P,f,n_psd);
			break;
		case 3: 
			emp_wrapper(nm,n_models,ts,dk,op);
			break;
		case 4:
			wls_wrapper(nm,n_models,ts,dk,op); 
			break;
	}	
	if (got_psdfile || op.est_method == 2) {
		free(P);
		free(f);
	}

	/***********************/
	/* Output results here */
	/***********************/

	fprintf(op.fpout, "+%s  MLE    :    %12.6f\n", comp_names[ts.c_id[c]], dk.MLE[0]); 
	for (j = 0; j < dk.n_par; j++) {
		fprintf(op.fpout, "+%s ", comp_names[ts.c_id[c]]);
		fprintf(op.fpout, " %s ", par_names[dk.name_id[j]]);
		fprintf(op.fpout, " %11.4f +- %11.4f\n", dk.params[j], sqrt(dk.covar[j + j * dk.n_par]));
	}
	for (j = 0; j < n_models; j++) print_model(nm[j],op,1.0);
		
	for (j = 0; j < n_models; j++) { 
		if (nm[j].n_pvec > 0) {
			free(nm[j].names);
			free(nm[j].pvec_flag);
			free(nm[j].pvec);
			free(nm[j].pvec_sigma);
		}
		free(nm[j].sigma);
		free(nm[j].dsigma);
		free(nm[j].sigma_flag);
	}

	free(nm);

	free(dk.A);
	free(dk.d);
	free(dk.dd);
	free(dk.name_id);
	free(dk.alt_id);
	free(dk.params);
	free(dk.covar);
	free(dk.MLE);

}

free(ts.index);
free(ts.off_code);
free(ts.c_id);
free(ts.t);
free(ts.data);
free(ts.formal_error);
free(ts.offsets);
free(period);
free(op.tol);

end_time = time(NULL);
time2 = clock();
loctime = localtime(&end_time);
fprintf(op.fpout, " End Time : %s\n", asctime(loctime));
fprintf(op.fpout, " Total Time : %lg\n", difftime(end_time, start_time) );


if (op.fpout != stdout) fclose(op.fpout);
if (got_psdfile) fclose(op.fpsd);

}