text:
/* eigen/gensymmv.c
 *
 * Copyright (C) 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 <config.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>

/*
 * This module computes the eigenvalues and eigenvectors of a real
 * generalized symmetric-definite eigensystem A x = \lambda B x, where
 * A and B are symmetric, and B is positive-definite.
 */

static void gensymmv_normalize_eigenvectors(gsl_matrix *evec);

/*
gsl_eigen_gensymmv_alloc()

Allocate a workspace for solving the generalized symmetric-definite
eigenvalue problem. The size of this workspace is O(4n).

Inputs: n - size of matrices

Return: pointer to workspace
*/

gsl_eigen_gensymmv_workspace *
gsl_eigen_gensymmv_alloc(const size_t n)
{
  gsl_eigen_gensymmv_workspace *w;

  if (n == 0)
    {
      GSL_ERROR_NULL ("matrix dimension must be positive integer",
                      GSL_EINVAL);
    }

  w = (gsl_eigen_gensymmv_workspace *) calloc (1, sizeof (gsl_eigen_gensymmv_workspace));

  if (w == 0)
    {
      GSL_ERROR_NULL ("failed to allocate space for workspace", GSL_ENOMEM);
    }

  w->size = n;

  w->symmv_workspace_p = gsl_eigen_symmv_alloc(n);
  if (!w->symmv_workspace_p)
    {
      gsl_eigen_gensymmv_free(w);
      GSL_ERROR_NULL("failed to allocate space for symmv workspace", GSL_ENOMEM);
    }

  return (w);
} /* gsl_eigen_gensymmv_alloc() */

/*
gsl_eigen_gensymmv_free()
  Free workspace w
*/

void
gsl_eigen_gensymmv_free (gsl_eigen_gensymmv_workspace * w)
{
  RETURN_IF_NULL (w);

  if (w->symmv_workspace_p)
    gsl_eigen_symmv_free(w->symmv_workspace_p);

  free(w);
} /* gsl_eigen_gensymmv_free() */

/*
gsl_eigen_gensymmv()

Solve the generalized symmetric-definite eigenvalue problem

A x = \lambda B x

for the eigenvalues \lambda and eigenvectors x.

Inputs: A    - real symmetric matrix
        B    - real symmetric and positive definite matrix
        eval - where to store eigenvalues
        evec - where to store eigenvectors
        w    - workspace

Return: success or error
*/

int
gsl_eigen_gensymmv (gsl_matrix * A, gsl_matrix * B, gsl_vector * eval,
                    gsl_matrix * evec, gsl_eigen_gensymmv_workspace * w)
{
  const size_t N = A->size1;

  /* check matrix and vector sizes */

  if (N != A->size2)
    {
      GSL_ERROR ("matrix must be square to compute eigenvalues", GSL_ENOTSQR);
    }
  else if ((N != B->size1) || (N != B->size2))
    {
      GSL_ERROR ("B matrix dimensions must match A", GSL_EBADLEN);
    }
  else if (eval->size != N)
    {
      GSL_ERROR ("eigenvalue vector must match matrix size", GSL_EBADLEN);
    }
  else if (evec->size1 != evec->size2)
    {
      GSL_ERROR ("eigenvector matrix must be square", GSL_ENOTSQR);
    }
  else if (evec->size1 != N)
    {
      GSL_ERROR ("eigenvector matrix has wrong size", GSL_EBADLEN);
    }
  else if (w->size != N)
    {
      GSL_ERROR ("matrix size does not match workspace", GSL_EBADLEN);
    }
  else
    {
      int s;

      /* compute Cholesky factorization of B */
      s = gsl_linalg_cholesky_decomp(B);
      if (s != GSL_SUCCESS)
        return s; /* B is not positive definite */

      /* transform to standard symmetric eigenvalue problem */
      gsl_eigen_gensymm_standardize(A, B);

      /* compute eigenvalues and eigenvectors */
      s = gsl_eigen_symmv(A, eval, evec, w->symmv_workspace_p);
      if (s != GSL_SUCCESS)
        return s;

      /* backtransform eigenvectors: evec -> L^{-T} evec */
      gsl_blas_dtrsm(CblasLeft,
                     CblasLower,
                     CblasTrans,
                     CblasNonUnit,
                     1.0,
                     B,
                     evec);

      /* the blas call destroyed the normalization - renormalize */
      gensymmv_normalize_eigenvectors(evec);

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

/********************************************
 *           INTERNAL ROUTINES              *
 ********************************************/

/*
gensymmv_normalize_eigenvectors()
  Normalize eigenvectors so that their Euclidean norm is 1

Inputs: evec - eigenvectors
*/

static void
gensymmv_normalize_eigenvectors(gsl_matrix *evec)
{
  const size_t N = evec->size1;
  size_t i;     /* looping */

  for (i = 0; i < N; ++i)
    {
      gsl_vector_view vi = gsl_matrix_column(evec, i);
      double scale = 1.0 / gsl_blas_dnrm2(&vi.vector);

      gsl_blas_dscal(scale, &vi.vector);
    }
} /* gensymmv_normalize_eigenvectors() */

max_stars_repo_name: ruslankuzmin/julia


text:
#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>

int
main (void)
{
  double data[] = { 1.0  , 1/2.0, 1/3.0, 1/4.0,
                    1/2.0, 1/3.0, 1/4.0, 1/5.0,
                    1/3.0, 1/4.0, 1/5.0, 1/6.0,
                    1/4.0, 1/5.0, 1/6.0, 1/7.0 };

  gsl_matrix_view m
    = gsl_matrix_view_array (data, 4, 4);

  gsl_vector *eval = gsl_vector_alloc (4);
  gsl_matrix *evec = gsl_matrix_alloc (4, 4);

  gsl_eigen_symmv_workspace * w =
    gsl_eigen_symmv_alloc (4);

  gsl_eigen_symmv (&m.matrix, eval, evec, w);

  gsl_eigen_symmv_free (w);

  gsl_eigen_symmv_sort (eval, evec,
                        GSL_EIGEN_SORT_ABS_ASC);

  {
    int i;

    for (i = 0; i < 4; i++)
      {
        double eval_i
           = gsl_vector_get (eval, i);
        gsl_vector_view evec_i
           = gsl_matrix_column (evec, i);

        printf ("eigenvalue = %g\n", eval_i);
        printf ("eigenvector = \n");
        gsl_vector_fprintf (stdout,
                            &evec_i.vector, "%g");
      }
  }

  gsl_vector_free (eval);
  gsl_matrix_free (evec);

  return 0;
}

ext: c

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


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

#ifndef __GSL_VECTOR_UINT_H__
#define __GSL_VECTOR_UINT_H__

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_inline.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_block_uint.h>

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

__BEGIN_DECLS

typedef struct
{
  size_t size;
  size_t stride;
  unsigned int *data;
  gsl_block_uint *block;
  int owner;
}
gsl_vector_uint;

typedef struct
{
  gsl_vector_uint vector;
} _gsl_vector_uint_view;

typedef _gsl_vector_uint_view gsl_vector_uint_view;

typedef struct
{
  gsl_vector_uint vector;
} _gsl_vector_uint_const_view;

typedef const _gsl_vector_uint_const_view gsl_vector_uint_const_view;

/* Allocation */

gsl_vector_uint *gsl_vector_uint_alloc (const size_t n);
gsl_vector_uint *gsl_vector_uint_calloc (const size_t n);

gsl_vector_uint *gsl_vector_uint_alloc_from_block (gsl_block_uint * b,
                                                     const size_t offset,
                                                     const size_t n,
                                                     const size_t stride);

gsl_vector_uint *gsl_vector_uint_alloc_from_vector (gsl_vector_uint * v,
                                                      const size_t offset,
                                                      const size_t n,
                                                      const size_t stride);

void gsl_vector_uint_free (gsl_vector_uint * v);

/* Views */

_gsl_vector_uint_view
gsl_vector_uint_view_array (unsigned int *v, size_t n);

_gsl_vector_uint_view
gsl_vector_uint_view_array_with_stride (unsigned int *base,
                                         size_t stride,
                                         size_t n);

_gsl_vector_uint_const_view
gsl_vector_uint_const_view_array (const unsigned int *v, size_t n);

_gsl_vector_uint_const_view
gsl_vector_uint_const_view_array_with_stride (const unsigned int *base,
                                               size_t stride,
                                               size_t n);

_gsl_vector_uint_view
gsl_vector_uint_subvector (gsl_vector_uint *v,
                            size_t i,
                            size_t n);

_gsl_vector_uint_view
gsl_vector_uint_subvector_with_stride (gsl_vector_uint *v,
                                        size_t i,
                                        size_t stride,
                                        size_t n);

_gsl_vector_uint_const_view
gsl_vector_uint_const_subvector (const gsl_vector_uint *v,
                                  size_t i,
                                  size_t n);

_gsl_vector_uint_const_view
gsl_vector_uint_const_subvector_with_stride (const gsl_vector_uint *v,
                                              size_t i,
                                              size_t stride,
                                              size_t n);

/* Operations */

void gsl_vector_uint_set_zero (gsl_vector_uint * v);
void gsl_vector_uint_set_all (gsl_vector_uint * v, unsigned int x);
int gsl_vector_uint_set_basis (gsl_vector_uint * v, size_t i);

int gsl_vector_uint_fread (FILE * stream, gsl_vector_uint * v);
int gsl_vector_uint_fwrite (FILE * stream, const gsl_vector_uint * v);
int gsl_vector_uint_fscanf (FILE * stream, gsl_vector_uint * v);
int gsl_vector_uint_fprintf (FILE * stream, const gsl_vector_uint * v,
                              const char *format);

int gsl_vector_uint_memcpy (gsl_vector_uint * dest, const gsl_vector_uint * src);

int gsl_vector_uint_reverse (gsl_vector_uint * v);

int gsl_vector_uint_swap (gsl_vector_uint * v, gsl_vector_uint * w);
int gsl_vector_uint_swap_elements (gsl_vector_uint * v, const size_t i, const size_t j);

unsigned int gsl_vector_uint_max (const gsl_vector_uint * v);
unsigned int gsl_vector_uint_min (const gsl_vector_uint * v);
void gsl_vector_uint_minmax (const gsl_vector_uint * v, unsigned int * min_out, unsigned int * max_out);

size_t gsl_vector_uint_max_index (const gsl_vector_uint * v);
size_t gsl_vector_uint_min_index (const gsl_vector_uint * v);
void gsl_vector_uint_minmax_index (const gsl_vector_uint * v, size_t * imin, size_t * imax);

int gsl_vector_uint_add (gsl_vector_uint * a, const gsl_vector_uint * b);
int gsl_vector_uint_sub (gsl_vector_uint * a, const gsl_vector_uint * b);
int gsl_vector_uint_mul (gsl_vector_uint * a, const gsl_vector_uint * b);
int gsl_vector_uint_div (gsl_vector_uint * a, const gsl_vector_uint * b);
int gsl_vector_uint_scale (gsl_vector_uint * a, const double x);
int gsl_vector_uint_add_constant (gsl_vector_uint * a, const double x);

int gsl_vector_uint_equal (const gsl_vector_uint * u,
                            const gsl_vector_uint * v);

int gsl_vector_uint_isnull (const gsl_vector_uint * v);
int gsl_vector_uint_ispos (const gsl_vector_uint * v);
int gsl_vector_uint_isneg (const gsl_vector_uint * v);
int gsl_vector_uint_isnonneg (const gsl_vector_uint * v);

INLINE_DECL unsigned int gsl_vector_uint_get (const gsl_vector_uint * v, const size_t i);
INLINE_DECL void gsl_vector_uint_set (gsl_vector_uint * v, const size_t i, unsigned int x);
INLINE_DECL unsigned int * gsl_vector_uint_ptr (gsl_vector_uint * v, const size_t i);
INLINE_DECL const unsigned int * gsl_vector_uint_const_ptr (const gsl_vector_uint * v, const size_t i);

#ifdef HAVE_INLINE

INLINE_FUN
unsigned int
gsl_vector_uint_get (const gsl_vector_uint * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_VAL ("index out of range", GSL_EINVAL, 0);
    }
#endif
  return v->data[i * v->stride];
}

INLINE_FUN
void
gsl_vector_uint_set (gsl_vector_uint * v, const size_t i, unsigned int x)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_VOID ("index out of range", GSL_EINVAL);
    }
