text:
#ifndef DBJ_FLOAT_VARIOUS_MATMULS_INC
#define DBJ_FLOAT_VARIOUS_MATMULS_INC

#include "dbj_matrix.h"

#include <stdlib.h>
#include <intrin.h> // #define __SSE__ 1
#ifdef HAVE_CBLAS
#include <cblas.h>
#endif

#undef DBJ_API
#define DBJ_API

//////////////////////////////////////////////////////////////////////
// all of the bellow uses:
typedef DBJ_MATRIX_STRUCT(float) float_matrix_struct;

//////////////////////////////////////////////////////////////////////

typedef void* (*mat_mul_function) (
	const unsigned, const unsigned, const unsigned,
	float a[][*], float b[][*], float m[][*]
	);

typedef struct {
	const char* description;
	mat_mul_function function;
} description_function_pair;

#define DF_PAIR(D_,F_) (description_function_pair){ .description = D_ , .function = F_ }

//////////////////////////////////////////////////////////////////////

DBJ_API float_matrix_struct* new_float_matrix(const unsigned /*n_rows*/, const unsigned  /*n_cols*/);
DBJ_API float_matrix_struct* make_random_float_matrix(
	const unsigned /*n_rows*/, const unsigned  /*n_cols*/, float (*)(void)
);
DBJ_API void* mat_transpose(
	const unsigned /*n_rows*/, const unsigned /*n_cols*/,
	float /*a*/[/*static n_rows*/][*/*n_cols*/], float /*m*/[/*static n_cols*/][*/*n_rows*/]
);

DBJ_API void* mat_mul_null(
	const unsigned n_a_rows, const unsigned n_a_cols, const unsigned n_b_cols,
	float a[static n_a_rows][n_a_cols], float b[static n_a_rows][n_b_cols], float m[static n_a_rows][n_b_cols]
) {
	// used for when algorithms are not implemented yet
	(void)a;	(void)b;	(void)m;	return NULL;
}

DBJ_API void* mat_mul_0(const unsigned /*n_a_rows*/, const unsigned /*n_a_cols*/, const unsigned /*n_b_cols*/,
	float a[ /*n_a_rows*/][* /*n_a_cols*/], float b[ /*n_a_rows*/][* /*n_b_cols*/], float m[ /*n_a_rows*/][* /*n_b_cols*/]
);

DBJ_API void* mat_mul_1(const unsigned /*n_a_rows*/, const unsigned /*n_a_cols*/, const unsigned /*n_b_cols*/,
	float a[ /*n_a_rows*/][* /*n_a_cols*/], float b[ /*n_a_rows*/][* /*n_b_cols*/], float m[ /*n_a_rows*/][* /*n_b_cols*/]
);

DBJ_API void* mat_mul_3(const unsigned /*n_a_rows*/, const unsigned /*n_a_cols*/, const unsigned /*n_b_cols*/,
	float a[ /*n_a_rows*/][* /*n_a_cols*/], float b[ /*n_a_rows*/][* /*n_b_cols*/], float m[ /*n_a_rows*/][* /*n_b_cols*/]
);

DBJ_API void* mat_mul_4(const unsigned /*n_a_rows*/, const unsigned /*n_a_cols*/, const unsigned /*n_b_cols*/,
	float a[ /*n_a_rows*/][* /*n_a_cols*/], float b[ /*n_a_rows*/][* /*n_b_cols*/], float m[ /*n_a_rows*/][* /*n_b_cols*/]
);

#if __SSE__

DBJ_API void* mat_mul_2(const unsigned /*n_a_rows*/, const unsigned /*n_a_cols*/, const unsigned /*n_b_cols*/,
	float a[ /*n_a_rows*/][* /*n_a_cols*/], float b[ /*n_a_rows*/][* /*n_b_cols*/], float m[ /*n_a_rows*/][* /*n_b_cols*/]
);

DBJ_API void* mat_mul_7(const unsigned /*n_a_rows*/, const unsigned /*n_a_cols*/, const unsigned /*n_b_cols*/,
	float a[ /*n_a_rows*/][* /*n_a_cols*/], float b[ /*n_a_rows*/][* /*n_b_cols*/], float m[ /*n_a_rows*/][* /*n_b_cols*/]
);

#endif // __SSE__

#ifdef HAVE_CBLAS
DBJ_API void* mat_mul_5(const unsigned /*n_a_rows*/, const unsigned /*n_a_cols*/, const unsigned /*n_b_cols*/,
	float a[ /*n_a_rows*/][* /*n_a_cols*/], float b[ /*n_a_rows*/][* /*n_b_cols*/], float m[ /*n_a_rows*/][* /*n_b_cols*/]
);

DBJ_API void* mat_mul_6(const unsigned /*n_a_rows*/, const unsigned /*n_a_cols*/, const unsigned /*n_b_cols*/,
	float a[ /*n_a_rows*/][* /*n_a_cols*/], float b[ /*n_a_rows*/][* /*n_b_cols*/], float m[ /*n_a_rows*/][* /*n_b_cols*/]
);
#endif // HAVE_CBLAS

const description_function_pair dbj_float_matmuls_algo_table[] =
{
DF_PAIR("0: naive - no optimization" , mat_mul_0) ,
DF_PAIR("1: transposing the second matrix" , mat_mul_1) ,
DF_PAIR("2: explicitly vectorized sdot() with SSE" , mat_mul_2) , /* requires __SSE__ */
DF_PAIR("3: implicitly vectorized sdot()" , mat_mul_3) ,
DF_PAIR("4: no vectorization hints", mat_mul_4) ,
#ifdef HAVE_CBLAS
DF_PAIR("5: with sdot() from an external CBLAS library" , mat_mul_5) ,
DF_PAIR("6: with sgemm() from an external CBLAS library", mat_mul_6) ,
#else
DF_PAIR("5: with sdot() from CBLAS library not implemented " , mat_mul_null) ,
DF_PAIR("6: with sgemm() from CBLAS library not implemented ", mat_mul_null) ,
#endif // ! HAVE_CBLAS
DF_PAIR("7: explicitly SSE sdot() plus loop tiling", mat_mul_7) /* requires __SSE__ */
};

static const unsigned dbj_float_matmuls_algo_table_size = (sizeof(dbj_float_matmuls_algo_table) / sizeof(dbj_float_matmuls_algo_table[0]));

//////////////////////////////////////////////////////////////////////
#ifdef DBJ_FLOAT_VARIOUS_MATMULS_IMPLEMENTATION
//////////////////////////////////////////////////////////////////////

DBJ_API float_matrix_struct* new_float_matrix(const unsigned n_rows, const unsigned  n_cols)
{
	// float_matrix_struct* retval = NULL; // must be init to NULL

	float_matrix_struct* retval = DBJ_MATRIX_ALLOC(n_rows, n_cols, DBJ_MATRIX_STRUCT_SIZE(float, n_rows, n_cols));

	if (retval) {
		retval->rows = n_rows;
		retval->cols = n_cols;
	}
	else {
		fprintf(stderr, "\n%s(%d) DBJ_MATRIX_NEW failed?\n", __FILE__, __LINE__);
		perror(" ");
		exit(1);
	}
	return retval;
}

/*
allocate new matrix and populate it with random float's
*/
DBJ_API float_matrix_struct* make_random_float_matrix(
	const unsigned n_rows, const unsigned  n_cols, float (*float_rand)(void)
)
{
	float_matrix_struct* fmt_pointer = new_float_matrix(n_rows, n_cols);

	DBJ_MATRIX_ALIAS(matrix, float, n_cols);
	DBJ_MATRIX_CAST(mx, matrix, fmt_pointer);

	for (unsigned i = 0; i < n_rows; ++i)
		for (unsigned j = 0; j < n_cols; ++j)
			mx[i][j] = float_rand(); //
	return fmt_pointer;
}

DBJ_API void* mat_transpose(
	const unsigned n_rows, const unsigned n_cols,
	float a[static n_rows][n_cols], float m[static n_cols][n_rows]
)
{
	for (unsigned i = 0; i < n_rows; ++i)
		for (unsigned j = 0; j < n_cols; ++j)
			m[j][i] = a[i][j];
	return m;
}

DBJ_API float sdot_1(int n, const float x[static n], const float y[static n])
{
	float s = 0.0f;
	for (int i = 0; i < n; ++i) s += x[i] * y[i];
	return s;
}

DBJ_API float sdot_8(int n, const float x[static n], const float y[static n])
{
	int i, n8 = n >> 3 << 3;
	float s = 0.0f, t[8] = { 0.0f };
	// t[0] = t[1] = t[2] = t[3] = t[4] = t[5] = t[6] = t[7] = 0.0f;
	for (i = 0; i < n8; i += 8) {
		t[0] += x[i + 0] * y[i + 0];
		t[1] += x[i + 1] * y[i + 1];
		t[2] += x[i + 2] * y[i + 2];
		t[3] += x[i + 3] * y[i + 3];
		t[4] += x[i + 4] * y[i + 4];
		t[5] += x[i + 5] * y[i + 5];
		t[6] += x[i + 6] * y[i + 6];
		t[7] += x[i + 7] * y[i + 7];
	}
	for (s = 0.0f; i < n; ++i) s += x[i] * y[i];
	s += t[0] + t[1] + t[2] + t[3] + t[4] + t[5] + t[6] + t[7];
	return s;
}

#ifdef __SSE__

DBJ_API float sdot_sse(int n, const float x[static n], const float y[static n])
{
	int i, n8 = n >> 3 << 3;
	__m128 vs1, vs2;
	float s, t[4];
	vs1 = _mm_setzero_ps();
	vs2 = _mm_setzero_ps();
	for (i = 0; i < n8; i += 8) {
		__m128 vx1, vx2, vy1, vy2;
		vx1 = _mm_loadu_ps(&x[i]);
		vx2 = _mm_loadu_ps(&x[i + 4]);
		vy1 = _mm_loadu_ps(&y[i]);
		vy2 = _mm_loadu_ps(&y[i + 4]);
		vs1 = _mm_add_ps(vs1, _mm_mul_ps(vx1, vy1));
		vs2 = _mm_add_ps(vs2, _mm_mul_ps(vx2, vy2));
	}
	for (s = 0.0f; i < n; ++i) s += x[i] * y[i];
	_mm_storeu_ps(t, vs1);
	s += t[0] + t[1] + t[2] + t[3];
	_mm_storeu_ps(t, vs2);
	s += t[0] + t[1] + t[2] + t[3];
	return s;
}
#endif // __SSE__

/**************************************************

 Various Matrix multiplication algorithms

 NOTE: check and choose for your RT Env
	   there are differences but are *very* dependant on
	   compiler, OS and hardware

 */

DBJ_API void* mat_mul_0(
	const unsigned n_a_rows, const unsigned n_a_cols, const unsigned n_b_cols,
	float a[static n_a_rows][n_a_cols], float b[static n_a_rows][n_b_cols], float m[static n_a_rows][n_b_cols]
)
{
	for (unsigned i = 0; i < n_a_rows; ++i) {
		for (unsigned j = 0; j < n_b_cols; ++j) {
			float t = 0.0;
			for (unsigned k = 0; k < n_a_cols; ++k)
				t += a[i][k] * b[k][j];
			m[i][j] = t;
		}
	}
	return m;
}

DBJ_API void* mat_mul_1(
	const unsigned n_a_rows, const unsigned n_a_cols, const unsigned n_b_cols,
	float a[static n_a_rows][n_a_cols], float b[static n_a_rows][n_b_cols], float m[static n_a_rows][n_b_cols]
)
{
	const unsigned n_b_rows = n_a_cols;
	// Temp rows and cols are inverted
	float_matrix_struct* Temp = new_float_matrix(n_b_cols, n_b_rows);
	DBJ_MATRIX_ALIAS(matrix, float, n_b_rows);
	DBJ_MATRIX_CAST(bT, matrix, Temp);
	(void)mat_transpose(n_b_rows, n_b_cols, b, bT);

	for (unsigned i = 0; i < n_a_rows; ++i) {
		const float* ai = a[i];
		float* mi = m[i];
		for (unsigned j = 0; j < n_b_cols; ++j) {
			float t = 0.0f, * bTj = bT[j];
			for (unsigned k = 0; k < n_a_cols; ++k)
				t += ai[k] * bTj[k];
			mi[j] = t;
		}
	}

	free(Temp);

	return m;
}

DBJ_API void* mat_mul_2(
	const unsigned n_a_rows, const unsigned n_a_cols, const unsigned n_b_cols,
	float a[static n_a_rows][n_a_cols], float b[static n_a_rows][n_b_cols], float m[static n_a_rows][n_b_cols]
)
{
#ifndef __SSE__
#error __SEE__ is required here
#endif

	const unsigned n_b_rows = n_a_cols;

	float_matrix_struct* Temp = new_float_matrix(n_b_cols, n_b_rows);
	DBJ_MATRIX_ALIAS(matrix, float, n_b_rows);
	DBJ_MATRIX_CAST(bT, matrix, Temp);
	(void)mat_transpose(n_b_rows, n_b_cols, b, bT);

	for (unsigned i = 0; i < n_a_rows; ++i)
		for (unsigned j = 0; j < n_b_cols; ++j)
			m[i][j] = sdot_sse(n_a_cols, a[i], bT[j]);
	free(Temp);
	return m;
}

DBJ_API void* mat_mul_7(
	const unsigned n_a_rows, const unsigned n_a_cols, const unsigned n_b_cols,
	float a[static n_a_rows][n_a_cols], float b[static n_a_rows][n_b_cols], float m[static n_a_rows][n_b_cols]
)
{
#ifndef __SSE__
#error __SEE__ is required here
#endif

	const unsigned x = 16, n_b_rows = n_a_cols;

	float_matrix_struct* Temp = new_float_matrix(n_b_cols, n_b_rows);
	DBJ_MATRIX_ALIAS(matrix, float, n_b_rows);
	DBJ_MATRIX_CAST(bT, matrix, Temp);
	(void)mat_transpose(n_b_rows, n_b_cols, b, bT);

	for (unsigned i = 0; i < n_a_rows; i += x) {
		for (unsigned j = 0; j < n_b_cols; j += x) {
			unsigned je = n_b_cols < j + x ? n_b_cols : j + x;
			unsigned ie = n_a_rows < i + x ? n_a_rows : i + x;
			for (unsigned ii = i; ii < ie; ++ii)
				for (unsigned jj = j; jj < je; ++jj)
					m[ii][jj] += sdot_sse(n_a_cols, a[ii], bT[jj]);
		}
	}
	free(Temp);
	return m;
}

/////////////////////////////////////////////////////////////////////////////////////

DBJ_API void* mat_mul_3(
	const unsigned n_a_rows, const unsigned n_a_cols, const unsigned n_b_cols,
	float a[static n_a_rows][n_a_cols], float b[static n_a_rows][n_b_cols], float m[static n_a_rows][n_b_cols]
)
{
	int  n_b_rows = n_a_cols;

	float_matrix_struct* Temp = new_float_matrix(n_b_cols, n_b_rows);
	DBJ_MATRIX_ALIAS(matrix, float, n_b_rows);
	DBJ_MATRIX_CAST(bT, matrix, Temp);
	(void)mat_transpose(n_b_rows, n_b_cols, b, bT);

	for (unsigned i = 0; i < n_a_rows; ++i)
		for (unsigned j = 0; j < n_b_cols; ++j)
			m[i][j] = sdot_8(n_a_cols, a[i], bT[j]);
	free(Temp);
	return m;
}

