text:
/* cdf/hypergeometric.c
 *
 * Copyright (C) 2004 Jason H. Stover.
 *
 * 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.
 */

/*
 * Computes the cumulative distribution function for a hypergeometric
 * random variable. A hypergeometric random variable X is the number
 * of elements of type 1 in a sample of size t, drawn from a population
 * of size n1 + n2, in which n1 are of type 1 and n2 are of type 2.
 *
 * This algorithm computes Pr( X <= k ) by summing the terms from
 * the mass function, Pr( X = k ).
 *
 * References:
 *
 * T. Wu. An accurate computation of the hypergeometric distribution
 * function. ACM Transactions on Mathematical Software. Volume 19, number 1,
 * March 1993.
 *  This algorithm is not used, since it requires factoring the
 *  numerator and denominator, then cancelling. It is more accurate
 *  than the algorithm used here, but the cancellation requires more
 *  time than the algorithm used here.
 *
 * W. Feller. An Introduction to Probability Theory and Its Applications,
 * third edition. 1968. Chapter 2, section 6.
 */

#include <config.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_randist.h>

#include "error.h"

static double
lower_tail (const unsigned int k, const unsigned int n1,
            const unsigned int n2, const unsigned int t)
{
  double relerr;
  int i = k;
  double s, P;

  s = gsl_ran_hypergeometric_pdf (i, n1, n2, t);
  P = s;

  while (i > 0)
    {
      double factor =
        (i / (n1 - i + 1.0)) * ((n2 + i - t) / (t - i + 1.0));
      s *= factor;
      P += s;
      relerr = s / P;
      if (relerr < GSL_DBL_EPSILON)
        break;
      i--;
    }

  return P;
}

static double
upper_tail (const unsigned int k, const unsigned int n1,
            const unsigned int n2, const unsigned int t)
{
  double relerr;
  unsigned int i = k + 1;
  double s, Q;

  s = gsl_ran_hypergeometric_pdf (i, n1, n2, t);
  Q = s;

  while (i < t)
    {
      double factor =
        ((n1 - i) / (i + 1.0)) * ((t - i) / (n2 + i + 1.0 - t));
      s *= factor;
      Q += s;
      relerr = s / Q;
      if (relerr < GSL_DBL_EPSILON)
        break;
      i++;
    }

  return Q;
}

/*
 * Pr (X <= k)
 */
double
gsl_cdf_hypergeometric_P (const unsigned int k,
                          const unsigned int n1,
                          const unsigned int n2, const unsigned int t)
{
  double P;

  if (t > (n1 + n2))
    {
      CDF_ERROR ("t larger than population size", GSL_EDOM);
    }
  else if (k >= n1 || k >= t)
    {
      P = 1.0;
    }
  else if (k < 0.0)
    {
      P = 0.0;
    }
  else
    {
      double midpoint = (int) (t * n1 / (n1 + n2));

      if (k >= midpoint)
        {
          P = 1 - upper_tail (k, n1, n2, t);
        }
      else
        {
          P = lower_tail (k, n1, n2, t);
        }
    }

  return P;
}

/*
 * Pr (X > k)
 */
double
gsl_cdf_hypergeometric_Q (const unsigned int k,
                          const unsigned int n1,
                          const unsigned int n2, const unsigned int t)
{
  double Q;

  if (t > (n1 + n2))
    {
      CDF_ERROR ("t larger than population size", GSL_EDOM);
    }
  else if (k >= n1 || k >= t)
    {
      Q = 0.0;
    }
  else if (k < 0.0)
    {
      Q = 1.0;
    }
  else
    {
      double midpoint = (int) (t * n1 / (n1 + n2));

      if (k < midpoint)
        {
          Q = 1 - lower_tail (k, n1, n2, t);
        }
      else
        {
          Q = upper_tail (k, n1, n2, t);
        }
    }

  return Q;
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
#include <assert.h>
#include <inttypes.h>
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

#include <gsl/gsl_math.h>

#include "imshift.h"
#include "conv1.h"

void showMatrix(double * I, size_t M, size_t N)
{
  for(size_t kk = 0; kk<N; kk++) {
    for(size_t ll = 0; ll<M; ll++)
    {
      printf("%5.2f ", I[ll*M+kk]);
    }
    printf("\n");
  }
}

int imshift3(double * V, size_t M, size_t N, size_t P,
    double dx, double dy, double dz, int method)
{

  int status = 0;

  double * Kx = NULL; size_t nKx;
  double * Ky = NULL; size_t nKy;
  double * Kz = NULL; size_t nKz;

  status = generateShift(&Kx, dx, &nKx, method);
  if(status != 0)
    return -1;

  status = generateShift(&Ky, dy, &nKy, method);
  if(status != 0) {
    free(Kx);
    return -2;
  }

  status = generateShift(&Kz, dz, &nKz, method);
  if(status != 0) {
    free(Kx);
    free(Ky);
    return -3;
  }

  status = conv1_3(V, M, N, P,
      Kx, nKx,
      Ky, nKy,
      Kz, nKz);

      free(Kz);
      free(Ky);
      free(Kx);

      return status;
}

int generateShift(double ** K, double delta, size_t * size, int order)
  /** Allocate space for a 1D shifting kernel, K,
   * shifted by delta in [-1,1]
   * returns the size of the kernel in size
   * order sets the polynomial order, 1 = linear, 2 = quadratic
   *
   * Errors:
   *  -1, K != NULL
   *  -2, no valid interpolation method selected
   *  -3, unable to allocate memory for the kernel
   *
   * */
{

  if(K[0]!=NULL)
    return -1;

  if(order == 1)
  {
    size_t N = 3;
    size[0] = N;
    K[0] = malloc(N*sizeof(double));
    if(K[0] == NULL)
      return(-3);

    for(size_t kk = 0; kk<N; kk++)
      K[0][kk] = GSL_MAX(0, 1-fabs(kk-(1-delta)));
    return 0;
  }

  if(order == 3)
  {
    size_t N = 5;
    size[0] = N;
    K[0] = malloc(N*sizeof(double));
    if(K[0] == NULL)
      return -3;

    // set up and solve polynomial equation here
    double a = -.5;

    for(size_t kk = 0; kk < N; kk++)
    {
      double x = kk-(2-delta);
      double ax = fabs(x);
      if(ax<=1)
        K[0][kk] = (a+2)*ax*ax*ax - (a+3)*ax*ax + 1;
      if(ax>1)
        K[0][kk] = a*ax*ax*ax -5*a*ax*ax  +8*a*ax -4*a;
      if(ax>2)
        K[0][kk] = 0;
    }
    return 0;
  }

  return -2;
}

int main(int argc, char ** argv)
{

  printf("Running %s\n", argv[0]);

  if(argc>1)
    printf("Warning: unused command line arguments\n");

  {
    printf("--> Testing kernels\n");
    double * K;
    double sum;

    for(double delta = -.5; delta < .6; delta+=.1) {
    for(int method = 1; method < 4; method = method +2)
    {
      K = NULL;
      size_t nK = 0;
      assert(generateShift(&K, delta, &nK, 3)==0);

      sum = 0;
      for(int kk = 0; kk<(int) nK; kk++)
        sum += K[kk];

      printf("delta: %f order: %d sum: %f\n", delta, method, sum);

     free(K);
    }
    }
  }

  {
    printf("--> Testing 1D\n");
    double delta = .3;
    size_t N = 10;
    double * V = malloc(N*sizeof(double));

    for(size_t kk=0; kk<N; kk++)
      V[kk] = kk%2;

    double * K = NULL;
    size_t nK = 0;
    assert(generateShift(&K, delta, &nK, 3)==0);

    printf("Kernel:\n");
    for(size_t kk = 0; kk<nK; kk++)
      printf("%.1f ", K[kk]);
    printf("\n");
    printf("Input:\n");
    for(size_t kk = 0; kk<N; kk++)
      printf("%.1f ", V[kk]);
    printf("\n");

    double * buff = calloc(N,sizeof(double));

    conv1(V, buff, N, K, nK, 1);
    printf("Output\n");
    for(size_t kk = 0; kk<N; kk++)
      printf("%.1f ", V[kk]);
    printf("\n");

    free(buff);
    free(K);
    free(V);
  }

  {
    printf("Testing 2D\n");
    size_t M = 11;
    size_t N = 11;

    double * Kx = NULL;
    double * Ky = NULL;
    size_t nKx;
    size_t nKy;

    double deltax = .1;
    double deltay = .6;

    int method = 1;
    assert(generateShift(&Kx, deltax, &nKx, method)==0);
    assert(generateShift(&Ky, deltay, &nKy, method)==0);

    printf("Kernels:\n");
double sum = 0;
    printf(" deltax: %f\n", deltax);
    for(size_t kk = 0; kk<nKx; kk++) {
      printf("%5.1f ", Kx[kk]);
      sum += Kx[kk];
    }
    printf("s: %.2f\n", sum);
    sum = 0;
    printf(" deltay: %f\n", deltay);
    for(size_t kk = 0; kk<nKy; kk++) {
      printf("%5.1f ", Ky[kk]);
      sum += Ky[kk];
    }
    printf("s: %.2f\n", sum);

    double * V = calloc(M*N, sizeof(double));
    V[(M*N-1)/2] = 1;

    printf("V:\n");
    showMatrix(V, M, N);
    double * C1 = malloc(M*N*sizeof(double));
    double * C2 = malloc(M*N*sizeof(double));
    memcpy(C1,V, M*N*sizeof(double));
    memcpy(C2,V, M*N*sizeof(double));

    conv1_3(C1,  M, N, 1, Kx, nKx, Ky, nKy, NULL, 0);
    conv1_3(C2,  M, N, 1, Ky, nKy, Kx, nKx, NULL, 0);

    printf("C1\n");
    showMatrix(C1, M, N);

    printf("C2\n");
    showMatrix(C2, M, N);

    printf("Testing symmetry, C1=C2' ... ");
    for(size_t kk = 0; kk<M; kk++)
      for(size_t ll = 0; ll<N; ll++)
        assert(C1[kk+M*ll] == C2[ll+M*kk]);
    printf("ok!\n");

    free(V);
    free(Ky);
    free(Kx);
    free(C1);
    free(C2);

  }

  {
    printf("Testing 3D\n");
    double delta = .3;
    size_t M = 11;
    size_t N = 11;
    size_t P = 11;
    double * V3 = calloc(M*N*P,sizeof(double));
    V3[(M*N*P-1)/2] = 1;
    double * K = NULL;
    size_t nK = 0;
    assert(generateShift(&K, delta, &nK, 1)==0);

    clock_t t1 = clock();
    conv1_3(V3, M, N,P, K,nK,K,nK,K,nK);
    clock_t t2 = clock();

    printf("V3(:,:,6):\n");
    showMatrix(V3 + M*N*6, M, N);

    printf("Time: %f ms\n", (double) 1000*(t2-t1)/CLOCKS_PER_SEC);

    free(V3);
    free(K);
  }

  return 0;
}

ext: c

max_stars_repo_name: elgw/dotter

max_issues_repo_name: elgw/dotter

max_forks_repo_name: elgw/dotter


text:
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_histogram.h>

int
main (int argc, char **argv)
{
  double a, b;
  size_t n;

  if (argc != 4)
    {
      printf ("Usage: gsl-histogram xmin xmax n\n"
              "Computes a histogram of the data "
              "on stdin using n bins from xmin "
              "to xmax\n");
      exit (0);
    }

  a = atof (argv[1]);
  b = atof (argv[2]);
  n = atoi (argv[3]);

  {
    double x;
    gsl_histogram * h = gsl_histogram_alloc (n);
    gsl_histogram_set_ranges_uniform (h, a, b);

    while (fscanf (stdin, "%lg", &x) == 1)
      {
        gsl_histogram_increment (h, x);
      }
    gsl_histogram_fprintf (stdout, h, "%g", "%g");
    gsl_histogram_free (h);
  }
  exit (0);
}

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
//Causal FIR filtering of each row or col of X according to dim.

//FIR impulse responses are given in matrix B.
//For dim=0, X is NxT and B is NxL
//For dim=1, X is TxN and B is LxN
//where N is the number of neurons, T the number of time points, and L the FIR filter order.

#include <stdio.h>
#include <cblas.h>

#ifdef __cplusplus
namespace codee {
extern "C" {
#endif

int fir_s (float *Y, const float *X, const float *B, const size_t N, const size_t T, const size_t L, const char iscolmajor, const size_t dim);
int fir_d (double *Y, const double *X, const double *B, const size_t N, const size_t T, const size_t L, const char iscolmajor, const size_t dim);
int fir_c (float *Y, const float *X, const float *B, const size_t N, const size_t T, const size_t L, const char iscolmajor, const size_t dim);
int fir_z (double *Y, const double *X, const double *B, const size_t N, const size_t T, const size_t L, const char iscolmajor, const size_t dim);

int fir_s (float *Y, const float *X, const float *B, const size_t N, const size_t T, const size_t L, const char iscolmajor, const size_t dim)
{
    const float z = 0.0f;

    //Initialize Y to 0
    cblas_scopy((int)(N*T),&z,0,Y,1);

    if (N==1u)
    {
        for (size_t l=0u; l<L; l++) { cblas_saxpy((int)(T-l),B[l],&X[0],1,&Y[l],1); }
    }
    else if (dim==0u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                for (size_t l=0u; l<L; l++) { cblas_saxpy((int)(T-l),B[n+l*N],&X[n],(int)N,&Y[n+l*N],(int)N); }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                for (size_t l=0u; l<L; l++) { cblas_saxpy((int)(T-l),B[n*L+l],&X[n*T],1,&Y[n*T+l],1); }
            }
        }
    }
    else if (dim==1u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                for (size_t l=0u; l<L; l++) { cblas_saxpy((int)(T-l),B[n*L+l],&X[n*T],1,&Y[n*T+l],1); }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                for (size_t l=0u; l<L; l++) { cblas_saxpy((int)(T-l),B[n+l*N],&X[n],(int)N,&Y[n+l*N],(int)N); }
            }
        }
    }
    else
    {
        fprintf(stderr,"error in fir_s: dim must be 0 or 1.\n"); return 1;
    }

    return 0;
}

