text:
// Computational routines for one-dimensional orthogonal polynomial transforms.

#ifndef FASTTRANSFORMS_H
#define FASTTRANSFORMS_H

#include <stdlib.h>
#include <math.h>
#include <immintrin.h>
#include <cblas.h>

#ifdef _OPENMP
    #include <omp.h>
#else
    #define omp_get_thread_num() 0
    #define omp_get_num_threads() 1
#endif

#define M_SQRT_PI    1.772453850905516027   /* sqrt(pi)       */
#define M_1_SQRT_PI  0.564189583547756287   /* 1/sqrt(pi)     */
#define M_SQRT_PI_2  0.886226925452758014   /* sqrt(pi)/2     */

#define MAX(a,b) ((a) > (b) ? a : b)
#define MIN(a,b) ((a) < (b) ? a : b)

#if __SSE2__
    #define VECTOR_SIZE_2 2
    typedef double double2 __attribute__ ((vector_size (VECTOR_SIZE_2*8)));
    #define vall2(x) ((double2) _mm_set1_pd(x))
    #define vload2(v) ((double2) _mm_loadu_pd(v))
    #define vstore2(u, v) (_mm_storeu_pd(u, v))
#endif
#if __AVX__
    #define VECTOR_SIZE_4 4
    typedef double double4 __attribute__ ((vector_size (VECTOR_SIZE_4*8)));
    #define vall4(x) ((double4) _mm256_set1_pd(x))
    #define vload4(v) ((double4) _mm256_loadu_pd(v))
    #define vstore4(u, v) (_mm256_storeu_pd(u, v))
#endif
#if __AVX512F__
    #define VECTOR_SIZE_8 8
    typedef double double8 __attribute__ ((vector_size (VECTOR_SIZE_8*8)));
    #define vall8(x) ((double8) _mm512_set1_pd(x))
    #define vload8(v) ((double8) _mm512_loadu_pd(v))
    #define vstore8(u, v) (_mm512_storeu_pd(u, v))
#endif

static inline double stirlingseries(const double z);
static inline double Aratio(const int n, const double alpha, const double beta);
static inline double Analphabeta(const int n, const double alpha, const double beta);
static inline double lambda(const double x);
static inline double lambda2(const double x, const double l1, const double l2);

double * plan_leg2cheb(const int normleg, const int normcheb, const int n);
double * plan_cheb2leg(const int normcheb, const int normleg, const int n);
double * plan_ultra2ultra(const int normultra1, const int normultra2, const int n, const double lambda1, const double lambda2);
double * plan_jac2jac(const int normjac1, const int normjac2, const int n, const double alpha, const double beta, const double gamma);

typedef struct {
    double * s;
    double * c;
    int n;
} RotationPlan;

void freeRotationPlan(RotationPlan * RP);

RotationPlan * plan_rotsphere(const int n);

void kernel_sph_hi2lo(const RotationPlan * RP, const int m, double * A);
void kernel_sph_lo2hi(const RotationPlan * RP, const int m, double * A);

void kernel_sph_hi2lo_SSE(const RotationPlan * RP, const int m, double * A);
void kernel_sph_lo2hi_SSE(const RotationPlan * RP, const int m, double * A);

void kernel_sph_hi2lo_AVX(const RotationPlan * RP, const int m, double * A);
void kernel_sph_lo2hi_AVX(const RotationPlan * RP, const int m, double * A);

void kernel_sph_hi2lo_AVX512(const RotationPlan * RP, const int m, double * A);
void kernel_sph_lo2hi_AVX512(const RotationPlan * RP, const int m, double * A);

RotationPlan * plan_rottriangle(const int n, const double alpha, const double beta, const double gamma);

void kernel_tri_hi2lo(const RotationPlan * RP, const int m, double * A);
void kernel_tri_lo2hi(const RotationPlan * RP, const int m, double * A);

void kernel_tri_hi2lo_SSE(const RotationPlan * RP, const int m, double * A);
void kernel_tri_lo2hi_SSE(const RotationPlan * RP, const int m, double * A);

void kernel_tri_hi2lo_AVX(const RotationPlan * RP, const int m, double * A);
void kernel_tri_lo2hi_AVX(const RotationPlan * RP, const int m, double * A);

void kernel_tri_hi2lo_AVX512(const RotationPlan * RP, const int m, double * A);
void kernel_tri_lo2hi_AVX512(const RotationPlan * RP, const int m, double * A);

RotationPlan * plan_rotdisk(const int n);

void kernel_disk_hi2lo(const RotationPlan * RP, const int m, double * A);
void kernel_disk_lo2hi(const RotationPlan * RP, const int m, double * A);

void kernel_disk_hi2lo_SSE(const RotationPlan * RP, const int m, double * A);
void kernel_disk_lo2hi_SSE(const RotationPlan * RP, const int m, double * A);

void kernel_disk_hi2lo_AVX(const RotationPlan * RP, const int m, double * A);
void kernel_disk_lo2hi_AVX(const RotationPlan * RP, const int m, double * A);

void kernel_disk_hi2lo_AVX512(const RotationPlan * RP, const int m, double * A);
void kernel_disk_lo2hi_AVX512(const RotationPlan * RP, const int m, double * A);

typedef struct {
    double * s1;
    double * c1;
    double * s2;
    double * c2;
    double * s3;
    double * c3;
    int n;
    int s;
} SpinRotationPlan;

void freeSpinRotationPlan(SpinRotationPlan * SRP);

SpinRotationPlan * plan_rotspinsphere(const int n, const int s);

void kernel_spinsph_hi2lo(const SpinRotationPlan * SRP, const int m, double * A);
void kernel_spinsph_lo2hi(const SpinRotationPlan * SRP, const int m, double * A);

void kernel_spinsph_hi2lo_SSE(const SpinRotationPlan * SRP, const int m, double * A);
void kernel_spinsph_lo2hi_SSE(const SpinRotationPlan * SRP, const int m, double * A);

void kernel_spinsph_hi2lo_AVX(const SpinRotationPlan * SRP, const int m, double * A);
void kernel_spinsph_lo2hi_AVX(const SpinRotationPlan * SRP, const int m, double * A);

void kernel_spinsph_hi2lo_AVX512(const SpinRotationPlan * SRP, const int m, double * A);
void kernel_spinsph_lo2hi_AVX512(const SpinRotationPlan * SRP, const int m, double * A);

static inline void apply_givens(const double S, const double C, double * X, double * Y);
static inline void apply_givens_t(const double S, const double C, double * X, double * Y);

static inline void apply_givens_SSE(const double S, const double C, double * X, double * Y);
static inline void apply_givens_t_SSE(const double S, const double C, double * X, double * Y);

static inline void apply_givens_AVX(const double S, const double C, double * X, double * Y);
static inline void apply_givens_t_AVX(const double S, const double C, double * X, double * Y);

static inline void apply_givens_AVX512(const double S, const double C, double * X, double * Y);
static inline void apply_givens_t_AVX512(const double S, const double C, double * X, double * Y);

void execute_sph_hi2lo(const RotationPlan * RP, double * A, const int M);
void execute_sph_lo2hi(const RotationPlan * RP, double * A, const int M);

void execute_sph_hi2lo_SSE(const RotationPlan * RP, double * A, double * B, const int M);
void execute_sph_lo2hi_SSE(const RotationPlan * RP, double * A, double * B, const int M);

void execute_sph_hi2lo_AVX(const RotationPlan * RP, double * A, double * B, const int M);
void execute_sph_lo2hi_AVX(const RotationPlan * RP, double * A, double * B, const int M);

void execute_sph_hi2lo_AVX512(const RotationPlan * RP, double * A, double * B, const int M);
void execute_sph_lo2hi_AVX512(const RotationPlan * RP, double * A, double * B, const int M);

void execute_tri_hi2lo(const RotationPlan * RP, double * A, const int M);
void execute_tri_lo2hi(const RotationPlan * RP, double * A, const int M);

void execute_tri_hi2lo_SSE(const RotationPlan * RP, double * A, double * B, const int M);
void execute_tri_lo2hi_SSE(const RotationPlan * RP, double * A, double * B, const int M);

void execute_tri_hi2lo_AVX(const RotationPlan * RP, double * A, double * B, const int M);
void execute_tri_lo2hi_AVX(const RotationPlan * RP, double * A, double * B, const int M);

void execute_tri_hi2lo_AVX512(const RotationPlan * RP, double * A, double * B, const int M);
void execute_tri_lo2hi_AVX512(const RotationPlan * RP, double * A, double * B, const int M);

void execute_disk_hi2lo(const RotationPlan * RP, double * A, const int M);
void execute_disk_lo2hi(const RotationPlan * RP, double * A, const int M);

void execute_disk_hi2lo_SSE(const RotationPlan * RP, double * A, double * B, const int M);
void execute_disk_lo2hi_SSE(const RotationPlan * RP, double * A, double * B, const int M);

void execute_disk_hi2lo_AVX(const RotationPlan * RP, double * A, double * B, const int M);
void execute_disk_lo2hi_AVX(const RotationPlan * RP, double * A, double * B, const int M);

void execute_disk_hi2lo_AVX512(const RotationPlan * RP, double * A, double * B, const int M);
void execute_disk_lo2hi_AVX512(const RotationPlan * RP, double * A, double * B, const int M);

void execute_spinsph_hi2lo(const SpinRotationPlan * SRP, double * A, const int M);
void execute_spinsph_lo2hi(const SpinRotationPlan * SRP, double * A, const int M);

void execute_spinsph_hi2lo_SSE(const SpinRotationPlan * SRP, double * A, double * B, const int M);
void execute_spinsph_lo2hi_SSE(const SpinRotationPlan * SRP, double * A, double * B, const int M);

void execute_spinsph_hi2lo_AVX(const SpinRotationPlan * SRP, double * A, double * B, const int M);
void execute_spinsph_lo2hi_AVX(const SpinRotationPlan * SRP, double * A, double * B, const int M);

void execute_spinsph_hi2lo_AVX512(const SpinRotationPlan * SRP, double * A, double * B, const int M);
void execute_spinsph_lo2hi_AVX512(const SpinRotationPlan * SRP, double * A, double * B, const int M);

typedef struct {
    RotationPlan * RP;
    double * B;
    double * P1;
    double * P2;
    double * P1inv;
    double * P2inv;
} SphericalHarmonicPlan;

void freeSphericalHarmonicPlan(SphericalHarmonicPlan * P);

SphericalHarmonicPlan * plan_sph2fourier(const int n);

void execute_sph2fourier(const SphericalHarmonicPlan * P, double * A, const int N, const int M);
void execute_fourier2sph(const SphericalHarmonicPlan * P, double * A, const int N, const int M);

typedef struct {
    RotationPlan * RP;
    double * B;
    double * P1;
    double * P2;
    double * P3;
    double * P4;
    double * P1inv;
    double * P2inv;
    double * P3inv;
    double * P4inv;
    double alpha;
    double beta;
    double gamma;
} TriangularHarmonicPlan;

void freeTriangularHarmonicPlan(TriangularHarmonicPlan * P);

TriangularHarmonicPlan * plan_tri2cheb(const int n, const double alpha, const double beta, const double gamma);

void execute_tri2cheb(const TriangularHarmonicPlan * P, double * A, const int N, const int M);
void execute_cheb2tri(const TriangularHarmonicPlan * P, double * A, const int N, const int M);

static void alternate_sign(double * A, const int N);

static void chebyshev_normalization(double * A, const int N, const int M);
static void chebyshev_normalization_t(double * A, const int N, const int M);

void permute(const double * A, double * B, const int N, const int M, const int L);
void permute_t(double * A, const double * B, const int N, const int M, const int L);

void permute_sph(const double * A, double * B, const int N, const int M, const int L);
void permute_t_sph(double * A, const double * B, const int N, const int M, const int L);

void permute_tri(const double * A, double * B, const int N, const int M, const int L);
void permute_t_tri(double * A, const double * B, const int N, const int M, const int L);