#endif
  v->data[i * v->stride] = x;
}

INLINE_FUN
unsigned int *
gsl_vector_uint_ptr (gsl_vector_uint * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (unsigned int *) (v->data + i * v->stride);
}

INLINE_FUN
const unsigned int *
gsl_vector_uint_const_ptr (const gsl_vector_uint * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (const unsigned int *) (v->data + i * v->stride);
}
#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_VECTOR_UINT_H__ */

ext: h

max_stars_repo_licenses: Net-SNMP Xnet

max_issues_repo_licenses: Net-SNMP Xnet

max_forks_repo_licenses: Net-SNMP Xnet


text:
#ifdef ESL_WITH_GSL
/* interface_gsl.h
 * Easel's interfaces to the GNU Scientific Library
 *
 * SRE, Tue Jul 13 15:36:48 2004
 * SVN $Id: interface_gsl.h 11 2005-01-06 11:44:17Z eddy $
 */
#ifndef ESL_INTERFACE_GSL_INCLUDED
#define ESL_INTERFACE_GSL_INCLUDED

#include <stdlib.h>
#include <easel/easel.h>
#include <easel/dmatrix.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_eigen.h>

extern int esl_GSL_MatrixInversion(ESL_DMATRIX *A, ESL_DMATRIX **ret_Ai);

#endif /*ESL_INTERFACE_GSL_INCLUDED*/
#endif /*ESL_WITH_GSL*/

max_stars_repo_name: juan-rodriguez-rivas/cointerfaces

max_issues_repo_name: juan-rodriguez-rivas/cointerfaces

max_forks_repo_name: juan-rodriguez-rivas/cointerfaces


text:
/*
 *  gsl_binomial_randomdev.h
 *
 *  This file is part of NEST.
 *
 *  Copyright (C) 2004 The NEST Initiative
 *
 *  NEST 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.
 *
 *  NEST 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 NEST.  If not, see <http://www.gnu.org/licenses/>.
 *
 */

// Generated includes:
#include "config.h"

#ifndef GSL_BINOMIAL_RANDOMDEV_H
#define GSL_BINOMIAL_RANDOMDEV_H

// Includes from libnestutil:
#include "lockptr.h"

// Includes from librandom:
#include "gslrandomgen.h"
#include "randomdev.h"
#include "randomgen.h"

// Includes from sli:
#include "dictdatum.h"

#ifdef HAVE_GSL

// External includes:
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>

/*BeginDocumentation
Name: rdevdict::gsl_binomial - GSL binomial random deviate generator
Description:
This function returns a random integer from the binomial distribution,
the number of successes in n independent trials with probability
p. The probability distribution for binomial variates is,

   p(k) = (n! / k!(n-k)!) p^k (1-p)^(n-k)  , 0<=k<=n, n>0

Please note that the RNG used to initialize gsl_binomial has to be
from the GSL (prefixed gsl_ in rngdict)

Parameters:
   p - probability of success in a single trial (double)
   n - number of trials (positive integer)

SeeAlso: CreateRDV, RandomArray, rdevdict
Author: Jochen Martin Eppler
*/

namespace librandom
{

/**
 Class GSL_BinomialRandomDev

 Generates an RNG which returns Binomial(k;p;n)
 distributed random numbers out of an RNG which returns
 binomially distributed random numbers:

    p(k) = (n! / k!(n-k)!) p^k (1-p)^(n-k)  , 0<=k<=n, n<0

 Arguments:
  - pointer to an RNG
  - parameter p (optional, default = 0.5)
  - parameter n (optional, default = 1)

 @see
 http://www.gnu.org/software/gsl/manual/html_node/The-Binomial-Distribution.html
 @ingroup RandomDeviateGenerators
*/

class GSL_BinomialRandomDev : public RandomDev
{
public:
  // accept only lockPTRs for initialization,
  // otherwise creation of a lock ptr would
  // occur as side effect---might be unhealthy
  GSL_BinomialRandomDev( RngPtr, double p_s = 0.5, unsigned int n_s = 1 );
  GSL_BinomialRandomDev( double p_s = 0.5, unsigned int n_s = 1 );

