text:
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <math.h>
#include "cconfigspace_internal.h"
#include "distribution_internal.h"

struct _ccs_distribution_roulette_data_s {
	_ccs_distribution_common_data_t  common_data;
	size_t                           num_areas;
	ccs_float_t                     *areas;
};
typedef struct _ccs_distribution_roulette_data_s _ccs_distribution_roulette_data_t;

static ccs_result_t
_ccs_distribution_del(ccs_object_t o) {
	(void)o;
	return CCS_SUCCESS;
}

static ccs_result_t
_ccs_distribution_roulette_get_bounds(_ccs_distribution_data_t *data,
                                      ccs_interval_t           *interval_ret);

static ccs_result_t
_ccs_distribution_roulette_samples(_ccs_distribution_data_t *data,
                                   ccs_rng_t                 rng,
                                   size_t                    num_values,
                                   ccs_numeric_t            *values);

static ccs_result_t
_ccs_distribution_roulette_strided_samples(_ccs_distribution_data_t *data,
                                           ccs_rng_t                 rng,
                                           size_t                    num_values,
                                           size_t                    stride,
                                           ccs_numeric_t            *values);

static ccs_result_t
_ccs_distribution_roulette_soa_samples(_ccs_distribution_data_t  *data,
                                       ccs_rng_t                  rng,
                                       size_t                     num_values,
                                       ccs_numeric_t            **values);

static _ccs_distribution_ops_t _ccs_distribution_roulette_ops = {
	{ &_ccs_distribution_del },
	&_ccs_distribution_roulette_samples,
	&_ccs_distribution_roulette_get_bounds,
	&_ccs_distribution_roulette_strided_samples,
	&_ccs_distribution_roulette_soa_samples
};

static ccs_result_t
_ccs_distribution_roulette_get_bounds(_ccs_distribution_data_t *data,
                                      ccs_interval_t           *interval_ret) {
	_ccs_distribution_roulette_data_t *d = (_ccs_distribution_roulette_data_t *)data;

	interval_ret->type = CCS_NUM_INTEGER;
	interval_ret->lower = CCSI(INT64_C(0));
	interval_ret->upper = CCSI((ccs_int_t)(d->num_areas));
	interval_ret->lower_included = CCS_TRUE;
	interval_ret->upper_included = CCS_FALSE;
	return CCS_SUCCESS;
}

static ccs_result_t
_ccs_distribution_roulette_samples(_ccs_distribution_data_t *data,
                                   ccs_rng_t                 rng,
                                   size_t                    num_values,
                                   ccs_numeric_t            *values) {
	_ccs_distribution_roulette_data_t *d = (_ccs_distribution_roulette_data_t *)data;

	gsl_rng *grng;
	CCS_VALIDATE(ccs_rng_get_gsl_rng(rng, &grng));

	for (size_t i = 0; i < num_values; i++) {
		ccs_float_t rnd = gsl_rng_uniform(grng);
		ccs_int_t index = ccs_dichotomic_search(d->num_areas, d->areas, rnd);
		values[i].i = index;
	}
	return CCS_SUCCESS;
}

static ccs_result_t
_ccs_distribution_roulette_strided_samples(_ccs_distribution_data_t *data,
                                           ccs_rng_t                 rng,
                                           size_t                    num_values,
                                           size_t                    stride,
                                           ccs_numeric_t            *values) {
	_ccs_distribution_roulette_data_t *d = (_ccs_distribution_roulette_data_t *)data;

	gsl_rng *grng;
	CCS_VALIDATE(ccs_rng_get_gsl_rng(rng, &grng));

	for (size_t i = 0; i < num_values; i++) {
		ccs_float_t rnd = gsl_rng_uniform(grng);
		ccs_int_t index = ccs_dichotomic_search(d->num_areas, d->areas, rnd);
		values[i*stride].i = index;
	}
	return CCS_SUCCESS;
}

static ccs_result_t
_ccs_distribution_roulette_soa_samples(_ccs_distribution_data_t  *data,
                                     ccs_rng_t                  rng,
                                     size_t                     num_values,
                                     ccs_numeric_t            **values) {
	if (*values)
		return _ccs_distribution_roulette_samples(data, rng, num_values, *values);
	return CCS_SUCCESS;
}

ccs_result_t
ccs_create_roulette_distribution(size_t              num_areas,
                                 ccs_float_t        *areas,
                                 ccs_distribution_t *distribution_ret) {
	CCS_CHECK_ARY(num_areas, areas);
	CCS_CHECK_PTR(distribution_ret);
	if (!num_areas || num_areas > INT64_MAX)
		return -CCS_INVALID_VALUE;
	ccs_float_t sum = 0.0;

	for(size_t i = 0; i < num_areas; i++) {
		if (areas[i] < 0.0)
			return -CCS_INVALID_VALUE;
		sum += areas[i];
	}
	if (sum == 0.0)
		return -CCS_INVALID_VALUE;
	ccs_float_t inv_sum = 1.0/sum;
	if (isnan(inv_sum) || !isfinite(inv_sum))
		return -CCS_INVALID_VALUE;

	uintptr_t mem = (uintptr_t)calloc(1, sizeof(struct _ccs_distribution_s) + sizeof(_ccs_distribution_roulette_data_t) + sizeof(ccs_float_t)*(num_areas + 1) + sizeof(ccs_numeric_type_t));

	if (!mem)
		return -CCS_OUT_OF_MEMORY;

	ccs_distribution_t distrib = (ccs_distribution_t)mem;
	_ccs_object_init(&(distrib->obj), CCS_DISTRIBUTION, (_ccs_object_ops_t *)&_ccs_distribution_roulette_ops);
	_ccs_distribution_roulette_data_t * distrib_data = (_ccs_distribution_roulette_data_t *)(mem + sizeof(struct _ccs_distribution_s));
	distrib_data->common_data.data_types    = (ccs_numeric_type_t *)(mem + sizeof(struct _ccs_distribution_s) + sizeof(_ccs_distribution_roulette_data_t) + sizeof(ccs_float_t)*(num_areas + 1));
	distrib_data->common_data.type          = CCS_ROULETTE;
	distrib_data->common_data.dimension     = 1;
	distrib_data->common_data.data_types[0] = CCS_NUM_INTEGER;
	distrib_data->num_areas                 = num_areas;
	distrib_data->areas                     = (ccs_float_t *)(mem + sizeof(struct _ccs_distribution_s) + sizeof(_ccs_distribution_roulette_data_t));

	distrib_data->areas[0] = 0.0;
	for(size_t i = 1; i <= num_areas; i++) {
		distrib_data->areas[i] = distrib_data->areas[i-1] + areas[i-1] * inv_sum;
	}
	distrib_data->areas[num_areas] = 1.0;
	if (distrib_data->areas[num_areas] < distrib_data->areas[num_areas-1]) {
		free((void *)mem);
		return -CCS_INVALID_VALUE;
	}
	distrib->data = (_ccs_distribution_data_t *)distrib_data;
	*distribution_ret = distrib;
	return CCS_SUCCESS;
}

ccs_result_t
ccs_roulette_distribution_get_num_areas(ccs_distribution_t  distribution,
                                        size_t             *num_areas_ret) {
	CCS_CHECK_DISTRIBUTION(distribution, CCS_ROULETTE);
	CCS_CHECK_PTR(num_areas_ret);
	_ccs_distribution_roulette_data_t * data = (_ccs_distribution_roulette_data_t *)distribution->data;
	*num_areas_ret = data->num_areas;
	return CCS_SUCCESS;
}

ccs_result_t
ccs_roulette_distribution_get_areas(ccs_distribution_t  distribution,
                                    size_t              num_areas,
                                    ccs_float_t        *areas,
                                    size_t             *num_areas_ret) {
	CCS_CHECK_DISTRIBUTION(distribution, CCS_ROULETTE);
	CCS_CHECK_ARY(num_areas, areas);
	if (!areas && !num_areas_ret)
		return -CCS_INVALID_VALUE;
	_ccs_distribution_roulette_data_t * data = (_ccs_distribution_roulette_data_t *)distribution->data;
	if (areas) {
		if (num_areas < data->num_areas)
			return -CCS_INVALID_VALUE;
		for (size_t i = 0; i < data->num_areas; i++)
			areas[i] = data->areas[i+1] - data->areas[i];
		for (size_t i = data->num_areas; i < num_areas; i++)
			areas[i] = 0.0;
	}
	if (num_areas_ret)
		*num_areas_ret = data->num_areas;
	return CCS_SUCCESS;
}

max_stars_repo_name: deephyper/CCS

max_issues_repo_name: deephyper/CCS

max_forks_repo_name: deephyper/CCS


text:
/* wavelet/gsl_wavelet.h
 *
 * Copyright (C) 2004 Ivo Alxneit
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#ifndef __GSL_WAVELET_H__
#define __GSL_WAVELET_H__
#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS          /* empty */
# define __END_DECLS            /* empty */
#endif

__BEGIN_DECLS

#ifndef GSL_DISABLE_DEPRECATED
typedef enum {
  forward = 1, backward = -1,
  gsl_wavelet_forward = 1, gsl_wavelet_backward = -1
}
gsl_wavelet_direction;
#else
typedef enum {
  gsl_wavelet_forward = 1, gsl_wavelet_backward = -1
}
gsl_wavelet_direction;
#endif

typedef struct
{
  const char *name;
  int (*init) (const double **h1, const double **g1,
               const double **h2, const double **g2, size_t * nc,
               size_t * offset, size_t member);
}
gsl_wavelet_type;

typedef struct
{
  const gsl_wavelet_type *type;
  const double *h1;
  const double *g1;
  const double *h2;
  const double *g2;
  size_t nc;
  size_t offset;
}
gsl_wavelet;

typedef struct
{
  double *scratch;
  size_t n;
}
gsl_wavelet_workspace;

GSL_VAR const gsl_wavelet_type *gsl_wavelet_daubechies;
GSL_VAR const gsl_wavelet_type *gsl_wavelet_daubechies_centered;
GSL_VAR const gsl_wavelet_type *gsl_wavelet_haar;
GSL_VAR const gsl_wavelet_type *gsl_wavelet_haar_centered;
GSL_VAR const gsl_wavelet_type *gsl_wavelet_bspline;
GSL_VAR const gsl_wavelet_type *gsl_wavelet_bspline_centered;

gsl_wavelet *gsl_wavelet_alloc (const gsl_wavelet_type * T, size_t k);
void gsl_wavelet_free (gsl_wavelet * w);
const char *gsl_wavelet_name (const gsl_wavelet * w);

gsl_wavelet_workspace *gsl_wavelet_workspace_alloc (size_t n);
void gsl_wavelet_workspace_free (gsl_wavelet_workspace * work);

int gsl_wavelet_transform (const gsl_wavelet * w,
                           double *data, size_t stride, size_t n,
                           gsl_wavelet_direction dir,
                           gsl_wavelet_workspace * work);

int gsl_wavelet_transform_forward (const gsl_wavelet * w,
                                   double *data, size_t stride, size_t n,
                                   gsl_wavelet_workspace * work);

int gsl_wavelet_transform_inverse (const gsl_wavelet * w,
                                    double *data, size_t stride, size_t n,
                                    gsl_wavelet_workspace * work);

__END_DECLS

#endif /* __GSL_WAVELET_H__ */

ext: h

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
#include <bindings.cmacros.h>
#include <gsl/gsl_version.h>

BC_GLOBALARRAY(GSL_VERSION,char)
BC_GLOBALARRAY(gsl_version,char)

ext: c


text:
/* These functions compute linear preprocessing for
the UE using LAPACKE and CBLAS modules of
LAPACK libraries.
MMSE and MMSE whitening filters are available.
Functions are using RowMajor storage of the
matrices, like in conventional C. Traditional
Fortran functions of LAPACK employ ColumnMajor
data storage. */

#include<stdio.h>
#include<math.h>
#include<complex.h>
#include <stdlib.h>
#include <cblas.h>
#include <string.h>
#include <linux/version.h>
#if RHEL_RELEASE_CODE >= 1796
#include <lapacke/lapacke_utils.h>
#include <lapacke/lapacke.h>
#else
#include <lapacke_utils.h>
#include <lapacke.h>
#endif
//#define DEBUG_PREPROC

