text:
/* linalg/gsl_linalg.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2006, 2007 Gerard Jungman, Brian Gough, 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.
 */

#ifndef __GSL_LINALG_H__
#define __GSL_LINALG_H__

#include <gsl/gsl_mode.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>

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

__BEGIN_DECLS

typedef enum
  {
    GSL_LINALG_MOD_NONE = 0,
    GSL_LINALG_MOD_TRANSPOSE = 1,
    GSL_LINALG_MOD_CONJUGATE = 2
  }
gsl_linalg_matrix_mod_t;

/* Note: You can now use the gsl_blas_dgemm function instead of matmult */

/* Simple implementation of matrix multiply.
 * Calculates C = A.B
 *
 * exceptions: GSL_EBADLEN
 */
int gsl_linalg_matmult (const gsl_matrix * A,
                        const gsl_matrix * B,
                        gsl_matrix * C);

/* Simple implementation of matrix multiply.
 * Allows transposition of either matrix, so it
 * can compute A.B or Trans(A).B or A.Trans(B) or Trans(A).Trans(B)
 *
 * exceptions: GSL_EBADLEN
 */
int gsl_linalg_matmult_mod (const gsl_matrix * A,
                            gsl_linalg_matrix_mod_t modA,
                            const gsl_matrix * B,
                            gsl_linalg_matrix_mod_t modB,
                            gsl_matrix * C);

/* Calculate the matrix exponential by the scaling and
 * squaring method described in Moler + Van Loan,
 * SIAM Rev 20, 801 (1978). The mode argument allows
 * choosing an optimal strategy, from the table
 * given in the paper, for a given precision.
 *
 * exceptions: GSL_ENOTSQR, GSL_EBADLEN
 */
int gsl_linalg_exponential_ss(
  const gsl_matrix * A,
  gsl_matrix * eA,
  gsl_mode_t mode
  );

/* Householder Transformations */

double gsl_linalg_householder_transform (gsl_vector * v);
gsl_complex gsl_linalg_complex_householder_transform (gsl_vector_complex * v);

int gsl_linalg_householder_hm (double tau,
                               const gsl_vector * v,
                               gsl_matrix * A);

int gsl_linalg_householder_mh (double tau,
                               const gsl_vector * v,
                               gsl_matrix * A);

int gsl_linalg_householder_hv (double tau,
                               const gsl_vector * v,
                               gsl_vector * w);

int gsl_linalg_householder_hm1 (double tau,
                                gsl_matrix * A);

int gsl_linalg_complex_householder_hm (gsl_complex tau,
                                       const gsl_vector_complex * v,
                                       gsl_matrix_complex * A);

int gsl_linalg_complex_householder_mh (gsl_complex tau,
                                       const gsl_vector_complex * v,
                                       gsl_matrix_complex * A);

int gsl_linalg_complex_householder_hv (gsl_complex tau,
                                       const gsl_vector_complex * v,
                                       gsl_vector_complex * w);

/* Hessenberg reduction */

int gsl_linalg_hessenberg_decomp(gsl_matrix *A, gsl_vector *tau);
int gsl_linalg_hessenberg_unpack(gsl_matrix * H, gsl_vector * tau,
                                 gsl_matrix * U);
int gsl_linalg_hessenberg_unpack_accum(gsl_matrix * H, gsl_vector * tau,
                                       gsl_matrix * U);
int gsl_linalg_hessenberg_set_zero(gsl_matrix * H);
int gsl_linalg_hessenberg_submatrix(gsl_matrix *M, gsl_matrix *A,
                                    size_t top, gsl_vector *tau);

/* To support gsl-1.9 interface: DEPRECATED */
int gsl_linalg_hessenberg(gsl_matrix *A, gsl_vector *tau);

/* Hessenberg-Triangular reduction */

int gsl_linalg_hesstri_decomp(gsl_matrix * A, gsl_matrix * B,
                              gsl_matrix * U, gsl_matrix * V,
                              gsl_vector * work);

/* Singular Value Decomposition

 * exceptions:
 */

int
gsl_linalg_SV_decomp (gsl_matrix * A,
                      gsl_matrix * V,
                      gsl_vector * S,
                      gsl_vector * work);

int
gsl_linalg_SV_decomp_mod (gsl_matrix * A,
                          gsl_matrix * X,
                          gsl_matrix * V,
                          gsl_vector * S,
                          gsl_vector * work);

int gsl_linalg_SV_decomp_jacobi (gsl_matrix * A,
                                 gsl_matrix * Q,
                                 gsl_vector * S);

int
gsl_linalg_SV_solve (const gsl_matrix * U,
                     const gsl_matrix * Q,
                     const gsl_vector * S,
                     const gsl_vector * b,
                     gsl_vector * x);

int gsl_linalg_SV_leverage(const gsl_matrix *U, gsl_vector *h);

/* LU Decomposition, Gaussian elimination with partial pivoting
 */

int gsl_linalg_LU_decomp (gsl_matrix * A, gsl_permutation * p, int *signum);

int gsl_linalg_LU_solve (const gsl_matrix * LU,
                         const gsl_permutation * p,
                         const gsl_vector * b,
                         gsl_vector * x);

int gsl_linalg_LU_svx (const gsl_matrix * LU,
                       const gsl_permutation * p,
                       gsl_vector * x);

int gsl_linalg_LU_refine (const gsl_matrix * A,
                          const gsl_matrix * LU,
                          const gsl_permutation * p,
                          const gsl_vector * b,
                          gsl_vector * x,
                          gsl_vector * residual);

int gsl_linalg_LU_invert (const gsl_matrix * LU,
                          const gsl_permutation * p,
                          gsl_matrix * inverse);

double gsl_linalg_LU_det (gsl_matrix * LU, int signum);
double gsl_linalg_LU_lndet (gsl_matrix * LU);
int gsl_linalg_LU_sgndet (gsl_matrix * lu, int signum);

/* Complex LU Decomposition */

int gsl_linalg_complex_LU_decomp (gsl_matrix_complex * A,
                                  gsl_permutation * p,
                                  int *signum);

int gsl_linalg_complex_LU_solve (const gsl_matrix_complex * LU,
                                 const gsl_permutation * p,
                                 const gsl_vector_complex * b,
                                 gsl_vector_complex * x);

int gsl_linalg_complex_LU_svx (const gsl_matrix_complex * LU,
                               const gsl_permutation * p,
                               gsl_vector_complex * x);

int gsl_linalg_complex_LU_refine (const gsl_matrix_complex * A,
                                  const gsl_matrix_complex * LU,
                                  const gsl_permutation * p,
                                  const gsl_vector_complex * b,
                                  gsl_vector_complex * x,
                                  gsl_vector_complex * residual);

int gsl_linalg_complex_LU_invert (const gsl_matrix_complex * LU,
                                  const gsl_permutation * p,
                                  gsl_matrix_complex * inverse);

gsl_complex gsl_linalg_complex_LU_det (gsl_matrix_complex * LU,
                                       int signum);

double gsl_linalg_complex_LU_lndet (gsl_matrix_complex * LU);

gsl_complex gsl_linalg_complex_LU_sgndet (gsl_matrix_complex * LU,
                                          int signum);

/* QR decomposition */

int gsl_linalg_QR_decomp (gsl_matrix * A,
                          gsl_vector * tau);

int gsl_linalg_QR_solve (const gsl_matrix * QR,
                         const gsl_vector * tau,
                         const gsl_vector * b,
                         gsl_vector * x);

int gsl_linalg_QR_svx (const gsl_matrix * QR,
                       const gsl_vector * tau,
                       gsl_vector * x);

int gsl_linalg_QR_lssolve (const gsl_matrix * QR,
                           const gsl_vector * tau,
                           const gsl_vector * b,
                           gsl_vector * x,
                           gsl_vector * residual);

int gsl_linalg_QR_QRsolve (gsl_matrix * Q,
                           gsl_matrix * R,
                           const gsl_vector * b,
                           gsl_vector * x);

int gsl_linalg_QR_Rsolve (const gsl_matrix * QR,
                          const gsl_vector * b,
                          gsl_vector * x);

int gsl_linalg_QR_Rsvx (const gsl_matrix * QR,
                        gsl_vector * x);

int gsl_linalg_QR_update (gsl_matrix * Q,
                          gsl_matrix * R,
                          gsl_vector * w,
                          const gsl_vector * v);

int gsl_linalg_QR_QTvec (const gsl_matrix * QR,
                         const gsl_vector * tau,
                         gsl_vector * v);

int gsl_linalg_QR_Qvec (const gsl_matrix * QR,
                        const gsl_vector * tau,
                        gsl_vector * v);

int gsl_linalg_QR_QTmat (const gsl_matrix * QR,
                         const gsl_vector * tau,
                         gsl_matrix * A);

int gsl_linalg_QR_unpack (const gsl_matrix * QR,
                          const gsl_vector * tau,
                          gsl_matrix * Q,
                          gsl_matrix * R);

int gsl_linalg_R_solve (const gsl_matrix * R,
                        const gsl_vector * b,
                        gsl_vector * x);

int gsl_linalg_R_svx (const gsl_matrix * R,
                      gsl_vector * x);

/* Q R P^T decomposition */

int gsl_linalg_QRPT_decomp (gsl_matrix * A,
                            gsl_vector * tau,
                            gsl_permutation * p,
                            int *signum,
                            gsl_vector * norm);

int gsl_linalg_QRPT_decomp2 (const gsl_matrix * A,
                             gsl_matrix * q, gsl_matrix * r,
                             gsl_vector * tau,
                             gsl_permutation * p,
                             int *signum,
                             gsl_vector * norm);

int gsl_linalg_QRPT_solve (const gsl_matrix * QR,
                           const gsl_vector * tau,
                           const gsl_permutation * p,
                           const gsl_vector * b,
                           gsl_vector * x);

int gsl_linalg_QRPT_svx (const gsl_matrix * QR,
                         const gsl_vector * tau,
                         const gsl_permutation * p,
                         gsl_vector * x);

int gsl_linalg_QRPT_QRsolve (const gsl_matrix * Q,
                             const gsl_matrix * R,
                             const gsl_permutation * p,
                             const gsl_vector * b,
                             gsl_vector * x);

int gsl_linalg_QRPT_Rsolve (const gsl_matrix * QR,
                             const gsl_permutation * p,
                             const gsl_vector * b,
                             gsl_vector * x);

int gsl_linalg_QRPT_Rsvx (const gsl_matrix * QR,
                           const gsl_permutation * p,
                           gsl_vector * x);

int gsl_linalg_QRPT_update (gsl_matrix * Q,
                            gsl_matrix * R,
                            const gsl_permutation * p,
                            gsl_vector * u,
                            const gsl_vector * v);

/* LQ decomposition */

int gsl_linalg_LQ_decomp (gsl_matrix * A, gsl_vector * tau);

int gsl_linalg_LQ_solve_T (const gsl_matrix * LQ, const gsl_vector * tau,
			 const gsl_vector * b, gsl_vector * x);

int gsl_linalg_LQ_svx_T (const gsl_matrix * LQ, const gsl_vector * tau,
                         gsl_vector * x);

int gsl_linalg_LQ_lssolve_T (const gsl_matrix * LQ, const gsl_vector * tau,
			   const gsl_vector * b, gsl_vector * x,
			   gsl_vector * residual);

int gsl_linalg_LQ_Lsolve_T (const gsl_matrix * LQ, const gsl_vector * b,
			  gsl_vector * x);

int gsl_linalg_LQ_Lsvx_T (const gsl_matrix * LQ, gsl_vector * x);