  /**
   * set parameters for p and n
   * @parameters
   * p - success probability for single trial
   * n - number of trials
   */
  void set_p_n( double, unsigned int );
  void set_p( double );       //!<set p
  void set_n( unsigned int ); //!<set n

  /**
   * Import sets of overloaded virtual functions.
   * We need to explicitly include sets of overloaded
   * virtual functions into the current scope.
   * According to the SUN C++ FAQ, this is the correct
   * way of doing things, although all other compilers
   * happily live without.
   */
  using RandomDev::operator();
  using RandomDev::ldev;

  long ldev();               //!< draw integer
  long ldev( RngPtr ) const; //!< draw integer, threaded
  bool
  has_ldev() const
  {
    return true;
  }

  double operator()( RngPtr ) const; //!< return as double, threaded

  //! set distribution parameters from SLI dict
  void set_status( const DictionaryDatum& );

  //! get distribution parameters from SLI dict
  void get_status( DictionaryDatum& ) const;

private:
  double p_;       //!<probability p of binomial distribution
  unsigned int n_; //!<parameter n in binomial distribution

  gsl_rng* rng_;
};

inline double GSL_BinomialRandomDev::operator()( RngPtr rthrd ) const
{
  return static_cast< double >( ldev( rthrd ) );
}
}

#endif

#endif

max_stars_repo_name: OpenHEC/SNN-simulator-on-PYNQcluster

max_issues_repo_name: zlchai/SNN-simulator-on-PYNQcluster

max_forks_repo_name: OpenHEC/SNN-simulator-on-PYNQcluster


text:
#include "pyactpol.h"
#include <gsl/gsl_linalg.h>

/* sync.c
 *
 * scan-synchronous mode fitting and removal.
 */

static int fit_poly(float *x, double *y, int n, int order, int samp,
                    double *fit);
static int fit_sine(double *t, double *y, int n, double f,
                  double *fit, double *coeff);

/* Find the amplitude and phase of the synchronous signal in a set of
 * detectors together with the phase of the azimuth scan.
 *
 * az should already have the mean removed.  ctime should start at 0
 * to minimize rounding errors. nwin should be the number of samples
 * to use for the amplitude fit; i.e. an round number of azimuth scans.
 */

static int mbSyncFit(double* az, double* ctime,
		     float* data, int ndata,
		     int *dets, int ndets,
		     float *sync_amp, float* sync_phase,
		     float *sync_rms,
		     float *phaseAz,
		     float f,
		     int order, int samp, int nT)
{
  int n, nwin;
  double sampRate, thetaAz;

  // Calculate size of data window to do analysis
  sampRate = (double)ndata / (ctime[ndata-1] - ctime[0]);
  nwin = floor( (double)nT / f * sampRate );
  if (nwin > ndata) {
    nwin = ndata;
  }

  double *fit_az = malloc(nwin * sizeof(double));
  double *coeff_az = malloc(2 * sizeof(double));

  // Obtain phase of azimuth scan
  if (fit_sine(ctime, az, nwin, f, fit_az, coeff_az) != 0)
       return 1;
  if (coeff_az[0] > 0) thetaAz = atan(coeff_az[1]/coeff_az[0]);
  else thetaAz = atan(coeff_az[1]/coeff_az[0]) + M_PI;
  *phaseAz = thetaAz;

  free(fit_az);
  free(coeff_az);

  float *tcopy = malloc(nwin * sizeof(float));
  for (int j=0; j<nwin; j++)
       tcopy[j] = ctime[j] - ctime[0];

  // Process all selected detectors
#pragma omp parallel private(n)
  {
    double *y = malloc(nwin * sizeof(double));
    double *fit_1f = malloc(nwin * sizeof(double));
    double *fit_syn = malloc(nwin * sizeof(double));
    double *coeff_syn = malloc(2 * sizeof(double));

#pragma omp for
    for (n = 0; n < ndets; n++) {
      int i;
      double A, theta, rms;

      float *det_data = data + ndata * dets[n];

      // Obtain copy of data
      for (i = 0; i < nwin; i++)
	   y[i] = det_data[i];

      // Fit and subtract polynomial
      fit_poly(tcopy, y, nwin, order, samp, fit_1f);
      for (i = 0; i < nwin; i++)
	y[i] -= fit_1f[i];

      // Fit sinusoidal to obtain amplitude and phase
      fit_sine(ctime, y, nwin, f, fit_syn, coeff_syn);
      A = sqrt(coeff_syn[0]*coeff_syn[0] + coeff_syn[1]*coeff_syn[1]);
      if (A == 0) theta = 0.0;
      else if (coeff_syn[0] > 0) theta = atan(coeff_syn[1]/coeff_syn[0]);
      else theta = atan(coeff_syn[1]/coeff_syn[0]) + M_PI;
      sync_amp[n] = A;
      sync_phase[n] = theta;

      // Find fit RMS error
      rms = 0.0;
      for (i = 0; i < nwin; i++)
	   rms += (y[i] - fit_syn[i])*(y[i] - fit_syn[i]);
      sync_rms[n] = sqrt(rms/nwin);
    } /* pragma omp for */

    free(y);
    free(fit_1f);
    free(fit_syn);
    free(coeff_syn);
  } /* pragma omp parallel */

  return 0;
}

/// Fit sinusoidal of a given frequency to data.
/// \param t      Independent variable (time).
/// \param y      Dependent variable (data).
/// \param n      Number of elements in data.
/// \param f      Frequency of the sinusoid.
/// \param fit    Vector of fitted data (sinusoid).
/// \param coeff  Coefficients of sinusoid (A*cos(th) + B*sin(th)).
static int fit_sine(double *t, double *y, int n, double f,
	          double *fit, double *coeff) {
  int k;

  double temp_c, temp_s;
  double cos2, sin2, sincos, ycos, ysin;

  cos2 = sin2 = sincos = ycos = ysin = 0.0;
  for (k = 0; k < n; k++) {
    temp_c = cos(2*M_PI*f*t[k]);
    temp_s = sin(2*M_PI*f*t[k]);
    cos2 += temp_c*temp_c;
    sin2 += temp_s*temp_s;
    sincos += temp_s*temp_c;
    ycos += y[k]*temp_c;
    ysin += y[k]*temp_s;
  }

  coeff[0] = (sin2*ycos - sincos*ysin) / (cos2*sin2 - sincos*sincos);
  coeff[1] = (cos2*ysin - sincos*ycos) / (cos2*sin2 - sincos*sincos);

  // Evaluate result
  for (k = 0; k < n; k++) {
    fit[k] = cos(2*M_PI*f*t[k])*coeff[0] + sin(2*M_PI*f*t[k])*coeff[1];
  }

  return 0;
}

/// Fit polynomial to data to remove common mode.
/// \param x       Independent variable.
/// \param y       Dependent variable (data).
/// \param n       Number of elements in data.
/// \param order   Order of polinomial to fit.
/// \param samp    Number of samples to take from data to do the fit.
/// \param fit     vector with evaluated polynomial