void transpose (int N, float complex *A, float complex *Result)
{
  // COnputes C := alpha*op(A)*op(B) + beta*C,
  enum CBLAS_TRANSPOSE transa = CblasTrans;
  enum CBLAS_TRANSPOSE transb = CblasNoTrans;
  int rows_opA = N; // number of rows in op(A) and in C
  int col_opB = N; //number of columns of op(B) and in C
  int col_opA = N; //number of columns in op(A) and rows in op(B)
  int col_B; //number of columns in B
  float complex alpha = 1.0+I*0;
  int lda  = rows_opA;
  float complex beta = 0.0+I*0;
  int ldc = rows_opA;
  int i;
  float complex* B;

  int ldb = col_opB;

  if (transb == CblasNoTrans) {
    B = (float complex*)calloc(ldb*col_opB,sizeof(float complex));
    col_B= col_opB;
  }
  else {
    B = (float complex*)calloc(ldb*col_opA, sizeof(float complex));
    col_B = col_opA;
  }
  float complex* C = (float complex*)malloc(ldc*col_opB*sizeof(float complex));

  for (i=0; i<lda*col_B; i+=N+1)
    B[i]=1.0+I*0;

  cblas_cgemm(CblasRowMajor, transa, transb, rows_opA, col_opB, col_opA, &alpha, A, lda, B, ldb, &beta, C, ldc);

  memcpy(Result, C, N*N*sizeof(float complex));

  free(B);
  free(C);
 }

void conjugate_transpose (int N, float complex *A, float complex *Result)
{
  // Computes C := alpha*op(A)*op(B) + beta*C,
  enum CBLAS_TRANSPOSE transa = CblasConjTrans;
  enum CBLAS_TRANSPOSE transb = CblasNoTrans;
  int rows_opA = N; // number of rows in op(A) and in C
  int col_opB = N; //number of columns of op(B) and in C
  int col_opA = N; //number of columns in op(A) and rows in op(B)
  int col_B; //number of columns in B
  float complex alpha = 1.0+I*0;
  int lda  = rows_opA;
  float complex beta = 0.0+I*0;
  int ldc = rows_opA;
  int i;
  float complex* B;
  int ldb = col_opB;

  if (transb == CblasNoTrans) {
    B = (float complex*)calloc(ldb*col_opB,sizeof(float complex));
    col_B= col_opB;
  }
  else {
    B = (float complex*)calloc(ldb*col_opA, sizeof(float complex));
    col_B = col_opA;
  }
  float complex* C = (float complex*)malloc(ldc*col_opB*sizeof(float complex));

  for (i=0; i<lda*col_B; i+=N+1)
    B[i]=1.0+I*0;

  cblas_cgemm(CblasRowMajor, transa, transb, rows_opA, col_opB, col_opA, &alpha, A, lda, B, ldb, &beta, C, ldc);

  memcpy(Result, C, N*N*sizeof(float complex));

  free(B);
  free(C);
 }

void H_hermH_plus_sigma2I (int N, int M, float complex *A, float sigma2, float complex *Result)
{
  //C := alpha*op(A)*op(B) + beta*C,
  enum CBLAS_TRANSPOSE transa = CblasConjTrans;
  enum CBLAS_TRANSPOSE transb = CblasNoTrans;
  int rows_opA = N; // number of rows in op(A) and in C
  int col_opB = N; //number of columns of op(B) and in C
  int col_opA = N; //number of columns in op(A) and rows in op(B)
  int col_C = N; //number of columns in B
  float complex alpha = 1.0+I*0;
  int lda  = col_opA;
  float complex beta = 1.0 + I*0;
  int ldc = col_opA;
  int i;

  float complex* C = (float complex*)calloc(ldc*col_opB, sizeof(float complex));

  for (i=0; i<lda*col_C; i+=N+1)
    C[i]=sigma2*(1.0+I*0);

  cblas_cgemm(CblasRowMajor, transa, transb, rows_opA, col_opB, col_opA, &alpha, A, lda, A, lda, &beta, C, ldc);

  memcpy(Result, C, N*M*sizeof(float complex));
  free(C);

 }

 void HH_herm_plus_sigma2I (int M, int N, float complex *A, float sigma2, float complex *Result)
{
  //C := alpha*op(A)*op(B) + beta*C,
  enum CBLAS_TRANSPOSE transa = CblasNoTrans;
  enum CBLAS_TRANSPOSE transb = CblasConjTrans;
  int k = N; //number of columns in op(A) and rows in op(B),k
  float complex alpha = 1.0+I*0;
  int lda  = N;
  int ldb  = N;
  int ldc = M;
  int i;

  float complex* C = (float complex*)calloc(M*M, sizeof(float complex));

  for (i=0; i<M*M; i+=M+1)
    C[i]=1.0+I*0;

  cblas_cgemm(CblasRowMajor, transa, transb, M, M, k, &alpha, A, lda, A, ldb, &sigma2, C, ldc);

  memcpy(Result, C, M*M*sizeof(float complex));
  free(C);

}

void eigen_vectors_values (int N, float complex *A, float complex *Vectors, float *Values_Matrix)
{
  // This function computes ORTHONORMAL eigenvectors and eigenvalues of matrix A,
  // where Values_Matrix is a diagonal matrix of eigenvalues.
  // A=Vectors*Values_Matrix*Vectors'
  char jobz = 'V';
  char uplo = 'U';
  int order_A = N;
  int lda = N;
  int i;
  float* Values = (float*)malloc(sizeof(float)*1*N);

  LAPACKE_cheev(LAPACK_ROW_MAJOR, jobz, uplo, order_A, A, lda, Values);

  memcpy(Vectors, A, N*N*sizeof(float complex));

  for (i=0; i<lda; i+=1)
    Values_Matrix[i*(lda+1)]=Values[i];

  free(Values);
}

 void lin_eq_solver (int N, float complex* A, float complex* B, float complex* Result)
{
  int n = N;
  int lda = N;
  int ldb = N;
  int nrhs = N;

  char transa = 'N';
  int* IPIV = malloc(N*N*sizeof(int));

  // Compute LU-factorization
  LAPACKE_cgetrf(LAPACK_ROW_MAJOR, n, nrhs, A, lda, IPIV);

  // Solve AX=B
  LAPACKE_cgetrs(LAPACK_ROW_MAJOR, transa, n, nrhs, A, lda, IPIV, B, ldb);

  // cgetrs( "N", N, 4, A, lda, IPIV, B, ldb, INFO )

  memcpy(Result, B, N*N*sizeof(float complex));

  free(IPIV);

}

void mutl_matrix_matrix_row_based(float complex* M0, float complex* M1, int rows_M0, int col_M0, int rows_M1, int col_M1, float complex* Result ){
  enum CBLAS_TRANSPOSE transa = CblasNoTrans;
  enum CBLAS_TRANSPOSE transb = CblasNoTrans;
  int rows_opA = rows_M0; // number of rows in op(A) and in C
  int col_opB = col_M1; //number of columns of op(B) and in C
  int col_opA = col_M0; //number of columns in op(A) and rows in op(B)
  float complex alpha =1.0;
  int lda  = col_M0;
  float complex beta = 0.0;
  int ldc = col_M1;
  int ldb = col_M1;

#ifdef DEBUG_PREPROC
  int i=0;
  printf("rows_M0 %d, col_M0 %d, rows_M1 %d, col_M1 %d\n", rows_M0, col_M0, rows_M1, col_M1);

  for(i=0; i<rows_M0*col_M0; ++i)
    printf(" rows_opA = %d, col_opB = %d, W_MMSE[%d] = (%f + i%f)\n", rows_opA, col_opB,  i , creal(M0[i]), cimag(M0[i]));

  for(i=0; i<rows_M1*col_M1; ++i)
    printf(" M1[%d] = (%f + i%f)\n", i , creal(M1[i]), cimag(M1[i]));
#endif

  cblas_cgemm(CblasRowMajor, transa, transb, rows_opA, col_opB, col_opA, &alpha, M0, lda, M1, ldb, &beta, Result, ldc);

#ifdef DEBUG_PREPROC
  for(i=0; i<rows_opA*col_opB; ++i)
    printf(" result[%d] = (%f + i%f)\n", i , creal(Result[i]), cimag(Result[i]));
#endif

}
void mutl_matrix_matrix_col_based(float complex* M0, float complex* M1, int rows_M0, int col_M0, int rows_M1, int col_M1, float complex* Result ){
  enum CBLAS_TRANSPOSE transa = CblasNoTrans;
  enum CBLAS_TRANSPOSE transb = CblasNoTrans;
  int rows_opA = rows_M0; // number of rows in op(A) and in C
  int col_opB = col_M1; //number of columns of op(B) and in C
  int col_opA = col_M0; //number of columns in op(A) and rows in op(B)
  float complex alpha =1.0;
  int lda  = col_M0;
  float complex beta = 0.0;
  int ldc = rows_M1;
  int ldb = rows_M1;

#ifdef DEBUG_PREPROC
  int i = 0;
  printf("rows_M0 %d, col_M0 %d, rows_M1 %d, col_M1 %d\n", rows_M0, col_M0, rows_M1, col_M1);

  for(i=0; i<rows_M0*col_M0; ++i)
    printf(" rows_opA = %d, col_opB = %d, W_MMSE[%d] = (%f + i%f)\n", rows_opA, col_opB,  i , creal(M0[i]), cimag(M0[i]));

  for(i=0; i<rows_M1*col_M1; ++i)
    printf(" M1[%d] = (%f + i%f)\n", i , creal(M1[i]), cimag(M1[i]));
#endif

  cblas_cgemm(CblasColMajor, transa, transb, rows_opA, col_opB, col_opA, &alpha, M0, lda, M1, ldb, &beta, Result, ldc);

#ifdef DEBUG_PREPROC
  for(i=0; i<rows_opA*col_opB; ++i)
    printf(" result[%d] = (%f + i%f)\n", i , creal(Result[i]), cimag(Result[i]));
#endif
}

/*FILTERS */
void compute_MMSE(float complex* H, int order_H, float sigma2, float complex* W_MMSE)
{
  int N = order_H;
  float complex* H_hermH_sigmaI = malloc(N*N*sizeof(float complex));
  float complex* H_herm =  malloc(N*N*sizeof(float complex));

  H_hermH_plus_sigma2I(N, N, H, sigma2, H_hermH_sigmaI);

#ifdef DEBUG_PREPROC
  int i =0;
  for(i=0;i<N*N;i++)
    printf(" H_hermH_sigmaI[%d] = (%f + i%f)\n", i , creal(H_hermH_sigmaI[i]), cimag(H_hermH_sigmaI[i]));
#endif

  conjugate_transpose (N, H, H_herm); //equals H_herm

#ifdef DEBUG_PREPROC
  for(i=0;i<N*N;i++)
    printf(" H_herm[%d] = (%f + i%f)\n", i , creal(H_herm[i]), cimag(H_herm[i]));
#endif

  lin_eq_solver(N, H_hermH_sigmaI, H_herm, W_MMSE);

#ifdef DEBUG_PREPROC
  for(i=0;i<N*N;i++)
    printf(" W_MMSE[%d] = (%f + i%f)\n", i , creal(W_MMSE[i]), cimag(W_MMSE[i]));
#endif

  free(H_hermH_sigmaI);
  free(H_herm);
}

