libtpcmodel.h

Go to the documentation of this file.

 
#ifndef _LIBTPCMODEL_H
#define _LIBTPCMODEL_H
/*****************************************************************************/
 
/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
//#include <omp.h>
/*****************************************************************************/
#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <string.h>
#include <strings.h>
#include <ctype.h>
#include <math.h>
#include <float.h>
#include <sys/time.h>
#include <time.h>
/*****************************************************************************/
#include "libtpcmisc.h"
/*****************************************************************************/
 
/*****************************************************************************/
#ifndef MAX_PARAMETERS
#define MAX_PARAMETERS 50
#endif
#ifndef MAX_PARAMS
#define MAX_PARAMS MAX_PARAMETERS
#endif
/*****************************************************************************/
 
/*****************************************************************************/
#define TPCCLIB_MERTWI_NN 312
#define TPCCLIB_MERTWI_A UINT64_C(0xB5026F5AA96619E9)
typedef struct MERTWI {
  unsigned int n;
  unsigned int m;
  uint64_t a;
  uint64_t um;
  uint64_t lm;
  uint64_t mt[TPCCLIB_MERTWI_NN]; 
  uint64_t mti;
} MERTWI;
/* mertwi */
extern uint32_t mertwiSeed32(void);
extern uint64_t mertwiSeed64(void);
extern void mertwiInit(MERTWI *mt);
extern void mertwiInitWithSeed64(MERTWI *mt, uint64_t seed);
extern void mertwiInitByArray64(MERTWI *mt, uint64_t init_key[], uint64_t key_length);
extern uint64_t mertwiRandomInt64(MERTWI *mt);
extern int64_t mertwiRandomInt63(MERTWI *mt);
extern double mertwiRandomDouble1(MERTWI *mt);
extern double mertwiRandomDouble2(MERTWI *mt);
extern double mertwiRandomDouble3(MERTWI *mt);
/*****************************************************************************/
 
/*****************************************************************************/
typedef enum {
   BOBYQA_INVALID_ARGS = -1,
   BOBYQA_OUT_OF_MEMORY = -2,
   BOBYQA_ROUNDOFF_LIMITED = -3,
   BOBYQA_FAIL = -4, /* generic fail code */
   BOBYQA_SUCCESS = 0, /* generic success code */
   BOBYQA_MINF_MAX_REACHED = 1,
   BOBYQA_FTOL_REACHED = 2,
   BOBYQA_XTOL_REACHED = 3,
   BOBYQA_MAXEVAL_REACHED = 4,
   BOBYQA_RELFTOL_REACHED = 5,
   BOBYQA_ABSFTOL_REACHED = 6
} bobyqa_result;
 
typedef double (*bobyqa_func)(int n, const double *x, void *func_data);
 
typedef struct {
  int n;
  int npt;
  double *x;
  int x_size;
  double *xscale;
  int xscale_size;
  int nfull;
  double *xfull;
  int *xplace; 
  double *xl;
  int xl_size;
  double *xu;
  int xu_size;
  double rhobeg;
  double rhoend; 
  double minf_max;
  double ftol_rel;
  double ftol_abs;
  int maxeval;
  int nevals;
 
  bobyqa_func objf;
  void *objf_data;
  double minf;
 
  double *wmptr;
  double *lwmptr;
  int *liwmptr;
 
  double *xbase; 
  int xbase_size;
  double *xpt;
  int xpt_size;
  double *fval;
  int fval_size;
  double *xopt;
  int xopt_size;
  double *gopt;
  int gopt_size;
  double *hq;
  int hq_size;
  double *pq;
  int pq_size;
  double *bmat;
  int bmat_size;
  double *zmat; 
  int zmat_size; 
  int ndim;
  double *sl;
  int sl_size;
  double *su;
  int su_size;
  double *xnew; 
  int xnew_size; 
  double *xalt;
  int xalt_size;
  double *dtrial;
  int dtrial_size;
  double *vlag;
  int vlag_size;
 
  double *w2npt;
  int w2npt_size;
  double *wndim;
  int wndim_size;
  double *wn;
  int wn_size;
  double *gnew;
  int gnew_size;
  double *xbdi;
  int xbdi_size;
  double *s;
  int s_size;
  double *hs;
  int hs_size;
  double *hred;
  int hred_size;
  double *glag;
  int glag_size;
  double *hcol;
  int hcol_size;
  double *ccstep;
  int ccstep_size;
 
  int verbose;
 
 
  double _crvmin;
 
  int ntrits;
  double rho;
  int nresc;
  double delta;
  double diffa;
  double diffb;
  double diffc;
  double ratio;
  int itest;
  int nfsav;
  int kopt;
  double fsave;
  double vquad;
  double fopt;
  double dsq;
  double xoptsq;
  int nptm;
  double alpha;
  double beta;
  double dnorm;
  int rc;
  double newf;
  int knew;
  int kbase;
  double denom;
  double delsq;
  double scaden;
  double biglsq;
  double distsq;
  double cauchy;
  double adelt;
  int prelim_nr;
  int rescue_nr;
  int altmov_nr;
  int trsbox_nr;
  int update_nr;
} bobyqa_data;
/*****************************************************************************/
 