static int fit_poly(float *x, double *y, int n, int order, int samp,
                    double *fit) {
  int i, j, k;
  int N;
  double *coeff;
  double **A, *pows;
  double temp;

  gsl_vector *g_Ab = gsl_vector_calloc(order+1);
  gsl_matrix *g_AA = gsl_matrix_alloc(order+1,order+1);
  gsl_vector *g_x = gsl_vector_alloc(order+1);
  gsl_permutation *g_p = gsl_permutation_alloc(order+1);
  int signum = 0;

  int status = 1; // failure.

  if (order <= 0) {
       print_error("Invalid 'order=%i' in fit_poly\n", order);
       return 1;
  }

  // Initializa Matrices and Vectors
  coeff = malloc((order+1) * sizeof(double));
  N = n/samp;
  if (N < order) {
       print_error("Fitting 'order=%i' with 'N=%i' data points in fit_poly\n", order, N);
       return 1;
  }
  A = (double **)malloc((order+1) * sizeof(double *));
  A[0] = (double *)malloc(N * (order+1) * sizeof(double));
  for (i = 0; i < order+1; i++){
    A[i] = A[0] + N * i;
    for( j = 0; j < N; j++ ) A[i][j] = 0.0;
  }

  pows = (double *)malloc((2*order+1) * sizeof(double));
  for (i = 0; i < 2*order+1; i++)
    pows[i] = 0.0;

  // Generate A matrix
  for (i = 0; i < order+1; i++)
    for (k = 0; k < N; k++) {
      A[i][k] = 1.0;
      for (j = 0; j < i; j++)
        A[i][k] *= x[k*samp + samp/2];
    }

  // Generate AA as transpose(A)*A
  for (i = 2*order; i > 0; i -= 2) {
    for (k = 0; k < N; k++) {
      pows[i] += A[i/2][k]*A[i/2][k];
      pows[i-1] += A[i/2][k]*A[i/2-1][k];
    }
  }
  for (k = 0; k < N; k++)
    pows[0]++;

  for (i = 0; i < order+1; i++)
    for (j = 0; j < order+1; j++)
      g_AA->data[i*g_AA->tda+j] = pows[i+j];

/*  for (i = 0; i < order+1; i++) {
	for (j = 0; j <= i; j++) {
	  for (k = 0; k < N; k++)
		AA[i][j] += A[j][k]*A[i][k];
	}
  }*/

  // Generate Ab as transpose(A)*y
  for (i = 0; i < order+1; i++) {
    for (k = 0; k < N; k++) {
      temp = 0.0;
      for (j = 0; j < samp; j++)
        temp += y[k*samp + j];
      g_Ab->data[i] += A[i][k] * temp / (float)samp;
    }
  }

  // Solve system
  if ((status = gsl_linalg_LU_decomp(g_AA, g_p, &signum))!=0)
      goto exit_now;
  if ((status = gsl_linalg_LU_solve(g_AA, g_p, g_Ab, g_x))!=0)
      goto exit_now;

  // Copy out...
  for (i = 0; i < order+1; i++)
      coeff[i] = g_x->data[i];

  // Evaluate result
  for (k = 0; k < n; k++) {
    fit[k] = 0.0;
    for (i = 0; i < order+1; i++) {
      temp = 1.0;
      for (j = 0; j < i; j++)
        temp *= x[k];
      fit[k] += temp * coeff[i];
    }
  }

/*
  for (k = 0; k < n; k++)
	fit[k] = 0.0;
  for (i = 0; i < order+1; i++)
	for (k = 0; k < n; k++)
	  fit[k] += A[i][k] * Ab[i];
*/

  status = 0;

exit_now:
  free(A[0]);

  free(A);
  free(pows);
  free(coeff);

  gsl_vector_free(g_Ab);
  gsl_vector_free(g_x);
  gsl_matrix_free(g_AA);
  gsl_permutation_free(g_p);

  return status;
}

PyDoc_STRVAR(get_sync_amps__doc__,
	     "get_sync_amps(data, dets, az, ctime)\n"
	     "\n"
	     "The arrays must be C-ordered with dimensions like:\n"
	     "     data [     *,n_data]    (float)\n"
	     "     dets [n_det]    	   (int)\n"
             "     az   [n_data]   	   (double)\n"
	     "     ctime[n_data]   	   (double)\n"
	     "\n"
	     "Returns (az_phase, amps_cos, amps_sin)."
     );

static PyObject *get_sync_amps(PyObject *self, PyObject *args)
{
    PyArrayObject *data_array;
    PyArrayObject *dets_array;
    PyArrayObject *az_array;
    PyArrayObject *ctime_array;
    double scan_freq;
    int order, samp, nT;

    if (!PyArg_ParseTuple(args, "O!O!O!O!diii",
                          &PyArray_Type, &data_array,
                          &PyArray_Type, &dets_array,
                          &PyArray_Type, &az_array,
                          &PyArray_Type, &ctime_array,
			  &scan_freq,
			  &order,
			  &samp,
			  &nT
            ))
        po_raise("invalid arguments.");

    // Types and ordering
    ASSERT_CARRAY_TYPE_NDIM(data_array, NPY_FLOAT32, 2);
    ASSERT_CARRAY_TYPE_NDIM(dets_array, NPY_INT32, 1);
    ASSERT_CARRAY_TYPE_NDIM(az_array, NPY_FLOAT64, 1);
    ASSERT_CARRAY_TYPE_NDIM(ctime_array, NPY_FLOAT64, 1);

    int ndata = PyArray_DIMS(data_array)[1];
    int ndet = PyArray_DIMS(dets_array)[0];

    po_assert(PyArray_DIMS(az_array)[0] == ndata);
    po_assert(PyArray_DIMS(ctime_array)[0] == ndata);

    float *data = PyArray_DATA(data_array);
    double *az = PyArray_DATA(az_array);
    double *ctime = PyArray_DATA(ctime_array);
    int *dets = PyArray_DATA(dets_array);

    // We've been so thorough, it would be a shame to segfault now.
    for (int i=0; i<ndet; i++)
	 po_assert(dets[i] >= 0 &&
		   dets[i] < PyArray_DIMS(data_array)[0]);

    // And, places for the results.
    npy_intp ndet_ = ndet;
    PyArrayObject *amp_array = (PyArrayObject*)
	 PyArray_SimpleNew(1, &ndet_, NPY_FLOAT32);
    PyArrayObject *phase_array = (PyArrayObject*)
	 PyArray_SimpleNew(1, &ndet_, NPY_FLOAT32);
    PyArrayObject *rms_array = (PyArrayObject*)
	 PyArray_SimpleNew(1, &ndet_, NPY_FLOAT32);

    float phaseAz = -1.;
    mbSyncFit(az, ctime,
	      data, ndata,
	      dets, ndet,
	      PyArray_DATA(amp_array),
	      PyArray_DATA(phase_array),
	      PyArray_DATA(rms_array),
	      &phaseAz,
	      scan_freq, order, samp, nT);

    return Py_BuildValue("NNNf",
			 amp_array,
			 phase_array,
			 rms_array,
			 phaseAz);
}

PyMethodDef pyactpol_sync_methods[] = {
    {"get_sync_amps", get_sync_amps, METH_VARARGS,
    get_sync_amps__doc__},
    {NULL, NULL, 0, NULL}        /* Sentinel */
};


text:
/*****************************************************************\
           __
          / /
		 / /                     __  __
		/ /______    _______    / / / / ________   __       __
	   / ______  \  /_____  \  / / / / / _____  | / /      / /
	  / /      | / _______| / / / / / / /____/ / / /      / /
	 / /      / / / _____  / / / / / / _______/ / /      / /
	/ /      / / / /____/ / / / / / / |______  / |______/ /
   /_/      /_/ |________/ / / / /  \_______/  \_______  /
                          /_/ /_/                     / /
			                                         / /
		       High Level Game Framework            /_/

  ---------------------------------------------------------------

  Copyright (c) 2007-2011 - Rodrigo Braz Monteiro.
  This file is subject to the terms of halley_license.txt.

\*****************************************************************/

#pragma once

#include <iostream>
#include <iomanip>
#include <sstream>
#include <cmath>
#include "halley/text/halleystring.h"
#include <gsl/gsl_assert>
#include <cstdint>

namespace Halley {
	// This whole class is TERRIBLE
	// There should be a type distinction between linear and gamma space, with conversions,
	// and the default Colour type should be Colour4c

	template <typename T>
	constexpr T colMinValue()
	{
		return 0;
	}

	template <typename T>
	constexpr T colMaxValue()
	{
		return 255;
	}

	template <>
	constexpr float colMaxValue<float>()
	{
		return 1.0f;
	}

	template <typename T, typename U>
	constexpr U convertColour(T x)
	{
		return U(float(x) * colMaxValue<U>() / colMaxValue<T>());
	}

	template <>
	constexpr uint8_t convertColour(int x)
	{
		return static_cast<uint8_t>(x);
	}

	template <typename T>
	class Colour4 {
	public:
		T r = colMinValue<T>();
		T g = colMinValue<T>();
		T b = colMinValue<T>();
		T a = colMaxValue<T>();