#if 0
void compute_white_filter(float complex* H_re,
                          int order_H,
                          float sigma2,
                          float complex* W_Wh_0_re,
                          float complex* W_Wh_1_re){

  int aatx, aarx, re;
  int i,j;
  int M =n_rx;
  int N = n_tx;
  int sigma2=noise_power;

  float complex *H0_re = malloc(n_rx*(n_tx>>2)*sizeof(float complex));
  float complex *H1_re = malloc(n_rx*(n_tx>>2)*sizeof(float complex));
  float complex *R_corr_col_n_0_re = malloc(n_rx*n_tx*sizeof(float complex));
  float complex *R_corr_col_n_1_re = malloc(n_rx*n_tx*sizeof(float complex));
  float complex *U_0_re = malloc(n_rx*n_tx*sizeof(float complex));
  float complex *U_1_re = malloc(n_rx*n_tx*sizeof(float complex));
  float complex *U_0_herm_re = malloc(n_rx*n_tx*sizeof(float complex));
  float complex *U_1_herm_re = malloc(n_rx*n_tx*sizeof(float complex));
  float complex *D_0_re = malloc(n_rx*n_tx*sizeof(float complex));
  float complex *D_1_re = malloc(n_rx*n_tx*sizeof(float complex));
  float complex *W_Wh_0_re = malloc(n_rx*n_tx*sizeof(float complex));
  float complex *W_Wh_1_re = malloc(n_rx*n_tx*sizeof(float complex));

  for (aatx=0; aatx<n_tx/2; aatx++){
    for (aarx=0; aarx<n_rx; aarx++) {
      H0_re[aatx*n_rx + aarx] = H_re[aatx*n_rx + aarx][re]; // H0 gets [0 1 2 3; 4,5,6,7].' coefficients of H
      H1_re[aatx*n_rx + aarx] = H_re[aatx*n_rx + aarx + 8][re]; // H1 gets [8 9 10 11; 12, 13, 14, 15].' coefficients of H
        if (re == 0)
        printf("ant %d, H_re = (%f + i%f) \n", aatx*n_rx + aarx, creal(H[aatx*n_rx + aarx][re]), cimag(H[aatx*n_rx + aarx][re]));
      }
   }

    //HH_herm_plus_sigma2I(n_rx, (n_tx>>2), H1_re, sigma2, R_corr_col_n_0_re);
    HH_herm_plus_sigma2I(n_rx, (n_tx>>2), H0_re, sigma2, R_corr_col_n_1_re);

    eigen_vectors_values(n_rx, R_corr_col_n_0_re, U_0_re, D_0_re);
    eigen_vectors_values(n_rx, R_corr_col_n_1_re, U_1_re, D_1_re);

    transpose (n_rx, U_0_re, U_0_herm_re);
    transpose (n_rx, U_1_re, U_1_herm_re);

    sigma = (float)(sqrt((double)(sigma2)));

    /*The inverse of a diagonal matrix is obtained by replacing each element in the diagonal with its reciprocal.
    A square root of a diagonal matrix is given by the diagonal matrix, whose diagonal entries are just the square
     roots of the original matrix.*/

    D_0_re_inv_sqrt[0] = sqrt_float(1/D_0_re_inv[0]);
    D_0_re_inv_sqrt[5] = sqrt_float(1/D_0_re_inv[5]);
    D_0_re_inv_sqrt[10] = sqrt_float(1/D_0_re_inv[10]);
    D_0_re_inv_sqrt[15] = sqrt_float(1/D_0_re_inv[15]);

    D_1_re_inv[0] = sqrt_float(1/D_1_re_inv[0]);
    D_1_re_inv[5] = sqrt_float(1/D_1_re_inv[5]);
    D_1_re_inv[10] = sqrt_float(1/D_1_re_inv[10]);
    D_1_re_inv[15] = sqrt_float(1/D_1_re_inv[15]);

    now only to multiply

  free(H0);
  free(H1);
  free(R_corr_col_n_0);
  free(R_corr_col_n_1);
}
#endif

float sqrt_float(float x, float sqrt_x)
{
  sqrt_x = (float)(sqrt((double)(x)));
  return sqrt_x;
}

max_stars_repo_name: shadansari/onos-cu-cp

max_issues_repo_name: shadansari/onos-cu-cp

max_forks_repo_name: shadansari/onos-cu-cp


text:
/*
 * Copyright (c) 2011-2013  University of Texas at Austin. All rights reserved.
 *
 * $COPYRIGHT$
 *
 * Additional copyrights may follow
 *
 * This file is part of PerfExpert.
 *
 * PerfExpert is free software: you can redistribute it and/or modify it under
 * the terms of the The University of Texas at Austin Research License
 *
 * PerfExpert is distributed in the hope that it will be useful, but WITHOUT ANY
 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
 * A PARTICULAR PURPOSE.
 *
 * Authors: Leonardo Fialho and Ashay Rane
 *
 * $HEADER$
 */

#ifndef ANALYSIS_DEFS_H_
#define ANALYSIS_DEFS_H_

#include <vector>
#include <gsl/gsl_histogram.h>

#include "tools/macpo/common/avl_tree.h"
#include "tools/macpo/common/generic_defs.h"
#include "tools/macpo/common/macpo_record.h"
#include "tools/macpo/common/macpo_record_cxx.h"

typedef gsl_histogram histogram_t;

typedef std::vector<mem_info_list_t> mem_info_bucket_t;
typedef std::vector<trace_info_list_t> trace_info_bucket_t;
typedef std::vector<vector_stride_info_list_t> vector_stride_info_bucket_t;

typedef std::vector<histogram_t*> histogram_list_t;
typedef std::vector<histogram_list_t> histogram_matrix_t;

typedef std::vector<avl_tree*> avl_tree_list_t;

typedef std::pair<size_t, size_t> pair_t;
typedef std::vector<pair_t> pair_list_t;

typedef std::vector<double> double_list_t;
typedef std::vector<int> int_list_t;

typedef struct {
    size_t size, line_size, associativity, count;
} cache_data_t;

typedef struct {
    cache_data_t l1_data, l2_data, l3_data;
    name_list_t stream_list;
    mem_info_bucket_t mem_info_bucket;
    trace_info_bucket_t trace_info_bucket;
    vector_stride_info_bucket_t vector_stride_info_bucket;
} global_data_t;

typedef struct {
    size_t count, l1_conflicts, l2_conflicts;
} stream_list_t;

/* Different kinds of analyses options. */
#define ANALYSIS_REUSE_DISTANCE     (1 << 0)
#define ANALYSIS_CACHE_CONFLICTS    (1 << 1)
#define ANALYSIS_PREFETCH_STREAMS   (1 << 2)
#define ANALYSIS_STRIDES            (1 << 3)
#define ANALYSIS_VECTOR_STRIDES     (1 << 4)

#define ANALYSIS_ALL                (~0)

/* Number of reuse distance entries to display for each stream. */
#define REUSE_DISTANCE_COUNT        3

#define ADDR_TO_CACHE_LINE(x)   (x >> 6)

#endif /* ANALYSIS_DEFS_H_ */

max_stars_repo_name: roystgnr/perfexpert

max_issues_repo_name: roystgnr/perfexpert

max_forks_repo_name: roystgnr/perfexpert


text:
#include <linux/types.h>
#include "sweeny_uf.h"
#include "extract_args.h"
#include "timeseries.h"
#include <stdio.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_rng.h> //Verwendung von Mersenne-Twister Pseudozufallszahlengenerator

#define TS_FILE_D "../time_series/uf/simulation_l%u_q%.4f_b%.4f_c%.4f_s%u.hdf5"
#define BUF_LEN 512
char fn[BUF_LEN];
char impl_title[] = "Union-Find implementation";

char verbose=0;
__u64 *four_cs_moment ;
__u64 *sec_cs_moment ;
__u32 *size_giant;
__u32 *num_bonds;
__u32 *num_cluster;
__u32 DX;
__u32 N;
__u32 seed=123456;
double q=2;
double coupling;
double beta;
double v;
double K;
__u32 steps;
__u32 cutoff;
int main(int argc, char *argv[])
{

    if(!extractArgs(argc, argv,impl_title))
      return EXIT_FAILURE;
    snprintf(fn,BUF_LEN,TS_FILE_D,DX,q,beta,coupling,seed);
    if(!init_observables())
        return EXIT_FAILURE;
    if(!init_sweeny_uf( q,DX,beta,coupling,cutoff,steps,seed,num_bonds,num_cluster,size_giant,sec_cs_moment,four_cs_moment))
        return EXIT_FAILURE;
    if(!simulate_sweeny_uf())
        return EXIT_FAILURE;
    destroy_sweeny_uf();
    save_timeseries();
    destroy_observables();
    return EXIT_SUCCESS;
}
/******************************************************************************
 *****************************************************************************/

max_stars_repo_name: ernmeel/sweeny

max_issues_repo_name: ernmeel/sweeny


text:
/* roots/utility.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Reid Priedhorsky, Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

/* utility.c -- various root finding utility routines */

/* config headers */
#include <config.h>

/* standard headers */
#include <stddef.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <float.h>

/* gsl headers */
#include <gsl/gsl_errno.h>
#include <gsl/gsl_roots.h>

/* roots headers */
#include "roots.h"

/* Validate arguments common to gsl_root_bisection and gsl_root_falsepos.
   Return GSL_SUCCESS if all arguments are okay, complain appropriately (i.e.
   call GSL_ERROR and return GSL_FAILURE) otherwise. */
int
_gsl_root_validate_bfp_args (double *root, double (*f) (double),
                             double *lower_bound,
                             double *upper_bound, double epsrel,
                             double epsabs, unsigned int max_iterations,
                             double max_deltay)
{
  /* Is the maximum delta-y too small? */
  if (max_deltay < GSL_ROOT_MIN_MAX_DELTAY)
    GSL_ERROR ("maximum delta-y negative, zero, or too small", GSL_EBADTOL);

  /* Did the user give a lower bound that not less than the upper bound? */
  if (*lower_bound >= *upper_bound)
    GSL_ERROR ("lower bound larger than upper_bound", GSL_EINVAL);

  /* The rest of the arguments are common. */
  return _gsl_root_validate_args (root, f, lower_bound, upper_bound,
                                  epsrel, epsabs, max_iterations);
}

/* Validate arguments commond to gsl_root_secant and gsl_root_newtons. Return
   GSL_SUCCESS if all arguments are okay, complain appropriately (i.e. call
   GSL_ERROR and return GSL_FAILURE) otherwise. */
int
_gsl_root_validate_sn_args (double *root, double (*f) (double),
                            double *guess1,
                            double *guess2, double epsrel,
                            double epsabs, unsigned int max_iterations,
                            double max_step_size)
{
  /* Is the maximum step size ridiculous? */
  if (max_step_size <= 0.0)
    GSL_ERROR ("maximum step size <= 0", GSL_EBADTOL);

  /* The rest of the arguments are common. */
  return _gsl_root_validate_args (root, f, guess1, guess2, epsrel,
                                  epsabs, max_iterations);
}

/* Validate the arguments common to all four low level functions. Return
   GSL_SUCCESS if all of the following arguments hold true, call GSL_ERROR and
   return GSL_FAILURE otherwise.

   * No pointer arguments are null.
   * The maximum number of iterations is non-zero.
   * Relative and absolute error are non-negative.
   * The relative error is not too small. */
int
_gsl_root_validate_args (double *root, double (*f) (double), double *where1,
                         double *where2, double epsrel,
                         double epsabs, unsigned int max_iterations)
{
  /* Are any pointers null? */
  if ((root == NULL) || (f == NULL) || (where1 == NULL)
      || (where2 == NULL))
    GSL_ERROR ("pointer argument null", GSL_EINVAL);
  /* Did the user tell us to do no iterations? */
  if (max_iterations == 0)
    GSL_ERROR ("maximum iterations 0", GSL_EINVAL);
  /* Did the user try to pawn a negative tolerance off on us? */
  if (epsrel < 0.0 || epsabs < 0.0)
    GSL_ERROR ("relative or absolute tolerance negative", GSL_EBADTOL);
  /* Is the relative error too small? */
  if (epsrel < GSL_DBL_EPSILON * GSL_ROOT_EPSILON_BUFFER)
    GSL_ERROR ("relative tolerance too small", GSL_EBADTOL);

  /* All is well. */
  return GSL_SUCCESS;
}

/* Verify that the supplied interval is guaranteed by the Intermediate Value
   Theorem to contain a root and complain appropriately if it is not. (It
   might actually be a discontinuity, but we check for that elsewhere.) Return
   GSL_SUCCESS if all is well, otherwise, call GSL_ERROR and return
   GSL_FAILURE. */
