diff --git a/include/gsl/gsl_deriv.h b/include/gsl/gsl_deriv.h new file mode 100644 index 0000000..7f4694f --- /dev/null +++ b/include/gsl/gsl_deriv.h @@ -0,0 +1,50 @@ +/* deriv/gsl_deriv.h + * + * 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 3 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +#ifndef __GSL_DERIV_H__ +#define __GSL_DERIV_H__ +#include + +#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 + +int gsl_deriv_central (const gsl_function *f, + double x, double h, + double *result, double *abserr); + +int gsl_deriv_backward (const gsl_function *f, + double x, double h, + double *result, double *abserr); + +int gsl_deriv_forward (const gsl_function *f, + double x, double h, + double *result, double *abserr); + +__END_DECLS + +#endif /* __GSL_DERIV_H__ */ diff --git a/include/gsl/gsl_fit.h b/include/gsl/gsl_fit.h new file mode 100644 index 0000000..de83a41 --- /dev/null +++ b/include/gsl/gsl_fit.h @@ -0,0 +1,85 @@ +/* fit/gsl_fit.h + * + * Copyright (C) 2000, 2007 Brian Gough + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 3 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +#ifndef __GSL_FIT_H__ +#define __GSL_FIT_H__ + +#include +#include + +#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 + +int gsl_fit_linear (const double * x, const size_t xstride, + const double * y, const size_t ystride, + const size_t n, + double * c0, double * c1, + double * cov00, double * cov01, double * cov11, + double * sumsq); + + +int gsl_fit_wlinear (const double * x, const size_t xstride, + const double * w, const size_t wstride, + const double * y, const size_t ystride, + const size_t n, + double * c0, double * c1, + double * cov00, double * cov01, double * cov11, + double * chisq); + +int +gsl_fit_linear_est (const double x, + const double c0, const double c1, + const double cov00, const double cov01, const double cov11, + double *y, double *y_err); + + +int gsl_fit_mul (const double * x, const size_t xstride, + const double * y, const size_t ystride, + const size_t n, + double * c1, + double * cov11, + double * sumsq); + +int gsl_fit_wmul (const double * x, const size_t xstride, + const double * w, const size_t wstride, + const double * y, const size_t ystride, + const size_t n, + double * c1, + double * cov11, + double * sumsq); + + +int +gsl_fit_mul_est (const double x, + const double c1, + const double cov11, + double *y, double *y_err); + +__END_DECLS + +#endif /* __GSL_FIT_H__ */ diff --git a/include/gsl/gsl_machine.h b/include/gsl/gsl_machine.h new file mode 100644 index 0000000..c44ffc2 --- /dev/null +++ b/include/gsl/gsl_machine.h @@ -0,0 +1,104 @@ +/* Author: B. Gough and G. Jungman */ +#ifndef __GSL_MACHINE_H__ +#define __GSL_MACHINE_H__ + +#include +#include + +/* magic constants; mostly for the benefit of the implementation */ + +/* -*-MACHINE CONSTANTS-*- + * + * PLATFORM: Whiz-O-Matic 9000 + * FP_PLATFORM: IEEE-Virtual + * HOSTNAME: nnn.lanl.gov + * DATE: Fri Nov 20 17:53:26 MST 1998 + */ +#define GSL_DBL_EPSILON 2.2204460492503131e-16 +#define GSL_SQRT_DBL_EPSILON 1.4901161193847656e-08 +#define GSL_ROOT3_DBL_EPSILON 6.0554544523933429e-06 +#define GSL_ROOT4_DBL_EPSILON 1.2207031250000000e-04 +#define GSL_ROOT5_DBL_EPSILON 7.4009597974140505e-04 +#define GSL_ROOT6_DBL_EPSILON 2.4607833005759251e-03 +#define GSL_LOG_DBL_EPSILON (-3.6043653389117154e+01) + +#define GSL_DBL_MIN 2.2250738585072014e-308 +#define GSL_SQRT_DBL_MIN 1.4916681462400413e-154 +#define GSL_ROOT3_DBL_MIN 2.8126442852362996e-103 +#define GSL_ROOT4_DBL_MIN 1.2213386697554620e-77 +#define GSL_ROOT5_DBL_MIN 2.9476022969691763e-62 +#define GSL_ROOT6_DBL_MIN 5.3034368905798218e-52 +#define GSL_LOG_DBL_MIN (-7.0839641853226408e+02) + +#define GSL_DBL_MAX 1.7976931348623157e+308 +#define GSL_SQRT_DBL_MAX 1.3407807929942596e+154 +#define GSL_ROOT3_DBL_MAX 5.6438030941222897e+102 +#define GSL_ROOT4_DBL_MAX 1.1579208923731620e+77 +#define GSL_ROOT5_DBL_MAX 4.4765466227572707e+61 +#define GSL_ROOT6_DBL_MAX 2.3756689782295612e+51 +#define GSL_LOG_DBL_MAX 7.0978271289338397e+02 + +#define GSL_FLT_EPSILON 1.1920928955078125e-07 +#define GSL_SQRT_FLT_EPSILON 3.4526698300124393e-04 +#define GSL_ROOT3_FLT_EPSILON 4.9215666011518501e-03 +#define GSL_ROOT4_FLT_EPSILON 1.8581361171917516e-02 +#define GSL_ROOT5_FLT_EPSILON 4.1234622211652937e-02 +#define GSL_ROOT6_FLT_EPSILON 7.0153878019335827e-02 +#define GSL_LOG_FLT_EPSILON (-1.5942385152878742e+01) + +#define GSL_FLT_MIN 1.1754943508222875e-38 +#define GSL_SQRT_FLT_MIN 1.0842021724855044e-19 +#define GSL_ROOT3_FLT_MIN 2.2737367544323241e-13 +#define GSL_ROOT4_FLT_MIN 3.2927225399135965e-10 +#define GSL_ROOT5_FLT_MIN 2.5944428542140822e-08 +#define GSL_ROOT6_FLT_MIN 4.7683715820312542e-07 +#define GSL_LOG_FLT_MIN (-8.7336544750553102e+01) + +#define GSL_FLT_MAX 3.4028234663852886e+38 +#define GSL_SQRT_FLT_MAX 1.8446743523953730e+19 +#define GSL_ROOT3_FLT_MAX 6.9814635196223242e+12 +#define GSL_ROOT4_FLT_MAX 4.2949672319999986e+09 +#define GSL_ROOT5_FLT_MAX 5.0859007855960041e+07 +#define GSL_ROOT6_FLT_MAX 2.6422459233807749e+06 +#define GSL_LOG_FLT_MAX 8.8722839052068352e+01 + +#define GSL_SFLT_EPSILON 4.8828125000000000e-04 +#define GSL_SQRT_SFLT_EPSILON 2.2097086912079612e-02 +#define GSL_ROOT3_SFLT_EPSILON 7.8745065618429588e-02 +#define GSL_ROOT4_SFLT_EPSILON 1.4865088937534013e-01 +#define GSL_ROOT5_SFLT_EPSILON 2.1763764082403100e-01 +#define GSL_ROOT6_SFLT_EPSILON 2.8061551207734325e-01 +#define GSL_LOG_SFLT_EPSILON (-7.6246189861593985e+00) + +/* !MACHINE CONSTANTS! */ + + +/* a little internal backwards compatibility */ +#define GSL_MACH_EPS GSL_DBL_EPSILON + + + +/* Here are the constants related to or derived from + * machine constants. These are not to be confused with + * the constants that define various precision levels + * for the precision/error system. + * + * This information is determined at configure time + * and is platform dependent. Edit at your own risk. + * + * PLATFORM: WHIZ-O-MATIC + * CONFIG-DATE: Thu Nov 19 19:27:18 MST 1998 + * CONFIG-HOST: nnn.lanl.gov + */ + +/* machine precision constants */ +/* #define GSL_MACH_EPS 1.0e-15 */ +#define GSL_SQRT_MACH_EPS 3.2e-08 +#define GSL_ROOT3_MACH_EPS 1.0e-05 +#define GSL_ROOT4_MACH_EPS 0.000178 +#define GSL_ROOT5_MACH_EPS 0.00100 +#define GSL_ROOT6_MACH_EPS 0.00316 +#define GSL_LOG_MACH_EPS (-34.54) + + +#endif /* __GSL_MACHINE_H__ */ diff --git a/include/gsl/gsl_math.h b/include/gsl/gsl_math.h new file mode 100644 index 0000000..9687464 --- /dev/null +++ b/include/gsl/gsl_math.h @@ -0,0 +1,164 @@ +/* gsl_math.h + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004, 2007 Gerard Jungman, Brian Gough + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 3 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +#ifndef __GSL_MATH_H__ +#define __GSL_MATH_H__ +#include +//#include +//#include +#include +//#include +//#include +//#include +//#include + +#ifndef M_E +#define M_E 2.71828182845904523536028747135 /* e */ +#endif + +#ifndef M_LOG2E +#define M_LOG2E 1.44269504088896340735992468100 /* log_2 (e) */ +#endif + +#ifndef M_LOG10E +#define M_LOG10E 0.43429448190325182765112891892 /* log_10 (e) */ +#endif + +#ifndef M_SQRT2 +#define M_SQRT2 1.41421356237309504880168872421 /* sqrt(2) */ +#endif + +#ifndef M_SQRT1_2 +#define M_SQRT1_2 0.70710678118654752440084436210 /* sqrt(1/2) */ +#endif + + +#ifndef M_SQRT3 +#define M_SQRT3 1.73205080756887729352744634151 /* sqrt(3) */ +#endif + +#ifndef M_PI +#define M_PI 3.14159265358979323846264338328 /* pi */ +#endif + +#ifndef M_PI_2 +#define M_PI_2 1.57079632679489661923132169164 /* pi/2 */ +#endif + +#ifndef M_PI_4 +#define M_PI_4 0.78539816339744830961566084582 /* pi/4 */ +#endif + +#ifndef M_SQRTPI +#define M_SQRTPI 1.77245385090551602729816748334 /* sqrt(pi) */ +#endif + +#ifndef M_2_SQRTPI +#define M_2_SQRTPI 1.12837916709551257389615890312 /* 2/sqrt(pi) */ +#endif + +#ifndef M_1_PI +#define M_1_PI 0.31830988618379067153776752675 /* 1/pi */ +#endif + +#ifndef M_2_PI +#define M_2_PI 0.63661977236758134307553505349 /* 2/pi */ +#endif + +#ifndef M_LN10 +#define M_LN10 2.30258509299404568401799145468 /* ln(10) */ +#endif + +#ifndef M_LN2 +#define M_LN2 0.69314718055994530941723212146 /* ln(2) */ +#endif + +#ifndef M_LNPI +#define M_LNPI 1.14472988584940017414342735135 /* ln(pi) */ +#endif + +#ifndef M_EULER +#define M_EULER 0.57721566490153286060651209008 /* Euler constant */ +#endif + +#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 + +/* other needlessly compulsive abstractions */ + +#define GSL_IS_ODD(n) ((n) & 1) +#define GSL_IS_EVEN(n) (!(GSL_IS_ODD(n))) +#define GSL_SIGN(x) ((x) >= 0.0 ? 1 : -1) + +/* Return nonzero if x is a real number, i.e. non NaN or infinite. */ +#define GSL_IS_REAL(x) (gsl_finite(x)) + +/* Definition of an arbitrary function with parameters */ + +struct gsl_function_struct +{ + double (* function) (double x, void * params); + void * params; +}; + +typedef struct gsl_function_struct gsl_function ; + +#define GSL_FN_EVAL(F,x) (*((F)->function))(x,(F)->params) + +/* Definition of an arbitrary function returning two values, r1, r2 */ + +struct gsl_function_fdf_struct +{ + double (* f) (double x, void * params); + double (* df) (double x, void * params); + void (* fdf) (double x, void * params, double * f, double * df); + void * params; +}; + +typedef struct gsl_function_fdf_struct gsl_function_fdf ; + +#define GSL_FN_FDF_EVAL_F(FDF,x) (*((FDF)->f))(x,(FDF)->params) +#define GSL_FN_FDF_EVAL_DF(FDF,x) (*((FDF)->df))(x,(FDF)->params) +#define GSL_FN_FDF_EVAL_F_DF(FDF,x,y,dy) (*((FDF)->fdf))(x,(FDF)->params,(y),(dy)) + + +/* Definition of an arbitrary vector-valued function with parameters */ + +struct gsl_function_vec_struct +{ + int (* function) (double x, double y[], void * params); + void * params; +}; + +typedef struct gsl_function_vec_struct gsl_function_vec ; + +#define GSL_FN_VEC_EVAL(F,x,y) (*((F)->function))(x,y,(F)->params) + +__END_DECLS + +#endif /* __GSL_MATH_H__ */ diff --git a/src/core/deriv.c b/src/core/deriv.c new file mode 100644 index 0000000..aa7e2e6 --- /dev/null +++ b/src/core/deriv.c @@ -0,0 +1,176 @@ +/* deriv/deriv.c + * + * Copyright (C) 2004, 2007 Brian Gough + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 3 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +#include +#include +#include + +static void +central_deriv (const gsl_function * f, double x, double h, + double *result, double *abserr_round, double *abserr_trunc) +{ + /* Compute the derivative using the 5-point rule (x-h, x-h/2, x, + x+h/2, x+h). Note that the central point is not used. + + Compute the error using the difference between the 5-point and + the 3-point rule (x-h,x,x+h). Again the central point is not + used. */ + + double fm1 = GSL_FN_EVAL (f, x - h); + double fp1 = GSL_FN_EVAL (f, x + h); + + double fmh = GSL_FN_EVAL (f, x - h / 2); + double fph = GSL_FN_EVAL (f, x + h / 2); + + double r3 = 0.5 * (fp1 - fm1); + double r5 = (4.0 / 3.0) * (fph - fmh) - (1.0 / 3.0) * r3; + + double e3 = (fabs (fp1) + fabs (fm1)) * GSL_DBL_EPSILON; + double e5 = 2.0 * (fabs (fph) + fabs (fmh)) * GSL_DBL_EPSILON + e3; + + /* The next term is due to finite precision in x+h = O (eps * x) */ + + double dy = max (fabs (r3 / h), fabs (r5 / h)) *(fabs (x) / h) * GSL_DBL_EPSILON; + + /* The truncation error in the r5 approximation itself is O(h^4). + However, for safety, we estimate the error from r5-r3, which is + O(h^2). By scaling h we will minimise this estimated error, not + the actual truncation error in r5. */ + + *result = r5 / h; + *abserr_trunc = fabs ((r5 - r3) / h); /* Estimated truncation error O(h^2) */ + *abserr_round = fabs (e5 / h) + dy; /* Rounding error (cancellations) */ +} + +int +gsl_deriv_central (const gsl_function * f, double x, double h, + double *result, double *abserr) +{ + double r_0, round, trunc, error; + central_deriv (f, x, h, &r_0, &round, &trunc); + error = round + trunc; + + if (round < trunc && (round > 0 && trunc > 0)) + { + double r_opt, round_opt, trunc_opt, error_opt; + + /* Compute an optimised stepsize to minimize the total error, + using the scaling of the truncation error (O(h^2)) and + rounding error (O(1/h)). */ + + double h_opt = h * pow (round / (2.0 * trunc), 1.0 / 3.0); + central_deriv (f, x, h_opt, &r_opt, &round_opt, &trunc_opt); + error_opt = round_opt + trunc_opt; + + /* Check that the new error is smaller, and that the new derivative + is consistent with the error bounds of the original estimate. */ + + if (error_opt < error && fabs (r_opt - r_0) < 4.0 * error) + { + r_0 = r_opt; + error = error_opt; + } + } + + *result = r_0; + *abserr = error; + + return 0; +} + + +static void +forward_deriv (const gsl_function * f, double x, double h, + double *result, double *abserr_round, double *abserr_trunc) +{ + /* Compute the derivative using the 4-point rule (x+h/4, x+h/2, + x+3h/4, x+h). + + Compute the error using the difference between the 4-point and + the 2-point rule (x+h/2,x+h). */ + + double f1 = GSL_FN_EVAL (f, x + h / 4.0); + double f2 = GSL_FN_EVAL (f, x + h / 2.0); + double f3 = GSL_FN_EVAL (f, x + (3.0 / 4.0) * h); + double f4 = GSL_FN_EVAL (f, x + h); + + double r2 = 2.0*(f4 - f2); + double r4 = (22.0 / 3.0) * (f4 - f3) - (62.0 / 3.0) * (f3 - f2) + + (52.0 / 3.0) * (f2 - f1); + + /* Estimate the rounding error for r4 */ + + double e4 = 2 * 20.67 * (fabs (f4) + fabs (f3) + fabs (f2) + fabs (f1)) * GSL_DBL_EPSILON; + + /* The next term is due to finite precision in x+h = O (eps * x) */ + + double dy = max (fabs (r2 / h), fabs (r4 / h)) * fabs (x / h) * GSL_DBL_EPSILON; + + /* The truncation error in the r4 approximation itself is O(h^3). + However, for safety, we estimate the error from r4-r2, which is + O(h). By scaling h we will minimise this estimated error, not + the actual truncation error in r4. */ + + *result = r4 / h; + *abserr_trunc = fabs ((r4 - r2) / h); /* Estimated truncation error O(h) */ + *abserr_round = fabs (e4 / h) + dy; +} + +int +gsl_deriv_forward (const gsl_function * f, double x, double h, + double *result, double *abserr) +{ + double r_0, round, trunc, error; + forward_deriv (f, x, h, &r_0, &round, &trunc); + error = round + trunc; + + if (round < trunc && (round > 0 && trunc > 0)) + { + double r_opt, round_opt, trunc_opt, error_opt; + + /* Compute an optimised stepsize to minimize the total error, + using the scaling of the estimated truncation error (O(h)) and + rounding error (O(1/h)). */ + + double h_opt = h * pow (round / (trunc), 1.0 / 2.0); + forward_deriv (f, x, h_opt, &r_opt, &round_opt, &trunc_opt); + error_opt = round_opt + trunc_opt; + + /* Check that the new error is smaller, and that the new derivative + is consistent with the error bounds of the original estimate. */ + + if (error_opt < error && fabs (r_opt - r_0) < 4.0 * error) + { + r_0 = r_opt; + error = error_opt; + } + } + + *result = r_0; + *abserr = error; + + return 0; +} + +int +gsl_deriv_backward (const gsl_function * f, double x, double h, + double *result, double *abserr) +{ + return gsl_deriv_forward (f, x, -h, result, abserr); +} diff --git a/src/core/drand48.cpp b/src/core/drand48.cpp new file mode 100644 index 0000000..26e795b --- /dev/null +++ b/src/core/drand48.cpp @@ -0,0 +1,8 @@ +#include + +static std::default_random_engine generator; +static std::uniform_real_distribution distribution(0.0,1.0); + +double drand48() { + return distribution(generator); +} \ No newline at end of file diff --git a/src/core/ffs.c b/src/core/ffs.c new file mode 100644 index 0000000..8a2d8a1 --- /dev/null +++ b/src/core/ffs.c @@ -0,0 +1,39 @@ +/* Copyright (C) 1991-2018 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Torbjorn Granlund (tege@sics.se). + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library 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 + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + . */ + +/* Find the first bit set in I. */ +int find_first_bit_set (int i) +{ + static const unsigned char table[] = + { + 0,1,2,2,3,3,3,3,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5, + 6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6, + 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7, + 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7, + 8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8, + 8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8, + 8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8, + 8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8 + }; + unsigned int a; + unsigned int x = i & -i; + + a = x <= 0xffff ? (x <= 0xff ? 0 : 8) : (x <= 0xffffff ? 16 : 24); + + return table[x >> a] + a; +} diff --git a/src/core/linear.c b/src/core/linear.c new file mode 100644 index 0000000..976209b --- /dev/null +++ b/src/core/linear.c @@ -0,0 +1,344 @@ +/* fit/linear.c + * + * Copyright (C) 2000, 2007 Brian Gough + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 3 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +#include + +/* Fit the data (x_i, y_i) to the linear relationship + + Y = c0 + c1 x + + returning, + + c0, c1 -- coefficients + cov00, cov01, cov11 -- variance-covariance matrix of c0 and c1, + sumsq -- sum of squares of residuals + + This fit can be used in the case where the errors for the data are + uknown, but assumed equal for all points. The resulting + variance-covariance matrix estimates the error in the coefficients + from the observed variance of the points around the best fit line. +*/ + +int +gsl_fit_linear (const double *x, const size_t xstride, + const double *y, const size_t ystride, + const size_t n, + double *c0, double *c1, + double *cov_00, double *cov_01, double *cov_11, double *sumsq) +{ + double m_x = 0, m_y = 0, m_dx2 = 0, m_dxdy = 0; + + size_t i; + + for (i = 0; i < n; i++) + { + m_x += (x[i * xstride] - m_x) / (i + 1.0); + m_y += (y[i * ystride] - m_y) / (i + 1.0); + } + + for (i = 0; i < n; i++) + { + const double dx = x[i * xstride] - m_x; + const double dy = y[i * ystride] - m_y; + + m_dx2 += (dx * dx - m_dx2) / (i + 1.0); + m_dxdy += (dx * dy - m_dxdy) / (i + 1.0); + } + + /* In terms of y = a + b x */ + + { + double s2 = 0, d2 = 0; + double b = m_dxdy / m_dx2; + double a = m_y - m_x * b; + + *c0 = a; + *c1 = b; + + /* Compute chi^2 = \sum (y_i - (a + b * x_i))^2 */ + + for (i = 0; i < n; i++) + { + const double dx = x[i * xstride] - m_x; + const double dy = y[i * ystride] - m_y; + const double d = dy - b * dx; + d2 += d * d; + } + + s2 = d2 / (n - 2.0); /* chisq per degree of freedom */ + + *cov_00 = s2 * (1.0 / n) * (1 + m_x * m_x / m_dx2); + *cov_11 = s2 * 1.0 / (n * m_dx2); + + *cov_01 = s2 * (-m_x) / (n * m_dx2); + + *sumsq = d2; + } + + return 0; +} + + +/* Fit the weighted data (x_i, w_i, y_i) to the linear relationship + + Y = c0 + c1 x + + returning, + + c0, c1 -- coefficients + s0, s1 -- the standard deviations of c0 and c1, + r -- the correlation coefficient between c0 and c1, + chisq -- weighted sum of squares of residuals */ + +int +gsl_fit_wlinear (const double *x, const size_t xstride, + const double *w, const size_t wstride, + const double *y, const size_t ystride, + const size_t n, + double *c0, double *c1, + double *cov_00, double *cov_01, double *cov_11, + double *chisq) +{ + + /* compute the weighted means and weighted deviations from the means */ + + /* wm denotes a "weighted mean", wm(f) = (sum_i w_i f_i) / (sum_i w_i) */ + + double W = 0, wm_x = 0, wm_y = 0, wm_dx2 = 0, wm_dxdy = 0; + + size_t i; + + for (i = 0; i < n; i++) + { + const double wi = w[i * wstride]; + + if (wi > 0) + { + W += wi; + wm_x += (x[i * xstride] - wm_x) * (wi / W); + wm_y += (y[i * ystride] - wm_y) * (wi / W); + } + } + + W = 0; /* reset the total weight */ + + for (i = 0; i < n; i++) + { + const double wi = w[i * wstride]; + + if (wi > 0) + { + const double dx = x[i * xstride] - wm_x; + const double dy = y[i * ystride] - wm_y; + + W += wi; + wm_dx2 += (dx * dx - wm_dx2) * (wi / W); + wm_dxdy += (dx * dy - wm_dxdy) * (wi / W); + } + } + + /* In terms of y = a + b x */ + + { + double d2 = 0; + double b = wm_dxdy / wm_dx2; + double a = wm_y - wm_x * b; + + *c0 = a; + *c1 = b; + + *cov_00 = (1 / W) * (1 + wm_x * wm_x / wm_dx2); + *cov_11 = 1 / (W * wm_dx2); + + *cov_01 = -wm_x / (W * wm_dx2); + + /* Compute chi^2 = \sum w_i (y_i - (a + b * x_i))^2 */ + + for (i = 0; i < n; i++) + { + const double wi = w[i * wstride]; + + if (wi > 0) + { + const double dx = x[i * xstride] - wm_x; + const double dy = y[i * ystride] - wm_y; + const double d = dy - b * dx; + d2 += wi * d * d; + } + } + + *chisq = d2; + } + + return 0; +} + + + +int +gsl_fit_linear_est (const double x, + const double c0, const double c1, + const double cov00, const double cov01, const double cov11, + double *y, double *y_err) +{ + *y = c0 + c1 * x; + *y_err = sqrt (cov00 + x * (2 * cov01 + cov11 * x)); + return 0; +} + + +int +gsl_fit_mul (const double *x, const size_t xstride, + const double *y, const size_t ystride, + const size_t n, + double *c1, double *cov_11, double *sumsq) +{ + double m_x = 0, m_y = 0, m_dx2 = 0, m_dxdy = 0; + + size_t i; + + for (i = 0; i < n; i++) + { + m_x += (x[i * xstride] - m_x) / (i + 1.0); + m_y += (y[i * ystride] - m_y) / (i + 1.0); + } + + for (i = 0; i < n; i++) + { + const double dx = x[i * xstride] - m_x; + const double dy = y[i * ystride] - m_y; + + m_dx2 += (dx * dx - m_dx2) / (i + 1.0); + m_dxdy += (dx * dy - m_dxdy) / (i + 1.0); + } + + /* In terms of y = b x */ + + { + double s2 = 0, d2 = 0; + double b = (m_x * m_y + m_dxdy) / (m_x * m_x + m_dx2); + + *c1 = b; + + /* Compute chi^2 = \sum (y_i - b * x_i)^2 */ + + for (i = 0; i < n; i++) + { + const double dx = x[i * xstride] - m_x; + const double dy = y[i * ystride] - m_y; + const double d = (m_y - b * m_x) + dy - b * dx; + d2 += d * d; + } + + s2 = d2 / (n - 1.0); /* chisq per degree of freedom */ + + *cov_11 = s2 * 1.0 / (n * (m_x * m_x + m_dx2)); + + *sumsq = d2; + } + + return 0; +} + + +int +gsl_fit_wmul (const double *x, const size_t xstride, + const double *w, const size_t wstride, + const double *y, const size_t ystride, + const size_t n, + double *c1, double *cov_11, double *chisq) +{ + + /* compute the weighted means and weighted deviations from the means */ + + /* wm denotes a "weighted mean", wm(f) = (sum_i w_i f_i) / (sum_i w_i) */ + + double W = 0, wm_x = 0, wm_y = 0, wm_dx2 = 0, wm_dxdy = 0; + + size_t i; + + for (i = 0; i < n; i++) + { + const double wi = w[i * wstride]; + + if (wi > 0) + { + W += wi; + wm_x += (x[i * xstride] - wm_x) * (wi / W); + wm_y += (y[i * ystride] - wm_y) * (wi / W); + } + } + + W = 0; /* reset the total weight */ + + for (i = 0; i < n; i++) + { + const double wi = w[i * wstride]; + + if (wi > 0) + { + const double dx = x[i * xstride] - wm_x; + const double dy = y[i * ystride] - wm_y; + + W += wi; + wm_dx2 += (dx * dx - wm_dx2) * (wi / W); + wm_dxdy += (dx * dy - wm_dxdy) * (wi / W); + } + } + + /* In terms of y = b x */ + + { + double d2 = 0; + double b = (wm_x * wm_y + wm_dxdy) / (wm_x * wm_x + wm_dx2); + + *c1 = b; + + *cov_11 = 1 / (W * (wm_x * wm_x + wm_dx2)); + + /* Compute chi^2 = \sum w_i (y_i - b * x_i)^2 */ + + for (i = 0; i < n; i++) + { + const double wi = w[i * wstride]; + + if (wi > 0) + { + const double dx = x[i * xstride] - wm_x; + const double dy = y[i * ystride] - wm_y; + const double d = (wm_y - b * wm_x) + (dy - b * dx); + d2 += wi * d * d; + } + } + + *chisq = d2; + } + + return 0; +} + +int +gsl_fit_mul_est (const double x, + const double c1, const double cov11, + double *y, double *y_err) +{ + *y = c1 * x; + *y_err = sqrt (cov11) * fabs (x); + return 0; +} diff --git a/test/core/CMakeLists.txt b/test/core/CMakeLists.txt new file mode 100644 index 0000000..5b43892 --- /dev/null +++ b/test/core/CMakeLists.txt @@ -0,0 +1,6 @@ +set(SOURCES + test_discrete_set.cpp + ) + +add_executable(test_discrete_set ${SOURCES}) +target_link_libraries(test_discrete_set ${PROJECT_NAME}_static) \ No newline at end of file