		constexpr Colour4() = default;

		constexpr Colour4(T luma)
			: r(luma)
			, g(luma)
			, b(luma)
		{}

		constexpr Colour4(T r, T g, T b, T a=colMaxValue<T>())
			: r(r)
			, g(g)
			, b(b)
			, a(a)
		{}

		Colour4(const String& str)
		{
			*this = fromString(str);
		}

		template <typename U>
		Colour4(const Colour4<U> c)
		{
			r = convertColour<U, T>(c.r);
			g = convertColour<U, T>(c.g);
			b = convertColour<U, T>(c.b);
			a = convertColour<U, T>(c.a);
		}

		String toString() const
		{
			std::stringstream ss;
			ss << "#" << std::hex;
			writeByte(ss, r);
			writeByte(ss, g);
			writeByte(ss, b);
			if (byteRep(a) != 255) writeByte(ss, a);
			ss.flush();
			return ss.str();
		}

		static Colour4 fromString(String str)
		{
			Colour4 col;
			size_t len = str.length();
			if (len >= 1 && str[0] == '#') {
				if (len >= 3) {
					col.r = parseHex(str.substr(1, 2));
				}
				if (len >= 5) {
					col.g = parseHex(str.substr(3, 2));
				}
				if (len >= 7) {
					col.b = parseHex(str.substr(5, 2));
				}
				if (len >= 9) {
					col.a = parseHex(str.substr(7, 2));
				}
			}
			return col;
		}

		static Colour4 fromHSV(float h, float s, float v)
		{
			float r = 0;
			float g = 0;
			float b = 0;
			if (s == 0) {
				r = v * 255;
				if (r < 0) r = 0;
				if (r > 255) r = 255;
				g = b = r;
			} else {
				h = float(std::fmod(h, 1.0f));
				int hi = int(h * 6.0f);
				float f = h*6 - hi;
				float p = v*(1-s);
				float q = v*(1-f*s);
				float t = v*(1-(1-f)*s);

				switch (hi) {
				case 0:	r = v; g = t; b = p; break;
				case 1:	r = q; g = v; b = p; break;
				case 2:	r = p; g = v; b = t; break;
				case 3:	r = p; g = q; b = v; break;
				case 4:	r = t; g = p; b = v; break;
				case 5:	r = v; g = p; b = q; break;
				//default: throw Exception("Invalid value!");
				}
			}

			return Colour4<T>(Colour4<float>(r, g, b));
		}

		constexpr Colour4 multiplyLuma(float t) const
		{
			return Colour4(r*t, g*t, b*t, a);
		}

		constexpr Colour4 inverseMultiplyLuma(float t) const
		{
			return Colour4(1.0f - ((1.0f - r) * t), 1.0f - ((1.0f - g) * t), 1.0f - ((1.0f - b) * t), a);
		}

		constexpr Colour4 multiplyAlpha(float t) const
		{
			return Colour4(r, g, b, a * t);
		}

		constexpr Colour4 operator+(const Colour4& c) const
		{
			return Colour4(r+c.r, g+c.g, b+c.b, a+c.a);
		}

		constexpr Colour4 operator*(const Colour4& c) const
		{
			return Colour4(r * c.r, g * c.g, b * c.b, a * c.a);
		}

		constexpr Colour4 operator*(float t) const
		{
			return Colour4(r*t, g*t, b*t, a*t);
		}

		constexpr bool operator==(const Colour4& c) const
		{
			return r == c.r && g == c.g && b == c.b && a == c.a;
		}

		constexpr bool operator!=(const Colour4& c) const
		{
			return r != c.r || g != c.g || b != c.b || a != c.a;
		}

	private:
		uint8_t byteRep(T v) const
		{
			return convertColour<T, uint8_t>(v);
		}

		static T parseHex(String str)
		{
			std::stringstream ss(str);
			int value;
			ss << std::hex;
			ss >> value;
			return convertColour<uint8_t, T>(static_cast<uint8_t>(value));
		}

		template <typename U>
		void writeByte(U& os, T component) const
		{
			unsigned int value = byteRep(component);
			if (value < 16) os << '0';
			os << value;
		}
	};

	typedef Colour4<uint8_t> Colour4c;
	typedef Colour4<float> Colour4f;
	typedef Colour4f Colour;
}

max_stars_repo_name: chidddy/halley

max_issues_repo_name: chidddy/halley

max_forks_repo_name: chidddy/halley


text:
/* randist/rayleigh.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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

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

/* The Rayleigh distribution has the form

   p(x) dx = (x / sigma^2) exp(-x^2/(2 sigma^2)) dx

   for x = 0 ... +infty */

double
gsl_ran_rayleigh (const gsl_rng * r, const double sigma)
{
  double u = gsl_rng_uniform_pos (r);

  return sigma * sqrt(-2.0 * log (u));
}

double
gsl_ran_rayleigh_pdf (const double x, const double sigma)
{
  if (x < 0)
    {
      return 0 ;
    }
  else
    {
      double u = x / sigma ;
      double p = (u / sigma) * exp(-u * u / 2.0) ;

      return p;
    }
}

/* The Rayleigh tail distribution has the form

   p(x) dx = (x / sigma^2) exp((a^2 - x^2)/(2 sigma^2)) dx

   for x = a ... +infty */

double
gsl_ran_rayleigh_tail (const gsl_rng * r, const double a, const double sigma)
{
  double u = gsl_rng_uniform_pos (r);

  return sqrt(a * a - 2.0 * sigma * sigma * log (u));
}

double
gsl_ran_rayleigh_tail_pdf (const double x, const double a, const double sigma)
{
  if (x < a)
    {
      return 0 ;
    }
  else
    {
      double u = x / sigma ;
      double v = a / sigma ;

      double p = (u / sigma) * exp((v + u) * (v - u) / 2.0) ;

      return p;
    }
}


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "compearth.h"
/*
#ifdef COMPEARTH_USE_MKL
#include <mkl_cblas.h>
#else
#include <cblas.h>
#endif
*/

/*!
 * @brief Converts a moment tensor, M, from input system defined by i1in
 *        to moment tensor, Mout, to output system defind by i2in.
 *
 * @param[in] nmt      Number of moment tensors.
 * @param[in] i1in     Coordinate system for M: \n
 *                       = CE_USE (1) -> up, south, east \n
 *                       = CE_NED (2) -> north, east, down \n
 *                       = CE_NWU (3) -> north west, up \n
 *                       = CE_ENU (4) -> east, north, up \n
 *                       = CE_SEU (5) -> south, east, up \n
 * @param[in] i2in     Coordinate system for Mout: \n
 *                       = CE_USE (1) -> up, south, east \n
 *                       = CE_NED (2) -> north, east, down \n
 *                       = CE_NWU (3) -> north west, up \n
 *                       = CE_ENU (4) -> east, north, up \n
 *                       = CE_SEU (5) -> south, east, up \n
 * @param[in] M        Input moment tensor in system i1in.  This is an
 *                     an array of dimension [6 x nmt].
 *                     The C indices {0,1,2,3,4,5} correspond to matrix
 *                     indices: {11, 22, 33, 12, 13, 23}.
 *
 * @param[out] Mout    Corresponding moment tensor now in system i2in.
 *                     This is an array of dimension [6 x nmt].
 *                     the C indices {0,1,2,3,4,5} correspond to matrix
 *                     indices: {11, 22, 33, 12, 13, 23}
 *
 * @result 0 indicates success.
 *
 * @date 2016 - Ben Baker converted Carl Tape's convert_MT.m to C
 *
 * @copyright MIT
 *
 */