/*****************************************************************************/
/* aic */
extern double aicSS(double ss, const int n, const int k);
extern int parFreeNr(const int n, double *pLower, double *pUpper);
extern int aicWeights(double *aic, double *w, int n);
extern double aicWeightedAvg(double *w, double *p, int n);
extern double aicModel(double *w, int n);
/*****************************************************************************/
 
/*****************************************************************************/
/* bobyqa */
extern bobyqa_result bobyqb(
  bobyqa_data *bdata
);
extern bobyqa_result bobyqa(
  int n,
  int npt,
  double *x,   
  const double *xl,
  const double *xu, 
  const double *dx,
  const double rhoend,
  double xtol_rel,
  double minf_max,
  double ftol_rel,
  double ftol_abs,
  int maxeval,
  int *nevals,
  double *minf,
  double (*f)(int n, double *x, void *objf_data),
  void *objf_data,
  double *working_space,
  int verbose
);
extern int bobyqa_minimize_single_parameter(
  bobyqa_data *bdata
);
extern char *bobyqa_rc(
  bobyqa_result rc
);
extern int fixed_params(
  int n, const double *lower, const double *upper, const double *delta
);
extern int bobyqa_working_memory_size(
  int n, int fitted_n, int npt, bobyqa_data *bdata
);
extern bobyqa_result bobyqa_set_memory(
  int n, int fitted_n, int npt, bobyqa_data *bdata, double *wm
);
extern bobyqa_result bobyqa_free_memory(
  bobyqa_data *bdata
);
extern bobyqa_result bobyqa_reset_memory(
  bobyqa_data *bdata
);
extern void bobyqa_print(
  bobyqa_data *bdata, int sw, FILE *fp
);
extern double bobyqa_x_funcval(
  bobyqa_data *bdata, double *x
);
extern void bobyqa_xfull(
  bobyqa_data *bdata
);
extern bobyqa_result bobyqa_set_optimization(
  int full_n,
  double *x,
  const double *dx,
  const double *xl,
  const double *xu,
  const double rhoend,
  double xtol_rel,
  double minf_max,
  double ftol_rel,
  double ftol_abs,
  int maxeval,
  double (*f)(int n, double *x, void *objf_data),
  void *objf_data,
  int verbose,
  bobyqa_data *bdata
);
/*****************************************************************************/
 
/*****************************************************************************/
/* bootstrap */
extern int bootstrap(
  int iterNr,
  double *cLim1, double *cLim2, double *SD, double *parameter,
  double *lowlim, double *uplim, int frameNr, double *origTac,
  double *fitTac, double *bsTAC, int parNr, double *weight,
  double (*objf)(int, double*, void*), char *status, int verbose
);
/*****************************************************************************/
 