#define permute_disk(A, B, N, M, L) permute_sph(A, B, N, M, L)
#define permute_t_disk(A, B, N, M, L) permute_t_sph(A, B, N, M, L)

#define permute_spinsph(A, B, N, M, L) permute_sph(A, B, N, M, L)
#define permute_t_spinsph(A, B, N, M, L) permute_t_sph(A, B, N, M, L)

void swap(double * A, double * B, const int N);
void warp(double * A, const int N, const int M, const int L);
void warp_t(double * A, const int N, const int M, const int L);

#endif //FASTTRANSFORMS_H

max_stars_repo_name: dawese/FastTransforms

max_issues_repo_name: dawese/FastTransforms

max_forks_repo_name: dawese/FastTransforms


text:
#pragma once

#include <gsl/string_span.h>

using zstr = gsl::zstring<>;
using czstr = gsl::czstring<>;
using zwstr = gsl::wzstring<>;
using czwstr = gsl::cwzstring<>;

ext: h

max_stars_repo_name: lcmftianci/licodeanalysis

max_issues_repo_name: lcmftianci/licodeanalysis

max_forks_repo_name: lcmftianci/licodeanalysis


text:
/*
   NAME:
      bovy_randvec
   PURPOSE:
      returns a (uniform) random vector with a maximum length
   CALLING SEQUENCE:
      bovy_randvec(gsl_vector * eps, int d, double length)
   INPUT:
      d      - dimension of the vector
      length - maximum length of the random vector
   OUTPUT:
      eps    - random vector
   REVISION HISTORY:
      2008-09-21
*/
#include <math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_rng.h>
#include <proj_gauss_mixtures.h>

void bovy_randvec(gsl_vector * eps, int d, double length){
  length /= sqrt((double)d);
  int dd;
  for (dd = 0; dd != d; ++dd)
    gsl_vector_set(eps,dd,(2.*gsl_rng_uniform(randgen)-1.)*length);

  return;
}

ext: c


text:
#ifndef TTT_UTIL_UTILITIES_H
#define TTT_UTIL_UTILITIES_H

#include <SDL_keyboard.h>

#include <gsl/util>

namespace ttt::util
{
	// Used to avoid having to cast std::size_t to gsl::index.
	template<typename Container>
	[[nodiscard]] inline gsl::index getSize(const Container &container) noexcept
	{
		return gsl::narrow_cast<gsl::index>(container.size());
	}
}

#endif

max_stars_repo_name: itsArtem/TicTacToe

max_issues_repo_name: itsArtem/TicTacToe

max_forks_repo_name: itsArtem/TicTacToe


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

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

#include <config.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>

#include "givens.c"

/*
 * This module contains routines related to the Hessenberg-Triangular
 * reduction of two general real matrices
 *
 * See Golub & Van Loan, "Matrix Computations", 3rd ed, sec 7.7.4
 */

/*
gsl_linalg_hesstri_decomp()
  Perform a reduction to generalized upper Hessenberg form.
Given A and B, this function overwrites A with an upper Hessenberg
matrix H = U^T A V and B with an upper triangular matrix R = U^T B V
with U and V orthogonal.

See Golub & Van Loan, "Matrix Computations" (3rd ed) algorithm 7.7.1

Inputs: A    - real square matrix
        B    - real square matrix
        U    - (output) if non-null, U is stored here on output
        V    - (output) if non-null, V is stored here on output
        work - workspace (length n)

Return: success or error
*/

int
gsl_linalg_hesstri_decomp(gsl_matrix * A, gsl_matrix * B, gsl_matrix * U,
                          gsl_matrix * V, gsl_vector * work)
{
  const size_t N = A->size1;

  if ((N != A->size2) || (N != B->size1) || (N != B->size2))
    {
      GSL_ERROR ("Hessenberg-triangular reduction requires square matrices",
                 GSL_ENOTSQR);
    }
  else if (N != work->size)
    {
      GSL_ERROR ("length of workspace must match matrix dimension",
                 GSL_EBADLEN);
    }
  else
    {
      double cs, sn;          /* rotation parameters */
      size_t i, j;            /* looping */
      gsl_vector_view xv, yv; /* temporary views */

      /* B -> Q^T B = R (upper triangular) */
      gsl_linalg_QR_decomp(B, work);

      /* A -> Q^T A */
      gsl_linalg_QR_QTmat(B, work, A);

      /* initialize U and V if desired */

      if (U)
        {
          gsl_linalg_QR_unpack(B, work, U, B);
        }
      else
        {
          /* zero out lower triangle of B */
          for (j = 0; j < N - 1; ++j)
            {
              for (i = j + 1; i < N; ++i)
                gsl_matrix_set(B, i, j, 0.0);
            }
        }

      if (V)
        gsl_matrix_set_identity(V);

      if (N < 3)
        return GSL_SUCCESS; /* nothing more to do */

      /* reduce A and B */
      for (j = 0; j < N - 2; ++j)
        {
          for (i = N - 1; i >= (j + 2); --i)
            {
              /* step 1: rotate rows i - 1, i to kill A(i,j) */

              /*
               * compute G = [ CS SN ] so that G^t [ A(i-1,j) ] = [ * ]
               *             [-SN CS ]             [ A(i, j)  ]   [ 0 ]
               */
              create_givens(gsl_matrix_get(A, i - 1, j),
                            gsl_matrix_get(A, i, j),
                            &cs,
                            &sn);
              /* invert so drot() works correctly (G -> G^t) */
              sn = -sn;

              /* compute G^t A(i-1:i, j:n) */
              xv = gsl_matrix_subrow(A, i - 1, j, N - j);
              yv = gsl_matrix_subrow(A, i, j, N - j);
              gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);

              /* compute G^t B(i-1:i, i-1:n) */
              xv = gsl_matrix_subrow(B, i - 1, i - 1, N - i + 1);
              yv = gsl_matrix_subrow(B, i, i - 1, N - i + 1);
              gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);

              if (U)
                {
                  /* accumulate U: U -> U G */
                  xv = gsl_matrix_column(U, i - 1);
                  yv = gsl_matrix_column(U, i);
                  gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);
                }

              /* step 2: rotate columns i, i - 1 to kill B(i, i - 1) */

              create_givens(-gsl_matrix_get(B, i, i),
                            gsl_matrix_get(B, i, i - 1),
                            &cs,
                            &sn);
              /* invert so drot() works correctly (G -> G^t) */
              sn = -sn;

              /* compute B(1:i, i-1:i) G */
              xv = gsl_matrix_subcolumn(B, i - 1, 0, i + 1);
              yv = gsl_matrix_subcolumn(B, i, 0, i + 1);
              gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);

              /* apply to A(1:n, i-1:i) */
              xv = gsl_matrix_column(A, i - 1);
              yv = gsl_matrix_column(A, i);
              gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);

              if (V)
                {
                  /* accumulate V: V -> V G */
                  xv = gsl_matrix_column(V, i - 1);
                  yv = gsl_matrix_column(V, i);
                  gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);
                }
            }
        }

      return GSL_SUCCESS;
    }
} /* gsl_linalg_hesstri_decomp() */

max_stars_repo_name: ruslankuzmin/julia


text:
/* rng/test.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 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 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>

void rng_test (const gsl_rng_type * T, unsigned long int seed, unsigned int n,
               unsigned long int result);
void rng_float_test (const gsl_rng_type * T);
void generic_rng_test (const gsl_rng_type * T);
void rng_state_test (const gsl_rng_type * T);
void rng_parallel_state_test (const gsl_rng_type * T);
void rng_read_write_test (const gsl_rng_type * T);
int rng_max_test (gsl_rng * r, unsigned long int *kmax, unsigned long int ran_max) ;
int rng_min_test (gsl_rng * r, unsigned long int *kmin, unsigned long int ran_min, unsigned long int ran_max) ;
int rng_sum_test (gsl_rng * r, double *sigma);
int rng_bin_test (gsl_rng * r, double *sigma);

#define N  10000
#define N2 200000

int
main (void)
{
  const gsl_rng_type ** rngs = gsl_rng_types_setup();  /* get all rng types */

  const gsl_rng_type ** r ;

  gsl_ieee_env_setup ();

  gsl_rng_env_setup ();

  /* specific tests of known results for 10000 iterations with seed = 1 */

  rng_test (gsl_rng_rand, 1, 10000, 1910041713);
  rng_test (gsl_rng_randu, 1, 10000, 1623524161);
  rng_test (gsl_rng_cmrg, 1, 10000, 719452880);
  rng_test (gsl_rng_minstd, 1, 10000, 1043618065);
  rng_test (gsl_rng_mrg, 1, 10000, 2064828650);
  rng_test (gsl_rng_taus, 1, 10000, 2733957125UL);
  rng_test (gsl_rng_taus2, 1, 10000, 2733957125UL);
  rng_test (gsl_rng_taus113, 1, 1000, 1925420673UL);
  rng_test (gsl_rng_transputer, 1, 10000, 1244127297UL);
  rng_test (gsl_rng_vax, 1, 10000, 3051034865UL);

  /* Borosh13 test value from PARI: (1812433253^10000)%(2^32) */
  rng_test (gsl_rng_borosh13, 1, 10000, 2513433025UL);

  /* Fishman18 test value from PARI: (62089911^10000)%(2^31-1) */
  rng_test (gsl_rng_fishman18, 1, 10000, 330402013UL);

  /* Fishman2x test value from PARI:
     ((48271^10000)%(2^31-1) - (40692^10000)%(2^31-249))%(2^31-1) */
  rng_test (gsl_rng_fishman2x, 1, 10000, 540133597UL);

  /* Knuthran2 test value from PARI:
     { xn1=1; xn2=1; for (n=1,10000,
            xn = (271828183*xn1 - 314159269*xn2)%(2^31-1);
            xn2=xn1; xn1=xn; print(xn); ) } */
  rng_test (gsl_rng_knuthran2, 1, 10000, 1084477620UL);

  /* Knuthran test value taken from p188 in Knuth Vol 2. 3rd Ed */
  rng_test (gsl_rng_knuthran, 310952, 1009 * 2009 + 1, 461390032);

  /* Knuthran improved test value from Knuth's source */
  rng_test (gsl_rng_knuthran2002, 310952, 1, 708622036);
  rng_test (gsl_rng_knuthran2002, 310952, 2, 1005450560);
  rng_test (gsl_rng_knuthran2002, 310952, 100 * 2009 + 1, 995235265);
  rng_test (gsl_rng_knuthran2002, 310952, 1009 * 2009 + 1, 704987132);

  /* Lecuyer21 test value from PARI: (40692^10000)%(2^31-249) */
  rng_test (gsl_rng_lecuyer21, 1, 10000, 2006618587UL);

  /* Waterman14 test value from PARI: (1566083941^10000)%(2^32) */
  rng_test (gsl_rng_waterman14, 1, 10000, 3776680385UL);

  /* specific tests of known results for 10000 iterations with seed = 6 */

  /* Coveyou test value from PARI:
     x=6; for(n=1,10000,x=(x*(x+1))%(2^32);print(x);) */

  rng_test (gsl_rng_coveyou, 6, 10000, 1416754246UL);

  /* Fishman20 test value from PARI: (6*48271^10000)%(2^31-1) */
  rng_test (gsl_rng_fishman20, 6, 10000, 248127575UL);

  /* FIXME: the ranlux tests below were made by running the fortran code and
     getting the expected value from that. An analytic calculation
     would be preferable. */

  rng_test (gsl_rng_ranlux, 314159265, 10000, 12077992);
  rng_test (gsl_rng_ranlux389, 314159265, 10000, 165942);

  rng_test (gsl_rng_ranlxs0, 1, 10000, 11904320);
  /* 0.709552764892578125 * ldexp(1.0,24) */

  rng_test (gsl_rng_ranlxs1, 1, 10000, 8734328);
  /* 0.520606517791748047 * ldexp(1.0,24) */

  rng_test (gsl_rng_ranlxs2, 1, 10000, 6843140);
  /* 0.407882928848266602 * ldexp(1.0,24) */

  rng_test (gsl_rng_ranlxd1, 1, 10000, 1998227290UL);
  /* 0.465248546261094020 * ldexp(1.0,32) */

  rng_test (gsl_rng_ranlxd2, 1, 10000, 3949287736UL);
  /* 0.919515205581550532 * ldexp(1.0,32) */

  /* FIXME: the tests below were made by running the original code in
     the ../random directory and getting the expected value from
     that. An analytic calculation would be preferable. */

  rng_test (gsl_rng_slatec, 1, 10000, 45776);
  rng_test (gsl_rng_uni, 1, 10000, 9214);
  rng_test (gsl_rng_uni32, 1, 10000, 1155229825);
  rng_test (gsl_rng_zuf, 1, 10000, 3970);

  /* The tests below were made by running the original code and
     getting the expected value from that. An analytic calculation
     would be preferable. */

  rng_test (gsl_rng_r250, 1, 10000, 1100653588);
  rng_test (gsl_rng_mt19937, 4357, 1000, 1186927261);
  rng_test (gsl_rng_mt19937_1999, 4357, 1000, 1030650439);
  rng_test (gsl_rng_mt19937_1998, 4357, 1000, 1309179303);
  rng_test (gsl_rng_tt800, 0, 10000, 2856609219UL);

  rng_test (gsl_rng_ran0, 0, 10000, 1115320064);
  rng_test (gsl_rng_ran1, 0, 10000, 1491066076);
  rng_test (gsl_rng_ran2, 0, 10000, 1701364455);
  rng_test (gsl_rng_ran3, 0, 10000, 186340785);

  rng_test (gsl_rng_ranmar, 1, 10000, 14428370);

  rng_test (gsl_rng_rand48, 0, 10000, 0xDE095043UL);
  rng_test (gsl_rng_rand48, 1, 10000, 0xEDA54977UL);

  rng_test (gsl_rng_random_glibc2, 0, 10000, 1908609430);
  rng_test (gsl_rng_random8_glibc2, 0, 10000, 1910041713);
  rng_test (gsl_rng_random32_glibc2, 0, 10000, 1587395585);
  rng_test (gsl_rng_random64_glibc2, 0, 10000, 52848624);
  rng_test (gsl_rng_random128_glibc2, 0, 10000, 1908609430);
  rng_test (gsl_rng_random256_glibc2, 0, 10000, 179943260);

  rng_test (gsl_rng_random_bsd, 0, 10000, 1457025928);
  rng_test (gsl_rng_random8_bsd, 0, 10000, 1910041713);
  rng_test (gsl_rng_random32_bsd, 0, 10000, 1663114331);
  rng_test (gsl_rng_random64_bsd, 0, 10000, 864469165);
  rng_test (gsl_rng_random128_bsd, 0, 10000, 1457025928);
  rng_test (gsl_rng_random256_bsd, 0, 10000, 1216357476);

  rng_test (gsl_rng_random_libc5, 0, 10000, 428084942);
  rng_test (gsl_rng_random8_libc5, 0, 10000, 1910041713);
  rng_test (gsl_rng_random32_libc5, 0, 10000, 1967452027);
  rng_test (gsl_rng_random64_libc5, 0, 10000, 2106639801);
  rng_test (gsl_rng_random128_libc5, 0, 10000, 428084942);
  rng_test (gsl_rng_random256_libc5, 0, 10000, 116367984);

  rng_test (gsl_rng_ranf, 0, 10000, 2152890433UL);
  rng_test (gsl_rng_ranf, 2, 10000, 339327233);

  /* Test constant relationship between int and double functions */

  for (r = rngs ; *r != 0; r++)
    rng_float_test (*r);

  /* Test save/restore functions */

  for (r = rngs ; *r != 0; r++)
    rng_state_test (*r);

  for (r = rngs ; *r != 0; r++)
    rng_parallel_state_test (*r);

  for (r = rngs ; *r != 0; r++)
    rng_read_write_test (*r);

  /* generic statistical tests (these are just to make sure that we
     don't get any crazy results back from the generator, i.e. they
     aren't a test of the algorithm, just the implementation) */

  for (r = rngs ; *r != 0; r++)
    generic_rng_test (*r);

  exit (gsl_test_summary ());
}

