libbpm: bpmnr/bpm_nr.h Source File

00001 
00013 #ifndef BPMNR_H__
00014 #define BPMNR_H__
00015 
00016 /* -----------------------------------------------------------------------------
00017    includes
00018    -------------------------------------------------------------------------- */
00019 #include <stdio.h>
00020 #include <stdlib.h>
00021 #include <math.h>
00022 #include <float.h>
00023 #include <string.h>
00024 
00025 #include <bpm/bpm_defs.h>
00026 
00027 /* -----------------------------------------------------------------------------
00028    macro definitions
00029    -------------------------------------------------------------------------- */
00030 #define GCF_ITMAX  100
00031 #define GCF_FPMIN  1.0e-30
00032 #define GCF_EPS    3.0e-7
00033 
00034 #define GSER_EPS   3.0e-7
00035 #define GSER_ITMAX 100
00036 
00037 #define RAN1_IA 16807
00038 #define RAN1_IM 2147483647
00039 #define RAN1_AM (1.0/RAN1_IM)
00040 #define RAN1_IQ 127773
00041 #define RAN1_IR 2836
00042 #define RAN1_NTAB 32
00043 #define RAN1_NDIV (1+(RAN1_IM-1)/RAN1_NTAB)
00044 #define RAN1_EPS 1.2e-7
00045 #define RAN1_RNMX (1.0-RAN1_EPS)
00046 
00055 #define __LM_BLOCKSZ__ 32
00056 #define __LM_BLOCKSZ__SQ    (__LM_BLOCKSZ__)*(__LM_BLOCKSZ__)
00057 
00058 #define LINSOLVERS_RETAIN_MEMORY
00059 #ifdef LINSOLVERS_RETAIN_MEMORY
00060 #define __LM_STATIC__ static
00061 #else
00062 #define __LM_STATIC__ /* empty */
00063 #endif /* LINSOLVERS_RETAIN_MEMORY */
00064 
00065 #define FABS(x) (((x)>=0.0)? (x) : -(x))
00066 #define CNST(x) (x)
00067 #define _LM_POW_ CNST(2.1)
00068 
00073 #define LM_DER_WORKSZ(npar, nmeas) (2*(nmeas) + 4*(npar) + (nmeas)*(npar) + (npar)*(npar))
00074 
00075 #define LM_DIF_WORKSZ(npar, nmeas) (3*(nmeas) + 4*(npar) + (nmeas)*(npar) + (npar)*(npar))
00076 
00077 #define LM_EPSILON       1E-12
00078 #define LM_ONE_THIRD     0.3333333334 /* 1.0/3.0 */
00079 
00080 #define LM_OPTS_SZ       5 /* max(4, 5) */
00081 #define LM_INFO_SZ       9
00082 #define LM_INIT_MU       1E-03
00083 #define LM_STOP_THRESH   1E-17
00084 #define LM_DIFF_DELTA    1E-06
00085 
00086 #define NR_FFTFORWARD     1  
00087 #define NR_FFTBACKWARD   -1  
00090 #define __LM_MEDIAN3(a, b, c) ( ((a) >= (b))?\
00091         ( ((c) >= (a))? (a) : ( ((c) <= (b))? (b) : (c) ) ) : \
00092         ( ((c) >= (b))? (b) : ( ((c) <= (a))? (a) : (c) ) ) )
00093 
00094 /* Some GSL defines */
00095 #define NULL_VECTOR {0, 0, 0, 0, 0}
00096 #define NULL_VECTOR_VIEW {{0, 0, 0, 0, 0}}
00097 #define NULL_MATRIX {0, 0, 0, 0, 0, 0}
00098 #define NULL_MATRIX_VIEW {{0, 0, 0, 0, 0, 0}}
00099 #define GSL_DBL_EPSILON        2.2204460492503131e-16
00100 #define OFFSET(N, incX) ((incX) > 0 ?  0 : ((N) - 1) * (-(incX)))
00101 #define GSL_MIN(a,b) ((a) < (b) ? (a) : (b))
00102 
00103 
00104 /* -----------------------------------------------------------------------------
00105    typedefs, struct, enums and other declarations
00106    -------------------------------------------------------------------------- */
00107 #ifdef __cplusplus
00108 extern "C" {
00109 #endif
00110 
00111   enum CBLAS_TRANSPOSE {CblasNoTrans=111, CblasTrans=112, CblasConjTrans=113};
00112   enum CBLAS_ORDER {CblasRowMajor=101, CblasColMajor=102};
00113   typedef  enum CBLAS_TRANSPOSE   CBLAS_TRANSPOSE_t;
00114   
00118   struct lm_fstate{ 
00119     int n, *nfev;
00120     double *hx, *x;
00121     void *adata;
00122   };
00123 
00124   /*
00125     Here are some GSL structures and the licence info. Please don't sue us!
00126     * 
00127     * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004 Gerard Jungman, Brian Gough
00128     * 
00129     * This program is free software; you can redistribute it and/or modify
00130     * it under the terms of the GNU General Public License as published by
00131     * the Free Software Foundation; either version 2 of the License, or (at
00132     * your option) any later version.
00133     * 
00134     * This program is distributed in the hope that it will be useful, but
00135     * WITHOUT ANY WARRANTY; without even the implied warranty of
00136     * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00137     * General Public License for more details.
00138     * 
00139     * You should have received a copy of the GNU General Public License
00140     * along with this program; if not, write to the Free Software
00141     * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
00142     */
00143   
00144 /* --------------------------------------
00145    Blocks */
00146   struct gsl_block_struct
00147   {
00148     size_t size;
00149     double *data;
00150   };
00151   
00152   typedef struct gsl_block_struct gsl_block;
00153   
00154 /* --------------------------------------
00155    Matrices */
00156   typedef struct 
00157   {
00158     size_t size1;
00159     size_t size2;
00160     size_t tda;
00161     double * data;
00162     gsl_block * block;
00163     int owner;
00164   } gsl_matrix;
00165   
00166   typedef struct
00167   {
00168     gsl_matrix matrix;
00169   } 
00170     _gsl_matrix_view;
00171   
00172   typedef _gsl_matrix_view gsl_matrix_view;
00173   
00174 /* --------------------------------------
00175    Vectors */
00176   typedef struct 
00177   {
00178     size_t size;
00179     size_t stride;
00180     double *data;
00181     gsl_block *block;
00182     int owner;
00183   } 
00184     gsl_vector;
00185   
00186   typedef struct
00187   {
00188     gsl_vector vector;
00189   }
00190     _gsl_vector_view;
00191   
00192   typedef _gsl_vector_view gsl_vector_view;
00193   
00194   typedef struct
00195   {
00196     gsl_vector vector;
00197   } 
00198     _gsl_vector_const_view;
00199   
00200   typedef const _gsl_vector_const_view gsl_vector_const_view;
00201   
00202   
00206   typedef struct {
00207     double re;
00208     double im;
00209   } complex_t;
00210 
00211 /* -----------------------------------------------------------------------------
00212    function prototypes and declarations
00213    -------------------------------------------------------------------------- */
00214   
00215   EXTERN double nr_gammln( double xx );
00216   EXTERN double nr_gammq( double a, double x );
00217   EXTERN int nr_gcf( double *gammcf, double a, double x, double *gln );
00218   EXTERN int nr_gser( double *gamser, double a, double x, double *gln );
00219   EXTERN int nr_fit ( double *x, double y[], int ndata, double sig[], 
00220                       int mwt, double *a, double *b, 
00221                       double *siga, double *sigb, double *chi2, double *q );
00222   EXTERN int nr_is_pow2( unsigned long n );
00223   EXTERN int nr_four1( double data[], unsigned long nn, int isign );
00224   EXTERN int nr_realft( double data[], unsigned long n, int isign );
00225   EXTERN double nr_ran1( long *idum );
00226   EXTERN int nr_seed( long seed );
00227   EXTERN double nr_ranuniform( double lower, double upper );
00228   EXTERN double nr_rangauss( double mean, double std_dev );
00229   EXTERN long bpm_rseed;
00230 
00231   /* the levenberg marquardt routines from nr_levmar.c */
00232   /* unconstrained minimization with known jacobian */
00233   EXTERN int nr_lmder( void (*func)(double *p, double *hx, int m, int n, void *adata),
00234                        void (*jacf)(double *p, double *j, int m, int n, void *adata),
00235                        double *p, double *x, int m, int n, int itmax, double *opts,
00236                        double *info, double *work, double *covar, void *adata );
00237   
00238   /* unconstrained minimization with unknown jacobian */
00239   EXTERN int nr_lmdif( void (*func)(double *p, double *hx, int m, int n, void *adata),
00240                        double *p, double *x, int m, int n, int itmax, double *opts,
00241                        double *info, double *work, double *covar, void *adata );
00242 
00243   /* box-constrained minimization with known jacobian */
00244   EXTERN int nr_lmder_bc( void (*func)(double *p, double *hx, int m, int n, void *adata),
00245                           void (*jacf)(double *p, double *j, int m, int n, void *adata),  
00246                           double *p, double *x, int m, int n, double *lb, double *ub, int itmax,
00247                           double *opts, double *info, double *work, double *covar, void *adata );
00248 
00249   /* box-constrained minimization with unknown jacobian */
00250   EXTERN int nr_lmdif_bc( void (*func)(double *p, double *hx, int m, int n, void *adata),
00251                           double *p, double *x, int m, int n, double *lb, double *ub, int itmax, 
00252                           double *opts, double *info, double *work, double *covar, void *adata );
00253 
00254   /* jacobian verification */
00255   EXTERN void nr_lmchkjac( void (*func)(double *p, double *hx, int m, int n, void *adata),
00256                            void (*jacf)(double *p, double *j, int m, int n, void *adata),
00257                            double *p, int m, int n, void *adata, double *err );
00258   
00259   /* covariance matrix */
00260   EXTERN int nr_lmcovar(double *JtJ, double *C, double sumsq, int m, int n );
00261 
00262   /* lm helper routine, LU solver for matrix inversion */
00263   EXTERN int nr_ax_eq_b_LU(double *A, double *B, double *x, int n);
00264   EXTERN void nr_trans_mat_mat_mult(double *a, double *b, int n, int m);
00265   EXTERN void nr_fdif_forw_jac_approx( void (*func)(double *p, double *hx, 
00266                                                     int m, int n, void *adata),
00267                                        double *p, double *hx, double *hxx, double delta,
00268                                        double *jac, int m, int n, void *adata );
00269   
00270   EXTERN void nr_fdif_cent_jac_approx(void (*func)(double *p, double *hx, 
00271                                                    int m, int n, void *adata),
00272                                       double *p, double *hxm, double *hxp, double delta,
00273                                       double *jac, int m, int n, void *adata );
00274 
00275   EXTERN double nr_median( int n, double *arr );
00276 
00277   EXTERN double nr_select( int k, int n, double *org_arr );
00278   
00279   /* 
00280      Here are a load of GNU Scientific Library functions. See above for the license info.
00281   */
00282 
00283   /* ----------------------------------------
00284      Matrix related */
00285   EXTERN gsl_matrix * gsl_matrix_calloc (const size_t n1, const size_t n2);
00286   EXTERN _gsl_vector_view gsl_matrix_column (gsl_matrix * m, const size_t i);
00287   EXTERN _gsl_matrix_view gsl_matrix_submatrix (gsl_matrix * m, 
00288                             const size_t i, const size_t j, 
00289                             const size_t n1, const size_t n2);
00290   EXTERN double gsl_matrix_get(const gsl_matrix * m, const size_t i, const size_t j); 
00291   EXTERN void gsl_matrix_set(gsl_matrix * m, const size_t i, const size_t j, const double x);
00292   EXTERN int gsl_matrix_swap_columns(gsl_matrix * m, const size_t i, const size_t j);
00293   EXTERN gsl_matrix * gsl_matrix_alloc (const size_t n1, const size_t n2);
00294   EXTERN _gsl_vector_const_view gsl_matrix_const_row (const gsl_matrix * m, const size_t i);
00295   EXTERN _gsl_vector_view gsl_matrix_row (gsl_matrix * m, const size_t i);
00296   EXTERN _gsl_vector_const_view gsl_matrix_const_column (const gsl_matrix * m, const size_t j);
00297   EXTERN void gsl_matrix_set_identity (gsl_matrix * m);
00298   
00299   /* ----------------------------------------
00300      Vector related */
00301   EXTERN gsl_vector *gsl_vector_calloc (const size_t n);
00302   EXTERN _gsl_vector_view gsl_vector_subvector (gsl_vector *v, size_t offset, size_t n);
00303   EXTERN double gsl_vector_get (const gsl_vector * v, const size_t i);
00304   EXTERN void gsl_vector_set (gsl_vector * v, const size_t i, double x);
00305   EXTERN int gsl_vector_swap_elements (gsl_vector * v, const size_t i, const size_t j);
00306   EXTERN _gsl_vector_const_view gsl_vector_const_subvector (const gsl_vector *v, size_t i, 
00307                                                             size_t n);
00308   EXTERN void gsl_vector_free (gsl_vector * v);
00309 
00310   /* ----------------------------------------
00311      linalg related */
00312   EXTERN int gsl_linalg_SV_solve (const gsl_matrix * U, const gsl_matrix * Q,
00313                                   const gsl_vector * S, const gsl_vector * b,
00314                                   gsl_vector * x);
00315   EXTERN int gsl_linalg_bidiag_unpack (const gsl_matrix * A, const gsl_vector * tau_U, 
00316                                        gsl_matrix * U, const gsl_vector * tau_V,
00317                                        gsl_matrix * V, gsl_vector * diag, 
00318                                        gsl_vector * superdiag);
00319   EXTERN int gsl_linalg_householder_hm (double tau, const gsl_vector * v, gsl_matrix * A);
00320   EXTERN int gsl_linalg_bidiag_unpack2 (gsl_matrix * A, gsl_vector * tau_U, 
00321                                         gsl_vector * tau_V, gsl_matrix * V);
00322   EXTERN int gsl_linalg_householder_hm1 (double tau, gsl_matrix * A);
00323   EXTERN void create_givens (const double a, const double b, double *c, double *s);
00324   EXTERN double gsl_linalg_householder_transform (gsl_vector * v);
00325   EXTERN int gsl_linalg_householder_mh (double tau, const gsl_vector * v, gsl_matrix * A);
00326   EXTERN int gsl_linalg_SV_solve (const gsl_matrix * U,
00327                                   const gsl_matrix * V,
00328                                   const gsl_vector * S,
00329                                   const gsl_vector * b, gsl_vector * x);
00330   EXTERN void chop_small_elements (gsl_vector * d, gsl_vector * f);
00331   EXTERN void qrstep (gsl_vector * d, gsl_vector * f, gsl_matrix * U, gsl_matrix * V);
00332   EXTERN double trailing_eigenvalue (const gsl_vector * d, const gsl_vector * f);
00333   EXTERN void create_schur (double d0, double f0, double d1, double * c, double * s);
00334   EXTERN void svd2 (gsl_vector * d, gsl_vector * f, gsl_matrix * U, gsl_matrix * V);
00335   EXTERN void chase_out_intermediate_zero (gsl_vector * d, gsl_vector * f, gsl_matrix * U, 
00336                                            size_t k0);
00337   EXTERN void chase_out_trailing_zero (gsl_vector * d, gsl_vector * f, gsl_matrix * V);
00338   EXTERN int gsl_isnan (const double x);
00339 
00340   /* ----------------------------------------
00341      blas/cblas related */
00342   EXTERN double gsl_blas_dnrm2 (const gsl_vector * X);
00343   EXTERN double cblas_dnrm2(const int N, const double *X, const int incX);
00344   EXTERN void gsl_blas_dscal  (double alpha, gsl_vector * X);
00345   EXTERN void cblas_dscal(const int N, const double alpha, double *X, const int incX);
00346   EXTERN void cblas_dgemv (const enum CBLAS_ORDER order, const enum CBLAS_TRANSPOSE TransA,
00347                            const int M, const int N, const double alpha, const double *A,
00348                            const int lda, const double *X, const int incX,
00349                            const double beta, double *Y, const int incY);
00350 
00351   /* ----------------------------------------
00352      block related */
00353   EXTERN gsl_block *gsl_block_alloc (const size_t n);
00354   EXTERN void gsl_block_free (gsl_block * b);
00355 
00356 
00357 
00358   /* Some complex functions */
00359   EXTERN complex_t complex( double re, double im );
00360   EXTERN double    c_real( complex_t z );
00361   EXTERN double    c_imag( complex_t z );
00362   EXTERN complex_t c_conj( complex_t z );
00363   EXTERN complex_t c_neg( complex_t z );
00364   EXTERN complex_t c_sum( complex_t z1, complex_t z2 );
00365   EXTERN complex_t c_diff( complex_t z1, complex_t z2 );
00366   EXTERN complex_t c_mult( complex_t z1, complex_t z2 );
00367   EXTERN complex_t c_div( complex_t z1, complex_t z2 );
00368   EXTERN complex_t c_scale( double r, complex_t z );
00369   EXTERN complex_t c_sqr( complex_t z );
00370   EXTERN complex_t c_sqrt( complex_t z );
00371   EXTERN double    c_norm2( complex_t z );
00372   EXTERN double    c_abs( complex_t z );
00373   EXTERN double    c_arg( complex_t z );
00374   EXTERN complex_t c_exp( complex_t z );
00375   EXTERN int       c_isequal( complex_t z1, complex_t z2 );
00376 
00377 
00382   EXTERN double nr_quadinterpol( double x, 
00383                                  double x1, double x2, double x3,
00384                                  double y1, double y2, double y3 );
00385   
00386   
00388   EXTERN double sinc( double x );
00389 
00391   EXTERN double lanczos( double x, int a );
00392 
00394   EXTERN double dround( double x );
00395 
00396 
00397 #ifdef __cplusplus
00398 }
00399 #endif
00400 
00401 #endif /* #ifndef BPMNR_H__ */
00402 
00403 /* ================================ end of file ============================= */