int gsl_linalg_L_solve_T (const gsl_matrix * L, const gsl_vector * b,
			gsl_vector * x);

int gsl_linalg_LQ_vecQ (const gsl_matrix * LQ, const gsl_vector * tau,
			gsl_vector * v);

int gsl_linalg_LQ_vecQT (const gsl_matrix * LQ, const gsl_vector * tau,
			 gsl_vector * v);

int gsl_linalg_LQ_unpack (const gsl_matrix * LQ, const gsl_vector * tau,
			  gsl_matrix * Q, gsl_matrix * L);

int gsl_linalg_LQ_update (gsl_matrix * Q, gsl_matrix * R,
			  const gsl_vector * v, gsl_vector * w);
int gsl_linalg_LQ_LQsolve (gsl_matrix * Q, gsl_matrix * L,
			   const gsl_vector * b, gsl_vector * x);

/* P^T L Q decomposition */

int gsl_linalg_PTLQ_decomp (gsl_matrix * A, gsl_vector * tau,
			    gsl_permutation * p, int *signum,
			    gsl_vector * norm);

int gsl_linalg_PTLQ_decomp2 (const gsl_matrix * A, gsl_matrix * q,
			     gsl_matrix * r, gsl_vector * tau,
			     gsl_permutation * p, int *signum,
			     gsl_vector * norm);

int gsl_linalg_PTLQ_solve_T (const gsl_matrix * QR,
			   const gsl_vector * tau,
			   const gsl_permutation * p,
			   const gsl_vector * b,
			   gsl_vector * x);

int gsl_linalg_PTLQ_svx_T (const gsl_matrix * LQ,
                           const gsl_vector * tau,
                           const gsl_permutation * p,
                           gsl_vector * x);

int gsl_linalg_PTLQ_LQsolve_T (const gsl_matrix * Q, const gsl_matrix * L,
			     const gsl_permutation * p,
			     const gsl_vector * b,
			     gsl_vector * x);

int gsl_linalg_PTLQ_Lsolve_T (const gsl_matrix * LQ,
			    const gsl_permutation * p,
			    const gsl_vector * b,
			    gsl_vector * x);

int gsl_linalg_PTLQ_Lsvx_T (const gsl_matrix * LQ,
			  const gsl_permutation * p,
			  gsl_vector * x);

int gsl_linalg_PTLQ_update (gsl_matrix * Q, gsl_matrix * L,
			    const gsl_permutation * p,
			    const gsl_vector * v, gsl_vector * w);

/* Cholesky Decomposition */

int gsl_linalg_cholesky_decomp (gsl_matrix * A);

int gsl_linalg_cholesky_solve (const gsl_matrix * cholesky,
                               const gsl_vector * b,
                               gsl_vector * x);

int gsl_linalg_cholesky_svx (const gsl_matrix * cholesky,
                             gsl_vector * x);

int gsl_linalg_cholesky_invert(gsl_matrix * cholesky);

/* Cholesky decomposition with unit-diagonal triangular parts.
 *   A = L D L^T, where diag(L) = (1,1,...,1).
 *   Upon exit, A contains L and L^T as for Cholesky, and
 *   the diagonal of A is (1,1,...,1). The vector Dis set
 *   to the diagonal elements of the diagonal matrix D.
 */
int gsl_linalg_cholesky_decomp_unit(gsl_matrix * A, gsl_vector * D);

/* Complex Cholesky Decomposition */

int gsl_linalg_complex_cholesky_decomp (gsl_matrix_complex * A);

int gsl_linalg_complex_cholesky_solve (const gsl_matrix_complex * cholesky,
                                       const gsl_vector_complex * b,
                                       gsl_vector_complex * x);

int gsl_linalg_complex_cholesky_svx (const gsl_matrix_complex * cholesky,
                                     gsl_vector_complex * x);

int gsl_linalg_complex_cholesky_invert(gsl_matrix_complex * cholesky);

/* Symmetric to symmetric tridiagonal decomposition */

int gsl_linalg_symmtd_decomp (gsl_matrix * A,
                              gsl_vector * tau);

int gsl_linalg_symmtd_unpack (const gsl_matrix * A,
                              const gsl_vector * tau,
                              gsl_matrix * Q,
                              gsl_vector * diag,
                              gsl_vector * subdiag);

int gsl_linalg_symmtd_unpack_T (const gsl_matrix * A,
                                gsl_vector * diag,
                                gsl_vector * subdiag);

/* Hermitian to symmetric tridiagonal decomposition */

int gsl_linalg_hermtd_decomp (gsl_matrix_complex * A,
                              gsl_vector_complex * tau);

int gsl_linalg_hermtd_unpack (const gsl_matrix_complex * A,
                              const gsl_vector_complex * tau,
                              gsl_matrix_complex * U,
                              gsl_vector * diag,
                              gsl_vector * sudiag);

int gsl_linalg_hermtd_unpack_T (const gsl_matrix_complex * A,
                                gsl_vector * diag,
                                gsl_vector * subdiag);

/* Linear Solve Using Householder Transformations

 * exceptions:
 */

int gsl_linalg_HH_solve (gsl_matrix * A, const gsl_vector * b, gsl_vector * x);
int gsl_linalg_HH_svx (gsl_matrix * A, gsl_vector * x);

/* Linear solve for a symmetric tridiagonal system.

 * The input vectors represent the NxN matrix as follows:
 *
 *     diag[0]  offdiag[0]             0    ...
 *  offdiag[0]     diag[1]    offdiag[1]    ...
 *           0  offdiag[1]       diag[2]    ...
 *           0           0    offdiag[2]    ...
 *         ...         ...           ...    ...
 */
int gsl_linalg_solve_symm_tridiag (const gsl_vector * diag,
                                   const gsl_vector * offdiag,
                                   const gsl_vector * b,
                                   gsl_vector * x);

/* Linear solve for a nonsymmetric tridiagonal system.

 * The input vectors represent the NxN matrix as follows:
 *
 *       diag[0]  abovediag[0]              0    ...
 *  belowdiag[0]       diag[1]   abovediag[1]    ...
 *             0  belowdiag[1]        diag[2]    ...
 *             0             0   belowdiag[2]    ...
 *           ...           ...            ...    ...
 */
int gsl_linalg_solve_tridiag (const gsl_vector * diag,
                                   const gsl_vector * abovediag,
                                   const gsl_vector * belowdiag,
                                   const gsl_vector * b,
                                   gsl_vector * x);

/* Linear solve for a symmetric cyclic tridiagonal system.

 * The input vectors represent the NxN matrix as follows:
 *
 *      diag[0]  offdiag[0]             0   .....  offdiag[N-1]
 *   offdiag[0]     diag[1]    offdiag[1]   .....
 *            0  offdiag[1]       diag[2]   .....
 *            0           0    offdiag[2]   .....
 *          ...         ...
 * offdiag[N-1]         ...
 */
int gsl_linalg_solve_symm_cyc_tridiag (const gsl_vector * diag,
                                       const gsl_vector * offdiag,
                                       const gsl_vector * b,
                                       gsl_vector * x);

/* Linear solve for a nonsymmetric cyclic tridiagonal system.

 * The input vectors represent the NxN matrix as follows:
 *
 *        diag[0]  abovediag[0]             0   .....  belowdiag[N-1]
 *   belowdiag[0]       diag[1]  abovediag[1]   .....
 *              0  belowdiag[1]       diag[2]
 *              0             0  belowdiag[2]   .....
 *            ...           ...
 * abovediag[N-1]           ...
 */
int gsl_linalg_solve_cyc_tridiag (const gsl_vector * diag,
                                  const gsl_vector * abovediag,
                                  const gsl_vector * belowdiag,
                                  const gsl_vector * b,
                                  gsl_vector * x);

/* Bidiagonal decomposition */

int gsl_linalg_bidiag_decomp (gsl_matrix * A,
                              gsl_vector * tau_U,
                              gsl_vector * tau_V);

int gsl_linalg_bidiag_unpack (const gsl_matrix * A,
                              const gsl_vector * tau_U,
                              gsl_matrix * U,
                              const gsl_vector * tau_V,
                              gsl_matrix * V,
                              gsl_vector * diag,
                              gsl_vector * superdiag);

int gsl_linalg_bidiag_unpack2 (gsl_matrix * A,
                               gsl_vector * tau_U,
                               gsl_vector * tau_V,
                               gsl_matrix * V);

int gsl_linalg_bidiag_unpack_B (const gsl_matrix * A,
                                gsl_vector * diag,
                                gsl_vector * superdiag);

/* Balancing */

int gsl_linalg_balance_matrix (gsl_matrix * A, gsl_vector * D);
int gsl_linalg_balance_accum (gsl_matrix * A, gsl_vector * D);
int gsl_linalg_balance_columns (gsl_matrix * A, gsl_vector * D);

__END_DECLS

#endif /* __GSL_LINALG_H__ */

max_stars_repo_licenses: Net-SNMP Xnet

max_issues_repo_licenses: Net-SNMP Xnet

max_forks_repo_licenses: Net-SNMP Xnet


text:
/**
 * Interface to flatfielding
 */

#ifndef _SPC_FLATFIELD_H
#define _SPC_FLATFIELD_H 1

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

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

#include "aXe_grism.h"
#include "aXe_utils.h"

typedef struct
{
  int poly_order;      /* order of the polynomial */
  gsl_matrix **coeffs; /* the list of gsl_matrix'es containing */
                       /* the polynomial coefficients of the FF */
  double wmin;         /* Minumum wavelength value for the coeff. normalization */
  double wmax;         /* Maximum wavelength value for the coeff. normalization */
}
poly_cube_flatfield;

/**
  Descriptor of flatfield with a single flatfield image and a
  polynomial flatfield dependence.
*/
typedef struct
{
  gsl_matrix *flatfield; /* The flatfield image, giving the relative
                            efficiency for each pixel in the aperture */
  gsl_matrix *errors;	 /* relative errors for flatfield */
  int poly_order;	 /* order of the polynomial */
  double *poly_coeffs;	 /* coefficients for the polynom */
}
polynom_flatfield;

/**
   Descriptor of flatfield with multiple flatfields at different
   wave lengths
*/
typedef struct
{
  int num_flats;                /* number of (lambda,flatfield) pairs */
  gsl_matrix **flatfields;      /* flatfield images, giving the relative
				   efficiency for each pixel in the aperture */
  gsl_matrix **errors;	        /* relative errors for flatfields */
  double *lambdas;              /* wave lengths for entries in flatfields */
}
multi_flatfield;

/**
  Descriptor of a flatfielding function, including the function
  to compute the flat field value of a given pixel, the bounding box
  we handle, and auxillary data (usually a flatfield image)
*/
typedef struct s_flatfield
{
  void (*func) (const double lambda, const int x, const int y,
		PIXEL_T * const val, PIXEL_T * const err,
		const struct s_flatfield * const self);
  int ll_x, ll_y, w, h;	/* bounding box we handle */
  union
  {
    polynom_flatfield *poly;
    multi_flatfield *multi;
  }
  data;
}
flatfield_d;

poly_cube_flatfield * load_flat_poly_cube(char *fname);
void free_flat_poly_cube(poly_cube_flatfield *poly_cube);

flatfield_d *make_poly_flatfield (const beam * const curbeam,
				  const char *const flat_name,
				  const int order,
				  const double *const coeffs);
void free_poly_flatfield (flatfield_d * const flat);
flatfield_d *make_multi_flatfield (const beam * const curbeam,
				   const int num_flats, const double lambda,
				   const char *const flat_name, ...);