int fir_d (double *Y, const double *X, const double *B, const size_t N, const size_t T, const size_t L, const char iscolmajor, const size_t dim)
{
    const double z = 0.0;

    //Initialize Y to 0
    cblas_dcopy((int)(N*T),&z,0,Y,1);

    if (N==1u)
    {
        for (size_t l=0u; l<L; l++) { cblas_daxpy((int)(T-l),B[l],&X[0],1,&Y[l],1); }
    }
    else if (dim==0u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                for (size_t l=0u; l<L; l++) { cblas_daxpy((int)(T-l),B[n+l*N],&X[n],(int)N,&Y[n+l*N],(int)N); }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                for (size_t l=0u; l<L; l++) { cblas_daxpy((int)(T-l),B[n*L+l],&X[n*T],1,&Y[n*T+l],1); }
            }
        }
    }
    else if (dim==1u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                for (size_t l=0u; l<L; l++) { cblas_daxpy((int)(T-l),B[n*L+l],&X[n*T],1,&Y[n*T+l],1); }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                for (size_t l=0u; l<L; l++) { cblas_daxpy((int)(T-l),B[n+l*N],&X[n],(int)N,&Y[n+l*N],(int)N); }
            }
        }
    }
    else
    {
        fprintf(stderr,"error in fir_d: dim must be 0 or 1.\n"); return 1;
    }

    return 0;
}

int fir_c (float *Y, const float *X, const float *B, const size_t N, const size_t T, const size_t L, const char iscolmajor, const size_t dim)
{
    const float z[2] = {0.0f,0.0f};

    //Initialize Y to 0
    cblas_ccopy((int)(N*T),z,0,Y,1);

    if (N==1u)
    {
        for (size_t l=0u; l<L; l++) { cblas_caxpy((int)(T-l),&B[2*l],&X[0],1,&Y[2*l],1); }
    }
    else if (dim==0u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                for (size_t l=0u; l<L; l++) { cblas_caxpy((int)(T-l),&B[2*(n+l*N)],&X[2*n],(int)N,&Y[2*(n+l*N)],(int)N); }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                for (size_t l=0u; l<L; l++) { cblas_caxpy((int)(T-l),&B[2*(n*L+l)],&X[2*n*T],1,&Y[2*(n*T+l)],1); }
            }
        }
    }
    else if (dim==1u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                for (size_t l=0u; l<L; l++) { cblas_caxpy((int)(T-l),&B[2*(n*L+l)],&X[2*n*T],1,&Y[2*(n*T+l)],1); }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                for (size_t l=0u; l<L; l++) { cblas_caxpy((int)(T-l),&B[2*(n+l*N)],&X[2*n],(int)N,&Y[2*(n+l*N)],(int)N); }
            }
        }
    }
    else
    {
        fprintf(stderr,"error in fir_c: dim must be 0 or 1.\n"); return 1;
    }

    return 0;
}

int fir_z (double *Y, const double *X, const double *B, const size_t N, const size_t T, const size_t L, const char iscolmajor, const size_t dim)
{
    const double z[2] = {0.0,0.0};

    //Initialize Y to 0
    cblas_zcopy((int)(N*T),z,0,Y,1);

    if (N==1u)
    {
        for (size_t l=0u; l<L; l++) { cblas_zaxpy((int)(T-l),&B[2*l],&X[0],1,&Y[2*l],1); }
    }
    else if (dim==0u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                for (size_t l=0u; l<L; l++) { cblas_zaxpy((int)(T-l),&B[2*(n+l*N)],&X[2*n],(int)N,&Y[2*(n+l*N)],(int)N); }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                for (size_t l=0u; l<L; l++) { cblas_zaxpy((int)(T-l),&B[2*(n*L+l)],&X[2*n*T],1,&Y[2*(n*T+l)],1); }
            }
        }
    }
    else if (dim==1u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                for (size_t l=0u; l<L; l++) { cblas_zaxpy((int)(T-l),&B[2*(n*L+l)],&X[2*n*T],1,&Y[2*(n*T+l)],1); }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                for (size_t l=0u; l<L; l++) { cblas_zaxpy((int)(T-l),&B[2*(n+l*N)],&X[2*n],(int)N,&Y[2*(n+l*N)],(int)N); }
            }
        }
    }
    else
    {
        fprintf(stderr,"error in fir_z: dim must be 0 or 1.\n"); return 1;
    }

    return 0;
}

#ifdef __cplusplus
}
}
#endif


text:
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <assert.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_sort.h>

#include "cosmocalc.h"

gsl_spline *cosmocalc_aexpn2age_spline = NULL;
gsl_interp_accel *cosmocalc_aexpn2age_acc = NULL;
double AGE = 0;

/* function for integration using gsl integration */
static double age_integ_funct(double a, void *p)
{
  return 1.0/a/hubble_noscale(a);
}

/* init function  - some help from Gadget-2 applied here */
void init_cosmocalc_age_table(void)
{
#define WORKSPACE_NUM 100000
#define ABSERR 0.0
#define RELERR 1e-8

  static int initFlag = 1;
  static int currCosmoNum;

  gsl_integration_workspace *workspace;
  gsl_function F;
  long i;
  double result,abserr,afact;
  double age_table[COSMOCALC_AGE_TABLE_LENGTH];
  double aexpn_table[COSMOCALC_AGE_TABLE_LENGTH];

  if(initFlag == 1 || currCosmoNum != cosmoData.cosmoNum)
    {
      initFlag = 0;
      currCosmoNum = cosmoData.cosmoNum;

      workspace = gsl_integration_workspace_alloc((size_t) WORKSPACE_NUM);

      aexpn_table[0] = 0.0;
      age_table[0] = 0.0;
      for(i=1;i<COSMOCALC_AGE_TABLE_LENGTH;++i)
	{
	  afact = (1.0 - 0.0)/(COSMOCALC_AGE_TABLE_LENGTH-1.0)*((double) i);
	  F.function = &age_integ_funct;
	  gsl_integration_qag(&F,1e-12,afact,ABSERR,RELERR,(size_t) WORKSPACE_NUM,GSL_INTEG_GAUSS51,workspace,&result,&abserr);
	  aexpn_table[i] = afact;
	  age_table[i] = result*TH/cosmoData.h;
	}

      gsl_integration_workspace_free(workspace);
#undef ABSERR
#undef RELERR
#undef WORKSPACE_NUM

      //init the spline and accelerators
      if(cosmocalc_aexpn2age_spline != NULL)
	gsl_spline_free(cosmocalc_aexpn2age_spline);
      cosmocalc_aexpn2age_spline = gsl_spline_alloc(GSL_SPLINE_TYPE,(size_t) (COSMOCALC_AGE_TABLE_LENGTH));
      gsl_spline_init(cosmocalc_aexpn2age_spline,aexpn_table,age_table,(size_t) (COSMOCALC_AGE_TABLE_LENGTH));
      if(cosmocalc_aexpn2age_acc != NULL)
	gsl_interp_accel_reset(cosmocalc_aexpn2age_acc);
      else
	cosmocalc_aexpn2age_acc = gsl_interp_accel_alloc();

      AGE = gsl_spline_eval(cosmocalc_aexpn2age_spline,1.0,cosmocalc_aexpn2age_acc);
    }
}

double age(double a)
{
  static int initFlag = 1;
  static int currCosmoNum;

  if(initFlag == 1 || currCosmoNum != cosmoData.cosmoNum)
    {
      initFlag = 0;
      currCosmoNum = cosmoData.cosmoNum;
      init_cosmocalc_age_table();
    }

  return gsl_spline_eval(cosmocalc_aexpn2age_spline,a,cosmocalc_aexpn2age_acc);
}

//call age first to set AGE
double lookback(double a)
{
  double age_a = age(a);
  return AGE - age_a;
}

max_stars_repo_name: beckermr/cosmocalc

max_issues_repo_name: beckermr/cosmocalc

max_forks_repo_name: beckermr/cosmocalc


text:
/* eigen/gsl_eigen.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2006, 2007 Gerard Jungman, Brian Gough, Patrick Alken
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#ifndef __GSL_EIGEN_H__
#define __GSL_EIGEN_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 <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>

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

__BEGIN_DECLS

typedef struct {
  size_t size;
  double * d;
  double * sd;
} gsl_eigen_symm_workspace;

GSL_FUN gsl_eigen_symm_workspace * gsl_eigen_symm_alloc (const size_t n);
GSL_FUN void gsl_eigen_symm_free (gsl_eigen_symm_workspace * w);
GSL_FUN int gsl_eigen_symm (gsl_matrix * A, gsl_vector * eval, gsl_eigen_symm_workspace * w);

typedef struct {
  size_t size;
  double * d;
  double * sd;
  double * gc;
  double * gs;
} gsl_eigen_symmv_workspace;

GSL_FUN gsl_eigen_symmv_workspace * gsl_eigen_symmv_alloc (const size_t n);
GSL_FUN void gsl_eigen_symmv_free (gsl_eigen_symmv_workspace * w);
GSL_FUN int gsl_eigen_symmv (gsl_matrix * A, gsl_vector * eval, gsl_matrix * evec, gsl_eigen_symmv_workspace * w);

typedef struct {
  size_t size;
  double * d;
  double * sd;
  double * tau;
} gsl_eigen_herm_workspace;

GSL_FUN gsl_eigen_herm_workspace * gsl_eigen_herm_alloc (const size_t n);
GSL_FUN void gsl_eigen_herm_free (gsl_eigen_herm_workspace * w);
GSL_FUN int gsl_eigen_herm (gsl_matrix_complex * A, gsl_vector * eval,
                         gsl_eigen_herm_workspace * w);

typedef struct {
  size_t size;
  double * d;
  double * sd;
  double * tau;
  double * gc;
  double * gs;
} gsl_eigen_hermv_workspace;

GSL_FUN gsl_eigen_hermv_workspace * gsl_eigen_hermv_alloc (const size_t n);
GSL_FUN void gsl_eigen_hermv_free (gsl_eigen_hermv_workspace * w);
GSL_FUN int gsl_eigen_hermv (gsl_matrix_complex * A, gsl_vector * eval,
                           gsl_matrix_complex * evec,
                           gsl_eigen_hermv_workspace * w);

typedef struct {
  size_t size;           /* matrix size */
  size_t max_iterations; /* max iterations since last eigenvalue found */
  size_t n_iter;         /* number of iterations since last eigenvalue found */
  size_t n_evals;        /* number of eigenvalues found so far */

  int compute_t;         /* compute Schur form T = Z^t A Z */

  gsl_matrix *H;         /* pointer to Hessenberg matrix */
  gsl_matrix *Z;         /* pointer to Schur vector matrix */
} gsl_eigen_francis_workspace;