DBJ_API void* mat_mul_4(
	const unsigned n_a_rows, const unsigned n_a_cols, const unsigned n_b_cols,
	float a[static n_a_rows][n_a_cols], float b[static n_a_rows][n_b_cols], float m[static n_a_rows][n_b_cols]
)
{
	int n_b_rows = n_a_cols;

	float_matrix_struct* Temp = new_float_matrix(n_b_cols, n_b_rows);
	DBJ_MATRIX_ALIAS(matrix, float, n_b_rows);
	DBJ_MATRIX_CAST(bT, matrix, Temp);
	(void)mat_transpose(n_b_rows, n_b_cols, b, bT);

	for (unsigned i = 0; i < n_a_rows; ++i)
		for (unsigned j = 0; j < n_b_cols; ++j)
			m[i][j] = sdot_1(n_a_cols, a[i], bT[j]);
	free(Temp);
	return m;
}

#ifdef HAVE_CBLAS

DBJ_API void* mat_mul_5(
	const unsigned n_a_rows, const unsigned n_a_cols, const unsigned n_b_cols,
	float a[static n_a_rows][n_a_cols], float b[static n_a_rows][n_b_cols], float m[static n_a_rows][n_b_cols]
)
{
	const unsigned n_b_rows = n_a_cols;

	float_matrix_struct* Temp = new_float_matrix(n_b_cols, n_b_rows);
	DBJ_MATRIX_ALIAS(matrix, float, n_b_rows);
	DBJ_MATRIX_CAST(bT, matrix, Temp);
	(void)mat_transpose(n_b_rows, n_b_cols, b, bT);

	for (unsigned i = 0; i < n_a_rows; ++i)
		for (unsigned j = 0; j < n_b_cols; ++j)
			m[i][j] = cblas_sdot(n_a_cols, a[i], 1, bT[j], 1); // clas_sdot() ???
	free(Temp);
	return m;
}

DBJ_API void* mat_mul_6(
	const unsigned n_a_rows, const unsigned n_a_cols, const unsigned n_b_cols,
	float a[static n_a_rows][n_a_cols], float b[static n_a_rows][n_b_cols], float m[static n_a_rows][n_b_cols]
)
{
	cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, n_a_rows, n_b_cols, n_a_cols, 1.0f, a[0], n_a_rows, b[0], n_b_rows, 0.0f, m[0], n_a_rows);
	return m;
}
#endif // HAVE_CBLAS

//////////////////////////////////////////////////////////////////////
#endif // DBJ_FLOAT_VARIOUS_MATMULS_IMPLEMENTATION
//////////////////////////////////////////////////////////////////////

#undef DF_PAIR

#endif //  DBJ_FLOAT_VARIOUS_MATMULS_INC


text:
/*! \file ExtremaRootFinding.h
 *  \brief header file for finding roots and extrema of a function
    \todo not implemented yet
 */

#ifndef MINMAXROOT_H
#define MINMAXROOT_H

#include <cmath>
#include <iostream>
#include <gsl/gsl_math.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>
#include <Precision.h>
#include <Function.h>
#include <Matrix.h>
#include <GMatrix.h>
#include <Interpolate.h>

using namespace std;
namespace Math
{
    //// Brent bracket search for extrema
    //Double_t Brent();
}

#endif

ext: h

max_stars_repo_name: ICRAR/NBodylib

max_issues_repo_name: ICRAR/NBodylib

max_forks_repo_name: ICRAR/NBodylib


text:
/*
** ECM - eigenvector centrality mapping, I/O interface
**
** G.Lohmann, MPI-KYB, updated Oct 2018
*/

#include <viaio/Vlib.h>
#include <viaio/VImage.h>
#include <viaio/mu.h>
#include <viaio/option.h>

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_math.h>

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

#ifdef _OPENMP
#include <omp.h>
#endif /*_OPENMP*/

#define SQR(x) ((x) * (x))
#define ReLU(x) ((x) > 0 ? (x) : 0)

VImage WriteOutput(VAttrList list,VImage map,float *eigvec,size_t n)
{
  size_t i;
  int b,r,c;
  int nrows=0,ncols=0,ntimesteps=0;
  int nslices = VAttrListNumImages(list);
  VImage *src = VAttrListGetImages(list,nslices);
  VImageDimensions(src,nslices,&ntimesteps,&nrows,&ncols);

  VImage dest = VCreateImage(nslices,nrows,ncols,VFloatRepn);
  VFillImage(dest,VAllBands,0);
  VCopyImageAttrs (src[0], dest);
  VSetAttr(VImageAttrList(dest),"modality",NULL,VStringRepn,"ECM");

  for (i=0; i<n; i++) {
    b = VPixel(map,0,0,i,VInteger);
    r = VPixel(map,0,1,i,VInteger);
    c = VPixel(map,0,2,i,VInteger);
    if (b < 0 || r < 0 || c < 0) VError("  WriteOutput: illegal address %d %d %d",b,r,c);
    VPixel(dest,b,r,c,VFloat) = eigvec[i];
  }
  return dest;
}

void NormalizeData(gsl_matrix_float *X,int nt)
{
  int i,j;
  double s1,s2,mean,sd,u,nx,tiny=1.0e-6;
  float *arr=NULL;
  nx = (double)nt;

  for (i=0; i<X->size1; i++) {
    s1=s2=0;
    arr = gsl_matrix_float_ptr(X,i,0);
    for (j=0; j<nt; j++) {
      u = (double)arr[j];
      s1 += u;
      s2 += u*u;
    }
    mean = s1/nx;
    sd = sqrt((s2 - nx * mean * mean) / (nx - 1.0));

    arr = gsl_matrix_float_ptr(X,i,0);
    for (j=0; j<nt; j++) {
      u = (double)arr[j];
      if (sd > tiny) arr[j] = (float)((u-mean)/sd);
      else arr[j] = 0;
    }
  }

  for (i=0; i<X->size1; i++) {
    for (j=0; j<X->size2; j++) {
      double u = (double)gsl_matrix_float_get(X,i,j);
      if (gsl_isnan(u) || gsl_isinf(u)) gsl_matrix_float_set(X,i,j,0);
    }
  }
}

VImage VoxelMap(VAttrList list,VImage mask,size_t nvox)
{
  int b,r,c;

  /* read image dims */
  int nrows=0,ncols=0,ntimesteps=0;
  int nslices = VAttrListNumImages(list);
  VImage *src = VAttrListGetImages(list,nslices);
  VImageDimensions(src,nslices,&ntimesteps,&nrows,&ncols);

  /*  voxel addresses */
  VImage map = VCreateImage(1,5,nvox,VIntegerRepn);
  if (map == NULL) VError(" error allocating addr map");
  VFillImage(map,VAllBands,0);
  VPixel(map,0,3,0,VInteger) = nslices;
  VPixel(map,0,3,1,VInteger) = nrows;
  VPixel(map,0,3,2,VInteger) = ncols;

  int i = 0;
  for (b=0; b<nslices; b++) {
    for (r=0; r<nrows; r++) {
      for (c=0; c<ncols; c++) {
	if (VGetPixel(mask,b,r,c) < 0.1) continue;
	VPixel(map,0,0,i,VInteger) = b;
	VPixel(map,0,1,i,VInteger) = r;
	VPixel(map,0,2,i,VInteger) = c;
	i++;
      }
    }
  }
  return map;
}

/* count number of voxels */
size_t NumVoxels(VImage mask)
{
  int b,r,c;
  size_t nvox = 0;
  for (b=0; b<VImageNBands(mask); b++) {
    for (r=0; r<VImageNRows(mask); r++) {
      for (c=0; c<VImageNColumns(mask); c++) {
	if (VGetPixel(mask,b,r,c) < 0.1) continue;
	nvox++;
      }
    }
  }
  fprintf(stderr," nvoxels: %ld\n",(long)nvox);
  return nvox;
}

gsl_matrix_float *ReadDataECM(VAttrList list,VImage mask,VImage map,VShort first,VShort length,int type,size_t nvox)
{
  size_t i,j;
  int b,r,c;

  /* read image dims */
  int nrows=0,ncols=0,ntimesteps=0,nt=0;
  int nslices = VAttrListNumImages(list);
  VImage *src = VAttrListGetImages(list,nslices);
  VImageDimensions(src,nslices,&ntimesteps,&nrows,&ncols);

  if (VImageNRows(mask) != nrows)
    VError(" inconsistent image dims, mask has %d rows, image has %d rows",VImageNRows(mask),nrows);
  if (VImageNColumns(mask) != ncols)
    VError(" inconsistent image dims, mask has %d columns, image has %d columns",VImageNColumns(mask),ncols);
  if (VImageNBands(mask) != nslices)
    VError(" inconsistent image dims, mask has %d slices, image has %d slices",VImageNBands(mask),nslices);

   /* get time steps to include */
  if (length < 1) length = ntimesteps-1;
  int last = first + length;
  if (last >= ntimesteps) last = ntimesteps-1;
  if (first < 0) first = 0;
  nt = last - first + 1;
  if (nt < 2) VError(" not enough timesteps, nt= %d",nt);
  fprintf(stderr," ntimesteps: %d, first= %d, last= %d, nt= %d\n",
	  (int)ntimesteps,(int)first,(int)last,(int)nt);

  /* alloc and fill  matrix X */
  size_t dim = nt;
  if (type == 0 || type == 7) dim = 2*nt;
  gsl_matrix_float *X = gsl_matrix_float_calloc(nvox,dim);
  for (i=0; i<nvox; i++) {

    b = VPixel(map,0,0,i,VInteger);
    r = VPixel(map,0,1,i,VInteger);
    c = VPixel(map,0,2,i,VInteger);

    float *ptr = gsl_matrix_float_ptr(X,i,0);
    int k;
    j = 0;
    for (k=first; k<=last; k++) {
      if (j >= X->size2) VError(" j= %d %d",j,X->size2);
      if (k >= VImageNBands(src[b])) VError(" k= %d %d",k, VImageNBands(src[b]));
      *ptr++ = (float) VGetPixel(src[b],k,r,c);
      j++;
    }
  }

  /* normalize */
  NormalizeData(X,nt);

  return X;
}

ext: c

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
#include <gsl/gsl_math.h>
#include <gsl/gsl_cblas.h>
#include "cblas.h"

void
cblas_sscal (const int N, const float alpha, float *X, const int incX)
{
#define BASE float
#include "source_scal_r.h"
#undef BASE
}

max_stars_repo_licenses: Net-SNMP Xnet

max_issues_repo_licenses: Net-SNMP Xnet

max_forks_repo_licenses: Net-SNMP Xnet


text:
// Sammlung verwendeter Funktionen
#include "global.h"
#include "steppmethods.h"
#include "targets.h"
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_sf_gamma.h>
#include <math.h>
#include <stdio.h>
#include "utility.h"
#include <string.h>
#include <time.h>
#include <omp.h>

double Norm_Diff(double *vec1,double *vec2, int dim, int p) //berechnet norm von || vec1-vec2 ||_p
{
	double result = 0.0d;
	double temp;
	if (p == 1)
	{
		for (int i = 0; i < DIM; i++)
		{
			result += fabs(vec1[i] - vec2[i]);
		}
		return(result);
	}
	if (p == 2)
	{
		for (int i = 0; i < DIM; i++)
		{
			result += pow(vec1[i] - vec2[i],2.0);
		}
		return(sqrt(result));
	}
	if (p == 0)
	{
		for (int i = 0; i < DIM; i++)
		{
			temp = fabs(vec1[i] - vec2[i]);
			if (result < temp)
			{
				result = temp;
			}
		}
		return(result);
	}
	return (0.0f);
}

void Lattice_Setup(double **positions,int  numbertoinf,int dim, double latticespacing) // puts coordinates into lattice matrix of dimension numbertoinf**dim, dim
{
	printf("setze Gitter auf \n");
	int m = -numbertoinf + 1;
	int l = -numbertoinf + 1;
	int count = 0;
	while(count<ipow(numbertoinf,dim))
	{
		if (m> (numbertoinf - 1) ) // if over line, restart
		{
			m = -numbertoinf +1;
			l = l + 2;
		}
		positions[count][0] = m * latticespacing/2;
		positions[count][1] = l * latticespacing/2;
		m = m + 2;
		count = count +1;
	}
	printf("Gitter aufgesetzt \n");
}

int ipow(int base, int exp) // I^p for integer i
{
	int result = 1;
	if(		(base < 0)\
		||	(exp < 0))
	{
		printf("wrong base or exponent, smaller zero!\n");
		exit;
	}
	while (exp)
	{
		if (exp & 1)
			result *= base;
		exp >>= 1;
		base *= base;
	}
	return result;
}

int Kahan_Sum(int N, double input[], double *ans){  // Add elements via Kahan summation, code is pretty much Wikipedia copy paste
	int i;
	double y,t,sum = 0.0d;
    double c = 0.0d;                 // A running compensation for lost low-order bits.
    for (i=0; i<N; i++){
        y = input[i] - c;	    // So far, so good: c is zero.
        t = sum + y;         	// Alas, sum is big, y small, so low-order digits of y are lost.
        c = (t - sum) - y;       // (t - sum) cancels the high-order part of y; subtracting y recovers negative (low part of y)
        sum = t;                 // Algebraically, c should always be zero. Beware overly-aggressive optimizing compilers!
    }                      		// Next time around, the lost low part will be added to y in a fresh attempt.
    *ans  = (sum);
    return 0;
}

int Sign_Sum(int N, double *input, double *ans){  // Add seperated by sign
	int i;
	double sum = 0.0, negsum = 0.0;
	for (i=0; i<N; i++){
		if(input[i]<0.0) negsum += input[i];
		if(input[i]>=0.0) sum += input[i];
	}
	*ans = sum + negsum;
	return 0;
}

int Neumaier_Sum(int N, double input[N], double *ans){  // Add elements of input by Neumaier summation
	int i;
	double t,sum = input[0];
    double c = 0.0;                 // A running compensation for lost low-order bits.
    for (i=1; i<N; i++){
    	t = sum + input[i];
    	if (fabs(sum)>=fabs(input[i]))
    	{
        	c += (sum - t) + input[i]; // If sum is bigger, low-order digits of input[i] are lost.
        } else
        {
            c += (input[i] - t) + sum; // Else low-order digits of sum are lost
        }
        sum = t;
    }                      // Next time around, the lost low part will be added to y in a fresh attempt.
    *ans = (sum + c);         // Correction only applied once in the very end
    return 0;
}

void deleteSpaces(char src[], char dst[]){      // deletes spaces from string
   // src is supposed to be zero ended
   // dst is supposed to be large enough to hold src
	int s, d=0;
	for (s=0; src[s] != 0; s++)
		if (src[s] != ' ')
		{
			dst[d] = src[s];
			if (dst[d] == ':')
			{
				dst[d] = '_';
			}
			d++;
		}

		dst[d] = 0;
	}