void
rng_test (const gsl_rng_type * T, unsigned long int seed, unsigned int n,
          unsigned long int result)
{
  gsl_rng *r = gsl_rng_alloc (T);
  unsigned int i;
  unsigned long int k = 0;
  int status;

  if (seed != 0)
    {
      gsl_rng_set (r, seed);
    }

  for (i = 0; i < n; i++)
    {
      k = gsl_rng_get (r);
    }

  status = (k != result);
  gsl_test (status, "%s, %u steps (%u observed vs %u expected)",
            gsl_rng_name (r), n, k, result);

  gsl_rng_free (r);
}

void
rng_float_test (const gsl_rng_type * T)
{
  gsl_rng *ri = gsl_rng_alloc (T);
  gsl_rng *rf = gsl_rng_alloc (T);

  double u, c ;
  unsigned int i;
  unsigned long int k = 0;
  int status = 0 ;

  do
    {
      k = gsl_rng_get (ri);
      u = gsl_rng_get (rf);
    }
  while (k == 0) ;

  c = k / u ;
  for (i = 0; i < N2; i++)
    {
      k = gsl_rng_get (ri);
      u = gsl_rng_get (rf);
      if (c*k != u)
        {
          status = 1 ;
          break ;
        }
    }

  gsl_test (status, "%s, ratio of int to double (%g observed vs %g expected)",
            gsl_rng_name (ri), c, k/u);

  gsl_rng_free (ri);
  gsl_rng_free (rf);
}

void
rng_state_test (const gsl_rng_type * T)
{
  unsigned long int test_a[N], test_b[N];

  int i;

  gsl_rng *r = gsl_rng_alloc (T);
  gsl_rng *r_save = gsl_rng_alloc (T);

  for (i = 0; i < N; ++i)
    {
      gsl_rng_get (r);  /* throw away N iterations */
    }

  gsl_rng_memcpy (r_save, r);   /* save the intermediate state */

  for (i = 0; i < N; ++i)
    {
      test_a[i] = gsl_rng_get (r);
    }

  gsl_rng_memcpy (r, r_save);   /* restore the intermediate state */
  gsl_rng_free (r_save);

  for (i = 0; i < N; ++i)
    {
      test_b[i] = gsl_rng_get (r);
    }

  {
    int status = 0;
    for (i = 0; i < N; ++i)
      {
        status |= (test_b[i] != test_a[i]);
      }
    gsl_test (status, "%s, random number state consistency",
              gsl_rng_name (r));
  }

  gsl_rng_free (r);
}

void
rng_parallel_state_test (const gsl_rng_type * T)
{
  unsigned long int test_a[N], test_b[N];
  unsigned long int test_c[N], test_d[N];
  double test_e[N], test_f[N];

  int i;

  gsl_rng *r1 = gsl_rng_alloc (T);
  gsl_rng *r2 = gsl_rng_alloc (T);

  for (i = 0; i < N; ++i)
    {
      gsl_rng_get (r1);         /* throw away N iterations */
    }

  gsl_rng_memcpy (r2, r1);              /* save the intermediate state */

  for (i = 0; i < N; ++i)
    {
      /* check that there is no hidden state intermixed between r1 and r2 */
      test_a[i] = gsl_rng_get (r1);
      test_b[i] = gsl_rng_get (r2);
      test_c[i] = gsl_rng_uniform_int (r1, 1234);
      test_d[i] = gsl_rng_uniform_int (r2, 1234);
      test_e[i] = gsl_rng_uniform (r1);
      test_f[i] = gsl_rng_uniform (r2);
    }

  {
    int status = 0;
    for (i = 0; i < N; ++i)
      {
        status |= (test_b[i] != test_a[i]);
        status |= (test_c[i] != test_d[i]);
        status |= (test_e[i] != test_f[i]);
      }
    gsl_test (status, "%s, parallel random number state consistency",
              gsl_rng_name (r1));
  }

  gsl_rng_free (r1);
  gsl_rng_free (r2);

}

void
rng_read_write_test (const gsl_rng_type * T)
{
  unsigned long int test_a[N], test_b[N];

  int i;

  gsl_rng *r = gsl_rng_alloc (T);

  for (i = 0; i < N; ++i)
    {
      gsl_rng_get (r);  /* throw away N iterations */
    }

  { /* save the state to a binary file */
    FILE *f = fopen("test.dat", "wb");
    gsl_rng_fwrite(f, r);
    fclose(f);
  }

  for (i = 0; i < N; ++i)
    {
      test_a[i] = gsl_rng_get (r);
    }

  { /* read the state from a binary file */
    FILE *f = fopen("test.dat", "rb");
    gsl_rng_fread(f, r);
    fclose(f);
  }

  for (i = 0; i < N; ++i)
    {
      test_b[i] = gsl_rng_get (r);
    }

  {
    int status = 0;
    for (i = 0; i < N; ++i)
      {
        status |= (test_b[i] != test_a[i]);
      }
    gsl_test (status, "%s, random number generator read and write",
              gsl_rng_name (r));
  }

  gsl_rng_free (r);
}

void
generic_rng_test (const gsl_rng_type * T)
{
  gsl_rng *r = gsl_rng_alloc (T);
  const char *name = gsl_rng_name (r);
  unsigned long int kmax = 0, kmin = 1000;
  double sigma = 0;
  const unsigned long int ran_max = gsl_rng_max (r);
  const unsigned long int ran_min = gsl_rng_min (r);

  int status = rng_max_test (r, &kmax, ran_max);

  gsl_test (status,
            "%s, observed vs theoretical maximum (%lu vs %lu)",
            name, kmax, ran_max);

  status = rng_min_test (r, &kmin, ran_min, ran_max);

  gsl_test (status,
            "%s, observed vs theoretical minimum (%lu vs %lu)",
            name, kmin, ran_min);

  status = rng_sum_test (r, &sigma);

  gsl_test (status,
            "%s, sum test within acceptable sigma (observed %.2g sigma)",
            name, sigma);

  status = rng_bin_test (r, &sigma);

  gsl_test (status,
            "%s, bin test within acceptable chisq (observed %.2g sigma)",
            name, sigma);

  gsl_rng_set (r, 1);   /* set seed to 1 */
  status = rng_max_test (r, &kmax, ran_max);

  gsl_rng_set (r, 1);   /* set seed to 1 */
  status |= rng_min_test (r, &kmin, ran_min, ran_max);

  gsl_rng_set (r, 1);   /* set seed to 1 */
  status |= rng_sum_test (r, &sigma);

  gsl_test (status, "%s, maximum and sum tests for seed=1", name);

  gsl_rng_set (r, 12345);       /* set seed to a "typical" value */
  status = rng_max_test (r, &kmax, ran_max);

  gsl_rng_set (r, 12345);       /* set seed to a "typical" value */
  status |= rng_min_test (r, &kmin, ran_min, ran_max);

  gsl_rng_set (r, 12345);       /* set seed to a "typical" value */
  status |= rng_sum_test (r, &sigma);

  gsl_test (status, "%s, maximum and sum tests for non-default seeds", name);

  gsl_rng_free (r);
}

int
rng_max_test (gsl_rng * r, unsigned long int *kmax, unsigned long int ran_max)
{
  unsigned long int actual_uncovered;
  double expect_uncovered;
  int status;
  unsigned long int max = 0;
  int i;

  for (i = 0; i < N2; ++i)
    {
      unsigned long int k = gsl_rng_get (r);
      if (k > max)
        max = k;
    }

  *kmax = max;

  actual_uncovered = ran_max - max;
  expect_uncovered = (double) ran_max / (double) N2;

  status = (max > ran_max) || (actual_uncovered > 7 * expect_uncovered) ;

  return status;
}