int compearth_convertMT(const int nmt,
                        const enum compearthCoordSystem_enum i1in,
                        const enum compearthCoordSystem_enum i2in,
                        const double *__restrict__ M,
                        double *__restrict__ Mout)
{
    int i, i1, i2;
    // Check the inputs to avoid seg faults
    if (nmt < 1 || M == NULL || Mout == NULL)
    {
        if (nmt < 1){fprintf(stderr, "%s: No moment tensors\n", __func__);}
        if (M == NULL){fprintf(stderr, "%s: Error M is NULL\n", __func__);}
        if (Mout == NULL)
        {
            fprintf(stderr, "%s: Error Mout is NULL\n", __func__);
        }
        return -1;
    }
    for (i=0; i<6*nmt; i++){Mout[i] = 0.0;}
    // Quick checks
    i1 = (int) i1in;
    i2 = (int) i2in;
    if (i1 < 1 || i1 > 5)
    {
        fprintf(stderr, "%s: Error unknown input coordinate system %d\n",
                __func__, i1);
        return -1;
    }
    if (i2 < 1 || i2 > 5)
    {
        fprintf(stderr, "%s: Error unkonwn output coordinate system %d\n",
                __func__, i2);
        return -1;
    }
    // Base case
    if (i1 == i2)
    {
        memcpy(Mout, M, (size_t) nmt*6*sizeof(double));
        //cblas_dcopy(6*nmt, M, 1, Mout, 1);
        return 0;
    }
    // Convert
    if (i1 == 1)
    {
        // up-south-east (GCMT) to north-east-down (AkiRichards)
        // (AR, 1980, p. 118)
        if (i2 == 2)
        {
            for (i=0; i<nmt; i++)
            {
                Mout[6*i+0] = M[6*i+1]; //  tt -> xx
                Mout[6*i+1] = M[6*i+2]; //  pp -> yy
                Mout[6*i+2] = M[6*i+0]; //  rr -> zz
                Mout[6*i+3] =-M[6*i+5]; // -tp -> xy
                Mout[6*i+4] = M[6*i+3]; //  rt -> xz
                Mout[6*i+5] =-M[6*i+4]; // -rp -> yz
            }
       }
       // up-south-east (GCMT) to north-west-up
       else if (i2 == 3)
       {
            for (i=0; i<nmt; i++)
            {
                Mout[6*i+0] = M[6*i+1];
                Mout[6*i+1] = M[6*i+2];
                Mout[6*i+2] = M[6*i+0];
                Mout[6*i+3] = M[6*i+5];
                Mout[6*i+4] =-M[6*i+3];
                Mout[6*i+5] =-M[6*i+4];
            }
       }
       // up-south-east (GCMT) to east-north-up
       else if (i2 == 4)
       {
            for (i=0; i<nmt; i++)
            {
                Mout[6*i+0] = M[6*i+2];
                Mout[6*i+1] = M[6*i+1];
                Mout[6*i+2] = M[6*i+0];
                Mout[6*i+3] =-M[6*i+5];
                Mout[6*i+4] = M[6*i+4];
                Mout[6*i+5] =-M[6*i+3];
            }
       }
       // up-south-east (GCMT) to south-east-up
       else if (i2 == 5)
       {
            for (i=0; i<nmt; i++)
            {
                Mout[6*i+0] = M[6*i+1];
                Mout[6*i+1] = M[6*i+2];
                Mout[6*i+2] = M[6*i+0];
                Mout[6*i+3] = M[6*i+5];
                Mout[6*i+4] = M[6*i+3];
                Mout[6*i+5] = M[6*i+4];
            }
        }
    }
    else if (i1 == 2)
    {
        // north-east-down (AkiRichards) to up-south-east (GCMT)
        // (AR, 1980, p. 118)
        if (i2 == 1)
        {
            for (i=0; i<nmt; i++)
            {
                Mout[6*i+0] = M[6*i+2]; //  zz -> rr
                Mout[6*i+1] = M[6*i+0]; //  xx -> tt
                Mout[6*i+2] = M[6*i+1]; //  yy -> pp
                Mout[6*i+3] = M[6*i+4]; //  xz -> rt
                Mout[6*i+4] =-M[6*i+5]; // -yz -> rp
                Mout[6*i+5] =-M[6*i+3]; // -xy -> tp
            }
        }
        // north-east-down (AkiRichards) to north-west-up
        else if (i2 == 3)
        {
            for (i=0; i<nmt; i++)
            {
                Mout[6*i+0] = M[6*i+0];
                Mout[6*i+1] = M[6*i+1];
                Mout[6*i+2] = M[6*i+2];
                Mout[6*i+3] =-M[6*i+3];
                Mout[6*i+4] =-M[6*i+4];
                Mout[6*i+5] = M[6*i+5];
            }
        }
        // north-east-down (AkiRichards) to east-north-up
        else if (i2 == 4)
        {
            for (i=0; i<nmt; i++)
            {
                Mout[6*i+0] = M[6*i+1];
                Mout[6*i+1] = M[6*i+0];
                Mout[6*i+2] = M[6*i+2];
                Mout[6*i+3] = M[6*i+3];
                Mout[6*i+4] =-M[6*i+5];
                Mout[6*i+5] =-M[6*i+4];
            }
        }
        // north-east-down (AkiRichards) to south-east-up
        else if (i2 == 5)
        {
            for (i=0; i<nmt; i++)
            {
                Mout[6*i+0] = M[6*i+0];
                Mout[6*i+1] = M[6*i+1];
                Mout[6*i+2] = M[6*i+2];
                Mout[6*i+3] =-M[6*i+3];
                Mout[6*i+4] = M[6*i+4];
                Mout[6*i+5] =-M[6*i+5];
            }
        }
    }
    else if (i1 == 3)
    {
        // north-west-up to up-south-east (GCMT)
        if (i2 == 1)
        {
            for (i=0; i<nmt; i++)
            {
                Mout[6*i+0] = M[6*i+2];
                Mout[6*i+1] = M[6*i+0];
                Mout[6*i+2] = M[6*i+1];
                Mout[6*i+3] =-M[6*i+4];
                Mout[6*i+4] =-M[6*i+5];
                Mout[6*i+5] = M[6*i+3];
            }
        }
        // north-west-up to north-east-down (AkiRichards)
        else if (i2 == 2)
        {
            for (i=0; i<nmt; i++)
            {
                Mout[6*i+0] = M[6*i+0];
                Mout[6*i+1] = M[6*i+1];
                Mout[6*i+2] = M[6*i+2];
                Mout[6*i+3] =-M[6*i+3];
                Mout[6*i+4] =-M[6*i+4];
                Mout[6*i+5] = M[6*i+5];
            }
        }
        // north-west-up to east-north-up
        else if (i2 == 4)
        {
            for (i=0; i<nmt; i++)
            {
                Mout[6*i+0] = M[6*i+1];
                Mout[6*i+1] = M[6*i+0];
                Mout[6*i+2] = M[6*i+2];
                Mout[6*i+3] =-M[6*i+3];
                Mout[6*i+4] =-M[6*i+5];
                Mout[6*i+5] = M[6*i+4];
            }
        }
        // north-west-up to south-east-up
        else if (i2 == 5)
        {
            for (i=0; i<nmt; i++)
            {
                Mout[6*i+0] = M[6*i+0];
                Mout[6*i+1] = M[6*i+1];
                Mout[6*i+2] = M[6*i+2];
                Mout[6*i+3] = M[6*i+3];
                Mout[6*i+4] =-M[6*i+4];
                Mout[6*i+5] =-M[6*i+5];
            }
        }
    }
    else if (i1 == 4)
    {
        // east-north-up to up-south-east (GCMT)
        if (i2 == 1)
        {
            for (i=0; i<nmt; i++)
            {
                Mout[6*i+0] = M[6*i+2];
                Mout[6*i+1] = M[6*i+1];
                Mout[6*i+2] = M[6*i+0];
                Mout[6*i+3] =-M[6*i+5];
                Mout[6*i+4] = M[6*i+4];
                Mout[6*i+5] =-M[6*i+3];
            }
        }
        // east-north-up to north-east-down (AkiRichards)
        else if (i2 == 2)
        {
            for (i=0; i<nmt; i++)
            {
                Mout[6*i+0] = M[6*i+1];
                Mout[6*i+1] = M[6*i+0];
                Mout[6*i+2] = M[6*i+2];
                Mout[6*i+3] = M[6*i+3];
                Mout[6*i+4] =-M[6*i+5];
                Mout[6*i+5] =-M[6*i+4];
            }
        }
        // east-north-up to north-west-up
        else if (i2 == 3)
        {
            for (i=0; i<nmt; i++)
            {
                Mout[6*i+0] = M[6*i+1];
                Mout[6*i+1] = M[6*i+0];
                Mout[6*i+2] = M[6*i+2];
                Mout[6*i+3] =-M[6*i+3];
                Mout[6*i+4] = M[6*i+5];
                Mout[6*i+5] =-M[6*i+4];
            }
        }
        // east-north-up to south-east-up
        else if (i2 == 5)
        {
            for (i=0; i<nmt; i++)
            {
                Mout[6*i+0] = M[6*i+1];
                Mout[6*i+1] = M[6*i+0];
                Mout[6*i+2] = M[6*i+2];
                Mout[6*i+3] =-M[6*i+3];
                Mout[6*i+4] =-M[6*i+5];
                Mout[6*i+5] = M[6*i+4];
            }
        }
    }
    else if (i1 == 5)
    {
        // south-east-up to up-south-east (GCMT)
        if (i2 == 1)
        {
            for (i=0; i<nmt; i++)
            {
                Mout[6*i+0] = M[6*i+2];
                Mout[6*i+1] = M[6*i+0];
                Mout[6*i+2] = M[6*i+1];
                Mout[6*i+3] = M[6*i+4];
                Mout[6*i+4] = M[6*i+5];
                Mout[6*i+5] = M[6*i+3];
            }
        }
        // south-east-up to north-east-down (AkiRichards)
        else if (i2 == 2)
        {
            for (i=0; i<nmt; i++)
            {
                Mout[6*i+0] = M[6*i+0];
                Mout[6*i+1] = M[6*i+1];
                Mout[6*i+2] = M[6*i+2];
                Mout[6*i+3] =-M[6*i+3];
                Mout[6*i+4] = M[6*i+4];
                Mout[6*i+5] =-M[6*i+5];
            }
        }
        // south-east-up to north-west-up
        else if (i2 == 3)
        {
            for (i=0; i<nmt; i++)
            {
                Mout[6*i+0] = M[6*i+0];
                Mout[6*i+1] = M[6*i+1];
                Mout[6*i+2] = M[6*i+2];
                Mout[6*i+3] = M[6*i+3];
                Mout[6*i+4] =-M[6*i+4];
                Mout[6*i+5] =-M[6*i+5];
            }
        }
        // south-east-up to east-north-up
        else if (i2 == 4)
        {
            for (i=0; i<nmt; i++)
            {
                Mout[6*i+0] = M[6*i+1];
                Mout[6*i+1] = M[6*i+0];
                Mout[6*i+2] = M[6*i+2];
                Mout[6*i+3] =-M[6*i+3];
                Mout[6*i+4] = M[6*i+5];
                Mout[6*i+5] =-M[6*i+4];
            }
        }
    }
    return 0;
}

