text:
/* randist/gsl_randist.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 James Theiler, 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.
 */

#ifndef __GSL_RANDIST_H__
#define __GSL_RANDIST_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif
#include <gsl/gsl_rng.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.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

GSL_FUN unsigned int gsl_ran_bernoulli (const gsl_rng * r, double p);
GSL_FUN double gsl_ran_bernoulli_pdf (const unsigned int k, double p);

GSL_FUN double gsl_ran_beta (const gsl_rng * r, const double a, const double b);
GSL_FUN double gsl_ran_beta_pdf (const double x, const double a, const double b);

GSL_FUN unsigned int gsl_ran_binomial (const gsl_rng * r, double p, unsigned int n);
GSL_FUN unsigned int gsl_ran_binomial_knuth (const gsl_rng * r, double p, unsigned int n);
GSL_FUN unsigned int gsl_ran_binomial_tpe (const gsl_rng * r, double p, unsigned int n);
GSL_FUN double gsl_ran_binomial_pdf (const unsigned int k, const double p, const unsigned int n);

GSL_FUN double gsl_ran_exponential (const gsl_rng * r, const double mu);
GSL_FUN double gsl_ran_exponential_pdf (const double x, const double mu);

GSL_FUN double gsl_ran_exppow (const gsl_rng * r, const double a, const double b);
GSL_FUN double gsl_ran_exppow_pdf (const double x, const double a, const double b);

GSL_FUN double gsl_ran_cauchy (const gsl_rng * r, const double a);
GSL_FUN double gsl_ran_cauchy_pdf (const double x, const double a);

GSL_FUN double gsl_ran_chisq (const gsl_rng * r, const double nu);
GSL_FUN double gsl_ran_chisq_pdf (const double x, const double nu);

GSL_FUN void gsl_ran_dirichlet (const gsl_rng * r, const size_t K, const double alpha[], double theta[]);
GSL_FUN double gsl_ran_dirichlet_pdf (const size_t K, const double alpha[], const double theta[]);
GSL_FUN double gsl_ran_dirichlet_lnpdf (const size_t K, const double alpha[], const double theta[]);

GSL_FUN double gsl_ran_erlang (const gsl_rng * r, const double a, const double n);
GSL_FUN double gsl_ran_erlang_pdf (const double x, const double a, const double n);

GSL_FUN double gsl_ran_fdist (const gsl_rng * r, const double nu1, const double nu2);
GSL_FUN double gsl_ran_fdist_pdf (const double x, const double nu1, const double nu2);

GSL_FUN double gsl_ran_flat (const gsl_rng * r, const double a, const double b);
GSL_FUN double gsl_ran_flat_pdf (double x, const double a, const double b);

GSL_FUN double gsl_ran_gamma (const gsl_rng * r, const double a, const double b);
GSL_FUN double gsl_ran_gamma_int (const gsl_rng * r, const unsigned int a);
GSL_FUN double gsl_ran_gamma_pdf (const double x, const double a, const double b);
GSL_FUN double gsl_ran_gamma_mt (const gsl_rng * r, const double a, const double b);
GSL_FUN double gsl_ran_gamma_knuth (const gsl_rng * r, const double a, const double b);

GSL_FUN double gsl_ran_gaussian (const gsl_rng * r, const double sigma);
GSL_FUN double gsl_ran_gaussian_ratio_method (const gsl_rng * r, const double sigma);
GSL_FUN double gsl_ran_gaussian_ziggurat (const gsl_rng * r, const double sigma);
GSL_FUN double gsl_ran_gaussian_pdf (const double x, const double sigma);

GSL_FUN double gsl_ran_ugaussian (const gsl_rng * r);
GSL_FUN double gsl_ran_ugaussian_ratio_method (const gsl_rng * r);
GSL_FUN double gsl_ran_ugaussian_pdf (const double x);

GSL_FUN double gsl_ran_gaussian_tail (const gsl_rng * r, const double a, const double sigma);
GSL_FUN double gsl_ran_gaussian_tail_pdf (const double x, const double a, const double sigma);

GSL_FUN double gsl_ran_ugaussian_tail (const gsl_rng * r, const double a);
GSL_FUN double gsl_ran_ugaussian_tail_pdf (const double x, const double a);

GSL_FUN void gsl_ran_bivariate_gaussian (const gsl_rng * r, double sigma_x, double sigma_y, double rho, double *x, double *y);
GSL_FUN double gsl_ran_bivariate_gaussian_pdf (const double x, const double y, const double sigma_x, const double sigma_y, const double rho);

GSL_FUN int gsl_ran_multivariate_gaussian (const gsl_rng * r, const gsl_vector * mu, const gsl_matrix * L, gsl_vector * result);
GSL_FUN int gsl_ran_multivariate_gaussian_log_pdf (const gsl_vector * x,
                                           const gsl_vector * mu,
                                           const gsl_matrix * L,
                                           double * result,
                                           gsl_vector * work);
GSL_FUN int gsl_ran_multivariate_gaussian_pdf (const gsl_vector * x,
                                       const gsl_vector * mu,
                                       const gsl_matrix * L,
                                       double * result,
                                       gsl_vector * work);
GSL_FUN int gsl_ran_multivariate_gaussian_mean (const gsl_matrix * X, gsl_vector * mu_hat);
GSL_FUN int gsl_ran_multivariate_gaussian_vcov (const gsl_matrix * X, gsl_matrix * sigma_hat);

GSL_FUN int gsl_ran_wishart (const gsl_rng * r,
                     const double df,
                     const gsl_matrix * L,
                     gsl_matrix * result,
                     gsl_matrix * work);
GSL_FUN int gsl_ran_wishart_log_pdf (const gsl_matrix * X,
                             const gsl_matrix * L_X,
                             const double df,
                             const gsl_matrix * L,
                             double * result,
                             gsl_matrix * work);
GSL_FUN int gsl_ran_wishart_pdf (const gsl_matrix * X,
                         const gsl_matrix * L_X,
                         const double df,
                         const gsl_matrix * L,
                         double * result,
                         gsl_matrix * work);

GSL_FUN double gsl_ran_landau (const gsl_rng * r);
GSL_FUN double gsl_ran_landau_pdf (const double x);

GSL_FUN unsigned int gsl_ran_geometric (const gsl_rng * r, const double p);
GSL_FUN double gsl_ran_geometric_pdf (const unsigned int k, const double p);

GSL_FUN unsigned int gsl_ran_hypergeometric (const gsl_rng * r, unsigned int n1, unsigned int n2, unsigned int t);
GSL_FUN double gsl_ran_hypergeometric_pdf (const unsigned int k, const unsigned int n1, const unsigned int n2, unsigned int t);

GSL_FUN double gsl_ran_gumbel1 (const gsl_rng * r, const double a, const double b);
GSL_FUN double gsl_ran_gumbel1_pdf (const double x, const double a, const double b);

GSL_FUN double gsl_ran_gumbel2 (const gsl_rng * r, const double a, const double b);
GSL_FUN double gsl_ran_gumbel2_pdf (const double x, const double a, const double b);

GSL_FUN double gsl_ran_logistic (const gsl_rng * r, const double a);
GSL_FUN double gsl_ran_logistic_pdf (const double x, const double a);

GSL_FUN double gsl_ran_lognormal (const gsl_rng * r, const double zeta, const double sigma);
GSL_FUN double gsl_ran_lognormal_pdf (const double x, const double zeta, const double sigma);

GSL_FUN unsigned int gsl_ran_logarithmic (const gsl_rng * r, const double p);
GSL_FUN double gsl_ran_logarithmic_pdf (const unsigned int k, const double p);

GSL_FUN void gsl_ran_multinomial (const gsl_rng * r, const size_t K,
                          const unsigned int N, const double p[],
                          unsigned int n[] );
GSL_FUN double gsl_ran_multinomial_pdf (const size_t K,
                                const double p[], const unsigned int n[] );
GSL_FUN double gsl_ran_multinomial_lnpdf (const size_t K,
                           const double p[], const unsigned int n[] );

GSL_FUN unsigned int gsl_ran_negative_binomial (const gsl_rng * r, double p, double n);
GSL_FUN double gsl_ran_negative_binomial_pdf (const unsigned int k, const double p, double n);

GSL_FUN unsigned int gsl_ran_pascal (const gsl_rng * r, double p, unsigned int n);
GSL_FUN double gsl_ran_pascal_pdf (const unsigned int k, const double p, unsigned int n);

GSL_FUN double gsl_ran_pareto (const gsl_rng * r, double a, const double b);
GSL_FUN double gsl_ran_pareto_pdf (const double x, const double a, const double b);

GSL_FUN unsigned int gsl_ran_poisson (const gsl_rng * r, double mu);
GSL_FUN void gsl_ran_poisson_array (const gsl_rng * r, size_t n, unsigned int array[],
                            double mu);
GSL_FUN double gsl_ran_poisson_pdf (const unsigned int k, const double mu);

GSL_FUN double gsl_ran_rayleigh (const gsl_rng * r, const double sigma);
GSL_FUN double gsl_ran_rayleigh_pdf (const double x, const double sigma);

GSL_FUN double gsl_ran_rayleigh_tail (const gsl_rng * r, const double a, const double sigma);
GSL_FUN double gsl_ran_rayleigh_tail_pdf (const double x, const double a, const double sigma);

GSL_FUN double gsl_ran_tdist (const gsl_rng * r, const double nu);
GSL_FUN double gsl_ran_tdist_pdf (const double x, const double nu);

GSL_FUN double gsl_ran_laplace (const gsl_rng * r, const double a);
GSL_FUN double gsl_ran_laplace_pdf (const double x, const double a);

GSL_FUN double gsl_ran_levy (const gsl_rng * r, const double c, const double alpha);
GSL_FUN double gsl_ran_levy_skew (const gsl_rng * r, const double c, const double alpha, const double beta);

GSL_FUN double gsl_ran_weibull (const gsl_rng * r, const double a, const double b);
GSL_FUN double gsl_ran_weibull_pdf (const double x, const double a, const double b);

GSL_FUN void gsl_ran_dir_2d (const gsl_rng * r, double * x, double * y);
GSL_FUN void gsl_ran_dir_2d_trig_method (const gsl_rng * r, double * x, double * y);
GSL_FUN void gsl_ran_dir_3d (const gsl_rng * r, double * x, double * y, double * z);
GSL_FUN void gsl_ran_dir_nd (const gsl_rng * r, size_t n, double * x);

GSL_FUN void gsl_ran_shuffle (const gsl_rng * r, void * base, size_t nmembm, size_t size);
GSL_FUN int gsl_ran_choose (const gsl_rng * r, void * dest, size_t k, void * src, size_t n, size_t size) ;
GSL_FUN void gsl_ran_sample (const gsl_rng * r, void * dest, size_t k, void * src, size_t n, size_t size) ;

typedef struct {                /* struct for Walker algorithm */
    size_t K;
    size_t *A;
    double *F;
} gsl_ran_discrete_t;

GSL_FUN gsl_ran_discrete_t * gsl_ran_discrete_preproc (size_t K, const double *P);
GSL_FUN void gsl_ran_discrete_free(gsl_ran_discrete_t *g);
GSL_FUN size_t gsl_ran_discrete (const gsl_rng *r, const gsl_ran_discrete_t *g);
GSL_FUN double gsl_ran_discrete_pdf (size_t k, const gsl_ran_discrete_t *g);

__END_DECLS

#endif /* __GSL_RANDIST_H__ */

max_stars_repo_name: mgreter/astrometrylib

max_issues_repo_name: zzpwahaha/VimbaCamJILA

max_forks_repo_name: zzpwahaha/VimbaCamJILA


text:
#include <cblas.h>
#ifdef LAPACKE_FOLDER
#include <lapacke/lapacke.h>
#include <lapacke/lapacke_utils.h>
#else
#include <lapacke.h>
#include <lapacke_utils.h>
#endif

#include "math.h"
#include "stdbool.h"
#include "stdio.h"
#include "string.h"
#include "sv_matrix.h"

#include <limits.h>
#include <stdarg.h>

#include "linmath.h"

#ifdef _WIN32
#define SURVIVE_LOCAL_ONLY
#include <malloc.h>
#define alloca _alloca
#else
#define SURVIVE_LOCAL_ONLY __attribute__((visibility("hidden")))
#endif

#define SV_Error(code, msg) assert(0 && msg); // cv::error( code, msg, SV_Func, __FILE__, __LINE__ )

const int DECOMP_SVD = 1;
const int DECOMP_LU = 2;

void print_mat(const SvMat *M);

static size_t mat_size_bytes(const SvMat *mat) { return (size_t)sizeof(FLT) * mat->cols * mat->rows; }