int
rng_min_test (gsl_rng * r, unsigned long int *kmin,
              unsigned long int ran_min, unsigned long int ran_max)
{
  unsigned long int actual_uncovered;
  double expect_uncovered;
  int status;
  unsigned long int min = 1000000000UL;
  int i;

  for (i = 0; i < N2; ++i)
    {
      unsigned long int k = gsl_rng_get (r);
      if (k < min)
        min = k;
    }

  *kmin = min;

  actual_uncovered = min - ran_min;
  expect_uncovered = (double) ran_max / (double) N2;

  status = (min < ran_min) || (actual_uncovered > 7 * expect_uncovered);

  return status;
}

int
rng_sum_test (gsl_rng * r, double *sigma)
{
  double sum = 0;
  int i, status;

  for (i = 0; i < N2; ++i)
    {
      double x = gsl_rng_uniform (r) - 0.5;
      sum += x;
    }

  sum /= N2;

  /* expect the average to have a variance of 1/(12 n) */

  *sigma = sum * sqrt (12.0 * N2);

  /* more than 3 sigma is an error (increased to 3.1 to avoid false alarms) */

  status = (fabs (*sigma) > 3.1 || fabs(*sigma) < 0.003);

  if (status) {
      fprintf(stderr,"sum=%g, sigma=%g\n",sum,*sigma);
  }

  return status;
}

