create_covariance.c

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

#include "timeseries.h"

void create_covariance(noise_model nm, time_series ts, int method, double *C, int is_unit) {

int j, k, l;
int n_pts, n_full, n_squared;

double dt, sec_per_year;
double alpha, gamma, s;
double fl, fh, td, e;

double *tv, *Q;

FILE *fp;

alpha = 1.0;
gamma = 0.0;
sec_per_year = 365.24219*24.0*3600.0;

dt = 1.0 / ts.fs / sec_per_year;

n_pts  = ts.n_data;
n_full = ts.n_full;
n_squared = n_pts * n_pts;

switch (nm.model) {
	case 'b':
		tv = (double *) calloc((size_t) n_full, sizeof(double));

		fl = nm.pvec[0] / nm.pvec[1];
		fh = nm.pvec[0] * nm.pvec[1];
		band_pass_tv(tv, n_full, fl, fh, nm.pvec[2]); 
		alpha =  is_unit ? 1.0  : nm.sigma[0] * nm.sigma[0];

		Q = (double *)calloc((size_t)(n_pts*n_full), sizeof(double));
		for (k = 0; k < n_pts; k++) {
			l = ts.index[k];	
			for (j = 0; j <= l; j++) {
				Q[k + j * n_pts] = tv[l-j];
			}
		}

        	dgemm_("N","T",&n_pts,&n_pts,&n_full,&alpha,Q,&n_pts,Q,&n_pts,&gamma,C,&n_pts);

		free(Q);
		free(tv);

		break;

	case 'f':
		alpha = is_unit ? 1.0  : nm.sigma[0] * nm.sigma[0];

		if (method == 1) {
			tv = (double *) calloc((size_t) n_full, sizeof(double));

			fogm_tv(tv, n_full, nm.pvec[0], dt);

			Q = (double *)calloc((size_t)(n_pts*n_full), sizeof(double));
			for (k = 0; k < n_pts; k++) {
				l = ts.index[k];	
				for (j = 0; j <= l; j++) {
					Q[k + j * n_pts] = tv[l-j];
				}
			}

        		dgemm_("N","T",&n_pts,&n_pts,&n_full,&alpha,Q,&n_pts,Q,&n_pts,&gamma,C,&n_pts);

			free(Q);
			free(tv);
		} else {
			generate_FOGM_cov(ts.t, ts.n_data, dt, alpha, nm.pvec[0], C); 
		}
		break;

	case 'p':
		alpha = is_unit ? pow(dt,-nm.pvec[0]/2.0) : nm.sigma[0] * nm.sigma[0] * pow(dt,-nm.pvec[0]/2.0);
	
		if (method == 1) {
	
			tv = (double *) calloc((size_t) n_full, sizeof(double));

			power_law_tv(tv, n_full, nm.pvec[0]/2.0); 

			Q = (double *)calloc((size_t)(n_pts*n_full), sizeof(double));
			for (k = 0; k < n_pts; k++) {
				l = ts.index[k];	
				for (j = 0; j <= l; j++) {
					Q[k + j * n_pts] = tv[l-j];
				}
			}

        		dgemm_("N","T",&n_pts,&n_pts,&n_full,&alpha,Q,&n_pts,Q,&n_pts,&gamma,C,&n_pts);

			free(Q);
			free(tv);

		} else {
			power_law_cov(ts.t, ts.n_data, dt,  alpha, nm.pvec[0], C);
		}
		break;
	case 's':

		alpha = is_unit ? 1.0  : nm.sigma[0] * nm.sigma[0];

		for (k = 0; k < n_squared; k++) C[k] = 0.0;

		for (k = 0; k < n_pts; k++) {
			if (ts.t[k] >= nm.pvec[0] && ts.t[k] <= nm.pvec[1]) {
				C[k + k * n_pts] = alpha;
			} else {
				C[k + k * n_pts] = 1e-10;
			}
		}

		break;

	case 't':
		alpha = is_unit ? 1.0  : nm.sigma[0] * nm.sigma[0];

		for (k = 0; k < n_squared; k++) C[k] = 0.0;

		for (k = 0; k < n_pts; k++) {
			td = ts.t[k] - ts.t[0];
			e = exp( nm.pvec[0] * td);
			C[k + k * n_pts] = alpha * e * e; 
		}

		break;
		

	case 'v':

		alpha = is_unit ? 1.0  : nm.sigma[0] * nm.sigma[0];
		l = ts.current_series;

		for (k = 0; k < n_squared; k++) C[k] = 0.0;

		for (k = 0; k < n_pts; k++) {
				s = ts.formal_error[(k*ts.n_series)+l];
				C[k + k * n_pts] = alpha * s * s;
		}

		break;
	case 'w':

		alpha = is_unit ? 1.0  : nm.sigma[0] * nm.sigma[0];

		for (k = 0; k < n_squared; k++) C[k] = 0.0;

		for (k = 0; k < n_pts; k++) C[k + k * n_pts] = alpha;

		break;
	case 'g':

		alpha = is_unit ? pow(dt,-nm.pvec[1]/2.0) : nm.sigma[0] * nm.sigma[0] * pow(dt,-nm.pvec[1]/2.0);

		tv = (double *) calloc((size_t) n_full, sizeof(double));

		fogm_tv(tv, n_full, nm.pvec[0], dt);

		Q = (double *)calloc((size_t)(n_pts*n_full), sizeof(double));

		for (k = 0; k < n_pts; k++) {
			l = ts.index[k];	
			for (j = 0; j <= l; j++) {
				Q[k + j * n_pts] = tv[l-j] / tv[0];
			}
		}

		free(tv);
		tv = (double *) calloc((size_t) n_full, sizeof(double));

		power_law_tv(tv, n_full, nm.pvec[1]/2.0); 

		for (k = 0; k < n_pts; k++) {
			l = ts.index[k];	
			for (j = 0; j <= l; j++) {
				Q[k + j * n_pts] *= tv[l-j];
			}
		}
        	dgemm_("N","T",&n_pts,&n_pts,&n_full,&alpha,Q,&n_pts,Q,&n_pts,&gamma,C,&n_pts);

		free(Q);
		free(tv);
		break;
		
	default:
		fprintf(stderr, " create_covariance : unknown model (%s) /method for covariance\n", nm.model);
		exit(EXIT_FAILURE);
		break;
}

}