void vec_zero(double *v,int  N){ // zero target vector of double
	int i;
	for (i= 0 ; i < N; i++)
	{
		v[i] = 0.0;
	}

}
void vec_zero_i(int *v,int  N){ // zero target vector of ints
	int i;
	for (i= 0 ; i < N; i++)
	{
		v[i] = 0;
	}

}
void zero_Impuls(double *v)
{ 								// zero target vector of momenta
	double P[OSSZI + 1];
	double sum;
	double OS = OSSZI * 1.0;
	for (int j = 0; j < DIM; j++)
	{
		for (int i = 0; i < OSSZI; i++)
		{
			P[i] = v[DIM + i * DIM + j];
		}
		P[OSSZI] = v[j] * mass;
		Sum_Method(OSSZI+1,P,&sum);
		for (int i = 0; i < OSSZI; i++)
		{
			v[DIM + i * DIM + j] -= sum/( OS  + 1.0);
		}
		v[j] -= sum/( OS  + 1.0) /mass;
	}
}

int Update_Lattice_Position(int * pos, const  double *vec, const double lattice_spacing){ // updatet Int-Vector pos when leaving square lattice cell

	int ret_val = 0;			// return zero if no update neccessary, -1 if out of bounds
	for (int j = 0; j< DIM; j++)
	{
		int trigger = 1;
		while(trigger)
		{
			if (vec[j]< pos[j]*lattice_spacing - lattice_spacing/2.0)
			{
				pos[j]-=1;
				trigger = 1;
				ret_val = 1;
			}else
			{
				trigger = 0;
			}
			if (vec[j]> pos[j]*lattice_spacing + lattice_spacing/2.0)
			{
				pos[j]+=1;
				trigger = 1;
				ret_val = 1;
			}else
			{
				trigger = 0;
			}
		}
	}
	if ((abs(pos[0])> NUMBER_TO_INF) || (abs(pos[1])> NUMBER_TO_INF))
	{
		printf(" \n particle out of bounds");
		FLAG_DUMP_VALUE = 1;											// error occured
		return(-1);
	}
	return(ret_val);
}

void VVerlet(int N, double *y, double **ans, double *t,							// Velocity Verlet for Start y, return ans[LengthT][ORDER] to times T
				void (*derivmethod) (double *y, double *ans, double t, int N)){
	int i,j;
	for(i=0; i< N; i++){
	ans[0][i] = y[i];  				// give initial vec to answer
	}
	double x[N/2], v[N/2], a[N/2], z[N];
	vec_zero(x,N/2);
	vec_zero(v,N/2);
	vec_zero(a,N/2);
	vec_zero(z,N);
	int current_index = 0;
	double current_t = t[current_index];
	derivmethod(y, z, current_t, N); 				// get first derivative
	for (j = 0; j < DIM; j++ )						// bring y-p  Form to x-v-a
	{
		x[j] = y[j];
		v[j] = y[DIM + j]/mass;
		a[j] = z[DIM + j]/mass;
		#pragma omp parallel for
		for (i= 0 ; i < OSSZI; i++)
		{
					x[DIM + j + i*DIM] = y[2*DIM + j + i * 2*DIM];
					v[DIM + j + i*DIM] = y[3*DIM + j + i * 2*DIM]/massq[i];
					a[DIM + j + i*DIM] = z[3*DIM + j + i * 2*DIM]/massq[i];
		}
	}

	while(current_index < LENGTH_T-1)
	{
		while(fabs(current_t - t[current_index + 1])> MAXSTEPSIZE)
		{ 	//integrate till space till next time smaller than stepsize
			Stepper_Method(N, x, v, a, &current_t, derivmethod);
		}
												// get to y-Form
		current_index = current_index +1;
		for (j = 0; j < DIM; j++ )
		{
			ans[current_index][DIM + j] = mass * v[j];
			ans[current_index][j] 	    = x[j];
			#pragma omp parallel for
			for (i= 0 ; i < OSSZI; i++)
			{
						ans[current_index][2*DIM + j + i * 2*DIM] = x[DIM +j+  i*DIM];
						ans[current_index][3*DIM + j + i * 2*DIM] = v[DIM +j+  i*DIM] * massq[i]; // übergib SChritt an antwort
			}
		}
	}

}

void VVerlet_parallel(const int N, const  double *y, double **ans, double *t,							// Velocity Verlet for Start y, return ans[LengthT][ORDER] to times T
				void (*derivmethod) (double *y, double *ans, double t, int N)){
	int i,j;
	for(i=0; i< N; i++){
	ans[0][i] = y[i];  				// give initial vec to answer
	}
	TARGET_FLAG = 1;
	double x[N/2], v[N/2], a[N/2], z[N];
	vec_zero(x,N/2);
	vec_zero(v,N/2);
	vec_zero(a,N/2);
	vec_zero(z,N);
	int current_index = 0;
	double current_t = t[current_index];
	derivmethod(y, z, current_t, N); 						// giv f'
	int steps_in_intervall = 0;						// counts steps done in intervall  t, t+ dt
	int max_steps_in_intervall = 100000;			// after X steps, assume error occured
	for (j = 0; j < DIM; j++ )						//  y-p  --> x-v-a form
	{
		x[j] = y[j];
		VIRTUAL_X[j] = x[j];
		v[j] = y[DIM + j]/mass;
		VIRTUAL_V[j] = v[j];
		a[j] = z[DIM + j]/mass;
	}

	#pragma omp parallel for
	for (i= 0 ; i < OSSZI; i++)
	{
		for (j = 0; j < DIM; j++ )
		{
			x[DIM + j + i*DIM] = y[2*DIM + j + i * 2*DIM];
			v[DIM + j + i*DIM] = y[3*DIM + j + i * 2*DIM]/massq[i];
			a[DIM + j + i*DIM] = z[3*DIM + j + i * 2*DIM]/massq[i];
		}
	}
	LASTSTEPSIZE = MAXSTEPSIZE;
	while(current_index < LENGTH_T-1)
	{
		while(fabs(current_t - t[current_index + 1])> MAXSTEPSIZE)			//integrate till space till next time smaller than stepsize
		{
			Stepper_Method(N, x, v, a, &current_t, derivmethod);
			steps_in_intervall++;
			if(steps_in_intervall > max_steps_in_intervall)
			{
				FLAG_DUMP_VALUE  = 1;										// too many steps, error occured
				printf("too many steps taken to advance one time point \n");
			}
			if (FLAG_DUMP_VALUE) 											// if error occured somewhere leave calculation
			{
				return;
			}
		}
		steps_in_intervall = 0; 				// reset stepcount
		current_index = current_index +1;		// get  back to y-form
		for (j = 0; j < DIM; j++ )
		{
			if(VIRTUAL_FLAG)
			{
				ans[current_index][DIM + j] = mass * VIRTUAL_V[j];
				ans[current_index][j] 	    = VIRTUAL_X[j];
			}else
			{
				ans[current_index][DIM + j] = mass * v[j];
				ans[current_index][j] 	    = x[j];
			}
		}
		#pragma omp parallel for
		for (i= 0 ; i < OSSZI; i++)
		{
			for (j = 0; j < DIM; j++ )
			{
				ans[current_index][2*DIM + j + i * 2*DIM] = x[DIM +j+  i*DIM];
				ans[current_index][3*DIM + j + i * 2*DIM] = v[DIM +j+  i*DIM] * massq[i] ; // give step to answer
			}
		}
	}
	if (TARGET_FLAG)		// if no contact, give simulation length as contact time
	{
		time_first_contact += TIME_END;
		TARGET_FLAG = 0;
	}
}

void VVerlet_parallel_burn(const int N, const  double *y, double **ans, double *t,							// Velocity Verlet für Start y, Ausgabe ans[LengthT][ORDER] zu Zeiten T
				void (*derivmethod) (double *y, double *ans, double t, int N))
// always involes free particle for burning samples
{
	int i,j;
	for(i=0; i< N; i++){
	ans[0][i] = y[i];  				// übergib startvektor an antwort
	}
	TARGET_FLAG = 1;
	double x[N/2], v[N/2], a[N/2], z[N];
	vec_zero(x,N/2);
	vec_zero(v,N/2);
	vec_zero(a,N/2);
	vec_zero(z,N);
	int current_index = 0;
	double current_t = t[current_index];
	derivmethod(y, z, current_t, N); 				// übergib erste Ableitung
	// bringe aus y-p y Form auf x-v-a
	for (j = 0; j < DIM; j++ )
	{
		x[j] = y[j];
		v[j] = y[DIM + j]/mass;
		a[j] = z[DIM + j]/mass;
	}
	#pragma omp parallel for
	for (i= 0 ; i < OSSZI; i++)
	{
		for (j = 0; j < DIM; j++ )
		{
			x[DIM + j + i*DIM] = y[2*DIM + j + i * 2*DIM];
			v[DIM + j + i*DIM] = y[3*DIM + j + i * 2*DIM]/massq[i];
			a[DIM + j + i*DIM] = z[3*DIM + j + i * 2*DIM]/massq[i];
		}
	}

	while(current_index < LENGTH_T-1)
	{
		while(fabs(current_t - t[current_index + 1])> MAXSTEPSIZE)
		{ 	//integriere bis abstand zu  nächstem t kleiner schrittweite
			VVerlet_Step_deriv(N, x, v, a, &current_t, derivmethod);
			//printf ("%fx      %ft  \n", x[1], current_t);
		}
												// rechne wieder auf y-Form
		current_index = current_index +1;
		for (j = 0; j < DIM; j++ )
		{
			ans[current_index][DIM + j] = mass * v[j];
			ans[current_index][j] 	    = x[j];
		}
		#pragma omp parallel for
		for (i= 0 ; i < OSSZI; i++)
		{
			for (j = 0; j < DIM; j++ )
			{
				ans[current_index][2*DIM + j + i * 2*DIM] = x[DIM +j+  i*DIM];
				ans[current_index][3*DIM + j + i * 2*DIM] = v[DIM +j+  i*DIM] * massq[i]; // übergib SChritt an antwort
			}
		}
	}
	if (TARGET_FLAG)		// FAlls kein erstkontakt addiere gesamtzeit
	{
		time_first_contact += TIME_END;
		TARGET_FLAG = 0;
	}
}

void Bath_Setup(double *y, char *inputlabel, gsl_rng * r){
	int a,b,i;
	double u,P[OSSZI], sum=0.0;

	if (strcmp(inputlabel,"ZWANZIG")==0)
	{
		for (b = 0; b < DIM; b++)
		{
			y[b] = 0.0;					//massives particle seperatly
			y[DIM+b] = 0.0;
			for (a = 0; a < OSSZI; a++)
			{
				u = gsl_ran_gaussian(r, 1.0);
				y[2*DIM + a * 2*DIM + b] = sqrt(KBOLTZ * TEMP / massq[a]) / ommega[a] * u; 	//q = (k*T/m)^0.5/ommega * N(0,1)
				u = gsl_ran_gaussian(r, 1.0);
				y[3*DIM + a * 2*DIM + b] = sqrt(KBOLTZ * TEMP * massq[a]) * u;			//p_q = (k*T*m)^0.5 * N(0,1)
			}
		}
		return;
	}
	if (strcmp(inputlabel,"ZWANZIG_THERMALIZED")==0)
	{
		for (b = 0; b < DIM; b++)
		{
			y[b] = 0.0;					//massives particle seperatly
			u = gsl_ran_gaussian(r, 1.0);
			y[DIM+b] =  sqrt(KBOLTZ * TEMP * mass) * u;			//p_q = (k*T*m)^0.5 * N(0,1)
			for (a = 0; a < OSSZI; a++)
			{
				u = gsl_ran_gaussian(r, 1.0);
				y[2*DIM + a * 2*DIM + b] = sqrt(KBOLTZ * TEMP / massq[a]) / ommega[a] * u; 	//q = (k*T/m)^0.5/ommega * N(0,1)
				u = gsl_ran_gaussian(r, 1.0);
				y[3*DIM + a * 2*DIM + b] = sqrt(KBOLTZ * TEMP * massq[a]) * u;			//p_q = (k*T*m)^0.5 * N(0,1)
			}
		}
		return;
	}
	printf("FALSCHES BADLABEL");
	exit(1);
}

void Ommega_Setup(gsl_rng * r)
{

 															// draw uniformly nonrandom
	double nr_osszillators = 1.0 * OSSZI; 					// cast Osszi to double
	double higher_cutoff = 1.0/(0.5 * TIME_STEPS); 			//cleanly resolve smallest steps
	double lower_cutoff = 1.0/(1.6 * TIME_END); 			//cleanly resolve till end time

	higher_cutoff = 5.0* pow( (GAMMA/mass), 1.0/(2.0-ALPHA));
	lower_cutoff = 5E-5;
	double ommega_step = (higher_cutoff - lower_cutoff)/(nr_osszillators);
	for (int i = 0; i < OSSZI; i++)
	{
		ommega[i] = lower_cutoff + i * ommega_step;
	}
	printf("\nnonrandom bath ");
	printf ("Ommega_Max =   %1.3e  ", ommega[OSSZI-1] );
	printf ("Ommega_Min =   %1.3e\n", ommega[0] );
}

void Gamma_Setup(double *gamma){ 						//Set coupling for Exponent Alpha so that MSD s^2 prop to
														//t^ALPHA,  ALPHA = 1 (classical diff) seperately
														// put constants in gamma
	int i;
	double TOL = 1e-3;
	double nr_osz = (double) OSSZI;
	double delta_ommega = ommega[1]-ommega[0];
	if (! (KUPF_BINARY))
	{
		if (fabs(ALPHA-1.0)<TOL){  // Alpha = 1
		for(i = 0; i < OSSZI; i++)
		{
			gamma[i] = sqrt(GAMMA) * ommega[i] * sqrt(delta_ommega);
		}
		if (OUTPUT_FLAG)
		{
			printf("ALPHA rounded to 1 \n");
		}
		}
		if (fabs(ALPHA-1.0)>TOL){
			for(i = 0; i < OSSZI; i++)
			{
				gamma[i] = sqrt(GAMMA) * pow(ommega[i], 0.5 + ALPHA * 0.5)\
				* sqrt( 2/sqrt(M_PI * 2) * sin(M_PI * ALPHA/2.0) * gsl_sf_gamma(1 - ALPHA)) * sqrt(delta_ommega);
			}
		}
		if (OUTPUT_FLAG)
		{
			printf("Konstanten aufgesetzt für Alpha =  %f \n " , ALPHA);
		}
	}else
	{
		if (fabs(ALPHA-1.0)<TOL)	// Alpha = 1
		{
			for(i = 0; i < OSSZI; i++)
			{

				gamma[i] = 1.0 / sqrt(2.0 * M_PI) * GAMMA * (delta_ommega);
			}
			printf("Konstanten aufgesetzt für KUPF_A rounded to 1  \n " );
		}else
		{
			for(i = 0; i < OSSZI; i++)
			{

				gamma[i] = 1.0 / sqrt(2.0 * M_PI) * GAMMA  * gsl_sf_gamma(1 - ALPHA)\
				* sin(M_PI * ALPHA/2.0) * pow(ommega[i], ALPHA - 1.0) * (delta_ommega);
			}
			printf("Konstanten aufgesetzt für KUPF_A = %2.2f  \n ", ALPHA );
		}

	}
}