#ifdef USE_FLOAT
#define cblas_gemm cblas_sgemm
#define cblas_symm cblas_ssymm
#define LAPACKE_getrs LAPACKE_sgetrs
#define LAPACKE_getrf LAPACKE_sgetrf
#define LAPACKE_getri LAPACKE_sgetri
#define LAPACKE_gelss LAPACKE_sgelss
#define LAPACKE_gesvd LAPACKE_sgesvd
#define LAPACKE_gesvd_work LAPACKE_sgesvd_work
#define LAPACKE_getri_work LAPACKE_sgetri_work
#define LAPACKE_ge_trans LAPACKE_sge_trans
#else
#define cblas_gemm cblas_dgemm
#define cblas_symm cblas_dsymm
#define LAPACKE_getrs LAPACKE_dgetrs
#define LAPACKE_getrf LAPACKE_dgetrf
#define LAPACKE_getri LAPACKE_dgetri
#define LAPACKE_gelss LAPACKE_dgelss
#define LAPACKE_gesvd LAPACKE_dgesvd
#define LAPACKE_gesvd_work LAPACKE_dgesvd_work
#define LAPACKE_getri_work LAPACKE_dgetri_work
#define LAPACKE_ge_trans LAPACKE_dge_trans
#endif

// dst = alpha * src1 * src2 + beta * src3 or dst = alpha * src2 * src1 + beta * src3 where src1 is symm
SURVIVE_LOCAL_ONLY void svSYMM(const SvMat *src1, const SvMat *src2, double alpha, const SvMat *src3, double beta,
							   SvMat *dst, bool src1First) {

	int rows1 = src1->rows;
	int cols1 = src1->cols;

	int rows2 = src2->rows;
	int cols2 = src2->cols;

	if (src3) {
		int rows3 = src3->rows;
		int cols3 = src3->cols;
		assert(rows3 == dst->rows);
		assert(cols3 == dst->cols);
	}

	// assert(src3 == 0 || beta != 0);
	assert(cols1 == rows2);
	assert(rows1 == dst->rows);
	assert(cols2 == dst->cols);

	lapack_int lda = src1->cols;
	lapack_int ldb = src2->cols;

	if (src3)
		svCopy(src3, dst, 0);
	else
		beta = 0;

	assert(SV_RAW_PTR(dst) != SV_RAW_PTR(src1));
	assert(SV_RAW_PTR(dst) != SV_RAW_PTR(src2));
	/*
		void cblas_dsymm(OPENBLAS_CONST enum CBLAS_ORDER Order,
						 OPENBLAS_CONST enum CBLAS_SIDE Side,
						 OPENBLAS_CONST enum CBLAS_UPLO Uplo,
						 OPENBLAS_CONST blasint M,
						 OPENBLAS_CONST blasint N,
						 OPENBLAS_CONST double alpha,
						 OPENBLAS_CONST double *A,
						 OPENBLAS_CONST blasint lda,
						 OPENBLAS_CONST double *B,
						 OPENBLAS_CONST blasint ldb,
						 OPENBLAS_CONST double beta,
						 double *C,
						 OPENBLAS_CONST blasint ldc);
	*/
	cblas_symm(CblasRowMajor, src1First ? CblasLeft : CblasRight, CblasUpper, dst->rows, dst->cols, alpha,
			   SV_RAW_PTR(src1), lda, SV_RAW_PTR(src2), ldb, beta, SV_RAW_PTR(dst), dst->cols);
}

// Special case dst = alpha * src2 * src1 * src2' + beta * src3
void mulBABt(const SvMat *src1, const SvMat *src2, double alpha, const SvMat *src3, double beta, SvMat *dst) {
	size_t dims = src2->rows;
	assert(src2->cols == src2->rows);
	SV_CREATE_STACK_MAT(tmp, dims, dims);

	// This has been profiled; and weirdly enough the SYMM version is slower for a 19x19 matrix. Guessing access order
	// or some other cache thing matters more than the additional 2x multiplications.
//#define USE_SYM
#ifdef USE_SYM
	svSYMM(src1, src2, 1, 0, 0, &tmp, false);
	svGEMM(&tmp, src2, alpha, src3, beta, dst, SV_GEMM_B_T);
#else
	svGEMM(src1, src2, 1, 0, 0, &tmp, SV_GEMM_FLAG_B_T);
	svGEMM(src2, &tmp, alpha, src3, beta, dst, 0);
#endif
	SV_FREE_STACK_MAT(tmp);
}
// dst = alpha * src1 * src2 + beta * src3
SURVIVE_LOCAL_ONLY void svGEMM(const SvMat *src1, const SvMat *src2, double alpha, const SvMat *src3, double beta,
							   SvMat *dst, enum svGEMMFlags tABC) {

	int rows1 = (tABC & SV_GEMM_FLAG_A_T) ? src1->cols : src1->rows;
	int cols1 = (tABC & SV_GEMM_FLAG_A_T) ? src1->rows : src1->cols;

	int rows2 = (tABC & SV_GEMM_FLAG_B_T) ? src2->cols : src2->rows;
	int cols2 = (tABC & SV_GEMM_FLAG_B_T) ? src2->rows : src2->cols;

	if (src3) {
		int rows3 = (tABC & SV_GEMM_FLAG_C_T) ? src3->cols : src3->rows;
		int cols3 = (tABC & SV_GEMM_FLAG_C_T) ? src3->rows : src3->cols;
		assert(rows3 == dst->rows);
		assert(cols3 == dst->cols);
	}

	// assert(src3 == 0 || beta != 0);
	assert(cols1 == rows2);
	assert(rows1 == dst->rows);
	assert(cols2 == dst->cols);

	lapack_int lda = src1->cols;
	lapack_int ldb = src2->cols;

	if (src3)
		svCopy(src3, dst, 0);
	else
		beta = 0;

	assert(SV_RAW_PTR(dst) != SV_RAW_PTR(src1));
	assert(SV_RAW_PTR(dst) != SV_RAW_PTR(src2));
	assert(dst->cols > 0);
	cblas_gemm(CblasRowMajor, (tABC & SV_GEMM_FLAG_A_T) ? CblasTrans : CblasNoTrans,
			   (tABC & SV_GEMM_FLAG_B_T) ? CblasTrans : CblasNoTrans, dst->rows, dst->cols, cols1, alpha,
			   SV_RAW_PTR(src1), lda, SV_RAW_PTR(src2), ldb, beta, SV_RAW_PTR(dst), dst->cols);
}

// dst = scale * src ^ t * src     iff order == 1
// dst = scale *     src * src ^ t iff order == 0
SURVIVE_LOCAL_ONLY void svMulTransposed(const SvMat *src, SvMat *dst, int order, const SvMat *delta, double scale) {
	lapack_int rows = src->rows;
	lapack_int cols = src->cols;

	lapack_int drows = order == 0 ? dst->rows : dst->cols;
	assert(drows == dst->cols);
	assert(order == 1 ? (dst->cols == src->cols) : (dst->cols == src->rows));
	assert(delta == 0 && "This isn't implemented yet");
	double beta = 0;

	bool isAT = order == 1;
	bool isBT = !isAT;

	lapack_int dstCols = dst->cols;
	assert(dstCols > 0);
	cblas_gemm(CblasRowMajor, isAT ? CblasTrans : CblasNoTrans, isBT ? CblasTrans : CblasNoTrans, dst->rows, dst->cols,
			   order == 1 ? src->rows : src->cols, scale, SV_RAW_PTR(src), src->cols, SV_RAW_PTR(src), src->cols, beta,
			   SV_RAW_PTR(dst), dstCols);
}

/* IEEE754 constants and macros */
#define SV_TOGGLE_FLT(x) ((x) ^ ((int)(x) < 0 ? 0x7fffffff : 0))
#define SV_TOGGLE_DBL(x) ((x) ^ ((int64)(x) < 0 ? SV_BIG_INT(0x7fffffffffffffff) : 0))

#define SV_DbgAssert assert

#define SV_CREATE_MAT_HEADER_ALLOCA(stack_mat, rows, cols)                                                             \
	SvMat *stack_mat = svInitMatHeader(SV_MATRIX_ALLOC(sizeof(SvMat)), rows, cols);

#define SV_CREATE_MAT_ALLOCA(stack_mat, height, width)                                                                 \
	SV_CREATE_MAT_HEADER_ALLOCA(stack_mat, height, width);                                                             \
	(stack_mat)->data = SV_MATRIX_ALLOC(mat_size_bytes(stack_mat));

#define SV_MAT_ALLOCA_FREE(stack_mat)                                                                                \
	{                                                                                \
		SV_MATRIX_FREE(stack_mat->data);                                                                               \
		SV_MATRIX_FREE(stack_mat);                                                                                \
	}