void free_multi_flatfield (flatfield_d * const flat);

void apply_flatfield (ap_pixel * const ap_p, const flatfield_d * const flat);
double poly_cube_flatfield_lambda (const double lambda, const int x, const int y,
				   poly_cube_flatfield *poly_cube);

#endif

ext: h

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe


text:
#include <stdio.h>
#define _USE_MATH_DEFINES
#include <math.h>
#include <gsl/gsl_sf.h>

#define HBARC 197.3269804 // MeV•fm
#define MASS_TRITON 2808.92113298 // MeV
#define MASS_DEUTERON 1875.61294257 // MeV
#define MASS_NEUTRON 939.565420 // MeV
#define MASS_ALPHA 3727.37940 // MeV
const double MU_D = MASS_TRITON*MASS_DEUTERON/(MASS_TRITON + MASS_DEUTERON);
const double MU_N = MASS_ALPHA*MASS_NEUTRON/(MASS_ALPHA + MASS_NEUTRON);
const double ALPHA = 1.0/137.035999084;
const double Q = 17.6;
const double J = 1.5;
const double J1 = 0.5;
const double J2 = 1.0;

double relative_momentum(double e, double mu)
{
  return sqrt(2*mu*e)/HBARC;
}

double sommerfeld(double e, double mu)
{
  return ALPHA * sqrt(mu/(2*e));
}

void coulomb_functions(double* F, double* Fp, double* G, double* Gp, int l, double eta,
    double x)
{
  int info;
  double exp_f, exp_g;
  gsl_sf_result f, fp, g, gp;
  info = gsl_sf_coulomb_wave_FG_e(eta, x, l, 0, &f, &fp, &g, &gp, &exp_f, &exp_g);
  // printf("%.4f  %.4f\n", exp_f, exp_g);
  if (info == GSL_EOVRFLW) printf("overflow or underflow encountered\n");
  *F = f.val*exp(exp_f);
  *Fp = fp.val*exp(exp_f);
  *G = g.val*exp(exp_g);
  *Gp = gp.val*exp(exp_g);
}

double penetration_factor(int l, double eta, double x)
{
  double F, Fp, G, Gp;
  coulomb_functions(&F, &Fp, &G, &Gp, l, eta, x);
  return x / (F*F + G*G);
}

double shift_factor(int l, double eta, double x)
{
  double F, Fp, G, Gp;
  coulomb_functions(&F, &Fp, &G, &Gp, l, eta, x);
  return x * (F*Fp + G*Gp) / (F*F + G*G);
}

double Gamma_c(double gamma_c_sq, int l, double eta, double x)
{
  return 2 * gamma_c_sq * penetration_factor(l, eta, x);
}

double Delta_c(double gamma_c_2, double k, double r, double Bc, int l, double eta)
{
	return -gamma_c_2 * (shift_factor(l, eta, k*r)-Bc);
}

double Sdn2(double energy, double e0, double Bd, double Bn, double gamma_d_2,
    double gamma_n_2, double kd, double ad, double kn, double an, double eta,
    double A)
{
	double Gamma_d = Gamma_c(gamma_d_2, 0, eta, kd*ad);
	double Gamma_n = Gamma_c(gamma_n_2, 2, 0, kn*an);
	double Gamma = Gamma_d + Gamma_n;
	double Delta_d = Delta_c(gamma_d_2, kd, ad, Bd, 0, eta);
	double Delta_n = Delta_c(gamma_n_2, kn, an, Bn, 2, 0);
	double Delta = Delta_d + Delta_n;
  // printf("Sdn2: %.2e  %.2e  %.2e\n", Gamma_d, Gamma_n, Delta);
	return Gamma_d*Gamma_n / ((e0+Delta-energy)*(e0+Delta-energy) + Gamma*Gamma/4) +
	       2.0/M_PI*A * penetration_factor(0, eta, kd*ad) * penetration_factor(0, 0.0, kn*an);
}

/*
Evaluates the cross section (barn = 100fm^2) at energy en (MeV) with R-matrix parameters:
		eB (MeV) where Bc = shift_factor(eB),
		gamma_d_2 (MeV) is the reduced deuteron width,
		gamma_n_2 (MeV) is the reduced neutron width,
		ad (fm) is the deuteron channel radius,
		an (fm) is the neutron channel radius,
and ue (MeV) is the electron screening potential.
*/
double cross_section(double en, double e0, double eB, double gamma_d_2,
		double gamma_n_2, double ad, double an, double ue, double A)
{
	if (en < 0.0 || e0 < 0.0 || eB < 0.0 || gamma_d_2 < 0.0 ||
			gamma_n_2 < 0.0 || ad < 0.0 || an < 0.0 || ue < 0.0) {
		return INFINITY;
	}
	double kd = relative_momentum(en, MU_D); // 1/fm
	double kn = relative_momentum(en+Q, MU_N); // 1/fm
	double kdB = relative_momentum(eB, MU_D); // 1/fm
	double knB = relative_momentum(eB+Q, MU_N); // 1/fm
	double eta = sommerfeld(en, MU_D); // dimensionless
	double etaB = sommerfeld(eB, MU_D); // dimensionless
  // printf("kdB = %.2e\n", kdB);
	double Bd = shift_factor(0, etaB, kdB*ad); // dimensionless
	double Bn = shift_factor(2, 0, knB*an); // dimensionless
  // printf("Bd = %.2e  kdB = %.2e  ad = %.2e\n", Bd, kdB, ad);
	return exp(M_PI*eta*ue/en) *
			M_PI/(kd*kd) * (2.0*J+1.0)/((2.0*J1+1.0)*(2.0*J2+1.0)) *
			Sdn2(en, e0, Bd, Bn, gamma_d_2, gamma_n_2, kd, ad, kn, an, eta, A) / 100;
}

/*
Evaluates the S-factor (MeV•barn = MeV•100fm^2) at energy en (MeV) with R-matrix parameters:
		eB (MeV) where Bc = shift_factor(eB),
		gamma_d_2 (MeV) is the reduced deuteron width,
		gamma_n_2 (MeV) is the reduced neutron width,
		ad (fm) is the deuteron channel radius,
		an (fm) is the neutron channel radius,
and ue (MeV) is the electron screening potential.
*/
double S_factor(double en, double e0, double eB, double gamma_d_2,
		double gamma_n_2, double ad, double an, double ue, double A)
{
	if (en < 0.0 || e0 < 0.0 || eB < 0.0 || gamma_d_2 < 0.0 ||
			gamma_n_2 < 0.0 || ad < 0.0 || an < 0.0 || ue < 0.0) {
		return INFINITY;
	}
	double kd = relative_momentum(en, MU_D); // 1/fm
	double kn = relative_momentum(en+Q, MU_N); // 1/fm
	double kdB = relative_momentum(eB, MU_D); // 1/fm
	double knB = relative_momentum(eB+Q, MU_N); // 1/fm
	double eta = sommerfeld(en, MU_D); // dimensionless
	double etaB = sommerfeld(eB, MU_D); // dimensionless
  // printf("kdB = %.2e\n", kdB);
	double Bd = shift_factor(0, etaB, kdB*ad); // dimensionless
	double Bn = shift_factor(2, 0, knB*an); // dimensionless
  // printf("Bd = %.2e  kdB = %.2e  ad = %.2e\n", Bd, kdB, ad);
	return exp(M_PI*eta*ue/en) * en * exp(2*M_PI*eta) *\
			M_PI/(kd*kd) * (2.0*J+1.0)/((2.0*J1+1.0)*(2.0*J2+1.0)) *\
			Sdn2(en, e0, Bd, Bn, gamma_d_2, gamma_n_2, kd, ad, kn, an, eta, A) / 100;
}

double S_factor_unitary_limit(double en, double ue)
{
	if (en < 0.0 || ue < 0.0) {
		return INFINITY;
	}
	double kd = relative_momentum(en, MU_D); // 1/fm
	double eta = sommerfeld(en, MU_D); // dimensionless
	return exp(M_PI*eta*ue/en) * en * exp(2*M_PI*eta) *\
			M_PI/(kd*kd) * (2.0*J+1.0)/((2.0*J1+1.0)*(2.0*J2+1.0)) * 1.0 / 100;
}

ext: c

max_stars_repo_name: odell/r-matrix-emulation

max_issues_repo_name: odell/r-matrix-emulation

max_forks_repo_name: odell/r-matrix-emulation


text:
/* integration/cquad.c
 *
 * Copyright (C) 2010 Pedro Gonnet
 *
 * 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 <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>

#include "cquad_const.c"

/* Allocates a workspace for the given maximum number of intervals.
    Note that if the workspace gets filled, the intervals with the
    lowest error estimates are dropped. The maximum number of
    intervals is therefore not the maximum number of intervals
    that will be computed, but merely the size of the buffer.
    */

gsl_integration_cquad_workspace *
gsl_integration_cquad_workspace_alloc (const size_t n)
{

  gsl_integration_cquad_workspace *w;

  /* Check inputs */
  if (n < 3)
    GSL_ERROR_VAL ("workspace size n must be at least 3", GSL_EDOM, 0);

  /* Allocate first the workspace struct */
  if ((w =
       (gsl_integration_cquad_workspace *)
       malloc (sizeof (gsl_integration_cquad_workspace))) == NULL)
    GSL_ERROR_VAL ("failed to allocate space for workspace struct",
		   GSL_ENOMEM, 0);

  /* Allocate the intervals */
  if ((w->ivals =
       (gsl_integration_cquad_ival *)
       malloc (sizeof (gsl_integration_cquad_ival) * n)) == NULL)
    {
      free (w);
      GSL_ERROR_VAL ("failed to allocate space for the intervals", GSL_ENOMEM,
		     0);
    }

  /* Allocate the max-heap indices */
  if ((w->heap = (size_t *) malloc (sizeof (size_t) * n)) == NULL)
    {
      free (w->ivals);
      free (w);
      GSL_ERROR_VAL ("failed to allocate space for the heap", GSL_ENOMEM, 0);
    }

  /* Remember the size of the workspace */
  w->size = n;

  /* Return the result */
  return w;

}

/* Liberates the workspace memory.
    */

void
gsl_integration_cquad_workspace_free (gsl_integration_cquad_workspace * w)
{

  /* Nothing to be done? */
  if (w == NULL)
    return;

  /* Free the intervals first */
  if (w->ivals != NULL)
    free (w->ivals);

  /* Free the heap */
  if (w->heap != NULL)
    free (w->heap);

  /* Free the structure */
  free (w);

}

/* Compute the product of the fx with one of the inverse
    Vandermonde-like matrices. */

static void
Vinvfx (const double *fx, double *c, const int d)
{

  int i, j;

  switch (d)
    {
    case 0:
      for (i = 0; i <= 4; i++)
	{
	  c[i] = 0.0;
	  for (j = 0; j <= 4; j++)
	    c[i] += V1inv[i * 5 + j] * fx[j * 8];
	}
      break;
    case 1:
      for (i = 0; i <= 8; i++)
	{
	  c[i] = 0.0;
	  for (j = 0; j <= 8; j++)
	    c[i] += V2inv[i * 9 + j] * fx[j * 4];
	}
      break;
    case 2:
      for (i = 0; i <= 16; i++)
	{
	  c[i] = 0.0;
	  for (j = 0; j <= 16; j++)
	    c[i] += V3inv[i * 17 + j] * fx[j * 2];
	}
      break;
    case 3:
      for (i = 0; i <= 32; i++)
	{
	  c[i] = 0.0;
	  for (j = 0; j <= 32; j++)
	    c[i] += V4inv[i * 33 + j] * fx[j];
	}
      break;
    }

}