void deriv(double *yin, double *ans, double t, int N){ 		// Hamiltonian derivative

	int i,j;
	double P[OSSZI], Q[OSSZI], temp_vec[OSSZI];
	double sum;

	for(j = 0; j < DIM; j++)
	{
		#pragma omp parallel for
		for(i = 0; i < OSSZI; i++)
		{
			Q[i] = yin[(2 + i * 2)  * DIM  + j ];  //get P und Q for easier indexing
			P[i] = yin[(3 + i * 2)  * DIM  + j ];
		}
		sum = 0.0;
		#pragma omp parallel for
		for (i = 0 ; i < OSSZI; i++)
		{ // Sum
			if(KUPF_BINARY)
			{
				temp_vec[i] = -coupling[i] * (Q[i] - yin[j]);

			}else
			{
				temp_vec[i] = coupling[i] * Q[i] - pow(coupling[i],2.0)/pow(ommega[i],2.0) * yin[j];
			}

		}
		Sum_Method(OSSZI, temp_vec, &sum);
		ans[j] = yin[DIM+ j]/mass;  // dx/dx = p/m
		ans[DIM+j] = sum  ;   	//p-deriv
		#pragma omp parallel for
		for (i = 0 ; i < OSSZI; i++)
		{
			ans[(2 + i * 2)  * DIM  + j ] = P[i]/massq[i];      		// dq/dt = p_q/m_q
			if(KUPF_BINARY)
			{
				ans[(3 + i * 2)  * DIM  + j ] =  coupling[i] * (Q[i] - yin[j]); 		// dp_q/dt = -w^2 * q+ gamma * x

			}else
			{
				ans[(3 + i * 2)  * DIM  + j ] = - ommega[i] * ommega[i] * Q[i]  + coupling[i] * yin[j]; 		// dp_q/dt = -w^2 * q+ gamma * x

			}

		}
	}
}

void deriv_parallel(const double *yin, double *ans, double t, int N)		// deriv with i-j loops reversed for parralelization
{ 		// Ableitungsfunktion
	int i,j;
	double temp_vec[OSSZI];
	double sum;

	#pragma omp parallel for 		// Values for Bath particles parralized
	for (i = 0 ; i < OSSZI; i++)
	{
		for(j = 0; j < DIM; j++)
		{
			ans[(2 + i * 2)  * DIM  + j ] = yin[(3 + i * 2)  * DIM  + j ]/massq[i];      		// dq/dt = p_q/m_q
			if(KUPF_BINARY)
			{
				ans[(3 + i * 2)  * DIM  + j ] =  -coupling[i] *  (yin[(2 + i * 2)  * DIM  + j] - yin[j]); 		// dp_q/dt = -w^2 * q+ gamma * x
			}else
			{
				ans[(3 + i * 2)  * DIM  + j ] = - ommega[i] * ommega[i] * yin[(2 + i * 2)  * DIM  + j]\
										  + coupling[i] * yin[j]; 		// dp_q/dt = -w^2 * q+ gamma * x
			}
		}
	}

	for(j = 0; j < DIM; j++)		// Values for heavy particle
	{

		sum = 0.0;
		#pragma omp parallel for
		for (i = 0 ; i < OSSZI; i++)
		{ // berechne die Summe

			if(KUPF_BINARY)
			{
				temp_vec[i] =  coupling[i] *  (yin[(2 + i * 2)  * DIM  + j] - yin[j]); 		// dp_q/dt = -w^2 * q+ gamma * x
			}else
			{
				temp_vec[i] = coupling[i] * yin[(2 + i * 2)  * DIM  + j]\
						 - pow(coupling[i],2.0)/pow(ommega[i],2.0) * yin[j];  // for Kahan summation : vec gamma^2/ommega^2
			}
		}
		Sum_Method(OSSZI, temp_vec, &sum);
		ans[j] = yin[DIM+ j]/mass;  // dx/dx = p/m
		ans[DIM+j] = sum  ;   	//p-deriv
	}
}


text:
#include <stdlib.h>
#include <stdio.h>
#include <unistd.h>
#include <string.h>
#include <stdint.h>
#include <math.h>
#include <assert.h>
#include <gsl/gsl_math.h>

#include "com3.h"

// #define debugmode 1
int unit_tests(void);

const int64_t comr = 2; // Radius of com filter in pixels
// The filter will be (2*comr+1)^3 pixels

int pointInDomain(int64_t * restrict Point, size_t * Domain, int nDim)
  // Check if a point is in the domain of a 3D image
{
  for(int kk = 0; kk<nDim; kk++)
    if(Point[kk]<0 || Point[kk]>= (int64_t) Domain[kk])
      return 0;

  return 1;
}

int checkBounds(int64_t * restrict D,
    const size_t M, const size_t N, const size_t P)
{
  if(D[0]>=comr && D[0]+comr<(int64_t) M &&
      D[1]>=comr && D[1]+comr<(int64_t) N &&
      D[2]>=comr && D[2]+comr<(int64_t) P)
    return 1; // ok

  return 0;
}

void com3_localw(double * restrict V, double * restrict W, const size_t M, const size_t N, const size_t P,
    int64_t * restrict D, double * restrict C)
{

#ifdef debugmode
  printf("D: [%u %u %u]\n", D[0], D[1], D[2]);
  printf("M: %lu N: %lu P: %lu\n", M, N, P);
#endif

  double sum = 0;
  double dx = 0;
  double dy = 0;
  double dz = 0;

#ifdef debugmode
  printf("sum: %f\n", sum);
#endif

  if(checkBounds(D, M, N, P))
  {

    double val = V[D[0] + D[1]*M + D[2]*M*N];

    for(int kk = -comr; kk<=comr; kk++) {
      for(int ll = -comr; ll<=comr; ll++) {
        for(int mm = -comr; mm<=comr; mm++) {

          size_t pos = (D[0]+kk) + (D[1]+ll)*M + (D[2]+mm)*M*N;
          //  printf("pos: %lu V[pos]: %f\n", pos, V[pos]);
          dx += kk*V[pos]*val/W[pos];
          dy += ll*V[pos]*val/W[pos];
          dz += mm*V[pos]*val/W[pos];

          sum += V[pos];
        } } }

    //printf("sum: %f\n", sum);
  }

  if(sum>0)
  {
#ifdef debugmode
    printf("sum: %f (%f, %f, %f)\n", sum, dx, dy, dz);
#endif
    C[0] = D[0] + dx/sum;
    C[1] = D[1] + dy/sum;
    C[2] = D[2] + dz/sum;
  }
  else
  {
#ifdef debugmode
    printf("No com calculation\n");
#endif
    C[0] = D[0];
    C[1] = D[1];
    C[2] = D[2];
  }
}

void com3_local(double * restrict V, const size_t M, const size_t N, const size_t P,
    int64_t * restrict D, double * restrict C)
{

#ifdef debugmode
  printf("D: [%u %u %u]\n", D[0], D[1], D[2]);
  printf("M: %lu N: %lu P: %lu\n", M, N, P);
#endif

  double sum = 0;
  double dx = 0;
  double dy = 0;
  double dz = 0;

#ifdef debugmode
  printf("sum: %f\n", sum);
#endif

  if(checkBounds(D, M, N, P))
  {

    for(int kk = -comr; kk<=comr; kk++) {
      for(int ll = -comr; ll<=comr; ll++) {
        for(int mm = -comr; mm<=comr; mm++) {

          size_t pos = (D[0]+kk) + (D[1]+ll)*M + (D[2]+mm)*M*N;
          //  printf("pos: %lu V[pos]: %f\n", pos, V[pos]);
          dx += kk*V[pos];
          dy += ll*V[pos];
          dz += mm*V[pos];

          sum += V[pos];
        } } }

    //printf("sum: %f\n", sum);
  }

  if(sum>0)
  {
#ifdef debugmode
    printf("sum: %f (%f, %f, %f)\n", sum, dx, dy, dz);
#endif
    C[0] = D[0] + dx/sum;
    C[1] = D[1] + dy/sum;
    C[2] = D[2] + dz/sum;
  }
  else
  {
#ifdef debugmode
    printf("No com calculation\n");
#endif
    C[0] = D[0];
    C[1] = D[1];
    C[2] = D[2];
  }
}

void setLmax(double *V, double * W, size_t * Domain, int64_t * Dot)
  // Adds V(dot) to all W within comr of dot
  // uses global value comr
{
  uint32_t m = Dot[0];
  uint32_t n = Dot[1];
  uint32_t p = Dot[2];

  size_t M = Domain[0];
  size_t N = Domain[1];
  size_t P = Domain[2];

  double val = V[m + n*M + p*M*N];

  for(uint32_t     pp = GSL_MAX(p-comr,0) ; pp<= GSL_MIN(p+comr, P-1) ; pp++) {
    for(uint32_t   nn = GSL_MAX(n-comr,0) ; nn<= GSL_MIN(n+comr, N-1) ; nn++) {
      for(uint32_t mm = GSL_MAX(m-comr,0) ; mm<= GSL_MIN(m+comr, M-1) ; mm++) {
        W[mm + nn*M + pp*M*N] += val;
      }
    }
  }
}

void setWeights(double * V, double * W, size_t M, size_t N, size_t P, double * D, size_t nD)
{
  size_t Domain[] = {M, N, P};
  int64_t Dot[] = {0,0,0};
  for(size_t kk = 0; kk<nD; kk++) {
    Dot[0] = D[kk*3 + 0];
    Dot[1] = D[kk*3 + 1];
    Dot[2] = D[kk*3 + 2];

    setLmax(V, W, Domain, Dot);
  }
}

void com3(double * V, size_t M, size_t N, size_t P,
    double * D, double * C, size_t L, int weighted) {
  /* Centre of mass
   *
   * V: MxNxP image
   * P 3xL list of dots
   * C 3xL list of fitted dots
   */

  int64_t Dround[] = {0,0,0};
  size_t Domain[] = {M, N, P};

  double * W;
  if(weighted == 1)
    // Set up the weighting matrix
  {
    W = (double *) calloc(M*N*P, sizeof(double));
    for(size_t kk = 0; kk<L; kk++)
    {
      Dround[0] = nearbyint(D[kk*3]-1);
      Dround[1] = nearbyint(D[kk*3+1]-1);
      Dround[2] = nearbyint(D[kk*3+2]-1);
      if(pointInDomain(Dround, Domain, 3))

        setLmax(V, W, Domain, Dround);
    }
  }

  if(weighted == 1)
    for(size_t kk = 0; kk<L; kk++)
    {
      Dround[0] = nearbyint(D[kk*3]-1);
      Dround[1] = nearbyint(D[kk*3+1]-1);
      Dround[2] = nearbyint(D[kk*3+2]-1);

      //    printf("%d %d %d\n", Dround[0], Dround[1], Dround[2]);

      com3_localw(V, W, M, N, P, Dround, C+kk*3);
    }

  if(weighted ==0)
    for(size_t kk = 0; kk<L; kk++)
    {
      Dround[0] = nearbyint(D[kk*3]-1);
      Dround[1] = nearbyint(D[kk*3+1]-1);
      Dround[2] = nearbyint(D[kk*3+2]-1);

      //    printf("%d %d %d\n", Dround[0], Dround[1], Dround[2]);

      com3_local(V, M, N, P, Dround, C+kk*3);
    }

  if(weighted == 1){

    // This was used for debugging
    //for(size_t xx = 0; xx<M*N*P; xx++)
    //  V[xx] = W[xx];

    free(W);
  }
}

#ifdef standalone

int unit_tests() {

  // Size of image
  uint32_t M = 100;
  uint32_t N = 100;
  uint32_t P = 100;
  uint32_t L = 2; // number of dots

  double * V = malloc(M*N*P*sizeof(double));
  double * D = malloc(L*3*sizeof(double));
  double * C = malloc(L*3*sizeof(double));

  memset(V, 0, M*N*P*sizeof(double));
  memset(C, 0, 3*L*sizeof(double));

  for(size_t kk = 0; kk<M*N*P; kk++)
    V[kk] = 0;

  D[0] = 11; D[1] = 12; D[2] = 13;
  D[3] = 0; D[4] = 15; D[5] = 16;

  size_t pos = D[0] + M*D[1] + M*N*D[2];
  V[pos] = 3;
  //  V[pos+1] = 1;
  //  V[pos + M] = 1;
  V[pos + M*N] = 1;

  // Ordinary
  com3(V, M, N, P,
      D, C, L, 0);

  // Weighted
  com3(V, M, N, P,
      D, C, L, 1);

  for(uint32_t kk=0; kk<L; kk++) {
    printf("%d [%f %f %f] -> ", kk, D[3*kk], D[3*kk+1], D[3*kk+2]);
    printf(" [%f %f %f]\n",         C[3*kk], C[3*kk+1], C[3*kk+2]);
  }

  free(C);
  free(D);
  free(V);

  return 0;
}

int main(int argc, char ** argv)
{
  printf("%s\n", argv[0]);
  if(argc == 1)
    return unit_tests();
}
#endif

max_stars_repo_name: elgw/dotter

max_issues_repo_name: elgw/dotter

max_forks_repo_name: elgw/dotter


text:
/* -----------------------------------------------------------------------------
 * Copyright 2021 Jonathan Haigh
 * SPDX-License-Identifier: MIT
 * ---------------------------------------------------------------------------*/

#ifndef SQ_INCLUDE_GUARD_parser_Parser_inl_h_
#define SQ_INCLUDE_GUARD_parser_Parser_inl_h_

#include <charconv>

#include "core/ASSERT.h"
#include "core/errors.h"

#include <gsl/gsl>

namespace sq::parser {

template <std::integral Int> std::optional<Int> Parser::parse_integer() {
  const auto opt_token = accept_token(TokenKind::Integer);
  if (!opt_token) {
    return std::nullopt;
  }
  const auto &token = opt_token.value();
  const auto str_view = token.view();
  ASSERT(!str_view.empty());
  const auto *begin = str_view.data();
  // We can't really avoid using pointer arithmetic when using
  // std::from_chars - it requires a const char* to indicate the end of the
  // string, but we can only reasonably make one using pointer arithmetic.
  // In particular, std::string, std::string_view, std::span etc. only
  // provide iterator versions of end(), but we need a pointer.
  // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-pointer-arithmetic)
  const auto *end = begin + str_view.size();
  Int value = 0;
  const auto [ptr, ec] = std::from_chars(begin, end, value, 10);
  if (ec == std::errc::result_out_of_range) {
    auto ss = std::ostringstream{};
    ss << "integer " << str_view << " does not fit in required type; "
       << "must be in the closed interval [" << std::numeric_limits<Int>::min()
       << ", " << std::numeric_limits<Int>::max() << "]";
    throw OutOfRangeError{token, ss.str()};
  }
  ASSERT(ec == std::errc{});
  ASSERT(ptr == end);
  return value;
}

} // namespace sq::parser

#endif // SQ_INCLUDE_GUARD_parser_Parser_inl_h_

max_stars_repo_name: jonathanhaigh/sq

max_issues_repo_name: jonathanhaigh/sq

max_forks_repo_name: jonathanhaigh/sq