#define CREATE_SV_STACK_MAT(name, rows, cols, type)                                                                    \
	FLT *_##name = alloca(rows * cols * sizeof(FLT));                                                                  \
	SvMat name = svMat(rows, cols, _##name);

static lapack_int LAPACKE_gesvd_static_alloc(int matrix_layout, char jobu, char jobvt, lapack_int m, lapack_int n,
											 FLT *a, lapack_int lda, FLT *s, FLT *u, lapack_int ldu, FLT *vt,
											 lapack_int ldvt, FLT *superb) {
	lapack_int info = 0;
	lapack_int lwork = -1;
	FLT *work = NULL;
	FLT work_query;
	lapack_int i;
	if (matrix_layout != LAPACK_COL_MAJOR && matrix_layout != LAPACK_ROW_MAJOR) {
		LAPACKE_xerbla("LAPACKE_dgesvd", -1);
		return -1;
	}

	/* Query optimal working array(s) size */
	info = LAPACKE_gesvd_work(matrix_layout, jobu, jobvt, m, n, a, lda, s, u, ldu, vt, ldvt, &work_query, lwork);
	if (info != 0) {
		goto exit_level_0;
	}
	lwork = (lapack_int)work_query;
	/* Allocate memory for work arrays */
	work = (FLT *)SV_MATRIX_ALLOC(sizeof(FLT) * lwork);
	memset(work, 0, sizeof(FLT) * lwork);

	if (work == NULL) {
		info = LAPACK_WORK_MEMORY_ERROR;
		goto exit_level_0;
	}
	/* Call middle-level interface */
	info = LAPACKE_gesvd_work(matrix_layout, jobu, jobvt, m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork);
	/* Backup significant data from working array(s) */
	for (i = 0; i < MIN(m, n) - 1; i++) {
		superb[i] = work[i + 1];
	}
	SV_MATRIX_FREE(work);

exit_level_0:
	if (info == LAPACK_WORK_MEMORY_ERROR) {
		LAPACKE_xerbla("LAPACKE_dgesvd", info);
	}
	return info;
}

static inline lapack_int LAPACKE_getri_static_alloc(int matrix_layout, lapack_int n, FLT *a, lapack_int lda,
													const lapack_int *ipiv) {
	lapack_int info = 0;
	lapack_int lwork = -1;
	FLT *work = NULL;
	FLT work_query;
	if (matrix_layout != LAPACK_COL_MAJOR && matrix_layout != LAPACK_ROW_MAJOR) {
		LAPACKE_xerbla("LAPACKE_dgetri", -1);
		return -1;
	}
	/* Query optimal working array(s) size */
	info = LAPACKE_getri_work(matrix_layout, n, a, lda, ipiv, &work_query, lwork);
	if (info != 0) {
		goto exit_level_0;
	}
	lwork = (lapack_int)work_query;
	/* Allocate memory for work arrays */
	work = (FLT *)alloca(sizeof(FLT) * lwork);
	memset(work, 0, sizeof(FLT) * lwork);
	if (work == NULL) {
		info = LAPACK_WORK_MEMORY_ERROR;
		goto exit_level_0;
	}
	/* Call middle-level interface */
	info = LAPACKE_getri_work(matrix_layout, n, a, lda, ipiv, work, lwork);
	/* Release memory and exit */

exit_level_0:
	if (info == LAPACK_WORK_MEMORY_ERROR) {
		LAPACKE_xerbla("LAPACKE_dgetri", info);
	}
	return info;
}

SURVIVE_LOCAL_ONLY double svInvert(const SvMat *srcarr, SvMat *dstarr, enum svInvertMethod method) {
	lapack_int inf;
	lapack_int rows = srcarr->rows;
	lapack_int cols = srcarr->cols;
	lapack_int lda = srcarr->cols;

	svCopy(srcarr, dstarr, 0);
	FLT *a = SV_RAW_PTR(dstarr);

#ifdef DEBUG_PRINT
	printf("a: \n");
	print_mat(srcarr);
#endif
	if (method == SV_INVERT_METHOD_LU) {
		lapack_int *ipiv = SV_MATRIX_ALLOC(sizeof(lapack_int) * MIN(srcarr->rows, srcarr->cols));

		lapack_int lda_t = MAX(1, rows);

		FLT *a_t = (FLT *)SV_MATRIX_ALLOC(sizeof(FLT) * lda_t * MAX(1, cols));
		LAPACKE_ge_trans(LAPACK_ROW_MAJOR, rows, cols, a, lda, a_t, lda_t);

		inf = LAPACKE_getrf(LAPACK_COL_MAJOR, rows, cols, a_t, lda, ipiv);
		assert(inf == 0);

		inf = LAPACKE_getri_static_alloc(LAPACK_COL_MAJOR, rows, a_t, lda, ipiv);
		assert(inf >= 0);

		LAPACKE_ge_trans(LAPACK_COL_MAJOR, rows, cols, a_t, lda, a, lda_t);

		if (inf > 0) {
			printf("Warning: Singular matrix: \n");
			// print_mat(srcarr);
		}

		SV_MATRIX_FREE(a_t);
		SV_MATRIX_FREE(ipiv);
		// free(ipiv);

	} else if (method == DECOMP_SVD) {
		// TODO: There is no way this needs this many allocations,
		// but in my defense I was very tired when I wrote this code
		SV_CREATE_STACK_MAT(w, 1, MIN(dstarr->rows, dstarr->cols));
		SV_CREATE_STACK_MAT(u, dstarr->cols, dstarr->cols);
		SV_CREATE_STACK_MAT(v, dstarr->rows, dstarr->rows);
		SV_CREATE_STACK_MAT(um, w.cols, w.cols);

		svSVD(dstarr, &w, &u, &v, 0);

		svSetZero(&um);
		for (int i = 0; i < w.cols; i++) {
			if (_w[i] != 0.0)
				svMatrixSet(&um, i, i, 1. / (_w)[i]);
		}

		SvMat *tmp = svCreateMat(dstarr->cols, dstarr->rows);
		svGEMM(&v, &um, 1, 0, 0, tmp, SV_GEMM_FLAG_A_T);
		svGEMM(tmp, &u, 1, 0, 0, dstarr, SV_GEMM_FLAG_B_T);

		svReleaseMat(&tmp);

		SV_FREE_STACK_MAT(um);
		SV_FREE_STACK_MAT(v);
		SV_FREE_STACK_MAT(u);
		SV_FREE_STACK_MAT(w);
	} else {
		assert(0 && "Bad argument");
		return -1;
	}

	return 0;
}

#define SV_CLONE_MAT_ALLOCA(stack_mat, mat)                                                                            \
	SV_CREATE_MAT_ALLOCA(stack_mat, mat->rows, mat->cols)                                                              \
	svCopy(mat, stack_mat, 0);

static int svSolve_LU(const SvMat *Aarr, const SvMat *Barr, SvMat *xarr) {
	lapack_int inf;
	lapack_int arows = Aarr->rows;
	lapack_int acols = Aarr->cols;
	lapack_int xcols = Barr->cols;
	lapack_int xrows = Barr->rows;
	lapack_int lda = acols; // Aarr->step / sizeof(double);

	assert(Aarr->cols == xarr->rows);
	assert(Barr->rows == Aarr->rows);
	assert(xarr->cols == Barr->cols);

	svCopy(Barr, xarr, 0);
	FLT *a_ws = SV_MATRIX_ALLOC(mat_size_bytes(Aarr));
	memcpy(a_ws, SV_RAW_PTR(Aarr), mat_size_bytes(Aarr));

	lapack_int brows = xarr->rows;
	lapack_int bcols = xarr->cols;
	lapack_int ldb = bcols; // Barr->step / sizeof(double);

	lapack_int *ipiv = SV_MATRIX_ALLOC(sizeof(lapack_int) * MIN(Aarr->rows, Aarr->cols));

	inf = LAPACKE_getrf(LAPACK_ROW_MAJOR, arows, acols, (a_ws), lda, ipiv);
	assert(inf >= 0);
	if (inf > 0) {
		printf("Warning: Singular matrix: \n");
		// print_mat(a_ws);
	}

#ifdef DEBUG_PRINT
	printf("Solve A * x = B:\n");
	// print_mat(a_ws);
	print_mat(Barr);
#endif

	inf = LAPACKE_getrs(LAPACK_ROW_MAJOR, CblasNoTrans, arows, bcols, (a_ws), lda, ipiv, SV_RAW_PTR(xarr), ldb);
	assert(inf == 0);

	SV_MATRIX_FREE(a_ws);
	SV_MATRIX_FREE(ipiv);
	return 0;
}

static inline int svSolve_SVD(const SvMat *Aarr, const SvMat *Barr, SvMat *xarr) {
	lapack_int arows = Aarr->rows;
	lapack_int acols = Aarr->cols;
	lapack_int xcols = Barr->cols;

	bool xLargerThanB = Barr->rows > acols;
	SvMat *xCpy = 0;
	if (xLargerThanB) {
		SV_CLONE_MAT_ALLOCA(xCpyStack, Barr);
		xCpy = xCpyStack;
	} else {
		xCpy = xarr;
		memcpy(SV_RAW_PTR(xarr), SV_RAW_PTR(Barr), mat_size_bytes(Barr));
	}

	// SvMat *aCpy = svCloneMat(Aarr);
	FLT *aCpy = SV_MATRIX_ALLOC(mat_size_bytes(Aarr));
	memcpy(aCpy, SV_RAW_PTR(Aarr), mat_size_bytes(Aarr));

	FLT *S = SV_MATRIX_ALLOC(sizeof(FLT) * MIN(arows, acols));
	// FLT *S = malloc(sizeof(FLT) * MIN(arows, acols));
	FLT rcond = -1;
	lapack_int *rank = SV_MATRIX_ALLOC(sizeof(lapack_int) * MIN(arows, acols));
	lapack_int inf =
		LAPACKE_gelss(LAPACK_ROW_MAJOR, arows, acols, xcols, (aCpy), acols, SV_RAW_PTR(xCpy), xcols, S, rcond, rank);
	assert(xarr->rows == acols);
	assert(xarr->cols == xCpy->cols);

	if (xLargerThanB) {
		xCpy->rows = acols;
		svCopy(xCpy, xarr, 0);
		// svReleaseMat(&xCpy);
	}

	SV_MATRIX_FREE(rank);
	SV_MATRIX_FREE(aCpy);
	SV_MATRIX_FREE(S);
	if (xLargerThanB) {
		SV_MAT_ALLOCA_FREE(xCpy);
	}

	assert(inf == 0);
	if (inf != 0)
		return inf;
	return 0;
}

SURVIVE_LOCAL_ONLY int svSolve(const SvMat *Aarr, const SvMat *Barr, SvMat *xarr, enum svInvertMethod method) {
	if (method == SV_INVERT_METHOD_LU) {
		return svSolve_LU(Aarr, Barr, xarr);
	} else if (method == SV_INVERT_METHOD_SVD || method == SV_INVERT_METHOD_QR) {
		return svSolve_SVD(Aarr, Barr, xarr);
	} else {
		assert("Unknown method to solve" && 0);
	}
	return -1;
}

SURVIVE_LOCAL_ONLY void svTranspose(const SvMat *M, SvMat *dst) {
	bool inPlace = M == dst || SV_RAW_PTR(M) == SV_RAW_PTR(dst);
	FLT *src = SV_RAW_PTR(M);

	if (inPlace) {
		src = alloca(mat_size_bytes(M));
		memcpy(src, SV_RAW_PTR(M), mat_size_bytes(M));
	} else {
		assert(M->rows == dst->cols);
		assert(M->cols == dst->rows);
	}

	for (unsigned i = 0; i < M->rows; i++) {
		for (unsigned j = 0; j < M->cols; j++) {
			SV_RAW_PTR(dst)[j * M->rows + i] = src[i * M->cols + j];
		}
	}
}

#if defined(__has_feature)
#if __has_feature(memory_sanitizer)
#define MEMORY_SANITIZER_IGNORE __attribute__((no_sanitize("memory")))
#endif
#endif
#ifndef MEMORY_SANITIZER_IGNORE
#define MEMORY_SANITIZER_IGNORE
#endif

#define CALLOCA(size) memset(alloca(size), 0, size)

SURVIVE_LOCAL_ONLY void svSVD(SvMat *aarr, SvMat *warr, SvMat *uarr, SvMat *varr, enum svSVDFlags flags) {
	char jobu = 'A';
	char jobvt = 'A';

	lapack_int inf;

	if ((flags & SV_SVD_MODIFY_A) == 0) {
		aarr = svCloneMat(aarr);
	}

	if (uarr == 0)
		jobu = 'N';
	if (varr == 0)
		jobvt = 'N';

	lapack_int arows = aarr->rows, acols = aarr->cols;

	FLT *pw = warr ? SV_RAW_PTR(warr) : (FLT *)CALLOCA(sizeof(FLT) * arows * acols);
	FLT *pu = uarr ? SV_RAW_PTR(uarr) : (FLT *)CALLOCA(sizeof(FLT) * arows * arows);
	FLT *pv = varr ? SV_RAW_PTR(varr) : (FLT *)CALLOCA(sizeof(FLT) * acols * acols);

	lapack_int ulda = uarr ? uarr->cols : acols;
	lapack_int plda = varr ? varr->cols : acols;

	FLT *superb = CALLOCA(sizeof(FLT) * MIN(arows, acols));
	inf = LAPACKE_gesvd_static_alloc(LAPACK_ROW_MAJOR, jobu, jobvt, arows, acols, SV_RAW_PTR(aarr), acols, pw, pu, ulda,
									 pv, plda, superb);

	switch (inf) {
	case -6:
		assert(false && "matrix has NaNs");
		break;
	case 0:
		break;
	default:
		assert(inf == 0);
	}

	if (uarr && (flags & SV_SVD_U_T)) {
		svTranspose(uarr, uarr);
	}

	if (varr && (flags & SV_SVD_V_T) == 0) {
		svTranspose(varr, varr);
	}

	if ((flags & SV_SVD_MODIFY_A) == 0) {
		svReleaseMat(&aarr);
	}
}

SURVIVE_LOCAL_ONLY double svDet(const SvMat *M) {
	assert(M->rows == M->cols);
	assert(M->rows <= 3 && "svDet unimplemented for matrices >3");

	FLT *m = SV_RAW_PTR(M);

	switch (M->rows) {
	case 1:
		return m[0];
	case 2: {
		return m[0] * m[3] - m[1] * m[2];
	}
	case 3: {
		FLT m00 = m[0], m01 = m[1], m02 = m[2], m10 = m[3], m11 = m[4], m12 = m[5], m20 = m[6], m21 = m[7], m22 = m[8];

		return m00 * (m11 * m22 - m12 * m21) - m01 * (m10 * m22 - m12 * m20) + m02 * (m10 * m21 - m11 * m20);
	}
	default:
		abort();
	}
}

ext: c

max_stars_repo_name: jdavidberger/libsurvive

max_issues_repo_name: jdavidberger/libsurvive

max_forks_repo_name: jdavidberger/libsurvive


text:
/**
 *
 * @file qwrapper_zgemv.c
 *
 *  PLASMA core_blas quark wrapper
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Mark Gates
 * @author Mathieu Faverge
 * @date 2010-11-15
 * @precisions normal z -> c d s
 *
 **/
#include <cblas.h>
#include "common.h"

/***************************************************************************//**
 *
 **/
void QUARK_CORE_zgemv(Quark *quark, Quark_Task_Flags *task_flags,
                      PLASMA_enum trans, int m, int n,
                      PLASMA_Complex64_t alpha, const PLASMA_Complex64_t *A, int lda,
                                                const PLASMA_Complex64_t *x, int incx,
                      PLASMA_Complex64_t beta,        PLASMA_Complex64_t *y, int incy)
{
    DAG_CORE_GEMV;
    QUARK_Insert_Task(quark, CORE_zgemv_quark, task_flags,
        sizeof(PLASMA_enum),             &trans,  VALUE,
        sizeof(int),                     &m,      VALUE,
        sizeof(int),                     &n,      VALUE,
        sizeof(PLASMA_Complex64_t),      &alpha,  VALUE,
        sizeof(PLASMA_Complex64_t)*m*n,  A,               INPUT,
        sizeof(int),                     &lda,    VALUE,
        sizeof(PLASMA_Complex64_t)*n,    x,               INPUT,
        sizeof(int),                     &incx,   VALUE,
        sizeof(PLASMA_Complex64_t),      &beta,   VALUE,
        sizeof(PLASMA_Complex64_t)*m,    y,               INOUT,
        sizeof(int),                     &incy,   VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_zgemv_quark = PCORE_zgemv_quark
#define CORE_zgemv_quark PCORE_zgemv_quark
#endif
void CORE_zgemv_quark(Quark *quark)
{
    PLASMA_enum trans;
    int m, n, lda, incx, incy;
    PLASMA_Complex64_t alpha, beta;
    const PLASMA_Complex64_t *A, *x;
    PLASMA_Complex64_t *y;

    quark_unpack_args_11( quark, trans, m, n, alpha, A, lda, x, incx, beta, y, incy );
    cblas_zgemv(
        CblasColMajor,
        (CBLAS_TRANSPOSE)trans,
        m, n,
        CBLAS_SADDR(alpha), A, lda,
                            x, incx,
        CBLAS_SADDR(beta),  y, incy);
}


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdbool.h>
#include <math.h>
#include <float.h>
#include "parmt_postProcess.h"
#include "beachball.h"
#ifdef PARMT_USE_INTEL
#include <mkl_cblas.h>
#else
#include <cblas.h>
#endif
#include "compearth.h"

static void argsort3(const double *__restrict__ x, int *__restrict__ iperm);
static void sortUL(const double *__restrict__ lam,
                   const double *__restrict__ U,
                   double *__restrict__ lamSort,
                   double *__restrict__ Usort);
static void eigenvector2PrincipalAxis(const double eig,
                                      const double *__restrict__ ev,
                                      double *__restrict__ paxis);
static void vecstrd(const double str, const double dip, double r[3]);
static double dot3(const double *__restrict__ x,
                   const double *__restrict__ y);
static void gemv3(const double *__restrict__ A, const double *__restrict__ x,
                  double *__restrict__ y);
static void crossp(const double *__restrict__ u,
                   const double *__restrict__ v,
                   double *__restrict__ n);
static void strdvec(double r[3], double *str, double *dip);

int postprocess_tnp2beachballPolarity(const int nxp,
                                      const double xc, const double yc,
                                      const double rad,
                                      const double *__restrict__ pAxis,
                                      const double *__restrict__ nAxis,
                                      const double *__restrict__ tAxis,
                                      double *__restrict__ xw1,
                                      double *__restrict__ yw1,
                                      int8_t *__restrict__ pn1)
{
    double ax1[3] __attribute__ ((aligned(64)));
    double ax2[3] __attribute__ ((aligned(64)));
    double ax3[3] __attribute__ ((aligned(64)));
    double rm[9] __attribute__ ((aligned(64)));
    double rmt[9] __attribute__ ((aligned(64)));
    double tm[9] __attribute__ ((aligned(64)));
    double pt[3] __attribute__ ((aligned(64)));
    double ptnew[3] __attribute__ ((aligned(64)));
    double ampMax, az, baz, baz1, bdp, bdp1, bev, bevI, dp,
           dx, dxc, dy, dyc, evI,
           fclvd, fIso, paz, paz1, pdp, pdp1, pev, pevI, plusMin,
           rad2, rpt, rpt2,
           taz, taz1, tdp, tdp1, tev, tevI, theta, x, x0, y, y0;
    int ix, iy, k, nx, ny;
    const double degradi = 180.0/M_PI;
    const double third = 1.0/3.0;
    const double sqrt2 = 1.4142135623730951;
    const double sqrt2i = 1.0/sqrt2;
    taz1 = tAxis[0];
    tdp1 = tAxis[1];
    tevI = tAxis[2];
    baz1 = nAxis[0];
    bdp1 = nAxis[1];
    bevI = nAxis[2];
    paz1 = pAxis[0];
    pdp1 = pAxis[1];
    pevI = pAxis[2];
    // Avoid overwriting memory
    taz = taz1;
    tdp = tdp1;
    baz = baz1;
    bdp = bdp1;
    paz = paz1;
    pdp = pdp1;
    // Step 1: set up axes:
    //  - D axis ('dominant' 1-axis - with largest absolute eigenvalue)
    //  - M axis ('minor'    3-axis)
    //  - B axis (           2-axis)
    evI = (tevI + bevI + pevI)*third;
    tev = tevI - evI;
    bev = bevI - evI;
    pev = pevI - evI;
    vecstrd(paz, pdp, ax1);
    vecstrd(taz, tdp, ax3);
    fclvd = fabs(bev/pev);
    fIso = evI/pev;
    if (fabs(tev) > fabs(pev))
    {
        vecstrd(taz, tdp, ax1);
        vecstrd(paz, pdp, ax3);
        fclvd = fabs(bev/tev);
        fIso = evI/tev;
    }
    //if (fabs(evI) > .03)
    //{
    //    write(6,901) evI
    //}
    crossp(ax1, ax3, ax2);
    strdvec(ax2, &baz, &bdp);
    cliffsNodes_mtensor(taz, tdp, paz, pdp, tevI, bevI, pevI,
                        tm, rm, rmt, &ampMax);
    nx = nxp;
    ny = nxp;
    memset(xw1, 0, (size_t) (nx*ny)*sizeof(double));
    memset(yw1, 0, (size_t) (nx*ny)*sizeof(double));
    memset(pn1, 0, (size_t) (nx*ny)*sizeof(int8_t));
    dx = 2.0*rad/(double) (nx - 1);
    dy = 2.0*rad/(double) (ny - 1);
    x0 = xc - rad;
    y0 = yc - rad;
    rad2 = rad*rad; // saves a sqrt computation when comparing distance
    // loop on pixel grid
    for (iy=0; iy<ny; iy++)
    {
        for (ix=0; ix<nx; ix++)
        {
            x = x0 + (double) ix*dx;
            y = y0 + (double) iy*dy;
            k = iy*nx + ix;
            pn1[k] = 0;
            xw1[k] = x;
            yw1[k] = y;
            // require (x,y) be in the focal sphere
            dxc = x - xc;
            dyc = y - yc;
            rpt2 = dxc*dxc + dyc*dyc;
            if (rpt2 < rad2)
            {
                rpt = sqrt(rpt2);
                // compute the azimuth
                theta = atan2(y - yc, x - xc);  // measure from x axis
                az = (M_PI_2 - theta)*degradi;  // az is measured from north
                if (az < 0.0){az = az + 360.0;} // convention is [0,360]
                // compute the dip
                dp = (M_PI_2 - 2.0*asin(rpt*sqrt2i))*degradi;
                vecstrd(az, dp, pt);
                gemv3(tm, pt, ptnew);
                plusMin = dot3(pt, ptnew);
                pn1[k] =-1;
                if (plusMin > 100.0*DBL_EPSILON*ampMax){pn1[k] = 1;}
            } // end check on location
        } // loop on x
    } // loop on y
    return 0;
}
/*!
 * @brief Converts a Tape and Tape 2015 tensor to tension, null, and
 *        plunge eigenvectors in USE coordinates for plotting with
 *        cliffsnodes
 * @param[in] beta    colatitude (radians)
 * @param[in] gamma   longitude (radians)
 * @param[in] kappa   strike angle (radians)
 * @param[in] sigma   slip angle (radians)
 * @param[in] theta   dip angle (radians)
 *
 * @param[out] pAxis  (azimuth, plunge, eigenvalue) of pressure axis [3]
 * @param[out] nAxis  (azimuth, plunge, eigenvalue) of null axis [3]
 * @param[out] tAxis  (azimuth, plunge, eigenvalue) of tension axis [3]
 *
 * @result 0 indicate success
 *
 * @author Ben Baker
 *
 * @copyright ISTI distributed under Apache 2
 *
 */
int postprocess_tt2tnp(const double beta, const double gamma,
                       const double kappa, const double sigma,
                       const double theta,
                       double *__restrict__ pAxis,
                       double *__restrict__ nAxis,
                       double *__restrict__ tAxis)
{
    const double R[9] = {
             0.7071067811865476, 0.0,                -0.7071067811865476,
            -0.4082482904638631, 0.8164965809277261, -0.4082482904638631,
             0.5773502691896258, 0.5773502691896258,  0.5773502691896258};
    const double Yrot[9] = {
             0.7071067811865475, 0.0, 0.7071067811865475,
             0.0, 1.0, 0.0,
            -0.7071067811865475, 0.0, 0.7071067811865475};
    // Convert NWU to USE
    const double Rot[9] = {0, -1, 0, 0, 0, -1, 1, 0, 0};
    //const double RotInv[9] = {0, 0, 1, -1, 0, 0, 0, -1, 0};
    double Zsigma[9], Xtheta[9], Zkappa[9], ZX[9], V[9],
           U[9], Uuse[9], Usort[9], Uw[9],
           lambda[3], v3[3], lamSort[3],
           cosb, cosg, sinb, sing, nKappaDeg, sigmaDeg, thetaDeg;
    const double toDeg = 180.0/M_PI;
    // Compute the eigenvalues from Eqn 7
    sinb = sin(beta);
    sing = sin(gamma);
    cosb = cos(beta);
    cosg = cos(gamma);
    v3[0] = sinb*cosg;
    v3[1] = sinb*sing;
    v3[2] = cosb;
    // Compute R*lam
    lambda[0] = R[0]*v3[0] + R[3]*v3[1] + R[6]*v3[2];
    lambda[1] = R[1]*v3[0] + R[4]*v3[1] + R[7]*v3[2];
    lambda[2] = R[2]*v3[0] + R[5]*v3[1] + R[8]*v3[2];
    // Eqn 9-10
    sigmaDeg  = sigma*toDeg;
    thetaDeg  = theta*toDeg;
    nKappaDeg =-kappa*toDeg;
    compearth_eulerUtil_rotmat(1, &sigmaDeg, 3, Zsigma);
    compearth_eulerUtil_rotmat(1, &thetaDeg, 1, Xtheta);
    compearth_eulerUtil_rotmat(1, &nKappaDeg, 3, Zkappa);
    //compearth_eulerUtil_rotmatGen(1,-45.0, 2, Yrot);
    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
                3, 3, 3, 1.0, Zkappa, 3, Xtheta, 3, 0.0, ZX, 3);
    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
                3, 3, 3, 1.0, ZX, 3, Zsigma, 3, 0.0, V, 3);
    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
                3, 3, 3, 1.0, V, 3, Yrot, 3, 0.0, U, 3);
    // M could be computed by computing M = U*Lambda*inv(U) hence
    // U corresponds to the eigenvectors and lambda the scaling.
    // However, the current basis is NWU which must be converted to USE.
    // If performing a proper rotation of M to M_{use} we would compute:
    //   R M inv(R) = R U Lambda inv(U) inv(R)
    // thus what we must be do is transform the eigenvectors into the
    // apropriate frame.  Intuitively it is the rotation of each
    // eigenvector (column) of the U matrix which is indeed the matrix
    // matrix multiplication of R*U.
    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
                3, 3, 3, 1.0, Rot, 3, U, 3, 0.0, Uuse, 3);
    // If I want to finish this I could compute:
    //   R U Lambda inv(U) inv(R) = R U Lambda U^T R^T
