#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);
}
موضوع: نحوه اجرا کردن یک برنامه که در زبان c نوشته شده (دفعات بازدید: 8770 بار)