/*****************************************************************************/
/* bvls */
extern int bvls(
  int key, const /*unsigned*/ int m, const /*unsigned*/ int n,
  double *a, double *b, double *bl, double *bu, double *x,
  double *w, double *act, double *zz, int *istate, int *iter,
  int verbose
);
extern int llsqWght(int N, int M, double **A, double *a, double *b, double *weight);
extern int llsqWghtSquared(int N, int M, double **A, double *a, double *b, double *weight);
/*****************************************************************************/
 
/*****************************************************************************/
/* constraints */
extern int modelCheckParameters(
  int par_nr, double *lower_p, double *upper_p, double *test_p,
  double *accept_p, double *penalty
);
extern int modelCheckLimits(
  int par_nr, double *lower_p, double *upper_p, double *test_p
);
extern int fitExpDecayNNLS(
  double *x, double *y, int n, double fittime, double kmin, double kmax,
  int pnr, double *a, double *k, int *fnr, int verbose
);
/*****************************************************************************/
 
/*****************************************************************************/
/* gaussdev */
 
extern long int GAUSSDEV_SEED;
 
extern unsigned int drandSeed(short int seed);
extern double gaussdev();
extern double gaussdev2();
extern void init_gaussdev();
extern double drand();
extern int rand_range(
  int nr, double *d, double low, double up, int type
);
/*****************************************************************************/
 
/*****************************************************************************/
/* hholder */
extern double householder_transform(double *v, int N);
extern int householder_hm(
  double tau, double *vector, double **matrix, int rowNr, int columnNr
);
extern int householder_hv(
  double tau, int size, double *v, double *w
);
extern double householder_norm(double *v, int size);
/*****************************************************************************/
 
/*****************************************************************************/
/* integr */
 
extern int interpolate(
  double *x, double *y, int nr,
  double *newx, double *newy, double *newyi, double *newyii, int newnr
);
extern int finterpolate(
  float *x, float *y, int nr,
  float *newx, float *newy, float *newyi, float *newyii, int newnr
);
extern int integrate(
  double *x, double *y, int nr,
  double *yi
);
extern int fintegrate(
  float *x, float *y, int nr,
  float *yi
);
extern int petintegrate(
  double *x1, double *x2, double *y, int nr, 
  double *newyi, double *newyii
);
extern int fpetintegrate(
  float *x1, float *x2, float *y, int nr, 
  float *newyi, float *newyii
);
extern int interpolate4pet(
  double *x, double *y, int nr,
  double *newx1, double *newx2, 
  double *newy, double *newyi, double *newyii, int newnr
);
extern int finterpolate4pet(
  float *x, float *y, int nr,
  float *newx1, float *newx2, 
  float *newy, float *newyi, float *newyii, int newnr
);
extern int petintegral(
  double *x1, double *x2, double *y, int nr,
  double *ie, double *iie
);
extern int fpetintegral(
  float *x1, float *x2, float *y, int nr,
  float *ie, float *iie
);
extern int petintegrate2fe(
  double *x1, double *x2, double *y, int nr,
  double *e, double *ie, double *iie
);
extern int fpetintegrate2fe(
  float *x1, float *x2, float *y, int nr,
  float *e, float *ie, float *iie
);
/*****************************************************************************/
 
/*****************************************************************************/
/* llsqwt */
 
extern int llsqwt(
  double *x, double *y, int n, double *wx, double *wy, double tol, double *w,
  double *ic, double *slope, double *nwss, double *sic, double *sslope,
  double *cx, double *cy
);
extern int best_llsqwt(
  double *x, double *y, double *wx, double *wy, int nr, int min_nr, int mode,
  double *slope, double *ic, double *nwss, double *sslope, double *sic,
  double *cx, double *cy, int *bnr
);
extern int llsqperp(double *x, double *y, int nr,
  double *slope, double *ic, double *ssd
);
extern int llsqperp3(double *x, double *y, int nr,
  double *slope, double *ic, double *ssd
);
extern int quadratic(double a, double b, double c, double *m1, double *m2);
extern int medianline(double *x, double *y, int nr, double *slope, double *ic);
/*****************************************************************************/
 