#define BINS 17
#define EXTRA 10
int
rng_bin_test (gsl_rng * r, double *sigma)
{
  int count[BINS+EXTRA];
  double chisq = 0;
  int i, status;

  for (i = 0; i < BINS+EXTRA; i++)
      count[i] = 0 ;

  for (i = 0; i < N2; i++)
    {
      int j = gsl_rng_uniform_int (r, BINS);
      count[j]++ ;
    }

  chisq = 0 ;
  for (i = 0; i < BINS; i++)
    {
      double x = (double)N2/(double)BINS ;
      double d = (count[i] - x) ;
      chisq += (d*d) / x;
    }

  *sigma = sqrt(chisq/BINS) ;

  /* more than 3 sigma is an error */

  status = (fabs (*sigma) > 3 || fabs(*sigma) < 0.003);

  for (i = BINS; i < BINS+EXTRA; i++)
    {
      if (count[i] != 0)
        {
          status = 1 ;
          gsl_test (status,
                    "%s, wrote outside range in bin test "
                    "(%d observed vs %d expected)",
                    gsl_rng_name(r), i, BINS - 1);
        }
    }

  return status;
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
// This is a plugin program that serves mainly as an example of how
// analysis programs for ssv should be written, and gives people
// something to cut and paste from.  It computes and prints histogram
// information about the tiles it analyzes.

#include <stdlib.h>

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

#include <float_image.h>

// As specified in the analysis interface provided by ssv, the first
// argument (after argv[0], the invocation name) is the number of
// tiles being passed into the analysis plugin.  The second and
// succeeding arguments are also as as described in the output
// produced by 'ssv --help'.
#define TILE_COUNT_ARGUMENT_POSITION 1

// This type describes the information that we get passed for each
// tile.
typedef struct {
  GString *base_name;
  size_t width, height;
  GString *file_name;
  FloatImage *tile;
} tile_spec_type;

int
main (int argc, char **argv)
{
  char *prog_name = argv[0];
  prog_name = prog_name;   // You might use this for reporting errors.

  // How many image tiles are we getting passed?
  int tile_count = atoi (argv[TILE_COUNT_ARGUMENT_POSITION]);

  // Print a string showing how we were invoked.
  GString *invocation = g_string_new ("");
  int ii;
  for ( ii = 0 ; ii < argc ; ii++ ) {
    g_string_append_printf (invocation, "%s ", argv[ii]);
  }
  g_print ("\nRunning analysis program '%s'...\n", invocation->str);
  g_string_free (invocation, TRUE);

  // Read the argument into an array of tile_spec_type structures.
  GPtrArray *tile_specs = g_ptr_array_new ();
  int current_arg = TILE_COUNT_ARGUMENT_POSITION + 1;
  for ( ii = 0 ; ii < tile_count ; ii++ ) {
    tile_spec_type *cs = g_new (tile_spec_type, 1); // Current spec.
    cs->base_name = g_string_new (argv[current_arg++]);
    cs->width = atoi (argv[current_arg++]);
    cs->height = atoi (argv[current_arg++]);
    cs->file_name = g_string_new (argv[current_arg++]);
    cs->tile = float_image_new_from_file (cs->width, cs->height,
					  cs->file_name->str, 0,
					  FLOAT_IMAGE_BYTE_ORDER_BIG_ENDIAN);
    g_ptr_array_add (tile_specs, cs);
  }

  g_print ("\nGot %d total structures to process.\n", tile_count);

  // Compute and print histograms and other information for each tile.
  for ( ii = 0 ; ii < tile_specs->len ; ii++ ) {
    // Current tile spec.
    tile_spec_type *cs = g_ptr_array_index (tile_specs, ii);
    g_print ("\nTile from image %s\n", cs->base_name->str);
    g_print ("Width: %llu\n", (long long unsigned int) cs->width);
    g_print ("Height: %llu\n", (long long unsigned int) cs->height);
    if ( cs->width == 0 || cs->height == 0 ) {
      g_print ("Tile of zero area, nothing to analyze.\n");
    }
    else {
      g_print ("Histogram:\n");
      const size_t bin_count = 20;
      gsl_histogram *hg = gsl_histogram_alloc (bin_count);
      // The test data I use is from an originally byte-valued image; if
      // you have a different sort of image a histogram covering only
      // this data region may be pretty dull.
      gsl_histogram_set_ranges_uniform (hg, 0.0, 255.0);
      guint jj, kk;
      for ( jj = 0 ; jj < cs->width ; jj++ ) {
	for ( kk = 0 ; kk < cs->height ; kk++ ) {
	  float pv = float_image_get_pixel (cs->tile, jj, kk);
	  gsl_histogram_increment (hg, pv);
	}
      }
      g_print ("  Range Start     Range End    Occurences\n");
      g_print ("-----------------------------------------\n");
      gsl_histogram_fprintf (stdout, hg, "%13.2f", "%13.0f");
      gsl_histogram_free (hg);
    }
    g_print ("\n");
  }

  // Free the tile specifications array.
  for ( ii = 0 ; ii < tile_specs->len ; ii++ ) {
    tile_spec_type *cs = g_ptr_array_index (tile_specs, ii);
    float_image_free (cs->tile);
    g_string_free (cs->file_name, TRUE);
    g_string_free (cs->base_name, TRUE);
    g_free (cs);
  }
  g_ptr_array_free (tile_specs, TRUE);

  exit (EXIT_SUCCESS);
}

max_stars_repo_name: glshort/MapReady

max_issues_repo_name: glshort/MapReady

max_forks_repo_name: glshort/MapReady


text:
#ifndef METHODS_HELLER_H
#define METHODS_HELLER_H

#include <gsl/gsl_multimin.h>

namespace method {
// use your method name to create a subspace for your
// implementation of details
namespace heller_cwa {
namespace details {

struct HellerParam {

  // exp(a[0] p^2 + a[1] p + a[2] p q + a[3] q + a[4] q^2 + a[5])
  // alpha = 0.5 / a[0]
  // delta = a[1]
  // epsilon = a[2]

  arma::cx_vec a;
  arma::mat points;
  math::Polynomial<cx_double> V_eff_0;
  double mass;
};

inline
cx_double heller_gaussian(const arma::cx_vec & a,
                          const double q,
                          const double p) {
  return std::exp(a(0) * std::pow(p, 2)
                  + a(1) * p
                  + a(2) * p * q
                  + a(3) * q
                  + a(4) * std::pow(q, 2)
                  + a(5));
}

inline
arma::cx_vec heller_gaussian(const arma::cx_vec & a,
                             const arma::mat & points) {
  if (points.n_rows != 2) {
    throw Error("heller_gaussian only supports one-dimensional phase space");
  }

  arma::cx_vec result(points.n_cols);
#pragma omp parallel for
  for (arma::uword i = 0; i < result.n_elem; i++) {
    result(i) = heller_gaussian(a, points(0, i), points(1, i));
  }

  return result;
}

inline
math::Gaussian<cx_double> heller_gaussian(const arma::cx_vec & a) {

  const arma::cx_mat binomial{{{-a(4) * 2.0, -a(2)}, {-a(2), -a(0) * 2.0}}};
  const arma::cx_mat covariance = arma::inv(binomial);
  const arma::cx_vec center = covariance * arma::cx_vec{a(3),a(1)};

  return math::Gaussian<cx_double>(covariance, center,
      std::exp(a(5)) * std::exp(0.5 * arma::dot(center, arma::cx_vec{a(3),a(1)})));

}

template<typename Potential>
arma::mat force(const Potential & potential,
                const arma::mat & positions) {
  arma::mat result = arma::mat(arma::size(positions));

#pragma omp parallel for
  for (arma::uword i = 0; i < positions.n_cols; i++) {
    const arma::vec position = positions.col(i);
    for (arma::uword j = 0; j < positions.n_rows; j++) {
      result(j, i) = -potential.derivative(j).at(position);
    }
  }

  return result;
}

template<typename Function>
inline
auto effective_potential(const Function & potential,
                         const cx_double alpha,
                         const cx_double delta,
                         const cx_double epsilon,
                         const arma::uword l,
                         const arma::uword cut_off) {

  if (potential.dim() > 1) {
    throw Error("Currently effective potential is only valid "
                "for one dimensional potential");
  }

  const auto functor = [alpha](const Function & function) -> Function {
    return function.derivative(0).derivative(0) / 8.0 / alpha;
  };

  const auto a = exp(functor, derivative(potential, arma::uvec{l}), cut_off) *
                 std::pow(-1, l + 1);
  const auto b = exp(functor, derivative(potential, arma::uvec{l}), cut_off);

  const std::vector<math::Polynomial<cx_double>>
      a_inner =
      {math::Polynomial<cx_double>(
          arma::cx_vec{1.0 + epsilon / 2.0,
                       cx_double{0.0, 1.0} / 2.0 / alpha,
                       0.5 * delta},
          lmat{{1, 0, 0},
               {0, 1, 0}})};

  const std::vector<math::Polynomial<cx_double>>
      b_inner =
      {math::Polynomial<cx_double>(
          arma::cx_vec{1.0 - epsilon / 2.0,
                       -cx_double{0.0, 1.0} / 2.0 / alpha,
                       -0.5 * delta},
          lmat{{1, 0, 0},
               {0, 1, 0}})};

  const auto combined = (a(a_inner) + b(b_inner)) * 0.5;

  return combined;
}

inline
math::Polynomial<cx_double> E_function(const arma::cx_vec & a,
                                       const arma::cx_vec & a_derivatives,
                                       const math::Polynomial<cx_double> & V_eff_0,
                                       const double mass) {

  if (a.n_elem != 6 || a_derivatives.n_elem != 6) {
    throw Error(
        "The number of elements for a's or their derivatives is invalid for E function");
  }

  const math::Polynomial<cx_double> without_derivatives(
      arma::cx_vec{-a(2) / mass, -2.0 * a(4) / mass, -a(3) / mass},
      lmat{{{0, 1, 0}, {2, 1, 1}}}
  );

  const math::Polynomial<cx_double> derivatives(
      arma::cx_vec{-a_derivatives(0),
                   -a_derivatives(1),
                   -a_derivatives(2),
                   -a_derivatives(3),
                   -a_derivatives(4),
                   -a_derivatives(5)},
      lmat{{{0, 0, 1, 1, 2, 0}, {2, 1, 1, 0, 0, 0}}}
  );

  return without_derivatives + derivatives - V_eff_0 * cx_double{0.0, 2.0};

}

inline
double I_function(
    const arma::cx_vec & a,
    const arma::cx_vec & a_derivatives,
    const arma::mat & points,
    const math::Polynomial<cx_double> & V_eff_0,
    const double mass) {

  return arma::sum(arma::pow(arma::abs(
      at(E_function(a, a_derivatives, V_eff_0, mass), points)
      % heller_gaussian(a, points)), 2));

}

inline
arma::vec I_function_derivative(
    const arma::cx_vec & a,
    const arma::cx_vec & a_derivatives,
    const arma::mat & points,
    const math::Polynomial<cx_double> & V_eff_0,
    const double mass
) {

  arma::vec result(12);

  const auto E_func = E_function(a, a_derivatives, V_eff_0, mass);

  const auto polynomial_term_list = lmat{{{0, 0, 1, 1, 2, 0}, {2, 1, 1, 0, 0, 0}}};

#pragma omp parallel for
  for (arma::uword i = 0; i < 6; i++) {

    const lvec E_derivative_term = polynomial_term_list.col(i);
    const auto E_derivative =
        math::polynomial::Term<double>(1.0, E_derivative_term);

    result(i) = -arma::sum(2.0 *
                           arma::real(at(E_func.conj() * E_derivative, points))
                           %
                           arma::pow(arma::abs(heller_gaussian(a, points)), 2));

    result(i + 6) = arma::sum(2.0 *
                              arma::imag(
                                  at(E_func.conj() * E_derivative, points))
                              %
                              arma::pow(arma::abs(heller_gaussian(a, points)),
                                        2));
  }

  return result;
}

inline
double I_function_gsl_wrapper(
    const gsl_vector * a_derivatives,
    void * param
) {
  const HellerParam heller_param = *(HellerParam *) param;
  const arma::cx_vec a = heller_param.a;
  const arma::vec a_derivatives_all = gsl::convert_vec(a_derivatives);
  const arma::cx_vec a_derivatives_arma =
      arma::cx_vec{a_derivatives_all.rows(0, 5),
                   a_derivatives_all.rows(6, 11)};

  return I_function(a,
                    a_derivatives_arma,
                    heller_param.points,
                    heller_param.V_eff_0,
                    heller_param.mass);

}

inline
void I_function_derivative_gsl_wrapper(
    const gsl_vector * a_derivatives,
    void * param,
    gsl_vector * g
) {
  const HellerParam heller_param = *(HellerParam *) param;
  const arma::cx_vec a = heller_param.a;
  const arma::vec a_derivatives_all = gsl::convert_vec(a_derivatives);

  const arma::cx_vec a_derivatives_arma =
      arma::cx_vec{a_derivatives_all.rows(0, 5),
                   a_derivatives_all.rows(6, 11)};

  const auto result_pointer =
      gsl::convert_vec(I_function_derivative(a,
                                             a_derivatives_arma,
                                             heller_param.points,
                                             heller_param.V_eff_0,
                                             heller_param.mass));

  gsl_vector_memcpy(g, result_pointer);

  gsl_vector_free(result_pointer);
}

inline
void I_function_fdf_gsl_wrapper(
    const gsl_vector * a_derivatives,
    void * param,
    double * f,
    gsl_vector * g
) {
  I_function_derivative_gsl_wrapper(a_derivatives, param, g);

  *f = I_function_gsl_wrapper(a_derivatives, param);
}

inline
arma::cx_vec a_derivative(HellerParam input,
                          const double initial_step_size,
                          const double tolerance,
                          const double gradient_tolerance,
                          const size_t total_steps) {

  /* allocate memory for minimization process */
  const auto minimizer_type = gsl_multimin_fdfminimizer_vector_bfgs2;
  auto minimizer_environment = gsl_multimin_fdfminimizer_alloc(minimizer_type,
                                                               12);

  /* assigning function to minimizer object */
  gsl_multimin_function_fdf minimizer_object;
  minimizer_object.f = &I_function_gsl_wrapper;
  minimizer_object.df = &I_function_derivative_gsl_wrapper;
  minimizer_object.fdf = &I_function_fdf_gsl_wrapper;
  minimizer_object.n = 12;
  minimizer_object.params = (void *) &input;

  /* starting point */
  const auto a_derivatives = gsl_vector_calloc(12);

  /* set environment */
  gsl_multimin_fdfminimizer_set(minimizer_environment,
                                &minimizer_object, a_derivatives,
                                initial_step_size, tolerance);

  size_t iter = 0;
  int status = GSL_CONTINUE;
  do {
    iter++;
    status = gsl_multimin_fdfminimizer_iterate(minimizer_environment);

    if (status) {
      throw Error(gsl_strerror(status));
    }

    status = gsl_multimin_test_gradient(minimizer_environment->gradient,
                                        gradient_tolerance);

    if (status == GSL_SUCCESS) {
      const arma::vec result = gsl::convert_vec(minimizer_environment->x);

      gsl_multimin_fdfminimizer_free(minimizer_environment);
      gsl_vector_free(a_derivatives);

      return arma::cx_vec{result.rows(arma::span(0, 5)),
                          result.rows(arma::span(6, 11))};
    }
  } while (status == GSL_CONTINUE && iter < total_steps);

  return arma::ones<arma::cx_vec>(6) * 2.8375;
}

} // namespace details

struct State {
public:
  arma::mat points;
  arma::vec weights;
  arma::vec masses;

  // Establish an easy way to construct your State
  template<typename PhaseSpaceDistribution>
  State(const PhaseSpaceDistribution & initial,
        const arma::uvec & grid,
        const arma::mat & range,
        const arma::vec & masses) :
      points(math::space::points_generate(grid, range)),
      weights(arma::real(at(initial, points))),
      masses(masses) {
    if (grid.n_rows != range.n_rows) {
      throw Error("Different dimension between the grid and the range");
    }
    if (grid.n_rows != 2 * masses.n_rows) {
      throw Error("Different dimension between the grid and the masses");
    }
  }

  template<typename PhaseSpaceDistribution>
  State(const PhaseSpaceDistribution & initial,
        const arma::uvec & grid,
        const arma::mat & range) :
      points(math::space::points_generate(grid, range)),
      weights(arma::real(at(initial, points))),
      masses(arma::ones<arma::vec>(grid.n_elem / 2)) {
    if (grid.n_elem % 2 != 0) {
      throw Error("Odd number of dimension - it is not likely a phase space");
    }
    if (grid.n_rows != range.n_rows) {
      throw Error("Different dimension between the grid and the range");
    }
    if (grid.n_rows != 2 * masses.n_rows) {
      throw Error("Different dimension between the grid and the masses");
    }
  }

  inline
  State(const arma::mat & points,
        const arma::vec & weights,
        const arma::vec & masses) :
      points(points),
      weights(weights),
      masses(masses) {
    if (points.n_cols != weights.n_elem) {
      throw Error("Different number of points and corresponding weights");
    }

    if (points.n_rows != 2 * masses.n_rows) {
      throw Error("Different dimension between the points and the masses");
    }
  }

  inline
  arma::uword dim() const {
    return points.n_rows / 2;
  }

  inline
  State normalise() const {
    return State(this->points, this->weights / arma::sum(this->weights),
                 this->masses);
  }

  template<typename Function>
  arma::vec expectation(const std::vector<Function> & function) const {
    arma::vec result(function.size());

#pragma omp parallel for
    for (arma::uword i = 0; i < result.n_elem; i++) {

      if (function[i].dim() != this->dim() * 2) {
        throw Error(
            "The dimension of the function is invalid for the calculation of expectation");
      }
      result(i) = arma::dot(at(function[i], this->points), weights) /
                  arma::sum(weights);
    }

    return result;
  }

  template<typename Function>
  double expectation(const Function & function) const {
    const arma::vec result = at(function, this->points);

    return arma::dot(result, weights) / arma::sum(weights);
  }

  inline
  arma::vec positional_expectation() const {
    arma::uword dim = this->dim();

    return this->points.rows(0, dim - 1) * this->weights /
           arma::sum(this->weights);
  }

  inline
  arma::vec momentum_expectation() const {
    arma::uword dim = this->dim();

    return this->points.rows(dim, 2 * dim - 1) * this->weights /
           arma::sum(this->weights);
  }

  State operator+(const State & B) const {
    if (!arma::approx_equal(this->weights, B.weights, "abs_diff", 1e-16) ||
        !arma::approx_equal(this->masses, B.masses, "abs_diff", 1e-16)) {
      throw Error("Different cwa states are being added");
    }

    return State(this->points + B.points, this->weights, this->masses);
  }

  State operator*(const double B) const {
    return State(this->points * B, this->weights, this->masses);
  }
};

template<typename Potential>
struct Operator {

private:
  PropagationType type = Classic;

public:
  Potential potential;

  Operator(const State & state,
           const Potential & potential) :
      potential(potential) {}

  inline
  PropagationType propagation_type() const {
    return Classic;
  }

  State operator()(const State & state) const {

    arma::mat p_submatrix = state.points.rows(state.dim(), 2 * state.dim() - 1);
    p_submatrix.each_col() /= state.masses;

    const arma::mat change_list =
        arma::join_cols(p_submatrix,
                        details::force(potential,
                                       state.points.rows(0, state.dim() - 1)));

    return State(change_list, state.weights, state.masses);
  }

};

} // namespace heller
}

#endif //METHODS_HELLER_H

max_stars_repo_name: Walter-Feng/Quartz

max_issues_repo_name: Walter-Feng/Quartz

max_forks_repo_name: Walter-Feng/Quartz


text:
/*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*
**                                                                                **
** This file forms part of the Underworld geophysics modelling application.         **
**                                                                                **
** For full license and copyright information, please refer to the LICENSE.md file  **
** located at the project root, or contact the authors.                             **
**                                                                                **
**~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*/
#ifndef __BSSCR_SUM_h__
#define __BSSCR_SUM_h__

#include <petsc.h>
#include <petscmat.h>
#include <petscvec.h>
#include <petscksp.h>
#include <petscpc.h>
#include <StGermain/libStGermain/src/StGermain.h>
#include <StgDomain/libStgDomain/src/StgDomain.h>
#include <StgFEM/libStgFEM/src/StgFEM.h>
#include <PICellerator/libPICellerator/src/PICellerator.h>
#include <Underworld/libUnderworld/src/Underworld.h>
#include "Solvers/KSPSolvers/src/KSPSolvers.h"

#include "common-driver-utils.h"
#include "BSSCR.h" /* includes StokesBlockKSPInterface.h */

//void bsscr_summary(KSP_BSSCR * bsscrp_self, KSP ksp_S, KSP ksp_inner, Mat K,Mat Korig,Mat K2,Mat D,Mat G,Mat C,Vec u,Vec p,Vec f,Vec h,Vec t,

void bsscr_summary(KSP_BSSCR * bsscrp_self, KSP ksp_S, KSP ksp_inner, Mat K,Mat K2,Mat D,Mat G,Mat C,Vec u,Vec p,Vec f,Vec h,Vec t,
		   double penaltyNumber,PetscTruth KisJustK,double mgSetupTime, double RHSSolveTime, double scrSolveTime,double a11SingleSolveTime);

#endif


text:
/* #################################################
 * # stochastic simulation of adapting populations #
 * # with varying constraint type                  #
 * # and exponential mutation kernel               #
 * #################################################
 *
 * 2012-2015, Lukas Geyrhofer, lukas.geyrhofer@ds.mpg.de
 *
 *
 * #################################################
 *
 * simplest usage (population size 1e7):
 * ./travelingwavepeak_exp -N 1e7 1> out.txt 2> conf.txt
 *
 * ################################################# */

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

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

#include <mcheck.h>

int space = 300;
int space0 = 100;
double dx = 1e-2;
int maxSteps = 1000;

int outputstep = 100;
int quiet = 0;

int noise = 0;		// noise can assume 4 different values:
			// value in code	cmdline	description
			// 0 (default)			populationsize is constrained to 1, no noise is added
			//				constraint could be any value
			// 1			N	fixedN, set by parameter -N POPSIZE
			// 2			u U	ustar, set by parameter -u UFILENAME
			//				if ustar-profile is in finer grid than simulation, use parameter -U LATTICERATIO
			// 3			n g	timetrace, population size is constrained using a (single column) textfile,
			//				where the entries are the current population size
			//				use also parameter -g STEPS to indicate the number of timesteps (steps! not time in generations)
			//				between two values in this textfile
			// 4			y U	ustar, use file from -y UFILENAME
			//				no stochastic dynamics, noise is not added, only the constraint is used
			//				if ustar-profile is in finer grid than simulation, use parameter -U LATTICERATIO
double populationsize = 1.;
double populationvariance;
double xg,varg,xcmax,xnose;
int correctformeanfitness =0;
int printhistotype = 1;
int printfixationprob = 0;

int allshifts = 0;
int shiftthreshold = 1;
double current_mean_fitness = 0.;

int read_from_file = 0;
int write_to_file = 0;
char c_infile[128],c_outfile[128];

int have_u_infile = 0;
char u_infile[128];
double *u_read,*u;
int dens_ustar_latticeratio = 1;
int u_space,u_space0;
double u_dx;
double wavespeed = 0.;
double speedprefactor;

double epsilon = 1e-2;
double twoepssqrt;
double *nn;
double *tmp;
double *x;

double mutationrate = 1e-5;
double *mutation_inflow;
double mutation_outflow;
double mutation_sigma = 1e-2;

gsl_rng* rg; /* gsl, global generator */
const gsl_rng_type* T;
unsigned long int randseed = 0;

int averagepopdens = 0;
int averagepopdens_center = 1;
int averagepopdens_resolution = 1;
int averagepopdens_havefile = 0;
double *averagepopdens_dens;
char averagepopdens_outputfile[128];
double averagepopdens_count = 0.;
double averagepopdens_dx;
int averagepopdens_space;
int averagepopdens_space0;
int averagepopdens_lower;
int averagepopdens_higher;

double *timetrace_popsize;
char timetrace_filename[128];
int timetrace_ps_interval = 100;
int timetrace_length;

int print_error(char *msg) {
  fprintf(stderr,"ERROR: %s\n",msg);
  exit(1);
}

// ************************************************************
// **   parameters
// ************************************************************

void parsecomamndline(int argn, char *argv[]) {
  char c;
  while((c = getopt(argn,argv,"s:z:d:S:e:D:O:qQi:o:R:N:CM:u:U:T:H:h:Pn:g:y:F")) != -1) {
    switch(c) {
      case 's':	space = atoi(optarg);
		break;
      case 'z':	space0 = atoi(optarg);
		break;
      case 'd':	dx = atof(optarg);
		break;
      case 'i':	strcpy(c_infile,optarg);
		read_from_file = 1;
		break;
      case 'o':	strcpy(c_outfile,optarg);
		write_to_file = 1;
		break;
      case 'u':	strcpy(u_infile,optarg);
		if(noise > 0)print_error("Only a single constraint-type can be used (either option -N POPSIZE or -u FILENAME or -n FILENAME_TIMETRACE)");
		noise = 2;
		break;
      case 'y':	strcpy(u_infile,optarg);
		if(noise > 0)print_error("Only a single constraint-type can be used (either option -N POPSIZE or -u FILENAME or -n FILENAME_TIMETRACE or -y FILENAME)");
		noise = 4;
		break;
      case 'U':	if((noise == 2)||(noise == 4)) {
		  dens_ustar_latticeratio = atoi(optarg);
		}else{
		  print_error("option -u FILENAME needed before option -U RATIO");
		}
		break;
      case 'S':	maxSteps = atoi(optarg);
		break;
      case 'e': epsilon = atof(optarg);
		break;
      case 'D':	mutationrate = atof(optarg);
		break;
      case 'O':	outputstep = atoi(optarg);
		break;
      case 'q':	quiet = 2;
		break;
      case 'Q':	quiet = 1;
		break;
      case 'R':	randseed =atoi(optarg);
		break;
      case 'N': populationsize = atof(optarg);
		if(noise > 0)print_error("Only a single constraint-type can be used (either option -N POPSIZE or -u FILENAME or -n FILENAME_TIMETRACE)");
		noise = 1;
		break;
      case 'C':	correctformeanfitness = 1;
		break;
      case 'M':	mutation_sigma = atof(optarg);
		break;
      case 'T':	shiftthreshold = atoi(optarg);
		break;
      case 'h':	strcpy(averagepopdens_outputfile,optarg);
		averagepopdens_havefile = 1;
		averagepopdens = 1;
		break;
      case 'H':	averagepopdens = 1;
		averagepopdens_resolution = atoi(optarg);
		if(averagepopdens_resolution < 0) {
		  averagepopdens_resolution *= -1;
		  averagepopdens_center = 0;
		}
		break;
      case 'P':	printhistotype = 0;
		break;
      case 'F': printfixationprob = 1;
		break;
      case 'n':	strcpy(timetrace_filename,optarg);
		if(noise > 0)print_error("Only a single constraint-type can be used (either option -N POPSIZE or -u FILENAME or -n FILENAME_TIMETRACE)");
		noise = 3;
		break;
      case 'g':	timetrace_ps_interval = atoi(optarg);
		break;
      }
  }
  if(randseed==0)randseed=time(NULL);
}

// ************************************************************
// **   input and output for configuration files
// ************************************************************

void read_popdens(int importparameters) {
  int i;
  int icount,dcount;
  int *ival;
  double *dval;
  int c_space,c_space0;
  double c_dx;
  FILE *fp;

  fp = fopen(c_infile,"rb");
  if(fp != NULL) {
    fread(&icount,sizeof(int),1,fp);
    fread(&dcount,sizeof(int),1,fp);
    fread(&c_dx,sizeof(double),1,fp);
    fread(&c_space,sizeof(int),1,fp);
    fread(&c_space0,sizeof(int),1,fp);

    if(importparameters) {
      space = c_space;
      space0 = c_space0;
      dx = c_dx;
    }

    if(icount>0) {
      ival = (int*)malloc(icount*sizeof(int));
      fread(ival,sizeof(int),icount,fp);
      free(ival);
    }

    if(dcount>0) {
      dval = (double*)malloc(dcount*sizeof(int));
      fread(dval,sizeof(double),dcount,fp);
      free(dval);
    }

    if(space != c_space)print_error("lattice does not match!");
    nn = (double*)malloc(space*sizeof(double));
    fread(nn,sizeof(double),space,fp);
    for(i=0;i<space;i++) {
      nn[i] *= dx;
    }

    fclose(fp);
  }else{
    print_error("could not open c-infile");
  }
}

void write_popdens() {
  int icount=0,dcount=0;
  int i;
  FILE *fpc;

  fpc=fopen(c_outfile,"wb");
  if(fpc != NULL) {
    fwrite(&icount,sizeof(int),1,fpc);
    fwrite(&dcount,sizeof(int),1,fpc);
    fwrite(&dx,sizeof(double),1,fpc);
    fwrite(&space,sizeof(int),1,fpc);
    fwrite(&space0,sizeof(int),1,fpc);
    for(i=0;i<space;i++)nn[i] /= dx;
    fwrite(nn,sizeof(double),space,fpc);
    fclose(fpc);
  }else{
    print_error("could not open c-outfile");
  }
}

void read_constraint(int importparameters) {
  int icount,dcount;
  int *ival;
  double *dval;
  FILE *fp;

  fp = fopen(u_infile,"rb");
  if(fp!=NULL) {
    fread(&icount,sizeof(int),1,fp);
    fread(&dcount,sizeof(int),1,fp);
    fread(&u_dx,sizeof(double),1,fp);
    fread(&u_space,sizeof(int),1,fp);
    fread(&u_space0,sizeof(int),1,fp);

    if(icount>0) {
      ival = (int*)malloc(icount*sizeof(int));
      fread(ival,sizeof(int),icount,fp);
      free(ival);
    }

    if(dcount>=3) {
      dval = (double*)malloc(dcount*sizeof(int));
      fread(dval,sizeof(double),dcount,fp);
      mutationrate = dval[0];
      wavespeed = dval[1];
      mutation_sigma = dval[2];
      free(dval);
    }else{
      print_error("not enough values in constraint file! need at least 3 double parameters: mutationrate, wavespeed, mutationsigma!");
    }

    if(importparameters) {
      space = u_space/dens_ustar_latticeratio;
      space0 = u_space0/dens_ustar_latticeratio;
      dx = u_dx*dens_ustar_latticeratio;
    }else{
      if((space*dens_ustar_latticeratio != u_space)||(space0*dens_ustar_latticeratio != u_space0))print_error("lattice does not match! u");
    }

    u_read = (double*)malloc(u_space*sizeof(double));
    fread(u_read,sizeof(double),u_space,fp);

    fclose(fp);
  }else{
    print_error("could not open ufile");
  }

  u = (double*)malloc(space*sizeof(double));

}

void flat_constraint(double size) {
  int i;
  u=(double*)malloc(space*sizeof(double));
  for(i=0;i<space;i++)u[i] = 1./size;
}

void initialize_with_gaussian_popdens() {
  int i;
  double startvariance;
  if((noise == 2)||(noise == 4)) {
    startvariance = wavespeed;
  }else{
    // Good et al., PNAS (2012)
    startvariance = mutation_sigma*mutation_sigma*2*log(mutation_sigma*populationsize)/(log(mutation_sigma/mutationrate)*log(mutation_sigma/mutationrate));
    // behaves weird for extremal parameter values.
    // catch those cases and set on standard value: density takes half of the box
    if(startvariance < 0) {
      startvariance = (0.25*space*space*dx*dx);
    }
  }
  for(i=0;i<space;i++) {
    nn[i] = exp(-(i-space0)*(i-space0)*dx*dx/(2.*startvariance));
  }
}

void print_u(int time) {
  int i;
  for(i=0;i<space;i++) {
    fprintf(stderr,"%10.6lf %10.6lf %16.10e\n",time*epsilon,(i-space0+allshifts)*dx,u[i]);
  }
  fprintf(stderr,"\n");
}

// ************************************************************
// **   timetrace
// ************************************************************

void read_timetrace() {
  int i;
  timetrace_length = (int)(maxSteps/timetrace_ps_interval)+1;
  FILE *fp;
  timetrace_popsize = (double*)malloc(timetrace_length*sizeof(double));
  fp = fopen(timetrace_filename,"r");
  if(fp != NULL) {
    for(i=0;i<timetrace_length;i++) {
      fscanf(fp,"%lf",&timetrace_popsize[i]);
    }
  }else{
    print_error("could not open timetrace file");
  }
  fclose(fp);
}

double get_timetrace_popsize(int step) {
  int i;
  if(step%timetrace_ps_interval == 0) {
    return timetrace_popsize[step/timetrace_ps_interval];
  }else{
    int lasttimepoint = (int)floor(step/timetrace_ps_interval);
    int nexttimepoint = lasttimepoint + 1;
    double fractimepoint = (1.*(step%timetrace_ps_interval))/(1.*timetrace_ps_interval);
    return fractimepoint*timetrace_popsize[nexttimepoint] + (1-fractimepoint)*timetrace_popsize[lasttimepoint];
  }
}

// ************************************************************
// **   screen output
// ************************************************************

void print_populationdensity(int timestep) {
  int i;
  double corr = allshifts*dx;
  if(correctformeanfitness)corr += current_mean_fitness;
  for(i=0;i<space;i++) {
    fprintf(stderr,"%lf %14.10lf %20.10e",timestep*epsilon,(i-space0)*dx+corr,nn[i]);
    if(printfixationprob)fprintf(stderr," %20.10e",u[i]);
    fprintf(stderr,"\n");

    if(printhistotype) {
      fprintf(stderr,"%lf %14.10lf %20.10e",timestep*epsilon,(i-space0+1)*dx+corr,nn[i]);
      if(printfixationprob)fprintf(stderr," %20.10e",u[i]);
      fprintf(stderr,"\n");
    }
  }
  fprintf(stderr,"\n");
}

// ************************************************************
// **   average popdens
// ************************************************************

void init_averagepopdens() {

  averagepopdens_dx = dx/(1.*averagepopdens_resolution);
  averagepopdens_space = space*averagepopdens_resolution;
  averagepopdens_space0 = space0*averagepopdens_resolution;

  averagepopdens_dens = (double*)calloc(averagepopdens_space,sizeof(double));

  averagepopdens_lower = (int)(-averagepopdens_resolution/2);
  averagepopdens_higher = (int)(averagepopdens_resolution/2);
  if(averagepopdens_higher - averagepopdens_lower < averagepopdens_resolution)averagepopdens_higher++;

}

void update_averagepopdens() {
  int i,j,idx;
  int offset=0;

  if(averagepopdens_center) {
    offset = (int)(current_mean_fitness/dx*averagepopdens_resolution);
  }

  for(i=0;i<space;i++) {
    for(j=averagepopdens_lower;j<averagepopdens_higher;j++) {
      idx = i*averagepopdens_resolution+j-offset;
      if((averagepopdens_space > idx) &&( idx >= 0))
	averagepopdens_dens[idx] += nn[i];
    }
  }

  averagepopdens_count += 1.;

}

void write_averagepopdens() {
  int i;
  FILE *fp;

  if(averagepopdens_havefile == 1) {
    fp = fopen(averagepopdens_outputfile,"w");
  }else{
    fp = stdout;
  }

  for(i=0;i<averagepopdens_space;i++) {
    fprintf(fp,"%lf %e\n",(i-averagepopdens_space0)*averagepopdens_dx,averagepopdens_dens[i]/averagepopdens_count);
  }

  if(averagepopdens_havefile == 1) {
    fclose(fp);
  }

  free(averagepopdens_dens);
}

// ************************************************************
// **   main algorithm
// ************************************************************

void update_u(int timestep) {
   // bug if shiftthreshold -T >1, use shiftthreshold only with fixedN
  int i,j;
  int baseshift;
  double fracshift;

  baseshift = (int)floor(current_mean_fitness/u_dx);
  fracshift = current_mean_fitness/u_dx - 1.*baseshift;
  i = 0;
  for(j=0;j<u_space;j++) {
    if ( j + baseshift + fracshift > (i+.5)*dens_ustar_latticeratio) {
      u[i] += (1.-fracshift)*u_read[j];
      u[i] /= (1.*dens_ustar_latticeratio);
      i++;
      if(i>=space)break;
      u[i] = fracshift*u_read[j];
    }else{
      u[i] += u_read[j];
    }
  }
}

void shift_population_backward(int step) {
  int i;
  for(i=0;i<space-step;i++)     nn[i] = nn[i+step];
  for(i=space-step;i<space;i++) nn[i] = 0.;
}

void shift_population_forward(int step) {
  int i;
  for(i=space-1;i>-step;i--) nn[i] = nn[i+step];
  for(i=-step;i>=0;i--)      nn[i] = 0.;
}

void shift_population(int timestep) {
  int shift = (int)floor(current_mean_fitness/dx);
  if(shift >= shiftthreshold) {shift_population_backward(shift);}
  if(shift <= -shiftthreshold) {shift_population_forward(-shift);}
  allshifts += shift;
  current_mean_fitness -= shift*dx;
}

void reproduce(int timestep) {
  int i,j;
  double tmpn;
  long int poiss;

  if((noise == 2)||(noise==4)) {
    current_mean_fitness = timestep*wavespeed*epsilon - allshifts*dx;
  }	// else:	current_mean_fitness was already updated at end of last timestep
	//		with function "populationconstraint" below and is still at the correct value

  if((current_mean_fitness/dx >= shiftthreshold) || (current_mean_fitness/dx <= -shiftthreshold)) {
    shift_population(timestep);
  }

  tmp[0] = 0.;
  tmp[space-1] = 0.;
  for(i=1;i<space-1;i++) {
    tmpn = nn[i]*(1.+(x[i]-current_mean_fitness)*epsilon-mutation_outflow);
    for(j = 1;j<i;j++) tmpn += mutation_inflow[j]*nn[i-j];
    if(tmpn<0) {
      tmpn = 0.;
    }else if((noise > 0)&&(noise != 4)) {
      if(tmpn < 1.e9) {	//	Poisson-RNG breaks down for parameters > 1e9. see GSL doc.
			//	use smaller bins if this occurs too often or population size too large
	tmp[i] = tmpn + twoepssqrt*(gsl_ran_poisson(rg,tmpn)-tmpn);
      }else{		//	neglect noise term if occupancy is too large,
	tmp[i] = tmpn;	//	dynamics can be assumed to be almost deterministic then
      }
    }else{
      tmp[i] = tmpn;
    }
  }

  memcpy(&nn[1],&tmp[1],(space-2)*sizeof(double));

}

double populationconstraint(int timestep) {
  int i;
  double constraint = 0., inv, sumxg = 0., sumxxg = 0.;

  double popdens_0thmom = 0.;
  double popdens_1stmom = 0.;
  double popdens_2ndmom = 0.;
  double cmax = 0.;
  int havenose = 0;
  if((noise == 2)||(noise == 4)) {
    update_u(timestep);
  }

  for(i=space-1;i>=0;i--) {
    if (nn[i] > cmax) {
      cmax=nn[i];
      xcmax = x[i];
    }
    if((havenose==0)&&(nn[i]>1)) {
      havenose = 1;
      xnose = x[i];
    }
    constraint += nn[i]*u[i];
    sumxg += x[i] * nn[i] * u[i];
    sumxxg += x[i] * x[i] * nn[i] * u[i];
    popdens_0thmom += nn[i];
    popdens_1stmom += x[i]*nn[i];
    popdens_2ndmom += x[i]*x[i]*nn[i];
  }

  current_mean_fitness = popdens_1stmom/popdens_0thmom;
  populationvariance = popdens_2ndmom/popdens_0thmom - current_mean_fitness*current_mean_fitness;
  populationsize = popdens_0thmom;
  xg = sumxg/constraint;
  varg = sumxxg/constraint - xg*xg;

  inv = 1./constraint;
  if(noise==3) {
    inv *= get_timetrace_popsize(timestep);
  }
  for(i=0;i<space;i++) nn[i] *= inv;
  populationsize *= inv;

  return inv;
}

// ************************************************************
// **   initialization
// ************************************************************

int initialize() {
  int i;

  if(read_from_file) {
    read_popdens(1);
    if(noise == 0)flat_constraint(1.);
    if(noise == 1)flat_constraint(populationsize);
    if((noise == 2)||(noise == 4)) {
      read_constraint(0);
      update_u(0);
    }
    if(noise == 3) {
      flat_constraint(1.);
      read_timetrace();
    }
  }else{
    if(noise == 0)flat_constraint(1.);
    if(noise == 1)flat_constraint(populationsize);
    if((noise == 2) || (noise == 4)) {
      read_constraint(1);
      update_u(0);
    }
    if(noise == 3) {
      flat_constraint(1.);
      read_timetrace();
    }
    nn = (double*)malloc(space*sizeof(double));
    initialize_with_gaussian_popdens();
  }
  tmp = (double*)malloc(space*sizeof(double));

  gsl_rng_env_setup();
  T = gsl_rng_default;
  rg = gsl_rng_alloc(T);
  gsl_rng_set(rg, randseed);

  twoepssqrt = sqrt(2.*epsilon);
  speedprefactor = wavespeed/dx*epsilon;

  x = (double*)malloc(space*sizeof(double));
  mutation_inflow = (double*)malloc(space*sizeof(double));
  mutation_outflow = 0.;
  mutation_inflow[0] = 0.;
  for(i=1;i<space;i++) {
    x[i] = (i-space0)*dx;
    mutation_inflow[i] = epsilon*mutationrate*exp(-i*dx/mutation_sigma)*(1.-exp(-dx/mutation_sigma));
    mutation_outflow += mutation_inflow[i];
  }

  if(quiet<2) {
    printf("#################################################################################\n");
    printf("# stochastic simulation of adapting population with exponential mutation kernel #\n");
    printf("#################################################################################\n");
    printf("#    mutationrate     = %e\n",mutationrate);
    printf("#    mutation_sigma   = %e\n",mutation_sigma);
    if(noise==0) {
      printf("#    constraint       = deterministic\n");
    }
    if(noise==1) {
      printf("#    constraint       = fixedN\n");
      printf("#    populationsize   = %e\n",populationsize);
    }
    if(noise==2) {
      printf("#    constraint       = ustar\n");
      printf("#    ufile            = %s\n",u_infile);
      printf("#    wavespeed        = %e\n",wavespeed);
    }
    if(noise==3) {
      printf("#    constraint       = fixedN, timetrace\n");
      printf("#    timetracefile    = %s\n",timetrace_filename);
    }
    if(noise==4) {
      printf("#    constraint       = ustar, non-stochastic\n");
      printf("#    ufile            = %s\n",u_infile);
      printf("#    wavespeed        = %e\n",wavespeed);
    }
    printf("#    (lattice) space  = %d\n",space);
    printf("#    (lattice) space0 = %d\n",space0);
    printf("#    (lattice) dx     = %e\n",dx);
    printf("#    randseed         = %d\n",randseed);
  }

}

// ************************************************************
// **   cleanup
// ************************************************************

void cleanup() {
  free(nn);
  free(tmp);
  free(x);
  free(mutation_inflow);
  free(u);
  if((noise == 2)||(noise == 4)) {
    free(u_read);
  }
  if(noise == 3) {
    free(timetrace_popsize);
  }
  gsl_rng_free(rg);
}

// ************************************************************
// **   main
// ************************************************************

int main(int argn, char *argv[]) {
  int i,j;
  double v,lambda = 1.;

  parsecomamndline(argn,argv);
  initialize();
  lambda = populationconstraint(0);
  if (quiet<2) fprintf(stdout,"%10.3lf %20.10e %.10e %.10e %.10e %.10e %.10e %.10e %.10e\n",0.0,current_mean_fitness,populationvariance,populationsize,lambda,xg,varg,allshifts*dx+xcmax,allshifts*dx+xnose);
  if(averagepopdens) {
    init_averagepopdens();
  }
  for(i=1;i<=maxSteps;i++) {
    reproduce(i);
    lambda = populationconstraint(i);

    if(i%outputstep == 0) {
      if (quiet<2)fprintf(stdout,"%10.3lf %20.10e %.10e %.10e %.10e %.10e %.10e %.10e %.10e\n",i*epsilon,allshifts*dx + current_mean_fitness,populationvariance,populationsize,lambda,xg,varg,allshifts*dx+xcmax,allshifts*dx+xnose);
      if (quiet==0) print_populationdensity(i);
      if (averagepopdens) update_averagepopdens();
    }
  }

  if(write_to_file)write_popdens();
  if(averagepopdens)write_averagepopdens();
  cleanup();
  return 0;
}

max_stars_repo_name: lukasgeyrhofer/adaptivewaves

max_issues_repo_name: lukasgeyrhofer/adaptivewaves

max_forks_repo_name: lukasgeyrhofer/adaptivewaves


text:
#if !defined(incompressibleFlow_h)
#define incompressibleFlow_h
#include <petsc.h>

// Define the test field number
typedef enum { VTEST, QTEST, WTEST } LowMachFlowTestFields;

typedef enum { STROUHAL, REYNOLDS, PECLET, MU, K, CP, TOTAL_INCOMPRESSIBLE_FLOW_PARAMETERS } IncompressibleFlowParametersTypes;
PETSC_EXTERN const char* incompressibleFlowParametersTypeNames[TOTAL_INCOMPRESSIBLE_FLOW_PARAMETERS + 1];

typedef enum { VEL, PRES, TEMP, TOTAL_INCOMPRESSIBLE_FLOW_FIELDS } IncompressibleFlowFields;

typedef enum { MOM, MASS, ENERGY, TOTAL_INCOMPRESSIBLE_SOURCE_FIELDS } IncompressibleSourceFields;

// Incompressible Kernels

PETSC_EXTERN void IncompressibleFlow_vIntegrandTestFunction(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[],
                                                            const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[],
                                                            const PetscScalar a_x[], PetscReal t, const PetscReal X[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]);

PETSC_EXTERN void IncompressibleFlow_vIntegrandTestGradientFunction(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[],
                                                                    const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[],
                                                                    const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal X[], PetscInt numConstants,
                                                                    const PetscScalar constants[], PetscScalar f1[]);

PETSC_EXTERN void IncompressibleFlow_wIntegrandTestFunction(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[],
                                                            const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[],
                                                            const PetscScalar a_x[], PetscReal t, const PetscReal X[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]);
PETSC_EXTERN void IncompressibleFlow_wIntegrandTestGradientFunction(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[],
                                                                    const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[],
                                                                    const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal X[], PetscInt numConstants,
                                                                    const PetscScalar constants[], PetscScalar f1[]);

PETSC_EXTERN void IncompressibleFlow_qIntegrandTestFunction(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[],
                                                            const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[],
                                                            const PetscScalar a_x[], PetscReal t, const PetscReal X[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]);
PETSC_EXTERN void IncompressibleFlow_g1_qu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[],
                                           const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
                                           PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]);
