#include "timeseries.h"double color_only_mle(double *sigma, int n_data, double *Eig_values, double *Eig_vectors, double *residuals, int get_sigma) {	int j, k, i;	double alpha, beta, N, a, P, out, sig;	double *work1, *work2, *work3;	alpha = 1.0;	beta  = 0.0;	N = (double) n_data;	work1 = (double *) calloc( (size_t)(n_data * n_data),  sizeof(double) );	work2 = (double *) calloc( (size_t) n_data,            sizeof(double) );	work3 = (double *) calloc( (size_t) n_data,            sizeof(double) );	for (a = 0.0, j = 0; j < n_data; j++) a += log(Eig_values[j]);	for (j = 0; j < n_data; j++) {		for (k = 0; k < n_data; k++) {			work1[j + k * n_data] = Eig_vectors[k + j * n_data] / Eig_values[j];		}	}	i = 1;	dgemm_("T","N",&i,&n_data,&n_data,&alpha,residuals,&n_data,Eig_vectors,&n_data,&beta,work2,&i);	dgemm_("N","N",&i,&n_data,&n_data,&alpha,work2,&i,work1,&n_data,&beta,work3,&i);	for (P = 0.0, j = 0; j < n_data; j++) P += work3[j] * residuals[j];	free(work1);	free(work2);	free(work3);		if (get_sigma) {		sig = sqrt(P / N);		*sigma = sig;		}	else sig = *sigma;	out = -N * log(2.0 * M_PI) / 2.0 - a / 2.0 - N * log(sig*sig) / 2.0 - P / 2.0 / sig / sig;	return(out);}