GSL_FUN gsl_eigen_francis_workspace * gsl_eigen_francis_alloc (void);
GSL_FUN void gsl_eigen_francis_free (gsl_eigen_francis_workspace * w);
GSL_FUN void gsl_eigen_francis_T (const int compute_t,
                          gsl_eigen_francis_workspace * w);
GSL_FUN int gsl_eigen_francis (gsl_matrix * H, gsl_vector_complex * eval,
                       gsl_eigen_francis_workspace * w);
GSL_FUN int gsl_eigen_francis_Z (gsl_matrix * H, gsl_vector_complex * eval,
                         gsl_matrix * Z,
                         gsl_eigen_francis_workspace * w);

typedef struct {
  size_t size;                 /* size of matrices */
  gsl_vector *diag;            /* diagonal matrix elements from balancing */
  gsl_vector *tau;             /* Householder coefficients */
  gsl_matrix *Z;               /* pointer to Z matrix */
  int do_balance;              /* perform balancing transformation? */
  size_t n_evals;              /* number of eigenvalues found */

  gsl_eigen_francis_workspace *francis_workspace_p;
} gsl_eigen_nonsymm_workspace;

GSL_FUN gsl_eigen_nonsymm_workspace * gsl_eigen_nonsymm_alloc (const size_t n);
GSL_FUN void gsl_eigen_nonsymm_free (gsl_eigen_nonsymm_workspace * w);
GSL_FUN void gsl_eigen_nonsymm_params (const int compute_t, const int balance,
                               gsl_eigen_nonsymm_workspace *w);
GSL_FUN int gsl_eigen_nonsymm (gsl_matrix * A, gsl_vector_complex * eval,
                       gsl_eigen_nonsymm_workspace * w);
GSL_FUN int gsl_eigen_nonsymm_Z (gsl_matrix * A, gsl_vector_complex * eval,
                         gsl_matrix * Z, gsl_eigen_nonsymm_workspace * w);

typedef struct {
  size_t size;                 /* size of matrices */
  gsl_vector *work;            /* scratch workspace */
  gsl_vector *work2;           /* scratch workspace */
  gsl_vector *work3;           /* scratch workspace */

  gsl_matrix *Z;               /* pointer to Schur vectors */

  gsl_eigen_nonsymm_workspace *nonsymm_workspace_p;
} gsl_eigen_nonsymmv_workspace;

GSL_FUN gsl_eigen_nonsymmv_workspace * gsl_eigen_nonsymmv_alloc (const size_t n);
GSL_FUN void gsl_eigen_nonsymmv_free (gsl_eigen_nonsymmv_workspace * w);
GSL_FUN void gsl_eigen_nonsymmv_params (const int balance,
                                gsl_eigen_nonsymmv_workspace *w);
GSL_FUN int gsl_eigen_nonsymmv (gsl_matrix * A, gsl_vector_complex * eval,
                        gsl_matrix_complex * evec,
                        gsl_eigen_nonsymmv_workspace * w);
GSL_FUN int gsl_eigen_nonsymmv_Z (gsl_matrix * A, gsl_vector_complex * eval,
                          gsl_matrix_complex * evec, gsl_matrix * Z,
                          gsl_eigen_nonsymmv_workspace * w);

typedef struct {
  size_t size;            /* size of matrices */
  gsl_eigen_symm_workspace *symm_workspace_p;
} gsl_eigen_gensymm_workspace;

GSL_FUN gsl_eigen_gensymm_workspace * gsl_eigen_gensymm_alloc (const size_t n);
GSL_FUN void gsl_eigen_gensymm_free (gsl_eigen_gensymm_workspace * w);
GSL_FUN int gsl_eigen_gensymm (gsl_matrix * A, gsl_matrix * B,
                       gsl_vector * eval, gsl_eigen_gensymm_workspace * w);
GSL_FUN int gsl_eigen_gensymm_standardize (gsl_matrix * A, const gsl_matrix * B);

typedef struct {
  size_t size;            /* size of matrices */
  gsl_eigen_symmv_workspace *symmv_workspace_p;
} gsl_eigen_gensymmv_workspace;

GSL_FUN gsl_eigen_gensymmv_workspace * gsl_eigen_gensymmv_alloc (const size_t n);
GSL_FUN void gsl_eigen_gensymmv_free (gsl_eigen_gensymmv_workspace * w);
GSL_FUN int gsl_eigen_gensymmv (gsl_matrix * A, gsl_matrix * B,
                        gsl_vector * eval, gsl_matrix * evec,
                        gsl_eigen_gensymmv_workspace * w);

typedef struct {
  size_t size;            /* size of matrices */
  gsl_eigen_herm_workspace *herm_workspace_p;
} gsl_eigen_genherm_workspace;

GSL_FUN gsl_eigen_genherm_workspace * gsl_eigen_genherm_alloc (const size_t n);
GSL_FUN void gsl_eigen_genherm_free (gsl_eigen_genherm_workspace * w);
GSL_FUN int gsl_eigen_genherm (gsl_matrix_complex * A, gsl_matrix_complex * B,
                       gsl_vector * eval, gsl_eigen_genherm_workspace * w);
GSL_FUN int gsl_eigen_genherm_standardize (gsl_matrix_complex * A,
                                   const gsl_matrix_complex * B);

typedef struct {
  size_t size;            /* size of matrices */
  gsl_eigen_hermv_workspace *hermv_workspace_p;
} gsl_eigen_genhermv_workspace;

GSL_FUN gsl_eigen_genhermv_workspace * gsl_eigen_genhermv_alloc (const size_t n);
GSL_FUN void gsl_eigen_genhermv_free (gsl_eigen_genhermv_workspace * w);
GSL_FUN int gsl_eigen_genhermv (gsl_matrix_complex * A, gsl_matrix_complex * B,
                        gsl_vector * eval, gsl_matrix_complex * evec,
                        gsl_eigen_genhermv_workspace * w);

typedef struct {
  size_t size;            /* size of matrices */
  gsl_vector *work;       /* scratch workspace */

  size_t n_evals;         /* number of eigenvalues found */
  size_t max_iterations;  /* maximum QZ iterations allowed */
  size_t n_iter;          /* number of iterations since last eigenvalue found */
  double eshift;          /* exceptional shift counter */

  int needtop;            /* need to compute top index? */

  double atol;            /* tolerance for splitting A matrix */
  double btol;            /* tolerance for splitting B matrix */

  double ascale;          /* scaling factor for shifts */
  double bscale;          /* scaling factor for shifts */

  gsl_matrix *H;          /* pointer to hessenberg matrix */
  gsl_matrix *R;          /* pointer to upper triangular matrix */

  int compute_s;          /* compute generalized Schur form S */
  int compute_t;          /* compute generalized Schur form T */

  gsl_matrix *Q;          /* pointer to left Schur vectors */
  gsl_matrix *Z;          /* pointer to right Schur vectors */
} gsl_eigen_gen_workspace;

GSL_FUN gsl_eigen_gen_workspace * gsl_eigen_gen_alloc (const size_t n);
GSL_FUN void gsl_eigen_gen_free (gsl_eigen_gen_workspace * w);
GSL_FUN void gsl_eigen_gen_params (const int compute_s, const int compute_t,
                           const int balance, gsl_eigen_gen_workspace * w);
GSL_FUN int gsl_eigen_gen (gsl_matrix * A, gsl_matrix * B,
                   gsl_vector_complex * alpha, gsl_vector * beta,
                   gsl_eigen_gen_workspace * w);
GSL_FUN int gsl_eigen_gen_QZ (gsl_matrix * A, gsl_matrix * B,
                      gsl_vector_complex * alpha, gsl_vector * beta,
                      gsl_matrix * Q, gsl_matrix * Z,
                      gsl_eigen_gen_workspace * w);

typedef struct {
  size_t size;            /* size of matrices */

  gsl_vector *work1;      /* 1-norm of columns of A */
  gsl_vector *work2;      /* 1-norm of columns of B */
  gsl_vector *work3;      /* real part of eigenvector */
  gsl_vector *work4;      /* imag part of eigenvector */
  gsl_vector *work5;      /* real part of back-transformed eigenvector */
  gsl_vector *work6;      /* imag part of back-transformed eigenvector */

  gsl_matrix *Q;          /* pointer to left Schur vectors */
  gsl_matrix *Z;          /* pointer to right Schur vectors */

  gsl_eigen_gen_workspace *gen_workspace_p;
} gsl_eigen_genv_workspace;

GSL_FUN gsl_eigen_genv_workspace * gsl_eigen_genv_alloc (const size_t n);
GSL_FUN void gsl_eigen_genv_free (gsl_eigen_genv_workspace * w);
GSL_FUN int gsl_eigen_genv (gsl_matrix * A, gsl_matrix * B,
                    gsl_vector_complex * alpha, gsl_vector * beta,
                    gsl_matrix_complex * evec,
                    gsl_eigen_genv_workspace * w);
GSL_FUN int gsl_eigen_genv_QZ (gsl_matrix * A, gsl_matrix * B,
                       gsl_vector_complex * alpha, gsl_vector * beta,
                       gsl_matrix_complex * evec,
                       gsl_matrix * Q, gsl_matrix * Z,
                       gsl_eigen_genv_workspace * w);

typedef enum {
  GSL_EIGEN_SORT_VAL_ASC,
  GSL_EIGEN_SORT_VAL_DESC,
  GSL_EIGEN_SORT_ABS_ASC,
  GSL_EIGEN_SORT_ABS_DESC
}
gsl_eigen_sort_t;

/* Sort eigensystem results based on eigenvalues.
 * Sorts in order of increasing value or increasing
 * absolute value.
 *
 * exceptions: GSL_EBADLEN
 */

GSL_FUN int gsl_eigen_symmv_sort(gsl_vector * eval, gsl_matrix * evec,
                         gsl_eigen_sort_t sort_type);

GSL_FUN int gsl_eigen_hermv_sort(gsl_vector * eval, gsl_matrix_complex * evec,
                         gsl_eigen_sort_t sort_type);

GSL_FUN int gsl_eigen_nonsymmv_sort(gsl_vector_complex * eval,
                            gsl_matrix_complex * evec,
                            gsl_eigen_sort_t sort_type);

GSL_FUN int gsl_eigen_gensymmv_sort (gsl_vector * eval, gsl_matrix * evec,
                             gsl_eigen_sort_t sort_type);

GSL_FUN int gsl_eigen_genhermv_sort (gsl_vector * eval, gsl_matrix_complex * evec,
                             gsl_eigen_sort_t sort_type);

GSL_FUN int gsl_eigen_genv_sort (gsl_vector_complex * alpha, gsl_vector * beta,
                         gsl_matrix_complex * evec,
                         gsl_eigen_sort_t sort_type);

/* Prototypes for the schur module */

GSL_FUN int gsl_schur_gen_eigvals(const gsl_matrix *A, const gsl_matrix *B,
                          double *wr1, double *wr2, double *wi,
                          double *scale1, double *scale2);

GSL_FUN int gsl_schur_solve_equation(double ca, const gsl_matrix *A, double z,
                             double d1, double d2, const gsl_vector *b,
                             gsl_vector *x, double *s, double *xnorm,
                             double smin);

GSL_FUN int gsl_schur_solve_equation_z(double ca, const gsl_matrix *A,
                               gsl_complex *z, double d1, double d2,
                               const gsl_vector_complex *b,
                               gsl_vector_complex *x, double *s,
                               double *xnorm, double smin);