/* Downdate the interpolation given by the n coefficients c
    by removing the nodes with indices in nans. */

static void
downdate (double *c, int n, int d, int *nans, int nnans)
{

  static const int bidx[4] = { 0, 6, 16, 34 };
  double b_new[34], alpha;
  int i, j;

  for (i = 0; i <= n + 1; i++)
    b_new[i] = bee[bidx[d] + i];
  for (i = 0; i < nnans; i++)
    {
      b_new[n + 1] = b_new[n + 1] / Lalpha[n];
      b_new[n] = (b_new[n] + xi[nans[i]] * b_new[n + 1]) / Lalpha[n - 1];
      for (j = n - 1; j > 0; j--)
	b_new[j] =
	  (b_new[j] + xi[nans[i]] * b_new[j + 1] -
	   Lgamma[j + 1] * b_new[j + 2]) / Lalpha[j - 1];
      for (j = 0; j <= n; j++)
	b_new[j] = b_new[j + 1];
      alpha = c[n] / b_new[n];
      for (j = 0; j < n; j++)
	c[j] -= alpha * b_new[j];
      c[n] = 0;
      n--;
    }

}

/* The actual integration routine.
    */

int
gsl_integration_cquad (const gsl_function * f, double a, double b,
		       double epsabs, double epsrel,
		       gsl_integration_cquad_workspace * ws,
		       double *result, double *abserr, size_t * nevals)
{

  /* Some constants that we will need. */
  static const int n[4] = { 4, 8, 16, 32 };
  static const int skip[4] = { 8, 4, 2, 1 };
  static const int idx[4] = { 0, 5, 14, 31 };
  static const double w = M_SQRT2 / 2;
  static const int ndiv_max = 20;

  /* Actual variables (as opposed to constants above). */
  double m, h, temp;
  double igral, err, igral_final, err_final, err_excess;
  int nivals, neval = 0;
  int i, j, d, split, t;
  int nnans, nans[32];
  gsl_integration_cquad_ival *iv, *ivl, *ivr;
  double nc, ncdiff;

  /* Check the input arguments. */
  if (f == NULL)
    GSL_ERROR ("function pointer shouldn't be NULL", GSL_EINVAL);
  if (result == NULL)
    GSL_ERROR ("result pointer shouldn't be NULL", GSL_EINVAL);
  if (ws == NULL)
    GSL_ERROR ("workspace pointer shouldn't be NULL", GSL_EINVAL);

  /* Check for unreasonable accuracy demands */
  if (epsabs < 0.0 || epsrel < 0.0)
    GSL_ERROR ("tolerances may not be negative", GSL_EBADTOL);
  if (epsabs <= 0 && epsrel < GSL_DBL_EPSILON)
    GSL_ERROR ("unreasonable accuracy requirement", GSL_EBADTOL);

  /* Create the first interval. */
  iv = &(ws->ivals[0]);
  m = (a + b) / 2;
  h = (b - a) / 2;
  nnans = 0;
  for (i = 0; i <= n[3]; i++)
    {
      iv->fx[i] = GSL_FN_EVAL (f, m + xi[i] * h);
      neval++;
      if (!finite (iv->fx[i]))
	{
	  nans[nnans++] = i;
	  iv->fx[i] = 0.0;
	}
    }
  Vinvfx (iv->fx, &(iv->c[idx[0]]), 0);
  Vinvfx (iv->fx, &(iv->c[idx[3]]), 3);
  Vinvfx (iv->fx, &(iv->c[idx[2]]), 2);
  for (i = 0; i < nnans; i++)
    iv->fx[nans[i]] = GSL_NAN;
  iv->a = a;
  iv->b = b;
  iv->depth = 3;
  iv->rdepth = 1;
  iv->ndiv = 0;
  iv->igral = 2 * h * iv->c[idx[3]] * w;
  nc = 0.0;
  for (i = n[2] + 1; i <= n[3]; i++)
    {
      temp = iv->c[idx[3] + i];
      nc += temp * temp;
    }
  ncdiff = nc;
  for (i = 0; i <= n[2]; i++)
    {
      temp = iv->c[idx[2] + i] - iv->c[idx[3] + i];
      ncdiff += temp * temp;
      nc += iv->c[idx[3] + i] * iv->c[idx[3] + i];
    }
  ncdiff = sqrt (ncdiff);
  nc = sqrt (nc);
  iv->err = ncdiff * 2 * h;
  if (ncdiff / nc > 0.1 && iv->err < 2 * h * nc)
    iv->err = 2 * h * nc;

  /* Initialize the heaps. */
  for (i = 0; i < ws->size; i++)
    ws->heap[i] = i;

  /* Initialize some global values. */
  igral = iv->igral;
  err = iv->err;
  nivals = 1;
  igral_final = 0.0;
  err_final = 0.0;
  err_excess = 0.0;

  /* Main loop. */
  while (nivals > 0 && err > 0.0 &&
	 !(err <= fabs (igral) * epsrel || err <= epsabs)
	 && !(err_final > fabs (igral) * epsrel
	      && err - err_final < fabs (igral) * epsrel)
	 && !(err_final > epsabs && err - err_final < epsabs))
    {

      /* Put our finger on the interval with the largest error. */
      iv = &(ws->ivals[ws->heap[0]]);
      m = (iv->a + iv->b) / 2;
      h = (iv->b - iv->a) / 2;

/*      printf
        ("cquad: processing ival %i (of %i) with [%e,%e] int=%e, err=%e, depth=%i\n",
         ws->heap[0], nivals, iv->a, iv->b, iv->igral, iv->err, iv->depth);
*/
      /* Should we try to increase the degree? */
      if (iv->depth < 3)
	{

	  /* Keep tabs on some variables. */
	  d = ++iv->depth;

	  /* Get the new (missing) function values */
	  for (i = skip[d]; i <= 32; i += 2 * skip[d])
	    {
	      iv->fx[i] = GSL_FN_EVAL (f, m + xi[i] * h);
	      neval++;
	    }
	  nnans = 0;
	  for (i = 0; i <= 32; i += skip[d])
	    {
	      if (!finite (iv->fx[i]))
		{
		  nans[nnans++] = i;
		  iv->fx[i] = 0.0;
		}
	    }

	  /* Compute the new coefficients. */
	  Vinvfx (iv->fx, &(iv->c[idx[d]]), d);

	  /* Downdate any NaNs. */
	  if (nnans > 0)
	    {
	      downdate (&(iv->c[idx[d]]), n[d], d, nans, nnans);
	      for (i = 0; i < nnans; i++)
		iv->fx[nans[i]] = GSL_NAN;
	    }

	  /* Compute the error estimate. */
	  nc = 0.0;
	  for (i = n[d - 1] + 1; i <= n[d]; i++)
	    {
	      temp = iv->c[idx[d] + i];
	      nc += temp * temp;
	    }
	  ncdiff = nc;
	  for (i = 0; i <= n[d - 1]; i++)
	    {
	      temp = iv->c[idx[d - 1] + i] - iv->c[idx[d] + i];
	      ncdiff += temp * temp;
	      nc += iv->c[idx[d] + i] * iv->c[idx[d] + i];
	    }
	  ncdiff = sqrt (ncdiff);
	  nc = sqrt (nc);
	  iv->err = ncdiff * 2 * h;

	  /* Compute the local integral. */
	  iv->igral = 2 * h * w * iv->c[idx[d]];

	  /* Split the interval prematurely? */
	  split = (nc > 0 && ncdiff / nc > 0.1);

	}

      /* Maximum degree reached, just split. */
      else
	{
	  split = 1;
	}

      /* Should we drop this interval? */
      if ((m + h * xi[0]) >= (m + h * xi[1])
	  || (m + h * xi[31]) >= (m + h * xi[32])
	  || iv->err < fabs (iv->igral) * GSL_DBL_EPSILON * 10)
	{

/*          printf
            ("cquad: dumping ival %i (of %i) with [%e,%e] int=%e, err=%e, depth=%i\n",
             ws->heap[0], nivals, iv->a, iv->b, iv->igral, iv->err,
             iv->depth);
*/
	  /* Keep this interval's contribution */
	  err_final += iv->err;
	  igral_final += iv->igral;

	  /* Swap with the last element on the heap */
	  t = ws->heap[nivals - 1];
	  ws->heap[nivals - 1] = ws->heap[0];
	  ws->heap[0] = t;
	  nivals--;

	  /* Fix up the heap */
	  i = 0;
	  while (2 * i + 1 < nivals)
	    {

	      /* Get the kids */
	      j = 2 * i + 1;

	      /* If the j+1st entry exists and is larger than the jth,
	         use it instead. */
	      if (j + 1 < nivals
		  && ws->ivals[ws->heap[j + 1]].err >=
		  ws->ivals[ws->heap[j]].err)
		j++;

	      /* Do we need to move the ith entry up? */
	      if (ws->ivals[ws->heap[j]].err <= ws->ivals[ws->heap[i]].err)
		break;
	      else
		{
		  t = ws->heap[j];
		  ws->heap[j] = ws->heap[i];
		  ws->heap[i] = t;
		  i = j;
		}
	    }

	}

      /* Do we need to split this interval? */
      else if (split)
	{

	  /* Some values we will need often... */
	  d = iv->depth;

	  /* Generate the interval on the left */
	  ivl = &(ws->ivals[ws->heap[nivals++]]);
	  ivl->a = iv->a;
	  ivl->b = m;
	  ivl->depth = 0;
	  ivl->rdepth = iv->rdepth + 1;
	  ivl->fx[0] = iv->fx[0];
	  ivl->fx[32] = iv->fx[16];
	  for (i = skip[0]; i < 32; i += skip[0])
	    {
	      ivl->fx[i] =
		GSL_FN_EVAL (f, (ivl->a + ivl->b) / 2 + xi[i] * h / 2);
	      neval++;
	    }
	  nnans = 0;
	  for (i = 0; i <= 32; i += skip[0])
	    {
	      if (!finite (ivl->fx[i]))
		{
		  nans[nnans++] = i;
		  ivl->fx[i] = 0.0;
		}
	    }
	  Vinvfx (ivl->fx, ivl->c, 0);
	  if (nnans > 0)
	    {
	      downdate (ivl->c, n[0], 0, nans, nnans);
	      for (i = 0; i < nnans; i++)
		ivl->fx[nans[i]] = GSL_NAN;
	    }
	  for (i = 0; i <= n[d]; i++)
	    {
	      ivl->c[idx[d] + i] = 0.0;
	      for (j = i; j <= n[d]; j++)
		ivl->c[idx[d] + i] += Tleft[i * 33 + j] * iv->c[idx[d] + j];
	    }
	  ncdiff = 0.0;
	  for (i = 0; i <= n[0]; i++)
	    {
	      temp = ivl->c[i] - ivl->c[idx[d] + i];
	      ncdiff += temp * temp;
	    }
	  for (i = n[0] + 1; i <= n[d]; i++)
	    {
	      temp = ivl->c[idx[d] + i];
	      ncdiff += temp * temp;
	    }
	  ncdiff = sqrt (ncdiff);
	  ivl->err = ncdiff * h;

	  /* Check for divergence. */
	  ivl->ndiv = iv->ndiv + (fabs (iv->c[0]) > 0
				  && ivl->c[0] / iv->c[0] > 2);
	  if (ivl->ndiv > ndiv_max && 2 * ivl->ndiv > ivl->rdepth)
	    {
              /* need copysign(INFINITY, igral) */
	      *result = (igral >= 0) ? GSL_POSINF : GSL_NEGINF;
	      if (nevals != NULL)
		*nevals = neval;
	      return GSL_EDIVERGE;
	    }

	  /* Compute the local integral. */
	  ivl->igral = h * w * ivl->c[0];

	  /* Generate the interval on the right */
	  ivr = &(ws->ivals[ws->heap[nivals++]]);
	  ivr->a = m;
	  ivr->b = iv->b;
	  ivr->depth = 0;
	  ivr->rdepth = iv->rdepth + 1;
	  ivr->fx[0] = iv->fx[16];
	  ivr->fx[32] = iv->fx[32];
	  for (i = skip[0]; i < 32; i += skip[0])
	    {
	      ivr->fx[i] =
		GSL_FN_EVAL (f, (ivr->a + ivr->b) / 2 + xi[i] * h / 2);
	      neval++;
	    }
	  nnans = 0;
	  for (i = 0; i <= 32; i += skip[0])
	    {
	      if (!finite (ivr->fx[i]))
		{
		  nans[nnans++] = i;
		  ivr->fx[i] = 0.0;
		}
	    }
	  Vinvfx (ivr->fx, ivr->c, 0);
	  if (nnans > 0)
	    {
	      downdate (ivr->c, n[0], 0, nans, nnans);
	      for (i = 0; i < nnans; i++)
		ivr->fx[nans[i]] = GSL_NAN;
	    }
	  for (i = 0; i <= n[d]; i++)
	    {
	      ivr->c[idx[d] + i] = 0.0;
	      for (j = i; j <= n[d]; j++)
		ivr->c[idx[d] + i] += Tright[i * 33 + j] * iv->c[idx[d] + j];
	    }
	  ncdiff = 0.0;
	  for (i = 0; i <= n[0]; i++)
	    {
	      temp = ivr->c[i] - ivr->c[idx[d] + i];
	      ncdiff += temp * temp;
	    }
	  for (i = n[0] + 1; i <= n[d]; i++)
	    {
	      temp = ivr->c[idx[d] + i];
	      ncdiff += temp * temp;
	    }
	  ncdiff = sqrt (ncdiff);
	  ivr->err = ncdiff * h;

	  /* Check for divergence. */
	  ivr->ndiv = iv->ndiv + (fabs (iv->c[0]) > 0
				  && ivr->c[0] / iv->c[0] > 2);
	  if (ivr->ndiv > ndiv_max && 2 * ivr->ndiv > ivr->rdepth)
	    {
              /* need copysign(INFINITY, igral) */
	      *result = (igral >= 0) ? GSL_POSINF : GSL_NEGINF;
	      if (nevals != NULL)
		*nevals = neval;
	      return GSL_EDIVERGE;
	    }

	  /* Compute the local integral. */
	  ivr->igral = h * w * ivr->c[0];

	  /* Fix-up the heap: we now have one interval on top
	     that we don't need any more and two new, unsorted
	     ones at the bottom. */

	  /* Flip the last interval to the top of the heap and
	     sift down. */
	  t = ws->heap[nivals - 1];
	  ws->heap[nivals - 1] = ws->heap[0];
	  ws->heap[0] = t;
	  nivals--;

	  /* Sift this interval back down the heap. */
	  i = 0;
	  while (2 * i + 1 < nivals - 1)
	    {
	      j = 2 * i + 1;
	      if (j + 1 < nivals - 1
		  && ws->ivals[ws->heap[j + 1]].err >=
		  ws->ivals[ws->heap[j]].err)
		j++;
	      if (ws->ivals[ws->heap[j]].err <= ws->ivals[ws->heap[i]].err)
		break;
	      else
		{
		  t = ws->heap[j];
		  ws->heap[j] = ws->heap[i];
		  ws->heap[i] = t;
		  i = j;
		}
	    }

	  /* Now grab the last interval and sift it up the heap. */
	  i = nivals - 1;
	  while (i > 0)
	    {
	      j = (i - 1) / 2;
	      if (ws->ivals[ws->heap[j]].err < ws->ivals[ws->heap[i]].err)
		{
		  t = ws->heap[j];
		  ws->heap[j] = ws->heap[i];
		  ws->heap[i] = t;
		  i = j;
		}
	      else
		break;
	    }

	}

      /* Otherwise, just fix-up the heap. */
      else
	{
	  i = 0;
	  while (2 * i + 1 < nivals)
	    {
	      j = 2 * i + 1;
	      if (j + 1 < nivals
		  && ws->ivals[ws->heap[j + 1]].err >=
		  ws->ivals[ws->heap[j]].err)
		j++;
	      if (ws->ivals[ws->heap[j]].err <= ws->ivals[ws->heap[i]].err)
		break;
	      else
		{
		  t = ws->heap[j];
		  ws->heap[j] = ws->heap[i];
		  ws->heap[i] = t;
		  i = j;
		}
	    }

	}

      /* If the heap is about to overflow, remove the last two
         intervals. */
      while (nivals > ws->size - 2)
	{
	  iv = &(ws->ivals[ws->heap[nivals - 1]]);

/*          printf
            ("cquad: dumping ival %i (of %i) with [%e,%e] int=%e, err=%e, depth=%i\n",
             ws->heap[0], nivals, iv->a, iv->b, iv->igral, iv->err,
             iv->depth);
*/
	  err_final += iv->err;
	  igral_final += iv->igral;
	  nivals--;
	}

      /* Collect the value of the integral and error. */
      igral = igral_final;
      err = err_final;
      for (i = 0; i < nivals; i++)
	{
	  igral += ws->ivals[ws->heap[i]].igral;
	  err += ws->ivals[ws->heap[i]].err;
	}

    }

  /* Dump the contents of the heap. */
/*  for (i = 0; i < nivals; i++)
    {
      iv = &(ws->ivals[ws->heap[i]]);
      printf
        ("cquad: ival %i (%i) with [%e,%e], int=%e, err=%e, depth=%i, rdepth=%i\n",
         i, ws->heap[i], iv->a, iv->b, iv->igral, iv->err, iv->depth,
         iv->rdepth);
    }
*/
  /* Clean up and present the results. */
  *result = igral;
  if (abserr != NULL)
    *abserr = err;
  if (nevals != NULL)
    *nevals = neval;

  /* All is well that ends well. */
  return GSL_SUCCESS;

}