/*****************************************************************************/
/* lms */
extern double least_median_of_squares(
  double *data, int n
);
/*****************************************************************************/
 
/*****************************************************************************/
/* lts */
#ifndef CHI2INV_1
#define CHI2INV_1 0.45493642311957
#endif
 
extern int least_trimmed_square(
  double data[], long int n, double *mean, double *variance
);
/*****************************************************************************/
 
/*****************************************************************************/
/* median */
extern double d_kth_smallest(double *data, int n, int k);
extern double dmedian(double *data, int n);
extern double dmean(double *data, int n, double *sd);
extern double dmean_nan(double *data, int n, double *sd, int *vn);
/*****************************************************************************/
 
/*****************************************************************************/
/* mestim */
extern double mEstim(
  double *data, int nr, int iterNr, double cutoff
);
extern double huber(
  double x, double b
);
/*****************************************************************************/
 
/*****************************************************************************/
/* mtga */
 
#define MTGA_BEST_MIN_NR 5
 
typedef enum {
  PEARSON, 
  PERP,
  LLSQWT,
  MEDIAN
} linefit_method;
 
typedef enum {
  PRESET, EXCLUDE_BEGIN, EXCLUDE_END
} linefit_range;
 
extern int patlak_data(
  int data_nr, double *i, double *ii, double *c, double *x, double *y
);
extern int logan_data(
  int data_nr, double *i, double *ii, double *c, double *ci,
  double k2, double *x, double *y
);
extern int mtga_best_perp(
  double *x, double *y, int nr,
  double *slope, double *ic, double *ssd, int *fnr
);
/*****************************************************************************/
 
/*****************************************************************************/
/* nnls */
extern int nnls(
  double **a, int m, int n, double *b, double *x,
  double *rnorm, double *w, double *zz, int *index
);
extern int nnlsWght(
  int N, int M, double **A, double *b, double *weight
);
extern int nnlsWghtSquared(
  int N, int M, double **A, double *b, double *sweight
);
/*****************************************************************************/
 
/*****************************************************************************/
/* normaldistr */
extern double ndtr(double a);
extern double normal_pvalue_2(double x);
extern double normal_pvalue_1(double x);
/*****************************************************************************/
 
/*****************************************************************************/
/* pearson */
 
extern int pearson(
  double *x, double *y, int nr,
  double *k, double *kSD, double *b, double *bSD, double *r, double *ySD
);
extern int pearson2(
  double *x, double *y, char *is, int nr,
  double *k, double *kSD, double *b, double *bSD, double *r, double *ySD
);
extern int pearson3(
  double *x, double *y, int nr,
  double *k, double *kSD, double *b, double *bSD, double *r, double *ySD
);
extern int pearson4(
  double *x, double *y, int nr, double start, double end,
  double *k, double *kSD, double *b, double *bSD, double *r, double *ySD
);
extern int best_pearson(
  double *x, double *y, int nr, int min_nr, int *first, int *last,
  double *k, double *kSD, double *b, double *bSD, double *r, double *ySD
);
extern int mean(
  double *x, double *y, int nr,
  double *xmean, double *xsd, double *ymean, double *ysd
);
extern int regr_line(
  double *x, double *y, int n, double *m, double *c
);
extern int highest_slope(
  double *x, double *y, int n, int slope_n,
  double *m, double *c, double *xi, double *xh
);
extern int highest_slope_after(
  double *x, double *y, int n, int slope_n,
  double x_start, double *m, double *c, double *xi, double *xh
);
/*****************************************************************************/
 
/*****************************************************************************/
/* powell */
extern int POWELL_LINMIN_MAXIT;
 
extern int powell(double *p, double *delta, int parNr, double ftol, int *iterNr,
  double *fret, double (*_fun)(int, double*, void*), void *fundata, int verbose
);
/*****************************************************************************/
 