/* The following functions are obsolete: */

/* Eigensolve by Jacobi Method
 *
 * The data in the matrix input is destroyed.
 *
 * exceptions:
 */
GSL_FUN int
gsl_eigen_jacobi(gsl_matrix * matrix,
                      gsl_vector * eval,
                      gsl_matrix * evec,
                      unsigned int max_rot,
                      unsigned int * nrot);

/* Invert by Jacobi Method
 *
 * exceptions:
 */
GSL_FUN int
gsl_eigen_invert_jacobi(const gsl_matrix * matrix,
                             gsl_matrix * ainv,
                             unsigned int max_rot);

__END_DECLS

#endif /* __GSL_EIGEN_H__ */

ext: h

max_stars_repo_name: mgreter/astrometrylib

max_issues_repo_name: zzpwahaha/VimbaCamJILA

max_forks_repo_name: zzpwahaha/VimbaCamJILA


text:
#ifndef NumericalApproximation_h
#define NumericalApproximation_h
/** This is free and unencumbered software released into the public domain.
The authors of ISIS do not claim copyright on the contents of this file.
For more details about the LICENSE terms and the AUTHORS, you will
find files of those names at the top level of this repository. **/

/* SPDX-License-Identifier: CC0-1.0 */

#include <string>
#include <vector>
#include <map>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include "IException.h"

using namespace std;