max_stars_repo_name: ruslankuzmin/julia

max_issues_repo_name: ruslankuzmin/julia

max_forks_repo_name: ruslankuzmin/julia


text:
#ifndef MATMUL_CMAT_MUL_H
#define MATMUL_CMAT_MUL_H

#include <cblas.h>
#include <pthread.h>

#include "cmatrix.h"

#define _AVX512
//#define _AVX2
//#define _NEON

#if defined(_NEON)
#include "arm_neon.h"
#elif defined(_AVX512) || defined(_AVX2)

#include <immintrin.h>

#endif

#define BLOCKSIZE 1
#define UNROLL 4

typedef enum {
    Naive,
    Blas,
    Order,
    Simd,
    OpenCv,
    Threads,
    OpenMp,
    Unroll,
    Block,
//    UnrollMp,
} kMatMulMode;

CMatrix *MatMul(ReadOnlyMat mat1, ReadOnlyMat mat2, kMatMulMode mode, int times);

CMatrix *MatMulO3(ReadOnlyMat mat1, ReadOnlyMat mat2);

void MatMulChangeOrd(ReadOnlyMat mat1, ReadOnlyMat mat2, int times);

CMatrix *MatMulBlas(ReadOnlyMat mat1, ReadOnlyMat mat2);

//CMatrix *MatMulStrassen(ReadOnlyMat mat1, ReadOnlyMat mat2, int edge);

CMatrix *MatMulCV(ReadOnlyMat mat1, ReadOnlyMat mat2);

CMatrix *MatMulMP(ReadOnlyMat mat1, ReadOnlyMat mat2);

CMatrix *MatMulSimd(ReadOnlyMat mat1, ReadOnlyMat mat2);

CMatrix *MatMulUnroll(ReadOnlyMat mat1, ReadOnlyMat mat2);

CMatrix *MatMulBlock(ReadOnlyMat mat1, ReadOnlyMat mat2);

CMatrix *MatMulThreads(ReadOnlyMat mat1, ReadOnlyMat mat2);

//CMatrix *MatMulUnrollMp(ReadOnlyMat mat1, ReadOnlyMat mat2);

#endif //MATMUL_CMAT_MUL_H

max_stars_repo_path: Project 3/src/cmat_mul.h

max_issues_repo_path: Project 3/src/cmat_mul.h

max_forks_repo_path: Project 3/src/cmat_mul.h


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <limits.h>
#include <getopt.h>
#include <iniparser.h>
#include "parmt_utils.h"
#include "iscl/array/array.h"
#include "iscl/memory/memory.h"
#include "iscl/os/os.h"
#ifdef PARMT_USE_INTEL
#include <mkl_cblas.h>
#else
#include <cblas.h>
#endif

#define PROGRAM_NAME "mergemt"

static void printUsage(void);
static int parseArguments(int argc, char *argv[],
                          char iniFile[PATH_MAX]);

struct mergemtParms_struct
{
    char **files;
    double *wts;
    char archiveFile[PATH_MAX];
    int nfiles;
};

int mergemt_readIni(const char *iniFile, struct mergemtParms_struct *parms);

int main(int argc, char *argv[])
{
    struct mergemtParms_struct parms;
    char iniFile[PATH_MAX], programNameIn[256];
    double *beta0, *beta1, *betaLoc, *dep0, *dep1, *depLoc,
           *gamma0, *gamma1, *gammaLoc, *kappa0, *kappa1, *kappaLoc,
           *M0, *M1, *M0loc, *phi0, *phi1,
           *phiLoc, *sigma0, *sigma1, *sigmaLoc, *theta0, *theta1,
           *thetaLoc, xsum;
    int i, ierr, il, im, indx, jndx, nb, nb0, nb1, ncopy, ng, ng0, ng1,
        nk, nk0, nk1, nlocs, nlocs0, nlocs1, nm, nm0, nm1, nmt, nmt0, nmt1,
        ns, ns0, ns1, nt, nt0, nt1;
    const int jm = 0;
    ierr = parseArguments(argc, argv, iniFile);
    if (ierr != 0)
    {
        if (ierr ==-2){return EXIT_SUCCESS;}
        printf("%s: Error parsing arguments\n", PROGRAM_NAME);
        return EXIT_FAILURE;
    }
    if (!os_path_isfile(iniFile))
    {
        printf("%s: Error ini file doesn't exist\n", PROGRAM_NAME);
        return EXIT_FAILURE;
    }
    ierr = mergemt_readIni(iniFile, &parms);
    if (ierr != 0)
    {
        printf("%s: Error reading ini file\n", PROGRAM_NAME);
        return EXIT_FAILURE;
    }
    xsum = array_sum64f(parms.nfiles, parms.wts, &ierr);
    if (fabs(xsum) < 1.e-14)
    {
        printf("%s: Division by zero - setting to 1\n", PROGRAM_NAME);
        xsum = 1.0;
    }
    cblas_dscal(parms.nfiles, 1.0/xsum, parms.wts, 1);
    // Load the files
    nlocs0 =-1;
    nlocs1 =-1;
    nm0 =-1;
    nm1 =-1;
    nb0 =-1;
    nb1 =-1;
    ng0 =-1;
    ng1 =-1;
    nk0 =-1;
    nk1 =-1;
    ns0 =-1;
    ns1 =-1;
    nt0 =-1;
    nt1 =-1;
    nmt0 =-1;
    nmt1 =-1;
    dep0 = NULL;
    dep1 = NULL;
    M0 = NULL;
    M1 = NULL;
    beta0 = NULL;
    beta1 = NULL;
    gamma0 = NULL;
    gamma1 = NULL;
    kappa0 = NULL;
    kappa1 = NULL;
    sigma0 = NULL;
    sigma1 = NULL;
    theta0 = NULL;
    theta1 = NULL;
    phi0 = NULL;
    phi1 = NULL;
    for (i=0; i<parms.nfiles; i++)
    {
        if (!os_path_isfile(parms.files[i])){continue;}
        printf("%s: Loading archive: %s\n", PROGRAM_NAME, parms.files[i]);
        ierr = parmt_io_readObjfnArchive64f(parms.files[i],
                                            programNameIn,
                                            &nlocs, &depLoc,
                                            &nm, &M0loc,
                                            &nb, &betaLoc,
                                            &ng, &gammaLoc,
                                            &nk, &kappaLoc,
                                            &ns, &sigmaLoc,
                                            &nt, &thetaLoc,
                                            &nmt, &phiLoc);
        if (ierr != 0)
        {
            printf("%s: Error loading %s; skipping...\n",
                   PROGRAM_NAME, parms.files[i]);
            continue;
        }
        if (strcasecmp("parmt", programNameIn) == 0)
        {
            if (nm0 ==-1)
            {
                nlocs0 = nlocs;
                nm0 = nm;
                nb0 = nb;
                ng0 = ng;
                nk0 = nk;
                ns0 = ns;
                nt0 = nt;
                nmt0 = nmt;
                phi0 = array_zeros64f(nmt, &ierr);
                dep0 = array_copy64f(nlocs, depLoc, &ierr);
                M0 = array_copy64f(nm, M0loc, &ierr);
                beta0 = array_copy64f(nb, betaLoc, &ierr);
                gamma0 = array_copy64f(ng, gammaLoc, &ierr);
                kappa0 = array_copy64f(nk, kappaLoc, &ierr);
                sigma0 = array_copy64f(ns, sigmaLoc, &ierr);
                theta0 = array_copy64f(nt, thetaLoc, &ierr);
            }
            if (nm0 != nm || nb0 != nb || ng0 != ng || nk0 != nk ||
                ns0 != ns || nt0 != nt || nmt0 != nmt || nlocs0 != nlocs)
            {
                printf("%s: Size mismatch in parmt\n", PROGRAM_NAME);
                continue;
            }
            cblas_daxpy(nmt, parms.wts[i], phiLoc, 1, phi0, 1);
        }
        else if (strcasecmp("polarmt", programNameIn) == 0)
        {
            if (nm1 ==-1)
            {
                nlocs1 = nlocs;
                nm1 = nm;
                nb1 = nb;
                ng1 = ng;
                nk1 = nk;
                ns1 = ns;
                nt1 = nt;
                nmt1 = nmt;
                phi1 = array_zeros64f(nmt, &ierr);
                dep1 = array_copy64f(nlocs, depLoc, &ierr);
                M1 = array_copy64f(nm, M0loc, &ierr);
                beta1 = array_copy64f(nb, betaLoc, &ierr);
                gamma1 = array_copy64f(ng, gammaLoc, &ierr);
                kappa1 = array_copy64f(nk, kappaLoc, &ierr);
                sigma1 = array_copy64f(ns, sigmaLoc, &ierr);
                theta1 = array_copy64f(nt, thetaLoc, &ierr);
            }
printf("%f\n",parms.wts[i]);
            cblas_daxpy(nmt, parms.wts[i], phiLoc, 1, phi1, 1);
        }
        else
        {
            printf("%s: Unkown origin: %s\n", PROGRAM_NAME, programNameIn);
            goto NEXT;
        }

NEXT:;
        memory_free64f(&M0loc);
        memory_free64f(&betaLoc);
        memory_free64f(&gammaLoc);
        memory_free64f(&kappaLoc);
        memory_free64f(&sigmaLoc);
        memory_free64f(&thetaLoc);
        memory_free64f(&phiLoc);
    }
    if (phi0 != NULL && phi1 != NULL)
    {
        if (nb0 != nb1 || ng0 != ng1 || nk0 != nk1 ||
            ns0 != ns1 || nt0 != nt1 || nlocs0 != nlocs)
        {
            printf("%s: Inconsistent grid search sizes\n", PROGRAM_NAME);
            return -1;
        }
        printf("%f\n", array_max64f(nmt0, phi0, &ierr));
        printf("%f\n", array_max64f(nmt1, phi1, &ierr));
        printf("%s: Stacking polarities into waveforms...\n", PROGRAM_NAME);
/*
        for (i=0; i<nm0; i++)
        {
            //printf("%d %d\n", nmt, nmt0);
            cblas_daxpy(nmt1, 1.0, phi1, 1, &phi0[nmt1*i], 1);
        }
*/
        for (il=0; il<nlocs0; il++)
        {
            for (im=0; im<nm0; im++)
            {
                indx = il*nm0*nb0*ng0*nk0*ns0*nt0
                     + im*nb0*ng0*nk0*ns0*nt0;
                jndx = il*nm1*nb1*ng1*nk1*ns1*nt1
                     + jm*nb1*ng1*nk1*ns1*nt1;
                ncopy = nb0*ng0*nk0*ns0*nt0;
                cblas_daxpy(ncopy,  1.0, &phi1[jndx], 1, &phi0[indx], 1);
            }
        }
printf("%f\n", array_max64f(nmt0, phi0, &ierr));
        printf("%s: Writing joint archive\n", PROGRAM_NAME);
        ierr = parmt_io_createObjfnArchive64f(PROGRAM_NAME, parms.archiveFile,
                                              1, //data.nobs, TODO - read and fix
                                              nlocs0, dep0,
                                              nm0, M0,
                                              nb0, beta0,
                                              ng0, gamma0,
                                              nk0, kappa0,
                                              ns0, sigma0,
                                              nt0, theta0);
        ierr = parmt_io_writeObjectiveFunction64f(parms.archiveFile,
                                                  nmt0, phi0);
    }
    else
    {
        if (phi0 != NULL)
        {
            printf("%s: Writing joint waveform archive\n", PROGRAM_NAME);
            ierr = parmt_io_createObjfnArchive64f(PROGRAM_NAME, parms.archiveFile,
                                              1, //data.nobs, TODO - read and fix
                                              nlocs0, dep0,
                                              nm0, M0,
                                              nb0, beta0,
                                              ng0, gamma0,
                                              nk0, kappa0,
                                              ns0, sigma0,
                                              nt0, theta0);
            ierr = parmt_io_writeObjectiveFunction64f(parms.archiveFile,
                                                      nmt0, phi0);
        }
        else if (phi1 != NULL)
        {
            printf("%s: Writing joint waveform archive\n", PROGRAM_NAME);
            ierr = parmt_io_createObjfnArchive64f(PROGRAM_NAME, parms.archiveFile,
                                              1, //data.nobs, TODO - read and fix
                                              nlocs1, dep1,
                                              nm1, M1,
                                              nb1, beta1,
                                              ng1, gamma1,
                                              nk1, kappa1,
                                              ns1, sigma1,
                                              nt1, theta1);
        }
    }
    memory_free64f(&dep0);
    memory_free64f(&dep1);
    memory_free64f(&M0);
    memory_free64f(&M1);
    memory_free64f(&beta0);
    memory_free64f(&beta1);
    memory_free64f(&gamma0);
    memory_free64f(&gamma1);
    memory_free64f(&kappa0);
    memory_free64f(&kappa1);
    memory_free64f(&sigma0);
    memory_free64f(&sigma1);
    memory_free64f(&theta0);
    memory_free64f(&theta1);
    memory_free64f(&phi0);
    memory_free64f(&phi1);
    return EXIT_SUCCESS;
}
//============================================================================//
int mergemt_readIni(const char *iniFile,
                    struct mergemtParms_struct *parms)
{
    const char *fcnm = "mergemt_readIni\0";
    const char *s;
    char vname[256];
    dictionary *ini;
    int i, ierr;
    memset(parms, 0, sizeof(struct mergemtParms_struct));
    if (!os_path_isfile(iniFile))
    {
        printf("%s: Error ini file doesn't exist\n", fcnm);
        return -1;
    }
    ini = iniparser_load(iniFile);

    parms->nfiles = iniparser_getint(ini, "general:nmerge\0", 0);
    if (parms->nfiles < 1)
    {
        printf("%s: No files to merge\n", fcnm);
        return -1;
    }

    parms->files = (char **) calloc((size_t) parms->nfiles, sizeof(char *));
    parms->wts = array_set64f(parms->nfiles, 1.0, &ierr);
    for (i=0; i<parms->nfiles; i++)
    {
        parms->files[i] = (char *) calloc(PATH_MAX, sizeof(char));
        memset(vname, 0, 256*sizeof(char));
        sprintf(vname, "general:mergeFile_%d", i+1);
        s = iniparser_getstring(ini, vname, NULL);
        if (!os_path_isfile(s))
        {
            printf("%s: File %s doesn't exist\n", fcnm, s);
            continue;
        }
        strcpy(parms->files[i], s);

        memset(vname, 0, 256*sizeof(char));
        sprintf(vname, "general:wtFile_%d", i+1);
        parms->wts[i] = iniparser_getdouble(ini, vname, 1.0);
        if (parms->wts[i] < 0.0)
        {
            printf("%s: Weight can't be negative: %f\n", fcnm, parms->wts[i]);
            return -1;
        }
    }

    memset(vname, 0, 256*sizeof(char));
    strcpy(vname, "general:mergeArchive");
    s = iniparser_getstring(ini, vname, "mergemt.h5");
    strcpy(parms->archiveFile, s);

    iniparser_freedict(ini);
    return 0;
}
//============================================================================//
static int parseArguments(int argc, char *argv[],
                          char iniFile[PATH_MAX])
{
    bool linFile;
    int prod;
    linFile = false;
    memset(iniFile, 0, PATH_MAX*sizeof(char));
    while (true)
    {
        static struct option longOptions[] =
        {
            {"help", no_argument, 0, '?'},
            {"help", no_argument, 0, 'h'},
            {"ini_file", required_argument, 0, 'i'},
            {0, 0, 0, 0}
        };
        int c, optionIndex;
        c = getopt_long(argc, argv, "?hi:m:l:",
                        longOptions, &optionIndex);
        if (c ==-1){break;}
        if (c == 'i')
        {
            strcpy(iniFile, (const char *) optarg);
            linFile = true;
        }
        else if (c == 'h' || c == '?')
        {
            printUsage();
            return -2;
        }
        else
        {
            printf("%s: Unknown options: %s\n",
                   PROGRAM_NAME, argv[optionIndex]);
        }
    }
    if (!linFile)
    {
        printf("%s: Error must specify ini file\n\n", PROGRAM_NAME);
        printUsage();
        return -1;
    }
    return 0;
}
//============================================================================//
static void printUsage(void)
{
    printf("Usage:\n   mergemt -i input_file\n\n");
    printf("Required arguments:\n");
    printf("   -i input_file specifies the initialization file\n");
    printf("\n");
    printf("Optional arguments:\n");
    printf("   -h displays this message\n");
    return;
}


text:
#include <math.h>
#include <stdio.h>
#include <stdarg.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_statistics_double.h>
#include "../include/type.h"
#include "../include/util.h"
#include "../include/myfunc.h"
#include "../include/global.h"
#include "../include/cthreadpool.h"