text:
/*
 *  emu_noh.c
 *
 *
 *  Created by Earl Lawrence on 9/17/09.
 *  Update 11/30/2012
 *
 *  This program was prepared by Los Alamos National Security, LLC at Los Alamos National Laboratory (LANL)
 *  under contract No. DE-AC52-06NA25396 with the U.S. Department of Energy (DOE). All rights in the program
 *  are reserved by the DOE and Los Alamos National Security, LLC.  Permission is granted to the public to
 *  copy and use this software without charge, provided that this Notice and any statement of authorship are
 *  reproduced on all copies.  Neither the U.S. Government nor LANS makes any warranty, express or implied,
 *  or assumes any liability or responsibility for the use of this software.
 *
 *
 */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_errno.h>

#include "pcbasis_noh.h"
#include "design_noh.h"
#include "pcweights_noh.h"
#include "corrlengths_noh.h"
#include "precisions_noh.h"
#include "meansd.h"
#include "kemu.h"
#include "ksim.h"

#ifndef M_PI
#define M_PI           3.14159265358979323846
#endif

// Sizes of stuff and number of redshifts
static int m=37, neta=5500, p=5, peta=6, rs=11, nsim=582;
// Kriging basis computed by emuInit, sizes should be peta and m
static double KrigBasis[6][37];

// Initialization function that computes the Kriging basis
void emuInit_noh() {
    int i,j,k,l;
    double cov;
    gsl_matrix *SigmaSim = gsl_matrix_alloc(m,m);
    gsl_vector *b = gsl_vector_alloc(m);

    // Do these one principal component at a time
    for(i=0; i<peta; i++) {
        // Fill in the covariance matrix for the principals components
        // Also make a gsl_vector with the weights.
        for(j=0; j<m; j++) {
            // Diagonal
            gsl_matrix_set(SigmaSim, j, j, (1.0/lamz[i]) + (1.0/lamws[i]));
            // Off-diagonals
            for(k=0; k<j; k++) {
                // Compute the covariance
                cov = 0.0;
                for(l=0; l<p; l++) {
                    cov -= beta[i][l]*pow(x[j][l]-x[k][l], 2.0);
                }
                cov = exp(cov)/lamz[i];
                gsl_matrix_set(SigmaSim, j, k, cov);
                gsl_matrix_set(SigmaSim, k, j, cov);
            } // for(k=0; k<j; k++)
            gsl_vector_set(b, j, w[i][j]);
        } // for(j=0; j<m; j++)

        // Cholesky and solve
        gsl_linalg_cholesky_decomp(SigmaSim);
        gsl_linalg_cholesky_svx(SigmaSim, b);

        // Copy into the Kriging Basis
        for(j=0; j<m; j++) {
            KrigBasis[i][j] = gsl_vector_get(b, j);
        }
    } // for(i=0; i<peta; i++)
    gsl_matrix_free(SigmaSim);
    gsl_vector_free(b);
}

// The actual emulation
// Cosmological parameters, placeholder for the output, type of output
void emu_noh(double *xstar, double *ystar, int *outtype) {
    static int inited=0;
    int i, j, k;
    double wstar[peta], Sigmastar[peta][m], ystaremu[neta], ystar_allz[rs*nsim], logc;
    double xstarstd[p];
    double zemu[rs], ybyz[rs];
    FILE *fp;
    // Interpolation stuff for k and then z
    gsl_spline *lininterp_k = gsl_spline_alloc(gsl_interp_linear, neta/rs);
    gsl_spline *lininterp_z = gsl_spline_alloc(gsl_interp_linear, rs);
    gsl_interp_accel *accel = gsl_interp_accel_alloc();

    // Iinitialize if necessary
    if(inited==0) {
        emuInit_noh();
        inited=1;
    }

    // Check the inputs to make sure we're interpolating.
    for(i=0; i<p; i++) {
        if((xstar[i] < xmin[i]) || (xstar[i] > xmin[i]+xrange[i])) {
            //printf("The inputs are outside the domain of the emulator.\n");
            switch(i) {
                case 0:
                    printf("omega_b must be between %f and %f.\n", xmin[i], xmin[i]+xrange[i]);
                    break;
                case 1:
                    printf("omega_m must be between %f and %f.\n", xmin[i], xmin[i]+xrange[i]);
                    break;
                case 2:
                    printf("n_s must be between %f and %f.\n", xmin[i], xmin[i]+xrange[i]);
                    break;
                case 3:
                    printf("w must be between %f and %f.\n", xmin[i], xmin[i]+xrange[i]);
                    break;
                case 4:
                    printf("sigma_8 must be between %f and %f.\n", xmin[i], xmin[i]+xrange[i]);
                    break;
            }
            exit(1);
        }
    } // for(i=0; i<p; i++)

    // Check redshift to make sure we're interpolating
    if((xstar[p] < 0) || (xstar[p] > 4)) {
        //printf("The inputs are outside the domain of the emulator.\n");
        printf("z must between 0 and 4\n");
        exit(1);
    }

    // Standardize the inputs
    for(i=0; i<p; i++) {
        xstarstd[i] = (xstar[i] - xmin[i]) / xrange[i];
    }

    // Compute the covariances between the new input and sims for all PCs
    for(i=0; i<peta; i++) {
        for(j=0; j<m; j++) {
            logc = 0.0;
            for(k=0; k<p; k++) {
                logc -= beta[i][k]*pow(x[j][k]-xstarstd[k], 2.0);
            }
            Sigmastar[i][j] = exp(logc)/lamz[i];
        }
    }

    // Compute wstar, the predicted PC weights for the new input
    for(i=0; i<peta; i++) {
        wstar[i]=0.0;
        for(j=0; j<m; j++) {
            wstar[i] += Sigmastar[i][j] * KrigBasis[i][j];
        }
    }

    // Compute ystar, the new output
    for(i=0; i<neta; i++) {
        ystaremu[i] = 0.0;
        for(j=0; j<peta; j++) {
            ystaremu[i] += K[i][j]*wstar[j];
        }
        ystaremu[i] = ystaremu[i]*sd + mean[i];
    }

    //printf("Emulated\n");

    // Interpolate the emulated output onto the original domain.
    for(i=0; i<rs; i++) {
        gsl_spline_init(lininterp_k, kemu, &ystaremu[i*neta/rs], neta/rs);
        for(j=0; j<nsim; j++) {
            ystar_allz[i*nsim+j] = gsl_spline_eval(lininterp_k, ksim[j], accel);
        }
        gsl_interp_accel_reset(accel);
    }

    //printf("k interped\n");

    // Fill in the k values for the final output
    for(i=0; i<nsim; i++) {
        ystar[i] = ksim[i];
    }

    // Interpolate on to the desired redshift
    // The order needs to be reversed here.
    for(i=0; i<rs; i++) {
        zemu[i] = (1/aemu[rs-i-1]) - 1.0;
    }
    for(i=0; i<nsim; i++) {
        // Build an array with the values of y for a given value of z
        // Reverse the order
        for(j=0; j<rs; j++) {
            ybyz[rs-j-1] = ystar_allz[j*nsim+i];
        }
        gsl_spline_init(lininterp_z, zemu, ybyz, rs);
        ystar[nsim+i] = gsl_spline_eval(lininterp_z, xstar[p], accel);
        gsl_interp_accel_reset(accel);
    }

    //printf("z interped\n");

    switch(outtype[0]) {
        default:
	   for(j=0; j<nsim; j++) {
	     ystar[nsim+j] = pow(10.0,ystar[nsim+j]);
	   }

            break;
        case 1:
            // Transform to Delta^2
            for(j=0; j<nsim; j++) {
                ystar[nsim+j] = ystar[nsim+j] + 1.5*log10(ksim[j]);
                ystar[nsim+j] = pow(10.0,ystar[nsim+j]);
            }
            break;
        case 2:
            // Transform to P(k)
            for(j=0; j<nsim; j++) {
                ystar[nsim+j] = ystar[nsim+j] - 1.5*log10(ksim[j]);
                ystar[nsim+j] = pow(10.0,ystar[nsim+j])*2.0*M_PI*M_PI;
            }
            break;
    }

    // Free some stuff.  I always forget this.
    // Thanks to Tim Eifler for discovering it (apparently the hard way).
    gsl_spline_free(lininterp_z);
    gsl_spline_free(lininterp_k);
    gsl_interp_accel_free(accel);
}

// Linker function for use with Fortran
void emu_noh_(double *x, double *y, int *outtype) {
    void emu();
    emu(x,y,outtype);
}


text:
/**
 *
 * @file core_strdalg1.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Azzam Haidar
 * @date 2013-07-04
 * @generated s Tue Jan  7 11:44:50 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 * @ingroup CORE_float
 *
 *  CORE_strdalg1 is a part of the tridiagonal reduction algorithm (bulgechasing)
 *  It correspond to a local driver of the kernels that should be executed on a
 *  single core.
 *
 *******************************************************************************
 *
 * @param[in] n
 *          The order of the matrix A. n >= 0.
 *
 * @param[in] nb
 *          The size of the Bandwidth of the matrix A,
 *          which correspond to the tile size. nb >= 0.
 *
 * @param[in,out] A
 *          float array, dimension (lda,n)
 *          On entry, the (nb+1)-by-n band lower hermetian matrix to
 *          be reduced to tridiagonal.
 *          On exit, the diagonal and first subdiagonal of A are over-
 *          written by the corresponding elements of the tridiagonal.
 *
 * @param[in] lda
 *          (input) INTEGER
 *          The leading dimension of the array A.  LDA >= max(1,nb+1).
 *
 * @param[out] V
 *          float array, dimension (n) if wantz=0
 *          or ldv*Vblksiz*blkcnt if wantz>0.
 *          The scalar elementary reflectors are written in this
 *          array.
 *
 * @param[out] TAU
 *          float array, dimension (n) if wantz=0
 *          or Vblksiz*Vblksiz*blkcnt if wantz>0.
 *          The scalar factors of the elementary reflectors are written
 *          in this array.
 *
 * @param[in] Vblksiz
 *          Local parameter to Plasma. It correspond to the local bloccking
 *          of the applyQ2 used to apply the orthogonal matrix Q2.
 *
 * @param[in] wantz
 *          integer tobe 0 or 1. if wantz=0 the V and TAU are not stored on
 *          only they are kept for next step then overwritten.
 *
 * @param[in] i
 *          Integer that refer to the current sweep. (outer loop).
 *
 * @param[in] sweepid
 *          Integer that refer to the sweep to chase.(inner loop).
 *
 * @param[in] m
 *          Integer that refer to a sweep step, to ensure order dependencies.
 *
 * @param[in] grsiz
 *          Integer that refer to the size of a group.
 *          group mean the number of kernel that should be executed sequentially
 *          on the same core.
 *          group size is a trade-off between locality (cache reuse) and parallelism.
 *          a small group size increase parallelism while a large group size increase
 *          cache reuse.
 *
 * @param[in] work
 *          Workspace of size nb. Used by the core_ssbtype[123]cb.
 *
 ******************************************************************************/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_strdalg1 = PCORE_strdalg1
#define CORE_strdalg1 PCORE_strdalg1
#endif
void CORE_strdalg1( int n, int nb,
                    float *A, int lda,
                    float *V,
                    float *TAU,
                    int Vblksiz, int wantz,
                    int i, int sweepid, int m, int grsiz,
                    float *work)
{
    int k, shift=3;
    int myid, colpt, stind, edind, blklastind, stepercol;

    k = shift / grsiz;
    stepercol = (k*grsiz == shift) ? k : k+1;
    for (k = 0; k < grsiz; k++){
        myid = (i-sweepid)*(stepercol*grsiz) +(m-1)*grsiz + k+1;
        if(myid%2 ==0) {
            colpt      = (myid/2) * nb + 1 + sweepid - 1;
            stind      = colpt - nb + 1;
            edind      = min(colpt, n);
            blklastind = colpt;
        } else {
            colpt      = ((myid+1)/2)*nb + 1 +sweepid -1 ;
            stind      = colpt-nb+1;
            edind      = min(colpt,n);
            if( (stind>=edind-1) && (edind==n) )
                blklastind = n;
            else
                blklastind = 0;
        }

        if( myid == 1 )
           CORE_ssbtype1cb(n, nb, A, lda, V, TAU, stind-1, edind-1, sweepid-1, Vblksiz, wantz, work);
        else if(myid%2 == 0)
           CORE_ssbtype2cb(n, nb, A, lda, V, TAU, stind-1, edind-1, sweepid-1, Vblksiz, wantz, work);
        else /*if(myid%2 == 1)*/
           CORE_ssbtype3cb(n, nb, A, lda, V, TAU, stind-1, edind-1, sweepid-1, Vblksiz, wantz, work);

        if(blklastind >= (n-1))  break;
    }
}


text:
static char help[] =
"Structured-grid minimal surface equation (MSE) in 2D.  Option prefix mse_.\n"
"Equation is\n"
"  - div ( (1 + |grad u|^2)^q grad u ) = 0\n"
"on the unit square [0,1]x[0,1] subject to Dirichlet boundary\n"
"conditions u = g(x,y).  Power q defaults to -1/2 (i.e. MSE) but is adjustable.\n"
"Catenoid and tent boundary conditions are implemented; catenoid is an exact\n"
"solution.  We re-use the Jacobian from the Poisson equation, but it is suitable\n"
"only for low-amplitude g, or as preconditioning material in -snes_mf_operator.\n"
"Options -snes_fd_color and -snes_grid_sequence K are recommended.\n"
"This code is multigrid (GMG) capable.\n\n";

#include <petsc.h>
#include "../ch6/poissonfunctions.h"
#include "../quadrature.h"

typedef struct {
    double    q,          // the exponent in the diffusivity;
                          // =-1/2 for minimal surface eqn; =0 for Laplace eqn
              tent_H,     // height of tent door along y=0 boundary
              catenoid_c; // parameter in catenoid formula
    int       quaddegree; // quadrature degree used in -mse_monitor
} MinimalCtx;

// Dirichlet boundary conditions
static double g_bdry_tent(double x, double y, double z, void *ctx) {
    PoissonCtx *user = (PoissonCtx*)ctx;
    MinimalCtx *mctx = (MinimalCtx*)(user->addctx);
    if (x < 1.0e-8) {
        return 2.0 * mctx->tent_H * (y < 0.5 ? y : 1.0 - y);
    } else
        return 0;
}

static double g_bdry_catenoid(double x, double y, double z, void *ctx) {
    PoissonCtx   *user = (PoissonCtx*)ctx;
    MinimalCtx   *mctx = (MinimalCtx*)(user->addctx);
    const double c = mctx->catenoid_c;
    return c * cosh(x/c) * sin(acos( (y/c) / cosh(x/c) ));
}

// the coefficient (diffusivity) of minimal surface equation, as a function
//   of  w = |grad u|^2
static double DD(double w, double q) {
    return pow(1.0 + w,q);
}

typedef enum {TENT, CATENOID} ProblemType;
static const char* ProblemTypes[] = {"tent","catenoid",
                                     "ProblemType", "", NULL};

extern PetscErrorCode FormExactFromG(DMDALocalInfo*, Vec, PoissonCtx*);
extern PetscErrorCode FormFunctionLocal(DMDALocalInfo*, double**,
                                        double **FF, PoissonCtx*);
extern PetscErrorCode MSEMonitor(SNES, int, double, void*);