namespace Isis {
  /**
   * @brief NumericalApproximation provides various numerical
   *        analysis methods of interpolation, extrapolation and
   *        approximation of a tabulated set of @a x, @a y data.
   *
   * This class contains a merged version of the Isis classes
   * @b DataInterp and @b NumericalMethods.  In addition,
   * some methods from @b AtmosModel were moved to this class, the
   * @a CubicNeighborhood interpolation type was
   * adapted from the IDL routine called @b interpol (using the
   * "/spline" keyword), and several differentiation and
   * integration approximation methods were created.
   *
   * NumericalApproximation contains 9 types of data interpolation
   * routines.  These routines act on x, y pairs of data points.
   * The following forms of data interpolation are supported:
   *
   *

   *

   *     Name
   *     Type (enum)
   *     MinPoints
   *     Description
   *

   *

   *     Linear
   *     @a Linear
   *     2
   *
   *        Linear interpolation approximates a curve by
   *        concatinating line segments between known data
   *        points.  This results in a continuous curve with a
   *        discontinuities of the derivative at the known data
   *        points.  This interpolation type uses the GSL
   *        library routines.
   *
   *
   *

   *

   *     Polynomial
   *     @a Polynomial
   *     3
   *
   *        Computes a polynomial interpolation.  This method
   *        should only be used for interpolating small numbers
   *        of points because polynomial interpolation introduces
   *        large oscillations, even for well-behaved datasets.
   *        The number of terms in the interpolating polynomial
   *        is equal to the number of points.  This interpolation
   *        type uses the GSL library routines.
   *
   *

   *

   *     Neville's Algorithm for Polynomial Interpolation
   *     @a PolynomialNeville
   *     3
   *
   *     Polynomial interpolation using Neville's algorithm. This
   *     first fits a polynomial of degree 0 through each point.
   *     On the second iteration, the polynomial of adjoining
   *     indices are combined to fit through pairs of points.
   *     This process repeats until a pyramid of approximations is
   *     reached. This is based on the Newton form of the
   *     interpolating polynomial and the recursion relation of
   *     the divided differences.
   *
   *

   *

   *     Natural Cubic Spline
   *     @a CubicNatural
   *     3
   *
   *       Cubic spline with natural boundary conditions.  The
   *       resulting curve is piecewise cubic on each interval
   *       with matching first and second derivatives at the
   *       supplied data points.  The second derivative is chosen
   *       to be zero at the endpoints, i.e. the first and last
   *       points of the data set. This is also refered to as
   *       free boundary condtions.  This interpolation type uses
   *       the GSL
   *        library routines.
   *
   *

   *

   *     Clamped Cubic Spline
   *     @a CubicClamped
   *     3
   *
   *     Cubic Spline interpolation with clamped boundary
   *     conditions.  The resulting curve is piecewise cubic on
   *     each interval.  For this type of boundary condition to
   *     hold, it is necessary to have either the values of the
   *     derivative at the endpoints or an accurate approximation
   *     to those values. In general, clamped boundary conditions
   *     lead to more accurate approximations since they include
   *     more information about the function.
   *
   *

   *

   *     Periodic Natural Cubic Spline
   *     @a CubicNatPeriodic
   *     2
   *
   *       Cubic spline with periodic boundary conditions.  The
   *       resulting curve is piecewise cubic on each interval
   *       with matching first and second derivatives at the
   *       supplied data points.  The derivatives at the first
   *       and last points are also matched.  Note that the last
   *       point in the data must have the same y-value as the
   *       first point, otherwise the resulting periodic
   *       interpolation will have a discontinuity at the
   *       boundary.  This interpolation type uses the GSL
   *       library routines.
   *
   *

   *

   *     Cubic Spline about 4-point Neighborhood
   *     @a CubicNeighborhood
   *     4
   *
   *       Cubic Spline interpolation using 4-pt Neighborhoods
   *       with natural boundary conditions. This interpolation
   *       type is adapted from the IDL interpol.pro application
   *       using the "/spline" keyword on an irregular grid.
   *       This type of cubic spline fits a natural cubic spline
   *       to the 4-point neighborhood of known data points
   *       surrounding the @a x value at which we wish to
   *       evaluate. For example, suppose {@a x0,
   *       @a x1, ..., @a xn} is the
   *       array of known domain values in the data set and @f$
   *       x_i \leq u < x_{i+1} @f$ for some @e i such that
   *       @f$ 0 \leq i \leq n @f$, then to evaluate the
   *       @a y value at @e u, this method of interpolation
   *       evaluates the ordinary cubic spline consisting of the
   *       data set {@a xi-1, @a xi,
   *       @a xi+1, @a xi+2} at
   *       @e u.
   *
   *

   *

   *     Cubic Spline using Hermite cubic polynomial
   *     @a CubicHermite
   *     2
   *
   *       Cubic Spline interpolation using the Hermite cubic
   *       polynomial (also called the unique Hermite
   *       interpolating fundamental polynomial of degree 3).
   *       This method requires the input of the slopes of the
   *       tangent lines of the curve (i.e. velocities, or first
   *       derivatives) at each of the data points.  This ensures
   *       that the approximation has the same "shape" as the
   *       curve at the known data points. This interpolation type
   *       uses piecewise Hermite cubic polynomials for each
   *       interval, (@a x0, @a x1), using
   *       the known data points, (@a
   *       x0,@a f(@a x0)) and
   *       (@a x1,@a f(@a x1)), and their
   *       derivatives,
   *       @a f '(@a x0) and @a f '(@a
   *       x1). The Hermite cubic polynomial is defined
   *       @f[ H_3(x) =
   *       f(x_0)h_0(x)+f(x_1)h_1(x)
   *       +f\prime(x_0)\hat{h}_0(x)+f\prime(x_1)\hat{h}_1(x)
   *       @f]
   *       where
   *       @f$ h_k(x) =
   *       [1-2(x-x_k)L\prime_k(x_k)](L_k(x))^2
   *       @f$
   *       and
   *       @f$ \hat{h}_k(x) =
   *       (x-x_k)(L_k(x))^2
   *       @f$
   *       for the kth-Lagrange coefficient polynomials of degree
   *       n = 1,
   *       @f$ L_0(x) = \frac{x- x_1}{x_0-x_1}@f$ and
   *       @f$ L_1(x) = \frac{x- x_0}{x_1-x_0}@f$.
   *
   *
   *

   *

   *     Akima
   *     @a Akima
   *     5
   *
   *       Non-rounded Akima spline with natural boundary
   *       conditions.  This method uses non-rounded corner
   *       algorithm of Wodicka. This interpolation type uses the
   *       GSL library routines.
   *
   *

   *

   *     Periodic Akima
   *     @a AkimaPeriodic
   *     5
   *
   *       Non-rounded Akima spline with periodic boundary
   *       conditions.  This method uses non-rounded corner
   *       algorithm of Wodicka.  This interpolation type uses
   *       the GSL library routines.
   *
   *

   *

   *
   *
   * Numerical analysis approximation methods for differentiating
   * and integrating unknown functions represented by a data set
   * are implemented by using the interpolations on the data set
   * created with a NumericalApproximation object. The Lagrange polynomial,
   *       @f[ L(x) =
   *       \sum_{k=0}^{n} \left(f(x_k) \prod_{i=0,i \neq k}^{n}
   *            \frac{x- x_i}{x_k-x_i}\right)
   *       @f]
   * is used to determine formulas for numerical differentiation
   * and integration.
   * @n
   * @b Numerical Differentiation
   * @n
   * The differentiation methods use difference formulas to
   * approximate first and second derivatives at a domain value
   * for a given a set of known (@a x,@a y) data points. Once all
   * of the data points are added, a Isis::NumericalApproximation
   * interpolation is used on the dataset to estimate the
   * function, @f$f:x \to y@f$, mapping @a x values to
   * corresponding @a y values. Then, the derivative of
   * @e f evaluated at @a x0 is approximated.
   * To do so, a uniform step size of @a h is chosen to
   * determine the distance between @a xi and @a
   * xi+1.
   *

   *
 Differentiation methods require the parameters
   *     @a a (domain value at which the derivative will be
   *     evaluated), @a n (number of points used in the difference
   *     formula), and @a h (step size of points). Note: @a a
   *     is denoted @a x0 in the above formulas.
   *

   *

   *     Numerical Differentiation Type
   *     Difference Formulas
   *     Methods Available
   *     Description
   *

   *

   *     Backward Difference
   *
   *

   *
 2-point backward difference.
   *           @f[
   *                f\prime(x_0) \approx \frac{1}{h}[f(x_0) - f(x_0 - h)]
   *           @f]
   *
 3-point backward difference.
   *           @f[
   *                f\prime(x_0) \approx \frac{1}{2h}[3f(x_0) - 4f(x_0 - h) +
   *                f(x_0 - 2h)]
   *           @f]
   *
 3-point backward second difference.
   *           @f[
   *                f\prime\prime(x_0) \approx \frac{1}{h^2}[f(x_0)
   *                - 2f(x_0 - h) + f(x_0 - 2h)]
   *           @f]
   *

   *
   *
   *

   *
BackwardFirstDifference()
   *
BackwardSecondDifference()
   *

   *
   *
   *       Backward difference formulas use a uniform step-size
   *       moving in the negative direction from the
   *       given @a x0. These formulas are derived
   *       by differentiating the Lagrange polynomials for
   *       @a xi = @a x0 - @a ih
   *       and evaluating them at@a x0.
   *
   *

   *

   *     Forward Difference
   *
   *

   *
 2-point forward difference.
   *           @f[
   *                f\prime(x_0) \approx \frac{1}{h}[f(x_0 + h) - f(x_0)]
   *           @f]
   *
 3-point forward difference.
   *           @f[
   *                f\prime(x_0) \approx \frac{1}{2h}[-f(x_0 + 2h) +
   *                4f(x_0 + h) - 3f(x_0)]
   *           @f]
   *
 3-point forward second difference.
   *           @f[
   *                f\prime\prime(x_0) \approx \frac{1}{h^2}[f(x_0 +
   *                2h) - 2f(x_0 + h) + f(x_0)]
   *           @f]
   *

   *
   *
   *

   *
ForwardFirstDifference()
   *
ForwardSecondDifference()
   *

   *
   *
   *       Forward difference formulas use a uniform step-size
   *       moving in the positive direction from
   *       @a x0. These formulas are derived by
   *       differentiating the Lagrange polynomials for
   *       @a xi = @a x0 + @a ih
   *       and evaluating them at@a x0.
   *
   *

   *

   *     Center Difference
   *
   *

   *
 3-point center difference.
   *           @f[
   *                f\prime(x_0) \approx \frac{1}{2h}[f(x_0 + h) -
   *                f(x_0 - h)]
   *           @f]
   *
 5-point center difference.
   *           @f[
   *                f\prime(x_0) \approx \frac{1}{12h}[-f(x_0 + 2h) +
   *                8f(x_0 +h) - 8f(x_0 - h) + f(x_0 - 2h)]
   *           @f]
   *
 3-point center second difference.
   *           @f[
   *                f\prime\prime(x_0) \approx \frac{1}{h^2}[f(x_0 + h)
   *                - 2f(x_0) + f(x_0 - h)]
   *           @f]
   *
 5-point center second difference.
   *           @f[
   *                f\prime\prime(x_0) \approx \frac{1}{12h^2}[-f(x_0 +
   *                2h) + 16f(x_0 +h) - 30f(x_0) + 16f(x_0 - h) - f(x_0
   *                - 2h)]
   *           @f]
   *

   *
   *
   *

   *
CenterFirstDifference()
   *
CenterSecondDifference()
   *

   *
   *
   *       Center difference formulas use a uniform step-size
   *       moving in both directions from @a x0 so this
   *       point stays centered. These formulas are derived by
   *       differentiating the Lagrange polynomials for
   *       @a xi = @a x0
   *       + (@e i - (@a n - 1)/2)@a h and evaluating them
   *       at @a x0.
   *
   *

   *

   *     GSL Differentiation
   *      Unknown
   *
   *

   *
GslFirstDerivative()
   *
GslSecondDerivative()
   *

   *
   *
   *       No documentation was found about the algorithms used
   *       for these methods.  They may only be called when
   *       interpolation type is one of the following: @a Linear,
   *       @a Polynomial, @a CubicNatural, @a CubicNatPeriodic, @a Akima, or
   *       @a AkimaPeriodic.
   *
   *

   *

   * @n
   * @b Numerical Integration
   * @n
   * The integration methods were derived using
   * @b Newton-Cotes, or @a quadrature, formulas for
   * approximating integrals given a set of known (@a x,@a y)
   * data points.  The @a x values may be irregularly spaced, but
   * must be unique and added to the dataset in ascending order.
   * Once all of the data points are added, a
   * Isis::NumericalApproximation interpolation is used on
   * the dataset to estimate the function, @f$f:x \to y@f$,
   * mapping @a x values to corresponding @a y values.
   * Then, the integral of @e f on the interval from @a a
   * to @a b is approximated. To do so, an algorithm for creating
   * a composite formula by applying the original formula to
   * segments that share a boundary point is used. The
   * NumericalApproximation::InterpType chosen for
   * interpolation computation seems to have little affect on the
   * error between the actual integral and the return values of
   * the integration methods. The BoolesRule() method varies the
   * most in error. Although errors are not high for any of the
   * interpolation types,
   * @a CubicNatural spline interpolations seem to
   * have the smallest error most often with BoolesRule(). For any
   * other numerical integration method, the errors appear to be
   * identical for any NumericalApproximation::InterpType except
   * @a Polynomial, which usually has a slightly larger
   * error than the other interpolation types. Note: A portion of
   * this algorithm is derived from the IDL function int_tabulated
   * for Boole's Method.
   *

   *
 Integration methods require the parameters @a a,
   *     @a b (interval of domain over which to integrate).
   *

   *

   *     Numerical Integration Type
   *     Integration Formulas
   *     Methods Available
   *     Description
   *

   *

   *     Trapezoidal Rule
   *
   *

   *
 2-point Newton-Cotes trapezoidal rule:
   *            @f[
   *           \int_{a}^b f(x)dx \approx \frac{h}{2}[f(a) + f(b)]
   *            @f]
   *         where @e h = @a b - @a a.
   *

   *
   *
   *

   *
TrapezoidalRule()
   *

   *
   *
   *       The 2-point closed rule, known as the trapezoidal
   *       rule, uses straight line segments between known data
   *       points to estimate the area under the curve (integral).
   *       This Newton-Cotes formula uses a uniform step-size of
   *       @e h = @a b - @a a and is derived by integrating the
   *       Lagrange polynomials over the closed interval [@a a,
   *       @a a+h] for @a xi =
   *       @a x0 + @e ih.
   *
   *

   *

   *     Simpson's Rule, or Simpson's 3-Point Rule
   *
   *

   *
 3-point Newton-Cotes, Simpson's Rule
   *           @f[
   *             \int_{a}^b f(x)dx \approx \frac{h}{3}[f(a) +
   *               4f(a+h) + f(a+2h)]
   *           @f]
   *           where @e h = (@a b - @a a)/2.
   *

   *
   *
   *

   *
Simpsons3PointRule()
   *

   *
   *
   *       The 3-point closed rule, known as Simpson's Rule or
   *       Simpson's 3-Point Rule, uses parabolic arcs between
   *       known data points to estimate the area under the curve
   *       (integral). This Newton-Cotes formula uses a uniform
   *       step-size of @e h = (@a b - @a a)/2 and is derived by
   *       integrating the Lagrange polynomials over the closed
   *       interval [@a a, @a a+2h] for @a xi =
   *       @a x0 + @e ih.
   *
   *

   *

   *     Simpson's 3/8 Rule, or Simpson's 4-Point Rule
   *
   *

   *
 4-point Newton-Cotes, Simpson's 3/8 Rule
   *           @f[
   *             \int_{a}^b f(x)dx \approx \frac{3h}{8}[f(a) +
   *                 3f(a+h) + 3f(a+2h) + f(a+3h)]
   *           @f]
   *           where @e h = (@a b - @a a)/3.
   *

   *
   *
   *

   *
Simpsons4PointRule()
   *

   *
   *
   *       The 4-point closed rule, known as Simpson's 3/8 Rule or
   *       Simpson's 4-Point Rule, uses cubic curves between
   *       known data points to estimate the area under the curve
   *       (integral). This Newton-Cotes formula uses a uniform
   *       step-size of @e h = (@a b - @a a)/3 and is derived by
   *       integrating the Lagrange polynomials over the closed
   *       interval [@a a, @a a+3h] for @a xi =
   *       @a x0 + @e ih.
   *
   *

   *

   *     Boole's Rule
   *
   *

   *
 5-point Newton-Cotes, Boole's Rule
   *           @f[
   *             \int_{a}^b f(x)dx \approx \frac{2h}{45}[7f(a) +
   *                 32f(a+h) + 12f(a+2h) + 32f(a+3h) + 7f(a+4h)]
   *           @f]
   *           where @e h = (@a b - @a a)/4.
   *

   *
   *
   *

   *
BoolesRule()
   *

   *
   *
   *       The 5-point closed rule, known as Boole's Rule, uses
   *       quartic curves between known data points to estimate
   *       the area under the curve (integral). This Newton-Cotes
   *       formula uses a uniform step-size of @e h = (@a b -
   *       @a a)/4 and is derived by integrating the Lagrange
   *       polynomials over the closed interval [@a a, @a a+4h]
   *       for @a xi = @a x0 +
   *       @a ih.
   *
   *

   *

   *     Refinements of Extended Trapezoidal Rule
   *
   *

   *
 Extended closed trapezoidal rule
   *           @f[
   *             \int_{a}^b f(x)dx \approx h[\frac12 f(x_0) +
   *                 f(x_1) + ... + f(x_{n-2} + \frac12
   *                 f(x_{n-1})]
   *           @f]
   *           where @e h = (@a b - @a a)/4, @a x0 =
   *           @a a, and @a xn-1 = @a b.
   *

   *
   *
   *

   *
RefineExtendedTrap()
   *

   *
   *
   *       The extended (or composite) trapezoidal rule can be
   *       used to with a series of refinements to approximate the
   *       integral.  The first stage of refinement returns the
   *       ordinary trapezoidal estimate of the integral.
   *       Subsequent stages will improve the accuracy, where, for
   *       the @a nth stage 2@a n-2  interior
   *       points are added.
   *
   *

   *

   *     Romberg's Method
   *
   *

   *
 Romberg's Method for Integration
   *

   *
   *

   *
RombergsMethod()
   *

   *
   *
   *       Romberg's Method is a potent numerical integration
   *       tool. It uses the extended trapezoidal rule and
   *       Neville's algorithm for polynomial interpolation.
   *
   *

   *

   *     GSL Integration
   *      Unknown
   *
   *

   *
GslIntegral()
   *

   *
   *
   *       No documentation was found about the algorithm used for
   *       this methods.  They may only be called when
   *       interpolation type is one of the following: @a Linear,
   *       @a Polynomial, @a CubicNatural, @a CubicNatPeriodic, @a Akima, or
   *       @a AkimaPeriodic.
   *
   *

   *

   *
   *
   * Below is an example demonstating the use of this class:
   * @code
   *    inline double f(double x) { return (x + cos ((x * x))); }
   *
   *    NumericalApproximation
   *    spline(NumericalApproximation::CubicClamped);
   *    for (int x = 0; x < 200 ; x++) {
   *      spline.addPoint(x, f(x));
   *    }
   *    spline.SetCubicClampedEndptDeriv(0,0);
   *    double yinterp = spline.Evaluate(50.7);
   *    double yextrap =
   *        spline.Evaluate(201,NumericalApproximation::Extrapolate);
   *    double deriv = spline.CenterFirstDifference(10);
   *    double integ =
   *        spline.RombergsMethod(spline.DomainMinimum(),spline.DomainMaximum());
   * @endcode
   *
   * To compute the same data set using the Akima spline, just use the following:
   * @code
   *    spline.Reset(NumericalApproximation::Akima);
   *    double yinterp = spline.Evaluate(50.7);
   *    double yextrap =
   *        spline.Evaluate(201,NumericalApproximation::NearestEndpoint);
   *    double deriv = spline.CenterFirstDifference(10);
   *    double integ =
   *        spline.RombergsMethod(spline.DomainMinimum(),spline.DomainMaximum());
   * @endcode
   *
   *
Caveats

   * When using this class, there are some important details that
   * require consideration.  Several interpolation types of this
   * class use the GSL library. In addition, the GSL library
   * default error handling scheme has implications when used in
   * C++ objects.  The default behavior is to terminate the
   * application when an error occurs.  Options for developers are
   * to turn off error trapping within the GSL, which means users
   * of the library must do their own error checking, or a
   * specific handler can be specified to replace the default
   * mode.  Neither option is well suited in object-oriented
   * implementations because they apply globally, to all users of
   * the library. The approach used here is to turn off error
   * handling and require users to handle their own error
   * checking.  This is not entirely safe as any user that does
   * not follow this strategy could change this behavior at
   * runtime.  Other options should be explored. An additional
   * option is that the default behavior can be altered prior to
   * building GSL but this is also problematic since this behavior
   * cannot be guaranteed universally.
   *
   * For these types, interpolation is strictly adhered to and
   * extrapolation is not handled well. If the input to be
   * evaluated is outside of the
   * minimum and maximum @a x values of the original data
   * points, then the @a y value corresponding to the
   * @a x value that is nearest to that input is returned.
   * However, for sufficiently close values, the clamped cubic and
   * polynomial Neville's types are fairly accurate in their
   * extrapolation techniques.  All differentiation and
   * integration methods throw an error if the value passed is
   * outside of the domain.
   *
   * @ingroup Math
   * @author 2008-11-05 Janet Barrett, Kris Becker, K Teal
   *         Thompson, Jeannie Walldren
   * @internal
   *   @history 1999-08-11 K Teal Thompson - Original version of
   *            NumericalMethods subroutines in Isis2.
   *   @history 2006-06-14 Kris Becker - Created DataInterp class
   *   @history 2007-02-20 Janet Barrett - Created NumericalMethods
   *            class from Isis2 subroutines
   *   @history 2007-02-20 Janet Barrett - Created AtmosModel class
   *            from Isis2 subroutines
   *   @history 2008-06-18 Christopher Austin - Fixed documentation
   *            for DataInterp
   *   @history 2008-06-18 Steven Koechle - Updated NumericalMethods
   *            unitTest.
   *   @history 2008-11-05 Jeannie Walldren - Merged DataInterp
   *            class, NumericalMethods class, and methods from
   *            AtmosModel.  Modified methods from various classes
   *            to be able to function similarly. Added InterpType
   *            @a CubicNeighborhood, difference formula methods, and
   *            integration methods. Created new unitTest.
   *   @history 2008-11-12 Jeannie Walldren - Fixed documentation.
   *   @history 2008-12-18 Jeannie Walldren - Added address
   *            operator (&) to input variables of type vector,
   *            NumericalApproximation::InterpType, and
   *            NumericalApproximation::ExtrapType in
   *            Constructors, AddData(), Evaluate(),
   *            ValueToExtrapolate(), EvaluateCubicNeighborhood(),
   *   @history 2009-01-26 Jeannie Walldren - Fixed error
   *            documentation.
   *   @history 2009-02-12 Jeannie Walldren - Fixed documentation
   *            to include weblinks.
   *   @history 2009-07-17 Steven Lambright - Added algorithm include,
   *            removed macro MAX
   *   @history 2009-06-10 Jeannie Walldren - Added interpolation
   *            type CubicHermite.  This involved modifying some
   *            methods, as well as adding variable p_fprimeOfx
   *            and methods AddCubicHermiteDeriv() and
   *            EvaluateCubicHermite().  Added Contains() method.
   *   @history 2009-07-02 Jeannie Walldren - Replaced Hermite
   *            interpolation algorithm with simpler version and
   *            added methods EvaluateCubicHermiteFirstDeriv and
   *            EvaluateCubicHermiteSecDeriv
   *   @history 2009-07-09 Debbie Cook - Finished Jeannie's
   *            modifications.
   *   @history 2009-08-03 Jeannie Walldren - Clean up code,
   *            documentation and check in changes from
   *            2009-06-10, 2009-07-02, 2009-07-09
   *   @history 2010-07-21 Sharmila Prasad - Remove doxygen documentation
   *            warnings
   *   @history 2010-12-06 Steven Lambright - Optimized AddData(vector, vector),
   *            AddCubicHermiteDeriv(vector), and Init() which causes a
   *            very significant increase in performance when constructing a lot
   *            or adding a lot of data.
   */