/*
    double lamMat[9], Temp[9], M[9], UtRt[9];
    memset(lamMat, 0, 9*sizeof(double));
    lamMat[0] = lambda[0]; lamMat[4] = lambda[1]; lamMat[8] = lambda[2];
    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
                3, 3, 3, 1.0, Uuse, 3, lamMat, 3, 0.0, Temp, 3);
    cblas_dgemm(CblasColMajor, CblasTrans, CblasTrans,
                3, 3, 3, 1.0, U, 3, Rot, 3, 0.0, UtRt, 3);
    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
                3, 3, 3, 1.0, Temp, 3, UtRt, 3, 0.0, M, 3);
    for (int i=0; i<3; i++)
    {
        for (int j=0; j<3; j++)
        {
            printf("%f ", M[3*j+i]);
        }
        printf("\n");
     }
*/
    sortUL(lambda, Uuse, lamSort, Usort);
//printf("%f %f %f\n", lamSort[0], lamSort[1], lamSort[2]);
    // Compute pressure, null, and tension axes
    eigenvector2PrincipalAxis(lamSort[0], &Usort[0], pAxis);
    eigenvector2PrincipalAxis(lamSort[1], &Usort[3], nAxis);
    eigenvector2PrincipalAxis(lamSort[2], &Usort[6], tAxis);
    // Explosion - simplification (all positive eigenvalues)
    if (lamSort[0] > 0.0)
    {
        pAxis[2] = 1.0;
        nAxis[2] = 1.0;
        tAxis[2] = 1.0;
    }
    // Implosion - simplification (all negative eigenvalues)
    if (lamSort[2] < 0.0)
    {
        pAxis[0] =-1.0;
        nAxis[1] =-1.0;
        tAxis[2] =-1.0;
    }
    //printf("%f %f %f\n", pAxis[0], pAxis[1], pAxis[2]);
    //printf("%f %f %f\n", nAxis[0], nAxis[1], nAxis[2]);
    //printf("%f %f %f\n", tAxis[0], tAxis[1], tAxis[2]);
    return 0;
}

/*!
 * @brief Convert pressure, null, and tension principal axes to azimuth,
 *        plunge, and length.
 */