int main(int argc,char **argv) {
    PetscErrorCode ierr;
    DM             da;
    SNES           snes;
    Vec            u_initial, u;
    PoissonCtx     user;
    MinimalCtx     mctx;
    PetscBool      monitor = PETSC_FALSE,
                   exact_init = PETSC_FALSE;
    DMDALocalInfo  info;
    ProblemType    problem = CATENOID;

    // defaults:
    mctx.q = -0.5;
    mctx.tent_H = 1.0;
    mctx.catenoid_c = 1.1;  // case shown in Figure in book
    mctx.quaddegree = 3;
    user.cx = 1.0;
    user.cy = 1.0;
    user.cz = 1.0;

    PetscInitialize(&argc,&argv,NULL,help);
    ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"mse_",
                             "minimal surface equation solver options",""); CHKERRQ(ierr);
    ierr = PetscOptionsReal("-catenoid_c",
                            "parameter for problem catenoid; c >= 1 required",
                            "minimal.c",mctx.catenoid_c,&(mctx.catenoid_c),NULL); CHKERRQ(ierr);
    ierr = PetscOptionsBool("-exact_init",
                            "initial Newton iterate = continuum exact solution; only for catenoid",
                            "minimal.c",exact_init,&(exact_init),NULL);CHKERRQ(ierr);
    ierr = PetscOptionsBool("-monitor",
                            "print surface area and diffusivity bounds at each SNES iteration",
                            "minimal.c",monitor,&(monitor),NULL);CHKERRQ(ierr);
    ierr = PetscOptionsReal("-q",
                            "power of (1+|grad u|^2) in diffusivity",
                            "minimal.c",mctx.q,&(mctx.q),NULL); CHKERRQ(ierr);
    ierr = PetscOptionsInt("-quaddegree",
                            "quadrature degree (=1,2,3) used in -mse_monitor",
                            "minimal.c",mctx.quaddegree,&(mctx.quaddegree),NULL); CHKERRQ(ierr);
    ierr = PetscOptionsEnum("-problem",
                            "problem type determines boundary conditions",
                            "minimal.c",ProblemTypes,(PetscEnum)problem,(PetscEnum*)&problem,
                            NULL); CHKERRQ(ierr);
    ierr = PetscOptionsReal("-tent_H",
                            "'door' height for problem tent",
                            "minimal.c",mctx.tent_H,&(mctx.tent_H),NULL); CHKERRQ(ierr);
    ierr = PetscOptionsEnd(); CHKERRQ(ierr);

    user.addctx = &mctx;   // attach MSE-specific parameters
    switch (problem) {
        case TENT:
            if (exact_init) {
                SETERRQ(PETSC_COMM_WORLD,2,
                    "initialization with exact solution only possible for -mse_problem catenoid\n");
            }
            user.g_bdry = &g_bdry_tent;
            break;
        case CATENOID:
            if (mctx.catenoid_c < 1.0) {
                SETERRQ(PETSC_COMM_WORLD,3,
                    "catenoid exact solution only valid if c >= 1\n");
            }
            if ((exact_init) && (mctx.q != -0.5)) {
                SETERRQ(PETSC_COMM_WORLD,4,
                    "initialization with catenoid exact solution only possible if q=-0.5\n");
            }
            user.g_bdry = &g_bdry_catenoid;
            break;
        default:
            SETERRQ(PETSC_COMM_WORLD,5,"unknown problem type\n");
    }

    ierr = DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,
                        DMDA_STENCIL_BOX,  // contrast with fish2
                        3,3,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da); CHKERRQ(ierr);
    ierr = DMSetApplicationContext(da,&user); CHKERRQ(ierr);
    ierr = DMSetFromOptions(da); CHKERRQ(ierr);
    ierr = DMSetUp(da); CHKERRQ(ierr);  // this must be called BEFORE SetUniformCoordinates
    ierr = DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0); CHKERRQ(ierr);

    ierr = SNESCreate(PETSC_COMM_WORLD,&snes); CHKERRQ(ierr);
    ierr = SNESSetDM(snes,da); CHKERRQ(ierr);
    ierr = DMDASNESSetFunctionLocal(da,INSERT_VALUES,
               (DMDASNESFunction)FormFunctionLocal,&user); CHKERRQ(ierr);
    // this is the Jacobian of the Poisson equation, thus ONLY APPROXIMATE;
    //     generally use -snes_fd_color or -snes_mf_operator
    ierr = DMDASNESSetJacobianLocal(da,
               (DMDASNESJacobian)Poisson2DJacobianLocal,&user); CHKERRQ(ierr);
    if (monitor) {
        ierr = SNESMonitorSet(snes,MSEMonitor,&user,NULL); CHKERRQ(ierr);
    }
    ierr = SNESSetFromOptions(snes); CHKERRQ(ierr);

    ierr = DMGetGlobalVector(da,&u_initial); CHKERRQ(ierr);
    if ((problem == CATENOID) && (mctx.q == -0.5) && (exact_init)) {
        ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);
        ierr = FormExactFromG(&info,u_initial,&user); CHKERRQ(ierr);
    } else {
        // initial iterate has u=g on boundary and u=0 in interior
        ierr = InitialState(da, ZEROS, PETSC_TRUE, u_initial, &user); CHKERRQ(ierr);
    }

//STARTSNESSOLVE
    ierr = SNESSolve(snes,NULL,u_initial); CHKERRQ(ierr);
    ierr = DMRestoreGlobalVector(da,&u_initial); CHKERRQ(ierr);
    ierr = DMDestroy(&da); CHKERRQ(ierr);
    ierr = SNESGetDM(snes,&da); CHKERRQ(ierr);
    ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"done on %d x %d grid and problem %s",
                       info.mx,info.my,ProblemTypes[problem]); CHKERRQ(ierr);
    ierr = SNESGetSolution(snes,&u); CHKERRQ(ierr);
//ENDSNESSOLVE

    // evaluate numerical error in exact solution case
    if ((problem == CATENOID) && (mctx.q == -0.5)) {
        Vec    u_exact;
        double errnorm;
        ierr = VecDuplicate(u,&u_exact); CHKERRQ(ierr);
        ierr = FormExactFromG(&info,u_exact,&user); CHKERRQ(ierr);
        ierr = VecAXPY(u,-1.0,u_exact); CHKERRQ(ierr);    // u <- u + (-1.0) uexact
        ierr = VecDestroy(&u_exact); CHKERRQ(ierr);
        ierr = VecNorm(u,NORM_INFINITY,&errnorm); CHKERRQ(ierr);
        ierr = PetscPrintf(PETSC_COMM_WORLD,
                           ":  error |u-uexact|_inf = %.5e\n",errnorm); CHKERRQ(ierr);
    } else {
        ierr = PetscPrintf(PETSC_COMM_WORLD," ...\n"); CHKERRQ(ierr);
    }

    ierr = SNESDestroy(&snes); CHKERRQ(ierr);
    return PetscFinalize();
}

PetscErrorCode FormExactFromG(DMDALocalInfo *info, Vec uexact,
                         PoissonCtx *user) {
    PetscErrorCode ierr;
    int     i, j;
    double  xymin[2], xymax[2], hx, hy, x, y, **auexact;
    ierr = DMDAGetBoundingBox(info->da,xymin,xymax); CHKERRQ(ierr);
    hx = (xymax[0] - xymin[0]) / (info->mx - 1);
    hy = (xymax[1] - xymin[1]) / (info->my - 1);
    ierr = DMDAVecGetArray(info->da,uexact,&auexact); CHKERRQ(ierr);
    for (j = info->ys; j < info->ys + info->ym; j++) {
        y = j * hy;
        for (i = info->xs; i < info->xs + info->xm; i++) {
            x = i * hx;
            auexact[j][i] = user->g_bdry(x,y,0.0,user);
        }
    }
    ierr = DMDAVecRestoreArray(info->da,uexact,&auexact); CHKERRQ(ierr);
    return 0;
}

PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, double **au,
                                 double **FF, PoissonCtx *user) {
    PetscErrorCode ierr;
    MinimalCtx *mctx = (MinimalCtx*)(user->addctx);
    int        i, j;
    double     xymin[2], xymax[2], hx, hy, hxhy, hyhx, x, y,
               ue, uw, un, us, une, use, unw, usw,
               dux, duy, De, Dw, Dn, Ds;
    ierr = DMDAGetBoundingBox(info->da,xymin,xymax); CHKERRQ(ierr);
    hx = (xymax[0] - xymin[0]) / (info->mx - 1);
    hy = (xymax[1] - xymin[1]) / (info->my - 1);
    hxhy = hx / hy;
    hyhx = hy / hx;
    for (j = info->ys; j < info->ys + info->ym; j++) {
        y = j * hy;
        for (i = info->xs; i < info->xs + info->xm; i++) {
            x = i * hx;
            if (j==0 || i==0 || i==info->mx-1 || j==info->my-1) {
                FF[j][i] = au[j][i] - user->g_bdry(x,y,0.0,user);
            } else {
                // assign neighbor values with either boundary condition or
                //     current u at that point (==> symmetric matrix)
                ue  = (i+1 == info->mx-1) ? user->g_bdry(x+hx,y,0.0,user)
                                          : au[j][i+1];
                uw  = (i-1 == 0)          ? user->g_bdry(x-hx,y,0.0,user)
                                          : au[j][i-1];
                un  = (j+1 == info->my-1) ? user->g_bdry(x,y+hy,0.0,user)
                                          : au[j+1][i];
                us  = (j-1 == 0)          ? user->g_bdry(x,y-hy,0.0,user)
                                          : au[j-1][i];
                if (i+1 == info->mx-1 || j+1 == info->my-1) {
                    une = user->g_bdry(x+hx,y+hy,0.0,user);
                } else {
                    une = au[j+1][i+1];
                }
                if (i-1 == 0 || j+1 == info->my-1) {
                    unw = user->g_bdry(x-hx,y+hy,0.0,user);
                } else {
                    unw = au[j+1][i-1];
                }
                if (i+1 == info->mx-1 || j-1 == 0) {
                    use = user->g_bdry(x+hx,y-hy,0.0,user);
                } else {
                    use = au[j-1][i+1];
                }
                if (i-1 == 0 || j-1 == 0) {
                    usw = user->g_bdry(x-hx,y-hy,0.0,user);
                } else {
                    usw = au[j-1][i-1];
                }
                // gradient  (dux,duy)   at east point  (i+1/2,j):
                dux = (ue - au[j][i]) / hx;
                duy = (un + une - us - use) / (4.0 * hy);
                De = DD(dux * dux + duy * duy, mctx->q);
                // ...                   at west point  (i-1/2,j):
                dux = (au[j][i] - uw) / hx;
                duy = (unw + un - usw - us) / (4.0 * hy);
                Dw = DD(dux * dux + duy * duy, mctx->q);
                // ...                  at north point  (i,j+1/2):
                dux = (ue + une - uw - unw) / (4.0 * hx);
                duy = (un - au[j][i]) / hy;
                Dn = DD(dux * dux + duy * duy, mctx->q);
                // ...                  at south point  (i,j-1/2):
                dux = (ue + use - uw - usw) / (4.0 * hx);
                duy = (au[j][i] - us) / hy;
                Ds = DD(dux * dux + duy * duy, mctx->q);
                // evaluate residual
                FF[j][i] = - hyhx * (De * (ue - au[j][i]) - Dw * (au[j][i] - uw))
                           - hxhy * (Dn * (un - au[j][i]) - Ds * (au[j][i] - us));
            }
        }
    }
    return 0;
}

// compute surface area and bounds on diffusivity using Q^1 elements and
// tensor product gaussian quadrature
PetscErrorCode MSEMonitor(SNES snes, int its, double norm, void *user) {
    PetscErrorCode ierr;
    PoissonCtx     *pctx = (PoissonCtx*)(user);
    MinimalCtx     *mctx = (MinimalCtx*)(pctx->addctx);
    DM             da;
    Vec            u, uloc;
    DMDALocalInfo  info;
    const Quad1D   q = gausslegendre[mctx->quaddegree-1];   // from ../quadrature.h
    double         xymin[2], xymax[2], hx, hy, **au, x_i, y_j, x, y,
                   ux, uy, W, D,
                   Dminloc = PETSC_INFINITY, Dmaxloc = 0.0, Dmin, Dmax,
                   arealoc = 0.0, area;
    int            i, j, r, s;
    MPI_Comm       comm;
    ierr = SNESGetDM(snes, &da); CHKERRQ(ierr);
    ierr = SNESGetSolution(snes, &u); CHKERRQ(ierr);
    ierr = DMGetLocalVector(da, &uloc); CHKERRQ(ierr);
    ierr = DMGlobalToLocalBegin(da, u, INSERT_VALUES, uloc); CHKERRQ(ierr);
    ierr = DMGlobalToLocalEnd(da, u, INSERT_VALUES, uloc); CHKERRQ(ierr);
    ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);
    ierr = DMDAGetBoundingBox(info.da,xymin,xymax); CHKERRQ(ierr);
    hx = (xymax[0] - xymin[0]) / (info.mx - 1);
    hy = (xymax[1] - xymin[1]) / (info.my - 1);
    ierr = DMDAVecGetArrayRead(da,uloc,&au); CHKERRQ(ierr);
//STARTMONITORLOOP
    // loop over rectangular cells in grid
    for (j = info.ys; j < info.ys + info.ym; j++) {
        if (j == 0)
            continue;
        y_j = j * hy;  // NE corner of cell is (x_i,y_j)
        for (i = info.xs; i < info.xs + info.xm; i++) {
            if (i == 0)
                continue;
            x_i = i * hx;
            // loop over quadrature points in cell
            for (r = 0; r < q.n; r++) {
                x = x_i - hx + hx * 0.5 * (q.xi[r] + 1);
                for (s = 0; s < q.n; s++) {
                    y = y_j - hy + hy * 0.5 * (q.xi[s] + 1);
                    // gradient of u(x,y) at a quadrature point
                    ux =   (au[j][i] - au[j][i-1])     * (y - (y_j - hy))
                         + (au[j-1][i] - au[j-1][i-1]) * (y_j - y);
                    ux /= hx * hy;
                    uy =   (au[j][i] - au[j-1][i])     * (x - (x_i - hx))
                         + (au[j][i-1] - au[j-1][i-1]) * (x_i - x);
                    uy /= hx * hy;
                    W = ux * ux + uy * uy;
                    // min and max of diffusivity at quadrature points
                    D = DD(W,mctx->q);
                    Dminloc = PetscMin(Dminloc,D);
                    Dmaxloc = PetscMax(Dmaxloc,D);
                    // apply quadrature in surface area formula
                    arealoc += q.w[r] * q.w[s] * PetscSqrtReal(1.0 + W);
                }
            }
        }
    }
//ENDMONITORLOOP
    ierr = DMDAVecRestoreArrayRead(da,uloc,&au); CHKERRQ(ierr);
    ierr = DMRestoreLocalVector(da, &uloc); CHKERRQ(ierr);

    // do global reductions and report
    ierr = PetscObjectGetComm((PetscObject)da,&comm); CHKERRQ(ierr);
    arealoc *= hx * hy / 4.0;  // from change of variables formula
    ierr = MPI_Allreduce(&arealoc,&area,1,MPI_DOUBLE,MPI_SUM,comm); CHKERRQ(ierr);
    ierr = MPI_Allreduce(&Dminloc,&Dmin,1,MPI_DOUBLE,MPI_MIN,comm); CHKERRQ(ierr);
    ierr = MPI_Allreduce(&Dmaxloc,&Dmax,1,MPI_DOUBLE,MPI_MAX,comm); CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"area = %.8f; %.4f <= D <= %.4f\n",
               area,Dmin,Dmax); CHKERRQ(ierr);
    return 0;
}