max_stars_repo_name: OUCyf/mtbeach

max_issues_repo_name: carltape/mtbeach

max_forks_repo_name: carltape/mtbeach


text:
#ifndef _AVA_COMMON_SUPPORT_STD_SPAN_H_
#define _AVA_COMMON_SUPPORT_STD_SPAN_H_
#pragma once

#include <gsl/span>
#include <gsl/span_ext>

namespace std {
using gsl::span;
}  // namespace std

#define EMPTY_CHAR_SPAN std::span<const char>()

#define STRING_AS_SPAN(STR_VAR) std::span<const char>((STR_VAR).data(), (STR_VAR).length())

#define VECTOR_AS_SPAN(VEC_VAR) std::span<const decltype(VEC_VAR)::value_type>((VEC_VAR).data(), (VEC_VAR).size())

#define VECTOR_AS_CHAR_SPAN(VEC_VAR)                                      \
  std::span<const char>(reinterpret_cast<const char *>((VEC_VAR).data()), \
                        sizeof(decltype(VEC_VAR)::value_type) * (VEC_VAR).size())

#endif  // _AVA_COMMON_SUPPORT_STD_SPAN_H_

max_stars_repo_name: anjuna-security/ava

max_issues_repo_name: anjuna-security/ava

max_forks_repo_name: anjuna-security/ava


text:
#ifndef OPENMC_TALLIES_FILTER_SURFACE_H
#define OPENMC_TALLIES_FILTER_SURFACE_H

#include <cstdint>
#include <unordered_map>

#include <gsl/gsl>

#include "openmc/tallies/filter.h"
#include "openmc/vector.h"

namespace openmc {

//==============================================================================
//! Specifies which surface particles are crossing
//==============================================================================

class SurfaceFilter : public Filter {
public:
  //----------------------------------------------------------------------------
  // Constructors, destructors

  ~SurfaceFilter() = default;

  //----------------------------------------------------------------------------
  // Methods

  std::string type() const override { return "surface"; }

  void from_xml(pugi::xml_node node) override;

  void get_all_bins(const Particle& p, TallyEstimator estimator,
    FilterMatch& match) const override;

  void to_statepoint(hid_t filter_group) const override;

  std::string text_label(int bin) const override;

  //----------------------------------------------------------------------------
  // Accessors

  void set_surfaces(gsl::span<int32_t> surfaces);

private:
  //----------------------------------------------------------------------------
  // Data members

  //! The indices of the surfaces binned by this filter.
  vector<int32_t> surfaces_;

  //! A map from surface indices to filter bin indices.
  std::unordered_map<int32_t, int> map_;
};

} // namespace openmc
#endif // OPENMC_TALLIES_FILTER_SURFACE_H


text:
#ifndef CWANNIER_BAND_ENERGY
#define CWANNIER_BAND_ENERGY

#include <stdbool.h>
#include <stdio.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_sort_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_eigen.h>
#include "HTightBinding.h"
#include "ctetra/sum.h"

double BandEnergy(double *E_Fermi, HTightBinding *Hrs, gsl_matrix *R, double num_electrons, int na, int nb, int nc, bool use_cache);

#endif // CWANNIER_BAND_ENERGY

max_stars_repo_name: tflovorn/cwannier

max_issues_repo_name: tflovorn/cwannier

max_forks_repo_name: tflovorn/cwannier


text:
/*
 * Copyright (c) 2018-2021 Aleksas Mazeliauskas, Stefan Floerchinger,
 *                    Eduardo Grossi, and Derek Teaney
 * All rights reserved.
 *
 * FastReso is distributed under MIT license;
 * see the LICENSE file that should be present in the root
 * of the source distribution, or alternately available at:
 * https://github.com/amazeliauskas/FastReso/
 */
#ifndef FASTFO_qag_params_h
#define FASTFO_qag_params_h
#include <gsl/gsl_integration.h>
namespace qag_params{
  const double fEpsRel = 1e-6;
  const double fEpsAbs = 1e-10;
  const int fKey = GSL_INTEG_GAUSS51;
  const int fLimit = 1000;
}
#endif

max_stars_repo_name: amazeliauskas/FastReso

max_issues_repo_name: amazeliauskas/FastReso

max_forks_repo_name: amazeliauskas/FastReso


text:
#pragma once

#include <gsl/span>

/**
* This header is mainly for shorthands when dealing with `gsl::span`s
* see the `gsl` or c++ core guidelines for more info on what a `gsl::span` is/does.
*/

namespace riner {

    template<class T, std::ptrdiff_t Extent = gsl::dynamic_extent>
    using span = gsl::span<T, Extent>;

    template<std::ptrdiff_t Extent = gsl::dynamic_extent>
    using ByteSpan = span<uint8_t, Extent>;

    template<std::ptrdiff_t Extent = gsl::dynamic_extent>
    using cByteSpan = span<const uint8_t, Extent>;
}

max_stars_repo_name: DavHau/Riner

max_issues_repo_name: DavHau/Riner

max_forks_repo_name: DavHau/Riner


text:
#ifndef _PC_PATCH_H
#define _PC_PATCH_H
#include <petsc.h>
PETSC_EXTERN PetscErrorCode PCPatchInitializePackage(void);
PETSC_EXTERN PetscErrorCode PCCreate_PATCH(PC);
PETSC_EXTERN PetscErrorCode PCPatchSetCellNumbering(PC, PetscSection);
PETSC_EXTERN PetscErrorCode PCPatchSetDiscretisationInfo(PC, PetscInt, DM *, PetscInt *, PetscInt *, const PetscInt **, const PetscInt *, PetscInt, const PetscInt *, PetscInt, const PetscInt *);
PETSC_EXTERN PetscErrorCode PCPatchSetComputeOperator(PC, PetscErrorCode (*)(PC,Mat,PetscInt,const PetscInt *,PetscInt,const PetscInt *,void *),
                                                      void *);

typedef enum {PC_PATCH_STAR, PC_PATCH_VANKA, PC_PATCH_USER, PC_PATCH_PYTHON} PCPatchConstructType;

PETSC_EXTERN const char *const PCPatchConstructTypes[];
#endif

max_stars_repo_name: wence-/ssc

max_issues_repo_name: wence-/ssc

max_forks_repo_name: wence-/ssc


text:
/* sys/pow_int.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */
#include <config.h>
#include <math.h>
#include <gsl/gsl_pow_int.h>

#ifndef HIDE_INLINE_STATIC
double gsl_pow_2(const double x) { return x*x;   }
double gsl_pow_3(const double x) { return x*x*x; }
double gsl_pow_4(const double x) { double x2 = x*x;   return x2*x2;    }
double gsl_pow_5(const double x) { double x2 = x*x;   return x2*x2*x;  }
double gsl_pow_6(const double x) { double x2 = x*x;   return x2*x2*x2; }
double gsl_pow_7(const double x) { double x3 = x*x*x; return x3*x3*x;  }
double gsl_pow_8(const double x) { double x2 = x*x;   double x4 = x2*x2; return x4*x4; }
double gsl_pow_9(const double x) { double x3 = x*x*x; return x3*x3*x3; }
#endif

double gsl_pow_int(double x, int n)
{
  double value = 1.0;

  if(n < 0) {
    x = 1.0/x;
    n = -n;
  }

  /* repeated squaring method
   * returns 0.0^0 = 1.0, so continuous in x
   */
  do {
     if(n & 1) value *= x;  /* for n odd */
     n >>= 1;
     x *= x;
  } while (n);

  return value;
}


text:
/* specfunc/gsl_sf_mathieu.h
 *
 * Copyright (C) 2002 Lowell Johnson
 *
 * 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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Author:  L. Johnson */

#ifndef __GSL_SF_MATHIEU_H__
#define __GSL_SF_MATHIEU_H__

#include <gsl/gsl_sf_result.h>
#include <gsl/gsl_eigen.h>

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

__BEGIN_DECLS

#define GSL_SF_MATHIEU_COEFF 100

typedef struct
{
  size_t size;
  size_t even_order;
  size_t odd_order;
  int extra_values;
  double qa;   /* allow for caching of results: not implemented yet */
  double qb;   /* allow for caching of results: not implemented yet */
  double *aa;
  double *bb;
  double *dd;
  double *ee;
  double *tt;
  double *e2;
  double *zz;
  gsl_vector *eval;
  gsl_matrix *evec;
  gsl_eigen_symmv_workspace *wmat;
} gsl_sf_mathieu_workspace;

/* Compute an array of characteristic (eigen) values from the recurrence
   matrices for the Mathieu equations. */
int gsl_sf_mathieu_a_array(int order_min, int order_max, double qq, gsl_sf_mathieu_workspace *work, double result_array[]);
int gsl_sf_mathieu_b_array(int order_min, int order_max, double qq,  gsl_sf_mathieu_workspace *work, double result_array[]);

/* Compute the characteristic value for a Mathieu function of order n and
   type ntype. */
int gsl_sf_mathieu_a(int order, double qq, gsl_sf_result *result);
int gsl_sf_mathieu_b(int order, double qq, gsl_sf_result *result);

/* Compute the Fourier coefficients for a Mathieu function. */
int gsl_sf_mathieu_a_coeff(int order, double qq, double aa, double coeff[]);
int gsl_sf_mathieu_b_coeff(int order, double qq, double aa, double coeff[]);

/* Allocate computational storage space for eigenvalue solution. */
gsl_sf_mathieu_workspace *gsl_sf_mathieu_alloc(const size_t nn,
                                               const double qq);
void gsl_sf_mathieu_free(gsl_sf_mathieu_workspace *workspace);

/* Compute an angular Mathieu function. */
int gsl_sf_mathieu_ce(int order, double qq, double zz, gsl_sf_result *result);
int gsl_sf_mathieu_se(int order, double qq, double zz, gsl_sf_result *result);
int gsl_sf_mathieu_ce_array(int nmin, int nmax, double qq, double zz,
                            gsl_sf_mathieu_workspace *work,
                            double result_array[]);
int gsl_sf_mathieu_se_array(int nmin, int nmax, double qq, double zz,
                            gsl_sf_mathieu_workspace *work,
                            double result_array[]);

/* Compute a radial Mathieu function. */
int gsl_sf_mathieu_Mc(int kind, int order, double qq, double zz,
                      gsl_sf_result *result);
int gsl_sf_mathieu_Ms(int kind, int order, double qq, double zz,
                      gsl_sf_result *result);
int gsl_sf_mathieu_Mc_array(int kind, int nmin, int nmax, double qq,
                            double zz, gsl_sf_mathieu_workspace *work,
                            double result_array[]);
int gsl_sf_mathieu_Ms_array(int kind, int nmin, int nmax, double qq,
                            double zz, gsl_sf_mathieu_workspace *work,
                            double result_array[]);

__END_DECLS

#endif /* !__GSL_SF_MATHIEU_H__ */

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: ielomariala/Hex-Game


text:
#ifndef BSPTN_H
#define BSPTN_H

#include <vector>
#include <gsl/gsl_sf_bessel.h>
#include "SpecialFunctions.h"

// NUMERICAL KERNELS /////

/// CREATION OF KERNEL ARRAYS ////

//Kernel array sizes

// 1st ORDER
extern double F1b_k1;
extern double G1b_k1;
extern double F1b_k2;
extern double G1b_k2;
extern double F1b_k3;
extern double G1b_k3;

extern double F1b_1;
extern double G1b_1;

extern double F1b_2;
extern double G1b_2;

extern double F1b_3;
extern double G1b_3;

// 2nd order
extern double F2b_k12;
extern double G2b_k12;
extern double F2b_k13;
extern double G2b_k13;
extern double F2b_k23;
extern double G2b_k23;

extern double F2b_p2mp;
extern double G2b_p2mp;

extern double F2b_p3mp;
extern double G2b_p3mp;

extern double F2b_12a;
extern double G2b_12a;

extern double F2b_13a;
extern double G2b_13a;

extern double F2b_23a;
extern double G2b_23a;

// 3rd order

extern double F3b_12mp;
extern double G3b_12mp;

extern double F3b_13mp;
extern double G3b_13mp;

extern double F3b_21mp;
extern double G3b_21mp;

extern double F3b_23mp;
extern double G3b_23mp;

extern double F3b_31mp;
extern double G3b_31mp;

extern double F3b_32mp;
extern double G3b_32mp;

extern double F3b_1pp;
extern double G3b_1pp;

extern double F3b_2pp;
extern double G3b_2pp;

extern double F3b_3pp;
extern double G3b_3pp;

// 4th order
extern double F4b_12pp;
extern double G4b_12pp;

extern double F4b_13pp;
extern double G4b_13pp;
extern double F4b_23pp;
extern double G4b_23pp;

const double myh0sqr = pow2(1./2997.92);

class BSPTN {
public:
  // tree level  equation solvers : dgp, f(R)
  void initnb0_dgp(double A, double k[], double x[], double kargs[], double omega0, double par1, double par2, double par3);
  void initnb0_fofr(double A, double k[], double x[], double kargs[], double omega0, double par1, double par2, double par3);

  // 1-loop level  equation solvers : dgp, f(R)
  void initnb1_dgp(double A, double k[], double x[], double kargs[], double omega0, double par1, double par2, double par3, double epars[]);
  void initnb1_fr(double A, double k[], double x[], double kargs[], double omega0, double par1, double par2, double par3, double epars[]);

};

#endif // BSPTN_H

max_stars_repo_name: PedroCarrilho/ReACT

max_issues_repo_name: PedroCarrilho/ReACT

max_forks_repo_name: PedroCarrilho/ReACT