static void eigenvector2PrincipalAxis(const double eig,
                                      const double *__restrict__ ev,
                                      double *__restrict__ paxis)
{
    const char *fcnm = "cmopad_Eigenvector2PrincipalAxis\0";
    double v[3], az, plunge;
    double twopi = 2.0*M_PI;
    double pi180i = 180.0/M_PI;
    int ierr;
    //------------------------------------------------------------------------//    //
    // Compute azimuth and plunge angles where ev is USE coordinates
    plunge = asin(-ev[0]);        //arcsin(-z/r)
    az     = atan2(ev[2],-ev[1]); //atan(x/y)
    if (plunge <= 0.0)
    {
        plunge =-plunge;
        az = az + M_PI;
    }
    if (az < 0.0){az = az + twopi;}   //Shift to [0,360]
    if (az > twopi){az = az - twopi;} //Shift to [0,360]
    // Convert to degrees
    plunge = plunge*pi180i;
    az = az*pi180i;
    // Copy back result
    paxis[0] = az;     //First value is azimuth (degrees)
    paxis[1] = plunge; //Second value is plunge (degrees)
    paxis[2] = eig;    //Final value is eigenvalue (distance)
    return;
}
//===========================================================================//
static void sortUL(const double *__restrict__ lam,
                   const double *__restrict__ U,
                   double *__restrict__ lamSort,
                   double *__restrict__ Usort)
{
    int perm[3], j;
    argsort3(lam, perm);
    lamSort[0] = lam[perm[0]];
    lamSort[1] = lam[perm[1]];
    lamSort[2] = lam[perm[2]];
    for (j=0; j<3; j++)
    {
        Usort[3*j+0] = U[3*perm[j]+0];
        Usort[3*j+1] = U[3*perm[j]+1];
        Usort[3*j+2] = U[3*perm[j]+2];
    }
    return;
}

//============================================================================//
/*!
 * @brief Permutation to sort length 3 array into ascending order
 */
static void argsort3(const double *__restrict__ x, int *__restrict__ iperm)
{
    const char *fcnm = "__cmopad_argsort3\0";
    int i, temp;
    const int a = 0;
    const int b = 1;
    const int c = 2;
    // Copy
    iperm[a] = a;
    iperm[b] = b;
    iperm[c] = c;
    if (x[iperm[a]] > x[iperm[c]])
    {
        temp = iperm[c];
        iperm[c] = iperm[a];
        iperm[a] = temp;
    }
    if (x[iperm[a]] > x[iperm[b]])
    {
        temp = iperm[b];
        iperm[b] = iperm[a];
        iperm[a] = temp;
    }
    //Now the smallest element is the first one. Just check the 2-nd and 3-rd
    if (x[iperm[b]] > x[iperm[c]])
    {
        temp = iperm[c];
        iperm[c] = iperm[b];
        iperm[b] = temp;
    }
    return;
}
//============================================================================//
/*!
 * @brief Find components of downward pointing unit vector pole from the
 *        strike and dip angles.  This is based on Cliff Frolich's vecstrd.
 *
 * @param[in] str    strike angle of pole (degrees)
 * @param[in] dip    dip angle of pole (degrees)
 *
 * @param[out] r     corresponding unit vector pole (3)
 *
 * @author Ben Baker
 *
 * @copyright BSD
 *
 */
static void vecstrd(const double str, const double dip, double r[3])
{
    double cosdip, cosstr, diprad, sindip, sinstr, strrad;
    const double degrad = M_PI/180.0;
    // convert to radians
    strrad = str*degrad;
    diprad = dip*degrad;
    // compute cosines and sines
    sinstr = sin(strrad);
    cosstr = cos(strrad);
    cosdip = cos(diprad);
    sindip = sin(diprad);
    // compute (x,y,z) from strike and dip for unit radius
    r[0] = sinstr*cosdip;
    r[1] = cosstr*cosdip;
    r[2] =-sindip;
    return;
}
/*!
 * @brief Utility function for computing dot product of length 3 vector
 */
static double dot3(const double *__restrict__ x,
                   const double *__restrict__ y)
{
    double dot;
    dot = x[0]*y[0] + x[1]*y[1] + x[2]*y[2];
    return dot;
}
/*!
 * @brief Utility function for computing y = Ax for 3 x 3 matrix.
 */
static void gemv3(const double *__restrict__ A, const double *__restrict__ x,
                  double *__restrict__ y)
{
    y[0] = A[0]*x[0] + A[3]*x[1] + A[6]*x[2];
    y[1] = A[1]*x[0] + A[4]*x[1] + A[7]*x[2];
    y[2] = A[2]*x[0] + A[5]*x[1] + A[8]*x[2];
    return;
}
/*!
 * @brief Cross product of 2 length 3 vectors u and v.  Returns normal vector n.
 *
 * @param[in] u    first vector in cross product u x v [3]
 * @param[in] v    second vector in cross product u x v [3]
 *
 * @param[out] n   normal vector n = u x v [3]
 *
 * @author Ben Baker
 *
 * @copyright BSD
 *
 */
static void crossp(const double *__restrict__ u,
                   const double *__restrict__ v,
                   double *__restrict__ n)
{
    n[0] = u[1]*v[2] - u[2]*v[1];
    n[1] = u[2]*v[0] - u[0]*v[2];
    n[2] = u[0]*v[1] - u[1]*v[0];
    return;
}
/*!
 * @brief Finds the strike and dip of downward pointing unit vector.  This
 *        is based on Cliff Frolich's strdvec.
 *
 * @param[in,out] r      on input contains the downward pointing unit
 *                       vector.
 *                       if r[2] is > 0 then the unit vector is pointing
 *                       updward.  in this case then on output r3 will
 *                       be negated.
 *
 * @param[out] str       strike angle (degrees) \f$ s.t. \phi \in [0,360] \f$.
 * @param[out] dip       dip angle (degrees)
 *
 * @author Ben Baker
 *
 * @copyright BSD
 *
 */
static void strdvec(double r[3], double *str, double *dip)
{
    double rs;
    const double degradi = 180.0/M_PI;
    if (r[2] > 0.0)
    {
        r[0] =-r[0];
        r[1] =-r[1];
        r[2] =-r[2];
    }
    *str = atan2(r[0], r[1])*degradi;
    if (*str < 0.){*str = *str + 360.;}
    rs = sqrt(r[0]*r[0] + r[1]*r[1]);
    *dip=atan2(-r[2], rs)*degradi;
    return;
}


text:
/* integration/qk.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

#include <config.h>
#include <float.h>
#include <math.h>
#include <gsl/gsl_integration.h>
#include "err.c"

void
gsl_integration_qk (const int n,
                    const double xgk[], const double wg[], const double wgk[],
                    double fv1[], double fv2[],
                    const gsl_function * f, double a, double b,
                    double *result, double *abserr,
                    double *resabs, double *resasc)
{

  const double center = 0.5 * (a + b);
  const double half_length = 0.5 * (b - a);
  const double abs_half_length = fabs (half_length);
  const double f_center = GSL_FN_EVAL (f, center);

  double result_gauss = 0;
  double result_kronrod = f_center * wgk[n - 1];

  double result_abs = fabs (result_kronrod);
  double result_asc = 0;
  double mean = 0, err = 0;

  int j;

  if (n % 2 == 0)
    {
      result_gauss = f_center * wg[n / 2 - 1];
    }

  for (j = 0; j < (n - 1) / 2; j++)
    {
      const int jtw = j * 2 + 1;	/* j=1,2,3 jtw=2,4,6 */
      const double abscissa = half_length * xgk[jtw];
      const double fval1 = GSL_FN_EVAL (f, center - abscissa);
      const double fval2 = GSL_FN_EVAL (f, center + abscissa);
      const double fsum = fval1 + fval2;
      fv1[jtw] = fval1;
      fv2[jtw] = fval2;
      result_gauss += wg[j] * fsum;
      result_kronrod += wgk[jtw] * fsum;
      result_abs += wgk[jtw] * (fabs (fval1) + fabs (fval2));
    }

  for (j = 0; j < n / 2; j++)
    {
      int jtwm1 = j * 2;
      const double abscissa = half_length * xgk[jtwm1];
      const double fval1 = GSL_FN_EVAL (f, center - abscissa);
      const double fval2 = GSL_FN_EVAL (f, center + abscissa);
      fv1[jtwm1] = fval1;
      fv2[jtwm1] = fval2;
      result_kronrod += wgk[jtwm1] * (fval1 + fval2);
      result_abs += wgk[jtwm1] * (fabs (fval1) + fabs (fval2));
    };

  mean = result_kronrod * 0.5;

  result_asc = wgk[n - 1] * fabs (f_center - mean);

  for (j = 0; j < n - 1; j++)
    {
      result_asc += wgk[j] * (fabs (fv1[j] - mean) + fabs (fv2[j] - mean));
    }

  /* scale by the width of the integration region */

  err = (result_kronrod - result_gauss) * half_length;

  result_kronrod *= half_length;
  result_abs *= abs_half_length;
  result_asc *= abs_half_length;

  *result = result_kronrod;
  *resabs = result_abs;
  *resasc = result_asc;
  *abserr = rescale_error (err, result_abs, result_asc);

}


text:
/* Implementation of the interface in orbital_state_vector.h.  */

#include <assert.h>
#include <math.h>
#include <stdlib.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_odeiv.h>

#include "basic_types.h"
#include "earth_constants.h"
#include "orbital_state_vector.h"

OrbitalStateVector *
orbital_state_vector_new (double x, double y, double z, double vx, double vy,
			  double vz)
{
  OrbitalStateVector *self = malloc (sizeof (OrbitalStateVector));

  self->position = vector_new (x, y, z);
  self->velocity = vector_new (vx, vy, vz);

  return self;
}

OrbitalStateVector *
orbital_state_vector_copy (OrbitalStateVector *a)
{
  OrbitalStateVector *new = malloc (sizeof (OrbitalStateVector));

  new->position = vector_copy (a->position);
  new->velocity = vector_copy (a->velocity);

  return new;
}

void
orbital_state_vector_get_keplerian_elements
(OrbitalStateVector *self, double *a, double *e, double *i,
 double *capital_omega, double *lower_case_omega, double *nu, double *E,
 double *M)
{
  /* Shorthand.  */
  static const double gm = EARTH_GRAVITATIONAL_CONSTANT;

  /* Fundamental vectors as described in "Satellite Geodesy".  */
  Vector *h = vector_cross (self->position, self->velocity);
  Vector *X = vector_new (1.0, 0.0, 0.0);
  Vector *Z = vector_new (0.0, 0.0, 1.0);
  Vector *n = vector_cross (Z, h);
  double r_dot_v = vector_dot (self->position, self->velocity);
  double r_mag = vector_magnitude (self->position);
  double v_mag = vector_magnitude (self->velocity);
  Vector *tmp1 = vector_copy (self->position);
  Vector *tmp2 = vector_copy (self->velocity);
  /* Fundamental vector 'e' ('_fv' is appended to dodge name clash).  */
  Vector *e_fv;
  double h_mag;
  double n_mag;

  vector_multiply (tmp1, pow (v_mag, 2) - gm / r_mag);
  vector_multiply (tmp2, r_dot_v);
  vector_subtract (tmp1, tmp2);
  vector_multiply (tmp1, (1.0 / gm));
  e_fv = vector_copy (tmp1);

  /* Compute keplerian elements.  */
  *e = vector_magnitude (e_fv);
  h_mag = vector_magnitude (h);
  *a = (pow (h_mag, 2) / gm) / (1 - pow (*e, 2));
  *i = acos (vector_dot (h, Z) / h_mag);
  n_mag = vector_magnitude (n);
  *capital_omega = acos (vector_dot (n, X) / n_mag);
  *lower_case_omega = acos (vector_dot (n, e_fv) / (n_mag * (*e)));
  *nu = acos (vector_dot (e_fv, self->position)
	      / ((*e) * vector_magnitude (self->position)));
  *E = acos ((*e + cos (*nu)) / (1.0 + (*e) * cos (*nu)));
  *M = *E - (*e) * sin (*E);

  /* Free fundamental and temporary vectors used in computation.  */
  free (h);
  free (Z);
  free (X);
  free (n);
  free (tmp1);
  free (tmp2);
  free (e_fv);
}

double
orbital_state_vector_height_above_ellipsoid (OrbitalStateVector *self)
{
  double r = vector_magnitude (self->position);

  /* Compute height of satellite above the WGS84 ellipsoid.  */
  static Vector z_unit = {0.0, 0.0, 1.0};
  /* Angle between position vector and equatorial plane.  */
  double eq_angle = M_PI / 2 - vector_angle (self->position, &z_unit);
  double ae = EARTH_SEMIMAJOR_AXIS;
  /* Earth semiminor axis.  */
  double be = sqrt (pow (ae, 2) * (1 - pow (EARTH_ECCENTRICITY, 2)));
  /* Radius of ellipsoid at point under satellite.  */
  double ellipsoid_radius = sqrt (pow (ae * cos (eq_angle), 2)
				  + pow (be * sin (eq_angle), 2));

  return r - ellipsoid_radius;
}

/* Compute derivatives f from functions y(t), t, and other
   miscellaneous parameters.  Return GSL_SUCCESS after successful
   computation.  This code mostly follows the example in the GNU
   Scientific Library documentation, which uses a y array for the
   function values.  The y[1] value is the only one that corresponds
   to the y coordinate, and the y[4] vale the only one that
   corresponds to the y velocity.  */