text:
/* randist/bigauss.c
 *
 * 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.
 */

#include <config.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

/* The Bivariate Gaussian probability distribution is

   p(x,y) dxdy = (1/(2 pi sigma_x sigma_y sqrt(c)))
    exp(-((x/sigma_x)^2 + (y/sigma_y)^2 - 2 r (x/sigma_x)(y/sigma_y))/2c) dxdy

   where c = 1-r^2
*/

void
gsl_ran_bivariate_gaussian (const gsl_rng * r,
                            double sigma_x, double sigma_y, double rho,
                            double *x, double *y)
{
  double u, v, r2, scale;

  do
    {
      /* choose x,y in uniform square (-1,-1) to (+1,+1) */

      u = -1 + 2 * gsl_rng_uniform (r);
      v = -1 + 2 * gsl_rng_uniform (r);

      /* see if it is in the unit circle */
      r2 = u * u + v * v;
    }
  while (r2 > 1.0 || r2 == 0);

  scale = sqrt (-2.0 * log (r2) / r2);

  *x = sigma_x * u * scale;
  *y = sigma_y * (rho * u + sqrt(1 - rho*rho) * v) * scale;
}

double
gsl_ran_bivariate_gaussian_pdf (const double x, const double y,
                                const double sigma_x, const double sigma_y,
                                const double rho)
{
  double u = x / sigma_x ;
  double v = y / sigma_y ;
  double c = 1 - rho*rho ;
  double p = (1 / (2 * M_PI * sigma_x * sigma_y * sqrt(c)))
    * exp (-(u * u - 2 * rho * u * v + v * v) / (2 * c));
  return p;
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
// Copyright © Matt Jones and Contributors. Licensed under the MIT Licence (MIT). See LICENCE.md in the repository root
// for more information.

#ifndef NOVELRT_UTILITIES_MISC_H
#define NOVELRT_UTILITIES_MISC_H

#include <filesystem>
#include <gsl/span>
#include <type_traits>

#if defined(NDEBUG)
#define unused(x) (void)(x)
#else
#define unused(x) (void)(0)
#endif

namespace NovelRT::Utilities
{
    class Misc
    {
    public:
        static inline const char* CONSOLE_LOG_GENERIC = "NovelRT";
        static inline const char* CONSOLE_LOG_APP = "Application";
        static inline const char* CONSOLE_LOG_DOTNET = ".NET";
        static inline const char* CONSOLE_LOG_GFX = "GFX";
        static inline const char* CONSOLE_LOG_STATS = "Statistics";
        static inline const char* CONSOLE_LOG_AUDIO = "Audio";
        static inline const char* CONSOLE_LOG_INPUT = "Input";
        static inline const char* CONSOLE_LOG_WINDOWING = "WindowManager";

        /**
         * @brief Gets the path to the executable.
         *
         * @return The path to the executable.
         */
        static std::filesystem::path getExecutablePath();

        /**
         * @brief Gets the path to the directory that contains the executable. <br/>
         * For example, `/home/stuff/game/best-game-executable` will return `/home/stuff/game`
         *
         * @return The path to the directory that contains the executable.
         */
        static std::filesystem::path getExecutableDirPath()
        {
            return getExecutablePath().parent_path();
        }

        [[nodiscard]] static std::vector<const char*> GetStringSpanAsCharPtrVector(
            const gsl::span<const std::string>& target) noexcept
        {
            size_t extensionLength = target.size();
            std::vector<const char*> targetPtrs{};
            targetPtrs.reserve(extensionLength);

            for (auto&& extension : target)
            {
                targetPtrs.emplace_back(extension.c_str());
            }

            return targetPtrs;
        }
    };
}

template<class T, class U = std::underlying_type_t<T>> constexpr T operator~(T a)
{
    return static_cast<T>(~static_cast<U>(a));
}

template<class T, class U = std::underlying_type_t<T>> constexpr T operator|(T a, T b)
{
    return static_cast<T>((static_cast<U>(a) | static_cast<U>(b)));
}

template<class T, class U = std::underlying_type_t<T>> constexpr T operator&(T a, T b)
{
    return static_cast<T>((static_cast<U>(a) & static_cast<U>(b)));
}

template<class T, class U = std::underlying_type_t<T>> constexpr T operator^(T a, T b)
{
    return static_cast<T>((static_cast<U>(a) ^ static_cast<U>(b)));
}

template<class T, class U = std::underlying_type_t<T>> constexpr T& operator|=(T& a, T b)
{
    return a = static_cast<T>((static_cast<U>(a) | static_cast<U>(b)));
}

template<class T, class U = std::underlying_type_t<T>> constexpr T& operator&=(T& a, T b)
{
    return a = static_cast<T>((static_cast<U>(a) & static_cast<U>(b)));
}

template<class T, class U = std::underlying_type_t<T>> constexpr T& operator^=(T& a, T b)
{
    return a = static_cast<T>((static_cast<U>(a) ^ static_cast<U>(b)));
}

#endif //! NOVELRT_UTILITIES_MISC_H

max_stars_repo_name: Exadon/NovelRT

max_issues_repo_name: BanalityOfSeeking/NovelRT

max_forks_repo_name: BanalityOfSeeking/NovelRT


text:
#include <stdio.h>
#include <math.h>

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_blas.h>

#include <lapacke.h>
#include "kjg_gsl.h"

void
kjg_gsl_matrix_fprintf (FILE * stream, gsl_matrix * m, const char *template)
{
  size_t i, j;
  for (i = 0; i < m->size1; i++)
    {
      fprintf (stream, template, gsl_matrix_get (m, i, 0));
      for (j = 1; j < m->size2; j++)
	{
	  fprintf (stream, "\t");
	  fprintf (stream, template, gsl_matrix_get (m, i, j));
	}
      fprintf (stream, "\n");
    }
}

void
kjg_gsl_matrix_fscanf (FILE * stream, gsl_matrix * m)
{
  size_t i, j;
  double x;
  for (i = 0; i < m->size1; i++)
    {
      for (j = 0; j < m->size2; j++)
	{
	  fscanf (stream, "%lg", &x);
	  gsl_matrix_set (m, i, j, x);
	}
    }
}

void
kjg_gsl_evec_fprintf (FILE * stream,
		      gsl_vector * eval,
		      gsl_matrix * evec, const char *template)
{
  size_t i, j;
  fprintf (stream, "#");
  fprintf (stream, template, gsl_vector_get (eval, 0));
  for (i = 1; i < eval->size; i++)
    {
      fprintf (stream, "\t");
      fprintf (stream, template, gsl_vector_get (eval, i));
    }
  fprintf (stream, "\n");
  kjg_gsl_matrix_fprintf (stream, evec, template);
}

int
kjg_gsl_evec_fscanf (FILE * stream, gsl_vector * eval, gsl_matrix * evec)
{
  size_t i, j;
  int r;
  double x;

  r = fscanf (stream, "#%lg", &x);
  if (r != 1)
    return (r);
  gsl_vector_set (eval, 0, x);

  for (i = 1; i < eval->size; i++)
    {
      r = fscanf (stream, "%lg", &x);
      if (r != 1)
	return (r);
      gsl_vector_set (eval, i, x);
    }

  for (i = 0; i < evec->size1; i++)
    {
      for (j = 0; j < evec->size2; j++)
	{
	  r = fscanf (stream, "%lg", &x);
	  if (r != 1)
	    return (r);
	  gsl_matrix_set (evec, i, j, x);
	}
    }

  return (0);
}

gsl_rng *
kjg_gsl_rng_init ()
{
  const gsl_rng_type *T;
  gsl_rng *r;
  extern long seed;

  gsl_rng_env_setup ();

  gsl_rng_default_seed = seed;

  T = gsl_rng_default;
  r = gsl_rng_alloc (T);

//  fprintf (stderr, "generator type: %s\n", gsl_rng_name (r));
//  fprintf (stderr, "seed = %lu\n", gsl_rng_default_seed);

  return (r);
}

int
kjg_gsl_matrix_frobenius_normalize (gsl_matrix * m)
{
  double s = kjg_gsl_dlange ('F', m);
  double d = m->size1 * m->size2;
  return (gsl_matrix_scale (m, d / s));
}

double
kjg_gsl_dlange (const char norm, const gsl_matrix * m)
{
  return (LAPACKE_dlange (LAPACK_ROW_MAJOR, norm, m->size1, m->size2, m->data,
			  m->tda));
}

int
kjg_gsl_dgeqrf (gsl_matrix * m, gsl_vector * tau)
{
  return (LAPACKE_dgeqrf (LAPACK_ROW_MAJOR, m->size1, m->size2, m->data,
			  m->tda, tau->data));
}

int
kjg_gsl_dorgqr (gsl_matrix * m, gsl_vector * tau)
{
  return (LAPACKE_dorgqr (LAPACK_ROW_MAJOR, m->size2, m->size2, m->size2,
			  m->data, m->tda, tau->data));
}

void
kjg_gsl_ran_ugaussian_pair (const gsl_rng * r, double x[2])
{
  double r2;

  do
    {
      /* choose x,y in uniform square (-1,-1) to (+1,+1) */
      x[0] = -1 + 2 * gsl_rng_uniform_pos (r);
      x[1] = -1 + 2 * gsl_rng_uniform_pos (r);

      /* see if it is in the unit circle */
      r2 = x[0] * x[0] + x[1] * x[1];
    }
  while (r2 > 1.0 || r2 == 0);

  r2 = sqrt (-2.0 * log (r2) / r2);

  x[0] *= r2;
  x[1] *= r2;
}

void
kjg_gsl_ran_ugaussian_matrix (const gsl_rng * r, gsl_matrix * m)
{
  size_t i, j;
  double *data;
  double x, y, r2;

  for (i = 0; i < m->size1; i++)
    {
      data = gsl_matrix_ptr (m, i, 0);

      for (j = 0; j < m->size2 - 1; j += 2)
	{
	  kjg_gsl_ran_ugaussian_pair (r, data);
	  data += 2;
	}

      if (m->size2 % 2)
	*data = gsl_rng_uniform_pos (r);
    }
}

void
kjg_gsl_matrix_QR (gsl_matrix * m)
{
  gsl_vector *tau = gsl_vector_alloc (m->size2);
  kjg_gsl_dgeqrf (m, tau);
  kjg_gsl_dorgqr (m, tau);
  gsl_vector_free (tau);
}

int
kjg_gsl_SVD (gsl_matrix * M, gsl_matrix * V, gsl_vector * S)
{
  size_t big_enough = M->size1 + V->size2;
  double *superb = malloc (big_enough * sizeof (double));
  double *U;
  int info = LAPACKE_dgesvd (LAPACK_ROW_MAJOR,	// row major
			     'O', 'S', M->size1, M->size2, M->data, M->tda,
			     S->data, U,
			     big_enough, V->data, V->tda, superb);
  free (superb);
  return (info);
}

max_stars_repo_name: nevrome/EIG

max_issues_repo_name: nevrome/EIG

max_forks_repo_name: nevrome/EIG


text:
/**
@file

Copyright John Reid 2006

*/

#ifndef BIOPSY_LOG_PROBS_H_
#define BIOPSY_LOG_PROBS_H_

#ifdef _MSC_VER
# pragma once
#endif //_MSC_VER

#include <gsl/gsl_sf_log.h>
#include <gsl/gsl_sf_exp.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_machine.h>

#include <boost/test/floating_point_comparison.hpp>

namespace biopsy {

inline double safe_log( double x )
{
	return x <= 0.0
		? std::numeric_limits< double >::quiet_NaN()
		: gsl_sf_log( x );
}

template<
	typename ParameterRange
>
void
normalise_discrete_probabilities(
	ParameterRange & parameters )
{
	using namespace boost;

	const double sum = std::accumulate( begin( parameters ), end( parameters ), 0.0 );
	BOOST_FOREACH( typename range_value< ParameterRange >::type & p, parameters )
	{
		p /= sum;
	}
}

template<
	typename log_vector,
	typename out_vector
>
void
probabilities_from_logs_of_choose(
	const log_vector & logs,
	out_vector & out )
{
	out.resize( logs.size() );

	// find the largest - as we can get underflow problems
	const double largest = *( std::max_element( logs.begin(), logs.end() ) );
	for( unsigned i = 0; logs.size() != i; ++i )
	{
		//scale by the largest
		const double v = logs[i] - largest;
		out[i] = (v <= GSL_LOG_DBL_MIN) ? 0.0 : gsl_sf_exp( v );
		BOOST_ASSERT( ! _isnan( out[i] ) );
	}

	normalise_discrete_probabilities( out );
}

template<
	typename ValueRange,
	typename ParameterRange
>
double
dirichlet_log_pdf( const ValueRange & values, const ParameterRange & parameters )
{
	using namespace boost;

	const double param_sum = std::accumulate( begin( parameters ), end( parameters ), 0.0 );

	double result = gsl_sf_lngamma( param_sum );

	typedef boost::tuple< double, double > param_value_tuple;
	BOOST_FOREACH(
		const param_value_tuple & t,
		make_iterator_range(
			make_zip_iterator( make_tuple( begin( parameters ), begin( values ) ) ),
			make_zip_iterator( make_tuple( end( parameters ), end( values ) ) ) ) )
	{
		const double p = t.get<0>();
		const double v = t.get<1>();

		result -= gsl_sf_lngamma( p );
		result += ( p - 1.0 ) * gsl_sf_log( v );
	}

	return result;
}

template<
	typename ParameterRange,
	typename OutputRange
>
void
draw_from_dirichlet( const ParameterRange & alpha, OutputRange & output )
{
	using namespace boost;

	BOOST_ASSERT( size( alpha ) == size( output ) );

	gsl_ran_dirichlet(
		get_gsl_rng(),
		size( alpha ),
		&*begin( alpha ),
		&*begin( output ) );
}

template<
	typename ParameterRange
>
unsigned
draw_from_multi( const ParameterRange & alpha )
{
	using namespace boost;
	using namespace boost::test_tools;

	BOOST_ASSERT(
		close_at_tolerance< double >( percent_tolerance_t< double >( 1.0 ) )(
			1.0,
			std::accumulate( begin( alpha ), end( alpha ), 0.0 ) ) );

	double p = gsl_rng_uniform( get_gsl_rng() );

	unsigned result = 0;
	BOOST_FOREACH( const double a, alpha )
	{
		if( p < a )
		{
			break;
		}
		p -= a;
		++result;
	}

	BOOST_ASSERT( result < size( alpha ) );

	return result;
}

} //namespace biopsy

#endif //BIOPSY_LOG_PROBS_H_

max_stars_repo_name: JohnReid/biopsy

max_issues_repo_name: JohnReid/biopsy

max_forks_repo_name: JohnReid/biopsy


text:
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_math.h>
#include "gsl_cblas.h"

#include "tests.h"

void
test_amax (void) {
  {
   int N = 1;
   float X[] = { -0.388f };
   int incX = -1;
   int expected = 0;
   int k;
   k = cblas_isamax(N, X, incX);
   gsl_test_int(k, expected, "samax(case 52)");
  };

  {
   int N = 1;
   double X[] = { 0.247 };
   int incX = -1;
   int expected = 0;
   int k;
   k = cblas_idamax(N, X, incX);
   gsl_test_int(k, expected, "damax(case 53)");
  };

  {
   int N = 1;
   float X[] = { 0.704f, 0.665f };
   int incX = -1;
   int expected = 0;
   int k;
   k = cblas_icamax(N, X, incX);
   gsl_test_int(k, expected, "camax(case 54)");
  };

  {
   int N = 1;
   double X[] = { -0.599, -0.758 };
   int incX = -1;
   int expected = 0;
   int k;
   k = cblas_izamax(N, X, incX);
   gsl_test_int(k, expected, "zamax(case 55)");
  };

  {
   int N = 2;
   float X[] = { 0.909f, 0.037f };
   int incX = 1;
   int expected = 0;
   int k;
   k = cblas_isamax(N, X, incX);
   gsl_test_int(k, expected, "samax(case 56)");
  };

  {
   int N = 2;
   double X[] = { 0.271, -0.426 };
   int incX = 1;
   int expected = 1;
   int k;
   k = cblas_idamax(N, X, incX);
   gsl_test_int(k, expected, "damax(case 57)");
  };

  {
   int N = 2;
   float X[] = { -0.648f, 0.317f, 0.62f, 0.392f };
   int incX = 1;
   int expected = 1;
   int k;
   k = cblas_icamax(N, X, incX);
   gsl_test_int(k, expected, "camax(case 58)");
  };

  {
   int N = 2;
   double X[] = { -0.789, 0.352, 0.562, 0.697 };
   int incX = 1;
   int expected = 1;
   int k;
   k = cblas_izamax(N, X, incX);
   gsl_test_int(k, expected, "zamax(case 59)");
  };

  {
   int N = 2;
   float X[] = { 0.487f, 0.918f };
   int incX = -1;
   int expected = 0;
   int k;
   k = cblas_isamax(N, X, incX);
   gsl_test_int(k, expected, "samax(case 60)");
  };

  {
   int N = 2;
   double X[] = { 0.537, 0.826 };
   int incX = -1;
   int expected = 0;
   int k;
   k = cblas_idamax(N, X, incX);
   gsl_test_int(k, expected, "damax(case 61)");
  };

  {
   int N = 2;
   float X[] = { 0.993f, 0.172f, -0.825f, 0.873f };
   int incX = -1;
   int expected = 0;
   int k;
   k = cblas_icamax(N, X, incX);
   gsl_test_int(k, expected, "camax(case 62)");
  };

  {
   int N = 2;
   double X[] = { 0.913, -0.436, -0.134, 0.129 };
   int incX = -1;
   int expected = 0;
   int k;
   k = cblas_izamax(N, X, incX);
   gsl_test_int(k, expected, "zamax(case 63)");
  };

}


text:
/*
 * BRAINS
 * (B)LR (R)everberation-mapping (A)nalysis (I)n AGNs with (N)ested (S)ampling
 * Yan-Rong Li, liyanrong@ihep.ac.cn
 * Thu, Aug 4, 2016
 */

/*!
 *  \file reconstruct_sa.c
 *  \brief reconstruct sa and RM data.
 */

#ifdef SA

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <float.h>
#include <gsl/gsl_interp.h>

#include "brains.h"

void *best_model_sa;      /*!< best model */
void *best_model_std_sa;  /*!< standard deviation of the best model */

/*!
 * postprocessing.
 */
void postprocess_sa()
{
  char posterior_sample_file[BRAINS_MAX_STR_LENGTH];
  int num_ps, i, j, k;
  double *pm, *pmstd;
  void *posterior_sample, *post_model;
  int size_of_modeltype = num_params * sizeof(double);

  best_model_sa = malloc(size_of_modeltype);
  best_model_std_sa = malloc(size_of_modeltype);

  if(thistask == roottask)
  {
    // initialize smoothing workspace
    char fname[200];
    FILE *fp, *fline, *fsa;

    sa_smooth_init(n_vel_sa_data, vel_sa_data, parset.sa_InstRes);

    // get number of lines in posterior sample file
    get_posterior_sample_file(dnest_options_file, posterior_sample_file);

    //file for posterior sample
    fp = fopen(posterior_sample_file, "r");
    if(fp == NULL)
    {
      fprintf(stderr, "# Error: Cannot open file %s.\n", posterior_sample_file);
      exit(0);
    }

    //file for line reconstruction
    sprintf(fname, "%s/%s", parset.file_dir, "data/sa_line_rec.txt");
    fline = fopen(fname, "w");
    if(fline == NULL)
    {
      fprintf(stderr, "# Error: Cannot open file %s.\n", fname);
      exit(0);
    }

    //file for phase
    sprintf(fname, "%s/%s", parset.file_dir, "data/sa_phase_rec.txt");
    fsa = fopen(fname, "w");
    if(fsa == NULL)
    {
      fprintf(stderr, "# Error: Cannot open file %s.\n", fname);
      exit(0);
    }

    // read number of lines
    if(fscanf(fp, "# %d", &num_ps) < 1)
    {
      fprintf(stderr, "# Error: Cannot read file %s.\n", posterior_sample_file);
      exit(0);
    }
    printf("# Number of points in posterior sample: %d\n", num_ps);

    post_model = malloc(size_of_modeltype);
    posterior_sample = malloc(num_ps * size_of_modeltype);

    force_update = 1;
    which_parameter_update = -1; // force to update the transfer function
    which_particle_update = 0;

    for(i=0; i<num_ps; i++)
    {
      // read lines
      for(j=0; j<num_params; j++)
      {
        if(fscanf(fp, "%lf", (double *)post_model + j) < 1)
        {
          fprintf(stderr, "# Error: Cannot read file %s.\n", posterior_sample_file);
          exit(0);
        }
      }
      fscanf(fp, "\n");

      //store model
      memcpy(posterior_sample+i*size_of_modeltype, post_model, size_of_modeltype);

      calculate_sa_from_blrmodel(post_model, 0);

      //if( i % (num_ps/10+1) == 0)
      {
        // output sa line
        for(j=0; j<n_vel_sa_data; j++)
        {
          fprintf(fline, "%e %e\n", wave_sa_data[j], Fline_sa[j]);
        }
        fprintf(fline, "\n");

        // output sa phase
        for(k=0; k<n_base_sa_data; k++)
        {
          for(j=0; j<n_vel_sa_data; j++)
          {
            fprintf(fsa, "%e %e\n", wave_sa_data[j], phase_sa[k*n_vel_sa_data + j]/(PhaseFactor * wave_sa_data[j]) );
          }
        }
        fprintf(fsa, "\n");
      }
    }

    fclose(fp);
    fclose(fline);
    fclose(fsa);

    sa_smooth_end();

    pm = (double *)best_model_sa;
    pmstd = (double *)best_model_std_sa;
    for(j=0; j<num_params; j++)
    {
      pm[j] = pmstd[j] = 0.0;
    }
    for(i=0; i<num_ps; i++)
    {
      for(j =0; j<num_params; j++)
        pm[j] += *((double *)posterior_sample + i*num_params + j );
    }

    for(j=0; j<num_params; j++)
      pm[j] /= num_ps;

    for(i=0; i<num_ps; i++)
    {
      for(j=0; j<num_params; j++)
        pmstd[j] += pow( *((double *)posterior_sample + i*num_params + j ) - pm[j], 2.0 );
    }

    for(j=0; j<num_params; j++)
    {
      if(num_ps > 1)
        pmstd[j] = sqrt(pmstd[j]/(num_ps-1.0));
      else
        pmstd[j] = 0.0;
    }

    for(j = 0; j<num_params_sa; j++)
      printf("Best params %d %f +- %f\n", j, *((double *)best_model_sa + j),
                                             *((double *)best_model_std_sa+j) );

    free(post_model);
    free(posterior_sample);
  }
  return;
}

/*!
 * this function run dnest sampleing, reconstruct light curves using the best estimates for parameters.
 */
void reconstruct_sa()
{
  int i, argc=0;
  char **argv;

  //configure restart of dnest
  argv = malloc(9*sizeof(char *));
  for(i=0; i<9; i++)
  {
    argv[i] = malloc(BRAINS_MAX_STR_LENGTH*sizeof(char));
  }
  //setup argc and argv
  strcpy(argv[argc++], "dnest");
  strcpy(argv[argc++], "-s");
  strcpy(argv[argc], parset.file_dir);
  strcat(argv[argc++], "/data/restartsa_dnest.txt");

  if(parset.flag_restart == 1)
  {
    strcpy(argv[argc++], "-r");
    strcpy(argv[argc], parset.file_dir);
    strcat(argv[argc], "/");
    strcat(argv[argc++], "data/restartsa_dnest.txt");
  }
  if(parset.flag_postprc == 1)
  {
    strcpy(argv[argc++], "-p");
  }
  if(parset.flag_temp == 1)
  {
    sprintf(argv[argc++], "-t%f", parset.temperature);
  }
  if(parset.flag_sample_info == 1)
  {
    strcpy(argv[argc++], "-c");
  }

  //level-dependent sampling
  {
    strcpy(argv[argc++], "-l");
  }

  reconstruct_sa_init();

  sa_smooth_init(n_vel_sa_data, vel_sa_data, parset.sa_InstRes);
  dnest_sa(argc, argv);
  sa_smooth_end();

  if(parset.flag_exam_prior != 1 && parset.flag_para_name != 1)
  {
    postprocess_sa();

    if(thistask == roottask)
    {
      FILE *fp;
      char fname[200];
      int j, k;

      force_update = 1;
      which_parameter_update = -1; // force to update the transfer function
      which_particle_update = 0;

      sa_smooth_init(n_vel_sa_data, vel_sa_data, parset.sa_InstRes);
      calculate_sa_from_blrmodel(best_model_sa, 1);
      sa_smooth_end();

      sprintf(fname, "%s/%s", parset.file_dir, "data/psa_line.txt");
      fp = fopen(fname, "w");
      if(fp == NULL)
      {
        fprintf(stderr, "# Error: Cannot open file %s.\n", fname);
        exit(0);
      }
      // output sa line
      for(j=0; j<n_vel_sa_data; j++)
      {
        fprintf(fp, "%e %e\n", wave_sa_data[j], Fline_sa[j]);
      }
      fclose(fp);

      //file for phase
      sprintf(fname, "%s/%s", parset.file_dir, "data/psa_phase.txt");
      fp = fopen(fname, "w");
      if(fp == NULL)
      {
        fprintf(stderr, "# Error: Cannot open file %s.\n", fname);
        exit(0);
      }
      for(k=0; k<n_base_sa_data; k++)
      {
        for(j=0; j<n_vel_sa_data; j++)
        {
          fprintf(fp, "%e %e\n", wave_sa_data[j], phase_sa[k*n_vel_sa_data + j]/(PhaseFactor * wave_sa_data[j]) );
        }
        fprintf(fp, "\n");
      }
      fclose(fp);
    }
  }

  reconstruct_sa_end();

  //clear up argv
  for(i=0; i<9; i++)
  {
    free(argv[i]);
  }
  free(argv);

  return;
}

void reconstruct_sa_init()
{
  sprintf(dnest_options_file, "%s/%s", parset.file_dir, "src/OPTIONSSA");
  if(thistask == roottask)
  {
    get_num_particles(dnest_options_file);
  }
  MPI_Bcast(&parset.num_particles, 1, MPI_INT, roottask, MPI_COMM_WORLD);

  phase_sa = malloc(n_base_sa_data * n_vel_sa_data * sizeof(double));
  Fline_sa = malloc(n_vel_sa_data * sizeof(double));

  /* cloud sample related */
  clouds_weight = malloc(parset.n_cloud_per_task * sizeof(double));
  clouds_alpha = malloc(parset.n_cloud_per_task * sizeof(double));
  clouds_beta = malloc(parset.n_cloud_per_task * sizeof(double));
  clouds_vel = malloc(parset.n_cloud_per_task * parset.n_vel_per_cloud * sizeof(double));

  workspace_phase = malloc( (3*n_vel_sa_data)* sizeof(double));

  if(parset.flag_save_clouds && thistask == roottask)
  {
    if(parset.n_cloud_per_task <= 1000)
      icr_cloud_save = 1;
    else
      icr_cloud_save = parset.n_cloud_per_task/1000;

    char fname[200];
    sprintf(fname, "%s/%s", parset.file_dir, parset.cloud_out_file);
    fcloud_out = fopen(fname, "w");
    if(fcloud_out == NULL)
    {
      fprintf(stderr, "# Error: Cannot open file %s\n", fname);
      exit(-1);
    }
  }
  return;
}

void reconstruct_sa_end()
{
  int i;

  free(phase_sa);
  free(Fline_sa);

  for(i=0; i<num_params; i++)
  {
    free(par_range_model[i]);
    free(par_prior_gaussian[i]);
  }
  free(par_range_model);
  free(par_prior_gaussian);
  free(par_prior_model);

  free(par_fix);
  free(par_fix_val);
  free(best_model_sa);
  free(best_model_std_sa);

  /* clouds sample related */
  free(clouds_weight);
  free(clouds_alpha);
  free(clouds_beta);
  free(clouds_vel);

  free(workspace_phase);

  if(parset.flag_save_clouds && thistask==roottask)
  {
    fclose(fcloud_out);
  }
  if(thistask == roottask)
  {
    printf("Ends reconstruct_sa.\n");
  }
  return;
}

/*!
 * this function calculate probability.
 *
 * At each MCMC step, only one parameter is updated, which only changes some values; thus,
 * optimization that reuses the unchanged values can improve computation efficiency.
 */
double prob_sa(const void *model)
{
  double prob_sa = 0.0, var2, dy;
  int i, j;

  calculate_sa_from_blrmodel(model, 0);

  for(j=0; j<n_epoch_sa_data; j++)
  {
    for(i=0; i<n_vel_sa_data; i++)
    {
      dy = Fline_sa[i] - Fline_sa_data[i + j*n_vel_sa_data];
      var2 = Flerrs_sa_data[i+j*n_vel_sa_data]*Flerrs_sa_data[i+j*n_vel_sa_data];
      prob_sa += -0.5 * (dy*dy)/var2 - 0.5*log(var2 * 2.0*PI);
    }
  }

  for(j=0; j<n_base_sa_data; j++)
  {
    for(i=0; i<n_vel_sa_data; i++)
    {
      dy = phase_sa[j*n_vel_sa_data + i] - phase_sa_data[j*n_vel_sa_data + i];
      var2 = pherrs_sa_data[j*n_vel_sa_data + i] * pherrs_sa_data[j*n_vel_sa_data + i];
      prob_sa += -0.5 * (dy*dy)/var2 - 0.5*log(var2 * 2.0*PI);
    }
  }

  return prob_sa;
}

#endif

max_stars_repo_name: yzxamos/BRAINS

max_issues_repo_name: yzxamos/BRAINS

max_forks_repo_name: yzxamos/BRAINS