PETSC_EXTERN void IncompressibleFlow_g0_vu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[],
                                           const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
                                           PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]);
PETSC_EXTERN void IncompressibleFlow_g1_vu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[],
                                           const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
                                           PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]);
PETSC_EXTERN void IncompressibleFlow_g2_vp(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[],
                                           const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
                                           PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]);
PETSC_EXTERN void IncompressibleFlow_g3_vu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[],
                                           const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
                                           PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]);
PETSC_EXTERN void IncompressibleFlow_g0_wT(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[],
                                           const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
                                           PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]);
PETSC_EXTERN void IncompressibleFlow_g0_wu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[],
                                           const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
                                           PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]);
PETSC_EXTERN void IncompressibleFlow_g1_wT(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[],
                                           const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
                                           PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]);
PETSC_EXTERN void IncompressibleFlow_g3_wT(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[],
                                           const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
                                           PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]);
#endif

max_stars_repo_name: mtmcgurn-buffalo/ablate

max_issues_repo_name: mtmcgurn-buffalo/ablate

max_forks_repo_name: mtmcgurn-buffalo/ablate


text:
/*
 ============================================================================
 Name        : mathutils.h
 Author      : Davide Malvezzi
 Version     : 1.0
 Copyright   :
 Description : Math utils function
 ============================================================================
 */