static int
func (double t, const double y[], double f[], void *params)
{
  /* Shorthand for some earth constants.  */
  static const double gm = EARTH_GRAVITATIONAL_CONSTANT;
  static const double ae = EARTH_SEMIMAJOR_AXIS;

  /* Create position and velocity vectors.  We use static vectors to
     avoid allocation overhead.  */
  static Vector p_st;
  static Vector v_st;
  static Vector *p = &p_st;
  static Vector *v = &v_st;

  double r;
  double j2;
  double a_j2_x, a_j2_y, a_j2_z;
  double ksx, ksy, ksz;

  vector_set (p, y[0], y[1], y[2]);
  vector_set (v, y[3], y[4], y[5]);

  /* Current range to satellife from center of earth.  */
  r = vector_magnitude (p);
  double r3 = r*r*r;

  /* The so-called first zonal harmonic is used to account for the
     effects of earth oblatedness (see below for references).  */
  j2 = 1082.63e-6;	/* First zonal harmonic coefficient.  */
  a_j2_x = (gm * p->x / r3) * j2 * (3.0 / 2) * (ae*ae / (r*r))
    * (5.0 * (p->z * p->z) / (r*r) - 1.0);
  a_j2_y = (p->y / p->x) * a_j2_x;
  a_j2_z = (-gm * p->z / r3) * j2 * (3.0 / 2)
    * (ae*ae / (r*r)) * (3.0 - 5.0 * (p->z * p->z) / (r*r));

  /* Total perturbations.  */
  ksx = a_j2_x;
  ksy = a_j2_y;
  ksz = a_j2_z;

  f[0] = v->x;
  f[1] = v->y;
  f[2] = v->z;
  f[3] = ksx - gm * (p->x / r3);
  f[4] = ksy - gm * (p->y / r3);
  f[5] = ksz - gm * (p->z / r3);

  return GSL_SUCCESS;
}

void
orbital_state_vector_propagate (OrbitalStateVector *self, double time)
{
  /* We need to solve the following six dimensional system of ordinary
     differential equations:

         dx/dt = vx
	 dy/dt = vy
	 dz/dt = vz

	 dvx/dt = ksx - (GM / r^3) * x
	 dvy/dt = ksy - (GM / r^3) * y
	 dvz/dt = ksz - (GM / r^3) * z

     where

         r = sqrt (x^2 + y^2 + z^2)

     with initial conditions given by the current state.  The values
     ksx, ksy, and ksz represent the perturbations due to earth
     oblatedness as described in section 3.2.2.3 respectively of
     "Satellite Geodesy", 2nd edition, by Gunter Seeber.  An ordinary
     differential equation solver from the GNU Scientific Library
     (GSL) is used to solve the system.  It is well described in the
     GSL documentation section "Ordinary Differential Equations".  */

  const int dimension = 6;
  /* Create some things the ODE solver uses.  */
  const gsl_odeiv_step_type *step_type = gsl_odeiv_step_rkf45;
  gsl_odeiv_step *ode_step = gsl_odeiv_step_alloc (step_type, dimension);
  /* Hold the absolute integration error for each coordinate for each
     step to within this value.  This value was chosen using the
     following sophisticated methodology: Keep tightening the
     tolerance until the answer no longer changes with further
     tightenings, for an short (~2.5 minute) propagation.  */
  double mae = 0.000001;
  gsl_odeiv_control *ode_control = gsl_odeiv_control_y_new (mae, 0.0);
  gsl_odeiv_evolve *ode_evolve = gsl_odeiv_evolve_alloc (dimension);
  gsl_odeiv_system ode_system;
  double t0 = 0.0, t1 = time;	/* Start and end times.  */
 /* Initial guess for step size.  */
  const double initial_step_size = GSL_SIGN (t1) * 1.0;
  double step_size = initial_step_size;
  double *y = malloc (dimension * sizeof (double));
  double t = t0;		/* Current time.  */

  ode_system.function = func;
  ode_system.jacobian = NULL;
  ode_system.dimension = dimension;
  ode_system.params = NULL;

  y[0] = self->position->x;
  y[1] = self->position->y;
  y[2] = self->position->z;
  y[3] = self->velocity->x;
  y[4] = self->velocity->y;
  y[5] = self->velocity->z;
  /* Here is the actual propagation.  */
  while ( fabs (t) < fabs (t1) ) {
    int status = gsl_odeiv_evolve_apply (ode_evolve, ode_control, ode_step,
					 &ode_system, &t, t1, &step_size, y);
    assert (status == GSL_SUCCESS);
  }
  self->position->x = y[0];
  self->position->y = y[1];
  self->position->z = y[2];
  self->velocity->x = y[3];
  self->velocity->y = y[4];
  self->velocity->z = y[5];

  free (y);
  gsl_odeiv_evolve_free (ode_evolve);
  gsl_odeiv_control_free (ode_control);
  gsl_odeiv_step_free (ode_step);
}

void
orbital_state_vector_propagate_interpolate (OrbitalStateVector *self,
					    double time,
					    OrbitalStateVector *other)
{
  /* Not implemented.  */
  assert (FALSE);
}

max_stars_repo_name: glshort/MapReady

max_issues_repo_name: glshort/MapReady

max_forks_repo_name: glshort/MapReady


text:
#ifndef CURVATURE_H_LSXONAFA
#define CURVATURE_H_LSXONAFA

#include <cmath>
#include <gsl/gsl>
#include <type_traits>

namespace sens_loc::math {

/// This function calculates the gaussian curvature for a function given its
/// derivatives.
/// \tparam Real precision of the calculation
/// \param f_u,f_v,f_uu,f_vv,f_uv partial derivatives of the depth values
/// \returns gaussian curvature for these derivatives
/// \sa conversion::depth_to_gaussian_curvature
/// \note check https://en.wikipedia.org/wiki/Gaussian_curvature for more info
template <typename Real>
inline Real gaussian_curvature(
    Real f_u, Real f_v, Real f_uu, Real f_vv, Real f_uv) noexcept {
    static_assert(std::is_floating_point_v<Real>);

    return (f_uu * f_vv - f_uv * f_uv) / (Real(1.) + f_u * f_u + f_v * f_v);
}

/// This function calculates the mean curvature for a function given its
/// derivatives.
/// \tparam Real precicions of the calculation
/// \param f_u,f_v,f_uu,f_vv,f_uv partial derivatives of the depth values
/// \returns the mean curvature for these derivatives
/// \sa conversion::depth_to_mean_curvature
/// \note check https://en.wikipedia.org/wiki/Mean_curvature for more info
template <typename Real>
inline Real
mean_curvature(Real f_u, Real f_v, Real f_uu, Real f_vv, Real f_uv) noexcept {
    static_assert(std::is_floating_point_v<Real>);

    using std::pow;
    using std::sqrt;
    // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers)
    return ((Real(1.) + f_v * f_v) * f_uu - Real(2.) * f_u * f_v * f_uv +
            (Real(1.) + f_u * f_u) * f_vv) /
           // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers)
           pow(Real(2.) * sqrt(Real(1.) + f_u * f_u + f_v * f_v), Real(3.));
}
}  // namespace sens_loc::math

#endif /* end of include guard: CURVATURE_H_LSXONAFA */

ext: h

max_stars_repo_name: JonasToth/depth-conversions

max_issues_repo_name: JonasToth/depth-conversions

max_forks_repo_name: JonasToth/depth-conversions


text:
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.h>
#include <float.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include "common.h"
#include "graph.h"
#include "utils.h"

void graph_init(graph *g, int vertices)
{
    g->num_v = vertices;
    g->num_e = (int*)malloc(vertices * sizeof(int));
    g->capacity = (int*)malloc(vertices * sizeof(int));
    g->adj_list = (int**)malloc(vertices * sizeof(int*));
    g->w_list = (double**)malloc(vertices * sizeof(double*));

    int i = 0;
    for (; i < vertices; ++i)
    {
        g->num_e[i] = 0;
        g->capacity[i] = GRAPH_INIT_CAPACITY;
        g->adj_list[i] = (int*)malloc(GRAPH_INIT_CAPACITY * sizeof(int));
        g->w_list[i] = (double*)malloc(GRAPH_INIT_CAPACITY * sizeof(double));
    }
}

ORIGIN_INLINE int graph_edges(const graph *g)
{
    const int num_v = g->num_v;
    int sum = 0;
    int v = 0;
    for (; v < num_v; ++v)
    {
        sum += g->num_e[v];
    }
    return sum;
}

ORIGIN_INLINE int graph_proper_edges(const graph *g)
{
    const int num_v = g->num_v;
    int sum = 0;
    int v = 0;
    for (; v < num_v; ++v)
    {
        int e = 0;
        for (; e < g->num_e[v]; ++e)
        {
            if (g->adj_list[v][e] != v)
            {
                ++sum;
            }
        }
    }
    return sum;
}

ORIGIN_INLINE int graph_loops(const graph *g)
{
    const int num_v = g->num_v;
    int sum = 0;
    int v = 0;
    for (; v < num_v; ++v)
    {
        int e = 0;
        for (; e < g->num_e[v]; ++e)
        {
            if (g->adj_list[v][e] == v)
            {
                ++sum;
            }
        }
    }
    return sum;
}

ORIGIN_INLINE int graph_outdegree(const graph *g, int u)
{
    return g->num_e[u];
}

ORIGIN_INLINE int graph_indegree(const graph *g, int u)
{
    const int num_v = g->num_v;
    int sum = 0;
    int v = 0;
    for (; v < num_v; ++v)
    {
        int e = 0;
        for (; e < g->num_e[v]; ++e)
        {
            if (g->adj_list[v][e] == u)
            {
                ++sum;
            }
        }
    }
    return sum;
}

ORIGIN_INLINE int graph_is_balanced(const graph *g)
{
    const int num_v = g->num_v;
    int u = 0;
    for (; u < num_v; ++u)
    {
        if (graph_indegree(g, u) != graph_outdegree(g, u))
        {
            return FALSE;
        }
    }
    return TRUE;
}

ORIGIN_INLINE void graph_add_edge(graph *g, int u, int v, double weight)
{
    if (g->num_e[u] == g->capacity[u])
    {
        graph_grow_lists(g, u);
    }
    g->adj_list[u][g->num_e[u]] = v;
    g->w_list[u][g->num_e[u]] = weight;
    g->num_e[u]++;
}

void graph_add_sym_edges(graph *g, int u, int v, double weight)
{
    graph_add_edge(g, u, v, weight);
    graph_add_edge(g, v, u, weight);
}

int graph_rmv_edge(graph *g, int u, int v)
{
    const int num_e = g->num_e[u];
    int e = 0;
    for (; e < num_e; ++e)
    {
        if (g->adj_list[u][e] == v)
        {
            g->num_e[u]--;
            int f = e;
            for (; f < num_e - 1; ++f)
            {
                g->adj_list[u][f] = g->adj_list[u][f+1];
                g->w_list[u][f] = g->w_list[u][f+1];
            }
            return TRUE;
        }
    }
    return FALSE;
}

int graph_rmv_sym_edges(graph *g, int u, int v)
{
    return (graph_rmv_edge(g, u, v) && graph_rmv_edge(g, v, u));
}

int graph_has_edge(graph *g, int u, int v)
{
    const int num_e = g->num_e[u];
    int e = 0;
    for (; e < num_e; ++e)
    {
        if (g->adj_list[u][e] == v)
        {
            return TRUE;
        }
    }
    return FALSE;
}

int graph_strongly_connected(const graph *g)
{
    const int num_v = g->num_v;
    int *group = (int*)malloc(num_v * sizeof(int));

    int u = 0;
    for (; u < num_v; ++u)
    {
        int v = 0;
        for (; v < num_v; ++v)
        {
            group[v] = FALSE;
        }

        group[u] = TRUE;
        graph_test_cc(g, group, u); // Call recursive function

        for (v = 0; v < num_v; ++v)
        {
            if (group[v] == FALSE)
            {
                free(group);
                return FALSE;
            }
        }
    }
    free(group);
    return TRUE;
}

void graph_svg(const graph *g, double *x, double *y, double size, double offset, FILE *out)
{
    fprintf(out, "<svg xmlns=\"http://www.w3.org/2000/svg\" version=\"1.1\">");
    const int num_v = g->num_v;
    for (int i = 0; i < num_v; ++i)
    {
        const int num_e = g->num_e[i];
        for (int j = 0; j < num_e; ++j)
        {
            const int e = g->adj_list[i][j];
            fprintf(out, "  <line x1=\"%f\" y1=\"%f\" x2=\"%f\" y2=\"%f\" style=\"stroke:rgb(0,0,0);stroke-width:2\"/>\n", offset + x[i] * size, offset + y[i] * size, offset + x[e] * size, offset + y[e] * size);
        }
    }
    for (int i = 0; i < num_v; ++i)
    {
        fprintf(out, "  <circle id=\"%d\" cx=\"%f\" cy=\"%f\" r=\"8\" fill=\"rgb(0,0,0)\"/>\n", i, offset + x[i] * size, offset + y[i] * size);
    }
    fprintf(out, "</svg>");
}