  class NumericalApproximation {
    public:
      /**
       * This enum defines the types of interpolation supported in this class
       */
      enum InterpType { Linear,               //!< Linear interpolation.
                        Polynomial,           //!< Polynomial interpolation.
                        PolynomialNeville,    //!< Polynomial interpolation using Neville's algorithm.
                        CubicNatural,         //!< Cubic Spline interpolation with natural boundary conditions.
                        CubicClamped,         //!< Cubic Spline interpolation with clamped boundary conditions.
                        CubicNatPeriodic,     //!< Cubic Spline interpolation with periodic boundary conditions.
                        CubicNeighborhood,    //!< Cubic Spline interpolation using 4-pt Neighborhoods with natural boundary conditions.
                        CubicHermite,         //!< Cubic Spline interpolation using the Hermite cubic polynomial.
                        Akima,                //!< Non-rounded Akima Spline interpolation with natural boundary conditions.
                        AkimaPeriodic         //!< Non-rounded Akima Spline interpolation with periodic boundary conditions.
                      };

      // CONSTRUCTORS
      NumericalApproximation(const NumericalApproximation::InterpType &itype = CubicNatural);
      NumericalApproximation(unsigned int n, double *x, double *y,
                             const NumericalApproximation::InterpType &itype = CubicNatural);
      NumericalApproximation(const vector <double> &x, const vector <double> &y,
                             const NumericalApproximation::InterpType &itype = CubicNatural);
      NumericalApproximation(const NumericalApproximation &dint);
      // ASSIGNMENT OPERATOR
      NumericalApproximation &operator=(const NumericalApproximation &numApMeth);
      // DESTRUCTOR
      virtual ~NumericalApproximation();

      // ACCESSOR METHODS FOR OBJECT PROPERTIES
      string Name() const;
      /**
       * @brief Returns the enumerated type of interpolation chosen
       *
       * This method can be selected after all the points are added in the
       * interpolation by using Compute() method.
       * Note that this prints out as an integer representaion of the
       * enumerated type:
       *

       *
 Linear             = 0
       *
 Polynomial         = 1
       *
 PolynomialNeville  = 2
       *
 CubicNatural       = 3
       *
 CubicClamped       = 4
       *
 CubicNatPeriodic   = 5
       *
 CubicNeighborhood  = 6
       *
 CubicHermite       = 7
       *
 Akima              = 8
       *
 AkimaPeriodic      = 9
       *

       *
       * @return NumericalApproximation::InterpType Currently assigned
       *         interpolation type
       */
      inline InterpType   InterpolationType() {
        return (p_itype);
      }
      int MinPoints();
      int MinPoints(NumericalApproximation::InterpType itype);

      // ADD DATA TO OBJECT
      void AddData(const double x, const double y);
      void AddData(unsigned int n, double *x, double *y);
      void AddData(const vector <double> &x, const vector <double> &y);
      void SetCubicClampedEndptDeriv(const double yp1, const double ypn);
      void AddCubicHermiteDeriv(unsigned int n, double *fprimeOfx);
      void AddCubicHermiteDeriv(const vector <double> &fprimeOfx);
      void AddCubicHermiteDeriv(const double fprimeOfx);

      //ACCESSOR METHODS AFTER DATA IS ENTERED
      double DomainMinimum();
      double DomainMaximum();
      bool Contains(double x);
      /**
       * Returns the number of the coordinates added to the data set.
       *
       * @return @b unsigned @b int Size of data set.
       */
      inline unsigned int Size() {
        return(p_x.size());
      }

      /**
      * This enum defines the manner in which a value outside of the
      * domain should be handled if passed to the Evaluate() method.
      */
      enum ExtrapType { ThrowError,     //!< Evaluate() throws an error if @a a is outside of the domain.
                        Extrapolate,    //!< Evaluate() attempts to extrapolate if @a a is outside of the domain.
                        //!< This is only valid for NumericalApproximation::InterpType @a CubicClamped or @a PolynomialNeville
                        //!< and the result will be accurate only if sufficiently close to the domain boundary.
                        NearestEndpoint //!< Evaluate() returns the y-value of the nearest endpoint if @a a is outside of the domain.
                      };
      // INTERPOLATION AND EXTRAPOLATION RESULTS
      double          Evaluate(const double          a, const ExtrapType &etype = ThrowError);
      vector <double> Evaluate(const vector <double> &a, const ExtrapType &etype = ThrowError);
      vector <double> PolynomialNevilleErrorEstimate();
      vector <double> CubicClampedSecondDerivatives();
      double EvaluateCubicHermiteFirstDeriv(const double a);
      double EvaluateCubicHermiteSecDeriv(const double a);

      // DIFFERENTIATION METHODS
      double GslFirstDerivative(const double a);
      double BackwardFirstDifference(const double a, const unsigned int n = 3,
                                     const double h = 0.1);
      double ForwardFirstDifference(const double a, const unsigned int n = 3,
                                    const double h = 0.1);
      double CenterFirstDifference(const double a, const unsigned int n = 5,
                                   const double h = 0.1);

      double GslSecondDerivative(const double a);
      double BackwardSecondDifference(const double a, const unsigned int n = 3,
                                      const double h = 0.1);
      double ForwardSecondDifference(const double a, const unsigned int n = 3,
                                     const double h = 0.1);
      double CenterSecondDifference(const double a, const unsigned int n = 5,
                                    const double h = 0.1);

      // INTERGRATION METHODS
      double GslIntegral(const double a, const double b);
      double TrapezoidalRule(const double a, const double b);
      double Simpsons3PointRule(const double a, const double b);
      double Simpsons4PointRule(const double a, const double b);
      double BoolesRule(const double a, const double b);
      double RefineExtendedTrap(double a, double b, double s, unsigned int n);
      double RombergsMethod(double a, double b);

      // ASSIGNMENT OPERATORS
      void Reset();
      void Reset(NumericalApproximation::InterpType itype);
      void SetInterpType(NumericalApproximation::InterpType itype);

    protected:
      // == CLASS VARIABLES
      // VARIABLES FOR ALL INTERP TYPES
      InterpType        p_itype;                            //!< Interpolation type
      vector<double>    p_x;                                //!< List of X values
      vector<double>    p_y;                                //!< List of Y values
      bool              p_dataValidated;                    //!< Flag variable to determine whether ValidateDataSet() has been called
      // GSL INTERP VARIABLES
      typedef const gsl_interp_type *InterpFunctor;         //!< GSL Interpolation specs
      typedef map<InterpType, InterpFunctor> FunctorList;   //!< Set up a std::map of GSL interpolator functors.  List of function types
      typedef FunctorList::const_iterator FunctorConstIter; //!< GSL Iterator
      gsl_interp_accel *p_acc;                              //!< Lookup accelorator
      gsl_spline       *p_interp;                           //!< Currently active interpolator
      static FunctorList p_interpFunctors;                   //!< Maintains list of interpolator options
      // CUBIC CLAMPED VARIABLES
      bool              p_clampedEndptsSet;                 //!< Flag variable to determine whether SetCubicClampedEndptDeriv() has been called after all data was added for @a CubicClamped interpolation.
      bool              p_clampedComputed;                  //!< Flag variable to determine whether ComputeCubicClamped() has been called.
      double            p_clampedDerivFirstPt;              //!< First derivative of first x-value, p_x[0].  This is only used for the @a CubicClamped interpolation type.
      double            p_clampedDerivLastPt;               //!< First derivative of last x-value, p_x[n-1].  This is only used for the @a CubicClamped interpolation type.
      vector<double>    p_clampedSecondDerivs;              //!< List of second derivatives evaluated at p_x values.  This is only used for the @a CubicClamped interpolation type.
      // POLYNOMIAL NEVILLE VARIABLES
      vector <double>   p_polyNevError;                     //!< Estimate of error for interpolation evaluated at x.  This is only used for the @a PolynomialNeville interpolation type. 91 taken from AtmosModel.
      // CUBIC HERMITE VARIABLES
      vector<double>    p_fprimeOfx;                        //!< List of first derivatives corresponding to each x value in the data set (i.e. each value in p_x)

      // == PROTECTED METHODS
      // CREATING, DESTROYING OBJECT
      void Init(NumericalApproximation::InterpType itype);
      bool GslInterpType(NumericalApproximation::InterpType itype) const;
      void GslAllocation(unsigned int npoints);
      void GslDeallocation();
      InterpFunctor GslFunctor(NumericalApproximation::InterpType itype) const;
      // VERIFICATION METHODS
      void GslIntegrityCheck(int gsl_status, const char *src,
                             int line);
      void ValidateDataSet();
      bool InsideDomain(const double a);
      // COMPUTATION AND EVALUATION METHODS
      bool   GslComputed() const;
      void   ComputeGsl();
      void   ComputeCubicClamped();
      double ValueToExtrapolate(const double a, const ExtrapType &etype);
      double          EvaluateCubicNeighborhood(const double a);
      vector <double> EvaluateCubicNeighborhood(const vector <double> &a, const ExtrapType &etype);
      double          EvaluateCubicClamped(const double a);
      double          EvaluateCubicHermite(const double a);
      double          EvaluatePolynomialNeville(const double a);
      vector <double> EvaluateForIntegration(const double a, const double b,
                                             const unsigned int n);
      int FindIntervalLowerIndex(const double a);
      // STANDARDIZE ERRORS
      void ReportException(IException::ErrorType type, const string &method,
                           const string &message, const char *filesrc,
                           int lineno) const;
  };
};

#endif


text:
#ifndef _SNPLIB_SRC_MATH_LIB_H_
#define _SNPLIB_SRC_MATH_LIB_H_

#ifdef USE_MKL
#include <mkl.h>
#define set_num_threads mkl_set_num_threads
#endif
#ifdef USE_OPENBLAS
#ifdef _MSC_VER
#include <complex>
#define lapack_complex_float std::complex<float>
#define lapack_complex_double std::complex<double>
#endif
#include <cblas.h>
#include <lapacke.h>
#define set_num_threads openblas_set_num_threads
#endif

#endif //_SNPLIB_SRC_MATH_LIB_H_

max_stars_repo_name: VasLem/SNPLIB

max_issues_repo_name: VasLem/SNPLIB

max_forks_repo_name: VasLem/SNPLIB


text:
/**
 *
 * @file core_cgetrf.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Mathieu Faverge
 * @date 2010-11-15
 * @generated c Tue Jan  7 11:44:48 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 * @ingroup CORE_PLASMA_Complex32_t
 *
 *  CORE_cgetrf - Computes an LU factorization of a general M-by-N matrix A
 *  using the tile LU algorithm with partial tile pivoting with row interchanges.
 *
 *******************************************************************************
 *
 * @param[in] m
 *          The number of rows of the matrix A. m >= 0.
 *
 * @param[in] n
 *          The number of columns of the matrix A. n >= 0.
 *
 * @param[in,out] A
 *          On entry, the M-by-N matrix to be factored.
 *          On exit, the tile factors L and U from the factorization.
 *
 * @param[in] lda
 *          The leading dimension of the array A. LDA >= max(1,M).
 *
 * @param[out] IPIV
 *          The pivot indices that define the permutations.
 *
 * @param[out] info
 *          - 0 on successful exit
 *          - <0 if -i, the i-th argument had an illegal value
 *          - >0 if i, U(i,i) is exactly zero. The factorization has been
 *            completed, but the factor U is exactly singular, and division by
 *            zero will occur if it is used to solve a system of equations.
 *
 *******************************************************************************
 *
 * @return
 *         \retval PLASMA_SUCCESS successful exit
 *
 ******************************************************************************/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_cgetrf = PCORE_cgetrf
#define CORE_cgetrf PCORE_cgetrf
#endif
int CORE_cgetrf(int m, int n,
                 PLASMA_Complex32_t *A, int lda,
                 int *IPIV, int *info)
{
    *info = LAPACKE_cgetrf_work(LAPACK_COL_MAJOR, m, n, A, lda, IPIV );
    return PLASMA_SUCCESS;
}


text:
/**********************************\
* Test benchmark                   *
*                                  *
* by                               *
*   Elliott Forney                 *
*   3.15.2010                      *
\**********************************/

/*
 *  Libraries
 */

#include <stdio.h>
#include <stdlib.h>
#include <stdio.h>
#include <stdbool.h>
#include <unistd.h>
#include <getopt.h>
#include <math.h>

#include <cblas.h>

#include "matrix.h"
#include "benchmark.h"
#include "errcheck.h"

/*
 *  Macros
 */

// print command line usage
#define print_usage() fprintf(stdout, "Usage: %s [-m rows] [-n cols] [-s]\n", arg[0])

// acceptable maximum relative error tolerance
#define TOLERANCE 0.01

#define DEBUG 0

/*
 *  Global variables
 */

// flag for simple, one line output
bool simple_out = false;

/*
 *  Function prototypes
 */

// parse command line arguments
void parse_args(int narg, char **arg,
                long long *m, long long *n);

/*
 *  Function bodies
 */

// setup network
int main(int narg, char **arg)
{
  // default matrix dimensions
  long long m = 6400;
  long long n = 0;

  // parse command line arguments
  parse_args(narg, arg, &m, &n);

  // if n or p == 0 then set to m
  if (n == 0)
    n = m;

  // load test matrix b
  matrix a;
  matrix_init(&a, m, n);
  matrix_load_runif(a, 0, 1);

  // create matrix to hold result
  matrix result;
  matrix_init(&result, m, n);

  // run warm up
  matrix_phi_prime(result, a);

  // wait for all kernels to finish
  matrix_wait();

  // create a new benchmark timer
  benchmark ben;
  benchmark_init(&ben);

  // start timer
  benchmark_start_timer(&ben);

  // run multiplication
  matrix_phi_prime(result, a);

  // wait for all kernels to finish
  matrix_wait();

  // stop timer
  benchmark_stop_timer(&ben);

  // figure giga floating point operatins per second
  benchmark_add_flop(&ben, 3ll*m*n);
  double gflops = benchmark_check_gflops(ben);
  double time   = benchmark_check_timer(ben);

  // print time and gbytess
  printf("Time: %f  GFlopS: %f\n",
         time, gflops);

  // clean up
  matrix_dest(&a);
  matrix_dest(&result);
  benchmark_dest(&ben);

  // Come back soon now ya'hear!
  return 0;
}

// parse command line arguments
void parse_args(int narg, char **arg,
                long long *m, long long *n)
{
  int opt; // getopt output

  // for each argument
  while ((opt = getopt(narg, arg, "sm:n:")) != -1)
  {
    if (opt == 's')
      simple_out = true;

    else if (opt == 'm')
      *m = (long long)atoi(arg[optind-1]);

    else if (opt == 'n')
      *n = (long long)atoi(arg[optind-1]);

    // print usage and quit on unknown option
    else
    {
      print_usage();
      exit(1);
    }
  }

  // assume last non dash arg is m
  if (optind < narg)
    *m = (long long)atoi(arg[optind]);
}

max_stars_repo_name: idfah/badger

max_issues_repo_name: idfah/badger

max_forks_repo_name: idfah/badger


text:
/* matrix/gsl_matrix_int.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 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 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#ifndef __GSL_MATRIX_INT_H__
#define __GSL_MATRIX_INT_H__

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector_int.h>

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

__BEGIN_DECLS

typedef struct
{
  size_t size1;
  size_t size2;
  size_t tda;
  int * data;
  gsl_block_int * block;
  int owner;
} gsl_matrix_int;

typedef struct
{
  gsl_matrix_int matrix;
} _gsl_matrix_int_view;

typedef _gsl_matrix_int_view gsl_matrix_int_view;

typedef struct
{
  gsl_matrix_int matrix;
} _gsl_matrix_int_const_view;

typedef const _gsl_matrix_int_const_view gsl_matrix_int_const_view;

/* Allocation */

GSL_EXPORT
gsl_matrix_int *
gsl_matrix_int_alloc (const size_t n1, const size_t n2);

GSL_EXPORT
gsl_matrix_int *
gsl_matrix_int_calloc (const size_t n1, const size_t n2);

GSL_EXPORT
gsl_matrix_int *
gsl_matrix_int_alloc_from_block (gsl_block_int * b,
                                   const size_t offset,
                                   const size_t n1,
                                   const size_t n2,
                                   const size_t d2);

GSL_EXPORT
gsl_matrix_int *
gsl_matrix_int_alloc_from_matrix (gsl_matrix_int * m,
                                    const size_t k1,
                                    const size_t k2,
                                    const size_t n1,
                                    const size_t n2);

GSL_EXPORT
gsl_vector_int *
gsl_vector_int_alloc_row_from_matrix (gsl_matrix_int * m,
                                        const size_t i);

GSL_EXPORT
gsl_vector_int *
gsl_vector_int_alloc_col_from_matrix (gsl_matrix_int * m,
                                        const size_t j);

GSL_EXPORT void gsl_matrix_int_free (gsl_matrix_int * m);

/* Views */

GSL_EXPORT
_gsl_matrix_int_view
gsl_matrix_int_submatrix (gsl_matrix_int * m,
                            const size_t i, const size_t j,
                            const size_t n1, const size_t n2);

GSL_EXPORT
_gsl_vector_int_view
gsl_matrix_int_row (gsl_matrix_int * m, const size_t i);

GSL_EXPORT
_gsl_vector_int_view
gsl_matrix_int_column (gsl_matrix_int * m, const size_t j);

GSL_EXPORT
_gsl_vector_int_view
gsl_matrix_int_diagonal (gsl_matrix_int * m);

GSL_EXPORT
_gsl_vector_int_view
gsl_matrix_int_subdiagonal (gsl_matrix_int * m, const size_t k);

GSL_EXPORT
_gsl_vector_int_view
gsl_matrix_int_superdiagonal (gsl_matrix_int * m, const size_t k);

GSL_EXPORT
_gsl_matrix_int_view
gsl_matrix_int_view_array (int * base,
                             const size_t n1,
                             const size_t n2);

GSL_EXPORT
_gsl_matrix_int_view
gsl_matrix_int_view_array_with_tda (int * base,
                                      const size_t n1,
                                      const size_t n2,
                                      const size_t tda);

GSL_EXPORT
_gsl_matrix_int_view
gsl_matrix_int_view_vector (gsl_vector_int * v,
                              const size_t n1,
                              const size_t n2);

GSL_EXPORT
_gsl_matrix_int_view
gsl_matrix_int_view_vector_with_tda (gsl_vector_int * v,
                                       const size_t n1,
                                       const size_t n2,
                                       const size_t tda);

GSL_EXPORT
_gsl_matrix_int_const_view
gsl_matrix_int_const_submatrix (const gsl_matrix_int * m,
                                  const size_t i, const size_t j,
                                  const size_t n1, const size_t n2);

GSL_EXPORT
_gsl_vector_int_const_view
gsl_matrix_int_const_row (const gsl_matrix_int * m,
                            const size_t i);

GSL_EXPORT
_gsl_vector_int_const_view
gsl_matrix_int_const_column (const gsl_matrix_int * m,
                               const size_t j);

GSL_EXPORT
_gsl_vector_int_const_view
gsl_matrix_int_const_diagonal (const gsl_matrix_int * m);

GSL_EXPORT
_gsl_vector_int_const_view
gsl_matrix_int_const_subdiagonal (const gsl_matrix_int * m,
                                    const size_t k);

GSL_EXPORT
_gsl_vector_int_const_view
gsl_matrix_int_const_superdiagonal (const gsl_matrix_int * m,
                                      const size_t k);

GSL_EXPORT
_gsl_matrix_int_const_view
gsl_matrix_int_const_view_array (const int * base,
                                   const size_t n1,
                                   const size_t n2);

GSL_EXPORT
_gsl_matrix_int_const_view
gsl_matrix_int_const_view_array_with_tda (const int * base,
                                            const size_t n1,
                                            const size_t n2,
                                            const size_t tda);

GSL_EXPORT
_gsl_matrix_int_const_view
gsl_matrix_int_const_view_vector (const gsl_vector_int * v,
                                    const size_t n1,
                                    const size_t n2);

GSL_EXPORT
_gsl_matrix_int_const_view
gsl_matrix_int_const_view_vector_with_tda (const gsl_vector_int * v,
                                             const size_t n1,
                                             const size_t n2,
                                             const size_t tda);

/* Operations */

GSL_EXPORT int   gsl_matrix_int_get(const gsl_matrix_int * m, const size_t i, const size_t j);
GSL_EXPORT void    gsl_matrix_int_set(gsl_matrix_int * m, const size_t i, const size_t j, const int x);

GSL_EXPORT int * gsl_matrix_int_ptr(gsl_matrix_int * m, const size_t i, const size_t j);
GSL_EXPORT const int * gsl_matrix_int_const_ptr(const gsl_matrix_int * m, const size_t i, const size_t j);

GSL_EXPORT void gsl_matrix_int_set_zero (gsl_matrix_int * m);
GSL_EXPORT void gsl_matrix_int_set_identity (gsl_matrix_int * m);
GSL_EXPORT void gsl_matrix_int_set_all (gsl_matrix_int * m, int x);

GSL_EXPORT int gsl_matrix_int_fread (FILE * stream, gsl_matrix_int * m) ;
GSL_EXPORT int gsl_matrix_int_fwrite (FILE * stream, const gsl_matrix_int * m) ;
GSL_EXPORT int gsl_matrix_int_fscanf (FILE * stream, gsl_matrix_int * m);
GSL_EXPORT int gsl_matrix_int_fprintf (FILE * stream, const gsl_matrix_int * m, const char * format);

GSL_EXPORT int gsl_matrix_int_memcpy(gsl_matrix_int * dest, const gsl_matrix_int * src);
GSL_EXPORT int gsl_matrix_int_swap(gsl_matrix_int * m1, gsl_matrix_int * m2);

GSL_EXPORT int gsl_matrix_int_swap_rows(gsl_matrix_int * m, const size_t i, const size_t j);
GSL_EXPORT int gsl_matrix_int_swap_columns(gsl_matrix_int * m, const size_t i, const size_t j);
GSL_EXPORT int gsl_matrix_int_swap_rowcol(gsl_matrix_int * m, const size_t i, const size_t j);
GSL_EXPORT int gsl_matrix_int_transpose (gsl_matrix_int * m);
GSL_EXPORT int gsl_matrix_int_transpose_memcpy (gsl_matrix_int * dest, const gsl_matrix_int * src);