extern double rho;
extern double cutoff;
extern clock_t dur;

#define ITER_CUTOFF 0.01

/**
 * @brief invoke this function before return from
 *        MLE_marginal_iteration and MLE_marginal_iteration_constrain.
 *
 * @param mle
 * @param sum
 * @param psi1
 * @param psi2
 * @param beta0
 * @param beta1
 * @param var1
 * @param var2
 */
void mle_result_set(mle_result* mle, double sum, gsl_vector* psi1, gsl_vector* psi2,
                    double beta0, double beta1, double var1, double var2) {
    mle->sum = sum;
    mle->params.psi1 = psi1;
    mle->params.psi2 = psi2;
    mle->params.beta0 = beta0;
    mle->params.beta1 = beta1;
    mle->params.var1 = var1;
    mle->params.var2 = var2;
    return;
}

double myfunc_multivar(const double x[], va_list argv) {
    gsl_vector *psi1 = va_arg(argv, gsl_vector*);
    gsl_vector *psi2 = va_arg(argv, gsl_vector*);
    double var1 = va_arg(argv, double), var2 = va_arg(argv, double);
    double sum1, sum2;

    // TODO  memcpy replaced.
    sum1 = cuscumsum(psi1, sum_for_multivar, 1, x[0]);
    sum2 = cuscumsum(psi2, sum_for_multivar, 1, x[1]);
    sum1 = sum1/(var1*2), sum2 = sum2/(var2*2);

    return sum1 + sum2 + 0.05 * pow(rho, 2) / (1 - pow(rho, 2)) * \
        (pow(gsl_cdf_gaussian_Pinv(x[0], 1), 2) + pow(gsl_cdf_gaussian_Pinv(x[1], 1), 2) - \
         2 * rho * gsl_cdf_gaussian_Pinv(x[0], 1) * gsl_cdf_gaussian_Pinv(x[1], 1));
}

void myfunc_multivar_der(const double x[], double res[], va_list argv) {
    gsl_vector *psi1 = va_arg(argv, gsl_vector*);
    gsl_vector *psi2 = va_arg(argv, gsl_vector*);
    double var1 = va_arg(argv, double), var2 = va_arg(argv, double);
    double sum1, sum2, tmp = pow(rho,2)/(1-pow(rho,2));
    double ppfx0 = gsl_cdf_gaussian_Pinv(x[0], 1);
    double ppfx1 = gsl_cdf_gaussian_Pinv(x[1], 1);

    sum1 = cuscumsum(psi1, sum_for_multivar_der, 1, x[0]);
    sum2 = cuscumsum(psi2, sum_for_multivar_der, 1, x[1]);
    sum1 = sum1/(var1*2), sum2 = sum2/(var2*2);

    // TODO are there any better way to convert cdf to pdf?
    res[0] = sum1 + tmp * 0.1*(ppfx0-rho*ppfx1)/gsl_ran_gaussian_pdf(ppfx0, 1);
    res[1] = sum2 + tmp * 0.1*(ppfx1-rho*ppfx0)/gsl_ran_gaussian_pdf(ppfx1, 1);

    return;
}

double myfunc_1_2(const double x[], va_list argv) {
    int flag = va_arg(argv, int);
    gsl_vector *psi1 = va_arg(argv, gsl_vector*);
    gsl_vector *psi2 = va_arg(argv, gsl_vector*);
    double var1 = va_arg(argv, double), var2 = va_arg(argv, double);
    double sum1, sum2;

    if (flag == 1) {
        sum1 = cuscumsum(psi1, sum_for_multivar, 1, x[0]+cutoff);
        sum2 = cuscumsum(psi2, sum_for_multivar, 1, x[0]);
    } else {
        sum1 = cuscumsum(psi1, sum_for_multivar, 1, x[0]);
        sum2 = cuscumsum(psi2, sum_for_multivar, 1, x[0]+cutoff);
    }
    sum1 = sum1/(var1*2), sum2 = sum2/(var2*2);

    return sum1 + sum2 + 0.05 * pow(rho, 2) / (1 - pow(rho, 2)) * \
        (pow(gsl_cdf_gaussian_Pinv(x[0]+cutoff, 1), 2) + pow(gsl_cdf_gaussian_Pinv(x[0], 1), 2) - \
         2 * rho * gsl_cdf_gaussian_Pinv(x[0]+cutoff, 1) * gsl_cdf_gaussian_Pinv(x[0], 1));
}

void myfunc_der_1_2(const double x[], double res[], va_list argv) {
    int flag = va_arg(argv, int);
    gsl_vector *psi1 = va_arg(argv, gsl_vector*);
    gsl_vector *psi2 = va_arg(argv, gsl_vector*);
    double var1 = va_arg(argv, double), var2 = va_arg(argv, double);
    double sum1, sum2, tmp = pow(rho,2)/(1-pow(rho,2));
    double ppfx0 = gsl_cdf_gaussian_Pinv(x[0], 1);
    double ppfxc = gsl_cdf_gaussian_Pinv(x[0]+cutoff, 1);

    if (flag == 1) {
        sum1 = cuscumsum(psi1, sum_for_multivar_der, 1, x[0]+cutoff);
        sum2 = cuscumsum(psi2, sum_for_multivar_der, 1, x[0]);
    } else {
        sum1 = cuscumsum(psi1, sum_for_multivar_der, 1, x[0]);
        sum2 = cuscumsum(psi2, sum_for_multivar_der, 1, x[0]+cutoff);
    }
    sum1 = sum1/(var1*2), sum2 = sum2/(var2*2);

    // TODO are there any better way to convert cdf to pdf?
    res[0] = sum1 + tmp * 0.1*(ppfxc-rho*ppfx0)/gsl_ran_gaussian_pdf(ppfxc, 1);
    res[1] = sum2 + tmp * 0.1*(ppfx0-rho*ppfxc)/gsl_ran_gaussian_pdf(ppfx0, 1);
    res[0] += res[1];

    return;
}

double myfunc_marginal(const double x[], va_list argv) {
    gsl_vector* I = va_arg(argv, gsl_vector*), *S = va_arg(argv, gsl_vector*);
    gsl_vector* psi = va_arg(argv, gsl_vector*);
    double var = va_arg(argv, double), sum = 0;
    int inclu_len = va_arg(argv, int), skip_len = va_arg(argv, int), idx = 0;

    sum = cuscumsum(psi, sum_for_marginal, 7, &idx, x[0], I, S, var, inclu_len, skip_len);

    return sum;
}

void myfunc_marginal_der(const double x[], double res[], va_list argv) {
    gsl_vector* I = va_arg(argv, gsl_vector*), *S = va_arg(argv, gsl_vector*);
    gsl_vector* psi = va_arg(argv, gsl_vector*);
    double var = va_arg(argv, double), sum;
    int inclu_len = va_arg(argv, int), skip_len = va_arg(argv, int), idx = 0;

    sum = cuscumsum(psi, sum_for_marginal_der, 7, &idx, x[0], I, S, var, inclu_len, skip_len);
    res[0] = sum;

    return;
}

double myfunc_marginal_1_2(const double x[], va_list argv) {
    int flag = va_arg(argv, int);
    double beta1, beta2;
    if (flag == 1) {
        beta1 = x[0] + cutoff, beta2 = x[0];
    } else {
        beta1 = x[0], beta2 = x[0] + cutoff;
    }
    // TODO parameter ordering
    return myfunc_marginal(&beta1, argv) + myfunc_marginal(&beta2, argv);
}

void myfunc_marginal_1_2_der(const double x[], double res[], va_list argv) {
    int flag = va_arg(argv, int);
    double beta1, beta2;
    double tmp1, tmp2;
    if (flag == 1) {
        beta1 = x[0] + cutoff, beta2 = x[0];
    } else {
        beta1 = x[0], beta2 = x[0] + cutoff;
    }
    // TODO parameter ordering
    myfunc_marginal_der(&beta1, &tmp1, argv);
    myfunc_marginal_der(&beta2, &tmp2, argv);
    res[0] = tmp1 + tmp2;

    return;
}

double myfunc_individual(const double x[], va_list argv) {
    double I = va_arg(argv, double), S = va_arg(argv, double);
    double beta = va_arg(argv, double), var = va_arg(argv, double);
    int inclu_len = va_arg(argv, int), skip_len = va_arg(argv, int);
    double new_psi = inclu_len * x[0]/(inclu_len * x[0] + skip_len * (1 - x[0]));

    // TODO This change the result.
    return pow((logit(x[0]) - logit(beta)), 2)/(2 * var) -
           (I * log(new_psi) + S * log(1 - new_psi) - log(x[0]) - log(1-x[0]));
    // return -(I * log(new_psi) + S * log(1 - new_psi) -
    //         pow((logit(x[0]) - logit(beta)), 2)/(2 * var) + log(x[0] - x[0]*x[0]));
    // return -(I * log(new_psi) + S * log(1 - new_psi) -
    //         pow((logit(x[0]) - logit(beta)), 2)/(2 * var) + log(x[0]) + log(1-x[0]));
}

void myfunc_individual_der(const double x[], double res[], va_list argv) {
    double I = va_arg(argv, double), S = va_arg(argv, double);
    double beta = va_arg(argv, double), var = va_arg(argv, double);
    int inclu_len = va_arg(argv, int), skip_len = va_arg(argv, int);
    double new_psi = inclu_len * x[0]/(inclu_len * x[0] + skip_len * (1 - x[0]));
    double new_psi_der = inclu_len * skip_len/pow(inclu_len * x[0] + skip_len * (1 - x[0]), 2);

    res[0] = 1/x[0] + S/(1 - new_psi) * new_psi_der + \
             (logit(x[0]) - logit(beta))/(var * x[0] * (1 - x[0])) - (1/(1-x[0]) + I/new_psi * new_psi_der);

    return;
}

#define individual_for_loop() \
    for (idx = 0; idx < psi1->size; ++idx) { \
        x[0] = gsl_vector_get(psi1, idx); \
        cur_sum += l_bfgs_b_wrapper(n, m, x, l, u, nbd, \
                                    myfunc_individual, myfunc_individual_der, \
                                    factr, pgtol, iprint, 15000, 15000, \
                                    6, gsl_vector_get(i1, idx), gsl_vector_get(s1, idx), \
                                    beta0, var1, inclu_len, skip_len); \
        gsl_vector_set(psi1, idx, x[0]); \
    } \
    for (idx = 0; idx < psi2->size; ++idx) { \
        x[0] = gsl_vector_get(psi2, idx); \
        cur_sum += l_bfgs_b_wrapper(n, m, x, l, u, nbd, \
                                    myfunc_individual, myfunc_individual_der, \
                                    factr, pgtol, iprint, 15000, 15000, \
                                    6, gsl_vector_get(i2, idx), gsl_vector_get(s2, idx), \
                                    beta1, var2, inclu_len, skip_len); \
        gsl_vector_set(psi2, idx, x[0]); \
    } \

