text:
/* ieee-utils/fp-tru64.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Tim Mooney
 *
 * 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.
 */

/*
 * Under Compaq's Unix with the silly name, read the man pages for read_rnd,
 * write_rnd, and ieee(3) for more information on the functions used here.
 *
 * Note that enabling control of dynamic rounding mode (via write_rnd) requires
 * that you pass a special flag to your C compiler.  For Compaq's C compiler
 * the flag is `-fprm d', for gcc it's `-mfp-rounding-mode=d'.
 *
 * Enabling the trap control (via ieee_set_fp_control) also requires a
 * flag be passed to the C compiler.  The flag for Compaq's C compiler
 * is `-ieee' and for gcc it's `-mieee'.

 * We have not implemented the `inexact' case, since it is rarely used
 * and requires the library being built with an additional compiler
 * flag that can degrade performance for everything else. If you need
 * to add support for `inexact' the relevant flag for Compaq's
 * compiler is `-ieee_with_inexact', and the flag for gcc is
 * `-mieee-with-inexact'.
 *
 * Problem have been reported with the "fixed" float.h installed with
 * gcc-2.95 lacking some of the definitions in the system float.h (the
 * symptoms are errors like: `FP_RND_RN' undeclared). To work around
 * this we can include the system float.h before the gcc version, e.g.
 *
 *  #include "/usr/include/float.h"
 *  #include <float.h>
 */

#include <float.h>

#ifndef FP_RND_RN
#  undef _FLOAT_H_
#  include "/usr/include/float.h"
#  undef _FLOAT_H_
#  include <float.h>
#endif

#include <machine/fpu.h>
#include <stdio.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_errno.h>

int
gsl_ieee_set_mode (int precision, int rounding, int exception_mask)
{
  unsigned long int mode = 0 ;
  unsigned int    rnd  = 0 ;

/* I'm actually not completely sure that the alpha only supports default
 * precisions rounding, but I couldn't find any information regarding this, so
 * it seems safe to assume this for now until it's proven otherwise.
 */

  switch (precision)
    {
    case GSL_IEEE_SINGLE_PRECISION:
      GSL_ERROR ("Tru64 Unix on the alpha only supports default precision rounding",
                 GSL_EUNSUP) ;
      break ;
    case GSL_IEEE_DOUBLE_PRECISION:
      GSL_ERROR ("Tru64 Unix on the alpha only supports default precision rounding",
                 GSL_EUNSUP) ;
      break ;
    case GSL_IEEE_EXTENDED_PRECISION:
      GSL_ERROR ("Tru64 Unix on the alpha only supports default precision rounding",
                 GSL_EUNSUP) ;
      break ;
    }

  switch (rounding)
    {
    case GSL_IEEE_ROUND_TO_NEAREST:
      rnd = FP_RND_RN ;
      write_rnd (rnd) ;
      break ;
    case GSL_IEEE_ROUND_DOWN:
      rnd = FP_RND_RM ;
      write_rnd (rnd) ;
      break ;
    case GSL_IEEE_ROUND_UP:
      rnd = FP_RND_RP ;
      write_rnd (rnd) ;
      break ;
    case GSL_IEEE_ROUND_TO_ZERO:
      rnd = FP_RND_RZ ;
      write_rnd (rnd) ;
      break ;
    default:
      rnd = FP_RND_RN ;
      write_rnd (rnd) ;
    }

  /* Turn on all the exceptions apart from 'inexact' */

  /* from the ieee(3) man page:
   * IEEE_TRAP_ENABLE_INV       ->      Invalid operation
   * IEEE_TRAP_ENABLE_DZE       ->      Divide by 0
   * IEEE_TRAP_ENABLE_OVF       ->      Overflow
   * IEEE_TRAP_ENABLE_UNF       ->      Underflow
   * IEEE_TRAP_ENABLE_INE       ->      Inexact (requires special option to C compiler)
   * IEEE_TRAP_ENABLE_DNO       ->      denormal operand
   * Note: IEEE_TRAP_ENABLE_DNO is not supported on OSF 3.x or Digital Unix
   * 4.0 - 4.0d(?).
   * IEEE_TRAP_ENABLE_MASK      ->      mask of all the trap enables
   * IEEE_MAP_DMZ                       ->      map denormal inputs to zero
   * IEEE_MAP_UMZ                       ->      map underflow results to zero
   */

  mode = IEEE_TRAP_ENABLE_INV | IEEE_TRAP_ENABLE_DZE | IEEE_TRAP_ENABLE_OVF
                | IEEE_TRAP_ENABLE_UNF ;

  if (exception_mask & GSL_IEEE_MASK_INVALID)
    mode &= ~ IEEE_TRAP_ENABLE_INV ;

  if (exception_mask & GSL_IEEE_MASK_DENORMALIZED)
    {
#ifdef IEEE_TRAP_ENABLE_DNO
     mode &= ~ IEEE_TRAP_ENABLE_DNO ;
#else
     GSL_ERROR ("Sorry, this version of Digital Unix does not support denormalized operands", GSL_EUNSUP) ;
#endif
    }

  if (exception_mask & GSL_IEEE_MASK_DIVISION_BY_ZERO)
    mode &= ~ IEEE_TRAP_ENABLE_DZE ;

  if (exception_mask & GSL_IEEE_MASK_OVERFLOW)
    mode &= ~ IEEE_TRAP_ENABLE_OVF ;

  if (exception_mask & GSL_IEEE_MASK_UNDERFLOW)
    mode &=  ~ IEEE_TRAP_ENABLE_UNF ;

  if (exception_mask & GSL_IEEE_TRAP_INEXACT)
    {
      /* To implement this would require a special flag to the C
       compiler which can cause degraded performance */

      GSL_ERROR ("Sorry, GSL does not implement trap-inexact for Tru64 Unix on the alpha - see fp-tru64.c for details", GSL_EUNSUP) ;

      /* In case you need to add it, the appropriate line would be
       *
       *  mode |= IEEE_TRAP_ENABLE_INE ;
       *
       */

    }
  else
    {
      mode &= ~ IEEE_TRAP_ENABLE_INE ;
    }

  ieee_set_fp_control (mode) ;

  return GSL_SUCCESS ;
}


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

double f1(double x, void *params) {
    (void) (params); /* avoid unused parameter warning */
    return pow(x, 2);
}

double fi1(double x) {
    return pow(x, 3) / 3;
}

double f2(double x, void *params) {
    (void) (params); /* avoid unused parameter warning */
    if (0 == x) {
        return 0;
    }
    return 1 / sqrt(x);
}

double fi2(double x) {
    return 2 * sqrt(x);
}

double pi(double (*fi)(double), double x1, double x2) {
    return fi(x2) - fi(x1);
}

double e(double p1, double p2) {
    return fabs(p1 - p2);
}

int main(int argc, char *argv[]) {
    double a, b, pg1, e1, pi1, pg2, e2, pi2;
    if (argc < 2) {
        errno = EINVAL;
        perror("");
        return errno;
    }
    a = strtod(argv[1], NULL);
    b = strtod(argv[2], NULL);
    if (0 != errno) {
        perror("");
        return errno;
    }
    gsl_integration_workspace *w1
            = gsl_integration_workspace_alloc(1000);
    gsl_function F1;
    F1.function = &f1;
    gsl_integration_qags(&F1, a, b, 0, 1e-7, 1000,
                         w1, &pg1, &e1);
    pi1 = pi(fi1, a, b);
    gsl_integration_workspace *w2
            = gsl_integration_workspace_alloc(1000);
    gsl_function F2;
    F2.function = &f2;
    gsl_integration_qags(&F2, a, b, 0, 1e-7, 1000,
                         w2, &pg2, &e2);
    pi2 = pi(fi2, a, b);
    printf("%f\t%f\n", pg1, pg2);
    printf("%f\t%f\n", pi1, pi2);
    printf("%e\t%e\n", e(pg1, pi1), e(pg2, pi2));
    printf("%zu\t%zu\n", w1->size, w2->size);
    gsl_integration_workspace_free(w1);
    gsl_integration_workspace_free(w2);
    return EXIT_SUCCESS;
}

ext: c

max_stars_repo_name: mistyfiky/agh-mownit

max_issues_repo_name: mistyfiky/agh-mownit

max_forks_repo_name: mistyfiky/agh-mownit


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

void
cblas_sger (const enum CBLAS_ORDER order, const int M, const int N,
            const float alpha, const float *X, const int incX, const float *Y,
            const int incY, float *A, const int lda)
{
#define BASE float
#include "source_ger.h"
#undef BASE
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
/* ============================================================ *
 * lensing.h							*
 *								*
 * Martin Kilbinger, Karim Benabed 2006-2012			*
 * ============================================================ */

#ifndef __LENSING_H
#define __LENSING_H

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

#include <gsl/gsl_linalg.h>
#include <gsl/gsl_monte_plain.h>
#include <gsl/gsl_monte_miser.h>
#include <gsl/gsl_monte_vegas.h>
#include <gsl/gsl_rng.h>

#include "errorlist.h"
#include "maths.h"
#include "io.h"
#include "mvdens.h"
#include "par.h"

#include "hod.h"
#include "halomodel.h"
#include "cosmo.h"
#include "nofz.h"
#include "decomp_eb.h"
#include "reduced_fit.h"

/* Dimensions of interpolation tables */
/* N_s was increased from 200 to 400, for linear tabulation of P_kappa */
#define N_s     400
#define N_theta 100

/* Ranges of interpolation table for convergence power spectrum. *
 * Power-law extrapolation outside these ranges.		 */
#define s_min     1.0e-2
#define s_max     1.0e6

/* Ranges of interpolation table for reduced-shear correction    *
 * power spectrum. No extrapolation outside these ranges.	 */
#define ELL_MIN_REDUCED            0.1
#define ELL_MAX_REDUCED            2.0e5
#define THETA_P_MIN_REDUCED        (0.1*arcmin)
#define THETA_M_MIN_REDUCED        (0.5*arcmin)
#define THETA_MAP_MIN_REDUCED      (0.2*arcmin)
#define THETA_GSQR_MIN_REDUCED     (0.1*arcmin)
#define THETA_MAPGAUSS_MIN_REDUCED (0.1*arcmin)
#define THETA_MAX_REDUCED          (1000.0*arcmin)

#define NELL_REDUCED 50

#define lensing_base                     -1400
#define lensing_inconsistent             -1 + lensing_base
#define lensing_baryon_fraction          -2 + lensing_base
#define lensing_tomoij			 -3 + lensing_base
#define lensing_initialised              -4 + lensing_base
#define lensing_unknown                  -5 + lensing_base
#define lensing_pm                       -6 + lensing_base
#define lensing_type                     -7 + lensing_base
#define lensing_fastxi                   -8 + lensing_base
#define lensing_nperm                    -9 + lensing_base
#define lensing_range                   -10 + lensing_base
#define lensing_cosebi_n_max            -11 + lensing_base
#define lensing_ia                      -12 + lensing_base
#define lensing_angle_format            -13 + lensing_base
#define lensing_nzbin                   -14 + lensing_base

/* If Ob/Oc > BARYON_FRAC, chi2 produces an error */
#define BARYON_FRAC 0.75

/* Intrinsic alignment, constant amplitude C_1 * rho_crit, *
 * with C_1 = 5e-14 h^2 Mpc^3/M_sol.                       */
#define ia_c1_rho_crit     0.0134

typedef enum {xipm, xip, xim, map2poly, map2gauss, gsqr, decomp_eb, nofz, pkappa, map3gauss,
	      map3gauss_diag, map2gauss_map3gauss_diag, map2gauss_map3gauss,
	      decomp_eb_map3gauss_diag, decomp_eb_map3gauss}
  lensdata_t;
#define slensdata_t(i) ( \
 i==xipm      ? "xipm" : \
 i==xip       ? "xip"  : \
 i==xim       ? "xim"  : \
 i==map2poly  ? "map2poly" : \
 i==map2gauss ? "map2gauss" : \
 i==gsqr      ? "gsqr" : \
 i==decomp_eb ? "decomp_eb" : \
 i==nofz      ? "nofz" : \
 i==pkappa    ? "pkappa" : \
 i==map3gauss ? "map3gauss" : \
 i==map3gauss_diag ? "map3gauss_diag" : \
 i==map2gauss_map3gauss_diag ? "map2gauss_map3gauss_diag" : \
 i==map2gauss_map3gauss ? "map2gauss_map3gauss" : \
 i==decomp_eb_map3gauss_diag ? "decomp_eb_map3gauss_diag" : \
 i==decomp_eb_map3gauss ? "decomp_eb_map3gauss" : \
 "")
#define Nlensdata_t 15