/*****************************************************************************/
/* qr */
extern int qrLSQ(
  double **mat, double *rhs, double *sol, 
  const unsigned int rows, const unsigned int cols, double *r2
);
extern int qr(
  double **A, int m, int n, double *B, double *X, double *rnorm,
  double *tau, double *res, double **wws, double *ws
);
extern int qr_decomp(
  double **a, int M, int N, double *tau, double **cchain, double *chain
);
extern int qr_solve(
  double **QR, int M, int N, double *tau, double *b, double *x,
  double *residual, double *resNorm, double **cchain, double *chain
);
extern int qr_weight(
  int N, int M, double **A, double *b, double *weight, double *ws
);
extern int qrLH(
  const unsigned int m, const unsigned int n,
  double *a, double *b, double *x, double *r2
);
/*****************************************************************************/
 
/*****************************************************************************/
/* runs_test */
extern int runs_test(
  double *data1, double *data2, int N, double alpha, double *p
);
extern int residuals(
  double *d1, double *d2, int Nr, int *Rnr, int *Nminus, int *Nplus
);
extern int mrl_between_tacs(
  double y1[], double y2[], int n
);
/*****************************************************************************/
 
/*****************************************************************************/
/* shuffle */
extern void random_shuffle(int *array, int n);
extern void randperm(int *array, int n, int a);
/*****************************************************************************/
 
/*****************************************************************************/
/* simplex */
 
extern double simplex(
  double (*_fun)(double*), int parNr, double *par, double *delta,
  double maxerr, int maxiter, int verbose
);
/*****************************************************************************/
 
/*****************************************************************************/
/* simulate */
/*****************************************************************************/
extern int simC3s(
  double *t, double *ca, int nr, double k1, double k2,
  double k3, double k4, double k5, double k6,
  double *ct, double *cta, double *ctb, double *ctc
);
extern int simC3p(
  double *t, double *ca, int nr, double k1, double k2,
  double k3, double k4, double k5, double k6,
  double *ct, double *cta, double *ctb, double *ctc
);
extern int simC3vs(
  double *t, double *ca, double *cb, int nr,
  double k1, double k2, double k3, double k4, double k5, double k6,
  double f, double vb, double fa,
  double *cpet, double *cta, double *ctb, double *ctc,
  double *ctab, double *ctvb
);
extern int simC3vp(
  double *t, double *ca, double *cb, int nr,
  double k1, double k2, double k3, double k4, double k5, double k6,
  double f, double vb, double fa,
  double *cpet, double *cta, double *ctb, double *ctc,
  double *ctab, double *ctvb
);
extern int simC2l(
  double *t, double *ca, int nr, double k1, double k2,
  double k3, double kLoss, double *ct, double *cta, double *ctb
);
extern int simC2vl(
  double *t, double *ca, double *cb, int nr,
  double k1, double k2, double k3, double kL,
  double f, double vb, double fa,
  double *cpet, double *cta, double *ctb,
  double *ctab, double *ctvb
);
extern int simC3vpKLoss(
  double *t, double *ca, double *cb, int nr,
  double k1, double k2, double k3, double k4, double k5, double k6,
  double kLoss, double f, double vb, double fa,
  double *cpet, double *cta, double *ctb, double *ctc,
  double *ctab, double *ctvb
);
extern int simRTCM(
  double *t, double *cr, int nr, double R1, double k2,
  double k3, double k4, double *ct, double *cta, double *ctb
);
extern int simSRTM(
  double *t, double *cr, int nr, double R1, double k2,
  double BP, double *ct
);
extern int simTRTM(
  double *t, double *cr, int nr, double R1, double k2,
  double k3, double *ct
);
extern int simHuangmet(
  double *t, double *ctot, int nr,
  double k01, double k12, double k21, double k03, double k34, double k43,
  double *c0, double *c1, double *c3 
);
extern int simTPCMOD0009c(
  double *t, double *ctot, int nr, double km, double k1m,
  double k2m, double k3m, double k4m, double *ca, double *cm
);
extern int simMBF(
  double *t, double *ci, int nr, double k1, double k2, double Vfit, double *ct
);
extern int simC1(
  double *t, double *ca, int nr, double k1, double k2, double *ct
);
extern int simC3DIvs(
  double *t, double *ca1, double *ca2, double *cb, int nr,
  double k1, double k2, double k3, double k4, double k5, double k6,
  double k1b, double k2b, double f, double vb, double fa,
  double *scpet, double *sct1, double *sct2, double *sct3, double *sct1b,
  double *sctab, double *sctvb
);
extern int simC4DIvp(
  double *t, double *ca1, double *ca2, double *cb, int nr,
  double k1, double k2, double k3, double k4, double k5, double k6, double k7,
  double km, double k1b, double k2b, double f, double vb, double fa,
  double *scpet, double *sct1, double *sct2, double *sct3,
  double *sct1b, double *sctab, double *sctvb,
  int verbose
);
extern int simC4DIvs(
  double *t, double *ca1, double *ca2, double *cb, int nr,
  double k1, double k2, double k3, double k4, double k5, double k6, double k7,
  double km, double k1b, double k2b, double f, double vb, double fa,
  double *scpet, double *sct1, double *sct2, double *sct3,
  double *sct1b, double *sctab, double *sctvb,
  int verbose
);
extern int simDispersion(
  double *x, double *y, int n,
  double tau1, double tau2, double *tmp
);
extern int simOxygen(
  double *t, double *ca1, double *ca2, double *ca1i, double *ca2i, const int n,
  const double k1a, const double k2a, const double km, 
  const double k1b, const double k2b, 
  const double vb, const double fa,
  double *scpet, double *sct1, double *sct2, double *sctab, 
  double *sctvb1, double *sctvb2, double *scvb1, double *scvb2,
  const int verbose
);
 
 
/* Deprecated functions. Please don't use these any more */
#define autointerpolateDFT dftAutointerpolate
#define c3sSIM simC3s
#define c3pSIM simC3p
#define c3vsSIM simC3vs
#define c3vpSIM simC3vp
#define c2lSIM simC2l
#define c2vlSIM simC2vl
#define rtcmSIM simRTCM
#define srtmSIM simSRTM
#define trtmSIM simTRTM
#define sim_dispersion simDispersion
/*****************************************************************************/
 