int MLE_marginal_iteration(gsl_vector* i1, gsl_vector* i2,
                           gsl_vector* s1, gsl_vector* s2,
                           const int inclu_len, const int skip_len,
                           mle_result* mle) {
    gsl_vector* psi1 = gsl_vector_alloc(i1->size); // TODO deallocation.
    gsl_vector* psi2 = gsl_vector_alloc(i2->size);
    vec2psi(psi1, i1, s1, inclu_len, skip_len);
    vec2psi(psi2, i2, s2, inclu_len, skip_len);
    int iter_max = 100, count = 0, i = 0;
    size_t idx = 0;
    double prev_sum = 0, cur_sum = 0;
    double var1 = 10 * gsl_stats_variance(psi1->data, 1, psi1->size) * (psi1->size-1)/(psi1->size);
    double var2 = 10 * gsl_stats_variance(psi2->data, 1, psi2->size) * (psi2->size-1)/(psi2->size);
    double beta0, beta1, iter_cutoff = 1;
    integer n, m = 10, nbd[nmax];
    doublereal x[nmax], l[nmax], u[nmax], factr=1.0e7, pgtol=1.0e-5, iprint=-1;

    if (var1 <= 0.01 || psi1->size == 1) {
        var1 = 0.01;
    }
    if (var2 <= 0.01 || psi2->size == 1) {
        var2 = 0.01;
    }

    // TODO & or && ?
    // According to original python code 'while((iter_cutoff>0.01)&(count<=iter_maxrun)):',
    // it's a bit arithmetic '&'. However,  it should be a logical 'and' in such senario.
    while((iter_cutoff > ITER_CUTOFF) && (count <= iter_max)) {
        ++count, n = 2;
        x[0] = gsl_stats_mean(psi1->data, 1, psi1->size);
        x[1] = gsl_stats_mean(psi2->data, 1, psi2->size);
        for (i = 0; i < n; ++i) {
            nbd[i] = 2;
            l[i] = 0.01;
            u[i] = 0.99;
        }
        l_bfgs_b_wrapper(n, m, x, l, u, nbd,
                         myfunc_multivar, myfunc_multivar_der, factr, pgtol,
                         iprint, 15000, 15000, 4, psi1, psi2, var1, var2);
        beta0 = x[0], beta1 = x[1];
        n = 1, cur_sum = 0;

        individual_for_loop();

        if (count > 1) {
            iter_cutoff = fabs(prev_sum - cur_sum);
        }
        prev_sum = cur_sum;
    }
    if (count > iter_max) {
        mle_result_set(mle, cur_sum, psi1, psi2, 0, 0, var1, var2);
        return 0;
    }

    iter_cutoff = 1, iter_max = 100, count = 0, prev_sum = 0;
    while((iter_cutoff > ITER_CUTOFF) && (count <= iter_max)) {
        ++count, n = 1;
        x[0] = beta0;
        l_bfgs_b_wrapper(n, m, x, l, u, nbd,
                         myfunc_marginal, myfunc_marginal_der,
                         factr, pgtol, iprint, 15000, 15000,
                         6, i1, s1, psi1, var1, inclu_len, skip_len);
        beta0 = x[0];
        x[0] = beta1;
        l_bfgs_b_wrapper(n, m, x, l, u, nbd,
                         myfunc_marginal, myfunc_marginal_der,
                         factr, pgtol, iprint, 15000, 15000,
                         6, i2, s2, psi2, var2, inclu_len, skip_len);
        beta1 = x[0];

        cur_sum = 0;

        individual_for_loop();

        if (count > 1) {
            iter_cutoff = fabs(prev_sum - cur_sum);
        }
        prev_sum = cur_sum;
    }

    gsl_vector_free(psi1);
    gsl_vector_free(psi2);

    if (count > iter_max) {
        mle_result_set(mle, cur_sum, psi1, psi2, 0, 0, var1, var2);
        return 0;
    }

    mle_result_set(mle, cur_sum, psi1, psi2, beta0, beta1, var1, var2);

    return 0;
}

int MLE_marginal_iteration_constrain(gsl_vector* i1, gsl_vector* i2,
                                     gsl_vector* s1, gsl_vector* s2,
                                     const int inclu_len, const int skip_len,
                                     mle_result* mle) {
    gsl_vector* psi1 = gsl_vector_alloc(i1->size); // TODO deallocation.
    gsl_vector* psi2 = gsl_vector_alloc(i2->size);
    vec2psi(psi1, i1, s1, inclu_len, skip_len);
    vec2psi(psi2, i2, s2, inclu_len, skip_len);
    int iter_max = 100, count = 0, i = 0;
    size_t idx = 0;
    double prev_sum = 0, cur_sum = 0, iter_cutoff = 1;
    double var1 = 10 * gsl_stats_variance(psi1->data, 1, psi1->size) * (psi1->size-1)/(psi1->size);
    double var2 = 10 * gsl_stats_variance(psi2->data, 1, psi2->size) * (psi2->size-1)/(psi2->size);
    double beta0, beta1;
    integer n = 1, m = 10, nbd[nmax];
    doublereal x[nmax], l[nmax], u[nmax], factr=1.0e7, pgtol=1.0e-5, iprint=-1;

    if (var1 <= 0.01 || psi1->size == 1) {
        var1 = 0.01;
    }
    if (var2 <= 0.01 || psi2->size == 1) {
        var2 = 0.01;
    }

    // TODO & or && ?
    // According to original python code 'while((iter_cutoff>0.01)&(count<=iter_maxrun)):',
    // it's a bit arithmetic '&'. However,  it should be a logical 'and' in such senario.
    while((iter_cutoff > ITER_CUTOFF) && (count <= iter_max)) {
        ++count;
        x[0] = gsl_stats_mean(psi1->data, 1, psi1->size);
        x[1] = gsl_stats_mean(psi2->data, 1, psi2->size);
        for (i = 0; i < n; ++i) {
            nbd[i] = 2;
            l[i] = 0.001;
            u[i] = 0.999-cutoff;
        }
        if (x[0] > x[1]) {
            l_bfgs_b_wrapper(n, m, x+1, l, u, nbd,
                             myfunc_1_2, myfunc_der_1_2, factr, pgtol,
                             iprint, 15000, 15000, 5, 1, psi1, psi2, var1, var2);
            beta1 = fmax(fmin(x[1], 1-cutoff), 0);
            beta0 = beta1 + cutoff;
        } else {
            l_bfgs_b_wrapper(n, m, x, l, u, nbd,
                             myfunc_1_2, myfunc_der_1_2, factr, pgtol,
                             iprint, 15000, 15000, 5, 2, psi1, psi2, var1, var2);
            beta0 = fmax(fmin(x[0], 1-cutoff), 0);
            beta1 = beta0 + cutoff;
        }
        cur_sum = 0;
        for (i = 0; i < n; ++i) {
            nbd[i] = 2;
            l[i] = 0.01;
            u[i] = 0.99;
        }

        individual_for_loop();

        if (count > 1) {
            iter_cutoff = fabs(prev_sum - cur_sum);
        }
        prev_sum = cur_sum;
    }
    if (count > iter_max) {
        mle_result_set(mle, cur_sum, psi1, psi2, 0, 0, var1, var2);
        return 0;
    }

    iter_cutoff = 1, iter_max = 100, count = 0, prev_sum = 0;
    while((iter_cutoff > ITER_CUTOFF) && (count <= iter_max)) {
        ++count;
        for (i = 0; i < n; ++i) {
            nbd[i] = 2;
            l[i] = 0.001;
            u[i] = 0.999-cutoff;
        }
        x[0] = beta0;
        x[1] = beta1;
        if (gsl_stats_mean(psi1->data, 1, psi1->size) > gsl_stats_mean(psi2->data, 1, psi2->size)) {
            l_bfgs_b_wrapper(n, m, x+1, l, u, nbd,
                             myfunc_marginal_1_2, myfunc_marginal_1_2_der, factr, pgtol,
                             iprint, 15000, 15000, 13, 1, i1, s1, psi1, var1, inclu_len, skip_len,
                             i2, s2, psi2, var2, inclu_len, skip_len);
            beta1 = fmax(fmin(x[1], 1-cutoff), 0);
            beta0 = beta1 + cutoff;
        } else {
            l_bfgs_b_wrapper(n, m, x, l, u, nbd,
                             myfunc_marginal_1_2, myfunc_marginal_1_2_der, factr, pgtol,
                             iprint, 15000, 15000, 13, 2, i1, s1, psi1, var1, inclu_len, skip_len,
                             i2, s2, psi2, var2, inclu_len, skip_len);
            beta0 = fmax(fmin(x[0], 1-cutoff), 0);
            beta1 = beta0 + cutoff;
        }
        cur_sum = 0;
        for (i = 0; i < n; ++i) {
            nbd[i] = 2;
            l[i] = 0.01;
            u[i] = 0.99;
        }

        individual_for_loop();

        if (count > 1) {
            iter_cutoff = fabs(prev_sum - cur_sum);
        }
        prev_sum = cur_sum;
        // printf("%.12f, %.12f, %.12f, %.12f\n", beta0, beta1, iter_cutoff, prev_sum);
    }

    gsl_vector_free(psi1);
    gsl_vector_free(psi2);

    if (count > iter_max) {
        mle_result_set(mle, cur_sum, psi1, psi2, 0, 0, var1, var2);
        return 0;
    }

    mle_result_set(mle, cur_sum, psi1, psi2, beta0, beta1, var1, var2);

    return 0;
}

void* thread_wrapper_for_LT(void* arg) {
    double *ret = (double*)malloc(sizeof(double));
    odiff *data = (odiff*)arg;
    *ret = likelihood_test(data->inc1, data->inc2,
                           data->skp1, data->skp2,
                           data->inclu_len, data->skip_len,
                           data->flag, data->id);
    return (void*)ret;
}

void* batch_wrapper_for_LT(void* arg) {
    batch_datum *args = (batch_datum*)arg;
    int batch_size = args->batch_size, i;
    double *ret = (double*)malloc(sizeof(double)*batch_size);
    odiff *datum = (odiff*)*args->datum;
    odiff data;
    for (i = 0; i < batch_size; ++i) {
        data = datum[i];
        ret[i] = likelihood_test(data.inc1, data.inc2,
                                 data.skp1, data.skp2,
                                 data.inclu_len, data.skip_len,
                                 data.flag, data.id);
    }
    return (void*)ret;
}

double likelihood_test(gsl_vector *i1, gsl_vector *i2, gsl_vector *s1, gsl_vector *s2,
                       int inclu_len, int skip_len, int flag, char* id) {
    printf("Testing %s\n", id);
    if (!flag) {
        printf("1 return from: %s\n", id);
        return 1;
    } else {
        mle_result mle;
        MLE_marginal_iteration(i1, i2, s1, s2, inclu_len, skip_len, &mle);
        if (fabs(mle.params.beta0 - mle.params.beta1) <= cutoff) {
            printf("2 return from: %s\n", id);
			return 1;
        } else {
            mle_result mle_constrain;
            MLE_marginal_iteration_constrain(i1, i2, s1, s2, inclu_len, skip_len, &mle_constrain);
            printf("3 return from: %s\n", id);
            return 1 - gsl_cdf_chisq_P(2 * (fabs(mle_constrain.sum - mle.sum)), 1);
        }
    }
}

int vec2psi(gsl_vector* psi, gsl_vector *inc, gsl_vector *skp,
            int inclu_len, int skip_len) {
    size_t idx;
    for(idx = 0; idx < inc->size; ++idx) {
        if (gsl_vector_get(inc, idx) + gsl_vector_get(skp, idx) == 0) {
            gsl_vector_set(psi, idx, 0.5);
        } else {
            gsl_vector_set(psi, idx, gsl_vector_get(inc, idx)/
                    (gsl_vector_get(inc, idx) + inclu_len * gsl_vector_get(skp, idx)/skip_len));
        }
    }
    return 0;
}

max_stars_repo_name: chunjie-sam-liu/rmats-turbo

max_issues_repo_name: chunjie-sam-liu/rmats-turbo

max_forks_repo_name: chunjie-sam-liu/rmats-turbo