#ifndef MATH_UTILS_H
#define MATH_UTILS_H

#include <float.h>
#include <limits.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <stdint.h>
#include <gsl/gsl_rng.h>

#ifndef M_PI
	#define M_PI	3.14159265359
#endif

#ifndef max
	#define max( a, b ) ( ((a) > (b)) ? (a) : (b) )
#endif

#ifndef min
	#define min( a, b ) ( ((a) < (b)) ? (a) : (b) )
#endif

#ifndef swapv
	#define swapv(x, y) do { __typeof__(x) SWAP = x; x = y; y = SWAP; } while (0)
#endif

extern gsl_rng* rng;

extern void randInit();
extern int randInt(int min, int max);
extern float randFloat(float min, float max);

extern int isZero(float val);

extern double fast_sin(double x);

extern double fast_cos(double x);

extern double fast_atan(double x);

#define sin_f	fast_sin
#define cos_f	fast_cos
#define atan_f	fast_atan

#endif

max_stars_repo_name: DavideMalvezzi/ChangePointGA

max_issues_repo_name: DavideMalvezzi/ChangePointGA

max_forks_repo_name: DavideMalvezzi/ChangePointGA


text:
#ifndef libceed_solids_examples_setup_libceed_h
#define libceed_solids_examples_setup_libceed_h