GSL_EXPORT int gsl_matrix_int_max (const gsl_matrix_int * m);
GSL_EXPORT int gsl_matrix_int_min (const gsl_matrix_int * m);
GSL_EXPORT void gsl_matrix_int_minmax (const gsl_matrix_int * m, int * min_out, int * max_out);

GSL_EXPORT void gsl_matrix_int_max_index (const gsl_matrix_int * m, size_t * imax, size_t *jmax);
GSL_EXPORT void gsl_matrix_int_min_index (const gsl_matrix_int * m, size_t * imin, size_t *jmin);
GSL_EXPORT void gsl_matrix_int_minmax_index (const gsl_matrix_int * m, size_t * imin, size_t * jmin, size_t * imax, size_t * jmax);

GSL_EXPORT int gsl_matrix_int_isnull (const gsl_matrix_int * m);

GSL_EXPORT int gsl_matrix_int_add (gsl_matrix_int * a, const gsl_matrix_int * b);
GSL_EXPORT int gsl_matrix_int_sub (gsl_matrix_int * a, const gsl_matrix_int * b);
GSL_EXPORT int gsl_matrix_int_mul_elements (gsl_matrix_int * a, const gsl_matrix_int * b);
GSL_EXPORT int gsl_matrix_int_div_elements (gsl_matrix_int * a, const gsl_matrix_int * b);
GSL_EXPORT int gsl_matrix_int_scale (gsl_matrix_int * a, const double x);
GSL_EXPORT int gsl_matrix_int_add_constant (gsl_matrix_int * a, const double x);
GSL_EXPORT int gsl_matrix_int_add_diagonal (gsl_matrix_int * a, const double x);

/***********************************************************************/
/* The functions below are obsolete                                    */
/***********************************************************************/
GSL_EXPORT int gsl_matrix_int_get_row(gsl_vector_int * v, const gsl_matrix_int * m, const size_t i);
GSL_EXPORT int gsl_matrix_int_get_col(gsl_vector_int * v, const gsl_matrix_int * m, const size_t j);
GSL_EXPORT int gsl_matrix_int_set_row(gsl_matrix_int * m, const size_t i, const gsl_vector_int * v);
GSL_EXPORT int gsl_matrix_int_set_col(gsl_matrix_int * m, const size_t j, const gsl_vector_int * v);

/* inline functions if you are using GCC */

#ifdef HAVE_INLINE
extern inline
int
gsl_matrix_int_get(const gsl_matrix_int * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_VAL("first index out of range", GSL_EINVAL, 0) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_VAL("second index out of range", GSL_EINVAL, 0) ;
    }
#endif
  return m->data[i * m->tda + j] ;
}

extern inline
void
gsl_matrix_int_set(gsl_matrix_int * m, const size_t i, const size_t j, const int x)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_VOID("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_VOID("second index out of range", GSL_EINVAL) ;
    }
#endif
  m->data[i * m->tda + j] = x ;
}

extern inline
int *
gsl_matrix_int_ptr(gsl_matrix_int * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
    }
#endif
  return (int *) (m->data + (i * m->tda + j)) ;
}

extern inline
const int *
gsl_matrix_int_const_ptr(const gsl_matrix_int * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
    }
#endif
  return (const int *) (m->data + (i * m->tda + j)) ;
}

#endif

__END_DECLS

#endif /* __GSL_MATRIX_INT_H__ */

max_stars_repo_name: dynaryu/vaws

max_issues_repo_name: dynaryu/vaws

max_forks_repo_name: dynaryu/vaws


text:
#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_vegas.h>
#include "integral.c"
#include "timer.c"
#include "timer.h"
#define N 20
double g (double *t,size_t dim, void *params);

int main (void)
{
	double res, err;
	// 6 dimensional integral
	size_t dim = 6;

	int dims=6;
	//evaluating all dimensions from 0 to 1
	double x1[] = { 0., 0., 0. ,0.,0.,0. };
	double xu[] = { 1., 1., 1. ,1.,1.,1. };
	FILE *file;
	file=fopen("data.txt","w");

	double rmin=1.001;
	double rmax=4;
	double rtot=rmin;
	double rstep= (rmax-rmin)/(N);

	gsl_rng *r = gsl_rng_alloc (gsl_rng_taus2);
	unsigned long seed = 1UL;
	gsl_rng_set (r, seed);

	// 10^6 function evaluations per iteration

	size_t calls = 1000000;
	gsl_monte_function G= { &g, dim, &rtot};
	gsl_monte_vegas_state *sv = gsl_monte_vegas_alloc (dim);
	gsl_monte_vegas_init (sv);

	//Vegas integration
	double vegas[N];
	double rval[N];
	double di[N];
	timer_start ();
	for (int i = 0; i<=N-1; i++)
	{
		gsl_monte_vegas_integrate (&G,x1,xu,dim,calls/5,r,sv,&res,&err);
		do
		{
			gsl_monte_vegas_integrate (&G, x1,xu,dim,calls,r,sv,&res,&err);
			fflush (stdout);
		}

		while (fabs (gsl_monte_vegas_chisq (sv)-1.0) >0.2);
		rtot+=rstep;
		vegas[i]=res;
		rval[i]=rtot;
		di[i] = -.2 / pow (rtot,3.);
	}
	timer_stop ();
	double tvegas=timer_stop();
	gsl_monte_vegas_free (sv);
	rtot=rmin;
	printf (" %10.6f\n",tvegas);
	//homemade
	rtot=rmin;
	double x[6];

   	r = gsl_rng_alloc (gsl_rng_taus2);
   	gsl_rng_set (r, 1UL);
	double thome;
	long i, j;
	double r2=0.;
	double home[20];
	for (j = 0; j < N; j++)
   	{
		r2=0.;
       		for (i = 0; i < 1.e6; i++)
        	{
            		for (int k = 0; k <= dims; k++)
            		{
                		x[k] = gsl_rng_uniform (r) ;
            		}
		r2+= g(x,dim, &rtot);

        }
       res = r2/1.e6;
	fflush (stdout);
	rtot+=rstep;
	home [j]=res;
    }
	timer_stop();
	rtot=rmin;
	thome=timer_stop();
	printf (" %10.6f\n",thome);
	for (int t=0;t<N;t++)
        {
                fprintf (file, "%10.6f     %10.6f    %10.6f  %10.6f \n " , rval[t],fabs(-vegas[t]),fabs(-di[t]),fabs(-home[t]));
                rtot+=rstep;
        }

	return 0;
}


text:
#pragma once

// STL and std libs
#include <map>
#include <set>
#include <vector>
#include <list>
#include <memory>
#include <functional>
#include <iostream>
#include <fstream>
#include <algorithm>
#include <chrono>
#include <type_traits>

// Guildline Support Library
#include <gsl.h>

//// boost and std extension
//#include <boost\signals2.hpp>
//#include <boost\any.hpp>
//#include <boost\range\iterator_range_core.hpp>
//#include <boost\range\adaptor\transformed.hpp>
//#include <boost\operators.hpp>
//#include <boost\format.hpp>
//#include <boost\filesystem.hpp>

#include <stride_range.h>
#include <tree.h>
#include <minmax>

#include "Math3D.h"
#include "SmartPointers.h"
#include "String.h"

//#if defined __AVX__
//#undef __AVX__ //#error Eigen have problem with AVX now
//#endif
//#define EIGEN_HAS_CXX11_MATH 1
//#define EIGEN_HAS_STD_RESULT_OF 1
//#define EIGEN_HAS_VARIADIC_TEMPLATES 1
//#include <Eigen\Dense>

namespace Causality
{
	using time_seconds = std::chrono::duration<double>;
	typedef uint64_t id_t;

	using std::string;
	using std::iterator_range;
	//using boost::sub_range;

	//namespace adaptors = boost::adaptors;

	using gsl::owner;
	using gsl::byte;
	using gsl::not_null;

	using stdx::tree_node;
	using stdx::foward_tree_node;
	using stdx::stride_range;
	using stdx::stride_iterator;

	using std::vector;
	using std::map;
	using std::function;
	using std::unique_ptr;
	using std::shared_ptr;
	using std::list;
	using std::weak_ptr;
}

max_stars_repo_name: ArcEarth/SrInspection

max_issues_repo_name: ArcEarth/SrInspection

max_forks_repo_name: ArcEarth/SrInspection


text:
#ifndef COMMON_H
#define COMMON_H

#include <assert.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_statistics.h>

#include <plplot/plplot.h>

#endif

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
/* integration/qk15.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 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 <config.h>
#include <gsl/gsl_integration.h>

/* Gauss quadrature weights and kronrod quadrature abscissae and
   weights as evaluated with 80 decimal digit arithmetic by
   L. W. Fullerton, Bell Labs, Nov. 1981. */

static const double xgk[8] =    /* abscissae of the 15-point kronrod rule */
{
  0.991455371120812639206854697526329,
  0.949107912342758524526189684047851,
  0.864864423359769072789712788640926,
  0.741531185599394439863864773280788,
  0.586087235467691130294144838258730,
  0.405845151377397166906606412076961,
  0.207784955007898467600689403773245,
  0.000000000000000000000000000000000
};

/* xgk[1], xgk[3], ... abscissae of the 7-point gauss rule.
   xgk[0], xgk[2], ... abscissae to optimally extend the 7-point gauss rule */

static const double wg[4] =     /* weights of the 7-point gauss rule */
{
  0.129484966168869693270611432679082,
  0.279705391489276667901467771423780,
  0.381830050505118944950369775488975,
  0.417959183673469387755102040816327
};

static const double wgk[8] =    /* weights of the 15-point kronrod rule */
{
  0.022935322010529224963732008058970,
  0.063092092629978553290700663189204,
  0.104790010322250183839876322541518,
  0.140653259715525918745189590510238,
  0.169004726639267902826583426598550,
  0.190350578064785409913256402421014,
  0.204432940075298892414161999234649,
  0.209482141084727828012999174891714
};

void
gsl_integration_qk15 (const gsl_function * f, double a, double b,
      double *result, double *abserr,
      double *resabs, double *resasc)
{
  double fv1[8], fv2[8];
  gsl_integration_qk (8, xgk, wg, wgk, fv1, fv2, f, a, b, result, abserr, resabs, resasc);
}

max_stars_repo_name: ielomariala/Hex-Game


text:
#ifndef ConditionalTest_SERIAL_H
#define ConditionalTest_SERIAL_H
#include <ace/core/core.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_statistics_double.h>
#include <gsl/gsl_randist.h>
#include "conditionaltest.h"

class ConditionalTest::Serial : public EAbstractAnalyticSerial
{
    Q_OBJECT
public:
    explicit Serial(ConditionalTest* parent);
    virtual std::unique_ptr<EAbstractAnalyticBlock> execute(const EAbstractAnalyticBlock* block) override final;

    // helper functions
    bool isEmpty(QVector<QVector<double>>& matrix);

    // Statistical Tests
    void hypergeom(
        const QVector<QString>& amx_column,
        const CCMatrix::Pair& ccmPair,
        int clusterIndex,
        int featureIndex,
        int labelIndex,
        double& results);

    void test_proportions(
        const QVector<QString>& amx_column,
        const CCMatrix::Pair& ccmPair,
        int clusterIndex,
        int featureIndex,
        int labelIndex,
        double& results);

    void regression(
        const QVector<QString>& amx_column,
        const CCMatrix::Pair& ccmPair,
        int clusterIndex,
        QVector<double>& results);

private:
    /*!
     * Pointer to the base analytic for this object.
     */
    ConditionalTest* _base;

    // Performs power analysis for multiple-linear regression.
    double pwr_f2_test(int u, int v, double f2, double sig_level);
};

#endif

max_stars_repo_name: SystemsGenetics/KINC

max_issues_repo_name: SystemsGenetics/KINC

max_forks_repo_name: SystemsGenetics/KINC