/*****************************************************************************/
/* o2 */
 
#define DEFAULT_SAO2 0.97
#define DEFAULT_P50HB 3.6
#define DEFAULT_P50MB 0.319
#define DEFAULT_NHB 2.7
#define DEFAULT_CHB 150.0
#define DEFAULT_CMB 4.7
 
extern double mo2k1k2(
  const double OER, const double SaO2, const double p50Hb, const double p50Mb,
  const double nHb, const double cHb, const double cMb,
  const int verbose
);
extern double mo2pO2(
  const double OER, const double K1k2, const double SaO2, const double p50Mb,
  const double cHb, const double cMb,
  const int verbose
);
/*****************************************************************************/
 
/*****************************************************************************/
/* tgo */
extern int TGO_SQUARED_TRANSF;
extern int TGO_LOCAL_INSIDE;
extern int TGO_LOCAL_OPT;
 
typedef struct {
  int topomin;
  double fvalue;
  double par[MAX_PARAMS];
  double delta[MAX_PARAMS]; // added by VO 2011-11-25
  double fvalrange; // added by VO 2011-11-25
} TGO_POINT;
 
extern int tgo(
  double *lowlim, double *uplim, double (*objf)(int, double*, void*),
  void *objfData, int dim, int neighNr, double *fmin, double *gmin,
  int samNr, int tgoNr, int verbose
);    
extern void tgoRandomParameters(
  TGO_POINT *p, int parNr, int sNr,
  double *low, double *up
);  
extern void tgoRandomParametersST(
  TGO_POINT *p, int parNr, int sNr,
  double *low, double *up
);
/*****************************************************************************/
 
/*****************************************************************************/
/* nlopt1d */
extern int nlopt1D(
  double (*_fun)(double, void*), void *_fundata,
  double x, double xl, double xu, double delta, double tol, const int maxeval,
  double *nx, double *nf, int verbose
);
/*****************************************************************************/
 
/*****************************************************************************/
#endif // _LIBTPCMODEL_H