int
_gsl_root_ivt_guar (double (*f) (double), double lower_bound,
                    double upper_bound)
{
  double fl, fu;

  _BARF_FPCALL (f, lower_bound, fl);
  _BARF_FPCALL (f, upper_bound, fu);

  if (fl * fu > 0.0)
    {
      GSL_ERROR ("interval not guaranteed to contain a root", GSL_EINVAL);
    }
  else
    {
      return GSL_SUCCESS;
    }
}

/* Check if the user has the root but doesn't know it. If lower_bound or
   upper_bound is a root of f, or the interval [upper_bound, lower_bound] is
   within tolerance, return 1 and set *root appropriately. Otherwise, return
   0. On error, call GSL_ERROR and return GSL_FAILURE. Only worry about
   max_deltay if it is greater than 0 (which implies that you should validate
   arguments _before_ calling this function). */
int
_gsl_root_silly_user (double *root, double (*f) (double), double lower_bound,
                      double upper_bound, double epsrel,
                      double epsabs, double max_deltay)
{
  double fl, fu;

  /* Is lower_bound the root? */
  _BARF_FPCALL (f, lower_bound, fl);
  if (fl == 0.0)
    {
      *root = lower_bound;
      return 1;
    }

  /* Is upper_bound the root? */
  _BARF_FPCALL (f, upper_bound, fu);
  if (fu == 0.0)
    {
      *root = upper_bound;
      return 1;
    }

  /* Are lower_bound and upper_bound within tolerance? */
  _BARF_TOLS (lower_bound, upper_bound, 2 * epsrel, 2 * epsabs);
  if (max_deltay > 0.0)
    _BARF_DELTAY (fl, fu, max_deltay);
  if (_WITHIN_TOL (lower_bound, upper_bound, 2 * epsrel,
                   2 * epsabs))
    {
      *root = (lower_bound + upper_bound) / 2.0;
      return 1;
    }

  /* No? Bummer. */
  return 0;
}

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
/**
 * \author Sylvain Marsat, University of Maryland - NASA GSFC
 *
 * \brief C header for functions windowing and computing FFT/IFFT of time/frequency series.
 *
 */

#ifndef __GENERATELLVFD_H__
#define __GENERATELLVFD_H__ 1

#ifdef __GNUC__
#define UNUSED __attribute__ ((unused))
#else
#define UNUSED
#endif

#define _XOPEN_SOURCE 500

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>
#include <time.h>
#include <unistd.h>
#include <getopt.h>
#include <stdbool.h>
#include <string.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_complex.h>

#include "constants.h"
#include "struct.h"
#include "waveform.h"
#include "fft.h"
#include "LLVgeometry.h"
#include "LLVFDresponse.h"

/********************************** Structures ******************************************/

typedef enum GenLLVFDtag {
  LLVFD,
  LLVhlm
} GenLLVFDtag;

/* Parameters for the generation of a LLV waveform (in the form of a list of modes) */
typedef struct tagGenLLVFDparams {
  double tRef;               /* reference time (s) - GPS time at the frequency representing coalescence */
  double phiRef;             /* reference phase (rad) - phase at the frequency representing coalescence (or at fRef if specified) */
  double fRef;               /* reference frequency at which phiRef is set (Hz, default 0 which is interpreted as Mf=0.14) */
  double m1;                 /* mass of companion 1 (solar masses, default 2e6) */
  double m2;                 /* mass of companion 2 (solar masses, default 1e6) */
  double distance;           /* distance of source (Mpc, default 1e3) */
  double inclination;        /* inclination of source (rad, default pi/3) */
  double ra;                 /* first angle for the position in the sky (rad, default 0) */
  double dec;                /* second angle for the position in the sky (rad, default 0) */
  double polarization;       /* polarization angle (rad, default 0) */

  int nbmode;                /* number of modes to generate (starting with 22) - defaults to 5 (all modes) */
  double minf;               /* Minimal frequency, ignore if 0 (Hz, default=0) - will use first frequency covered by the ROM if higher */
  double maxf;               /* Maximal frequency, ignore if 0 (Hz, default=0) - will use last frequency covered by the ROM if lower */
  double deltaf;             /* When generating frequency series from the mode contributions, deltaf for the output (0 to set automatically at 1/2*1/(2T)) */
  int tagnetwork;            /* Tag selecting the desired output format */
  int taggenwave;            /* Tag selecting the desired output format */
  int fromLLVtdfile;         /* Tag for loading time series for LLV detectors and FFTing */
  int nsamplesinfile;        /* Number of lines of input file */
  int binaryin;              /* Tag for loading the data in gsl binary form instead of text (default false) */
  int binaryout;             /* Tag for outputting the data in gsl binary form instead of text (default false) */
  char indir[256];           /* Path for the input directory */
  char infile[256];          /* Path for the input file */
  char outdir[256];          /* Path for the output directory */
  char outfileprefix[256];   /* Path for the output file */
} GenLLVFDparams;

#if 0
{ /* so that editors will match succeeding brace */
#elif defined(__cplusplus)
}
#endif

#endif /* _GENERATELLVFD_H */

max_stars_repo_name: titodalcanton/flare

max_issues_repo_name: titodalcanton/flare

max_forks_repo_name: titodalcanton/flare


text:
#ifndef BN
#define BN

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <stdbool.h>
#include <getopt.h>
#include <assert.h>

#include <gsl/gsl_randist.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_matrix.h>

#include <omp.h>
#include <math.h>

#include "matrix.h"
#include "model.h"
#include "data.h"
#include "vector.h"
#include "a-queue.h"

bool verbose = false;
bool debug = false;

int** GENOTYPE; // GENOTYPE[i] stores the binary expansion of i.
int NUM_SAMPLES_READ; // total number of samples read in from dataset.

// set initial parameter values (seeimingly arbitrary, but makes it easier to tweak them since they're all here)
void initialize_MH_params(unsigned int* seed, int* number_samples, int* record_ith, int* burn_in, char* reg_scheme) {
	*seed = time(NULL); // random seed using time
	*number_samples = 100000;
	*record_ith = 1;
	*burn_in = 0;
	strcpy(reg_scheme, "bic");
}

/*
 * Check file inputs to bn. If no .poset file supplied (f flag), prints error message and exits.
 * Flags:
 *
 *
 *
 */
void check_file_inputs(int argc, char** argv, char* filestem, char* output, int* number_samples, unsigned int* seed, char* reg_scheme) {
	int c = 0;
	bool input_flag = false;
	bool output_flag = false;
	bool help_flag = false;
	static struct option long_options[] =
        {
          {"dataset",     		required_argument, 0, 'd'},
          {"output",  				required_argument, 0, 'o'},
          {"number_samples",	required_argument, 0, 'n'},
          {"seed", 						required_argument, 0, 's'},
          {"reg", 						required_argument, 0, 'r'},
          {"help", 						no_argument, 			 0, 'h'},
          {NULL, 0, NULL,0,}
        };
	while ((c = getopt_long(argc, argv, "d:o:n:s:r:hvt:e:", long_options, NULL)) != -1 ) {
		switch(c) {
			case 'd': // get input filename, satisfy boolean check for provided input.
				strcpy(filestem, optarg);
				input_flag = true;
				break;
			case 'o':
				strcpy(output, optarg);
				output_flag = true;
				break;
			case 'n':
				*number_samples = atoi(optarg);
				break;
			case 's':
				*seed = atoi(optarg);
				break;
			case 'r':
				strcpy(reg_scheme, optarg);
				break;
			case 'h':
				help_flag = true;
				break;
			case 'v':
				verbose = true;
				break;
		}
	}

	if (help_flag) {
		printf("Usage: ./cbn (-(option) (argument) )*\n");
		printf("where dataset input option (-d or --dataset) is required.\n");
		printf("Options include:\n");
		printf("-d | --dataset \tPath to input file with sample data. Please include file extension as well. ex: test1.txt\n");
		printf("-o | --output \tPath to output file to write results to. If not specified, then results not written to file.\n");
		printf("-n | --number_samples \tSet the number of iterations for Metropolis-Hastings to be run. Default value set to 500,000.\n");
		printf("-s | --seed \tSet the seed for random number generators in the program. Time used as default seed.\n");
		printf("-r | --reg \tSet the likelihood regularization scheme for the MCMC algorithm. Valid choices are bic, aic, and loglik.\n");
		printf("-h | --help \tPrint help information for program.\n");
		exit(1);
	}

	if (!input_flag) {
		fprintf(stderr, "Error: No input file specified. Please provide the filename (without filetype) of a file containing formatted sample data using the \"-d\" or \"--dataset\" options.\n");
		exit(1);
	}
	if (!output_flag) {
		strcpy(output, "");
	}
	if (strcmp(reg_scheme, "bic") &&  strcmp(reg_scheme, "aic") && strcmp(reg_scheme, "loglik")) {
		fprintf(stderr, "Error: Invalid regularization scheme. Please choose from bic (default), aic, or loglik.\n");
		exit(1);
	}
}

/*
 * Standard open_file helper. Throws error if error while opening file.
 */
FILE* open_file(char* filename, char* option) {
	FILE* input = fopen(filename, option);
	if (input == NULL) {
		fprintf(stderr, "Error: Could not read %s\n", filename);
		exit(1);
	}
	return input;
}

/*
 * Base 2 power function. Computes 2^n, for n >= 0
 */
static inline int pow2(int n) {
	assert(n >= 0);
	return 1 << n;
}

// Basic power function implementation. Only accepts integral power values.
double power(double base, int power) {
	double result = 1.0;
	for (int i = 0; i < power; i++) {
		result *= base;
	}
	return result;
}

/*
 * Generates a binary compression of pattern data, which is
 * used as a means of quickly determining whether two samples have
 * the same pattern. Presumably, this is faster (as opposed to
 * comparisons without such a compression scheme), since we compress
 * all samples once, and then compare the "hash" values.
 * Parallels index_of in TS code.
 */
int binary_compression(int* pat_i, int num_events) {
	int index = 0;
	for (int i = 0; i < num_events; i++) {
		if (pat_i[i] ==1)
			index += pow2(num_events - 1 - i);
	}
	return index;
}

// Used to generate binary expansions of genotypes, for easier access later on.
void precompute_binary(const int n) {
	int m = pow2(n);
	GENOTYPE = get_int_matrix(m, n);
	for (int i = 0; i < m; i++) {
		int count = 0;
		for (int j = 0; j < n; j++) {
			if ( ((1 << j) & i) > 0) {
				GENOTYPE[i][n - 1 - j] = 1;
				count++;
			} else {
				GENOTYPE[i][n - 1 - j] = 0;
			}
		}
	}
}

int count_edges(int** poset, int n) {
	int count = 0;
	for (int i=0; i< n; i++) {
		for (int j = 0; j< n; j++) {
			if (poset[i][j])
				count++;
		}
	}
	return count;
}

// TODO: reconsider the use of int* ptr vector here. We've changed the idea, and we're no longer looking
// at the probability of combined events. Instead, using the current poset to determine how children
// are valid or not in find_children_extended.

/*
 * Helper function for process_patterns function. Reads pattern data
 * from .txt file. Only 0/1 expected, but places negative values where
 * any deviations occur. Exits if any error is encountered.
 */
int** read_patterns(FILE *input, model* M, char* filename, int* num_samples, int* num_events) {
	size_t temp = 0;
	int n;
	char* line = NULL;
	int len = getline(&line, &temp, input);
	if (len != -1) {
		n = len / 2;
		free(line);
		line = NULL;
		rewind(input);
	} else {
		fprintf(stderr, "Error reading sample data.\n");
		exit(1);
	}

	M->n = n;
	*num_events = n;

	ptr_vector* v = new_ptr_vector(); // vector takes ints, hopefully won't get mutilated?
	while (getline(&line, &temp, input) != -1) {
		int line_index = 0;
		int* genotype = get_int_array(n);
		for (int i = 0; i< n; i++) { // fill in each event with genotype.
			if (line[line_index] == '0') genotype[i] = 0;
			if (line[line_index] == '1') genotype[i] = 1;
			line_index += 2;
		}
		push_back_ptr(v, genotype);
		free(line);
		line = NULL;
		temp = 0;
	}
	free(line);
	*num_samples = v->size;

	if (debug) {
		printf("Finished reading patterns.\nPatterns matrix: \n");
		print_int_matrix(v->elems, v->size, n);
	}

	return free_ptr_vector(v, false);
}