void graph_svg_abun(const graph *g, double *x, double *y, double size, double offset, double *abun, int color, FILE *out)
{
    fprintf(out, "<svg xmlns=\"http://www.w3.org/2000/svg\" version=\"1.1\">");
    const int num_v = g->num_v;
    for (int i = 0; i < num_v; ++i)
    {
        const int num_e = g->num_e[i];
        for (int j = 0; j < num_e; ++j)
        {
            const int e = g->adj_list[i][j];
            fprintf(out, "  <line x1=\"%f\" y1=\"%f\" x2=\"%f\" y2=\"%f\" style=\"stroke:rgb(0,0,0);stroke-width:2\"/>\n", offset + x[i] * size, offset + y[i] * size, offset + x[e] * size, offset + y[e] * size);
        }
    }
    int clr[3] = {0, 0, 0};
    for (int i = 0; i < num_v; ++i)
    {
        clr[color] = (int)(abun[i] * 255);
        fprintf(out, "  <circle id=\"%d\" cx=\"%f\" cy=\"%f\" r=\"8\" stroke=\"black\" stroke-width=\"1\" fill=\"rgb(%d,%d,%d)\"/>\n", i, offset + x[i] * size, offset + y[i] * size, clr[0], clr[1], clr[2]);
    }
    fprintf(out, "</svg>");
}

void graph_graphml(const graph *g, FILE *out, unsigned int id)
{
    const int num_v = g->num_v;
    fprintf(out, "<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n");
    fprintf(out, "<graphml xmlns=\"http://graphml.graphdrawing.org/xmlns\"");
    fprintf(out, " xmlns:xsi=\"http://www.w3.org/2001/XMLSchema-instance\"");
    fprintf(out, " xsi:schemaLocation=\"http://graphml.graphdrawing.org/xmlns http://graphml.graphdrawing.org/xmlns/1.0/graphml.xsd\">\n");

    fprintf(out, "  <graph id=\"%u\" edgedefault=\"directed\">\n", id);
    for (int i = 0; i < num_v; ++i)
    {
        fprintf(out, "    <node id=\"v%d\"/>\n", i);
    }
    for (int i = 0; i < num_v; ++i)
    {
        for (int e = 0; e < g->num_e[i]; ++e)
        {
            fprintf(out, "    <edge source=\"v%d\" target=\"v%d\"/>\n", i, g->adj_list[i][e]);
        }
    }
    fprintf(out, "  </graph>\n");
    fprintf(out, "</graphml>\n");
}

void graph_print(const graph *g, FILE *out)
{
    for (int i = 0; i < g->num_v; ++i)
    {
        fprintf(out, "%5d -> ", i);
        for (int e = 0; e < g->num_e[i]; ++e)
        {
            fprintf(out, "%d ", g->adj_list[i][e]);
        }
        fprintf(out, "\n");
    }
}

void graph_free(graph *g)
{
    const int num_v = g->num_v;
    for (int i = 0; i < num_v; ++i)
    {
        free(g->adj_list[i]);
        free(g->w_list[i]);
    }
    free(g->adj_list);
    g->adj_list = NULL;
    free(g->w_list);
    g->w_list = NULL;
    free(g->num_e);
    g->num_e = NULL;
    free(g->capacity);
    g->capacity = NULL;
}

void graph_get_rgg(graph *g, int vertices, double r, double *x, double *y, gsl_rng *rng)
{
    graph_init(g, vertices);

    for (int i = 0; i < vertices; ++i)
    {
        x[i] = gsl_rng_uniform(rng);
        y[i] = gsl_rng_uniform(rng);
    }
    double d;
    for (int i = 0; i < vertices; ++i)
    {
        for (int j = 0; j < vertices; ++j)
        {
            const double a = x[i] - x[j];
            const double b = y[i] - y[j];
            d = hypot(a, b);
            if (d < r)
            {
                graph_add_edge(g, i, j, r - d);
            }
        }
    }
}

void graph_get_crgg(graph *g, int vertices, double r, double *x, double *y, gsl_rng *rng)
{
    graph_init(g, vertices);
    graph_get_rgg(g, vertices, r, x, y, rng);

    while (graph_strongly_connected(g) == FALSE)
    {
        graph_free(g);
        graph_get_rgg(g, vertices, r, x, y, rng);
    }
}

void graph_get_rec_rgg(graph *g, int vertices, double width, double r, double *x, double *y, gsl_rng *rng)
{
    graph_init(g, vertices);

    const double length = 1.0 / width; // A = l * w so l = 1 / w

    for (int i = 0; i < vertices; ++i)
    {
        x[i] = gsl_rng_uniform(rng) * length;
        y[i] = gsl_rng_uniform(rng) * width;
    }
    double d;
    for (int i = 0; i < vertices; ++i)
    {
        for (int j = 0; j < vertices; ++j)
        {
            const double a = x[i] - x[j];
            const double b = y[i] - y[j];
            d = hypot(a, b);
            if (d < r)
            {
                graph_add_edge(g, i, j, r - d);
            }
        }
    }
}

void graph_get_rec_crgg(graph *g, int vertices, double width, double r, double *x, double *y, gsl_rng *rng)
{
    graph_init(g, vertices);
    graph_get_rec_rgg(g, vertices, width, r, x, y, rng);

    while (graph_strongly_connected(g) == FALSE)
    {
        graph_free(g);
        graph_get_rec_rgg(g, vertices, width, r, x, y, rng);
    }
}

void graph_get_complete(graph *g, int vertices)
{
    graph_init(g, vertices);

    for (int i = 0; i < vertices; ++i)
    {
        for (int j = 0; j < vertices; ++j)
        {
            graph_add_edge(g, i, j, 1.0);
        }
    }
}

void graph_get_circle(graph *g, int vertices)
{
    graph_init(g, vertices);

    graph_add_edge(g, vertices - 1, vertices - 1, 1.0);
    graph_add_sym_edges(g, 0, vertices - 1, 1.0);

    for (int u = 0; u < vertices - 1; ++u)
    {
        graph_add_edge(g, u, u, 1.0);
        graph_add_sym_edges(g, u, u + 1, 1.0);
    }
}

void graph_get_star(graph *g, int vertices)
{
    graph_init(g, vertices);

    for (int u = 0; u < vertices; ++u)
    {
        graph_add_edge(g, u, u, 1.0);
        graph_add_sym_edges(g, u, 0, 1.0);
    }
}

///////////////////////////////////////////////////////////////
// 'Private' functions

void graph_test_cc(const graph *g, int *group, int u)
{
    const int num_e = g->num_e[u];
    int e = 0;
    for (; e < num_e; ++e)
    {
        const int head = g->adj_list[u][e];
        if (head != u && group[head] == FALSE)
        {
            group[head] = TRUE;
            graph_test_cc(g, group, head);
        }
    }
}

// Grow the edge list of some vertex v.
void graph_grow_lists(graph *g, int u)
{
    g->capacity[u] <<= 1;

    int *tmp_adj = (int*)malloc(g->capacity[u] * sizeof(int));
    double *tmp_w = (double*)malloc(g->capacity[u] * sizeof(double));

    const int num_e = g->num_e[u];
    int e = 0;
    for (; e < num_e; ++e)
    {
        tmp_adj[e] = g->adj_list[u][e];
        tmp_w[e] = g->w_list[u][e];
    }
    int *swap_adj = g->adj_list[u];
    double *swap_w = g->w_list[u];

    g->adj_list[u] = tmp_adj;
    g->w_list[u] = tmp_w;

    free(swap_adj);
    free(swap_w);
}

max_stars_repo_name: PhDP/Origin

max_issues_repo_name: PhDP/Origin

max_forks_repo_name: PhDP/Origin


text:
/* File: set_cosmo_dist.c */
/*
  This file is a part of the Corrfunc package
  Copyright (C) 2015-- Manodeep Sinha (manodeep@gmail.com)
  License: MIT LICENSE. See LICENSE file under the top-level
  directory at https://github.com/manodeep/Corrfunc/
*/

/* function set_cosmo_dist (modelled after cosmodist.c)

   --- Computes the comoving distance as a function of redshift for
   a given cosmological model.
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_integration.h>

#include "set_cosmo_dist.h"
#include "cosmology_params.h"

#define SQR(x) ((x)*(x))
#define CUBE(x) ((x)*(x)*(x))
#define epsilon 1e-10

int set_cosmo_dist(const double zmax,const int max_size,double *zc,double *dc,const int lasdamas_cosmology)
{
    /* First, initialize cosmology*/
    int status = init_cosmology(lasdamas_cosmology);
    if(status != EXIT_SUCCESS) {
        return -1;
    }

    int i = 0;
    double Omegak;
    double smallh = 1.0;//Andreas pointed out that I don't need the real value of LITTLE_H
    Omegak = 1.0 - OMEGA_M - OMEGA_L;
    const double Dh = SPEED_OF_LIGHT*0.01/smallh ;// c/(100) -> in units of little h^-1 Mpc

    const double Deltaz = 1.0/max_size;
    const double dz = 1e-2*Deltaz;

    double Eint = 0.0;
    double E2 = 1.0 ;
    double z2 = Deltaz ;
    for(double z=2.0*dz;z<zmax;z+=2.0*dz) {
        double E0 = E2 ;

        double E1 = 1.0/sqrt(OMEGA_M * CUBE(1+z-dz) + Omegak *(1+z-dz) + OMEGA_L);
        E2 = 1.0/sqrt(OMEGA_M * CUBE(1+z) + Omegak *(1+z) + OMEGA_L);
        Eint += dz*(E0 + 4.*E1 + E2)/3. ;

        if(z>(z2-epsilon) && z<(z2+epsilon)) {
            if( i >= max_size )
                break ;

            const double Dc = Eint*Dh ;
            zc[i] = z ;
            dc[i] = Dc ;
            z2 += Deltaz ;
            i++ ;
        }
    }

#ifdef DEBUG
    fprintf(stderr,"%s> (Omega_m, Omega_L, Omega_k, h, zmax) = (%4.2f, %4.2f, %4.2f, %4.2f,%4.2lf)\n", __FUNCTION__, OMEGA_M,OMEGA_L,Omegak,LITTLE_H,zmax) ;
    fprintf(stderr,"%s> tabulated redshift: %g to %g  (distance: %g to %g Mpc)\n", __FUNCTION__, zc[0], zc[i-1], dc[0], dc[i-1]) ;
#endif
    return i ;
}

max_stars_repo_name: laperezNYC/Corrfunc

max_issues_repo_name: laperezNYC/Corrfunc

max_forks_repo_name: laperezNYC/Corrfunc


text:
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_rstat.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_sort.h>

