cats_sopt.c.old

最大似然估计算法

#include "timeseries.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;
int n_models, white;
int filetype, no_elem;
int estimation_method;
int n_used_series;
int fno;
int got_white, sm, n_psd;
int estimate_trend;

char period_type;

char *outfile = NULL;
char *psdfile = 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, time1, time2, end_time;

struct tm *loctime;

struct utsname uts;
struct passwd *pw;

FILE *fpsd;
/*****************************/
/* Initialize some variables */
/*****************************/

estimation_method = 1;

fno            = 0;
filetype       = 1;
white          = 1;
n_models       = 0;
error          = 0;
got_outfile    = 0;
got_psdfile    = 0;
estimate_trend = 1;
columns        = -1;

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

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.003;
op.tol[2] = 0.00002;
op.tol[3] = 0.000001;
op.delta  = 1.1; 
op.fpout  = NULL;

scale_factor = 1.0;

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


for (i = 1; i < argc; i++) {
	if (argv[i][0] == '-') {
		switch(argv[i][1]) {
		case 'm':
		case 'M':
			model_string[n_models++] = &argv[i][2];
			break;
		case 'n':
		case 'N':
			white = 0;
			break;
		case 'b':
		case 'B':
			switch (argv[i][2]) {
				case 'm':
				case 'M':
					estimation_method = 1;
					break;
				case 'e':
				case 'E':
					estimation_method = 2;
					break;
				case 's':
				case 'S':
					estimation_method = 3;
					break;
				case 'w':
				case 'W':
					estimation_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(&argv[i][2], "%d", &columns);
			if (n != 1) error++;
			break;
		case 's':
		case 'S':
			n = sscanf(&argv[i][2], "%lf", &scale_factor);
			if (n != 1) error++;
			break;
		case 'a':
		case 'A':
			n = sscanf(&argv[i][2], "%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':
		case 'V':
			op.verbose = 1;
			break;
		case 'x':
		case 'X':
			n = sscanf(&argv[i][2], "%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(&argv[i][2], "%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(&argv[i][2], "%lf", &op.delta);
			if (n != 1) error++;
			if (op.delta <= 1.0) error++;
			break;
		case 't':
		case 'T':
			n = sscanf (&argv[i][2], "%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 = &argv[i][2];
			got_outfile = 1;
			break;
		case 'k':
		case 'K':
			/* Currently not used */
			break;
		case 'f':
		case 'F':
			if (strncmp(&argv[i][2],"cats",4) == 0) {
				filetype = 1;
			} else if (strncmp(&argv[i][2],"psmsl",4) == 0) {
				filetype = 2;
			} else {
				fprintf(stderr, " filetype %s must be cats or psmsl.\n", optarg);
				error++;
			}
			break;
		case 'p':
		case 'P':
			psdfile = &argv[i][2];
			got_psdfile = 1;
			break;
		case 'e':
		case 'E':
			estimate_trend = 0;
			break;
		default:
			error = 1;
			break;
		}
	} else {
		if (fno == 0) fno = i;
	}
}
	
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 (fno > 0) {
	if (filetype == 1) {
		ts = read_cats_series(argv[fno],op,scale_factor);
	} else if (filetype == 2) {
		ts = read_rlrdata(argv[fno],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, " Print manual here\n");
	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");
}

/**********************************************************/
/* Sort out which columns of data are actually being used */
/**********************************************************/

if (columns > 0) {
	for (k = 0; k < ts.n_series; k++) {
		ts.series_flag[k] = 1;
		bf = (int) pow(2.0, (double) ts.n_series-k-1);
		if ((columns & bf) == 0) ts.series_flag[k] = 0;
		if (op.verbose) fprintf(op.fpout, " Series[%d] = %d\n", k, ts.series_flag[k]);
	}
}

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

for (n_used_series = 0, k = 0; k < ts.n_series; k++) n_used_series += ts.series_flag[k];

if (n_used_series == 0) {
	fprintf(stderr, " No series in file to be used! exiting...\n");
	exit(EXIT_FAILURE);
}
fprintf(op.fpout, " Number of series to process : %d\n", n_used_series);


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


for (got_white = 0, j = 0; j < n_models; j++) if (model_string[j][0] == 'w') got_white = 1;
if (white && !got_white) model_string[n_models++] = "wh:";
if (got_white && !white) {
	fprintf(stderr, " You have specified white noise as a model but requested not to estimate a white noise component!\n");
	exit(EXIT_FAILURE);
}

if (n_models - white > 1 && estimation_method != 4) {
	fprintf(stderr, " cats can only cope with one stochastic model at a time at present.\n");
	exit(EXIT_FAILURE);
}

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],n_used_series);
}

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

uname(&uts);
fprintf(op.fpout, " Data from file : %s\n", argv[fno]);
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");
}

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

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

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

	if (ts.series_flag[c] == 0) continue;
	ts.current_series = c;

	dk = create_A_and_d(ts, period, n_periods, c, estimate_trend);

	/* error checking on options */

	error = 0;
	for (j = 0; j < n_models; j++) {
		nm[j].current_field = c;
		if (nm[j].model == 'f' && estimation_method == 2) error++; 
		if (nm[j].model == 'b' && (estimation_method == 2 || estimation_method == 3)) error++; 
		if (estimation_method == 4 && !nm[j].got_sigma) error++; 
	}
	if (error) continue;

	if (got_psdfile || estimation_method == 3) {
		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);
		fflush(op.fpout);
	}
	if (got_psdfile) {
		for (j = 0; j < n_psd; j++) {
			fprintf(fpsd, "%d %g %g\n", c, f[j], P[j]);
		}
	}

	switch(estimation_method) {
		case 1:
			mle_wrapper(nm,n_models,ts,dk,op);	
			break;
		case 2:
			/* emp_wrapper(nm,n_models,ts,dk,op); */
			break;
		case 3: 
			fprintf(stderr, "Calling psd_wrapper\n");
			psd_wrapper(nm,n_models,ts,dk,op,P,f,n_psd); 
			fprintf(stderr, "Calling psd_wrapper\n");
			break;
		case 4:
			/* wls_wrapper(nm,n_models,ts,dk,op); */
			break;
	}	
	if (got_psdfile || estimation_method == 3) {
		free(P);
		free(f);
	}
	free(dk.A);
	free(dk.d);
	free(dk.name_id);
	free(dk.alt_id);

}

if (got_psdfile) fclose(fpsd);

end_time = time(NULL);
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) ); 

}