// counts uniq patterns among samples.
// It seems that pat_idx counts the number of unique samples
// up to that last unique sample.
// ----- variables -----
// index stores the "hashed" patterns
// count stores the number of samples with the same pattern
// Removed pat_idx necessity, since the same information is stored in index.

// Note: The binary compression places a limit on the number of loci that can be observed (32 for int, 64 for long int).
// For a completely scalable implementation, should re-work this scheme, use a hashset for speed.
void count_uniq_patterns(int* index, int* count, int num_samples, int num_events, int* N_u, int** pat) {

	//binary compression that "hashes" genotypes
	for (int i = 0; i< num_samples; i++)
		index[i] = binary_compression(pat[i], num_events);

	// count unique patterns. Runs in worst case O(n^2) time. can be optimized to O(nlgn) using sort -> iterate, or linear time using a hashtable/set.
	bool skips; // true if another sample has same pattern ("index" value)
	int num_unique = 0; // counts number of unique patterns
	for (int i = 0; i < num_samples; i++) {
		skips = false;
		for (int j = 0; j < i; j++) {
			if (index[i] == index[j]) { //samples have the same pattern
				count[j]++; //increment the count at first index
				skips = true;
				break;
			}
		}
		if (!skips) { // unique pattern --> NOT previously observed
			count[i]++;
			num_unique++;
		}
	}
	*N_u = num_unique;

	if (verbose) {
		printf("Number of samples processed = %d\n", num_samples);
		printf("Number of gene loci considered = %d\n", num_events);
		printf("Number of unique genotypes observed = %d\n", num_unique);
	}
	if (debug) {
		printf("num_samples = %d\n", num_samples);
		printf("count array: \n");
		print_int_array(count, num_samples);
		printf("index array: \n");
		print_int_array(index, num_samples);
	}
}

data* construct_data_set(int num_unique, int* count, int* index, int num_samples, int num_events, int** pat) {
	data* D = calloc(num_unique, sizeof(data));
	if (D == NULL) {
		fprintf(stderr, "Error: calloc failure - out of memory.\n");
		exit(1);
	}
	int d_index = 0;
	for (int i = 0; i < num_samples; i++) {
		if (count[i] > 0) {//if was not a repeat
			D[d_index].count = count[i];
			D[d_index].g = get_int_array(num_events);
			for (int j = 0; j < num_events; j++)
				D[d_index].g[j] = pat[i][j];
			d_index++;
		}
	}

	if (debug) {
		print_data_array(D, num_unique, num_events);
	}
	return D;
}

// N_u = number of unique genotypes in dataset
// n  = number of events
int process_patterns(char* filename, model* M, int* N_u, data** D) {

	FILE* input = open_file(filename, "r");
	int num_samples, num_events;
	int** pat = read_patterns(input, M, filename, &num_samples, &num_events);
	fclose(input);

	//Count unique patterns, construct data set.
	int* index = get_int_array(num_samples);
	int* count = get_int_array(num_samples);
	count_uniq_patterns(index, count, num_samples, num_events, N_u, pat);
	*D = construct_data_set(*N_u, count, index, num_samples, num_events, pat);

	NUM_SAMPLES_READ = num_samples;
	free_int_matrix(pat, num_samples);
	free_int_array(index);
	free_int_array(count);

	return 0;
}

// Suppes condition
// Confirm that the proposed modification (move) is supported by event probabilities in the data set, taking into account
// current theta types.
bool cumulative_probs_maintained(model* M, data* D, int* theta_types, int num_unique) {
	int n = M->n;
	for (int i = 0; i< n; i++) { // loop through each event
		int num_parents = M->num_pa[i];
		int* parents = M->pa[i];
		int type = theta_types[i];
		int count_parents = 0;
		int count_child = 0;
		int count_parents_child = 0;
		for (int j = 0; j< num_unique; j++) { // loop through all observed genotypes
			int* g = D[j].g;
			int found = 0;

			// get child present bool value
			bool child_present = false;
			if (g[i]) child_present = true;

			// get parent set satisfied bool value
			bool parents_present = true;
			switch (type) { // 0 or 1 parents, conjunctive
				case -1 ... 0:
					for (int k = 0; k < num_parents; k++) {
						if (!g[parents[k]]) parents_present = false;
					}
					break;
				case 1: // or
					for (int k = 0; k < num_parents; k++) {
						if (g[parents[k]]) found++;
					}
					if (!found) parents_present = false;
					break;
				case 2: // xor
					for (int k = 0; k < num_parents; k++) {
						if (g[parents[k]]) found++;
					}
					if (found != 1) parents_present = false;
					break;
			}

			if (child_present && parents_present) {
				count_parents_child += D[j].count;
			}

			if (child_present) {
				count_child += D[j].count;
			}

			if (parents_present) {
				count_parents += D[j].count;
			}
		}
		// check that the checks and additions are correct.
		if (num_parents > 0) {
			if (count_parents <= count_child || count_parents_child * NUM_SAMPLES_READ <= count_parents * count_child) {
				return false;
			}
		}
	}
	return true;
}

// Initializes the parent fields in the model so that they can be used to find children more efficiently.
void set_theta_types(model* M, int* theta_types, int* theta_types_p) {
	int num_events = M->n;
	int** poset = M->P;
	free_parents(M);
	M->num_pa = get_int_array(num_events); // create an array that stores the sizes of parents.
	M->pa = (int**) malloc(num_events * sizeof(int*)); // create an array to hold arrays of integers

	// for all events, create a vector of all other events that relate to them, according to the current poset.
	// Then, determine theta type if an event has > 1 parent events.
	for (int i = 0; i< num_events; i++) {
		vector* v = new_vector();
		for (int j = 0; j< num_events; j++) {
			if (poset[j][i]) {
				push_back(v, j);
			}
		}
		int num_parents = v->size;
		M->num_pa[i] = num_parents;
		M->pa[i] = free_vector(v, false);

		if (num_parents <= 1) { // if 1 parent at most, no need for special care.
			theta_types_p[i] = -1;
		} else {
			if (theta_types[i] < 0) { // event goes from 0/1 parents to 2+ parents
				theta_types_p[i] = rand() % 3; // randomly pick between 0, 1, and 2.
			} else { // carry over formula from previous iteration.
				theta_types_p[i] = theta_types[i];
			}
		}
	}
}

// imposes a total ordering on the events. Use Fisher-Yates shuffling algorithm to randomize ordering
void new_total_ordering(model* M, int* theta_types, data* D, int num_unique) {
	int num_events = M->n;
	// initialize new array to use.
	int rand_arr[num_events];
	for (int i = 0; i < num_events; i++) {
		rand_arr[i] = i;
	}

	// establish total ordering
	for (int i = num_events - 1; i >= 0; i--) {
		int randy = rand() % (i + 1);
		int temp = rand_arr[i];
		rand_arr[i] = rand_arr[randy];
		rand_arr[randy] = temp;
	}

	int** initial_poset = get_int_matrix(num_events, num_events);
	// mark in poset.
	for (int i = 0; i < num_events - 1; i++) {
		initial_poset[rand_arr[i]][rand_arr[i+1]] = 1;
	}

	// transitive closure & theta types
	int** T = get_int_matrix(num_events, num_events);
	transitive_closure(initial_poset, T, num_events);

	for (int i = 0; i < num_events; i++) {
		int num_parents = 0;
		for (int j = 0; j < num_events; j++) {
			if (T[j][i]) num_parents++;
		}
		if (num_parents == 0) continue; // don't bother checking if there are no parents, since we're only going to remove excess edges.
		int reduce = num_parents - (rand() % (num_parents + 1));
		for (int j = 0; j < reduce; j++) { // delete a random relation reduce number of times.
			int ith_event = rand() % num_parents;
			int counted = 0;
			for (int k = 0; k < num_events; k++) {
				if (counted == ith_event && T[k][i]) {
					T[k][i] = 0; // delete event
					num_parents--; // decrement number of parents.
					break;
				}
				if (T[k][i]) counted++;
			}
		}
	}

	free_int_matrix(initial_poset, num_events);
	free_int_matrix(M->P, num_events);
	M->P = T;

	int* theta_types_p = get_int_array(num_events);
	for (int i = 0; i < num_events; i++) {
		theta_types_p[i] = -1;
	}
	set_theta_types(M, theta_types_p, theta_types);

	model M_p;
	M_p.n = num_events;
	M_p.pa = M->pa;
	M_p.P = T;
	M_p.num_pa = get_int_array(num_events); // initializes each to 0.
	// do suppes checking for each events' parent set.
	for (int i = 0; i < num_events; i++) {
		// determine suppes prob condition for this event and its parent set.
		// copy the parent set of this event into a new poset.
		theta_types_p[i] = theta_types[i];
		M_p.num_pa[i] = M->num_pa[i];
		if (!cumulative_probs_maintained(&M_p, D, theta_types, num_unique)) {// reject entire parent set.
			for (int j = 0; j < num_events; j++) {
				M->P[j][i] = 0;
			}
			M->num_pa[i] = 0;
			theta_types[i] = -1;
		}
		theta_types_p[i] = -1; // push to conjunctive so that it doesn't go to 'or' and 'xor' cases.
		M_p.num_pa[i] = 0; // set back to 0.
	}
	free_int_array(M_p.num_pa);
	free_int_array(theta_types_p);
}

// returns the total number of events associated with a given genotype.
// ex. For n= 4, 0101 returns 2.
int total_events(const int genotype, const int n) {
	int* binary = GENOTYPE[genotype];
	int num_events = 0;
	for (int i=0; i< n; i++)
		num_events += binary[i];
	return num_events;
}

// Given two integers, returns the number of bitwise differences between the two.
// If hamming distance > 1, returns index of last (largest) bit.
int hamming_distance(unsigned int x, unsigned int y, const int num_events, int* diff_index) {
	unsigned int dist = 0;
	unsigned int diff = x ^ y;
	for (int i= 0; i < num_events; i++) {
		if ( (1 << i) & diff ) {
			dist++;
			if (diff_index != NULL)	*diff_index = num_events - 1 - i;
		}
	}
	return dist;
}

// Checks that the parent conditions of the ith event are satisfied in provided genotype g.
// Returns true if parent conditions (according to poset and theta_types) are fulfilled, false otherwise.
bool check_parent_set(model* M, int* theta_types, int* g, int i) {
	int num_parents = M->num_pa[i];
	int type = theta_types[i]; // relationship between parents.
	bool compatible = false;
	int found = 0;

	switch (type) { // will only go into one case for each, so this may be faster.
		case -1 ... 0: // 0/1 parents, or and case if > 1 parents. loop through all parents, confirm that all are present.
			for (int j = 0; j< num_parents; j++) {
				if (!g[M->pa[i][j]]) return false;
			}
			break;
		case 1: // or.
			for (int j = 0; j< num_parents; j++) {
				if (g[M->pa[i][j]]) {
					compatible = true;
					break;
				}
			}
			if (!compatible) return false;
			break;
		case 2: // xor. loop through all parents, confirm that only one present at a time.
			for (int j = 0; j< num_parents; j++) {
				if (g[M->pa[i][j]]) { // another event observed that isn't the current
					found++;
				}
			}
			if (found != 1) return false;
			break;
	}
	return true;
}

// Theta types gives the relationship between "parent events", if an event has more than one direct parent.
// If check_parent_set returns true for all events, then returns true.
bool compatible_with_poset(int genotype, model* M, int* theta_types) {
	int num_events = M->n;
	int* g = GENOTYPE[genotype];
	for (int i = 0; i< num_events; i++) { // for each event
		if (g[i]) { // if event is observed
			if (!check_parent_set(M, theta_types, g, i)) return false;
		}
	}
	return true;
}

