text:
/* specfunc/gsl_sf_bessel.h
 *
 * 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 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.
 */

/* Author:  G. Jungman */

#ifndef __GSL_SF_BESSEL_H__
#define __GSL_SF_BESSEL_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <stdlib.h>
#include <gsl/gsl_mode.h>
#include <gsl/gsl_precision.h>
#include <gsl/gsl_sf_result.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

/* Regular Bessel Function J_0(x)
 *
 * exceptions: none
 */
GSL_FUN int gsl_sf_bessel_J0_e(const double x,  gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_J0(const double x);

/* Regular Bessel Function J_1(x)
 *
 * exceptions: GSL_EUNDRFLW
 */
GSL_FUN int gsl_sf_bessel_J1_e(const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_J1(const double x);

/* Regular Bessel Function J_n(x)
 *
 * exceptions: GSL_EUNDRFLW
 */
GSL_FUN int gsl_sf_bessel_Jn_e(int n, double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_Jn(const int n, const double x);

/* Regular Bessel Function J_n(x),  nmin <= n <= nmax
 *
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
GSL_FUN int gsl_sf_bessel_Jn_array(int nmin, int nmax, double x, double * result_array);

/* Irregular Bessel function Y_0(x)
 *
 * x > 0.0
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
GSL_FUN int gsl_sf_bessel_Y0_e(const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_Y0(const double x);

/* Irregular Bessel function Y_1(x)
 *
 * x > 0.0
 * exceptions: GSL_EDOM, GSL_EOVRFLW, GSL_EUNDRFLW
 */
GSL_FUN int gsl_sf_bessel_Y1_e(const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_Y1(const double x);

/* Irregular Bessel function Y_n(x)
 *
 * x > 0.0
 * exceptions: GSL_EDOM, GSL_EOVRFLW, GSL_EUNDRFLW
 */
GSL_FUN int gsl_sf_bessel_Yn_e(int n,const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_Yn(const int n,const double x);

/* Irregular Bessel function Y_n(x), nmin <= n <= nmax
 *
 * x > 0.0
 * exceptions: GSL_EDOM, GSL_EOVRFLW, GSL_EUNDRFLW
 */
GSL_FUN int gsl_sf_bessel_Yn_array(const int nmin, const int nmax, const double x, double * result_array);

/* Regular modified Bessel function I_0(x)
 *
 * exceptions: GSL_EOVRFLW
 */
GSL_FUN int gsl_sf_bessel_I0_e(const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_I0(const double x);

/* Regular modified Bessel function I_1(x)
 *
 * exceptions: GSL_EOVRFLW, GSL_EUNDRFLW
 */
GSL_FUN int gsl_sf_bessel_I1_e(const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_I1(const double x);

/* Regular modified Bessel function I_n(x)
 *
 * exceptions: GSL_EOVRFLW, GSL_EUNDRFLW
 */
GSL_FUN int gsl_sf_bessel_In_e(const int n, const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_In(const int n, const double x);

/* Regular modified Bessel function  I_n(x) for n=nmin,...,nmax
 *
 * nmin >=0, nmax >= nmin
 * exceptions: GSL_EDOM, GSL_EOVRFLW, GSL_EUNDRFLW
 */
GSL_FUN int gsl_sf_bessel_In_array(const int nmin, const int nmax, const double x, double * result_array);

/* Scaled regular modified Bessel function
 *  exp(-|x|) I_0(x)
 *
 * exceptions: none
 */
GSL_FUN int gsl_sf_bessel_I0_scaled_e(const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_I0_scaled(const double x);

/* Scaled regular modified Bessel function
 *  exp(-|x|) I_1(x)
 *
 * exceptions: GSL_EUNDRFLW
 */
GSL_FUN int gsl_sf_bessel_I1_scaled_e(const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_I1_scaled(const double x);

/* Scaled regular modified Bessel function
 *  exp(-|x|) I_n(x)
 *
 * exceptions: GSL_EUNDRFLW
 */
GSL_FUN int gsl_sf_bessel_In_scaled_e(int n, const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_In_scaled(const int n, const double x);

/* Scaled regular modified Bessel function
 *  exp(-|x|) I_n(x)  for n=nmin,...,nmax
 *
 * nmin >=0, nmax >= nmin
 * exceptions: GSL_EUNDRFLW
 */
GSL_FUN int gsl_sf_bessel_In_scaled_array(const int nmin, const int nmax, const double x, double * result_array);

/* Irregular modified Bessel function K_0(x)
 *
 * x > 0.0
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
GSL_FUN int gsl_sf_bessel_K0_e(const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_K0(const double x);

/* Irregular modified Bessel function K_1(x)
 *
 * x > 0.0
 * exceptions: GSL_EDOM, GSL_EOVRFLW, GSL_EUNDRFLW
 */
GSL_FUN int gsl_sf_bessel_K1_e(const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_K1(const double x);

/* Irregular modified Bessel function K_n(x)
 *
 * x > 0.0
 * exceptions: GSL_EDOM, GSL_EOVRFLW, GSL_EUNDRFLW
 */
GSL_FUN int gsl_sf_bessel_Kn_e(const int n, const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_Kn(const int n, const double x);

/* Irregular modified Bessel function  K_n(x)  for n=nmin,...,nmax
 *
 * x > 0.0, nmin >=0, nmax >= nmin
 * exceptions: GSL_EDOM, GSL_EOVRFLW, GSL_EUNDRFLW
 */
GSL_FUN int gsl_sf_bessel_Kn_array(const int nmin, const int nmax, const double x, double * result_array);

/* Scaled irregular modified Bessel function
 *  exp(x) K_0(x)
 *
 * x > 0.0
 * exceptions: GSL_EDOM
 */
GSL_FUN int gsl_sf_bessel_K0_scaled_e(const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_K0_scaled(const double x);

/* Scaled irregular modified Bessel function
 *  exp(x) K_1(x)
 *
 * x > 0.0
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
GSL_FUN int gsl_sf_bessel_K1_scaled_e(const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_K1_scaled(const double x);

/* Scaled irregular modified Bessel function
 *  exp(x) K_n(x)
 *
 * x > 0.0
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
GSL_FUN int gsl_sf_bessel_Kn_scaled_e(int n, const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_Kn_scaled(const int n, const double x);

/* Scaled irregular modified Bessel function  exp(x) K_n(x)  for n=nmin,...,nmax
 *
 * x > 0.0, nmin >=0, nmax >= nmin
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
GSL_FUN int gsl_sf_bessel_Kn_scaled_array(const int nmin, const int nmax, const double x, double * result_array);

/* Regular spherical Bessel function j_0(x) = sin(x)/x
 *
 * exceptions: none
 */
GSL_FUN int gsl_sf_bessel_j0_e(const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_j0(const double x);

/* Regular spherical Bessel function j_1(x) = (sin(x)/x - cos(x))/x
 *
 * exceptions: GSL_EUNDRFLW
 */
GSL_FUN int gsl_sf_bessel_j1_e(const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_j1(const double x);

/* Regular spherical Bessel function j_2(x) = ((3/x^2 - 1)sin(x) - 3cos(x)/x)/x
 *
 * exceptions: GSL_EUNDRFLW
 */
GSL_FUN int gsl_sf_bessel_j2_e(const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_j2(const double x);

/* Regular spherical Bessel function j_l(x)
 *
 * l >= 0, x >= 0.0
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
GSL_FUN int gsl_sf_bessel_jl_e(const int l, const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_jl(const int l, const double x);

/* Regular spherical Bessel function j_l(x) for l=0,1,...,lmax
 *
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
GSL_FUN int gsl_sf_bessel_jl_array(const int lmax, const double x, double * result_array);

/* Regular spherical Bessel function j_l(x) for l=0,1,...,lmax
 * Uses Steed's method.
 *
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
GSL_FUN int gsl_sf_bessel_jl_steed_array(const int lmax, const double x, double * jl_x_array);

/* Irregular spherical Bessel function y_0(x)
 *
 * exceptions: none
 */
GSL_FUN int gsl_sf_bessel_y0_e(const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_y0(const double x);

/* Irregular spherical Bessel function y_1(x)
 *
 * exceptions: GSL_EUNDRFLW
 */
GSL_FUN int gsl_sf_bessel_y1_e(const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_y1(const double x);

/* Irregular spherical Bessel function y_2(x)
 *
 * exceptions: GSL_EUNDRFLW
 */
GSL_FUN int gsl_sf_bessel_y2_e(const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_y2(const double x);

/* Irregular spherical Bessel function y_l(x)
 *
 * exceptions: GSL_EUNDRFLW
 */
GSL_FUN int gsl_sf_bessel_yl_e(int l, const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_yl(const int l, const double x);

/* Irregular spherical Bessel function y_l(x) for l=0,1,...,lmax
 *
 * exceptions: GSL_EUNDRFLW
 */
GSL_FUN int gsl_sf_bessel_yl_array(const int lmax, const double x, double * result_array);

/* Regular scaled modified spherical Bessel function
 *
 * Exp[-|x|] i_0(x)
 *
 * exceptions: none
 */
GSL_FUN int gsl_sf_bessel_i0_scaled_e(const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_i0_scaled(const double x);

/* Regular scaled modified spherical Bessel function
 *
 * Exp[-|x|] i_1(x)
 *
 * exceptions: GSL_EUNDRFLW
 */
GSL_FUN int gsl_sf_bessel_i1_scaled_e(const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_i1_scaled(const double x);

/* Regular scaled modified spherical Bessel function
 *
 * Exp[-|x|] i_2(x)
 *
 * exceptions: GSL_EUNDRFLW
 */
GSL_FUN int gsl_sf_bessel_i2_scaled_e(const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_i2_scaled(const double x);

/* Regular scaled modified spherical Bessel functions
 *
 * Exp[-|x|] i_l(x)
 *
 * i_l(x) = Sqrt[Pi/(2x)] BesselI[l+1/2,x]
 *
 * l >= 0
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
GSL_FUN int gsl_sf_bessel_il_scaled_e(const int l, double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_il_scaled(const int l, const double x);

/* Regular scaled modified spherical Bessel functions
 *
 * Exp[-|x|] i_l(x)
 * for l=0,1,...,lmax
 *
 * exceptions: GSL_EUNDRFLW
 */
GSL_FUN int gsl_sf_bessel_il_scaled_array(const int lmax, const double x, double * result_array);

/* Irregular scaled modified spherical Bessel function
 * Exp[x] k_0(x)
 *
 * x > 0.0
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
GSL_FUN int gsl_sf_bessel_k0_scaled_e(const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_k0_scaled(const double x);

/* Irregular modified spherical Bessel function
 * Exp[x] k_1(x)
 *
 * x > 0.0
 * exceptions: GSL_EDOM, GSL_EUNDRFLW, GSL_EOVRFLW
 */
GSL_FUN int gsl_sf_bessel_k1_scaled_e(const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_k1_scaled(const double x);

/* Irregular modified spherical Bessel function
 * Exp[x] k_2(x)
 *
 * x > 0.0
 * exceptions: GSL_EDOM, GSL_EUNDRFLW, GSL_EOVRFLW
 */
GSL_FUN int gsl_sf_bessel_k2_scaled_e(const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_k2_scaled(const double x);

/* Irregular modified spherical Bessel function
 * Exp[x] k_l[x]
 *
 * k_l(x) = Sqrt[Pi/(2x)] BesselK[l+1/2,x]
 *
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
GSL_FUN int gsl_sf_bessel_kl_scaled_e(int l, const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_kl_scaled(const int l, const double x);

/* Irregular scaled modified spherical Bessel function
 * Exp[x] k_l(x)
 *
 * for l=0,1,...,lmax
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
GSL_FUN int gsl_sf_bessel_kl_scaled_array(const int lmax, const double x, double * result_array);

/* Regular cylindrical Bessel function J_nu(x)
 *
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
GSL_FUN int gsl_sf_bessel_Jnu_e(const double nu, const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_Jnu(const double nu, const double x);

/* Irregular cylindrical Bessel function Y_nu(x)
 *
 * exceptions:
 */
GSL_FUN int gsl_sf_bessel_Ynu_e(double nu, double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_Ynu(const double nu, const double x);

/* Regular cylindrical Bessel function J_nu(x)
 * evaluated at a series of x values. The array
 * contains the x values. They are assumed to be
 * strictly ordered and positive. The array is
 * over-written with the values of J_nu(x_i).
 *
 * exceptions: GSL_EDOM, GSL_EINVAL
 */
GSL_FUN int gsl_sf_bessel_sequence_Jnu_e(double nu, gsl_mode_t mode, size_t size, double * v);

/* Scaled modified cylindrical Bessel functions
 *
 * Exp[-|x|] BesselI[nu, x]
 * x >= 0, nu >= 0
 *
 * exceptions: GSL_EDOM
 */
GSL_FUN int gsl_sf_bessel_Inu_scaled_e(double nu, double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_Inu_scaled(double nu, double x);

/* Modified cylindrical Bessel functions
 *
 * BesselI[nu, x]
 * x >= 0, nu >= 0
 *
 * exceptions: GSL_EDOM, GSL_EOVRFLW
 */
GSL_FUN int gsl_sf_bessel_Inu_e(double nu, double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_Inu(double nu, double x);

/* Scaled modified cylindrical Bessel functions
 *
 * Exp[+|x|] BesselK[nu, x]
 * x > 0, nu >= 0
 *
 * exceptions: GSL_EDOM
 */
GSL_FUN int gsl_sf_bessel_Knu_scaled_e(const double nu, const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_Knu_scaled(const double nu, const double x);

GSL_FUN int gsl_sf_bessel_Knu_scaled_e10_e(const double nu, const double x, gsl_sf_result_e10 * result);

/* Modified cylindrical Bessel functions
 *
 * BesselK[nu, x]
 * x > 0, nu >= 0
 *
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
GSL_FUN int gsl_sf_bessel_Knu_e(const double nu, const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_Knu(const double nu, const double x);

/* Logarithm of modified cylindrical Bessel functions.
 *
 * Log[BesselK[nu, x]]
 * x > 0, nu >= 0
 *
 * exceptions: GSL_EDOM
 */
GSL_FUN int gsl_sf_bessel_lnKnu_e(const double nu, const double x, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_lnKnu(const double nu, const double x);

/* s'th positive zero of the Bessel function J_0(x).
 *
 * exceptions:
 */
GSL_FUN int gsl_sf_bessel_zero_J0_e(unsigned int s, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_zero_J0(unsigned int s);

/* s'th positive zero of the Bessel function J_1(x).
 *
 * exceptions:
 */
GSL_FUN int gsl_sf_bessel_zero_J1_e(unsigned int s, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_zero_J1(unsigned int s);

/* s'th positive zero of the Bessel function J_nu(x).
 *
 * exceptions:
 */
GSL_FUN int gsl_sf_bessel_zero_Jnu_e(double nu, unsigned int s, gsl_sf_result * result);
GSL_FUN double gsl_sf_bessel_zero_Jnu(double nu, unsigned int s);

__END_DECLS

#endif /* __GSL_SF_BESSEL_H__ */

max_stars_repo_name: mgreter/astrometrylib

max_issues_repo_name: zzpwahaha/VimbaCamJILA

max_forks_repo_name: zzpwahaha/VimbaCamJILA


text:
/*
 * Copyright (c) 1997-1999 Massachusetts Institute of Technology
 *
 * 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
 *
 */

/* This file was automatically generated --- DO NOT EDIT */
/* Generated on Sun Nov  7 20:43:56 EST 1999 */

#include <fftw-int.h>
#include <fftw.h>

/* Generated by: ./genfft -magic-alignment-check -magic-twiddle-load-all -magic-variables 4 -magic-loopi -real2hc 13 */

/*
 * This function contains 76 FP additions, 34 FP multiplications,
 * (or, 57 additions, 15 multiplications, 19 fused multiply/add),
 * 34 stack variables, and 26 memory accesses
 */
static const fftw_real K083333333 = FFTW_KONST(+0.083333333333333333333333333333333333333333333);
static const fftw_real K075902986 = FFTW_KONST(+0.075902986037193865983102897245103540356428373);
static const fftw_real K251768516 = FFTW_KONST(+0.251768516431883313623436926934233488546674281);
static const fftw_real K503537032 = FFTW_KONST(+0.503537032863766627246873853868466977093348562);
static const fftw_real K113854479 = FFTW_KONST(+0.113854479055790798974654345867655310534642560);
static const fftw_real K265966249 = FFTW_KONST(+0.265966249214837287587521063842185948798330267);
static const fftw_real K387390585 = FFTW_KONST(+0.387390585467617292130675966426762851778775217);
static const fftw_real K300462606 = FFTW_KONST(+0.300462606288665774426601772289207995520941381);
static const fftw_real K258260390 = FFTW_KONST(+0.258260390311744861420450644284508567852516811);
static const fftw_real K132983124 = FFTW_KONST(+0.132983124607418643793760531921092974399165133);
static const fftw_real K2_000000000 = FFTW_KONST(+2.000000000000000000000000000000000000000000000);
static const fftw_real K1_732050807 = FFTW_KONST(+1.732050807568877293527446341505872366942805254);
static const fftw_real K156891391 = FFTW_KONST(+0.156891391051584611046832726756003269660212636);
static const fftw_real K256247671 = FFTW_KONST(+0.256247671582936600958684654061725059144125175);
static const fftw_real K011599105 = FFTW_KONST(+0.011599105605768290721655456654083252189827041);
static const fftw_real K300238635 = FFTW_KONST(+0.300238635966332641462884626667381504676006424);
static const fftw_real K174138601 = FFTW_KONST(+0.174138601152135905005660794929264742616964676);
static const fftw_real K575140729 = FFTW_KONST(+0.575140729474003121368385547455453388461001608);
static const fftw_real K866025403 = FFTW_KONST(+0.866025403784438646763723170752936183471402627);
static const fftw_real K500000000 = FFTW_KONST(+0.500000000000000000000000000000000000000000000);

/*
 * Generator Id's :
 * $Id: exprdag.ml,v 1.41 1999/05/26 15:44:14 fftw Exp $
 * $Id: fft.ml,v 1.43 1999/05/17 19:44:18 fftw Exp $
 * $Id: to_c.ml,v 1.25 1999/10/26 21:41:32 stevenj Exp $
 */

void fftw_real2hc_13(const fftw_real *input, fftw_real *real_output, fftw_real *imag_output, int istride, int real_ostride, int imag_ostride)
{
     fftw_real tmp65;
     fftw_real tmp11;
     fftw_real tmp28;
     fftw_real tmp37;
     fftw_real tmp51;
     fftw_real tmp62;
     fftw_real tmp22;
     fftw_real tmp58;
     fftw_real tmp59;
     fftw_real tmp66;
     fftw_real tmp56;
     fftw_real tmp63;
     fftw_real tmp35;
     fftw_real tmp38;
     ASSERT_ALIGNED_DOUBLE;
     tmp65 = input[0];
     {
	  fftw_real tmp3;
	  fftw_real tmp53;
	  fftw_real tmp21;
	  fftw_real tmp30;
	  fftw_real tmp26;
	  fftw_real tmp16;
	  fftw_real tmp29;
	  fftw_real tmp25;
	  fftw_real tmp6;
	  fftw_real tmp32;
	  fftw_real tmp9;
	  fftw_real tmp33;
	  fftw_real tmp10;
	  fftw_real tmp54;
	  fftw_real tmp1;
	  fftw_real tmp2;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp1 = input[8 * istride];
	  tmp2 = input[5 * istride];
	  tmp3 = tmp1 - tmp2;
	  tmp53 = tmp1 + tmp2;
	  {
	       fftw_real tmp17;
	       fftw_real tmp18;
	       fftw_real tmp19;
	       fftw_real tmp20;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp17 = input[12 * istride];
	       tmp18 = input[4 * istride];
	       tmp19 = input[10 * istride];
	       tmp20 = tmp18 + tmp19;
	       tmp21 = tmp17 + tmp20;
	       tmp30 = tmp17 - (K500000000 * tmp20);
	       tmp26 = tmp18 - tmp19;
	  }
	  {
	       fftw_real tmp12;
	       fftw_real tmp13;
	       fftw_real tmp14;
	       fftw_real tmp15;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp12 = input[istride];
	       tmp13 = input[3 * istride];
	       tmp14 = input[9 * istride];
	       tmp15 = tmp13 + tmp14;
	       tmp16 = tmp12 + tmp15;
	       tmp29 = tmp12 - (K500000000 * tmp15);
	       tmp25 = tmp13 - tmp14;
	  }
	  {
	       fftw_real tmp4;
	       fftw_real tmp5;
	       fftw_real tmp7;
	       fftw_real tmp8;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp4 = input[6 * istride];
	       tmp5 = input[11 * istride];
	       tmp6 = tmp4 - tmp5;
	       tmp32 = tmp4 + tmp5;
	       tmp7 = input[2 * istride];
	       tmp8 = input[7 * istride];
	       tmp9 = tmp7 - tmp8;
	       tmp33 = tmp7 + tmp8;
	  }
	  tmp10 = tmp6 + tmp9;
	  tmp54 = tmp32 + tmp33;
	  tmp11 = tmp3 - tmp10;
	  {
	       fftw_real tmp24;
	       fftw_real tmp27;
	       fftw_real tmp49;
	       fftw_real tmp50;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp24 = tmp3 + (K500000000 * tmp10);
	       tmp27 = K866025403 * (tmp25 + tmp26);
	       tmp28 = tmp24 - tmp27;
	       tmp37 = tmp27 + tmp24;
	       tmp49 = tmp9 - tmp6;
	       tmp50 = tmp25 - tmp26;
	       tmp51 = tmp49 - tmp50;
	       tmp62 = tmp50 + tmp49;
	  }
	  tmp22 = tmp16 - tmp21;
	  tmp58 = tmp16 + tmp21;
	  tmp59 = tmp53 + tmp54;
	  tmp66 = tmp58 + tmp59;
	  {
	       fftw_real tmp52;
	       fftw_real tmp55;
	       fftw_real tmp31;
	       fftw_real tmp34;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp52 = tmp29 + tmp30;
	       tmp55 = tmp53 - (K500000000 * tmp54);
	       tmp56 = tmp52 - tmp55;
	       tmp63 = tmp52 + tmp55;
	       tmp31 = tmp29 - tmp30;
	       tmp34 = K866025403 * (tmp32 - tmp33);
	       tmp35 = tmp31 + tmp34;
	       tmp38 = tmp31 - tmp34;
	  }
     }
     real_output[0] = tmp65 + tmp66;
     {
	  fftw_real tmp23;
	  fftw_real tmp45;
	  fftw_real tmp40;
	  fftw_real tmp48;
	  fftw_real tmp44;
	  fftw_real tmp46;
	  fftw_real tmp41;
	  fftw_real tmp47;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp23 = (K575140729 * tmp11) - (K174138601 * tmp22);
	  tmp45 = (K575140729 * tmp22) + (K174138601 * tmp11);
	  {
	       fftw_real tmp36;
	       fftw_real tmp39;
	       fftw_real tmp42;
	       fftw_real tmp43;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp36 = (K300238635 * tmp28) + (K011599105 * tmp35);
	       tmp39 = (K256247671 * tmp37) + (K156891391 * tmp38);
	       tmp40 = tmp36 - tmp39;
	       tmp48 = K1_732050807 * (tmp39 + tmp36);
	       tmp42 = (K300238635 * tmp35) - (K011599105 * tmp28);
	       tmp43 = (K156891391 * tmp37) - (K256247671 * tmp38);
	       tmp44 = K1_732050807 * (tmp42 - tmp43);
	       tmp46 = tmp43 + tmp42;
	  }
	  imag_output[imag_ostride] = tmp23 + (K2_000000000 * tmp40);
	  tmp41 = tmp23 - tmp40;
	  imag_output[3 * imag_ostride] = tmp41 - tmp44;
	  imag_output[4 * imag_ostride] = -(tmp41 + tmp44);
	  imag_output[5 * imag_ostride] = -(tmp45 + (K2_000000000 * tmp46));
	  tmp47 = tmp46 - tmp45;
	  imag_output[2 * imag_ostride] = tmp47 - tmp48;
	  imag_output[6 * imag_ostride] = tmp48 + tmp47;
     }
     {
	  fftw_real tmp61;
	  fftw_real tmp70;
	  fftw_real tmp74;
	  fftw_real tmp76;
	  fftw_real tmp68;
	  fftw_real tmp69;
	  fftw_real tmp57;
	  fftw_real tmp60;
	  fftw_real tmp71;
	  fftw_real tmp75;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp57 = (K132983124 * tmp51) + (K258260390 * tmp56);
	  tmp60 = K300462606 * (tmp58 - tmp59);
	  tmp61 = (K2_000000000 * tmp57) + tmp60;
	  tmp70 = tmp60 - tmp57;
	  {
	       fftw_real tmp72;
	       fftw_real tmp73;
	       fftw_real tmp64;
	       fftw_real tmp67;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp72 = (K387390585 * tmp51) - (K265966249 * tmp56);
	       tmp73 = (K113854479 * tmp62) - (K503537032 * tmp63);
	       tmp74 = tmp72 + tmp73;
	       tmp76 = tmp73 - tmp72;
	       tmp64 = (K251768516 * tmp62) + (K075902986 * tmp63);
	       tmp67 = tmp65 - (K083333333 * tmp66);
	       tmp68 = (K2_000000000 * tmp64) + tmp67;
	       tmp69 = tmp67 - tmp64;
	  }
	  real_output[real_ostride] = tmp61 + tmp68;
	  real_output[5 * real_ostride] = tmp68 - tmp61;
	  tmp71 = tmp69 - tmp70;
	  real_output[2 * real_ostride] = tmp71 - tmp74;
	  real_output[6 * real_ostride] = tmp74 + tmp71;
	  tmp75 = tmp70 + tmp69;
	  real_output[3 * real_ostride] = tmp75 - tmp76;
	  real_output[4 * real_ostride] = tmp76 + tmp75;
     }
}

fftw_codelet_desc fftw_real2hc_13_desc =
{
     "fftw_real2hc_13",
     (void (*)()) fftw_real2hc_13,
     13,
     FFTW_FORWARD,
     FFTW_REAL2HC,
     288,
     0,
     (const int *) 0,
};

ext: c

max_stars_repo_name: albertsgrc/ftdock-opt

max_issues_repo_name: albertsgrc/ftdock-opt

max_forks_repo_name: albertsgrc/ftdock-opt


text:
// --------------------------------------------------------------------------
//                   OpenMS -- Open-Source Mass Spectrometry
// --------------------------------------------------------------------------
// Copyright The OpenMS Team -- Eberhard Karls University Tuebingen,
// ETH Zurich, and Freie Universitaet Berlin 2002-2013.
//
// This software is released under a three-clause BSD license:
//  * Redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//  * Redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
//  * Neither the name of any author or any participating institution
//    may be used to endorse or promote products derived from this software
//    without specific prior written permission.
// For a full list of authors, refer to the file AUTHORS.
// --------------------------------------------------------------------------
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL ANY OF THE AUTHORS OR THE CONTRIBUTING
// INSTITUTIONS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
// OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
// WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
// ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// --------------------------------------------------------------------------
// $Maintainer: Clemens Groepl $
// $Authors: $
// --------------------------------------------------------------------------

#ifndef OPENMS_MATH_STATISTICS_LINEARREGRESSION_H
#define OPENMS_MATH_STATISTICS_LINEARREGRESSION_H

#include <OpenMS/CONCEPT/Types.h>
#include <OpenMS/CONCEPT/Exception.h>

#include <vector>
#include <gsl/gsl_fit.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_cdf.h>

namespace OpenMS
{
  namespace Math
  {
    /**
      @brief This class offers functions to perform least-squares fits to a straight line model, \f$ Y(c,x) = c_0 + c_1 x \f$.

      It encapsulates the GSL methods for a weighted and an unweighted linear regression.

      Next to the intercept with the y-axis and the slope of the fitted line, this class computes the:
        - squared pearson coefficient
        - value of the t-distribution
        - standard deviation of the residuals
        - standard error of the slope
        - intercept with the x-axis (useful for additive series experiments)
        - lower border of confidence interval
        - higher border of confidence interval
        - chi squared value
        - x mean

      @ingroup Math
    */
    class OPENMS_DLLAPI LinearRegression
    {
public:

      /// Constructor
      LinearRegression() :
        intercept_(0),
        slope_(0),
        x_intercept_(0),
        lower_(0),
        upper_(0),
        t_star_(0),
        r_squared_(0),
        stand_dev_residuals_(0),
        mean_residuals_(0),
        stand_error_slope_(0),
        chi_squared_(0),
        rsd_(0)
      {
      }

      /// Destructor
      virtual ~LinearRegression()
      {
      }

      /**
          @brief This function computes the best-fit linear regression coefficients \f$ (c_0,c_1) \f$
          of the model \f$ Y = c_0 + c_1 X \f$ for the dataset \f$ (x, y) \f$.

          The values in x-dimension of the dataset \f$ (x,y) \f$ are given by the iterator range [x_begin,x_end)
          and the corresponding y-values start at position y_begin.

          For a  "x %" Confidence Interval use confidence_interval_P = x/100.
          For example the 95% Confidence Interval is supposed to be an interval that has a 95% chance of
          containing the true value of the parameter.

          @return If an error occurred during the fit.

          @exception Exception::UnableToFit is thrown if fitting cannot be performed
      */
      template <typename Iterator>
      void computeRegression(double confidence_interval_P, Iterator x_begin, Iterator x_end, Iterator y_begin);

      /**
          @brief This function computes the best-fit linear regression coefficient \f$ (c_0) \f$
          of the model \f$ Y = c_1 X \f$ for the dataset \f$ (x, y) \f$.

          The values in x-dimension of the dataset \f$ (x,y) \f$ are given by the iterator range [x_begin,x_end)
          and the corresponding y-values start at position y_begin.

          For a  "x %" Confidence Interval use confidence_interval_P = x/100.
          For example the 95% Confidence Interval is supposed to be an interval that has a 95% chance of
          containing the true value of the parameter.

          @return If an error occurred during the fit.

          @exception Exception::UnableToFit is thrown if fitting cannot be performed
      */
      template <typename Iterator>
      void computeRegressionNoIntercept(double confidence_interval_P, Iterator x_begin, Iterator x_end, Iterator y_begin);

      /**
          @brief This function computes the best-fit linear regression coefficients \f$ (c_0,c_1) \f$
          of the model \f$ Y = c_0 + c_1 X \f$ for the weighted dataset \f$ (x, y) \f$.

          The values in x-dimension of the dataset \f$ (x, y) \f$ are given by the iterator range [x_begin,x_end)
          and the corresponding y-values start at position y_begin. They will be weighted by the
          values starting at w_begin.

          For a  "x %" Confidence Interval use confidence_interval_P = x/100.
          For example the 95% Confidence Interval is supposed to be an interval that has a 95% chance of
          containing the true value of the parameter.

          @return If an error occurred during the fit.

          @exception Exception::UnableToFit is thrown if fitting cannot be performed
      */
      template <typename Iterator>
      void computeRegressionWeighted(double confidence_interval_P, Iterator x_begin, Iterator x_end, Iterator y_begin, Iterator w_begin);

      /// Non-mutable access to the y-intercept of the straight line
      DoubleReal getIntercept() const;
      /// Non-mutable access to the slope of the straight line
      DoubleReal getSlope() const;
      /// Non-mutable access to the x-intercept of the straight line
      DoubleReal getXIntercept() const;
      /// Non-mutable access to the lower border of confidence interval
      DoubleReal getLower() const;
      /// Non-mutable access to the upper border of confidence interval
      DoubleReal getUpper() const;
      /// Non-mutable access to the value of the t-distribution
      DoubleReal getTValue() const;
      /// Non-mutable access to the squared pearson coefficient
      DoubleReal getRSquared() const;
      /// Non-mutable access to the standard deviation of the residuals
      DoubleReal getStandDevRes() const;
      /// Non-mutable access to the residual mean
      DoubleReal getMeanRes() const;
      /// Non-mutable access to the standard error of the slope
      DoubleReal getStandErrSlope() const;
      /// Non-mutable access to the chi squared value
      DoubleReal getChiSquared() const;
      /// Non-mutable access to relative standard deviation
      DoubleReal getRSD() const;

protected:

      /// The intercept of the fitted line with the y-axis
      double intercept_;
      /// The slope of the fitted line
      double slope_;
      /// The intercept of the fitted line with the x-axis
      double x_intercept_;
      /// The lower bound of the confidence interval
      double lower_;
      /// The upper bound of the confidence interval
      double upper_;
      /// The value of the t-statistic
      double t_star_;
      /// The squared correlation coefficient (Pearson)
      double r_squared_;
      /// The standard deviation of the residuals
      double stand_dev_residuals_;
      /// Mean of residuals
      double mean_residuals_;
      /// The standard error of the slope
      double stand_error_slope_;
      /// The value of the Chi Squared statistic
      double chi_squared_;
      /// the relative standard deviation
      double rsd_;

      /// Computes the goodness of the fitted regression line
      void computeGoodness_(double * X, double * Y, int N, double confidence_interval_P);

      /// Copies the distance(x_begin,x_end) elements starting at x_begin and y_begin into the arrays x_array and y_array
      template <typename Iterator>
      void iteratorRange2Arrays_(Iterator x_begin, Iterator x_end, Iterator y_begin, double * x_array, double * y_array);

      /// Copy the distance(x_begin,x_end) elements starting at  x_begin, y_begin and w_begin into the arrays x_array, y_array and w_array
      template <typename Iterator>
      void iteratorRange3Arrays_(Iterator x_begin, Iterator x_end, Iterator y_begin, Iterator w_begin, double * x_array, double * y_array, double * w_array);

private:

      /// Not implemented
      LinearRegression(const LinearRegression & arg);
      /// Not implemented
      LinearRegression & operator=(const LinearRegression & arg);
    };

    template <typename Iterator>
    void LinearRegression::computeRegression(double confidence_interval_P, Iterator x_begin, Iterator x_end, Iterator y_begin)
    {
      int N = int(distance(x_begin, x_end));

      double * X = new double[N];
      double * Y = new double[N];
      iteratorRange2Arrays_(x_begin, x_end, y_begin, X, Y);

      double cov00, cov01, cov11;

      // Compute the unweighted linear fit.
      // Get the intercept and the slope of the regression Y_hat=intercept_+slope_*X
      // and the value of Chi squared, the covariances of the intercept and the slope
      int error = gsl_fit_linear(X, 1, Y, 1, N, &intercept_, &slope_, &cov00, &cov01, &cov11, &chi_squared_);

      if (!error)
      {
        computeGoodness_(X, Y, N, confidence_interval_P);
      }

      delete[] X;
      delete[] Y;

      if (error)
      {
        throw Exception::UnableToFit(__FILE__, __LINE__, __PRETTY_FUNCTION__, "UnableToFit-LinearRegression", "Could not fit a linear model to the data");
      }
    }

    template <typename Iterator>
    void LinearRegression::computeRegressionNoIntercept(double confidence_interval_P, Iterator x_begin, Iterator x_end, Iterator y_begin)
    {
      int N = int(distance(x_begin, x_end));

      double * X = new double[N];
      double * Y = new double[N];
      iteratorRange2Arrays_(x_begin, x_end, y_begin, X, Y);

      double cov;

      // Compute the linear fit.
      // Get the intercept and the slope of the regression Y_hat=intercept_+slope_*X
      // and the value of Chi squared, the covariances of the intercept and the slope
      int error = gsl_fit_mul(X, 1, Y, 1, N, &slope_, &cov, &chi_squared_);
      intercept_ = 0.0;

      if (!error)
      {
        computeGoodness_(X, Y, N, confidence_interval_P);
      }

      delete[] X;
      delete[] Y;

      if (error)
      {
        throw Exception::UnableToFit(__FILE__, __LINE__, __PRETTY_FUNCTION__, "UnableToFit-LinearRegression", "Could not fit a linear model to the data");
      }
    }

    template <typename Iterator>
    void LinearRegression::computeRegressionWeighted(double confidence_interval_P, Iterator x_begin, Iterator x_end, Iterator y_begin, Iterator w_begin)
    {
      int N = int(distance(x_begin, x_end));

      double * X = new double[N];
      double * Y = new double[N];
      double * W = new double[N];
      iteratorRange3Arrays_(x_begin, x_end, y_begin, w_begin, X, Y, W);

      double cov00, cov01, cov11;

      // Compute the weighted linear fit.
      // Get the intercept and the slope of the regression Y_hat=intercept_+slope_*X
      // and the value of Chi squared, the covariances of the intercept and the slope
      int error = gsl_fit_wlinear(X, 1, W, 1, Y, 1, N, &intercept_, &slope_, &cov00, &cov01, &cov11, &chi_squared_);

      if (!error)
      {
        computeGoodness_(X, Y, N, confidence_interval_P);
      }

      delete[] X;
      delete[] Y;
      delete[] W;

      if (error)
      {
        throw Exception::UnableToFit(__FILE__, __LINE__, __PRETTY_FUNCTION__, "UnableToFit-LinearRegression", "Could not fit a linear model to the data");
      }
    }

    template <typename Iterator>
    void LinearRegression::iteratorRange2Arrays_(Iterator x_begin, Iterator x_end, Iterator y_begin, double * x_array, double * y_array)
    {
      int i = 0;
      while (x_begin < x_end)
      {
        x_array[i] = *x_begin;
        y_array[i] = *y_begin;
        ++x_begin;
        ++y_begin;
        ++i;
      }
    }

    template <typename Iterator>
    void LinearRegression::iteratorRange3Arrays_(Iterator x_begin, Iterator x_end, Iterator y_begin, Iterator w_begin, double * x_array, double * y_array, double * w_array)
    {
      int i = 0;
      while (x_begin < x_end)
      {
        x_array[i] = *x_begin;
        y_array[i] = *y_begin;
        w_array[i] = *w_begin;
        ++x_begin;
        ++y_begin;
        ++w_begin;
        ++i;
      }
    }

  } // namespace Math
} // namespace OpenMS

#endif

ext: h

max_stars_repo_name: kreinert/OpenMS

max_issues_repo_name: kreinert/OpenMS

max_forks_repo_name: kreinert/OpenMS


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

#include <config.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_diff.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_test.h>

double
f1 (double x, void *params)
{
  return exp (x);
}

double
df1 (double x, void *params)
{
  return exp (x);
}

double
f2 (double x, void *params)
{
  if (x >= 0.0)
    {
      return x * sqrt (x);
    }
  else
    {
      return 0.0;
    }
}

double
df2 (double x, void *params)
{
  if (x >= 0.0)
    {
      return 1.5 * sqrt (x);
    }
  else
    {
      return 0.0;
    }
}

double
f3 (double x, void *params)
{
  if (x != 0.0)
    {
      return sin (1 / x);
    }
  else
    {
      return 0.0;
    }
}

double
df3 (double x, void *params)
{
  if (x != 0.0)
    {
      return -cos (1 / x) / (x * x);
    }
  else
    {
      return 0.0;
    }
}

double
f4 (double x, void *params)
{
  return exp (-x * x);
}

double
df4 (double x, void *params)
{
  return -2.0 * x * exp (-x * x);
}

double
f5 (double x, void *params)
{
  return x * x;
}

double
df5 (double x, void *params)
{
  return 2.0 * x;
}

double
f6 (double x, void *params)
{
  return 1.0 / x;
}

double
df6 (double x, void *params)
{
  return -1.0 / (x * x);
}

typedef int (diff_fn) (const gsl_function * f, double x, double * res, double *abserr);

void
test (diff_fn * diff, gsl_function * f, gsl_function * df, double x,
      const char * desc)
{
  double result, abserr;
  double expected = GSL_FN_EVAL (df, x);
  (*diff) (f, x, &result, &abserr);
  gsl_test_abs (result, expected, abserr, desc);
  gsl_test (fabs(result-expected) >  abserr, "%s, valid error estimate", desc);
}

int
main ()
{
  gsl_function F1, DF1, F2, DF2, F3, DF3, F4, DF4, F5, DF5, F6, DF6;

  F1.function = &f1;
  DF1.function = &df1;

  F2.function = &f2;
  DF2.function = &df2;

  F3.function = &f3;
  DF3.function = &df3;

  F4.function = &f4;
  DF4.function = &df4;

  F5.function = &f5;
  DF5.function = &df5;

  F6.function = &f6;
  DF6.function = &df6;

  test (&gsl_diff_central, &F1, &DF1, 1.0, "exp(x), x=1, central diff");
  test (&gsl_diff_forward, &F1, &DF1, 1.0, "exp(x), x=1, forward diff");
  test (&gsl_diff_backward, &F1, &DF1, 1.0, "exp(x), x=1, backward diff");

  test (&gsl_diff_central, &F2, &DF2, 0.1, "x^(3/2), x=0.1, central diff");
  test (&gsl_diff_forward, &F2, &DF2, 0.1, "x^(3/2), x=0.1, forward diff");
  test (&gsl_diff_backward, &F2, &DF2, 0.1, "x^(3/2), x=0.1, backward diff");

  test (&gsl_diff_central, &F3, &DF3, 0.45, "sin(1/x), x=0.45, central diff");
  test (&gsl_diff_forward, &F3, &DF3, 0.45, "sin(1/x), x=0.45, forward diff");
  test (&gsl_diff_backward, &F3, &DF3, 0.45, "sin(1/x), x=0.45, backward diff");

  test (&gsl_diff_central, &F4, &DF4, 0.5, "exp(-x^2), x=0.5, central diff");
  test (&gsl_diff_forward, &F4, &DF4, 0.5, "exp(-x^2), x=0.5, forward diff");
  test (&gsl_diff_backward, &F4, &DF4, 0.5, "exp(-x^2), x=0.5, backward diff");

  test (&gsl_diff_central, &F5, &DF5, 0.0, "x^2, x=0, central diff");
  test (&gsl_diff_forward, &F5, &DF5, 0.0, "x^2, x=0, forward diff");
  test (&gsl_diff_backward, &F5, &DF5, 0.0, "x^2, x=0, backward diff");

  test (&gsl_diff_central, &F6, &DF6, 10.0, "1/x, x=10, central diff");
  test (&gsl_diff_forward, &F6, &DF6, 10.0, "1/x, x=10, forward diff");
  test (&gsl_diff_backward, &F6, &DF6, 10.0, "1/x, x=10, backward diff");

  exit (gsl_test_summary ());
}


text:
/*
 * $Source: /Users/mariko/lib/source/sde_wa/sde_wa_manual/sde_wa_manual_progra
 * $Revision: 1.7 $
 * $Author: mariko $
 * $Date: 2011/11/29 10:48:07 $
 */
#include <stdlib.h> /* for malloc */
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_qrng.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_statistics.h>
#include <sde_wa.h>
#include "european_option_blacksholes.h"
#include "mlmc.h"

#define L_MAX 20

int main(void){
	enum ALG alg = E_M;
	int m_moment = 7; /* Romberg */
	int M = 4;
	int n;
	double x0=1.0;
	double K = 1.00;
	struct EBS_params ebs_params;
	double mu=0.05;
	double sigma=0.2;
	double T=1.0;
	ebs_params.mu=mu;
	ebs_params.sigma=sigma;
	ebs_params.T=T;
	ebs_params.K=K;
	MLMC_L **mlmc_est;
	{
		SDE_WA_SYSTEM *sde;
		SDE_WA_SLTN *sl;
		SDE_WA_PRICER *pricer;
		double sum;
		int dNl;
		double hl;
		double eps=0.0001;
		double epss[] = {0.001, 0.0005, 0.0002, 0.0001, 0.00005};
		int l;
		int L;
		int i;
		double con;
		int converged_flag=0;
		int extraporation_flag=0;
		sde=alloc_SDE_WA_SYSTEM(1, 1, &ebs_params);
		sde->sde_type=STR;

		sde->V[0]=ebs_str_V_0;
		sde->V[1]=ebs_V_1;

		sde->drift_corrector[0]=NULL;
		sde->drift_corrector[1]=diff_ebs_V_1; /* for EM */

		sde->exp_sV[0]=exp_ebs_sVy_0; /* for NV */
		sde->exp_sV[1]=exp_ebs_sVy_1; /* for NV */

		sl=alloc_SDE_WA_SLTN(alg, m_moment, sde);

		pricer=alloc_SDE_WA_PRICER(sl);
		pricer->payoff_function=european_option_blacksholes_call_payoff;
		pricer->init_y[0] = x0;
		pricer->T = 1.0;

		//gsl_qrng *q;
		gsl_rng_type *rng_type = (gsl_rng_type *)gsl_rng_default;
		gsl_rng *rng = gsl_rng_alloc(rng_type);

		if (pricer->gen_unif_rand == NULL) {
			/* sobol seq. => seq.dim.<=40*/
		//	q=gsl_qrng_alloc(gsl_qrng_sobol, n*sde->dim_BM);
		//	pricer->gen_unif_rand = q;
		//	pricer->next_rand = gen_unif_rand_number_gsl_quasi;

			pricer->gen_unif_rand = rng;
			pricer->next_rand = gen_unif_rand_number_gsl_pseudo;
		}

		mlmc_est = (MLMC_L **)malloc(sizeof(MLMC_L *)*L_MAX);
		for (l=0; l<L_MAX; l++) {
			mlmc_est[l] = alloc_level_l_estimator();
		}

		for (l=0; l<=-4; l++) {
			L=l;
			mlmc_est[L]->N = 4000000;
			mlmc_level_l_estimation(pricer, mlmc_est[L], L, M, mlmc_est[L]->N);
			mlmc_est[L]->V = mlmc_est[L]->sumV/mlmc_est[L]->N - (mlmc_est[L]->sumY/mlmc_est[L]->N)*(mlmc_est[L]->sumY/mlmc_est[L]->N);
			printf("L:%d M:%d\n", L, M);
			printf("N:%d, Y:%12e, V:%12.12lf, sumV:%lf\n", mlmc_est[L]->N, mlmc_est[L]->sumY/mlmc_est[L]->N, mlmc_est[L]->V, mlmc_est[L]->sumV);
		}
		//exit(-1);
		for (i=0; i<5; i++) {
			L=-1;
			eps = epss[i];
			for (l=0; l<L_MAX; l++) {
				init_level_l_estimator(mlmc_est[l]);
			}
			converged_flag=0;
			printf("eps:%lf \n", eps);
			while (!converged_flag) {
				L++;
				n = (int)pow(M,L);

				mlmc_level_l_estimation(pricer, mlmc_est[L], L, M, mlmc_est[L]->N);
				printf("N:%d, Y:%lf, sumV:%lf\n", mlmc_est[L]->N, mlmc_est[L]->sumY/mlmc_est[L]->N, mlmc_est[L]->sumV);
				//estimate Vl
				sum = 0.0;
				for (l=0; l<=L; l++) {
					mlmc_est[l]->V = mlmc_est[l]->sumV/mlmc_est[l]->N -
						(mlmc_est[l]->sumY/mlmc_est[l]->N)*(mlmc_est[l]->sumY/mlmc_est[l]->N);
					hl = 1.0/pow(M,l);
					sum += sqrt(mlmc_est[l]->V/hl);
				}

				//define optimal Nl
				for (l=0; l<=L; l++) {
					mlmc_est[l]->newN = ceil(2.0 * sqrt(mlmc_est[l]->V * hl) * sum / (eps * eps));
					printf("newN[%d]:%d\n", l, mlmc_est[l]->newN);
				}

				//update sample sum
				for (l=0; l<=L; l++) {
					dNl = mlmc_est[l]->newN - mlmc_est[l]->N;
					printf("dN[%d]:%d\n", l, dNl);
					if (dNl > 0) {
						mlmc_est[l]->N += dNl;
						mlmc_level_l_estimation(pricer, mlmc_est[l], l, M, dNl);
					}
				}

				//test for convergence
				if (extraporation_flag) {
					if (L > 1) {
						con = mlmc_est[L]->sumY/mlmc_est[L]->N - (1.0/(double)n) * mlmc_est[L-1]->sumY/mlmc_est[L-1]->N;
						if (fabs(con) < 1.0*(M*M - 1.0)*eps/sqrt(2.0)) {
							converged_flag = 1;
						}
					}
				} else {
					if (L > 1) {
						con = fmax(fabs(mlmc_est[L]->sumY/mlmc_est[L]->N), fabs(mlmc_est[L-1]->sumY/mlmc_est[L-1]->N/M));
						if (con < (M-1)*eps/sqrt(2.0)) {
							converged_flag = 1;
						}
					}
				}
				printf("N:%d, Y:%lf, V:%12.12lf, sumV:%lf\n", mlmc_est[L]->N, mlmc_est[L]->sumY/mlmc_est[L]->N, mlmc_est[L]->V, mlmc_est[L]->sumV);
				printf("L:%d, con:%lf, test:%lf \n\n", L, con, (M-1)*eps/sqrt(2.0));
			}
			sum=0.0;
			for (l=0; l<=L; l++) {
				sum += mlmc_est[l]->sumY / mlmc_est[l]->N;
			}
			if (extraporation_flag) {
				sum += mlmc_est[L]->sumY/mlmc_est[L]->N/(M-1);
			}
			printf("P:%.12e\n\n", sum);
		}

		for (l=0; l<L_MAX; l++) {
			free(mlmc_est[l]);
		}
		free_SDE_WA_PRICER(pricer);
		free_SDE_WA_SLTN(sl);
		free_SDE_WA_SYSTEM(sde);
	}

	return SDE_WA_SUCCESS;
}

/**
 * MLMC_L allocation
 */
MLMC_L *alloc_level_l_estimator(void) {
	MLMC_L *mlmc_l;

	mlmc_l = (MLMC_L *)malloc(sizeof(MLMC_L));
	mlmc_l->sumY = 0.0;
	mlmc_l->sumV = 0.0;
	mlmc_l->V = 0.0;
	mlmc_l->N = 10000;
	mlmc_l->newN = 0;

	return mlmc_l;
}

void init_level_l_estimator(MLMC_L *mlmc_l) {
	mlmc_l->sumY = 0.0;
	mlmc_l->sumV = 0.0;
	mlmc_l->V = 0.0;
	mlmc_l->N = 10000;
	mlmc_l->newN = 0;
}

void gen_unif_rand_number_gsl_pseudo(void *gen, double *rand, int n) {
	gsl_rng *rng = gen;
	int j=0;

	for (j=0; j<n; j++) {
		rand[j] = gsl_rng_uniform(rng);
	}
}

void gen_unif_rand_number_gsl_quasi(void *gen, double *rand, int n) {
	gsl_qrng *q = gen;

	gsl_qrng_get(q, rand);
}

/**
 * transform uniform random number to normal random number.
 * transform uniform random number to normal random number.
 * @param unif_rand n dimension uniform random number.
 * @param normal_rand n dimension normal random number.
 * @param n dimension of randm number. n must be even number.
 * @return transformed noraml random number
 */
void trans_rand_unif_to_normal(double *unif_rand, double *normal_rand, int n) {
	int i;
	double *u_seq1, *u_seq2;
	double *n_seq1, *n_seq2;

	if ((n % 2) != 0) {
		printf("%d error\n", n);
		return;
	}
	u_seq1 = unif_rand;
	u_seq2 = u_seq1 + n/2;
	n_seq1 = normal_rand;
	n_seq2 = n_seq1 + n/2;
	for (i=0; i<n/2; i++) {
		n_seq1[i] = sqrt(-2.0*log(u_seq1[i]))*cos(2.0*M_PIl*u_seq2[i]);
		n_seq2[i] = sqrt(-2.0*log(u_seq1[i]))*sin(2.0*M_PIl*u_seq2[i]);
	}
}

/**
 * level l estimaton
 * level l esitmator for MLMC
 * @param sl
 * @param sde
 * @param l level
 * @param M base of number of partitions
 * @param N num of samles for level l estimation
 * @return estimated expectation
 */
int mlmc_level_l_estimation(SDE_WA_PRICER *pricer, MLMC_L *estimator, int l, int M, int N) {
	SDE_WA_SLTN *sl = pricer->sltn;
	SDE_WA_SYSTEM *sde = sl->sde;
	int i,j,k,m;
	double sum=0.0;
	int n = pow(M, l);
	double dt1, dt2;
	double *u_seq;
	double *n_seq;
	double *sp, *sp2;
	double *tmp_pt;
	double tmp;
	double **xf;
	double **xc;

	u_seq=(double *)malloc(sizeof(double)*n*sde->dim_BM);
	n_seq=(double *)malloc(sizeof(double)*n*sde->dim_BM);

	xf=(double **)malloc(sizeof(double *)*2);
	for (i=0; i<2; i++) xf[i]=(double *)malloc(sizeof(double)*sde->dim_y);

	xc=(double **)malloc(sizeof(double *)*2);
	for (i=0; i<2; i++) xc[i]=(double *)malloc(sizeof(double)*sde->dim_y);

	if (l == 0) {
		sp=(double *)malloc(sizeof(double)*sde->dim_BM*2);
		sp2=(double *)malloc(sizeof(double)*sde->dim_BM*2);
		for (sum=0.0, i=0; i < N; i++){
			for (j=0; j<sde->dim_y; j++) {
				xf[0][j] = pricer->init_y[j];
			}
			pricer->next_rand(pricer->gen_unif_rand, sp2, sde->dim_BM*2);
			trans_rand_unif_to_normal(sp2, sp, sde->dim_BM*2);

			next_SDE_WA(sl, pricer->T, xf[0], xf[1], sp);
			tmp_pt=xf[0];
			xf[0]=xf[1];
			xf[1]=tmp_pt;

			estimator->sumY += pricer->payoff_function(xf[0], sde->params);
			estimator->sumV += pricer->payoff_function(xf[0], sde->params) * pricer->payoff_function(xf[0], sde->params);
		} /* for i */
	} else {
		sp=(double *)malloc(sizeof(double)*sde->dim_BM);
		sp2=(double *)malloc(sizeof(double)*sde->dim_BM);
		dt1=pricer->T/(double)n;
		dt2=M/(double)n;
		for (sum=0.0, i=0; i < N; i++){
			pricer->next_rand(pricer->gen_unif_rand, u_seq, n*sde->dim_BM);
			trans_rand_unif_to_normal(u_seq, n_seq, n*sde->dim_BM);

			for (j=0; j<sde->dim_y; j++) {
				xf[0][j] = pricer->init_y[j];
				xc[0][j] = pricer->init_y[j];
			}
			for (k=0; k<n/M; k++){
				for (m=0; m<sde->dim_BM; m++) {
					sp2[m] = 0.0;
				}
				for (j=0; j<M; j++){
					for (m=0; m<sde->dim_BM; m++) {
						sp[m] = n_seq[k*M + m + sde->dim_BM*j];
						sp2[m] += sp[m];
					}
					next_SDE_WA(sl, dt1, xf[0], xf[1], sp);
					tmp_pt=xf[0];
					xf[0]=xf[1];
					xf[1]=tmp_pt;
				}
				for (m=0; m<sde->dim_BM; m++) {
					sp2[m] = sp2[m]/sqrt(M);
				}
				next_SDE_WA(sl, dt2, xc[0], xc[1], sp2);
				tmp_pt=xc[0];
				xc[0]=xc[1];
				xc[1]=tmp_pt;
			}
			tmp = pricer->payoff_function(xf[0], sde->params) - pricer->payoff_function(xc[0], sde->params);
			estimator->sumY += tmp;
			estimator->sumV += tmp*tmp;
		} /* for i */
	}
	return sum;
}


text:
/**
@file

Copyright John Reid 2007

*/

#ifndef MYRRH_MATH_H_
#define MYRRH_MATH_H_

#ifdef _MSC_VER
# pragma once
#endif //_MSC_VER

#include "myrrh/defs.h"

#include <boost/iterator.hpp>
#include <boost/range.hpp>
#include <boost/tuple/tuple.hpp>
#include <boost/foreach.hpp>
#include <boost/shared_ptr.hpp>
#include <boost/math/special_functions/fpclassify.hpp>

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

#include <numeric>
#include <iostream>
#include <stdexcept>
#include <vector>

#include <math.h>

#include <boost/version.hpp>
#include <boost/test/floating_point_comparison.hpp>

#if BOOST_VERSION < 104400
# define MYRRH_FPC_NS ::boost::test_tools ///< Floating point comparison namespace
#elif BOOST_VERSION < 104500
# define MYRRH_FPC_NS ::boost::math::fpc ///< Floating point comparison namespace
#elif BOOST_VERSION < 104600
# define MYRRH_FPC_NS ::boost::test_tools ///< Floating point comparison namespace
#elif BOOST_VERSION < 104900
# define MYRRH_FPC_NS ::boost::math::fpc ///< Floating point comparison namespace
#elif BOOST_VERSION >= 105900
# include <boost/test/tools/old/impl.hpp>
# include <boost/test/tools/floating_point_comparison.hpp>
# define MYRRH_FPC_NS ::boost::math::fpc ///< Floating point comparison namespace
#else
# define MYRRH_FPC_NS ::boost::test_tools ///< Floating point comparison namespace
#endif

#ifdef _MSC_VER
# define MYRRH_ISNAN( x ) _isnan( x )
# define MYRRH_ISFINITE( x ) _finite( x )
#else
# define MYRRH_ISNAN( x ) ::boost::math::isnan( x )
# define MYRRH_ISFINITE( x ) ::boost::math::isfinite( x )
#endif

#ifdef max
#undef max
#endif //max

#ifdef min
#undef min
#endif //min

namespace myrrh {

template< typename T >
T int_power( T x, unsigned exponent )
{
    T result = T( 1.0 );
    for( unsigned i = 0; exponent != i; ++i ) result *= x;
    return result;
}

inline
int
nearest_integer( double a )
{
    return static_cast<int>( a < 0.0 ? a - 0.5 : a + 0.5 );
}

inline bool within_closed_range( double x, double lower, double higher )
{
    return MYRRH_ISFINITE( x ) && lower <= x && x <= higher;
}

inline bool is_probability( double x )
{
    return within_closed_range( x, 0.0, 1.0 );
}

template< typename T >
bool
is_close( T t1, T t2, T percent_tolerance )
{
    return boost::test_tools::check_is_close( t1, t2, MYRRH_FPC_NS::percent_tolerance( percent_tolerance ) );
}

inline
double
safe_exp( double x )
{
    return ( x <= GSL_LOG_DBL_MIN ) ? 0.0 : gsl_sf_exp( x );
}

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

/**
Returns the log of the sum of 2 values, when the values are passed as their logarithms.

I.e. we have log(p) and log(q), we want to calculate log(p+q).
*/
inline
double
log_sum_of_logs(
    double log_p,
    double log_q )
{
    //make sure log_p is the larger of the two, so that exp( log(q) - log(p) ) is small
    if( log_q > log_p ) std::swap( log_q, log_p );
    BOOST_ASSERT( log_p >= log_q );

    return log_p + gsl_sf_log_1plusx( safe_exp( log_q - log_p ) );
}

/** Ensures that the sum of the elements in the range adds to the given value. Returns the scale factor used. */
template< typename Range >
double
scale(
    Range & range,
    typename boost::range_value< Range >::type desired_sum = 1.0 )
{
    using boost::begin;
    using boost::end;
    using boost::range_value;

    typedef typename range_value< Range >::type value_type;

    const value_type current_sum =
        std::accumulate(
            begin( range ),
            end( range ),
            0.0 );

    const double scale_factor = desired_sum / current_sum;
    BOOST_FOREACH( value_type & v, range ) v *= scale_factor;

    return scale_factor;
}

namespace gsl {
namespace impl {

inline
void
error_handler(
    const char * reason,
    const char * file,
    int line,
    int gsl_errno)
{
    using namespace std;

    const string msg = MYRRH_MAKE_STRING( "GSL error(" << gsl_errno << "): " << file << "(" << line << "): " << reason );
    //cout << msg << endl;
    throw logic_error( msg );
}

} //namespace impl

inline
void init()
{
    static bool inited = false;
    if( ! inited )
    {
        gsl_set_error_handler( impl::error_handler );
        inited = true;
    }
}

inline
gsl_rng *
get_rng()
{
    static const gsl_rng_type * T = 0;
    static gsl_rng * r = 0;

    if( 0 == T || 0 == r )
    {
        init();

        gsl_rng_env_setup();

        T = gsl_rng_default;
        if( ! T ) throw std::logic_error( "Could not find default GSL random number generator type" );

        r = gsl_rng_alloc( T );
        if( ! r ) throw std::logic_error( "Could not allocate GSL random number generator" );
    }
    return r;
}

template< typename Range >
unsigned
one_sample_from_multinomial(
    const Range & dist,
    double weight = 1.0,
    gsl_rng * rng = gsl::get_rng() )
{
    double uniform_sample = gsl_rng_uniform( rng ) * weight;
    unsigned result = 0;
    BOOST_FOREACH( double p, dist )
    {
        uniform_sample -= p;
        if( uniform_sample < 0.0 ) return result;
        ++result;
    }
    throw std::logic_error( "Multinomial parameters did not sum to 1 (or the supplied weight)" );
}

template<
    typename ParameterRange,
    typename OutputRange
>
void
draw_from_dirichlet(
    const ParameterRange & alpha,
    OutputRange & output,
    gsl_rng * rng = gsl::get_rng() )
{
    using boost::size;
    using boost::begin;

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

    gsl_ran_dirichlet(
        rng,
        size( alpha ),
        &*boost::begin( alpha ),
        &*boost::begin( output ) );
}

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

    typedef typename range_value< ParameterRange >::type param_t;
    typedef typename range_value< ValueRange >::type value_t;
    typedef boost::tuple< param_t, value_t > tuple;

    BOOST_ASSERT( size( alpha ) == size( values ) );

    double result = 0.0;
    param_t alpha_sum = 0.0;

    BOOST_FOREACH(
        const tuple & t,
        make_iterator_range(
            make_zip_iterator(
                make_tuple(
                    const_begin( alpha ),
                    const_begin( values ) ) ),
            make_zip_iterator(
                make_tuple(
                    const_end( alpha ),
                    const_end( values )) ) ) )
    {
        const param_t a = t.get< 0 >();
        const value_t v = t.get< 1 >();
        alpha_sum += a;
        result += ( a - param_t( 1.0 ) ) * gsl_sf_log( v );
        std::cout << a << "\n";
        result -= gsl_sf_lngamma( a );
    }

    return result + gsl_sf_lngamma( alpha_sum );
}

/**
Pre-processed sampling from general discrete distributions.

See http://www.gnu.org/software/gsl/manual/html_node/General-Discrete-Distributions.html
*/
struct ran_discrete : boost::noncopyable
{
    typedef boost::shared_ptr< ran_discrete > ptr;

    gsl_ran_discrete_t * _pre_processed;

    inline ran_discrete(size_t K, const double * P)
        : _pre_processed( gsl_ran_discrete_preproc( K, P ) )
    {
    }

    template< typename Range >
    ran_discrete( const Range & P_range )
    {
        using namespace boost;
        const size_t K = size( P_range );
        std::vector< double > P( K );
        std::copy( begin( P_range ), end( P_range ), P.begin() );
        _pre_processed = gsl_ran_discrete_preproc( K, &(*(P.begin())) );
    }

    inline ~ran_discrete()
    {
        gsl_ran_discrete_free( _pre_processed );
    }

    size_t sample(gsl_rng * rng = get_rng()) const
    {
        return gsl_ran_discrete( rng, _pre_processed );
    }
};

} //namespace gsl
} //namespace myrrh

#endif //MYRRH_MATH_H_

max_stars_repo_name: JohnReid/myrrh

max_issues_repo_name: JohnReid/myrrh

max_forks_repo_name: JohnReid/myrrh


text:
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <float.h>
#include <math.h>
#include <time.h>
#include <sys/time.h>
#include <cblas.h>

/* Your function must have the following signature: */

void sgemm( int m, int n, float *A, float *C );

/* The benchmarking program */

int main( int argc, char **argv )
{
  srand(time(NULL));

  int n;
  double sum= 0;
  for(n = 32; n < 101; n+=54)

  /* Try different m */
  for( int m = 1000; m <10001; m+=511 )
  {
    /* Allocate and fill 2 random matrices A, C */
    float *A = (float*) malloc( m * n * sizeof(float) );
    float *C = (float*) malloc( m * m * sizeof(float) );

    for( int i = 0; i < m*n; i++ ) A[i] = 2 * drand48() - 1;
    for( int i = 0; i < m*m; i++ ) C[i] = 2 * drand48() - 1;

    /* measure Gflop/s rate; time a sufficiently long sequence of calls to eliminate noise */
    double Gflop_s, seconds = -1.0;
    for( int n_iterations = 1; seconds < 0.1; n_iterations *= 2 )
    {
      /* warm-up */
      sgemm( m, n, A, C );
      /* measure time */
      struct timeval start, end;
      gettimeofday( &start, NULL );
      for( int i = 0; i < n_iterations; i++ )
	sgemm( m,n, A, C );
      gettimeofday( &end, NULL );
      seconds = (end.tv_sec - start.tv_sec) + 1.0e-6 * (end.tv_usec - start.tv_usec);

      /* compute Gflop/s rate */
      Gflop_s = 2e-9 * n_iterations * m * m * n / seconds;
	 }
    sum += Gflop_s;

    printf( "%d by %d matrix \t %g Gflop/s\n", m, n, Gflop_s );
    //printf( "%g\n", Gflop_s );

    /* Ensure that error does not exceed the theoretical error bound */

    /* Set initial C to 0 and do matrix multiply of A*B */
    memset( C, 0, sizeof( float ) * m * m );
    sgemm( m,n, A, C );

    /* Subtract A*B from C using standard sgemm (note that this should be 0 to within machine roundoff) */
    cblas_sgemm( CblasColMajor,CblasNoTrans,CblasTrans, m,m,n, -1, A,m, A,m, 1, C,m );

    /* Subtract the maximum allowed roundoff from each element of C */
    for( int i = 0; i < m*n; i++ ) A[i] = fabs( A[i] );
    for( int i = 0; i < m*m; i++ ) C[i] = fabs( C[i] );
    cblas_sgemm( CblasColMajor,CblasNoTrans,CblasTrans, m,m,n, -3.0*FLT_EPSILON*n, A,m, A,m, 1, C,m );

    /* After this test if any element in C is still positive something went wrong in square_sgemm */
    for( int i = 0; i < m * m; i++ )
      if( C[i] > 0 ) {
	printf( "FAILURE: error in matrix multiply exceeds an acceptable margin\n" );
	return -1;
      }

    /* release memory */
    free( C );
    free( A );
  }

  return 0;
}

max_stars_repo_name: jasonaibrahim/AcademicWorks

max_issues_repo_name: jasonaibrahim/AcademicWorks

max_forks_repo_name: jasonaibrahim/AcademicWorks


text:
#pragma once

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

size_t nn_utils_sizeof_gsl_matrix(void);
size_t nn_utils_sizeof_gsl_vector(void);

max_stars_repo_name: mygulamali/neural-net

max_issues_repo_name: mygulamali/neural-net

max_forks_repo_name: mygulamali/neural-net


text:
#include <cerrno>
#include <cstdio>
#include <cstdlib>
#include <sys/stat.h>
#include <cmath>
#include <vector>
#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <time.h>

#include <gsl/gsl_odeiv2.h>
#include <gsl/gsl_math.h>
#include <math.h>       /* pow */
#include <gsl/gsl_sf.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_spline.h>

// functions particular to bh_clean.cpp

/* Time as a function of scale factor in seconds */
// integrand
double timeofa_int(double ap, void * params) {
  double integrand = 1./(ap*hubble(om,orad,ap));
  return integrand;
}

double timeofa(double a){
	gsl_integration_workspace * w  = gsl_integration_workspace_alloc (1000);
	double result, error;
	gsl_function F;
	F.function = &timeofa_int;
	gsl_integration_qags (&F, ai, a, 0, 1e-7, 1000, w, &result, &error);
  gsl_integration_workspace_free (w);
	return result;
}

/* scale factor as a function of time */
// fill array to then spline
void aoftime(arrays_T xxyy){
		for(int i=0; i<1001; i++){
			double a = ai*exp(i*log(0.9999/ai)/(1000.));
			(*xxyy).xx[i] = timeofa(a);
			(*xxyy).yy[i] = a;
		}
		// set final time
		(*xxyy).xx[1001] = 1./hubble(om,orad,1.);
		(*xxyy).yy[1001] = 1.;
}

/* Log10[Mass] as a function of scale factor  --  not a step function*/
// M = 10^lgm ,   a is the scale factor , rem= false for no planck mass remnants
double bhmasslg(double lgm, double a, bool rem){
	// If M_i^3 <= 3 K M_p^4 t  then the black hole should have  0 mass ---> Log10[M_i] <= 1/3 (Log10[ 3 K M_p^4] + Log10[t]) = 1/3 ( 17.803 + Log10[t])
	double logdect = 1./3. * (17.803 + log10(timeofa(a)));
	if (lgm>logdect) {
		return log10(pow(pow(10., lgm*3.) - decfac*timeofa(a), 1.0/3.));
	}
	else{
		if (!rem) {
			return -100.;
		}
		else{
			return lgmp;
		}
	}
}

// Initial mass from final mass. Any bh that has decayed completely earlier than a will be mapped to decfac*timeofa(a)^(1/3)
// (There's no way of knowing with certainty what its initial mass is if its final mass is 0)
double bhmasslgi(double lgm, double a, bool rem){
  double logdect = 1./3. * (17.803 + log10(timeofa(a)));
  double mymasslg = log10(pow(10., lgm*3.) + decfac*timeofa(a))/3.;
  if(lgm<logdect ){
    return -100;
  }
  else if(mymasslg>lgmax){
		return -100;
    }
  else{
    return mymasslg;
  }
}

/* Calculate normalisation for psibroad */
// normalisation integrand --- see equation 3.4 for example
// Params are peak mass (as proxy for M_s)
// Normalisation = 1/ Integrate[ 10^psilg dM ] =   1/ Integrate[ 10^(psilg + lgm) x Ln[10] dlgm ]

double normlg_int(double lgm, void * params) {
	myparam_type pars = *(myparam_type *)(params);
	double peakm = pars.peakm;
	double expt = lgm + psibroadlg(lgm,peakm);
	double integrand =pow(10.,expt);
  return integrand;
}

// Normalisation = 10^normlg
double normlg(double peakm){
	struct myparam_type pars = {peakm, 1., ai, false};

	gsl_integration_workspace * w  = gsl_integration_workspace_alloc (1000);
	double result, error;
	gsl_function F;
	F.function = &normlg_int;
	F.params = &pars;
	gsl_integration_qags (&F, lgmp, lgmax, 0., 1e-5, 1000, w, &result, &error);
  gsl_integration_workspace_free (w);
	return  log10(1./(log(10.)*result));
}

/* Calculate normalisation for psi_mono */
// normalisation integrand --- see equation 3.4 for example
// Params are peak mass
// Normalisation = 1/ Integrate[ 10^psilg dM ] =   1/ Integrate[ 10^(psilg + lgm) x Ln[10] dlgm ]

double normlg_int_mono(double lgm, void * params) {
	myparam_type pars = *(myparam_type *)(params);
	double peakm = pars.peakm;
	double expt = lgm + psilowlg(lgm,peakm);
	double integrand =pow(10.,expt);
  return integrand;
}

// Normalisation = 10^normlg
double normlg_mono(double peakm){
	struct myparam_type pars = {peakm, 1., ai, false};

	gsl_integration_workspace * w  = gsl_integration_workspace_alloc (1000);
	double result, error;
	gsl_function F;
	F.function = &normlg_int_mono;
	F.params = &pars;
	gsl_integration_qags (&F, lgmp, lgmax, 0., 1e-5, 1000, w, &result, &error);
  gsl_integration_workspace_free (w);
	return  log10(1./(log(10.)*result));
}

max_stars_repo_name: nebblu/PBH

max_issues_repo_name: nebblu/PBH

max_forks_repo_name: nebblu/PBH


text:
/* The MIT License

   Copyright (c) 2013-2015 Genome Research Ltd.

   Author: Petr Danecek <pd3@sanger.ac.uk>

   Permission is hereby granted, free of charge, to any person obtaining a copy
   of this software and associated documentation files (the "Software"), to deal
   in the Software without restriction, including without limitation the rights
   to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
   copies of the Software, and to permit persons to whom the Software is
   furnished to do so, subject to the following conditions:

   The above copyright notice and this permission notice shall be included in
   all copies or substantial portions of the Software.

   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
   IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
   AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
   LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
   OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
   THE SOFTWARE.

 */

#include <stdio.h>
#include <stdlib.h>
#include <getopt.h>
#include <sys/stat.h>
#include <sys/types.h>
#include <unistd.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multifit_nlin.h>
#include <htslib/vcf.h>
#include <htslib/synced_bcf_reader.h>
#include "bcftools.h"
#include "peakfit.h"

typedef struct
{
    int nvals;          // all values, including RR,AA peaks
    double *xvals;      // pointer to args_t.xvals
    double *yvals;
    int copy_number;    // heuristics to skip futile CN1 fits when no het peak is detected
    int irr, ira, iaa;  // chop off RR and AA peaks
    char *chr;
}
dist_t;

typedef struct
{
    int ndist, nbins, ra_rr_scaling, smooth;
    double *xvals;
    dist_t *dist;
    char **argv, *output_dir;
    double fit_th, peak_symmetry, cn_penalty, min_peak_size, min_fraction;
    int argc, plot, verbose, regions_is_file, targets_is_file, include_aa, force_cn;
    char *dat_fname, *fname, *regions_list, *targets_list, *sample;
    FILE *dat_fp;
}
args_t;

FILE *open_file(char **fname, const char *mode, const char *fmt, ...);

static void init_dist(args_t *args, dist_t *dist, int verbose)
{
    // isolate RR and AA peaks and rescale so that they are comparable to hets
    int i, irr, iaa, n = dist->nvals;

    // smooth the distribution, this is just to find the peaks
    double *tmp = (double*) malloc(sizeof(double)*n);
    int win  = args->smooth ? fabs(args->smooth)*2 + 1 : 7;   // must be an odd number
    int hwin = win/2;
    double avg = tmp[0] = dist->yvals[0];
    for (i=1; i<hwin; i++)
    {
        avg += dist->yvals[2*i-1];
        tmp[i] = avg/(2*i+1);
    }
    avg = 0;
    for (i=0; i<n; i++)
    {
        avg += dist->yvals[i];
        if ( i>=win-1 )
        {
            tmp[i-hwin] = avg/win;
            avg -= dist->yvals[i-win+1];
        }
    }
    for (i=n-hwin; i<n; i++)
    {
        avg -= dist->yvals[i-hwin];
        hwin--;
        tmp[i] = avg/(2*hwin+1);
        avg -= dist->yvals[i-hwin];
    }

    // find the extremes; first a simple approach: find a gap
    for (irr=0,i=0; i<n/2; i++) if ( tmp[i] < tmp[irr] ) irr = i;
    for (iaa=n-1,i=n-1; i>=n/2; i--) if ( tmp[i] < tmp[iaa] ) iaa = i;
    irr += win*0.5;
    iaa += win*0.5;
    if ( iaa>=n ) iaa = n-1;
    if ( irr>=iaa ) error("FIXME: oops, dist normalization failed for %s: %d vs %d\n", dist->chr,irr,iaa); // we may need to be smarter
    if ( args->smooth>0 ) for (i=0; i<n; i++) dist->yvals[i] = tmp[i];
    free(tmp);

    // clean the data: the AA peak is occasionally not centered at 1.0 but is closer to the center, chop off
    int imax_aa = iaa;
    for (i=iaa; i<n; i++)
        if ( dist->yvals[imax_aa] < dist->yvals[i] ) imax_aa = i;
    dist->nvals = imax_aa+1;
    if ( iaa>=dist->nvals ) iaa = dist->nvals-1;

    // find the maximum and scale the peaks (first draft: no attempt to join the segments smootly)
    double max_rr = 0, max_aa = 0, max_ra = 0, srr = 0, saa = 0, sra = 0;
    for (i=0; i<irr; i++)
    {
        srr += dist->yvals[i];
        if ( max_rr < dist->yvals[i] ) max_rr = dist->yvals[i];
    }
    for (i=irr; i<=iaa; i++)
    {
        sra += dist->yvals[i];
        if ( max_ra < dist->yvals[i] ) max_ra = dist->yvals[i];
    }
    for (i=iaa+1; i<n; i++)
    {
        saa += dist->yvals[i];
        if ( max_aa < dist->yvals[i] ) max_aa = dist->yvals[i];
    }

    // Does the het peak exist at all? Usually the numbers are as follows:
    //  1:  cn=0     ra/rr=0.205730      aa/ra=0.674922      nra=7066
    //  20: cn=0     ra/rr=0.258019      aa/ra=0.548929      nra=2381
    //  X:  cn=0     ra/rr=0.005976      aa/ra=44.116667     nra=60
    //  Y:  cn=0     ra/rr=0.008316      aa/ra=7.250000      nra=12
    //  MT: cn=0     ra/rr=0.013699      aa/ra=0.666667      nra=3

    if ( !args->ra_rr_scaling ) max_ra = max_aa = max_rr;
    if ( !sra || (sra/srr<0.1 && saa/sra>1.0) ) // too few hets, CN1
    {
        max_ra = max_aa;
        dist->copy_number = 1;
    }
    else if ( sra/srr<0.1 || saa/sra>1.0 )
    {
        max_ra = max_aa;
        dist->copy_number = -1;     // unknown copy number
    }
    if ( max_rr ) for (i=0; i<irr; i++) dist->yvals[i] /= max_rr;
    if ( max_ra ) for (i=irr; i<=iaa; i++) dist->yvals[i] /= max_ra;
    if ( max_aa ) for (i=iaa+1; i<n; i++) dist->yvals[i] /= max_aa;

    dist->irr = irr;
    dist->iaa = iaa;
    dist->ira = n*0.5;

    if ( verbose )
        fprintf(stderr,"%s:\t irr,ira,iaa=%.2f,%.2f,%.2f \t cn=%2d \t ra/rr=%f \t aa/ra=%f \t nra=%d\n",
            dist->chr, dist->xvals[irr],dist->xvals[dist->ira],dist->xvals[iaa],
            dist->copy_number,sra/srr,saa/sra, (int)sra);
}

static void init_data(args_t *args)
{
    bcf_srs_t *files = bcf_sr_init();
    if ( args->regions_list )
    {
        if ( bcf_sr_set_regions(files, args->regions_list, args->regions_is_file)<0 )
            error("Failed to read the regions: %s\n", args->regions_list);
    }
    if ( args->targets_list )
    {
        if ( bcf_sr_set_targets(files, args->targets_list, args->targets_is_file, 0)<0 )
            error("Failed to read the targets: %s\n", args->targets_list);
    }
    if ( !bcf_sr_add_reader(files, args->fname) ) error("Failed to open %s: %s\n", args->fname,bcf_sr_strerror(files->errnum));
    bcf_hdr_t *hdr = files->readers[0].header;
    if ( !args->sample )
    {
        if ( bcf_hdr_nsamples(hdr)>1 ) error("Missing the option -s, --sample\n");
        args->sample = hdr->samples[0];
    }
    else if ( bcf_hdr_id2int(hdr,BCF_DT_SAMPLE,args->sample)<0 ) error("No such sample: %s\n", args->sample);
    int ret = bcf_hdr_set_samples(hdr, args->sample, 0);
    if ( ret<0 ) error("Error setting the sample: %s\n", args->sample);

    if ( !bcf_hdr_idinfo_exists(hdr,BCF_HL_FMT,bcf_hdr_id2int(hdr,BCF_DT_ID,"BAF")) )
        error("The tag FORMAT/BAF is not present in the VCF: %s\n", args->fname);

    int i;
    args->xvals = (double*) calloc(args->nbins,sizeof(double));
    for (i=0; i<args->nbins; i++) args->xvals[i] = 1.0*i/(args->nbins-1);

    // collect BAF distributions for all chromosomes
    int idist = -1, nbaf = 0, nprocessed = 0, ntotal = 0, prev_chr = -1;
    float *baf = NULL;
    while ( bcf_sr_next_line(files) )
    {
        ntotal++;

        bcf1_t *line = bcf_sr_get_line(files,0);
        if ( bcf_get_format_float(hdr,line,"BAF",&baf,&nbaf) != 1 ) continue;
        if ( bcf_float_is_missing(baf[0]) ) continue;

        nprocessed++;

        if ( prev_chr==-1 || prev_chr!=line->rid )
        {
            // new chromosome
            idist = args->ndist++;
            args->dist = (dist_t*) realloc(args->dist, sizeof(dist_t)*args->ndist);
            memset(&args->dist[idist],0,sizeof(dist_t));
            args->dist[idist].chr   = strdup(bcf_seqname(hdr,line));
            args->dist[idist].yvals = (double*) calloc(args->nbins,sizeof(double));
            args->dist[idist].xvals = args->xvals;
            args->dist[idist].nvals = args->nbins;
            prev_chr = line->rid;
        }
        int bin = baf[0]*(args->nbins-1);
        args->dist[idist].yvals[bin]++;   // the distribution
    }
    free(baf);
    bcf_sr_destroy(files);

    for (idist=0; idist<args->ndist; idist++)
    {
        #if 0
            int j;
            for (j=0; j<args->nbins; j++)
            {
                double x = args->dist[idist].xvals[j];
                args->dist[idist].yvals[j] = exp(-(x-0.5)*(x-0.5)/1e-3);
            }
        #endif
        init_dist(args, &args->dist[idist],args->verbose);
    }

    args->dat_fp = open_file(&args->dat_fname,"w","%s/dist.dat", args->output_dir);
    fprintf(args->dat_fp, "# This file was produced by: bcftools polysomy(%s+htslib-%s), the command line was:\n", bcftools_version(),hts_version());
    fprintf(args->dat_fp, "# \t bcftools %s ", args->argv[0]);
    for (i=1; i<args->argc; i++)
        fprintf(args->dat_fp, " %s",args->argv[i]);
    fprintf(args->dat_fp,"\n#\n");
    fprintf(args->dat_fp,"# DIST\t[2]Chrom\t[3]BAF\t[4]Normalized Count\n");
    fprintf(args->dat_fp,"# FIT\t[2]Goodness of Fit\t[3]iFrom\t[4]iTo\t[5]The Fitted Function\n");
    fprintf(args->dat_fp,"# CN\t[2]Chrom\t[3]Estimated Copy Number\t[4]Absolute fit deviation\n");

    char *fname = NULL;
    FILE *fp = open_file(&fname,"w","%s/dist.py", args->output_dir);
//-------- matplotlib script --------------
    fprintf(fp,
        "#!/usr/bin/env python\n"
        "#\n"
        "import matplotlib as mpl\n"
        "mpl.use('Agg')\n"
        "import matplotlib.pyplot as plt\n"
        "import csv,sys,argparse\n"
        "from math import exp\n"
        "\n"
        "outdir = '%s'\n"
        "\n"
        "def read_dat(dat,fit,cn):\n"
        "   csv.register_dialect('tab', delimiter='\t', quoting=csv.QUOTE_NONE)\n"
        "   with open(outdir+'/dist.dat', 'rb') as f:\n"
        "      reader = csv.reader(f, 'tab')\n"
        "      for row in reader:\n"
        "          if row[0][0]=='#': continue\n"
        "          type = row[0]\n"
        "          chr  = row[1]\n"
        "          if type=='DIST':\n"
        "              if chr not in dat: dat[chr] = []\n"
        "              dat[chr].append(row)\n"
        "          elif type=='FIT':\n"
        "              if chr not in fit: fit[chr] = []\n"
        "              fit[chr].append(row)\n"
        "          elif type=='CN':\n"
        "              cn[chr] = row[2]\n"
        "\n"
        "def plot_dist(dat,fit,chr):\n"
        "   fig, ax = plt.subplots(1, 1, figsize=(7,5))\n"
        "   ax.plot([x[2] for x in dat[chr]],[x[3] for x in dat[chr]],'k-',label='Distribution')\n"
        "   if chr in fit:\n"
        "       for i in range(len(fit[chr])):\n"
        "           pfit = fit[chr][i]\n"
        "           exec('def xfit(x): return '+pfit[5])\n"
        "           istart = int(pfit[3])\n"
        "           iend   = int(pfit[4])+1\n"
        "           vals   = dat[chr][istart:iend]\n"
        "           args   = {}\n"
        "           if i==0: args = {'label':'Target to Fit'}\n"
        "           ax.plot([x[2] for x in vals],[x[3] for x in vals],'r-',**args)\n"
        "           if i==0: args = {'label':'Best Fit'}\n"
        "           ax.plot([x[2] for x in vals],[xfit(float(x[2])) for x in vals],'g-',**args)\n"
        "   ax.set_title('BAF distribution, chr'+chr)\n"
        "   ax.set_xlabel('BAF')\n"
        "   ax.set_ylabel('Frequency')\n"
        "   ax.legend(loc='best',prop={'size':7},frameon=False)\n"
        "   plt.savefig(outdir+'/dist.chr'+chr+'.png')\n"
        "   plt.close()\n"
        "\n"
        "def plot_copy_number(cn):\n"
        "   fig, ax = plt.subplots(1, 1, figsize=(7,5))\n"
        "   xlabels = sorted(cn.keys())\n"
        "   xvals = range(len(xlabels))\n"
        "   yvals = [float(cn[x]) for x in xlabels]\n"
        "   ax.plot(xvals,yvals,'o',color='red')\n"
        "   for i in range(len(xvals)):\n"
        "       if yvals[i]==-1: ax.annotate('?', xy=(xvals[i],0.5),va='center',ha='center',color='red',fontweight='bold')\n"
        "   ax.tick_params(axis='both', which='major', labelsize=9)\n"
        "   ax.set_xticks(xvals)\n"
        "   ax.set_xticklabels(xlabels,rotation=45)\n"
        "   ax.set_xlim(-1,len(xlabels))\n"
        "   ax.set_ylim(0,5.0)\n"
        "   ax.set_yticks([1.0,2.0,3.0,4.0])\n"
        "   ax.set_xlabel('Chromosome')\n"
        "   ax.set_ylabel('Copy Number')\n"
        "   plt.savefig(outdir+'/copy-number.png')\n"
        "   plt.close()\n"
        "\n"
        "class myParser(argparse.ArgumentParser):\n"
        "   def error(self, message):\n"
        "       self.print_help()\n"
        "       sys.stderr.write('error: %%s\\n' %% message)\n"
        "       sys.exit(2)\n"
        "\n"
        "def main():\n"
        "   parser = myParser()\n"
        "   parser.add_argument('-a', '--all', action='store_true', help='Create all plots')\n"
        "   parser.add_argument('-c', '--copy-number', action='store_true', help='Create copy-number plot')\n"
        "   parser.add_argument('-d', '--distrib', metavar='CHR', help='Plot BAF distribution of a single chromosome')\n"
        "   args = parser.parse_args()\n"
        "   dat = {}; fit = {}; cn = {}\n"
        "   read_dat(dat,fit,cn)\n"
        "   if args.distrib!=None:\n"
        "       plot_dist(dat,fit,args.distrib)\n"
        "   if args.all:\n"
        "       for chr in dat: plot_dist(dat,fit,chr)\n"
        "       plot_copy_number(cn)\n"
        "   elif args.copy_number:\n"
        "       plot_copy_number(cn)\n"
        "   else:\n"
        "       for chr in dat: plot_dist(dat,fit,chr)\n"
        "\n"
        "if __name__ == '__main__':\n"
        "   main()\n",
        args->output_dir);
//---------------------------------------
    chmod(fname, S_IWUSR|S_IRUSR|S_IRGRP|S_IROTH|S_IXUSR|S_IXGRP|S_IXOTH);
    free(fname);
    fclose(fp);
}

static void destroy_data(args_t *args)
{
    int i;
    for (i=0; i<args->ndist; i++)
    {
        free(args->dist[i].chr);
        free(args->dist[i].yvals);
    }
    free(args->dist);
    free(args->xvals);
    free(args->dat_fname);
    fclose(args->dat_fp);
}

static void save_dist(args_t *args, dist_t *dist)
{
    int i;
    for (i=0; i<args->nbins; i++)
        fprintf(args->dat_fp,"DIST\t%s\t%f\t%f\n",dist->chr,dist->xvals[i],dist->yvals[i]);
}
static void fit_curves(args_t *args)
{
    peakfit_t *pkf = peakfit_init();
    peakfit_verbose(pkf,args->verbose);

    int i, nmc = 50;
    for (i=0; i<args->ndist; i++)
    {
        dist_t *dist = &args->dist[i];
        save_dist(args, &args->dist[i]);

        if ( dist->copy_number!=0 )
        {
            fprintf(args->dat_fp,"CN\t%s\t%.2f\n", dist->chr,(float)dist->copy_number);
            continue;
        }

        if ( args->verbose )
            fprintf(stderr,"%s:\n", dist->chr);

        int nrr_aa  = dist->iaa - dist->irr + 1;
        int nrr_ra  = dist->ira - dist->irr + 1;
        int naa_max = dist->nvals - dist->iaa;
        double xrr  = dist->xvals[dist->irr], *xrr_vals = &dist->xvals[dist->irr], *yrr_vals = &dist->yvals[dist->irr];
        double xaa  = dist->xvals[dist->iaa], *xaa_vals = &dist->xvals[dist->iaa], *yaa_vals = &dist->yvals[dist->iaa];
        double xra  = dist->xvals[dist->ira];
        double xmax = dist->xvals[dist->nvals-1];

        // CN2
        double cn2aa_fit = 0, cn2ra_fit, cn2_fit;
        char *cn2aa_func = 0, *cn2ra_func;
        double cn2aa_params[3] = {1,1,1} ,cn2ra_params[3];
        if ( args->include_aa )
        {
            peakfit_reset(pkf);
            peakfit_add_exp(pkf, 1.0,1.0,0.2, 5);
            peakfit_set_mc(pkf, 0.01,0.3,2,nmc);
            peakfit_set_mc(pkf, 0.05,1.0,0,nmc);
            cn2aa_fit  = peakfit_run(pkf, naa_max, xaa_vals, yaa_vals);
            cn2aa_func = strdup(peakfit_sprint_func(pkf));
            peakfit_get_params(pkf,0,cn2aa_params,3);
        }
        peakfit_reset(pkf);
        peakfit_add_bounded_gaussian(pkf, 1.0,0.5,0.03, 0.45,0.55, 7);
        peakfit_set_mc(pkf, 0.01,0.3,2,nmc);
        peakfit_set_mc(pkf, 0.05,1.0,0,nmc);
        cn2ra_fit  = peakfit_run(pkf, nrr_aa,xrr_vals,yrr_vals);
        cn2ra_func = strdup(peakfit_sprint_func(pkf));
        cn2_fit    = cn2ra_fit + cn2aa_fit;
        peakfit_get_params(pkf,0,cn2ra_params,3);

        // CN3: fit two peaks, then enforce the symmetry and fit again
        double cn3rra_params[5], cn3raa_params[5], *cn3aa_params = cn2aa_params;
        double cn3aa_fit = cn2aa_fit, cn3ra_fit;
        char *cn3aa_func = cn2aa_func, *cn3ra_func;
        double min_dx3   = 0.5 - 1./(args->min_fraction+2);
        peakfit_reset(pkf);
        peakfit_add_bounded_gaussian(pkf, 1.0,1/3.,0.03, xrr,xra-min_dx3, 7);
        peakfit_set_mc(pkf, xrr,xra-min_dx3, 1,nmc);
        peakfit_add_bounded_gaussian(pkf, 1.0,2/3.,0.03, xra+min_dx3,xaa, 7);
        peakfit_set_mc(pkf, xra+min_dx3,xaa, 1,nmc);
        peakfit_run(pkf, nrr_aa, xrr_vals, yrr_vals);
        // force symmetry around x=0.5
        peakfit_get_params(pkf,0,cn3rra_params,5);
        peakfit_get_params(pkf,1,cn3raa_params,5);
        double cn3_dx = (0.5-cn3rra_params[1] + cn3raa_params[1]-0.5)*0.5;
        if ( cn3_dx > 0.5/3 ) cn3_dx = 0.5/3;   // CN3 peaks should not be separated by more than 1/3
        peakfit_reset(pkf);
        peakfit_add_gaussian(pkf, cn3rra_params[0],0.5-cn3_dx,cn3rra_params[2], 5);
        peakfit_add_gaussian(pkf, cn3raa_params[0],0.5+cn3_dx,cn3raa_params[2], 5);
        cn3ra_fit  = peakfit_run(pkf, nrr_aa, xrr_vals, yrr_vals);
        cn3ra_func = strdup(peakfit_sprint_func(pkf));
        // compare peak sizes
        peakfit_get_params(pkf,0,cn3rra_params,3);
        peakfit_get_params(pkf,1,cn3raa_params,3);
        double cn3rra_size = cn3rra_params[0]*cn3rra_params[0];
        double cn3raa_size = cn3raa_params[0]*cn3raa_params[0];
        double cn3_dy      = cn3rra_size > cn3raa_size ? cn3raa_size/cn3rra_size : cn3rra_size/cn3raa_size;
        double cn3_frac    = (1 - 2*cn3rra_params[1]) / cn3rra_params[1];
        double cn3_fit     = cn3ra_fit + cn3aa_fit;
        // A very reasonable heuristics: check if the peak's width converged, exclude far too broad or far too narrow peaks
        if ( cn3rra_params[2]>0.3  || cn3raa_params[2]>0.3 ) cn3_fit = HUGE_VAL;
        if ( cn3rra_params[2]<1e-2 || cn3raa_params[2]<1e-2 ) cn3_fit = HUGE_VAL;

        // CN4 (contaminations)
        // - first fit only the [0,0.5] part of the data, then enforce the symmetry and fit again
        // - min_frac=1 (resp. 0.5) is interpreted as 50:50% (rep. 75:25%) contamination
        double cn4AAaa_params[3] = {1,1,1} ,cn4AAra_params[3] = {1,1,1}, cn4RAra_params[3], cn4RArr_params[5], cn4RAaa_params[5];
        double cn4aa_fit = 0, cn4ra_fit;
        char *cn4aa_func = 0, *cn4ra_func;
        double min_dx4   = 0.25*args->min_fraction;
        if ( args->include_aa )
        {
            peakfit_reset(pkf);
            peakfit_add_exp(pkf, 0.5,1.0,0.2, 5);
            peakfit_set_mc(pkf, 0.01,0.3,2,nmc);
            peakfit_add_bounded_gaussian(pkf, 0.4,(xaa+xmax)*0.5,2e-2, xaa,xmax, 7);
            peakfit_set_mc(pkf, xaa,xmax, 1,nmc);
            cn4aa_fit  = peakfit_run(pkf, naa_max, xaa_vals,yaa_vals);
            cn4aa_func = strdup(peakfit_sprint_func(pkf));
            peakfit_get_params(pkf,0,cn4AAaa_params,3);
            peakfit_get_params(pkf,1,cn4AAra_params,5);
        }
        peakfit_reset(pkf);
        // first fit only the [0,0.5] part of the data
        peakfit_add_gaussian(pkf, 1.0,0.5,0.03, 5);
        peakfit_add_bounded_gaussian(pkf, 0.6,0.3,0.03, xrr,xra-min_dx4, 7);
        peakfit_set_mc(pkf, xrr,xra-min_dx4,2,nmc);
        peakfit_run(pkf, nrr_ra , xrr_vals, yrr_vals);
        // now forcet symmetry around x=0.5
        peakfit_get_params(pkf,0,cn4RAra_params,3);
        peakfit_get_params(pkf,1,cn4RArr_params,5);
        double cn4_dx = 0.5-cn4RArr_params[1];
        if ( cn4_dx > 0.25 ) cn4_dx = 0.25;   // CN4 peaks should not be separated by more than 0.5
        peakfit_reset(pkf);
        peakfit_add_gaussian(pkf, cn4RAra_params[0],0.5,cn4RAra_params[2], 5);
        peakfit_add_gaussian(pkf, cn4RArr_params[0],0.5-cn4_dx,cn4RArr_params[2], 5);
        peakfit_add_gaussian(pkf, cn4RArr_params[0],0.5+cn4_dx,cn4RArr_params[2], 5);
        peakfit_set_mc(pkf, 0.1,cn4RAra_params[0],0,nmc);
        peakfit_set_mc(pkf, 0.01,0.1,2,nmc);
        cn4ra_fit  = peakfit_run(pkf, nrr_aa , xrr_vals, yrr_vals);
        cn4ra_func = strdup(peakfit_sprint_func(pkf));
        peakfit_get_params(pkf,0,cn4RAra_params,3);
        peakfit_get_params(pkf,1,cn4RArr_params,3);
        peakfit_get_params(pkf,2,cn4RAaa_params,3);
        double cn4RAra_size = cn4RAra_params[0]==0 ? HUGE_VAL : cn4RAra_params[0]*cn4RAra_params[0];
        double cn4RArr_size = cn4RArr_params[0]*cn4RArr_params[0];
        double cn4RAaa_size = cn4RAaa_params[0]*cn4RAaa_params[0];
        double cn4RArr_dy   = cn4RArr_size < cn4RAra_size ? cn4RArr_size/cn4RAra_size : cn4RAra_size/cn4RArr_size;
        double cn4RAaa_dy   = cn4RAaa_size < cn4RAra_size ? cn4RAaa_size/cn4RAra_size : cn4RAra_size/cn4RAaa_size;
        double cn4_dy       = cn4RArr_dy < cn4RAaa_dy ? cn4RArr_dy/cn4RAaa_dy : cn4RAaa_dy/cn4RArr_dy;
        double cn4_ymin     = cn4RArr_size < cn4RAaa_size ? cn4RArr_size/cn4RAra_size : cn4RAaa_size/cn4RAra_size;
        cn4_dx              = (cn4RAaa_params[1]-0.5) - (0.5-cn4RArr_params[1]);
        double cn4_frac     = cn4RAaa_params[1] - cn4RArr_params[1];
        double cn4_fit      = cn4ra_fit + cn4aa_fit;
        // A very reasonable heuristics: check if the peak's width converged, exclude far too broad or far too narrow peaks
        if ( cn4RAra_params[2]>0.3 || cn4RArr_params[2]>0.3 || cn4RAaa_params[2]>0.3 ) cn4_fit = HUGE_VAL;
        if ( cn4RAra_params[2]<1e-2 || cn4RArr_params[2]<1e-2 || cn4RAaa_params[2]<1e-2 ) cn4_fit = HUGE_VAL;

        // Choose the best match
        char cn2_fail = '*', cn3_fail = '*', cn4_fail = '*';
        if ( cn2_fit > args->fit_th ) cn2_fail = 'f';

        if ( cn3_fit > args->fit_th ) cn3_fail = 'f';
        else if ( cn3_dy < args->peak_symmetry ) cn3_fail = 'y';    // size difference is too big

        if ( cn4_fit > args->fit_th ) cn4_fail = 'f';
        else if ( cn4_ymin < args->min_peak_size ) cn4_fail = 'y';      // side peak is too small
        else if ( cn4_dy < args->peak_symmetry ) cn4_fail = 'Y';    // size difference is too big
        else if ( cn4_dx > 0.1 ) cn4_fail = 'x';                    // side peaks placed assymetrically

        double cn = -1, fit = cn2_fit;
        if ( cn2_fail == '*' ) { cn = 2; fit = cn2_fit; }
        if ( cn3_fail == '*' )
        {
            // Use cn_penalty as a tiebreaker. If set to 0.3, cn3_fit must be 30% smaller than cn2_fit.
            if ( cn<0 || cn3_fit < (1-args->cn_penalty) * fit )
            {
                cn = 2 + cn3_frac;
                fit = cn3_fit;
                if ( cn2_fail=='*' ) cn2_fail = 'p';
            }
            else cn3_fail = 'p';
        }
        if ( cn4_fail == '*' )
        {
            if ( cn<0 || cn4_fit < (1-args->cn_penalty) * fit )
            {
                cn = 3 + cn4_frac;
                fit = cn4_fit;
                if ( cn2_fail=='*' ) cn2_fail = 'p';
                if ( cn3_fail=='*' ) cn3_fail = 'p';
            }
            else cn4_fail = 'p';
        }

        if ( args->verbose )
        {
            fprintf(stderr,"\tcn2 %c fit=%e\n", cn2_fail, cn2_fit);
            fprintf(stderr,"\t       .. %e\n", cn2ra_fit);
            fprintf(stderr,"\t            RA:   %f %f %f\n", cn2ra_params[0],cn2ra_params[1],cn2ra_params[2]);
            fprintf(stderr,"\t       .. %e\n", cn2aa_fit);
            fprintf(stderr,"\t            AA:   %f %f %f\n", cn2aa_params[0],cn2aa_params[1],cn2aa_params[2]);
            fprintf(stderr,"\t      func:\n");
            fprintf(stderr,"\t            %s\n", cn2ra_func);
            fprintf(stderr,"\t            %s\n", cn2aa_func);
            fprintf(stderr,"\n");
            fprintf(stderr,"\tcn3 %c fit=%e  frac=%f  symmetry=%f\n", cn3_fail, cn3_fit, cn3_frac, cn3_dy);
            fprintf(stderr,"\t       .. %e\n", cn3ra_fit);
            fprintf(stderr,"\t            RRA:  %f %f %f\n", cn3rra_params[0],cn3rra_params[1],cn3rra_params[2]);
            fprintf(stderr,"\t            RAA:  %f %f %f\n", cn3raa_params[0],cn3raa_params[1],cn3raa_params[2]);
            fprintf(stderr,"\t       .. %e\n", cn3aa_fit);
            fprintf(stderr,"\t            AAA:  %f %f %f\n", cn3aa_params[0],cn3aa_params[1],cn3aa_params[2]);
            fprintf(stderr,"\t      func:\n");
            fprintf(stderr,"\t            %s\n", cn3ra_func);
            fprintf(stderr,"\t            %s\n", cn3aa_func);
            fprintf(stderr,"\n");
            fprintf(stderr,"\tcn4 %c fit=%e  frac=%f  symmetry=%f ymin=%f\n", cn4_fail, cn4_fit, cn4_frac, cn4_dy, cn4_ymin);
            fprintf(stderr,"\t       .. %e\n", cn4ra_fit);
            fprintf(stderr,"\t            RArr:  %f %f %f\n", cn4RArr_params[0],cn4RArr_params[1],cn4RArr_params[2]);
            fprintf(stderr,"\t            RAra:  %f %f %f\n", cn4RAra_params[0],cn4RAra_params[1],cn4RAra_params[2]);
            fprintf(stderr,"\t            RAaa:  %f %f %f\n", cn4RAaa_params[0],cn4RAaa_params[1],cn4RAaa_params[2]);
            fprintf(stderr,"\t       .. %e\n", cn4aa_fit);
            fprintf(stderr,"\t            AAaa:  %f %f %f\n", cn4AAaa_params[0],cn4AAaa_params[1],cn4AAaa_params[2]);
            fprintf(stderr,"\t      func:\n");
            fprintf(stderr,"\t            %s\n", cn4ra_func);
            fprintf(stderr,"\t            %s\n", cn4aa_func);
            fprintf(stderr,"\n");
        }

        if ( args->force_cn==2 || cn2_fail == '*' )
        {
            fprintf(args->dat_fp,"FIT\t%s\t%e\t%d\t%d\t%s\n", dist->chr,cn2ra_fit,dist->irr,dist->iaa,cn2ra_func);
            if ( cn2aa_func ) fprintf(args->dat_fp,"FIT\t%s\t%e\t%d\t%d\t%s\n", dist->chr,cn2aa_fit,dist->iaa,dist->nvals-1,cn2aa_func);
        }
        if ( args->force_cn==3 || cn3_fail == '*' )
        {
            fprintf(args->dat_fp,"FIT\t%s\t%e\t%d\t%d\t%s\n", dist->chr,cn3ra_fit,dist->irr,dist->iaa,cn3ra_func);
            if ( cn3aa_func ) fprintf(args->dat_fp,"FIT\t%s\t%e\t%d\t%d\t%s\n", dist->chr,cn3aa_fit,dist->iaa,dist->nvals-1,cn3aa_func);
        }
        if ( args->force_cn==4 || cn4_fail == '*' )
        {
            fprintf(args->dat_fp,"FIT\t%s\t%e\t%d\t%d\t%s\n", dist->chr,cn4ra_fit,dist->irr,dist->iaa,cn4ra_func);
            if ( cn4aa_func ) fprintf(args->dat_fp,"FIT\t%s\t%e\t%d\t%d\t%s\n", dist->chr,cn4aa_fit,dist->iaa,dist->nvals-1,cn4aa_func);
        }
        fprintf(args->dat_fp,"CN\t%s\t%.2f\t%f\n", dist->chr, cn, fit);

        free(cn2aa_func);
        free(cn2ra_func);
        free(cn3ra_func);
        free(cn4ra_func);
        free(cn4aa_func);
    }

    peakfit_destroy(pkf);
}

static void usage(args_t *args)
{
    fprintf(stderr, "\n");
    fprintf(stderr, "About:   Detect number of chromosomal copies from Illumina's B-allele frequency (BAF)\n");
    fprintf(stderr, "Usage:   bcftools polysomy [OPTIONS] <file.vcf>\n");
    fprintf(stderr, "\n");
    fprintf(stderr, "General options:\n");
    fprintf(stderr, "    -o, --output-dir <path>        \n");
    fprintf(stderr, "    -r, --regions <region>         restrict to comma-separated list of regions\n");
    fprintf(stderr, "    -R, --regions-file <file>      restrict to regions listed in a file\n");
    fprintf(stderr, "    -s, --sample <name>            sample to analyze\n");
    fprintf(stderr, "    -t, --targets <region>         similar to -r but streams rather than index-jumps\n");
    fprintf(stderr, "    -T, --targets-file <file>      similar to -R but streams rather than index-jumps\n");
    fprintf(stderr, "    -v, --verbose                  \n");
    fprintf(stderr, "\n");
    fprintf(stderr, "Algorithm options:\n");
    fprintf(stderr, "    -b, --peak-size <float>        minimum peak size (0-1, larger is stricter) [0.1]\n");
    fprintf(stderr, "    -c, --cn-penalty <float>       penalty for increasing CN (0-1, larger is stricter) [0.7]\n");
    fprintf(stderr, "    -f, --fit-th <float>           goodness of fit threshold (>0, smaller is stricter) [3.3]\n");
    fprintf(stderr, "    -i, --include-aa               include the AA peak in CN2 and CN3 evaluation\n");
    fprintf(stderr, "    -m, --min-fraction <float>     minimum distinguishable fraction of aberrant cells [0.1]\n");
    fprintf(stderr, "    -p, --peak-symmetry <float>    peak symmetry threshold (0-1, larger is stricter) [0.5]\n");
    fprintf(stderr, "\n");
    exit(1);
}

int main_polysomy(int argc, char *argv[])
{
    args_t *args = (args_t*) calloc(1,sizeof(args_t));
    args->argc   = argc; args->argv = argv;
    args->nbins  = 150;
    args->fit_th = 3.3;
    args->cn_penalty = 0.7;
    args->peak_symmetry = 0.5;
    args->min_peak_size = 0.1;
    args->ra_rr_scaling = 1;
    args->min_fraction = 0.1;
    args->smooth = -3;

    static struct option loptions[] =
    {
        {"ra-rr-scaling",0,0,1},    // hidden option
        {"force-cn",1,0,2},         // hidden option
        {"smooth",1,0,'S'},         // hidden option
        {"nbins",1,0,'n'},          // hidden option
        {"include-aa",0,0,'i'},
        {"peak-size",1,0,'b'},
        {"min-fraction",1,0,'m'},
        {"verbose",0,0,'v'},
        {"fit-th",1,0,'f'},
        {"cn-penalty",1,0,'c'},
        {"peak-symmetry",1,0,'p'},
        {"output-dir",1,0,'o'},
        {"sample",1,0,'s'},
        {"targets",1,0,'t'},
        {"targets-file",1,0,'T'},
        {"regions",1,0,'r'},
        {"regions-file",1,0,'R'},
        {0,0,0,0}
    };
    char c, *tmp;
    while ((c = getopt_long(argc, argv, "h?o:vt:T:r:R:s:f:p:c:im:b:n:S:",loptions,NULL)) >= 0)
    {
        switch (c)
        {
            case  1 : args->ra_rr_scaling = 0; break;
            case  2 : args->force_cn = atoi(optarg); break;
            case 'n': args->nbins = atoi(optarg); break;
            case 'S': args->smooth = atoi(optarg); break;
            case 'i': args->include_aa = 1; break;
            case 'b':
                args->min_peak_size = strtod(optarg,&tmp);
                if ( *tmp ) error("Could not parse: -b %s\n", optarg);
                if ( args->min_peak_size<0 || args->min_peak_size>1 ) error("Range error: -b %s\n", optarg);
                break;
            case 'm':
                args->min_fraction = strtod(optarg,&tmp);
                if ( *tmp ) error("Could not parse: -n %s\n", optarg);
                if ( args->min_fraction<0 || args->min_fraction>1 ) error("Range error: -n %s\n", optarg);
                break;
            case 'f':
                args->fit_th = strtod(optarg,&tmp);
                if ( *tmp ) error("Could not parse: -f %s\n", optarg);
                break;
            case 'p':
                args->peak_symmetry = strtod(optarg,&tmp);
                if ( *tmp ) error("Could not parse: -p %s\n", optarg);
                break;
            case 'c':
                args->cn_penalty = strtod(optarg,&tmp);
                if ( *tmp ) error("Could not parse: -c %s\n", optarg);
                break;
            case 's': args->sample = optarg; break;
            case 't': args->targets_list = optarg; break;
            case 'T': args->targets_list = optarg; args->targets_is_file = 1; break;
            case 'r': args->regions_list = optarg; break;
            case 'R': args->regions_list = optarg; args->regions_is_file = 1; break;
            case 'o': args->output_dir = optarg; break;
            case 'v': args->verbose++; break;
            default: usage(args); break;
        }
    }
    if ( optind>=argc )
    {
        if ( !isatty(fileno((FILE *)stdin)) ) args->fname = "-";
    }
    else args->fname = argv[optind];
    if ( !args->fname ) usage(args);
    if ( !args->output_dir ) error("Missing the -o option\n");

    init_data(args);
    fit_curves(args);
    destroy_data(args);
    free(args);

    return 0;
}