int
main(void)
{
  const size_t N = 10000;
  double *data = malloc(N * sizeof(double));
  gsl_rstat_quantile_workspace *work_25 = gsl_rstat_quantile_alloc(0.25);
  gsl_rstat_quantile_workspace *work_50 = gsl_rstat_quantile_alloc(0.5);
  gsl_rstat_quantile_workspace *work_75 = gsl_rstat_quantile_alloc(0.75);
  gsl_rng *r = gsl_rng_alloc(gsl_rng_default);
  double exact_p25, exact_p50, exact_p75;
  double val_p25, val_p50, val_p75;
  size_t i;

  /* add data to quantile accumulators; also store data for exact
   * comparisons */
  for (i = 0; i < N; ++i)
    {
      data[i] = gsl_ran_rayleigh(r, 1.0);
      gsl_rstat_quantile_add(data[i], work_25);
      gsl_rstat_quantile_add(data[i], work_50);
      gsl_rstat_quantile_add(data[i], work_75);
    }

  /* exact values */
  gsl_sort(data, 1, N);
  exact_p25 = gsl_stats_quantile_from_sorted_data(data, 1, N, 0.25);
  exact_p50 = gsl_stats_quantile_from_sorted_data(data, 1, N, 0.5);
  exact_p75 = gsl_stats_quantile_from_sorted_data(data, 1, N, 0.75);

  /* estimated values */
  val_p25 = gsl_rstat_quantile_get(work_25);
  val_p50 = gsl_rstat_quantile_get(work_50);
  val_p75 = gsl_rstat_quantile_get(work_75);

  printf ("The dataset is %g, %g, %g, %g, %g, ...\n",
         data[0], data[1], data[2], data[3], data[4]);

  printf ("0.25 quartile: exact = %.5f, estimated = %.5f, error = %.6e\n",
          exact_p25, val_p25, (val_p25 - exact_p25) / exact_p25);
  printf ("0.50 quartile: exact = %.5f, estimated = %.5f, error = %.6e\n",
          exact_p50, val_p50, (val_p50 - exact_p50) / exact_p50);
  printf ("0.75 quartile: exact = %.5f, estimated = %.5f, error = %.6e\n",
          exact_p75, val_p75, (val_p75 - exact_p75) / exact_p75);

  gsl_rstat_quantile_free(work_25);
  gsl_rstat_quantile_free(work_50);
  gsl_rstat_quantile_free(work_75);
  gsl_rng_free(r);
  free(data);

  return 0;
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#pragma once

#include <gsl\gsl>
#include <winrt\Windows.Foundation.h>
#include <d3d11.h>
#include "DrawableGameComponent.h"
#include "MatrixHelper.h"
#include "PointLight.h"
#include "DepthMap.h"
#include "Rectangle.h"
#include "ProjectiveTextureMappingMaterial.h"
#include "RenderStateHelper.h"

namespace DirectX
{
	class SpriteBatch;
}

namespace Library
{
	class ProxyModel;
	class RenderableFrustum;
	class DepthMapMaterial;
}

namespace Rendering
{
	class ProjectiveTextureMappingDemo final : public Library::DrawableGameComponent
	{
	public:
		ProjectiveTextureMappingDemo(Library::Game& game, const std::shared_ptr<Library::Camera>& camera);
		ProjectiveTextureMappingDemo(const ProjectiveTextureMappingDemo&) = delete;
		ProjectiveTextureMappingDemo(ProjectiveTextureMappingDemo&&) = default;
		ProjectiveTextureMappingDemo& operator=(const ProjectiveTextureMappingDemo&) = default;
		ProjectiveTextureMappingDemo& operator=(ProjectiveTextureMappingDemo&&) = default;
		~ProjectiveTextureMappingDemo();

		ProjectiveTextureMappingDrawModes DrawMode() const;
		const std::string& DrawModeString() const;
		void SetDrawMode(ProjectiveTextureMappingDrawModes drawMode);

		float AmbientLightIntensity() const;
		void SetAmbientLightIntensity(float intensity);

		float PointLightIntensity() const;
		void SetPointLightIntensity(float intensity);

		float PointLightRadius() const;
		void SetPointLightRadius(float radius);

		const DirectX::XMFLOAT3& ProjectorPosition() const;
		const DirectX::XMVECTOR ProjectorPositionVector() const;
		void SetProjectorPosition(const DirectX::XMFLOAT3& position);
		void SetProjectorPosition(DirectX::FXMVECTOR position);

		const DirectX::XMFLOAT3& ProjectorDirection() const;
		void RotateProjector(const DirectX::XMFLOAT2& amount);

		virtual void Initialize() override;
		virtual void Update(const Library::GameTime& gameTime) override;
		virtual void Draw(const Library::GameTime& gameTime) override;

	private:
		inline static const std::uint32_t DepthMapWidth{ 1024 };
		inline static const std::uint32_t DepthMapHeight{ 1024 };
		static inline const RECT DepthMapDestinationRectangle{ 0, 512, 256, 768 };

		void DrawWithDepthMap();
		void DrawWithoutDepthMap();
		void UpdateTransforms(ProjectiveTextureMappingMaterial::VertexCBufferPerObject& transforms, DirectX::FXMMATRIX worldViewProjectionMatrix, DirectX::CXMMATRIX worldMatrix, DirectX::CXMMATRIX projectiveTextureMatrix);
		void InitializeProjectedTextureScalingMatrix(uint32_t textureWidth, uint32_t textureHeight);

		ProjectiveTextureMappingMaterial::VertexCBufferPerObject mPlaneTransforms;
		ProjectiveTextureMappingMaterial::VertexCBufferPerObject mTeapotTransforms;
		DirectX::XMFLOAT4X4 mPlaneWorldMatrix{ Library::MatrixHelper::Identity };
		DirectX::XMFLOAT4X4 mTeapotWorldMatrix{ Library::MatrixHelper::Identity };
		DirectX::XMFLOAT4X4 mProjectedTextureScalingMatrix{ Library::MatrixHelper::Zero };
		Library::PointLight mPointLight;
		Library::DepthMap mDepthMap;
		Library::RenderStateHelper mRenderStateHelper;
		std::shared_ptr<ProjectiveTextureMappingMaterial> mMaterial;
		std::shared_ptr<Library::DepthMapMaterial> mDepthMapMaterial;
		winrt::com_ptr<ID3D11Buffer> mPlaneVertexBuffer;
		winrt::com_ptr<ID3D11Buffer> mTeapotVertexBuffer;
		winrt::com_ptr<ID3D11Buffer> mTeapotPositionOnlyVertexBuffer;
		winrt::com_ptr<ID3D11Buffer> mTeapotIndexBuffer;
		std::uint32_t mPlaneVertexCount{ 0 };
		std::uint32_t mTeapotIndexCount{ 0 };
		std::unique_ptr<Library::ProxyModel> mProxyModel;
		std::unique_ptr<Library::Camera> mProjector;
		std::unique_ptr<Library::RenderableFrustum> mRenderableProjectorFrustum;
		std::unique_ptr<DirectX::SpriteBatch> mSpriteBatch;
		bool mUpdateMaterial{ true };
	};
}


text:
/* movstat/test_Sn.c
 *
 * Copyright (C) 2018 Patrick Alken
 *
 * 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.
 */

#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_sort_vector.h>
#include <gsl/gsl_test.h>

/* calculate S_n statistic for input vector using slow/naive algorithm */
static int
slow_movSn(const gsl_movstat_end_t etype, const gsl_vector * x, gsl_vector * y,
           const int H, const int J)
{
  const size_t n = x->size;
  const int K = H + J + 1;
  double *window = malloc(K * sizeof(double));
  double *work = malloc(K * sizeof(double));
  size_t i;

  for (i = 0; i < n; ++i)
    {
      size_t wsize = gsl_movstat_fill(etype, x, i, H, J, window);
      double Sn;

      gsl_sort(window, 1, wsize);
      Sn = gsl_stats_Sn_from_sorted_data(window, 1, wsize, work);
      gsl_vector_set(y, i, Sn);
    }

  free(window);
  free(work);

  return GSL_SUCCESS;
}

static double
func_Sn(const size_t n, double x[], void * params)
{
  double *work = malloc(n * sizeof(double));
  double Sn;

  (void) params;

  gsl_sort(x, 1, n);
  Sn = gsl_stats_Sn_from_sorted_data(x, 1, n, work);

  free(work);

  return Sn;
}

static void
test_Sn_proc(const double tol, const size_t n, const size_t H, const size_t J,
             const gsl_movstat_end_t etype, gsl_rng *rng_p)
{
  gsl_movstat_workspace *w;
  gsl_vector *x = gsl_vector_alloc(n);
  gsl_vector *y = gsl_vector_alloc(n);
  gsl_vector *z = gsl_vector_alloc(n);
  gsl_movstat_function F;
  char buf[2048];

  F.function = func_Sn;
  F.params = NULL;

  if (H == J)
    w = gsl_movstat_alloc(2*H + 1);
  else
    w = gsl_movstat_alloc2(H, J);

  /* test moving median with random input */
  random_vector(x, rng_p);

  /* y = S_n(x) with slow brute force algorithm */
  slow_movSn(etype, x, y, H, J);

  /* z = S_n(x) */
  gsl_movstat_Sn(etype, x, z, w);

  /* test y = z */
  sprintf(buf, "n=%zu H=%zu J=%zu endtype=%u Sn random", n, H, J, etype);
  compare_vectors(tol, z, y, buf);

  /* z = S_n(x) in-place */
  gsl_vector_memcpy(z, x);
  gsl_movstat_Sn(etype, z, z, w);

  sprintf(buf, "n=%zu H=%zu J=%zu endtype=%u Sn random in-place", n, H, J, etype);
  compare_vectors(tol, z, y, buf);

  /* z = S_n(x) with user-defined function */
  gsl_movstat_apply(etype, &F, x, z, w);

  sprintf(buf, "n=%zu H=%zu J=%zu endtype=%u Sn user", n, H, J, etype);
  compare_vectors(tol, z, y, buf);

  gsl_vector_free(x);
  gsl_vector_free(y);
  gsl_vector_free(z);
  gsl_movstat_free(w);
}

static void
test_Sn(gsl_rng * rng_p)
{
  test_Sn_proc(GSL_DBL_EPSILON, 1000, 0, 0, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_Sn_proc(GSL_DBL_EPSILON, 1000, 5, 5, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_Sn_proc(GSL_DBL_EPSILON, 1000, 5, 2, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_Sn_proc(GSL_DBL_EPSILON, 1000, 2, 5, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_Sn_proc(GSL_DBL_EPSILON, 2000, 50, 0, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_Sn_proc(GSL_DBL_EPSILON, 2000, 0, 50, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_Sn_proc(GSL_DBL_EPSILON, 20, 50, 50, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_Sn_proc(GSL_DBL_EPSILON, 20, 1, 50, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_Sn_proc(GSL_DBL_EPSILON, 20, 50, 1, GSL_MOVSTAT_END_PADZERO, rng_p);

  test_Sn_proc(GSL_DBL_EPSILON, 1000, 0, 0, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_Sn_proc(GSL_DBL_EPSILON, 1000, 5, 5, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_Sn_proc(GSL_DBL_EPSILON, 1000, 5, 2, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_Sn_proc(GSL_DBL_EPSILON, 1000, 2, 5, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_Sn_proc(GSL_DBL_EPSILON, 2000, 50, 0, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_Sn_proc(GSL_DBL_EPSILON, 2000, 0, 50, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_Sn_proc(GSL_DBL_EPSILON, 20, 50, 50, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_Sn_proc(GSL_DBL_EPSILON, 20, 1, 50, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_Sn_proc(GSL_DBL_EPSILON, 20, 50, 1, GSL_MOVSTAT_END_PADVALUE, rng_p);

  test_Sn_proc(GSL_DBL_EPSILON, 1000, 0, 0, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_Sn_proc(GSL_DBL_EPSILON, 1000, 5, 5, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_Sn_proc(GSL_DBL_EPSILON, 1000, 5, 2, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_Sn_proc(GSL_DBL_EPSILON, 1000, 2, 5, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_Sn_proc(GSL_DBL_EPSILON, 2000, 50, 0, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_Sn_proc(GSL_DBL_EPSILON, 2000, 0, 50, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_Sn_proc(GSL_DBL_EPSILON, 20, 50, 50, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_Sn_proc(GSL_DBL_EPSILON, 20, 1, 50, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_Sn_proc(GSL_DBL_EPSILON, 20, 50, 1, GSL_MOVSTAT_END_TRUNCATE, rng_p);
}

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
/* Author: G. Jungman
 */
#include <config.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_errno.h>
#include "gsl_qrng.h"

gsl_qrng *
gsl_qrng_alloc (const gsl_qrng_type * T, unsigned int dimension)
{

  gsl_qrng * r = (gsl_qrng *) malloc (sizeof (gsl_qrng));

  if (r == 0)
    {
      GSL_ERROR_VAL ("allocation failed for qrng struct",
			GSL_ENOMEM, 0);
    };

  r->dimension = dimension;
  r->state_size = T->state_size(dimension);
  r->state = malloc (r->state_size);

  if (r->state == 0)
    {
      free (r);
      GSL_ERROR_VAL ("allocation failed for qrng state",
			GSL_ENOMEM, 0);
    };

  r->type = T;

  T->init_state(r->state, r->dimension);

  return r;
}

int
gsl_qrng_memcpy (gsl_qrng * dest, const gsl_qrng * src)
{
  if (dest->type != src->type)
    {
      GSL_ERROR ("generators must be of the same type", GSL_EINVAL);
    }

  dest->dimension = src->dimension;
  dest->state_size = src->state_size;
  memcpy (dest->state, src->state, src->state_size);

  return GSL_SUCCESS;
}

gsl_qrng *
gsl_qrng_clone (const gsl_qrng * q)
{
  gsl_qrng * r = (gsl_qrng *) malloc (sizeof (gsl_qrng));

  if (r == 0)
    {
      GSL_ERROR_VAL ("failed to allocate space for rng struct",
			GSL_ENOMEM, 0);
    };

  r->dimension = q->dimension;
  r->state_size = q->state_size;
  r->state = malloc (r->state_size);

  if (r->state == 0)
    {
      free (r);
      GSL_ERROR_VAL ("failed to allocate space for rng state",
			GSL_ENOMEM, 0);
    };

  r->type = q->type;

  memcpy (r->state, q->state, q->state_size);

  return r;
}

#ifndef HIDE_INLINE_STATIC
int
gsl_qrng_get (const gsl_qrng * r, double x[])
{
  return (r->type->get) (r->state, r->dimension, x);
}
#endif

const char *
gsl_qrng_name (const gsl_qrng * r)
{
  return r->type->name;
}

size_t
gsl_qrng_size (const gsl_qrng * r)
{
  return r->state_size;
}

void *
gsl_qrng_state (const gsl_qrng * r)
{
  return r->state;
}

void
gsl_qrng_free (gsl_qrng * r)
{
  if(r != 0) {
    if(r->state != 0) free (r->state);
    free (r);
  }
}