#endif

max_stars_repo_name: BIMIB-DISCo/PMCE

max_issues_repo_name: BIMIB-DISCo/PMCE

max_forks_repo_name: BIMIB-DISCo/PMCE


text:
#ifndef RAYS_H
#define RAYS_H

#include <complex.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
#include <malloc.h>
#include <stdlib.h>
#include <string.h>
#include <sys/time.h>
#include <time.h>

#include "tsvread.h"

// A ray.
typedef struct {
  double p0[3];           // position
  double e0[3];           // direction
  complex double pol[3];  // polarization
  double OPL;             // optical path length
  double lambda;          // wavelength
  double intensity;       // intensity of this paricular ray
  double g;               // apodisation factor
} ray_t;

ray_t *readrays(char *filename, unsigned int *nrays);
void movefocus(ray_t *ray, unsigned int nrays);

#endif  // RAYS_H

max_stars_repo_name: AaronWebster/directdebye

max_issues_repo_name: AaronWebster/directdebye

max_forks_repo_name: AaronWebster/directdebye


text:
/*
  Template for applications of single array GSL functions
*/

#include <stdio.h>

#ifdef SINGLE_ARG
# include <stdlib.h>
#endif

#include <gsl/gsl_statistics_double.h>

#ifdef SORT_DATA
# include <gsl/gsl_sort.h>
#endif

#include "y.tab.h"
#include "double_matrix.h"

#define SINGLE_STEP 1

Double_Matrix *data = NULL;

int
main(int argc, char *argv[]) {
	size_t i = 0;
#ifdef SINGLE_ARG
	char *endptr = NULL;
	double arg;
#endif

#ifdef SINGLE_ARG
	if (2 != argc) {
		fprintf(stderr, "Argument required.\n");
		return 1;
	}

	arg = strtod(argv[1], &endptr);

	if (0 == arg && endptr == argv[ 1 ]) {
		fprintf(stderr, "Argument (%s) is not a valid number\n", argv[1]);
		return 1;
	}
#endif

	data = double_matrix_new();

	yyparse();

	for (i = 0; i < double_matrix_width(data); i++) {
#ifdef SORT_DATA
		gsl_sort(double_matrix_col_mut_data(data, i), SINGLE_STEP, double_matrix_depth(data));
#endif

		if (0 != i)
			printf( " " );

		printf("%.10g",
			FUNC_NAME(double_matrix_col_data(data, i), SINGLE_STEP, double_matrix_depth(data)
#ifdef SINGLE_ARG
		    , arg
#endif
			)
		);
	}

	printf("\n");

	double_matrix_free(data);

	return 0;
}


text:
/***************************************************************************
                          integration.h  -  description
                             -------------------
    copyright            : (C) 2005 by MOUSSA
    email                : mmousa@liris.cnrs.fr
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/

#ifndef INTEGRATIONFILE
#define INTEGRATIONFILE

#include "config.h"
#include "misc.h"

#include <gsl/gsl_math.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_plain.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>

void Integrate(int , int , double , Tetrahedron_3 , gsl_complex *, gsl_complex *);

#endif

max_stars_repo_name: kyeonghopark/jgt-code

max_issues_repo_name: kyeonghopark/jgt-code

max_forks_repo_name: kyeonghopark/jgt-code


text:
/* cdf/fdistinv.c
 *
 * Copyright (C) 2002 Przemyslaw Sliwa and Jason H. Stover.
 * Copyright (C) 2006 Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307, USA.
 */

#include <config.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>

#include "error.h"

double
gsl_cdf_fdist_Pinv (const double P, const double nu1, const double nu2)
{
  double result;
  double y;

  if (P < 0.0)
    {
      CDF_ERROR ("P < 0.0", GSL_EDOM);
    }
  if (P > 1.0)
    {
      CDF_ERROR ("P > 1.0", GSL_EDOM);
    }
  if (nu1 < 1.0)
    {
      CDF_ERROR ("nu1 < 1", GSL_EDOM);
    }
  if (nu2 < 1.0)
    {
      CDF_ERROR ("nu2 < 1", GSL_EDOM);
    }

  if (P < 0.5)
    {
      y = gsl_cdf_beta_Pinv (P, nu1 / 2.0, nu2 / 2.0);

      result = nu2 * y / (nu1 * (1.0 - y));
    }
  else
    {
      y = gsl_cdf_beta_Qinv (P, nu2 / 2.0, nu1 / 2.0);

      result = nu2 * (1 - y) / (nu1 * y);
    }

  return result;
}

double
gsl_cdf_fdist_Qinv (const double Q, const double nu1, const double nu2)
{
  double result;
  double y;

  if (Q < 0.0)
    {
      CDF_ERROR ("Q < 0.0", GSL_EDOM);
    }
  if (Q > 1.0)
    {
      CDF_ERROR ("Q > 1.0", GSL_EDOM);
    }
  if (nu1 < 1.0)
    {
      CDF_ERROR ("nu1 < 1", GSL_EDOM);
    }
  if (nu2 < 1.0)
    {
      CDF_ERROR ("nu2 < 1", GSL_EDOM);
    }

  if (Q > 0.5)
    {
      y = gsl_cdf_beta_Qinv (Q, nu1 / 2.0, nu2 / 2.0);

      result = nu2 * y / (nu1 * (1.0 - y));
    }
  else
    {
      y = gsl_cdf_beta_Pinv (Q, nu2 / 2.0, nu1 / 2.0);

      result = nu2 * (1 - y) / (nu1 * y);
    }

  return result;
}


text:
/*
** use median of trials
**
** G.Lohmann, May 2016
*/
#include <viaio/Vlib.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_statistics.h>

#include <gsl/gsl_sort.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_sort_vector.h>

extern void VNormalize(float *data,int nt,VBoolean stddev);
extern void GetRank(float *data,gsl_vector *v,gsl_permutation *perm,gsl_permutation *rank);

void GetMedian(gsl_matrix_float **X1,gsl_matrix_float **X2,int *table,int n,gsl_matrix_float *SNR,int metric)
{
  long j,nvox=X1[0]->size1;
  long k,s,nt=X1[0]->size2;

  double *data = (double *) VCalloc(n,sizeof(double));
  double median = 0;

  gsl_vector *vec = NULL;
  gsl_permutation *perm = NULL;
  gsl_permutation *rank = NULL;
  if (metric == 1) {
    vec = gsl_vector_calloc(nt);
    perm = gsl_permutation_alloc(nt);
    rank = gsl_permutation_alloc(nt);
  }

  /* for each voxel ... */
  for (j=0; j<nvox; j++) {

    for (k=0; k<nt; k++) {
      for (s=0; s<n; s++) {
	if (table[s] == 0)
	  data[s] = (double)gsl_matrix_float_get(X1[s],j,k);
	else
	  data[s] = (double)gsl_matrix_float_get(X2[s],j,k);
      }
      gsl_sort(data,1,n);
      median = gsl_stats_median_from_sorted_data (data,1,(size_t) n);
      gsl_matrix_float_set(SNR,j,k,(float)median);
    }
    if (metric == 0) {  /* pearson correlation */
      float *qq = gsl_matrix_float_ptr(SNR,j,0);
      VNormalize(qq,nt,TRUE);
    }
    if (metric == 1) {  /* spearman ranks      */
      float *qq = gsl_matrix_float_ptr(SNR,j,0);
      GetRank(qq,vec,perm,rank);
    }
  }
  VFree(data);

  if (metric == 1) {
    gsl_permutation_free (perm);
    gsl_permutation_free (rank);
    gsl_vector_free(vec);
  }
}

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
#pragma once

#include "movie_writer.h"

#include <boost/math/tools/roots.hpp>
#include <cxxopts.hpp>
#include <glm/glm.hpp>
#include <glm/gtc/matrix_transform.hpp>
#include <glm/gtx/rotate_vector.hpp>
#include <glm/gtx/string_cast.hpp>
#include <gsl/gsl_integration.h>
#include <omp.h>
#include <opencv2/opencv.hpp>

#include <QApplication>
#include <QCheckBox>
#include <QGraphicsView>
#include <QGuiApplication>
#include <QHBoxLayout>
#include <QLineEdit>
#include <QMainWindow>
#include <QMenuBar>
#include <QPushButton>
#include <QVBoxLayout>
#include <QFileDialog>
#include <QGraphicsPixmapItem>
#include <QPixmap>

#include <array>
#include <chrono>
#include <cmath>
#include <filesystem>
#include <iostream>
#include <random>
#include <thread>
#include <vector>

max_stars_repo_name: evopen/blackness

max_issues_repo_name: evopen/blackness

max_forks_repo_name: evopen/blackness


text:
/*
 * Copyright (c) 2016-2021 lymastee, All rights reserved.
 * Contact: lymastee@hotmail.com
 *
 * This file is part of the gslib project.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */

#pragma once

#ifndef file_d50c3538_0c5e_49a3_bc6f_5d6bc6144520_h
#define file_d50c3538_0c5e_49a3_bc6f_5d6bc6144520_h

#include <stdio.h>
#include <gslib/config.h>
#include <gslib/string.h>

__gslib_begin__

typedef FILE*   fileptr;
typedef fpos_t  fpos;

class file
{
protected:
    fileptr     _file;

public:
    file(): _file(0) {}
    file(const gchar* name, const gchar* mode): _file(0) { open(name, mode); }
    ~file() { close(); }

#ifndef _UNICODE

    void open(const gchar* name, const gchar* mode) { fopen_s(&_file, name, mode); }
    void open_text(const gchar* name, bool readonly)
    {
        readonly ? fopen_s(&_file, name, _t("r")) :
            fopen_s(&_file, name, _t("w"));
    }
    int read(gchar buf[], int len)
    {
        assert(buf && len);
        return fgets(buf, len, _file) ? strlen(buf) : 0;
    }
    int write(const gchar str[])
    {
        assert(str);
        return (int)fputs(str, _file);
    }

#else

    void open(const gchar* name, const gchar* mode) { _wfopen_s(&_file, name, mode); }
    void open_text(const gchar* name, bool readonly)
    {
        readonly ? _wfopen_s(&_file, name, _t("r,ccs=UNICODE")) :
            _wfopen_s(&_file, name, _t("w,ccs=UNICODE"));
    }
    int read(gchar buf[], int len)
    {
        assert(buf && len);
        return fgetws(buf, len, _file) ? wcslen(buf) : 0;
    }
    int write(const gchar str[])
    {
        assert(str);
        return fputws(str, _file);
    }

#endif

    bool is_valid() const { return _file != 0; }
    void flush() { fflush(_file); }
    void rewind() { ::rewind(_file); }
    int current() const { return ftell(_file); }
    int seek(int offset, int sps) { return fseek(_file, offset, sps); }
    int size()
    {
        int old = current();
        seek(0, SEEK_END);
        int pos = current();
        seek(old, SEEK_SET);
        return pos;
    }
    void close()
    {
        if (_file) {
            fclose(_file);
            _file = 0;
        }
    }
    int read_fillup(gchar buf[], int len)
    {
        assert(buf && len > 0);
        int acc = 0;
        while(acc < len) {
            int r = read(buf + acc, len - acc);
            if(r <= 0)
                break;
            acc += r;
        }
        return acc;
    }
    int read(string& buf, int start, int len)
    {
        int lento = start + len;
        if(lento && buf.size() < lento)
            buf.resize(lento);
        return read(&buf.front() + start, len);
    }
    void read_all(string& buf)
    {
        int total = size() + 1;     /* a weird problem about fgetws. */
        if(buf.size() < total)
            buf.resize(total);
        for(int i = 0, j; j = read(buf, i, total-i); i += j);
        /* right trim */
        int p = buf.find_last_not_of(_t('\0'));
        if(p != string::npos) {
            p += 1;
            if(p <= buf.length())
                buf.resize(p);
        }
    }
    int write(const string& buf) { return write(&buf.front()); }
    int get(byte buf[], int len) { return (int)fread_s(buf, len, 1, len, _file); }
    int get(word buf[], int len) { return (int)fread_s(buf, len * 2, 2, len, _file); }
    int get(dword buf[], int len) { return (int)fread_s(buf, len * 4, 4, len, _file); }
    int get(qword buf[], int len) { return (int)fread_s(buf, len * 8, 8, len, _file); }
    int put(const byte buf[], int len) { return (int)fwrite(buf, 1, len, _file); }
    int put(const word buf[], int len) { return (int)fwrite(buf, 2, len, _file); }
    int put(const dword buf[], int len) { return (int)fwrite(buf, 4, len, _file); }
    int put(const qword buf[], int len) { return (int)fwrite(buf, 8, len, _file); }

public:
    template<int _szelem>
    int put_elem(const void* buf, int len) {}
    template<>
    int put_elem<1>(const void* buf, int len) { return put((const byte*)buf, len); }
    template<>
    int put_elem<2>(const void* buf, int len) { return put((const word*)buf, len); }
    template<>
    int put_elem<4>(const void* buf, int len) { return put((const dword*)buf, len); }
    template<>
    int put_elem<8>(const void* buf, int len) { return put((const qword*)buf, len); }
};

__gslib_end__

#endif

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


text:
#include <assert.h>

#include <gsl/gsl_statistics_double.h>
#include <hdf5.h>
#include <hdf5_hl.h>

#include <qdm.h>

qdm_evaluation *
qdm_evaluation_new(const qdm_parameters *parameters)
{
  qdm_evaluation *e = malloc(sizeof(qdm_evaluation));

  e->parameters = *parameters;

  e->rng = gsl_rng_alloc(gsl_rng_default);

  e->years_min = 0;
  e->years_max = 0;

  e->lower_bound = 0;
  e->upper_bound = 0;

  e->waic = 0;
  e->pwaic = 0;
  e->dic = 0;
  e->pd = 0;

  e->t = NULL;

  e->mcmc = NULL;

  e->theta_bar = NULL;
  e->theta_star_bar = NULL;

  e->xi_cov = NULL;
  e->theta_star_cov = NULL;

  e->xi_low_bar = 0;
  e->xi_high_bar = 0;

  e->theta_acc = NULL;
  e->xi_acc = NULL;

  e->m_knots = NULL;

  e->rng_seed = 0;

  e->elapsed = 0;

  return e;
}

void
qdm_evaluation_free(qdm_evaluation *e)
{
  if (e == NULL) {
    return;
  }

  e->rng_seed = 0;

  gsl_vector_free(e->m_knots);
  e->m_knots = NULL;

  gsl_vector_free(e->xi_acc);
  e->xi_acc = NULL;

  gsl_matrix_free(e->theta_acc);
  e->theta_acc = NULL;

  gsl_matrix_free(e->theta_star_cov);
  e->theta_bar = NULL;

  gsl_matrix_free(e->xi_cov);
  e->theta_bar = NULL;

  gsl_matrix_free(e->theta_star_bar);
  e->theta_bar = NULL;

  gsl_matrix_free(e->theta_bar);
  e->theta_bar = NULL;

  qdm_mcmc_free(e->mcmc);
  e->mcmc = NULL;

  qdm_tau_free(e->t);
  e->t = NULL;

  gsl_rng_free(e->rng);
  e->rng = NULL;

  free(e);
}

void
qdm_evaluation_fprint(
    FILE *f,
    const qdm_evaluation *e
)
{
  const char *prefix = "  ";

  fprintf(f, "%sparameters:\n", prefix);
  qdm_parameters_fprint(f, &e->parameters);

  fprintf(f, "%syears_min: %f\n", prefix, e->years_min);
  fprintf(f, "%syears_max: %f\n", prefix, e->years_max);

  fprintf(f, "%slower_bound: %f\n", prefix, e->lower_bound);
  fprintf(f, "%supper_bound: %f\n", prefix, e->upper_bound);

  fprintf(f, "%swaic: %f\n", prefix, e->waic);
  fprintf(f, "%spwaic: %f\n", prefix, e->pwaic);
  fprintf(f, "%sdic: %f\n", prefix, e->dic);
  fprintf(f, "%spd: %f\n", prefix, e->pd);

  fprintf(f, "%smcmc r.s: %zu\n", prefix, e->mcmc->r.s);

  fprintf(f, "%stheta_bar:\n", prefix);
  qdm_matrix_csv_fwrite(f, e->theta_bar);

  fprintf(f, "%stheta_star_bar:\n", prefix);
  qdm_matrix_csv_fwrite(f, e->theta_star_bar);

  fprintf(f, "%sxi_cov:\n", prefix);
  qdm_matrix_csv_fwrite(f, e->xi_cov);

  fprintf(f, "%stheta_star_cov:\n", prefix);
  qdm_matrix_csv_fwrite(f, e->theta_star_cov);

  fprintf(f, "%sxi_low_bar: %f\n", prefix, e->xi_low_bar);
  fprintf(f, "%sxi_high_bar: %f\n", prefix, e->xi_high_bar);

  fprintf(f, "%stheta_acc:\n", prefix);
  qdm_matrix_csv_fwrite(f, e->theta_acc);

  fprintf(f, "%sxi_acc:\n", prefix);
  qdm_vector_csv_fwrite(f, e->xi_acc);

  fprintf(f, "%sm_knots:\n", prefix);
  qdm_vector_csv_fwrite(f, e->m_knots);

  fprintf(f, "%srng_seed: %lu\n", prefix, e->rng_seed);

  fprintf(f, "%selapsed: %f\n", prefix, e->elapsed);
}

int
qdm_evaluation_read(
    hid_t id,
    qdm_evaluation **e
)
{
  int status = 0;

  hid_t parameters_group = -1;
  hid_t mcmc_group = -1;

  parameters_group = H5Gopen(id, "parameters", H5P_DEFAULT);
  if (parameters_group < 0) {
    status = parameters_group;

    goto cleanup;
  }

  qdm_parameters parameters = {0};
  status = qdm_parameters_read(parameters_group, &parameters);
  if (status != 0) {
    goto cleanup;
  }

  *e = qdm_evaluation_new(&parameters);

#define READ(t, n) \
  status = H5LTread_dataset_##t(id, #n, &(*e)->n); \
  if (status < 0) { \
    goto cleanup; \
  }

  READ(double, years_min);
  READ(double, years_max);

  READ(double, lower_bound);
  READ(double, upper_bound);

  READ(double, waic);
  READ(double, pwaic);
  READ(double, dic);
  READ(double, pd);

  mcmc_group = H5Gopen(id, "mcmc", H5P_DEFAULT);
  if (mcmc_group < 0) {
    status = mcmc_group;

    goto cleanup;
  }

  status = qdm_mcmc_read(mcmc_group, &(*e)->mcmc);
  if (status != 0) {
    goto cleanup;
  }

  status = qdm_matrix_hd5_read(id, "theta_bar", &(*e)->theta_bar);
  if (status != 0) {
    goto cleanup;
  }

  status = qdm_matrix_hd5_read(id, "theta_star_bar", &(*e)->theta_star_bar);
  if (status != 0) {
    goto cleanup;
  }

  status = qdm_matrix_hd5_read(id, "xi_cov", &(*e)->xi_cov);
  if (status != 0) {
    goto cleanup;
  }

  status = qdm_matrix_hd5_read(id, "theta_star_cov", &(*e)->theta_star_cov);
  if (status != 0) {
    goto cleanup;
  }

  READ(double, xi_high_bar);
  READ(double, xi_low_bar);

  status = qdm_matrix_hd5_read(id, "theta_acc", &(*e)->theta_acc);
  if (status != 0) {
    goto cleanup;
  }

  status = qdm_vector_hd5_read(id, "xi_acc", &(*e)->xi_acc);
  if (status != 0) {
    goto cleanup;
  }

  status = qdm_vector_hd5_read(id, "m_knots", &(*e)->m_knots);
  if (status != 0) {
    goto cleanup;
  }

  status = H5LTread_dataset(id, "rng_seed", H5T_NATIVE_ULONG, &(*e)->rng_seed);
  if (status < 0) {
    goto cleanup;
  }

  READ(double, elapsed);

#undef READ

  (*e)->t = qdm_tau_alloc(
      (*e)->parameters.tau_table,
      (*e)->parameters.tau_low,
      (*e)->parameters.tau_high,
      (*e)->parameters.spline_df,
      (*e)->m_knots
  );

cleanup:
  if (parameters_group >= 0) {
    H5Gclose(parameters_group);
  }

  if (mcmc_group >= 0) {
    H5Gclose(mcmc_group);
  }

  return status;
}

int
qdm_evaluation_write(
    hid_t id,
    const qdm_evaluation *e
)
{
  int status = 0;

  hid_t parameters_group = -1;
  hid_t mcmc_group = -1;

#define WRITE_DOUBLE(n) \
  status = qdm_double_write(id, #n, e->n); \
  if (status != 0) { \
    goto cleanup; \
  }

  parameters_group = qdm_data_create_group(id, "parameters");
  if (parameters_group < 0) {
    status = parameters_group;

    goto cleanup;
  }

  status = qdm_parameters_write(parameters_group, &e->parameters);
  if (status != 0) {
    goto cleanup;
  }

  WRITE_DOUBLE(years_min);
  WRITE_DOUBLE(years_max);

  WRITE_DOUBLE(lower_bound);
  WRITE_DOUBLE(upper_bound);

  WRITE_DOUBLE(waic);
  WRITE_DOUBLE(pwaic);
  WRITE_DOUBLE(dic);
  WRITE_DOUBLE(pd);

  mcmc_group = qdm_data_create_group(id, "mcmc");
  if (mcmc_group < 0) {
    status = mcmc_group;

    goto cleanup;
  }

  status = qdm_mcmc_write(mcmc_group, e->mcmc);
  if (status != 0) {
    goto cleanup;
  }

  status = qdm_matrix_hd5_write(id, "theta_bar", e->theta_bar);
  if (status != 0) {
    goto cleanup;
  }

  status = qdm_matrix_hd5_write(id, "theta_star_bar", e->theta_star_bar);
  if (status != 0) {
    goto cleanup;
  }

  status = qdm_matrix_hd5_write(id, "xi_cov", e->xi_cov);
  if (status != 0) {
    goto cleanup;
  }

  status = qdm_matrix_hd5_write(id, "theta_star_cov", e->theta_star_cov);
  if (status != 0) {
    goto cleanup;
  }

  WRITE_DOUBLE(xi_high_bar);
  WRITE_DOUBLE(xi_low_bar);

  status = qdm_matrix_hd5_write(id, "theta_acc", e->theta_acc);
  if (status != 0) {
    goto cleanup;
  }

  status = qdm_vector_hd5_write(id, "xi_acc", e->xi_acc);
  if (status != 0) {
    goto cleanup;
  }

  status = qdm_vector_hd5_write(id, "m_knots", e->m_knots);
  if (status != 0) {
    goto cleanup;
  }

  {
    hsize_t dims[1] = {1};
    unsigned long int data[1] = {e->rng_seed};

    status = H5LTmake_dataset(id, "rng_seed", 1, dims, H5T_NATIVE_ULONG, data);
    if (status < 0) {
      goto cleanup;
    }
  }

  WRITE_DOUBLE(elapsed);

#undef WRITE_DOUBLE

cleanup:
  if (parameters_group >= 0) {
    H5Gclose(parameters_group);
  }

  if (mcmc_group >= 0) {
    H5Gclose(mcmc_group);
  }

  return status;
}

int
qdm_evaluation_run(
    qdm_evaluation *e,
    const gsl_matrix *data,
    size_t data_year_idx,
    size_t data_month_idx,
    size_t data_value_idx
)
{
  int status = 0;

  clock_t started_at = clock();

  gsl_vector *years = NULL;
  gsl_vector *years_norm = NULL;

  gsl_vector *values = NULL;
  gsl_vector *values_sorted = NULL;

  gsl_vector *interior_knots = NULL;

  gsl_vector *m_knots = NULL;
  gsl_matrix *theta = NULL;
  gsl_vector *ll = NULL;
  gsl_vector *tau = NULL;

  gsl_matrix *theta_tune = NULL;
  gsl_matrix *theta_acc = NULL;
  gsl_vector *xi_tune = NULL;
  gsl_vector *xi_acc = NULL;

  gsl_vector *xi = gsl_vector_alloc(2);
  gsl_vector_set(xi, 0, e->parameters.xi_low);
  gsl_vector_set(xi, 1, e->parameters.xi_high);

  /* Initialize RNG */

  gsl_rng_set(e->rng, e->parameters.rng_seed);

  /* Normalize Years */

  {
    years = qdm_matrix_filter(data, data_month_idx, e->parameters.month, data_year_idx);
    years_norm = qdm_matrix_filter(data, data_month_idx, e->parameters.month, data_year_idx);

    gsl_vector_minmax(years, &e->years_min, &e->years_max);

    double years_range = e->years_max - e->years_min;

    status = gsl_vector_add_constant(years_norm, -e->years_min);
    if (status != 0) {
      goto cleanup;
    }

    status = gsl_vector_scale(years_norm, 1 / years_range);
    if (status != 0) {
      goto cleanup;
    }
  }

  /* Extract Values */

  {
    values = qdm_matrix_filter(data, data_month_idx, e->parameters.month, data_value_idx);
    values_sorted = qdm_vector_sorted(values);

    e->lower_bound = gsl_vector_get(values_sorted, 0) - e->parameters.bound;
    e->upper_bound = gsl_vector_get(values_sorted, values_sorted->size - 1) + e->parameters.bound;
  }

  /* Optimize Knots */

  {
    gsl_vector *middle = qdm_vector_seq(e->parameters.tau_low, e->parameters.tau_high, 0.005);

    double possible_low = e->parameters.tau_low + 0.1;
    double possible_high = e->parameters.tau_high - 0.1;
    double possible_delta = (possible_high - possible_low) / 19;

    gsl_vector *possible_knots = qdm_vector_seq(possible_low, possible_high, possible_delta);

    interior_knots = gsl_vector_calloc(e->parameters.knot);

    status = qdm_knots_optimize(
        interior_knots,

        e->rng,

        values_sorted,
        middle,
        possible_knots,

        e->parameters.knot_try,
        e->parameters.spline_df
    );
    if (status != 0) {
      goto cleanup;
    }

    gsl_vector_free(possible_knots);
    gsl_vector_free(middle);
  }

  /* Optimize Theta */

  {
    size_t m = e->parameters.spline_df + e->parameters.knot;

    m_knots = qdm_knots_vector(e->parameters.spline_df, interior_knots);
    e->t = qdm_tau_alloc(
        e->parameters.tau_table,
        e->parameters.tau_low,
        e->parameters.tau_high,
        e->parameters.spline_df,
        m_knots
    );

    gsl_vector *middle = qdm_vector_seq(e->parameters.tau_low, e->parameters.tau_high, 0.01);

    gsl_matrix *ix = gsl_matrix_alloc(middle->size, m+1);
    qdm_ispline_matrix(ix, middle, e->parameters.spline_df, m_knots);

    theta = gsl_matrix_alloc(2, ix->size2);
    gsl_vector_view theta0 = gsl_matrix_row(theta, 0);
    gsl_vector_view theta1 = gsl_matrix_row(theta, 1);

    size_t theta_pivot = qdm_vector_greater_than(years_norm, 0.3);
    gsl_vector_view values0 = gsl_vector_subvector(values, 0, theta_pivot);
    gsl_vector_view values1 = gsl_vector_subvector(values, theta_pivot, values->size - theta_pivot);
    assert(values0.vector.size + values1.vector.size == values->size);

    gsl_vector *eq0 = qdm_vector_quantile(&values0.vector, middle);
    status = qdm_theta_optimize(&theta0.vector, eq0, ix);
    if (status != 0) {
      goto cleanup;
    }

    gsl_vector *eq1 = NULL;

    if (e->parameters.truncate) {
      gsl_vector_set_zero(&theta1.vector);
    } else {
      eq1 = qdm_vector_quantile(&values1.vector, middle);
      status = qdm_theta_optimize(&theta1.vector, eq1, ix);
      if (status != 0) {
        goto cleanup;
      }

      status = gsl_vector_sub(&theta1.vector, &theta0.vector);
      if (status != 0) {
        goto cleanup;
      }
    }

    qdm_theta_matrix_constrain(theta, e->parameters.theta_min);

    ll = gsl_vector_alloc(values->size);
    tau = gsl_vector_alloc(values->size);

    for (size_t i = 0; i < values->size; i++) {
      qdm_logl_3(
          &ll->data[i],
          &tau->data[i],

          gsl_vector_get(years_norm, i),
          gsl_vector_get(values, i),

          e->t,
          xi,

          theta
      );
    }

    gsl_vector_free(eq1);
    gsl_vector_free(eq0);
    gsl_matrix_free(ix);
    gsl_vector_free(middle);
  }

  /* MCMC */

  {
    qdm_mcmc_parameters mcmc_p = {
      .rng = e->rng,

      .burn = e->parameters.burn,
      .iter = e->parameters.iter,
      .thin = e->parameters.thin,

      .acc_check = e->parameters.acc_check,
      .spline_df = e->parameters.spline_df,

      .xi    = xi,
      .ll    = ll,
      .tau   = tau,
      .theta = theta,

      .xi_prior_mean = e->parameters.xi_prior_mean,
      .xi_prior_var  = e->parameters.xi_prior_var,
      .xi_tune_sd    = e->parameters.xi_tune_sd,

      .theta_min     = e->parameters.theta_min,
      .theta_tune_sd = e->parameters.theta_tune_sd,

      .x = years_norm,
      .y = values,

      .t = e->t,
      .m_knots = m_knots,

      .truncate = e->parameters.truncate,
    };

    e->mcmc = qdm_mcmc_alloc(mcmc_p);

    status = qdm_mcmc_run(e->mcmc);
    if (status != 0) {
      goto cleanup;
    }
  }

  /* Model Diagnostics */

  {
    // WAIC

    gsl_vector *ll_exp = gsl_vector_calloc(e->mcmc->r.s);

    gsl_vector *lppd = gsl_vector_calloc(values->size);
    gsl_vector *vlog = gsl_vector_calloc(values->size);

    for (size_t j = 0; j < values->size; j++) {
      gsl_vector_view ll_k = qdm_ijk_get_k(e->mcmc->r.ll, 0, j);

      for (size_t i = 0; i < ll_k.vector.size; i++) {
        gsl_vector_set(ll_exp, i, exp(gsl_vector_get(&ll_k.vector, i)));
      }

      gsl_vector_set(lppd, j, log(gsl_stats_mean(ll_exp->data, ll_exp->stride, ll_exp->size)));
      gsl_vector_set(vlog, j, gsl_stats_variance(ll_k.vector.data, ll_k.vector.stride, ll_k.vector.size));
    }

    status = gsl_vector_sub(lppd, vlog);
    if (status != 0) {
      goto cleanup;
    }

    e->waic = -2 * qdm_vector_sum(lppd);
    e->pwaic = qdm_vector_sum(vlog);

    // DIC

    gsl_vector *ll_sum = gsl_vector_alloc(e->mcmc->r.s);

    for (size_t k = 0; k < e->mcmc->r.s; k++) {
      gsl_matrix_view ll_ij_m = qdm_ijk_get_ij(e->mcmc->r.ll, k);
      gsl_vector_view ll_ij = gsl_matrix_row(&ll_ij_m.matrix, 0);

      gsl_vector_set(ll_sum, k, qdm_vector_sum(&ll_ij.vector));
    }

    status = gsl_vector_scale(ll_sum, -2);
    if (status != 0) {
      goto cleanup;
    }

    double d_bar = gsl_stats_mean(ll_sum->data, ll_sum->stride, ll_sum->size);
    double d_theta_bar = 0;

    gsl_vector_view xi_low = qdm_ijk_get_k(e->mcmc->r.xi, 0, 0);
    gsl_vector_view xi_high = qdm_ijk_get_k(e->mcmc->r.xi, 0, 1);

    double xi_low_bar = gsl_stats_mean(xi_low.vector.data, xi_low.vector.stride, xi_low.vector.size);
    double xi_high_bar = gsl_stats_mean(xi_high.vector.data, xi_high.vector.stride, xi_high.vector.size);

    gsl_vector *xi_bar = gsl_vector_alloc(2);
    gsl_vector_set(xi_bar, 0, xi_low_bar);
    gsl_vector_set(xi_bar, 1, xi_high_bar);

    gsl_matrix *theta_bar = gsl_matrix_alloc(theta->size1, theta->size2);

    for (size_t i = 0; i < e->mcmc->r.theta->size1; i++) {
      for (size_t j = 0; j < e->mcmc->r.theta->size2; j++) {
        gsl_vector_view theta_k = qdm_ijk_get_k(e->mcmc->r.theta, i, j);

        gsl_matrix_set(theta_bar, i, j, gsl_stats_mean(theta_k.vector.data, theta_k.vector.stride, theta_k.vector.size));
      }
    }

    gsl_matrix *theta_star_bar = gsl_matrix_alloc(theta->size1, theta->size2);

    for (size_t i = 0; i < e->mcmc->r.theta_star->size1; i++) {
      for (size_t j = 0; j < e->mcmc->r.theta_star->size2; j++) {
        gsl_vector_view theta_k = qdm_ijk_get_k(e->mcmc->r.theta_star, i, j);

        gsl_matrix_set(theta_star_bar, i, j, gsl_stats_mean(theta_k.vector.data, theta_k.vector.stride, theta_k.vector.size));
      }
    }

    for (size_t i = 0; i < values->size; i++) {
      double theta_bar_ll = 0;
      double tau_tmp = 0;

      qdm_logl_3(
          &theta_bar_ll,
          &tau_tmp,

          gsl_vector_get(years_norm, i),
          gsl_vector_get(values, i),

          e->t,
          xi_bar,

          theta_bar
      );

      d_theta_bar = d_theta_bar + -2 * theta_bar_ll;
    }

    e->pd = d_bar - d_theta_bar;
    e->dic = e->pd + d_bar;

    e->theta_bar = qdm_matrix_copy(theta_bar);
    e->theta_star_bar = qdm_matrix_copy(theta_star_bar);

    e->xi_cov = qdm_ijk_cov(e->mcmc->r.xi);
    e->theta_star_cov = qdm_ijk_cov(e->mcmc->r.theta_star);

    e->xi_low_bar = xi_low_bar;
    e->xi_high_bar = xi_high_bar;

    gsl_vector_free(xi_bar);
    gsl_matrix_free(theta_star_bar);
    gsl_matrix_free(theta_bar);
    gsl_vector_free(ll_sum);

    gsl_vector_free(vlog);
    gsl_vector_free(lppd);
    gsl_vector_free(ll_exp);
  }

  /* Results */
  {
    e->theta_acc = qdm_matrix_copy(e->mcmc->w.theta_acc);
    e->xi_acc = qdm_vector_copy(e->mcmc->w.xi_acc);

    e->m_knots = qdm_vector_copy(m_knots);

    e->rng_seed = gsl_rng_get(e->rng);
  }

  gsl_vector_free(years);
  gsl_vector_free(years_norm);
  gsl_vector_free(values);
  gsl_vector_free(values_sorted);
  gsl_vector_free(interior_knots);
  gsl_vector_free(m_knots);
  gsl_matrix_free(theta);
  gsl_vector_free(ll);
  gsl_vector_free(tau);

  gsl_matrix_free(theta_tune);
  gsl_matrix_free(theta_acc);
  gsl_vector_free(xi_tune);
  gsl_vector_free(xi_acc);

cleanup:
  e->elapsed = ((double)(clock() - started_at)) / CLOCKS_PER_SEC;

  return status;
}

max_stars_repo_name: calebcase/qdm

max_issues_repo_name: calebcase/qdm

max_forks_repo_name: calebcase/qdm