typedef enum {decomp_eb_none, FK10_SN, FK10_FoM_eta10, FK10_FoM_eta50, COSEBIs_log} decomp_eb_filter_t;
#define sdecomp_eb_filter_t(i) (		\
 i==decomp_eb_none ? "none" : \
 i==FK10_SN        ? "FK10_SN" : \
 i==FK10_FoM_eta10 ? "FK10_FoM_eta10" : \
 i==FK10_FoM_eta50 ? "FK10_FoM_eta50" : \
 i==COSEBIs_log    ? "COSEBIs_log" : \
 "")
#define Ndecomp_eb_filter_t 5

/* The following arrays are defined in decomp_eb.c */
extern const double a_FK10_SN[], a_FK10_FoM_eta10[], a_FK10_FoM_eta50[];
// r_COSEB[];

typedef enum {angle_center, angle_mean, angle_wlinear, angle_wquadr} lensformat_t;
#define slensformat_t(i) ( \
 i==angle_center  ? "angle_center" : \
 i==angle_mean    ? "angle_mean" : \
 i==angle_wlinear ? "angle_wlinear" : \
 i==angle_wquadr  ? "angle_wquadr" : \
 "")
#define Nlensformat_t 4

typedef enum {cov_const, cov_ESH09} cov_scaling_t;
#define scov_scaling_t(i) ( \
  i==cov_const ? "cov_const" : \
  i==cov_ESH09 ? "cov_ESH09" : \
  "__undef__")
#define Ncov_scaling_t 2

typedef enum {reduced_none, reduced_K10} reduced_t;
#define sreduced_t(i) ( \
  i==reduced_none ? "none" : \
  i==reduced_K10  ? "K10" : \
  "")
#define Nreduced_t 2

/* Intrinsic alignment model */
typedef enum {ia_none, ia_HS04} ia_t;
#define sia_t(i) ( \
  i==ia_none ? "none" :	\
  i==ia_HS04 ? "HS04" : \
  "")
#define Nia_t 2

/* Bit-coded IA terms */
typedef enum {ia_undef, ia_GI_II, ia_only_GI, ia_only_II} ia_terms_t;
#define sia_terms_t(i) (      \
  i==ia_undef   ? "undef" :   \
  i==ia_GI_II   ? "GI_II" :   \
  i==ia_only_GI ? "only_GI" : \
  i==ia_only_II ? "only_II" : \
  "")
#define Nia_terms_t 4

typedef enum {second_order=2, third_order=3} order_t;

typedef enum {tomo_all, tomo_auto_only, tomo_cross_only} tomo_t;
#define stomo_t(i) ( \
 i==tomo_all        ? "tomo_all" : \
 i==tomo_auto_only  ? "tomo_auto_only" : \
 i==tomo_cross_only ? "tomo_cross_only" : \
 "")
#define Ntomo_t 3

typedef struct {
  int n_max;
  double th_min, th_max;
  char path[1024];
} cosebi_info_t;

typedef struct {

  /* Basic cosmology */
  cosmo *cosmo;

  /* Redshift distribution(s) */
  redshift_t *redshift;

  /* Tomography type */
  tomo_t tomo;

  /* Reduced-shear correction */
  reduced_t reduced;
  double q_mag_size;  /* q_mag_size = 2(alpha+beta-1),         *
		       * alpha, beta: slopes of number density *
		       * with flux (alpha), size (beta)        */

  /* Intrinsic aligmnent */
  ia_t ia;
  double A_ia;           /* IA amplitude */
  ia_terms_t ia_terms;  /* Bit-coded terms, GG=1, GI=2, II=4 */

  /* Halomodel stuff (only initialised if cosmo->nonlinear=halodm) */
  cosmo_hm *hm;

  /* ============================================================ *
   * Precomputed stuff.						  *
   * ============================================================ */

  interTable **g_source;
  interTable **Pshear, **Pg1;

  /* Shear second-order functions */
  interTable **xiP, **xiM, **gamma, **map_gauss, **map_poly;
  double *c_cosebi, psimin_cosebi, psimax_cosebi;
  int N_cosebi;

} cosmo_lens;

typedef struct {
  double r;
  cosmo_lens* self;
} cosmo_lensANDdouble;

typedef struct {
  int i;
  double r;
  cosmo_lens *self;
} cosmo_lensANDintANDdouble;

typedef struct {
  int i, j;
  double r;
  cosmo_lens *self;
} cosmo_lensANDiid;

typedef struct {
  int i, j, t;
  double r;
  cosmo_lens *self;
} cosmo_lensANDiiid;

typedef struct {
  int i_bin, j_bin, pm, n;
  const double *c;
  double thmin;
  cosmo_lens *self;
  error **err;
} cosmo_lensANDextra;

typedef struct {
   int Ntheta, Nzbin;  /* Number of angular and redshift bins */
   int Ntheta2;        /* For combined 2nd and 3rd-order */
   int Nzcorr;         /* Number of z-correlations, Nzcorr=Nzbin*(Nzbin+1)/2 */
   int n;              /* Number of total entries in data vector, n=Ntheta*Nzcorr */
   double *theta;      /* n-dimensional vector of angular scales */
   double *theta2;     /* For angle_range lensformats: (theta,theta2) = (lower,upper) bin limits */
   double *data;       /* n-dimensional data vector */
   double *var;        /* n-dimensional vector with variance */
   double *cov[3];     /* Maximum three nxn-dimensional covariance matrix */
   double a1, a2;      /* Coefficients for 'angle_wquadr' */
   double lndetC;
   int usecov;
   lensdata_t type;
   lensformat_t format;
   order_t order;
   decomp_eb_filter_t decomp_eb_filter;
   cov_scaling_t cov_scaling;
   cosmo_lens *fiducial;   /* Needed for ESH09 cov scaling */
} datcov;

/* ============================================================ *
 * Initialisation.						*
 * ============================================================ */

cosmo_lens *init_parameters_lens(double OMEGAM, double OMEGAV, double W0_DE, double W1_DE,
				 double *W_POLY_DE, int N_POLY_DE,
				 double H100, double OMEGAB, double OMEGANUMASS,
				 double NEFFNUMASS, double NORM, double NSPEC,
				 int Nzbin, const int *Nnz, const nofz_t *nofz, double *par_nz,
				 nonlinear_t NONLINEAR, transfer_t TRANSFER,
				 growth_t GROWTH, de_param_t DEPARAM,
				 norm_t normmode, tomo_t TOMO, reduced_t REDUCED, double Q_MAG_SIZE,
				 ia_t IA, ia_terms_t ia_terms, double A_IA, error **err);

void consistency_parameters_lens(const cosmo_lens *self, error **err);
cosmo_lens* copy_parameters_lens_only(cosmo_lens* source, error **err);
cosmo_lens* copy_parameters_lens(cosmo_lens* source, sm2_error **err);
void updateFrom_lens(cosmo_lens* avant, cosmo_lens* apres, error **err);
void copy_parameters_lenshm_cosmo(cosmo_lens *model, error **err);
void read_cosmological_parameters_lens(cosmo_lens **self, FILE *F, error **err);
cosmo_lens* set_cosmological_parameters_to_default_lens(error **err);
void free_parameters_lens(cosmo_lens** self);
void dump_param_lens(cosmo_lens* self, FILE *F, int wnofz, error **err);

/* ============================================================ *
 * Lensing functions.						*
 * ============================================================ */

/* Projection */
double int_for_g(double aprime, void *intpar, error **err);
double g_source(cosmo_lens*, double a, int n_bin, error **err);
double G(cosmo_lens* self, double a, int n_bin, error **err);
double int_for_p_2(double a, void *intpar,error **err);
double P_NL_tot(cosmo_lens *self, double a, double k, error **err);
double Pshear(cosmo_lens *self, double a, int i_bin, int j_bin, error **err);
double P_projected_kappa(void *self, double l, int i_bin, int j_bin, error **err);
double int_over_P_kappa(cosmo_lens *self, funcwithpars int_for_p, void *intpar, error **err);

double int_for_p_GI(double a, void *intpar, error **err);
double int_for_p_II(double a, void *intpar, error **err);

/* Reduced-shear correction (K10) */
extern const int parameter[M_PAR];
cosmo *set_cosmological_parameters_to_WMAP7(const redshift_t *nofz, tomo_t tomo, error **err);
double *par_to_pointer(cosmo *self, par_t par, error **err);
void fill_dpar(cosmo *model, cosmo *wmap7, double *dpar, error **err);
double Fbar(cosmo_lens *self, double a, int m_bin, int n_bin, error **err);
void fill_Fbar_array(cosmo_lens *self, double *fbar, int m_bin, int n_bin, double amin, int N_a,
		     double da, error **err);
void fill_dFbar_dp_array(cosmo_lens *self, par_t par, double *dfbar_dp, int m_bin, int n_bin, double amin,
			 int N_a, double da, error **err);
double Pg1(cosmo_lens *self, double s, int i_bin, int j_bin, error **err);

/* Second-order shear functions */
double xi(cosmo_lens*, int pm, double theta, int i_bin, int j_bin, error **err);
double gamma2(cosmo_lens*, double theta, int i_bin, int j_bin, error **err);
double map2_poly(cosmo_lens*, double theta, int i_bin, int j_bin, error **err);
double map2_gauss(cosmo_lens*, double theta, int i_bin, int j_bin, error **err);
double RR(cosmo_lens *lens, double THETA_MIN, double THETA_MAX, const double *a, int N,
	  poly_t poly, int pm, error **err);
double E_cosebi(cosmo_lens *lens, int n, double Psimin, double Psimax, int i_bin, int j_bin,
		const char *path, double *B_cosebi, error **err);
double RR_cosebi(cosmo_lens *lens, double THETA_MIN, double THETA_MAX, int i_bin, int j_bin,
		 int n, int pm, error **err);
double dRR_cosebi_dz_MC(double *z, int ndim, void *intpar);
double dRR_cosebi_dz(double z, void *intpar, error **err);
double int_for_map2_slow(double ell, void *intpar, error **err);
double map2_slow(cosmo_lens *self, double theta, tpstat_t tpstat, int i_bin, int j_bin, error **err);

/* Reading data files */
datcov *init_data_cov_tomo(char* dataname, char *dataname2, char** covname_ptr, lensdata_t type,
			   decomp_eb_filter_t decomp_eb_filter,
			   lensformat_t format, double corr_invcov,
			   double a1, double a2, order_t order,
			   cov_scaling_t cov_scaling, error **err);
datcov *init_datcov_for_cov_only(int Nzbin, int Ntheta, error **err);

void del_data_cov(datcov** dc);
void read_data_tomo(datcov *dc, char data_name[], int Nzbin, order_t order, error **err);
void read_cov_tomo(datcov* dc, char cov_name[], int icov, error **err);
void datcov2xipm(const datcov *dc, int i_bin, int j_bin, double **xip, double **xim, double **theta,
		 double **theta2, int *N, error **err);

//void read_cov(datcov* dc, char cov_name[], error **err);
void read_cov_col(datcov *dc, char cov_name[], error **err);

int get_pm(lensdata_t type, int i, int Ntheta, error **err);
int find_bin(double x, const double *list, int N, int prev, error **err);
void scale_cosmic_variance_ESH09(cosmo_lens *model, gsl_matrix *cov, const datcov *dc, error **err);
void scale_mixed_ESH09(const cosmo_lens *model, gsl_matrix *cov, const datcov *dc, error **err);
double lensing_signal(cosmo_lens *model, double theta, int i_bin, int j_bin, lensdata_t type,
		      decomp_eb_filter_t decomp_eb_filter, const cosebi_info_t *cosebi_info, error **err);
double chi2_lensing(cosmo_lens* csm, datcov* dc, int return_model, double **model_array, int *Nmodel,
		    const cosebi_info_t *cosebi_info, error **err);

/* Some third-order stuff which is called from lensing.c */
int Nperm_to_Ntheta(int Nperm, error **err);
void read_data_3rd(datcov *dc, char data_name[], int Nzbin, error **err);
void read_data_2nd_3rd(datcov *res, char *dataname, char *dataname2, error **err);

#define CHANGE(fct) int change_##fct(cosmo_lens*, cosmo_lens*)
CHANGE(g_source);
CHANGE(Pshear);
CHANGE(xi);
CHANGE(gamma2);
CHANGE(map2);
#undef CHANGE

#endif /* __LENSING_H */

ext: h

max_stars_repo_name: danielgruen/ccv

max_issues_repo_name: danielgruen/ccv

max_forks_repo_name: danielgruen/ccv


text:
/* linalg/test_lu_band.c
 *
 * Copyright (C) 2020 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 <config.h>
#include <stdlib.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_permute_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_permutation.h>

static int
test_LU_band_decomp_eps(const size_t p, const size_t q, const gsl_matrix * A, const double eps, const char * desc)
{
  int s = 0;
  const size_t M = A->size1;
  const size_t N = A->size2;
  const size_t minMN = GSL_MIN(M, N);
  size_t i, j;

  gsl_matrix * AB = gsl_matrix_alloc(N, 2*p + q + 1);
  gsl_vector_uint * piv = gsl_vector_uint_alloc(minMN);
  gsl_matrix * L = gsl_matrix_alloc(M, minMN);
  gsl_matrix * U = gsl_matrix_alloc(minMN, N);
  gsl_matrix * B = gsl_matrix_alloc(M, N);

  /* convert A to packed banded format */
  gen2band_matrix(p, q, A, AB);

  s += gsl_linalg_LU_band_decomp(M, p, q, AB, piv);
  s += gsl_linalg_LU_band_unpack(M, p, q, AB, piv, L, U);

  /* compute B = L U */
  gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, L, U, 0.0, B);

  for (i = 0; i < M; ++i)
    {
      for (j = 0; j < N; ++j)
        {
          double aij = gsl_matrix_get(A, i, j);
          double bij = gsl_matrix_get(B, i, j);

          gsl_test_rel(bij, aij, eps, "%s (M=%lu,N=%lu)(p=%lu,q=%lu)[%lu,%lu]: %22.18g   %22.18g\n",
                       desc, M, N, p, q, i, j, aij, bij);
        }
    }

  gsl_matrix_free(AB);
  gsl_vector_uint_free(piv);
  gsl_matrix_free(L);
  gsl_matrix_free(U);
  gsl_matrix_free(B);

  return s;
}

static int
test_LU_band_decomp(gsl_rng * r)
{
  int s = 0;
  const size_t N_max = 20;
  gsl_matrix * A = gsl_matrix_alloc(N_max, N_max);
  size_t M, N, p, q;

  for (M = 1; M <= N_max; ++M)
    {
      for (N = 1; N <= N_max; ++N)
        {
          gsl_matrix_view m = gsl_matrix_submatrix(A, 0, 0, M, N);

          for (p = 0; p < GSL_MIN(M, 10); ++p)
            {
              for (q = 0; q < GSL_MIN(N, 10); ++q)
                {
                  create_band_matrix(p, q, &m.matrix, r);
                  s += test_LU_band_decomp_eps(p, q, &m.matrix, 1.0e5 * GSL_MAX(M,N) * GSL_DBL_EPSILON, "LU_band_decomp random");
                }
            }
        }
    }

  gsl_matrix_free(A);

  return s;
}

static int
test_LU_band_solve_eps(const size_t p, const size_t q, const gsl_matrix * A,
                       const gsl_vector * rhs, const gsl_vector * sol,
                       const double eps, const char * desc)
{
  int s = 0;
  const size_t N = A->size2;
  size_t i;

  gsl_matrix * AB = gsl_matrix_alloc(N, 2*p + q + 1);
  gsl_vector_uint * piv = gsl_vector_uint_alloc(N);
  gsl_matrix * L = gsl_matrix_alloc(N, N);
  gsl_matrix * U = gsl_matrix_alloc(N, N);
  gsl_vector * x = gsl_vector_alloc(N);

  /* convert A to packed banded format */
  gen2band_matrix(p, q, A, AB);

  s += gsl_linalg_LU_band_decomp(N, p, q, AB, piv);
  s += gsl_linalg_LU_band_solve(p, q, AB, piv, rhs, x);

  for (i = 0; i < N; i++)
    {
      double xi = gsl_vector_get(x, i);
      double yi = gsl_vector_get(sol, i);

      gsl_test_rel(xi, yi, eps, "%s: %3lu[%lu]: %22.18g    %22.18g\n",
                   desc, N, i, xi, yi);
    }

  gsl_matrix_free(AB);
  gsl_vector_uint_free(piv);
  gsl_matrix_free(L);
  gsl_matrix_free(U);
  gsl_vector_free(x);

  return s;
}

static int
test_LU_band_solve(gsl_rng * r)
{
  int s = 0;
  const size_t N_max = 20;
  gsl_matrix * A = gsl_matrix_alloc(N_max, N_max);
  gsl_vector * b = gsl_vector_alloc(N_max);
  gsl_vector * x = gsl_vector_alloc(N_max);
  size_t N, p, q;

  for (N = 1; N <= N_max; ++N)
    {
      gsl_matrix_view m = gsl_matrix_submatrix(A, 0, 0, N, N);
      gsl_vector_view sol = gsl_vector_subvector(x, 0, N);
      gsl_vector_view rhs = gsl_vector_subvector(b, 0, N);

      for (p = 0; p < GSL_MIN(N, 10); ++p)
        {
          for (q = 0; q < GSL_MIN(N, 10); ++q)
            {
              create_band_matrix(p, q, &m.matrix, r);
              create_random_vector(&sol.vector, r);
              gsl_blas_dgemv(CblasNoTrans, 1.0, &m.matrix, &sol.vector, 0.0, &rhs.vector);

              s += test_LU_band_solve_eps(p, q, &m.matrix, &rhs.vector, &sol.vector,
                                          1.0e8 * N * GSL_DBL_EPSILON, "LU_band_solve random");
            }
        }
    }

  gsl_matrix_free(A);
  gsl_vector_free(b);
  gsl_vector_free(x);

  return s;
}

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
/*
 *  ht_neuron.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/>.
 *
 */

#ifndef HT_NEURON_H
#define HT_NEURON_H

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

#ifdef HAVE_GSL

// C++ includes:
#include <string>
#include <vector>

// C includes:
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>

// Includes from nestkernel:
#include "archiving_node.h"
#include "connection.h"
#include "recordables_map.h"
#include "ring_buffer.h"
#include "universal_data_logger.h"

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

/* BeginDocumentation
   Name: ht_neuron - Neuron model after Hill & Tononi (2005).

   Description:
   This model neuron implements a slightly modified version of the
   neuron model described in [1]. The most important properties are:

   - Integrate-and-fire with threshold adaptive threshold.
   - Repolarizing potassium current instead of hard reset.
   - AMPA, NMDA, GABA_A, and GABA_B conductance-based synapses with
     beta-function (difference of exponentials) time course.
   - Voltage-dependent NMDA with instantaneous or two-stage unblocking [1, 2].
   - Intrinsic currents I_h, I_T, I_Na(p), and I_KNa.
   - Synaptic "minis" are not implemented.

   Documentation and Examples:
   - docs/model_details/HillTononiModels.ipynb
   - pynest/examples/intrinsic_currents_spiking.py
   - pynest/examples/intrinsic_currents_subthreshold.py

   Parameters:
   V_m            - membrane potential
   tau_m          - membrane time constant applying to all currents except
                    repolarizing K-current (see [1], p 1677)
   t_ref          - refractory time and duration of post-spike repolarizing
                    potassium current (t_spike in [1])
   tau_spike      - membrane time constant for post-spike repolarizing
                    potassium current
   voltage_clamp  - if true, clamp voltage to value at beginning of simulation
                    (default: false, mainly for testing)
   theta, theta_eq, tau_theta - threshold, equilibrium value, time constant
   g_KL, E_K, g_NaL, E_Na     - conductances and reversal potentials for K and
                                Na leak currents
   {E_rev,g_peak,tau_rise,tau_decay}_{AMPA,NMDA,GABA_A,GABA_B}
                                - reversal potentials, peak conductances and
                                  time constants for synapses (tau_rise/
                                  tau_decay correspond to tau_1/tau_2 in the
                                  paper)
   V_act_NMDA, S_act_NMDA, tau_Mg_{fast, slow}_NMDA
                                - parameters for voltage dependence of NMDA-
                                  conductance, see above
   instant_unblock_NMDA         - instantaneous NMDA unblocking (default: false)
   {E_rev,g_peak}_{h,T,NaP,KNa} - reversal potential and peak conductance for
                                  intrinsic currents
   tau_D_KNa                    - relaxation time constant for I_KNa
   receptor_types               - dictionary mapping synapse names to ports on
                                  neuron model
   recordables                  - list of recordable quantities
   equilibrate                  - if given and true, time-dependent activation
                                  and inactivation state variables (h, m) of
                                  intrinsic currents and NMDA channels are set
                                  to their equilibrium values during this
                                  SetStatus call; otherwise they retain their
                                  present values.

   Note: Conductances are unitless in this model and currents are in mV.

   Author: Hans Ekkehard Plesser

   Sends: SpikeEvent
   Receives: SpikeEvent, CurrentEvent, DataLoggingRequest

   FirstVersion: October 2009; full revision November 2016

   References:
   [1] S Hill and G Tononi (2005). J Neurophysiol 93:1671-1698.
   [2] M Vargas-Caballero HPC Robinson (2003). J Neurophysiol 89:2778-2783.

   SeeAlso: ht_synapse
*/

namespace nest
{
/**
 * Function computing right-hand side of ODE for GSL solver.
 * @note Must be declared here so we can befriend it in class.
 * @note Must have C-linkage for passing to GSL. Internally, it is
 *       a first-class C++ function, but cannot be a member function
 *       because of the C-linkage.
 * @note No point in declaring it inline, since it is called
 *       through a function pointer.
 * @param void* Pointer to model neuron instance.
 */
extern "C" int ht_neuron_dynamics( double, const double*, double*, void* );

class ht_neuron : public Archiving_Node
{
public:
  ht_neuron();
  ht_neuron( const ht_neuron& );
  ~ht_neuron();

  /**
   * Import sets of overloaded virtual functions.
   * @see Technical Issues / Virtual Functions: Overriding, Overloading, and
   * Hiding
   */
  using Node::handle;
  using Node::handles_test_event;

  port send_test_event( Node&, rport, synindex, bool );

  void handle( SpikeEvent& e );
  void handle( CurrentEvent& e );
  void handle( DataLoggingRequest& );

  port handles_test_event( SpikeEvent&, rport );
  port handles_test_event( CurrentEvent&, rport );
  port handles_test_event( DataLoggingRequest&, rport );

  void get_status( DictionaryDatum& ) const;
  void set_status( const DictionaryDatum& );

private:
  /**
   * Synapse types to connect to
   * @note Excluded upper and lower bounds are defined as INF_, SUP_.
   *       Excluding port 0 avoids accidental connections.
   */
  enum SynapseTypes
  {
    INF_SPIKE_RECEPTOR = 0,
    AMPA,
    NMDA,
    GABA_A,
    GABA_B,
    SUP_SPIKE_RECEPTOR
  };

  void init_state_( const Node& proto );
  void init_buffers_();
  void calibrate();

  void update( Time const&, const long, const long );

  double get_synapse_constant( double, double, double );

  // END Boilerplate function declarations ----------------------------

  // Friends --------------------------------------------------------

  // make dynamics function quasi-member
  friend int ht_neuron_dynamics( double, const double*, double*, void* );

  // ----------------------------------------------------------------

  /**
   * Independent parameters of the model.
   */
  struct Parameters_
  {
    Parameters_();

    void get( DictionaryDatum& ) const; //!< Store current values in dictionary
    void set( const DictionaryDatum& ); //!< Set values from dicitonary

    // Note: Conductances are unitless
    // Leaks
    double E_Na; // mV
    double E_K;  // mV
    double g_NaL;
    double g_KL;
    double tau_m; // ms

    // Dynamic threshold
    double theta_eq;  // mV
    double tau_theta; // ms

    // Post-spike potassium current
    double tau_spike; // ms, membrane time constant for this current
    double t_ref;     // ms, refractory time

    // Parameters for synapse of type AMPA, GABA_A, GABA_B and NMDA
    double g_peak_AMPA;
    double tau_rise_AMPA;  // ms
    double tau_decay_AMPA; // ms
    double E_rev_AMPA;     // mV

    double g_peak_NMDA;
    double tau_rise_NMDA;  // ms
    double tau_decay_NMDA; // ms
    double E_rev_NMDA;     // mV
    double V_act_NMDA;     // mV, inactive for V << Vact, inflection of sigmoid
    double S_act_NMDA;     // mV, scale of inactivation
    double tau_Mg_slow_NMDA; // ms
    double tau_Mg_fast_NMDA; // ms
    bool instant_unblock_NMDA;

    double g_peak_GABA_A;
    double tau_rise_GABA_A;  // ms
    double tau_decay_GABA_A; // ms
    double E_rev_GABA_A;     // mV

    double g_peak_GABA_B;
    double tau_rise_GABA_B;  // ms
    double tau_decay_GABA_B; // ms
    double E_rev_GABA_B;     // mV

    // parameters for intrinsic currents
    double g_peak_NaP;
    double E_rev_NaP; // mV

    double g_peak_KNa;
    double E_rev_KNa; // mV
    double tau_D_KNa; // ms

    double g_peak_T;
    double E_rev_T; // mV

    double g_peak_h;
    double E_rev_h; // mV

    bool voltage_clamp;
  };

  // ----------------------------------------------------------------

  /**
   * State variables of the model.
   */
public:
  struct State_
  {

    // y_ = [V, theta, Synapses]
    enum StateVecElems_
    {
      V_M = 0,
      THETA,
      DG_AMPA,
      G_AMPA,
      DG_NMDA_TIMECOURSE,
      G_NMDA_TIMECOURSE,
      DG_GABA_A,
      G_GABA_A,
      DG_GABA_B,
      G_GABA_B, // DO NOT INSERT ANYTHING UP TO HERE, WILL MIX UP
                // SPIKE DELIVERY
      m_fast_NMDA,
      m_slow_NMDA,
      m_Ih,
      D_IKNa,
      m_IT,
      h_IT,
      STATE_VEC_SIZE
    };

    //! neuron state, must be C-array for GSL solver
    double y_[ STATE_VEC_SIZE ];

    /** Timer (counter) for spike-activated repolarizing potassium current.
     * Neuron is absolutely refractory during this period.
     */
    long ref_steps_;

    double I_NaP_; //!< Persistent Na current; member only to allow recording
    double I_KNa_; //!< Depol act. K current; member only to allow recording
    double I_T_;   //!< Low-thresh Ca current; member only to allow recording
    double I_h_;   //!< Pacemaker current; member only to allow recording

    State_( const ht_neuron&, const Parameters_& p );
    State_( const State_& s );
    ~State_();

    State_& operator=( const State_& s );

    void get( DictionaryDatum& ) const;
    void set( const DictionaryDatum&, const ht_neuron& );
  };

private:
  // These friend declarations must be precisely here.
  friend class RecordablesMap< ht_neuron >;
  friend class UniversalDataLogger< ht_neuron >;

  // ----------------------------------------------------------------

  /**
   * Buffers of the model.
   */
  struct Buffers_
  {
    Buffers_( ht_neuron& );
    Buffers_( const Buffers_&, ht_neuron& );

    UniversalDataLogger< ht_neuron > logger_;

    /** buffers and sums up incoming spikes/currents */
    std::vector< RingBuffer > spike_inputs_;
    RingBuffer currents_;

    /** GSL ODE stuff */
    gsl_odeiv_step* s_;    //!< stepping function
    gsl_odeiv_control* c_; //!< adaptive stepsize control function
    gsl_odeiv_evolve* e_;  //!< evolution function
    gsl_odeiv_system sys_; //!< struct describing system

    // IntergrationStep_ should be reset with the neuron on ResetNetwork,
    // but remain unchanged during calibration. Since it is initialized with
    // step_, and the resolution cannot change after nodes have been created,
    // it is safe to place both here.
    double step_;             //!< step size in ms
    double integration_step_; //!< current integration time step, updated by GSL

    /**
     * Input current injected by CurrentEvent.
     * This variable is used to transport the current applied into the
     * _dynamics function computing the derivative of the state vector.
     * It must be a part of Buffers_, since it is initialized once before
     * the first simulation, but not modified before later Simulate calls.
     */
    double I_stim_;
  };

  // ----------------------------------------------------------------

  /**
   * Internal variables of the model.
   */
  struct Variables_
  {
    //! size of conductance steps for arriving spikes
    std::vector< double > cond_steps_;

    //! Duration of potassium current.
    int PotassiumRefractoryCounts_;

    //! Voltage at beginning of simulation, for clamping
    double V_clamp_;
  };

  // readout functions, can use template for vector elements
  template < State_::StateVecElems_ elem >
  double
  get_y_elem_() const
  {
    return S_.y_[ elem ];
  }
  double
  get_I_NaP_() const
  {
    return S_.I_NaP_;
  }
  double
  get_I_KNa_() const
  {
    return S_.I_KNa_;
  }
  double
  get_I_T_() const
  {
    return S_.I_T_;
  }
  double
  get_I_h_() const
  {
    return S_.I_h_;
  }

  double get_g_NMDA_() const;

  /**
   * NMDA activation for given parameters
   * Needs to take parameter values explicitly since it is called from
   * _dynamics.
   */
  double m_NMDA_( double V, double m_eq, double m_fast, double m_slow ) const;

  /**
   * Return equilibrium value of I_h activation
   *
   * @param V Membrane potential for which to evaluate
   *        (may differ from y_[V_M] when clamping)
   */
  double m_eq_h_( double V ) const;

  /**
   * Return equilibrium value of I_T activation
   *
   * @param V Membrane potential for which to evaluate
   *        (may differ from y_[V_M] when clamping)
   */
  double m_eq_T_( double V ) const;

  /**
   * Return equilibrium value of I_T inactivation
   *
   * @param V Membrane potential for which to evaluate
   *        (may differ from y_[V_M] when clamping)
   */
  double h_eq_T_( double V ) const;

  /**
   * Return steady-state magnesium unblock ratio.
   *
   * Receives V_m as argument since it is called from ht_neuron_dyamics
   * with temporary state values.
   */
  double m_eq_NMDA_( double V ) const;

  /**
   * Steady-state "D" value for given voltage.
   */
  double D_eq_KNa_( double V ) const;

  static RecordablesMap< ht_neuron > recordablesMap_;

  Parameters_ P_;
  State_ S_;
  Variables_ V_;
  Buffers_ B_;
};

inline port
ht_neuron::send_test_event( Node& target, rport receptor_type, synindex, bool )
{
  SpikeEvent e;
  e.set_sender( *this );

  return target.handles_test_event( e, receptor_type );
}

inline port
ht_neuron::handles_test_event( SpikeEvent&, rport receptor_type )
{
  assert( B_.spike_inputs_.size() == 4 );

  if ( not( INF_SPIKE_RECEPTOR < receptor_type
         && receptor_type < SUP_SPIKE_RECEPTOR ) )
  {
    throw UnknownReceptorType( receptor_type, get_name() );
    return 0;
  }
  else
  {
    return receptor_type - 1;
  }

  /*
if (receptor_type != 0)
{
  throw UnknownReceptorType(receptor_type, get_name());
}
return 0;*/
}

inline port
ht_neuron::handles_test_event( CurrentEvent&, rport receptor_type )
{
  if ( receptor_type != 0 )
  {
    throw UnknownReceptorType( receptor_type, get_name() );
  }
  return 0;
}

inline port
ht_neuron::handles_test_event( DataLoggingRequest& dlr, rport receptor_type )
{
  if ( receptor_type != 0 )
  {
    throw UnknownReceptorType( receptor_type, get_name() );
  }
  return B_.logger_.connect_logging_device( dlr, recordablesMap_ );
}
}

#endif // HAVE_GSL
#endif // HT_NEURON_H

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:
// @(#)root/matrix:$Id$
// Authors: Fons Rademakers, Eddy Offermann   Nov 2003

/*************************************************************************
 * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers.               *
 * All rights reserved.                                                  *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/LICENSE.                         *
 * For the list of contributors see $ROOTSYS/README/CREDITS.             *
 *************************************************************************/

#ifndef ROOT_TMatrixT
#define ROOT_TMatrixT

//////////////////////////////////////////////////////////////////////////
//                                                                      //
// TMatrixT                                                             //
//                                                                      //
// Template class of a general matrix in the linear algebra package     //
//                                                                      //
//////////////////////////////////////////////////////////////////////////

#include "TMatrixTBase.h"
#include "TMatrixTUtils.h"

#ifdef CBLAS
#include <vecLib/vBLAS.h>
//#include <cblas.h>
#endif

#include "Rtypes.h"
#include "TError.h"

template<class Element> class TMatrixTSym;
template<class Element> class TMatrixTSparse;
template<class Element> class TMatrixTLazy;

template<class Element> class TMatrixT : public TMatrixTBase<Element> {

protected:

   Element  fDataStack[TMatrixTBase<Element>::kSizeMax]; //! data container
   Element *fElements;                                   //[fNelems] elements themselves

   Element *New_m   (Int_t size);
   void     Delete_m(Int_t size,Element*&);
   Int_t    Memcpy_m(Element *newp,const Element *oldp,Int_t copySize,
                      Int_t newSize,Int_t oldSize);
   void     Allocate(Int_t nrows,Int_t ncols,Int_t row_lwb = 0,Int_t col_lwb = 0,Int_t init = 0,
                     Int_t /*nr_nonzeros*/ = -1);

public:

   enum {kWorkMax = 100};
   enum EMatrixCreatorsOp1 { kZero,kUnit,kTransposed,kInverted,kAtA };
   enum EMatrixCreatorsOp2 { kMult,kTransposeMult,kInvMult,kMultTranspose,kPlus,kMinus };

   TMatrixT(): fDataStack(), fElements(0) { }
   TMatrixT(Int_t nrows,Int_t ncols);
   TMatrixT(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb);
   TMatrixT(Int_t nrows,Int_t ncols,const Element *data,Option_t *option="");
   TMatrixT(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,const Element *data,Option_t *option="");
   TMatrixT(const TMatrixT      <Element> &another);
   TMatrixT(const TMatrixTSym   <Element> &another);
   TMatrixT(const TMatrixTSparse<Element> &another);
   template <class Element2> TMatrixT(const TMatrixT<Element2> &another): fElements(0)
   {
      R__ASSERT(another.IsValid());
      Allocate(another.GetNrows(),another.GetNcols(),another.GetRowLwb(),another.GetColLwb());
      *this = another;
   }

   TMatrixT(EMatrixCreatorsOp1 op,const TMatrixT<Element> &prototype);
   TMatrixT(const TMatrixT    <Element> &a,EMatrixCreatorsOp2 op,const TMatrixT   <Element> &b);
   TMatrixT(const TMatrixT    <Element> &a,EMatrixCreatorsOp2 op,const TMatrixTSym<Element> &b);
   TMatrixT(const TMatrixTSym <Element> &a,EMatrixCreatorsOp2 op,const TMatrixT   <Element> &b);
   TMatrixT(const TMatrixTSym <Element> &a,EMatrixCreatorsOp2 op,const TMatrixTSym<Element> &b);
   TMatrixT(const TMatrixTLazy<Element> &lazy_constructor);

   virtual ~TMatrixT() { Clear(); }

   // Elementary constructors

   void Plus (const TMatrixT   <Element> &a,const TMatrixT   <Element> &b);
   void Plus (const TMatrixT   <Element> &a,const TMatrixTSym<Element> &b);
   void Plus (const TMatrixTSym<Element> &a,const TMatrixT   <Element> &b) { Plus(b,a); }

   void Minus(const TMatrixT   <Element> &a,const TMatrixT   <Element> &b);
   void Minus(const TMatrixT   <Element> &a,const TMatrixTSym<Element> &b);
   void Minus(const TMatrixTSym<Element> &a,const TMatrixT   <Element> &b) { Minus(b,a); }

   void Mult (const TMatrixT   <Element> &a,const TMatrixT   <Element> &b);
   void Mult (const TMatrixT   <Element> &a,const TMatrixTSym<Element> &b);
   void Mult (const TMatrixTSym<Element> &a,const TMatrixT   <Element> &b);
   void Mult (const TMatrixTSym<Element> &a,const TMatrixTSym<Element> &b);

   void TMult(const TMatrixT   <Element> &a,const TMatrixT   <Element> &b);
   void TMult(const TMatrixT   <Element> &a,const TMatrixTSym<Element> &b);
   void TMult(const TMatrixTSym<Element> &a,const TMatrixT   <Element> &b) { Mult(a,b); }
   void TMult(const TMatrixTSym<Element> &a,const TMatrixTSym<Element> &b) { Mult(a,b); }

   void MultT(const TMatrixT   <Element> &a,const TMatrixT   <Element> &b);
   void MultT(const TMatrixT   <Element> &a,const TMatrixTSym<Element> &b) { Mult(a,b); }
   void MultT(const TMatrixTSym<Element> &a,const TMatrixT   <Element> &b);
   void MultT(const TMatrixTSym<Element> &a,const TMatrixTSym<Element> &b) { Mult(a,b); }

   virtual const Element *GetMatrixArray  () const;
   virtual       Element *GetMatrixArray  ();
   virtual const Int_t   *GetRowIndexArray() const { return 0; }
   virtual       Int_t   *GetRowIndexArray()       { return 0; }
   virtual const Int_t   *GetColIndexArray() const { return 0; }
   virtual       Int_t   *GetColIndexArray()       { return 0; }

   virtual       TMatrixTBase<Element> &SetRowIndexArray(Int_t * /*data*/) { MayNotUse("SetRowIndexArray(Int_t *)"); return *this; }
   virtual       TMatrixTBase<Element> &SetColIndexArray(Int_t * /*data*/) { MayNotUse("SetColIndexArray(Int_t *)"); return *this; }

   virtual void Clear(Option_t * /*option*/ ="") { if (this->fIsOwner) Delete_m(this->fNelems,fElements);
                                                   else fElements = 0;
                                                   this->fNelems = 0; }

           TMatrixT    <Element> &Use     (Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,Element *data);
   const   TMatrixT    <Element> &Use     (Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,const Element *data) const
                                            { return (const TMatrixT<Element>&)
                                                     ((const_cast<TMatrixT<Element> *>(this))->Use(row_lwb,row_upb,col_lwb,col_upb, const_cast<Element *>(data))); }
           TMatrixT    <Element> &Use     (Int_t nrows,Int_t ncols,Element *data);
   const   TMatrixT    <Element> &Use     (Int_t nrows,Int_t ncols,const Element *data) const;
           TMatrixT    <Element> &Use     (TMatrixT<Element> &a);
   const   TMatrixT    <Element> &Use     (const TMatrixT<Element> &a) const;

   virtual TMatrixTBase<Element> &GetSub  (Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,
                                           TMatrixTBase<Element> &target,Option_t *option="S") const;
           TMatrixT    <Element>  GetSub  (Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,Option_t *option="S") const;
   virtual TMatrixTBase<Element> &SetSub  (Int_t row_lwb,Int_t col_lwb,const TMatrixTBase<Element> &source);

   virtual TMatrixTBase<Element> &ResizeTo(Int_t nrows,Int_t ncols,Int_t /*nr_nonzeros*/ =-1);
   virtual TMatrixTBase<Element> &ResizeTo(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,Int_t /*nr_nonzeros*/ =-1);
   inline  TMatrixTBase<Element> &ResizeTo(const TMatrixT<Element> &m) {
                                            return ResizeTo(m.GetRowLwb(),m.GetRowUpb(),m.GetColLwb(),m.GetColUpb());
                                 }

   virtual Double_t Determinant  () const;
   virtual void     Determinant  (Double_t &d1,Double_t &d2) const;

           TMatrixT<Element> &Invert      (Double_t *det=0);
           TMatrixT<Element> &InvertFast  (Double_t *det=0);
           TMatrixT<Element> &Transpose   (const TMatrixT<Element> &source);
   inline  TMatrixT<Element> &T           () { return this->Transpose(*this); }
           TMatrixT<Element> &Rank1Update (const TVectorT<Element> &v,Element alpha=1.0);
           TMatrixT<Element> &Rank1Update (const TVectorT<Element> &v1,const TVectorT<Element> &v2,Element alpha=1.0);
           Element            Similarity  (const TVectorT<Element> &v) const;

   TMatrixT<Element> &NormByColumn(const TVectorT<Element> &v,Option_t *option="D");
   TMatrixT<Element> &NormByRow   (const TVectorT<Element> &v,Option_t *option="D");

   // Either access a_ij as a(i,j)
   inline       Element                     operator()(Int_t rown,Int_t coln) const;
   inline       Element                    &operator()(Int_t rown,Int_t coln);

   // or as a[i][j]
   inline const TMatrixTRow_const<Element>  operator[](Int_t rown) const { return TMatrixTRow_const<Element>(*this,rown); }
   inline       TMatrixTRow      <Element>  operator[](Int_t rown)       { return TMatrixTRow      <Element>(*this,rown); }

   TMatrixT<Element> &operator= (const TMatrixT      <Element> &source);
   TMatrixT<Element> &operator= (const TMatrixTSym   <Element> &source);
   TMatrixT<Element> &operator= (const TMatrixTSparse<Element> &source);
   TMatrixT<Element> &operator= (const TMatrixTLazy  <Element> &source);
   template <class Element2> TMatrixT<Element> &operator= (const TMatrixT<Element2> &source)
   {
      if (!AreCompatible(*this,source)) {
         Error("operator=(const TMatrixT2 &)","matrices not compatible");
         return *this;
      }

     TObject::operator=(source);
     const Element2 * const ps = source.GetMatrixArray();
           Element  * const pt = this->GetMatrixArray();
     for (Int_t i = 0; i < this->fNelems; i++)
        pt[i] = ps[i];
     this->fTol = source.GetTol();
     return *this;
   }

   TMatrixT<Element> &operator= (Element val);
   TMatrixT<Element> &operator-=(Element val);
   TMatrixT<Element> &operator+=(Element val);
   TMatrixT<Element> &operator*=(Element val);

   TMatrixT<Element> &operator+=(const TMatrixT   <Element> &source);
   TMatrixT<Element> &operator+=(const TMatrixTSym<Element> &source);
   TMatrixT<Element> &operator-=(const TMatrixT   <Element> &source);
   TMatrixT<Element> &operator-=(const TMatrixTSym<Element> &source);

   TMatrixT<Element> &operator*=(const TMatrixT            <Element> &source);
   TMatrixT<Element> &operator*=(const TMatrixTSym         <Element> &source);
   TMatrixT<Element> &operator*=(const TMatrixTDiag_const  <Element> &diag);
   TMatrixT<Element> &operator/=(const TMatrixTDiag_const  <Element> &diag);
   TMatrixT<Element> &operator*=(const TMatrixTRow_const   <Element> &row);
   TMatrixT<Element> &operator/=(const TMatrixTRow_const   <Element> &row);
   TMatrixT<Element> &operator*=(const TMatrixTColumn_const<Element> &col);
   TMatrixT<Element> &operator/=(const TMatrixTColumn_const<Element> &col);

   const TMatrixT<Element> EigenVectors(TVectorT<Element> &eigenValues) const;

   ClassDef(TMatrixT,4) // Template of General Matrix class
};

#ifndef __CINT__
// When building with -fmodules, it instantiates all pending instantiations,
// instead of delaying them until the end of the translation unit.
// We 'got away with' probably because the use and the definition of the
// explicit specialization do not occur in the same TU.
//
// In case we are building with -fmodules, we need to forward declare the
// specialization in order to compile the dictionary G__Matrix.cxx.
template <> TClass *TMatrixT<double>::Class();
#endif // __CINT__

template <class Element> inline const Element           *TMatrixT<Element>::GetMatrixArray() const { return fElements; }
template <class Element> inline       Element           *TMatrixT<Element>::GetMatrixArray()       { return fElements; }

template <class Element> inline       TMatrixT<Element> &TMatrixT<Element>::Use           (Int_t nrows,Int_t ncols,Element *data)
                                                                                { return Use(0,nrows-1,0,ncols-1,data); }
template <class Element> inline const TMatrixT<Element> &TMatrixT<Element>::Use           (Int_t nrows,Int_t ncols,const Element *data) const
                                                                                { return Use(0,nrows-1,0,ncols-1,data); }
template <class Element> inline       TMatrixT<Element> &TMatrixT<Element>::Use           (TMatrixT &a)
                                                                                {
                                                                                R__ASSERT(a.IsValid());
                                                                                return Use(a.GetRowLwb(),a.GetRowUpb(),
                                                                                a.GetColLwb(),a.GetColUpb(),a.GetMatrixArray());
                                                                                }
template <class Element> inline const TMatrixT<Element> &TMatrixT<Element>::Use           (const TMatrixT &a) const
                                                                                {
                                                                                R__ASSERT(a.IsValid());
                                                                                return Use(a.GetRowLwb(),a.GetRowUpb(),
                                                                                a.GetColLwb(),a.GetColUpb(),a.GetMatrixArray());
                                                                                }

template <class Element> inline       TMatrixT<Element>  TMatrixT<Element>::GetSub        (Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,
                                                                                Option_t *option) const
                                                                                {
                                                                                TMatrixT tmp;
                                                                                this->GetSub(row_lwb,row_upb,col_lwb,col_upb,tmp,option);
                                                                                return tmp;
                                                                                }

template <class Element> inline Element TMatrixT<Element>::operator()(Int_t rown,Int_t coln) const
{
   R__ASSERT(this->IsValid());
   const Int_t arown = rown-this->fRowLwb;
   const Int_t acoln = coln-this->fColLwb;
   if (arown >= this->fNrows || arown < 0) {
      Error("operator()","Request row(%d) outside matrix range of %d - %d",rown,this->fRowLwb,this->fRowLwb+this->fNrows);
      return TMatrixTBase<Element>::NaNValue();
   }
   if (acoln >= this->fNcols || acoln < 0) {
      Error("operator()","Request column(%d) outside matrix range of %d - %d",coln,this->fColLwb,this->fColLwb+this->fNcols);
      return TMatrixTBase<Element>::NaNValue();

   }
   return (fElements[arown*this->fNcols+acoln]);
}

template <class Element> inline Element &TMatrixT<Element>::operator()(Int_t rown,Int_t coln)
{
   R__ASSERT(this->IsValid());
   const Int_t arown = rown-this->fRowLwb;
   const Int_t acoln = coln-this->fColLwb;
   if (arown >= this->fNrows || arown < 0) {
      Error("operator()","Request row(%d) outside matrix range of %d - %d",rown,this->fRowLwb,this->fRowLwb+this->fNrows);
      return TMatrixTBase<Element>::NaNValue();
   }
   if (acoln >= this->fNcols || acoln < 0) {
      Error("operator()","Request column(%d) outside matrix range of %d - %d",coln,this->fColLwb,this->fColLwb+this->fNcols);
      return TMatrixTBase<Element>::NaNValue();
   }
   return (fElements[arown*this->fNcols+acoln]);
}

template <class Element> TMatrixT<Element>  operator+  (const TMatrixT   <Element> &source1,const TMatrixT   <Element> &source2);
template <class Element> TMatrixT<Element>  operator+  (const TMatrixT   <Element> &source1,const TMatrixTSym<Element> &source2);
template <class Element> TMatrixT<Element>  operator+  (const TMatrixTSym<Element> &source1,const TMatrixT   <Element> &source2);
template <class Element> TMatrixT<Element>  operator+  (const TMatrixT   <Element> &source ,      Element               val    );
template <class Element> TMatrixT<Element>  operator+  (      Element               val    ,const TMatrixT   <Element> &source );
template <class Element> TMatrixT<Element>  operator-  (const TMatrixT   <Element> &source1,const TMatrixT   <Element> &source2);
template <class Element> TMatrixT<Element>  operator-  (const TMatrixT   <Element> &source1,const TMatrixTSym<Element> &source2);
template <class Element> TMatrixT<Element>  operator-  (const TMatrixTSym<Element> &source1,const TMatrixT   <Element> &source2);
template <class Element> TMatrixT<Element>  operator-  (const TMatrixT   <Element> &source ,      Element               val    );
template <class Element> TMatrixT<Element>  operator-  (      Element               val    ,const TMatrixT   <Element> &source );
template <class Element> TMatrixT<Element>  operator*  (      Element               val    ,const TMatrixT   <Element> &source );
template <class Element> TMatrixT<Element>  operator*  (const TMatrixT   <Element> &source ,      Element               val    );
template <class Element> TMatrixT<Element>  operator*  (const TMatrixT   <Element> &source1,const TMatrixT   <Element> &source2);
template <class Element> TMatrixT<Element>  operator*  (const TMatrixT   <Element> &source1,const TMatrixTSym<Element> &source2);
template <class Element> TMatrixT<Element>  operator*  (const TMatrixTSym<Element> &source1,const TMatrixT   <Element> &source2);
template <class Element> TMatrixT<Element>  operator*  (const TMatrixTSym<Element> &source1,const TMatrixTSym<Element> &source2);
// Preventing warnings with -Weffc++ in GCC since overloading the || and && operators was a design choice.
#if (__GNUC__ * 10000 + __GNUC_MINOR__ * 100 + __GNUC_PATCHLEVEL__) >= 40600
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Weffc++"
#endif
template <class Element> TMatrixT<Element>  operator&& (const TMatrixT   <Element> &source1,const TMatrixT   <Element> &source2);
template <class Element> TMatrixT<Element>  operator&& (const TMatrixT   <Element> &source1,const TMatrixTSym<Element> &source2);
template <class Element> TMatrixT<Element>  operator&& (const TMatrixTSym<Element> &source1,const TMatrixT   <Element> &source2);
template <class Element> TMatrixT<Element>  operator|| (const TMatrixT   <Element> &source1,const TMatrixT   <Element> &source2);
template <class Element> TMatrixT<Element>  operator|| (const TMatrixT   <Element> &source1,const TMatrixTSym<Element> &source2);
template <class Element> TMatrixT<Element>  operator|| (const TMatrixTSym<Element> &source1,const TMatrixT   <Element> &source2);
#if (__GNUC__ * 10000 + __GNUC_MINOR__ * 100 + __GNUC_PATCHLEVEL__) >= 40600
#pragma GCC diagnostic pop
#endif
template <class Element> TMatrixT<Element>  operator>  (const TMatrixT   <Element> &source1,const TMatrixT   <Element> &source2);
template <class Element> TMatrixT<Element>  operator>  (const TMatrixT   <Element> &source1,const TMatrixTSym<Element> &source2);
template <class Element> TMatrixT<Element>  operator>  (const TMatrixTSym<Element> &source1,const TMatrixT   <Element> &source2);
template <class Element> TMatrixT<Element>  operator>= (const TMatrixT   <Element> &source1,const TMatrixT   <Element> &source2);
template <class Element> TMatrixT<Element>  operator>= (const TMatrixT   <Element> &source1,const TMatrixTSym<Element> &source2);
template <class Element> TMatrixT<Element>  operator>= (const TMatrixTSym<Element> &source1,const TMatrixT   <Element> &source2);
template <class Element> TMatrixT<Element>  operator<= (const TMatrixT   <Element> &source1,const TMatrixT   <Element> &source2);
template <class Element> TMatrixT<Element>  operator<= (const TMatrixT   <Element> &source1,const TMatrixTSym<Element> &source2);
template <class Element> TMatrixT<Element>  operator<= (const TMatrixTSym<Element> &source1,const TMatrixT   <Element> &source2);
template <class Element> TMatrixT<Element>  operator<  (const TMatrixT   <Element> &source1,const TMatrixT   <Element> &source2);
template <class Element> TMatrixT<Element>  operator<  (const TMatrixT   <Element> &source1,const TMatrixTSym<Element> &source2);
template <class Element> TMatrixT<Element>  operator<  (const TMatrixTSym<Element> &source1,const TMatrixT   <Element> &source2);
template <class Element> TMatrixT<Element>  operator!= (const TMatrixT   <Element> &source1,const TMatrixT   <Element> &source2);
template <class Element> TMatrixT<Element>  operator!= (const TMatrixT   <Element> &source1,const TMatrixTSym<Element> &source2);
template <class Element> TMatrixT<Element>  operator!= (const TMatrixTSym<Element> &source1,const TMatrixT   <Element> &source2);

template <class Element> TMatrixT<Element> &Add        (TMatrixT<Element> &target,      Element               scalar,const TMatrixT   <Element> &source);
template <class Element> TMatrixT<Element> &Add        (TMatrixT<Element> &target,      Element               scalar,const TMatrixTSym<Element> &source);
template <class Element> TMatrixT<Element> &ElementMult(TMatrixT<Element> &target,const TMatrixT   <Element> &source);
template <class Element> TMatrixT<Element> &ElementMult(TMatrixT<Element> &target,const TMatrixTSym<Element> &source);
template <class Element> TMatrixT<Element> &ElementDiv (TMatrixT<Element> &target,const TMatrixT   <Element> &source);
template <class Element> TMatrixT<Element> &ElementDiv (TMatrixT<Element> &target,const TMatrixTSym<Element> &source);

template <class Element> void AMultB (const Element * const ap,Int_t na,Int_t ncolsa,
                                      const Element * const bp,Int_t nb,Int_t ncolsb,Element *cp);
template <class Element> void AtMultB(const Element * const ap,Int_t ncolsa,
                                      const Element * const bp,Int_t nb,Int_t ncolsb,Element *cp);
template <class Element> void AMultBt(const Element * const ap,Int_t na,Int_t ncolsa,
                                      const Element * const bp,Int_t nb,Int_t ncolsb,Element *cp);

#endif


text:
/**
 * \author Sylvain Marsat, University of Maryland - NASA GSFC
 *
 * \brief C code for structures EOBNRv2HM reduced order model (non-spinning version).
 * See CQG 31 195010, 2014, arXiv:1402.4146 for details on the reduced order method.
 * See arXiv:1106.1021 for the EOBNRv2HM model.
 *
 * Borrows from the SEOBNR ROM LAL code written by Michael Puerrer and John Veitch.
 *
 * Put the untared data in the directory designated by the environment variable ROM_DATA_PATH.
 *
 * Parameter range:
 *   q = 1-12 (almost)
 *   No spin
 *   Mtot >= 10Msun for fstart=8Hz
 *
 */

#define _XOPEN_SOURCE 500

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

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

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

#include "constants.h"
#include "EOBNRv2HMROMstruct.h"
#include "EOBNRv2HMROM.h"

/************************************************************************/
/********************* Functions for list structures ********************/

/***************** Functions for the SplineList structure ****************/

/* Prepend a node to a linked list of splines, or create a new head */
SplineList* SplineList_AddElementNoCopy(
	   SplineList* appended,  /* List structure to prepend to */
	   gsl_spline* spline,  /* spline to contain */
           gsl_interp_accel* accel,  /* accelerator to contain */
	   int i /* index in the list */)
{
    SplineList* splinelist;
    /* Check if the node with this index already exists */
    splinelist = appended;
    while( splinelist ){
      if( i == splinelist->i ){
	break;
      }
      splinelist = splinelist->next;
    }
    if( splinelist ){ /* We don't allow for the case where the index already exists*/
      printf("Error: Tried to add an already existing index to a SplineList");
      return(NULL);
    } else { /* In that case, we do NOT COPY the input spline, which therefore can't be
                 used anywhere else; this will be acceptable as these operations will only be done
                 when initializing the data */
      splinelist = malloc( sizeof(SplineList) );
    }
    splinelist->i = i;
    if( spline ){
      splinelist->spline = spline;
    } else {
      splinelist->spline = NULL;
    }
    if( accel ){
      splinelist->accel = accel;
    } else {
      splinelist->accel = NULL;
    }
    if( appended ){
      splinelist->next = appended;
    } else {
        splinelist->next = NULL;
    }
    return splinelist;
}
/* Get the element of a SplineList with a given index */
SplineList* SplineList_GetElement(
	      SplineList* const splinelist,  /* List structure to get element from */
              const int i ) /* Index looked for */
{
    if( !splinelist ) return NULL;

    SplineList* itr = splinelist;
    while( itr->i != i ){
        itr = itr->next;
        if( !itr ) return NULL;
    }
    return itr; /* The element returned is itself a pointer to a SplineList */
}
/* Delete list from given pointer to the end of the list */
void SplineList_Destroy( SplineList* splinelist ) /* Head of linked list to destroy */
{
  SplineList* pop;
  while( (pop = splinelist) ){
    if( pop->spline ){ /* Internal spline and accelerator are freed */
      gsl_spline_free( pop->spline );
    }
    if( pop->accel ){
      gsl_interp_accel_free( pop->accel );
    }
    /* Notice that the index i is not freed, like in SphHarmTimeSeries struct indices l and m */
    splinelist = pop->next;
    free( pop );
  }
}

/***************** Functions for the EOBNRHMROMdata structure ****************/
ListmodesEOBNRHMROMdata* ListmodesEOBNRHMROMdata_AddModeNoCopy(
	   ListmodesEOBNRHMROMdata* appended,  /* List structure to prepend to */
	   EOBNRHMROMdata* data,  /* data to contain */
	   const int l, /* major mode number */
	   const int m  /* minor mode number */)
{
    ListmodesEOBNRHMROMdata* list;
    /* Check if the node with this mode already exists */
    list = appended;
    while( list ){
      if( l == list->l && m == list->m ){
	break;
      }
      list = list->next;
    }
    if( list ){ /* We don't allow for the case where the mode already exists in the list*/
      printf("Error: Tried to add an already existing mode to a ListmodesEOBNRHMROMdata ");
      return(NULL);
    } else { /* In that case, we do NOT COPY the input interpolated data, which therefore can't be
		used anywhere else; this will be acceptable as these operations will only be done
		when interpolating the initialization data */
      list = malloc( sizeof(ListmodesEOBNRHMROMdata) );
    }
    list->l = l;
    list->m = m;
    if( data ){
      list->data = data;
    } else {
      list->data = NULL;
    }
    if( appended ){
      list->next = appended;
    } else {
        list->next = NULL;
    }
    return list;
}
/* Get the element of a ListmodesEOBNRHMROMdata with a given index */
ListmodesEOBNRHMROMdata* ListmodesEOBNRHMROMdata_GetMode(
	   ListmodesEOBNRHMROMdata* const list,  /* List structure to get a particular mode from */
	   int l, /*< major mode number */
	   int m  /*< minor mode number */ )
{
    if( !list ) return NULL;

    ListmodesEOBNRHMROMdata *itr = list;
    while( itr->l != l || itr->m != m ){
        itr = itr->next;
        if( !itr ) return NULL;
    }
    return itr; /* The element returned is itself a pointer to a ListmodesEOBNRHMROMdata */
}
void ListmodesEOBNRHMROMdata_Destroy(
	   ListmodesEOBNRHMROMdata* list  /* List structure to destroy; notice that the data is destroyed too */
)
{
  ListmodesEOBNRHMROMdata* pop;
  while( (pop = list) ){
    if( pop->data ){ /* Destroying the EOBNRHMROMdata data */
      EOBNRHMROMdata_Cleanup( pop->data );
    }
    /* Notice that the mode indices l and m are not freed, like in SphHarmTimeSeries struct indices l and m */
    list = pop->next;
    free( pop );
  }
}

/***************** Functions for the EOBNRHMROMdata_interp structure ****************/
ListmodesEOBNRHMROMdata_interp* ListmodesEOBNRHMROMdata_interp_AddModeNoCopy(
	   ListmodesEOBNRHMROMdata_interp* appended,  /* List structure to prepend to */
	   EOBNRHMROMdata_interp* data_interp,  /* data to contain */
	   int l, /* major mode number */
	   int m  /* minor mode number */)
{
    ListmodesEOBNRHMROMdata_interp* list;
    /* Check if the node with this mode already exists */
    list = appended;
    while( list ){
      if( l == list->l && m == list->m ){
	break;
      }
      list = list->next;
    }
    if( list ){ /* We don't allow for the case where the mode already exists in the list*/
      printf("Error: Tried to add an already existing mode to a ListmodesEOBNRHMROMdata_interp ");
      return(NULL);
    } else { /* In that case, we do NOT COPY the input interpolated data, which therefore can't be
		used anywhere else; this will be acceptable as these operations will only be done
		when interpolating the initialization data */
      list = malloc( sizeof(ListmodesEOBNRHMROMdata_interp) );
    }
    list->l = l;
    list->m = m;
    if( data_interp ){
      list->data_interp = data_interp;
    } else {
      list->data_interp = NULL;
    }
    if( appended ){
      list->next = appended;
    } else {
        list->next = NULL;
    }
    return list;
}
/* Get the element of a ListmodesEOBNRHMROMdata with a given index */
ListmodesEOBNRHMROMdata_interp* ListmodesEOBNRHMROMdata_interp_GetMode(
	   ListmodesEOBNRHMROMdata_interp* const list,  /* List structure to get a particular mode from */
	   int l, /*< major mode number */
	   int m  /*< minor mode number */ )
{
    if( !list ) return NULL;

    ListmodesEOBNRHMROMdata_interp *itr = list;
    while( itr->l != l || itr->m != m ){
        itr = itr->next;
        if( !itr ) return NULL;
    }
    return itr; /* The element returned is itself a pointer to a ListmodesEOBNRHMROMdata_interp */
}
void ListmodesEOBNRHMROMdata_interp_Destroy(
	   ListmodesEOBNRHMROMdata_interp* list  /* List structure to destroy; notice that the data is destroyed too */
)
{
  ListmodesEOBNRHMROMdata_interp* pop;
  while( (pop = list) ){
    if( pop->data_interp ){ /* Destroying the EOBNRHMROMdata_interp data */
      EOBNRHMROMdata_interp_Cleanup( pop->data_interp );
    }
    /* Notice that the mode indices l and m are not freed, like in SphHarmTimeSeries struct indices l and m */
    list = pop->next;
    free( pop );
  }
}

/***************** Functions for the EOBNRHMROMdata_coeff structure ****************/
ListmodesEOBNRHMROMdata_coeff* ListmodesEOBNRHMROMdata_coeff_AddModeNoCopy(
	   ListmodesEOBNRHMROMdata_coeff* appended,  /* List structure to prepend to */
	   EOBNRHMROMdata_coeff* data_coeff,  /* data to contain */
	   int l, /* major mode number */
	   int m  /* minor mode number */)
{
    ListmodesEOBNRHMROMdata_coeff* list;
    /* Check if the node with this mode already exists */
    list = appended;
    while( list ){
      if( l == list->l && m == list->m ){
	break;
      }
      list = list->next;
    }
    if( list ){ /* We don't allow for the case where the mode already exists in the list*/
      printf("Error: Tried to add an already existing mode to a ListmodesEOBNRHMROMdata_coeff ");
      return(NULL);
    } else { /* In that case, we do NOT COPY the input interpolated data, which therefore can't be
		used anywhere else; this will be acceptable as these operations will only be done
		when interpolating the initialization data */
      list = malloc( sizeof(ListmodesEOBNRHMROMdata_coeff) );
    }
    list->l = l;
    list->m = m;
    if( data_coeff ){
      list->data_coeff = data_coeff;
    } else {
      list->data_coeff = NULL;
    }
    if( appended ){
      list->next = appended;
    } else {
        list->next = NULL;
    }
    return list;
}
/* Get the element of a ListmodesEOBNRHMROMdata_coeff with a given index */
ListmodesEOBNRHMROMdata_coeff* ListmodesEOBNRHMROMdata_coeff_GetMode(
	   ListmodesEOBNRHMROMdata_coeff* const list,  /* List structure to get a particular mode from */
	   int l, /*< major mode number */
	   int m  /*< minor mode number */ )
{
    if( !list ) return NULL;

    ListmodesEOBNRHMROMdata_coeff *itr = list;
    while( itr->l != l || itr->m != m ){
        itr = itr->next;
        if( !itr ) return NULL;
    }
    return itr; /* The element returned is itself a pointer to a ListmodesEOBNRHMROMdata_coeff */
}
void ListmodesEOBNRHMROMdata_coeff_Destroy(
	   ListmodesEOBNRHMROMdata_coeff* list  /* List structure to destroy; notice that the data is destroyed too */
)
{
  ListmodesEOBNRHMROMdata_coeff* pop;
  while( (pop = list) ){
    if( pop->data_coeff ){ /* Destroying the EOBNRHMROMdata_coeff data */
      EOBNRHMROMdata_coeff_Cleanup( pop->data_coeff );
    }
    /* Notice that the mode indices l and m are not freed, like in SphHarmTimeSeries struct indices l and m */
    list = pop->next;
    free( pop );
  }
}

max_stars_repo_name: titodalcanton/flare

max_issues_repo_name: titodalcanton/flare

max_forks_repo_name: titodalcanton/flare


text:
#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

int main(int argc, char *argv[])
{
  gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);

  int i, n;
  double gauss, gamma;

  n = atoi(argv[1]);

  for (i=0; i<n; i++){
    gauss = gsl_ran_gaussian(r, 2.0);
    gamma = gsl_ran_gamma(r, 2.0, 3.0);
    printf("%2.4f %2.4f\n", gauss, gamma);
  }
  return 0;

}

max_stars_repo_name: Bingwen-Hu/hackaway

max_issues_repo_name: Bingwen-Hu/hackaway

max_forks_repo_name: Bingwen-Hu/hackaway


text:
/* histogram/stat2d.c
 * Copyright (C) 2002  Achim Gaedke
 *
 * This library 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 library; if not, write to the Free Software Foundation, Inc.,
 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

/***************************************************************
 *
 * File histogram/stat2d.c:
 * Routine to return statistical values of the content of a 2D hisogram.
 *
 * Contains the routines:
 * gsl_histogram2d_sum sum up all bin values
 * gsl_histogram2d_xmean determine mean of x values
 * gsl_histogram2d_ymean determine mean of y values
 *
 * Author: Achim Gaedke Achim.Gaedke@zpr.uni-koeln.de
 * Jan. 2002
 *
 ***************************************************************/

#include <config.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_histogram2d.h>

/*
  sum up all bins of histogram2d
 */

double
gsl_histogram2d_sum (const gsl_histogram2d * h)
{
  const size_t n = h->nx * h->ny;
  double sum = 0;
  size_t i = 0;

  while (i < n)
    sum += h->bin[i++];

  return sum;
}

double
gsl_histogram2d_xmean (const gsl_histogram2d * h)
{
  const size_t nx = h->nx;
  const size_t ny = h->ny;
  size_t i;
  size_t j;

  /* Compute the bin-weighted arithmetic mean M of a histogram using the
     recurrence relation

     M(n) = M(n-1) + (x[n] - M(n-1)) (w(n)/(W(n-1) + w(n)))
     W(n) = W(n-1) + w(n)

   */

  long double wmean = 0;
  long double W = 0;

  for (i = 0; i < nx; i++)
    {
      double xi = (h->xrange[i + 1] + h->xrange[i]) / 2.0;
      double wi = 0;

      for (j = 0; j < ny; j++)
        {
          double wij = h->bin[i * ny + j];
          if (wij > 0)
            wi += wij;
        }
      if (wi > 0)
        {
          W += wi;
          wmean += (xi - wmean) * (wi / W);
        }
    }

  return wmean;
}

double
gsl_histogram2d_ymean (const gsl_histogram2d * h)
{
  const size_t nx = h->nx;
  const size_t ny = h->ny;
  size_t i;
  size_t j;

  /* Compute the bin-weighted arithmetic mean M of a histogram using the
     recurrence relation

     M(n) = M(n-1) + (x[n] - M(n-1)) (w(n)/(W(n-1) + w(n)))
     W(n) = W(n-1) + w(n)

   */

  long double wmean = 0;
  long double W = 0;

  for (j = 0; j < ny; j++)
    {
      double yj = (h->yrange[j + 1] + h->yrange[j]) / 2.0;
      double wj = 0;

      for (i = 0; i < nx; i++)
        {
          double wij = h->bin[i * ny + j];
          if (wij > 0)
            wj += wij;
        }

      if (wj > 0)
        {
          W += wj;
          wmean += (yj - wmean) * (wj / W);
        }
    }

  return wmean;
}

double
gsl_histogram2d_xsigma (const gsl_histogram2d * h)
{
  const double xmean = gsl_histogram2d_xmean (h);
  const size_t nx = h->nx;
  const size_t ny = h->ny;
  size_t i;
  size_t j;

  /* Compute the bin-weighted arithmetic mean M of a histogram using the
     recurrence relation

     M(n) = M(n-1) + (x[n] - M(n-1)) (w(n)/(W(n-1) + w(n)))
     W(n) = W(n-1) + w(n)

   */

  long double wvariance = 0;
  long double W = 0;

  for (i = 0; i < nx; i++)
    {
      double xi = (h->xrange[i + 1] + h->xrange[i]) / 2 - xmean;
      double wi = 0;

      for (j = 0; j < ny; j++)
        {
          double wij = h->bin[i * ny + j];
          if (wij > 0)
            wi += wij;
        }

      if (wi > 0)
        {
          W += wi;
          wvariance += ((xi * xi) - wvariance) * (wi / W);
        }
    }

  {
    double xsigma = sqrt (wvariance);
    return xsigma;
  }
}

double
gsl_histogram2d_ysigma (const gsl_histogram2d * h)
{
  const double ymean = gsl_histogram2d_ymean (h);
  const size_t nx = h->nx;
  const size_t ny = h->ny;
  size_t i;
  size_t j;

  /* Compute the bin-weighted arithmetic mean M of a histogram using the
     recurrence relation

     M(n) = M(n-1) + (x[n] - M(n-1)) (w(n)/(W(n-1) + w(n)))
     W(n) = W(n-1) + w(n)

   */

  long double wvariance = 0;
  long double W = 0;

  for (j = 0; j < ny; j++)
    {
      double yj = (h->yrange[j + 1] + h->yrange[j]) / 2.0 - ymean;
      double wj = 0;

      for (i = 0; i < nx; i++)
        {
          double wij = h->bin[i * ny + j];
          if (wij > 0)
            wj += wij;
        }
      if (wj > 0)
        {
          W += wj;
          wvariance += ((yj * yj) - wvariance) * (wj / W);
        }
    }

  {
    double ysigma = sqrt (wvariance);
    return ysigma;
  }
}

double
gsl_histogram2d_cov (const gsl_histogram2d * h)
{
  const double xmean = gsl_histogram2d_xmean (h);
  const double ymean = gsl_histogram2d_ymean (h);
  const size_t nx = h->nx;
  const size_t ny = h->ny;
  size_t i;
  size_t j;

  /* Compute the bin-weighted arithmetic mean M of a histogram using the
     recurrence relation

     M(n) = M(n-1) + (x[n] - M(n-1)) (w(n)/(W(n-1) + w(n)))
     W(n) = W(n-1) + w(n)

   */

  long double wcovariance = 0;
  long double W = 0;

  for (j = 0; j < ny; j++)
    {
      for (i = 0; i < nx; i++)
        {
          double xi = (h->xrange[i + 1] + h->xrange[i]) / 2.0 - xmean;
          double yj = (h->yrange[j + 1] + h->yrange[j]) / 2.0 - ymean;
          double wij = h->bin[i * ny + j];

          if (wij > 0)
            {
              W += wij;
              wcovariance += ((xi * yj) - wcovariance) * (wij / W);
            }
        }
    }

  return wcovariance;

}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/*
/
/ adkGSL.h
/
/ homebrewed addons to gsl */

#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_sort_vector.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <math.h>
#include <assert.h>

int gsl_vector_includes(gsl_vector *aVec, double aValue);
double gsl_vector_sum(gsl_vector *aVec, int aVecSize);
double gsl_matrix_row_sum(gsl_matrix *mat, int row, int rowSize);
double gsl_matrix_row_max(gsl_matrix *mat, int iRow);

double log_sum(gsl_vector *vec);
double  *gsl_matrix_2_lapack(gsl_matrix *m);
double gsl_vector_dot_product(gsl_vector *vec1, gsl_vector *vec2);
void gsl_vector_outer_product(gsl_vector *vec1, gsl_vector *vec2, gsl_matrix *result);

gsl_matrix *gsl_matrix_power(gsl_matrix *A, int k);
gsl_matrix *gsl_matrix_power_logs(gsl_matrix *A, int k);
void gsl_matrix_lower_tri_copy(gsl_matrix *src, gsl_matrix *dest);
void gsl_matrix_upper_tri_copy(gsl_matrix *src, gsl_matrix *dest);
void gsl_matrix_lower_tri(gsl_matrix *x);
void gsl_matrix_upper_tri(gsl_matrix *x);
double gsl_matrix_trace(gsl_matrix *x);
void gsl_matrix_covariance(gsl_matrix *data, gsl_matrix *cov);
void gsl_matrix_bootstrap(gsl_matrix *orig, gsl_matrix *boot, gsl_rng *rng);

void gsl_matrix_prettyPrint(gsl_matrix *m);

void fillMatrixFromArray(double *numbers, gsl_matrix *dest, int nrow, int ncol);
void fillMatrixFromVector(gsl_vector *numbers, gsl_matrix *dest, int nrow, int ncol);
void fillArrayFromMatrix(gsl_matrix *src, gsl_vector *dest);

void fillCholArrayFromMatrix(gsl_matrix *src, gsl_vector *dest);
void fillMatrixFromCholArray(double *numbers, gsl_matrix *dest, int nrow, int ncol);

int gsl_matrix_invert_lapack(gsl_matrix *m);

void gsl_matrix_prettyPrint(gsl_matrix *m);
double gsl_matrix_sum(gsl_matrix *mat);

//lapack wrappers
//int dgeev(char jobvl, char jobvr, int n,  double *a, int lda, double *wr, double *wi, double *vl,
//	  int ldvl,double *vr, int ldvr, double *work, int lwork);
//int dhseqr(char job, char compz, int n, int ilo, int ihi,double *h, int ldh, double *wr, double *wi, double *z, int ldz, double *work, int lwork);

max_stars_repo_name: dortegadelv/IMaDNA

max_issues_repo_name: dortegadelv/IMaDNA


text:
/* specfunc/bessel_amp_phase.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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

/* Author:  G. Jungman */

#include <config.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_bessel.h>
#include "bessel_amp_phase.h"

/* chebyshev expansions for amplitude and phase
   functions used in bessel evaluations

   These are the same for J0,Y0 and for J1,Y1, so
   they sit outside those functions.
*/

static double bm0_data[21] = {
   0.09284961637381644,
  -0.00142987707403484,
   0.00002830579271257,
  -0.00000143300611424,
   0.00000012028628046,
  -0.00000001397113013,
   0.00000000204076188,
  -0.00000000035399669,
   0.00000000007024759,
  -0.00000000001554107,
   0.00000000000376226,
  -0.00000000000098282,
   0.00000000000027408,
  -0.00000000000008091,
   0.00000000000002511,
  -0.00000000000000814,
   0.00000000000000275,
  -0.00000000000000096,
   0.00000000000000034,
  -0.00000000000000012,
   0.00000000000000004
};
const cheb_series _gsl_sf_bessel_amp_phase_bm0_cs = {
  bm0_data,
  20,
  -1, 1,
  10
};

static double bth0_data[24] = {
  -0.24639163774300119,
   0.001737098307508963,
  -0.000062183633402968,
   0.000004368050165742,
  -0.000000456093019869,
   0.000000062197400101,
  -0.000000010300442889,
   0.000000001979526776,
  -0.000000000428198396,
   0.000000000102035840,
  -0.000000000026363898,
   0.000000000007297935,
  -0.000000000002144188,
   0.000000000000663693,
  -0.000000000000215126,
   0.000000000000072659,
  -0.000000000000025465,
   0.000000000000009229,
  -0.000000000000003448,
   0.000000000000001325,
  -0.000000000000000522,
   0.000000000000000210,
  -0.000000000000000087,
   0.000000000000000036
};
const cheb_series _gsl_sf_bessel_amp_phase_bth0_cs = {
  bth0_data,
  23,
  -1, 1,
  12
};

static double bm1_data[21] = {
   0.1047362510931285,
   0.00442443893702345,
  -0.00005661639504035,
   0.00000231349417339,
  -0.00000017377182007,
   0.00000001893209930,
  -0.00000000265416023,
   0.00000000044740209,
  -0.00000000008691795,
   0.00000000001891492,
  -0.00000000000451884,
   0.00000000000116765,
  -0.00000000000032265,
   0.00000000000009450,
  -0.00000000000002913,
   0.00000000000000939,
  -0.00000000000000315,
   0.00000000000000109,
  -0.00000000000000039,
   0.00000000000000014,
  -0.00000000000000005,
};
const cheb_series _gsl_sf_bessel_amp_phase_bm1_cs = {
  bm1_data,
  20,
  -1, 1,
  10
};

static double bth1_data[24] = {
   0.74060141026313850,
  -0.004571755659637690,
   0.000119818510964326,
  -0.000006964561891648,
   0.000000655495621447,
  -0.000000084066228945,
   0.000000013376886564,
  -0.000000002499565654,
   0.000000000529495100,
  -0.000000000124135944,
   0.000000000031656485,
  -0.000000000008668640,
   0.000000000002523758,
  -0.000000000000775085,
   0.000000000000249527,
  -0.000000000000083773,
   0.000000000000029205,
  -0.000000000000010534,
   0.000000000000003919,
  -0.000000000000001500,
   0.000000000000000589,
  -0.000000000000000237,
   0.000000000000000097,
  -0.000000000000000040,
};
const cheb_series _gsl_sf_bessel_amp_phase_bth1_cs = {
  bth1_data,
  23,
  -1, 1,
  12
};

int
gsl_sf_bessel_asymp_Mnu_e(const double nu, const double x, double * result)
{
  const double r  = 2.0*nu/x;
  const double r2 = r*r;
  const double x2 = x*x;
  const double term1 = (r2-1.0/x2)/8.0;
  const double term2 = (r2-1.0/x2)*(r2-9.0/x2)*3.0/128.0;
  const double Mnu2_c = 2.0/(M_PI) * (1.0 + term1 + term2);
  *result = sqrt(Mnu2_c)/sqrt(x); /* will never underflow this way */
  return GSL_SUCCESS;
}

int
gsl_sf_bessel_asymp_thetanu_corr_e(const double nu, const double x, double * result)
{
  const double r  = 2.0*nu/x;
  const double r2 = r*r;
  const double x2 = x*x;
  const double term1 = x*(r2 - 1.0/x2)/8.0;
  const double term2 = x*(r2 - 1.0/x2)*(r2 - 25.0/x2)/384.0;
  *result = (-0.25*M_PI + term1 + term2);
  return GSL_SUCCESS;
}


text:
#pragma once

// Standard
#include <exception>
#include <cassert>
#include <string>
#include <iostream>
#include <sstream>
#include <fstream>
#include <memory>
#include <vector>
#include <map>
#include <stack>
#include <cstdint>
#include <iomanip>
#include <codecvt>
#include <algorithm>
#include <functional>
#include <limits>
#include <filesystem>
#include <tuple>
#include <iterator>
#include <random>
#include <cmath>

#if defined(DEBUG) || defined(_DEBUG)
#define _CRTDBG_MAP_ALLOC
#include <stdlib.h>
#include <crtdbg.h>
#endif

// Guidelines Support Library
#include <gsl\gsl>

// Windows
#include <windows.h>
#include <winrt\Windows.Foundation.h>

// DirectX
#include <d3d11_3.h>
#include <dxgi1_2.h>
#include <DirectXMath.h>
#include <DirectXPackedVector.h>
#include <DirectXColors.h>
#include <DirectXTK\DDSTextureLoader.h>
#include <DirectXTK\WICTextureLoader.h>
#include <DirectXTK\SpriteBatch.h>
#include <DirectXTK\SpriteFont.h>
#include <DirectXTK\GamePad.h>
#include <DirectXTK\Keyboard.h>
#include <DirectXTK\Mouse.h>


text:
// @file kn2row_conv.c
//
//  \date Created on: Sep 23, 2017
//  \author Gopalakrishna Hegde
//
//   Description:
//
//
//
#include <assert.h>
#include <stdbool.h>
#include <stdlib.h>
#include <string.h>
#include <cblas.h>
#include "common_types.h"
#include "data_reshape.h"
#include "utils.h"

//
// col_shift : +ve --> shift left overlap mat , -ve --> shift right overlap mat
// or shift left base mat and keep overlap mat as it is.
//
//
// row_shift : +ve (coeff is down the center coeff) --> shift up overlap mat ,
// -ve --> shift down overlap mat or shift up the base mat.
void MatrixShiftAdd(float *base_mat,
                     int base_no_rows, int base_no_cols,
                     float *overlap_mat,
                     int ov_no_rows, int ov_no_cols,
                     int row_shift, int col_shift) {
  if (row_shift == 0 && col_shift == 0 && (base_no_rows == ov_no_rows) &&
      (base_no_cols == ov_no_cols)) {
    // normal matrix add
    cblas_saxpy(base_no_rows * base_no_cols, 1.0, overlap_mat, 1, base_mat, 1);
    return;
  }
  int rows_to_add, cols_to_add;
  int base_row_start, base_col_start;
  int ov_row_start, ov_col_start;
  // without padding case
  if (ov_no_rows > base_no_rows) {
    rows_to_add = base_no_rows;
    cols_to_add = base_no_cols;
    base_row_start = 0;
    base_col_start = 0;
    ov_row_start = row_shift < 0? -row_shift : 0;
    ov_col_start = col_shift < 0? -col_shift : 0;

  } else {
    rows_to_add = ov_no_rows - abs(row_shift);
    cols_to_add = ov_no_cols - abs(col_shift);

    ov_col_start = col_shift > 0? col_shift : 0;
    ov_row_start = row_shift > 0? row_shift : 0;
    base_row_start = row_shift < 0? -row_shift : 0;
    base_col_start = col_shift < 0? -col_shift : 0;
  }

  for (int r = 0; r < rows_to_add; ++r) {
    int base_mat_offset = (r + base_row_start) * base_no_cols + base_col_start;
    int overlap_mat_offset = (r + ov_row_start) * ov_no_cols + ov_col_start;
    cblas_saxpy(cols_to_add, 1.0, overlap_mat + overlap_mat_offset, 1,
                base_mat + base_mat_offset, 1);
  }
}

/* Ker2Row convolution implementations.
 *
 * Assumptions:
 * 1. in_data is in NCHW format.
 * 2. filters are in MCKK format where M is the no of output maps.
 * 3. Stride will always be 1.
 * 4. pad will be zero or kernel_size / 2
 *
 * Output will be in NCHW format.
 */
bool Kn2RowConvLayer(const float *in_data, const float *filters,
                         const float *bias, TensorDim in_dim,
                         TensorDim filt_dim, int stride, int pad, int group,
                         float *output) {
  // Currently we have limited support.
  assert(group == 1);
  assert((pad == 0) || (pad == filt_dim.w / 2));
  assert(in_dim.n == 1);
  assert(filt_dim.h == filt_dim.w);
  assert(stride == 1);

  // Output dimensions.
  TensorDim out_dim;
  out_dim.w = (in_dim.w + (pad + pad) - filt_dim.w) / stride + 1;
  out_dim.h = (in_dim.h + (pad + pad) - filt_dim.h) / stride + 1;
  out_dim.c = filt_dim.n;
  out_dim.n = in_dim.n;

  // Re-arrange filters in the  k x k x no_out_maps x no_in_maps.
  // We can avoid this if the filters are already reshaped in this format.
  float *kkmc_filters = malloc(filt_dim.n * filt_dim.c * filt_dim.h *
                               filt_dim.w * sizeof(float));
  NCHW2HWNC(filters, filt_dim.n, filt_dim.c, filt_dim.h, filt_dim.w,
            kkmc_filters);

  // Just for convenience
  int H = in_dim.h;
  int W = in_dim.w;
  float alpha = 1.0;
  float beta = 0.0;

  // We need separate buffer because GEMM output will have width = H*W even
  // if there is no padding (pad = 0).
  float *gemm_output = malloc(out_dim.c * H * W * sizeof(float));

  // Prefill output buffer with bias if present else set to zero.
  if (bias) {
    for (int m = 0; m < out_dim.c; ++m) {
      for (int a = 0; a < out_dim.h * out_dim.w; ++a) {
        output[m * out_dim.h * out_dim.w + a] = bias[m];
      }
      // For batch size > 1
      for (int b = 1; b < out_dim.n; ++b) {
        memcpy(output + b * out_dim.c * out_dim.h * out_dim.w,
               output, out_dim.c * out_dim.h * out_dim.w * sizeof(float));
      }
    }
  } else {
    memset(output, 0, out_dim.n * out_dim.c * out_dim.h * out_dim.w *
           sizeof(float));
  }

  for (int kr = 0; kr < filt_dim.h; kr++) {
    int row_shift = kr - filt_dim.h / 2;
    for (int kc = 0; kc < filt_dim.w; kc++) {
      int group_no = kr * filt_dim.w + kc;
      int col_shift = kc - filt_dim.w / 2;
      // Matrix dimensions - A -> mxk B -> kxn  C --> mxn
      int m = filt_dim.n;
      int k = filt_dim.c;
      int n = in_dim.h * in_dim.w;
      // This is just 1x1 convolution
      cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
                  m, n, k, alpha, kkmc_filters + group_no * m * k,
                  k, in_data, n, beta, gemm_output, n);
      // Slide the resulting matrix which has contribution from one of the
      // KxK kernel coefficients and add to the output.
      for (int omap = 0; omap < filt_dim.n; omap++) {
        MatrixShiftAdd(output + omap * out_dim.h * out_dim.w,
                        out_dim.h, out_dim.w,
                        gemm_output + omap * H * W,
                        H, W, row_shift, col_shift);
      }
    }
  }
  free(kkmc_filters);
  free(gemm_output);
  return true;
}

max_stars_repo_name: gplhegde/convolution-flavors

max_issues_repo_name: chayitw/convolution-flavors

max_forks_repo_name: chayitw/convolution-flavors


text:
/***************************************************************************
 *   Copyright (C) 2008 by Mikhail Zaslavskiy   *
 *   mikhail.zaslavskiy@ensmp.fr   *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   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.             *
 ***************************************************************************/
#ifndef ALGORITHM_H
#define ALGORITHM_H
#define EPSILON 1e-100
#include "rpc.h"
#include "graph.h"
#include <math.h>
#include "hungarian.h"

#include <gsl/gsl_blas.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_eigen.h>

#include <vector>
#include <iostream>

/**
Class of graph matching results
*/
class match_result
{
public:
	match_result(){gm_P=NULL;gm_P_exact=NULL;salgo="";}
	std::vector<double> vd_trace;
	int inum_iteration;
	double dres;
	gsl_matrix* gm_P;
	gsl_matrix* gm_P_exact;
	double dtime;
	double dfvalue;
	double dfvalue_exact;
	std::string salgo;
	~match_result(){};

};
/**
Parent class for all graph matching algorithms

	@author Mikhail Zaslavskiy <mikhail.zaslavskiy@ensmp.fr>
*/
class algorithm : public rpc
{
public:
     algorithm(std::string );
     algorithm();
     match_result gmatch(graph& g, graph& h,gsl_matrix* gm_P_i=NULL, gsl_matrix* gm_ldh=NULL,double dalpha_ldh=-1);//common stuff,
     virtual match_result match(graph& g, graph& h, gsl_matrix* gm_P_i=NULL, gsl_matrix* gm_ldh=NULL,double dalpha_ldh=-1)=0;//particular method implementation
     double graph_dist(graph &g,graph &h,gsl_matrix* gm_P,char cscore_matrix);
     double graph_dist(graph &g, graph &h,char cscore_matrix);
      ~algorithm();

     const gsl_matrix* get_ldhmatrix(){return gm_ldh;};
     void  set_ldhmatrix(const gsl_matrix* _gm_A);

protected:
    gsl_matrix *gm_ldh;
    double dalpha_ldh;
    void update_C_hungarian(gsl_matrix* gm_C,int isign=1, bool bback=false);
    double f_qcv(gsl_matrix *gm_Ag_d,gsl_matrix *gm_Ah_d,gsl_matrix* gm_P,gsl_matrix * gm_temp,bool bqcv=false);
    char cdesc_matrix,cscore_matrix;
    parameter pdebug,pdebug_f;
    bool bverbose;
    std::string sverbfile;
    std::ofstream fverbose;
    long long N;
    double df_norm;
    bool bnosymm;

};

#endif