#include <ceed.h>
#include <petsc.h>
#include "../include/structs.h"

// -----------------------------------------------------------------------------
// libCEED Functions
// -----------------------------------------------------------------------------
// Destroy libCEED objects
PetscErrorCode CeedDataDestroy(CeedInt level, CeedData data);

// Utility function - essential BC dofs are encoded in closure indices as -(i+1)
PetscInt Involute(PetscInt i);

// Utility function to create local CEED restriction from DMPlex
PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt height,
    DMLabel domain_label, CeedInt value, CeedElemRestriction *elem_restr);

// Utility function to get Ceed Restriction for each domain
PetscErrorCode GetRestrictionForDomain(Ceed ceed, DM dm, CeedInt height,
                                       DMLabel domain_label, PetscInt value,
                                       CeedInt Q, CeedInt q_data_size,
                                       CeedElemRestriction *elem_restr_q,
                                       CeedElemRestriction *elem_restr_x,
                                       CeedElemRestriction *elem_restr_qd_i);

// Set up libCEED for a given degree
PetscErrorCode SetupLibceedFineLevel(DM dm, DM dm_energy, DM dm_diagnostic,
                                     Ceed ceed, AppCtx app_ctx,
                                     CeedQFunctionContext phys_ctx,
                                     ProblemData problem_data,
                                     PetscInt fine_level, PetscInt num_comp_u,
                                     PetscInt U_g_size, PetscInt U_loc_size,
                                     CeedVector force_ceed,
                                     CeedVector neumann_ceed, CeedData *data);

// Set up libCEED multigrid level for a given degree
PetscErrorCode SetupLibceedLevel(DM dm, Ceed ceed, AppCtx app_ctx,
                                 ProblemData problem_data, PetscInt level,
                                 PetscInt num_comp_u, PetscInt U_g_size,
                                 PetscInt U_loc_size, CeedVector fine_mult,
                                 CeedData *data);

#endif // libceed_solids_examples_setup_libceed_h

max_stars_repo_name: wence-/libCEED

max_issues_repo_name: wence-/libCEED

max_forks_repo_name: wence-/libCEED


text:
/* cdf/geometric.c
 *
 * Copyright (C) 2004 Jason H. Stover.
 * Copyright (C) 2010 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307, USA.
 */

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

#include "error.h"

/* Pr (X <= k), i.e., the probability of n or fewer failures until the
   first success. */

double
gsl_cdf_geometric_P (const unsigned int k, const double p)
{
  double P, a, q;

  if (p > 1.0 || p < 0.0)
    {
      CDF_ERROR ("p < 0 or p > 1", GSL_EDOM);
    }

  if (k < 1)
    {
      return 0.0;
    }

  q = 1.0 - p;
  a = (double) k;

  if (p < 0.5)
    {
      P = -expm1 (a * log1p (-p));
    }
  else
    {
      P = 1.0 - pow (q, a);
    }

  return P;
}

double
gsl_cdf_geometric_Q (const unsigned int k, const double p)
{
  double Q, a, q;

  if (p > 1.0 || p < 0.0)
    {
      CDF_ERROR ("p < 0 or p > 1", GSL_EDOM);
    }

  if (k < 1)
    {
      return 1.0;
    }

  q = 1.0 - p;
  a = (double) k;

  if (p < 0.5)
    {
      Q = exp (a * log1p (-p));
    }
  else
    {
      Q = pow (q, a);
    }

  return Q;
}


text:
/*
 *
 * Author : Pierre Schnizer <schnizer@users.sourceforge.net>
 * Date   : 5. October 2003
 * Only used to test the error handling of pygsl.
 */
#include <Python.h>
#include <gsl/gsl_errno.h>
#include <pygsl/error_helpers.h>

static char trigger_doc [] = "Calls gsl_error with the passed error number";
static PyObject *module;

static PyObject*
trigger(PyObject *self, PyObject *args)
{
	int gsl_errno = GSL_SUCCESS;

	FUNC_MESS_BEGIN();
	if (0 == PyArg_ParseTuple(args, "i", &gsl_errno)){
		PyGSL_add_traceback(module, __FILE__, __FUNCTION__, __LINE__ - 1);
		return NULL;
	}
	/*
	 *
	 * 25. October 2008
	 * set the internal gsl_error_handler to off
	 *
	 */
	pygsl_error("Just a test to see what pygsl is doing!",
		  __FILE__, __LINE__, gsl_errno);
	if (PyGSL_ERROR_FLAG(gsl_errno) != GSL_SUCCESS){
		FUNC_MESS_FAILED();
		return NULL;
	}
	FUNC_MESS_END();
	Py_INCREF(Py_None);
	return (Py_None);
}

static PyMethodDef errortestMethods[] = {
     /*densities*/
     {"trigger", trigger, METH_VARARGS, trigger_doc},
     {NULL, NULL, 0, NULL}
};

DL_EXPORT(void) initerrortest(void)
{
     PyObject  *m=NULL;

     m = Py_InitModule("errortest", errortestMethods);
     assert(m);
     module = m;

     init_pygsl();

     return;
}

max_stars_repo_name: juhnowski/FishingRod

max_issues_repo_name: juhnowski/FishingRod

max_forks_repo_name: juhnowski/FishingRod


text:
#ifndef rpsr_f7f329f9_ae61_4fcf_b0c9_8a96216ed0e0_h
#define rpsr_f7f329f9_ae61_4fcf_b0c9_8a96216ed0e0_h

#include <rathen\config.h>

/*
 * This version of API wrappers were used by the flex & bison due to the cpp generator was still unstable.
 * Any kind of cpp style announcement should not appear in the header file.
 */

#ifndef __cplusplus
#ifdef _UNICODE
typedef wchar_t rchar;
#else
typedef char rchar;
#endif
#else
#include <gslib\string.h>
typedef gs::gchar rchar;
#endif

#ifndef __cplusplus
typedef void rpswrap;
typedef void rpsobj;
typedef void rpsparser;
#else
#include <rathen\parser.h>
typedef gs::rathen::parser::wrapper rpswrap;
typedef gs::rathen::ps_obj rpsobj;
typedef gs::rathen::parser rpsparser;
#endif

/* file: parser.h, line: 14 */
enum rpobj_tag
{
    rpt_root,
    rpt_block,
    rpt_value,
    rpt_expression,
    rpt_operator,
    rpt_function,
    rpt_calling,
    rpt_if_statement,
    rpt_loop_statement,
    rpt_go_statement,
    rpt_stop_statement,
    rpt_exit_statement,
    rpt_variable,
    rpt_constant,
};

enum rpv_tag
{
    rvt_string,
    rvt_integer,
    rvt_float,
};

enum rps_encoding
{
    rpe_ascii,
    rpe_utf8,
    rpe_utf16,
    // ...
};

typedef struct
{
    rps_encoding    encoding;
    const char*     file;
    char*           source;
    char*           current;
    int             length;
    int             line;
    int             position;
}
rgcontext;

typedef struct
{
    rps_encoding    encoding;
    const rchar*    file;
    char*           source;
    char*           current;
    int             length;
    int             line;
    int             position;
    const rchar*    keyword;
    rpsparser*      parser;
    rpswrap*        wrapper;
}
rpscontext;

#ifdef __cplusplus
extern "C" {
#endif

extern void rathen_yy_setup(rgcontext* context, const char* file, rps_encoding encoding = rpe_utf8);
extern void rathen_yy_ssetup(rgcontext* context, const char* file, char* source, int length, rps_encoding encoding = rpe_utf8);
extern void rathen_yy_release(rgcontext* context);
extern rpscontext* rathen_yy_select(rpscontext* context);
extern void rathen_yy_parse(rpsparser* ps, const char* file);
extern void rathen_yy_sparse(rpsparser* ps, const char* file, char* source, int length);
extern int rathen_yy_read(char* buff, int* bytes, int cap);

extern bool rathen_is_type(const rchar* str, int len = -1);
extern bool rathen_is_keyword(const rchar* str, int len = -1);
extern rpsparser* rathen_create_parser(const rchar* name);
extern void rathen_destroy_parser(rpsparser* rps);
extern rpswrap* rathen_get_root(rpsparser* rps);
extern rpswrap* rathen_create_value(rpscontext* context, const rchar* name, const rchar* tyname, bool ref);
extern rpswrap* rathen_create_function(rpscontext* context, const rchar* name, const rchar* retyname, bool retref);
extern rpswrap* rathen_create_constant(rpscontext* context, const rchar* name, rpv_tag tag, const rchar* vstr);
//extern rpswrap* rathen_create();

#ifdef __cplusplus
};
#endif

#endif

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib