text:
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_multifit_nlinear.h>

/* parameters to model */
struct model_params
{
  double a1;
  double a2;
  double a3;
  double a4;
  double a5;
};

/* Branin function */
int
func_f (const gsl_vector * x, void *params, gsl_vector * f)
{
  struct model_params *par = (struct model_params *) params;
  double x1 = gsl_vector_get(x, 0);
  double x2 = gsl_vector_get(x, 1);
  double f1 = x2 + par->a1 * x1 * x1 + par->a2 * x1 + par->a3;
  double f2 = sqrt(par->a4) * sqrt(1.0 + (1.0 - par->a5) * cos(x1));

  gsl_vector_set(f, 0, f1);
  gsl_vector_set(f, 1, f2);

  return GSL_SUCCESS;
}

int
func_df (const gsl_vector * x, void *params, gsl_matrix * J)
{
  struct model_params *par = (struct model_params *) params;
  double x1 = gsl_vector_get(x, 0);
  double f2 = sqrt(par->a4) * sqrt(1.0 + (1.0 - par->a5) * cos(x1));

  gsl_matrix_set(J, 0, 0, 2.0 * par->a1 * x1 + par->a2);
  gsl_matrix_set(J, 0, 1, 1.0);

  gsl_matrix_set(J, 1, 0, -0.5 * par->a4 / f2 * (1.0 - par->a5) * sin(x1));
  gsl_matrix_set(J, 1, 1, 0.0);

  return GSL_SUCCESS;
}

int
func_fvv (const gsl_vector * x, const gsl_vector * v,
          void *params, gsl_vector * fvv)
{
  struct model_params *par = (struct model_params *) params;
  double x1 = gsl_vector_get(x, 0);
  double v1 = gsl_vector_get(v, 0);
  double c = cos(x1);
  double s = sin(x1);
  double f2 = sqrt(par->a4) * sqrt(1.0 + (1.0 - par->a5) * c);
  double t = 0.5 * par->a4 * (1.0 - par->a5) / f2;

  gsl_vector_set(fvv, 0, 2.0 * par->a1 * v1 * v1);
  gsl_vector_set(fvv, 1, -t * (c + s*s/f2) * v1 * v1);

  return GSL_SUCCESS;
}

void
callback(const size_t iter, void *params,
         const gsl_multifit_nlinear_workspace *w)
{
  gsl_vector * x = gsl_multifit_nlinear_position(w);
  double x1 = gsl_vector_get(x, 0);
  double x2 = gsl_vector_get(x, 1);

  /* print out current location */
  printf("%f %f\n", x1, x2);
}

void
solve_system(gsl_vector *x0, gsl_multifit_nlinear_fdf *fdf,
             gsl_multifit_nlinear_parameters *params)
{
  const gsl_multifit_nlinear_type *T = gsl_multifit_nlinear_trust;
  const size_t max_iter = 200;
  const double xtol = 1.0e-8;
  const double gtol = 1.0e-8;
  const double ftol = 1.0e-8;
  const size_t n = fdf->n;
  const size_t p = fdf->p;
  gsl_multifit_nlinear_workspace *work =
    gsl_multifit_nlinear_alloc(T, params, n, p);
  gsl_vector * f = gsl_multifit_nlinear_residual(work);
  gsl_vector * x = gsl_multifit_nlinear_position(work);
  int info;
  double chisq0, chisq, rcond;

  printf("# %s/%s\n",
         gsl_multifit_nlinear_name(work),
         gsl_multifit_nlinear_trs_name(work));

  /* initialize solver */
  gsl_multifit_nlinear_init(x0, fdf, work);

  /* store initial cost */
  gsl_blas_ddot(f, f, &chisq0);

  /* iterate until convergence */
  gsl_multifit_nlinear_driver(max_iter, xtol, gtol, ftol,
                              callback, NULL, &info, work);

  /* store final cost */
  gsl_blas_ddot(f, f, &chisq);

  /* store cond(J(x)) */
  gsl_multifit_nlinear_rcond(&rcond, work);

  /* print summary */
  fprintf(stderr, "%-25s %-6zu %-5zu %-5zu %-13.4e %-12.4e %-13.4e (%.2e, %.2e)\n",
          gsl_multifit_nlinear_trs_name(work),
          gsl_multifit_nlinear_niter(work),
          fdf->nevalf,
          fdf->nevaldf,
          chisq0,
          chisq,
          1.0 / rcond,
          gsl_vector_get(x, 0),
          gsl_vector_get(x, 1));

  printf("\n\n");

  gsl_multifit_nlinear_free(work);
}

int
main (void)
{
  const size_t n = 2;
  const size_t p = 2;
  gsl_vector *f = gsl_vector_alloc(n);
  gsl_vector *x = gsl_vector_alloc(p);
  gsl_multifit_nlinear_fdf fdf;
  gsl_multifit_nlinear_parameters fdf_params =
    gsl_multifit_nlinear_default_parameters();
  struct model_params params;

  params.a1 = -5.1 / (4.0 * M_PI * M_PI);
  params.a2 = 5.0 / M_PI;
  params.a3 = -6.0;
  params.a4 = 10.0;
  params.a5 = 1.0 / (8.0 * M_PI);

  /* print map of Phi(x1, x2) */
  {
    double x1, x2, chisq;

    for (x1 = -5.0; x1 < 15.0; x1 += 0.1)
      {
        for (x2 = -5.0; x2 < 15.0; x2 += 0.1)
          {
            gsl_vector_set(x, 0, x1);
            gsl_vector_set(x, 1, x2);
            func_f(x, &params, f);

            gsl_blas_ddot(f, f, &chisq);

            printf("%f %f %f\n", x1, x2, chisq);
          }
        printf("\n");
      }
    printf("\n\n");
  }

  /* define function to be minimized */
  fdf.f = func_f;
  fdf.df = func_df;
  fdf.fvv = func_fvv;
  fdf.n = n;
  fdf.p = p;
  fdf.params = &params;

  /* starting point */
  gsl_vector_set(x, 0, 6.0);
  gsl_vector_set(x, 1, 14.5);

  fprintf(stderr, "%-25s %-6s %-5s %-5s %-13s %-12s %-13s %-15s\n",
          "Method", "NITER", "NFEV", "NJEV", "Initial Cost",
          "Final cost", "Final cond(J)", "Final x");

  fdf_params.trs = gsl_multifit_nlinear_trs_lm;
  solve_system(x, &fdf, &fdf_params);

  fdf_params.trs = gsl_multifit_nlinear_trs_lmaccel;
  solve_system(x, &fdf, &fdf_params);

  fdf_params.trs = gsl_multifit_nlinear_trs_dogleg;
  solve_system(x, &fdf, &fdf_params);

  fdf_params.trs = gsl_multifit_nlinear_trs_ddogleg;
  solve_system(x, &fdf, &fdf_params);

  fdf_params.trs = gsl_multifit_nlinear_trs_subspace2D;
  solve_system(x, &fdf, &fdf_params);

  gsl_vector_free(f);
  gsl_vector_free(x);

  return 0;
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_linalg.h>

/* number of data points to fit */
#define N        6

/* number of fit coefficients */
#define NCOEFFS  6

/* Spline order */
#define K  3

/* nbreak = ncoeffs + 2 - k */
#define NBREAK   (NCOEFFS + 2 - K)

void gsl_solve( gsl_matrix *A,  const size_t n,
                double xmax, double xmin,
                double *B )
{

  printf("\n *** Solve with GSL *** \n");
  gsl_vector *d; d = gsl_vector_alloc(n);
  gsl_vector *e; e = gsl_vector_alloc(n-1);
  gsl_vector *f; f = gsl_vector_alloc(n-1);
  gsl_vector *b; b = gsl_vector_alloc(n);
  gsl_vector *x; x = gsl_vector_alloc(n);

  for (int i = 0; i < n; i++)
    {
      gsl_vector_set(d, i, gsl_matrix_get(A,i,i));
      if (i>0)
        gsl_vector_set(f, i-1, gsl_matrix_get(A,i,i-1));
      if (i<n-1)
        gsl_vector_set(e, i, gsl_matrix_get(A,i,i+1));
    }

  for (int i=0; i<n; i++)
    gsl_vector_set(b,i,B[i]);

  gsl_linalg_solve_tridiag ( d, e, f, b, x);

  for (int i=0; i<n; i++)
    printf ("x[%d] = %f \n", i, gsl_vector_get(x,i));

  gsl_vector_free(e); e = gsl_vector_alloc(n);
  gsl_vector_free(f); f = gsl_vector_alloc(n);

  double alpha = 0.5;
  double beta  = 0.5;
  printf("\n *** A Matrix *** \n");
  for (int i = 0; i < n; i++)
    { for (int j = 0; j < n; j++)
        {
          printf ("%6.3f", gsl_matrix_get (A, i, j));
        }
      printf("\n");
    }

  for (int i = 0; i < n; i++)
    {
      gsl_vector_set(d, i, gsl_matrix_get(A,i,i));
      if (i>0)
        gsl_vector_set(f, i-1, gsl_matrix_get(A,i,i-1));
      if (i<n-1)
        gsl_vector_set(e, i, gsl_matrix_get(A,i,i+1));
    }

  gsl_vector_set(e,n-1,alpha);
  gsl_vector_set(f,n-1,beta);

  for (int i=0; i<n; i++) gsl_vector_set(b,i,B[i]);

  // A = ( d_0 e_0  0  f_3 )
  //     ( f_0 d_1 e_1  0  )
  //     (  0  f_1 d_2 e_2 )
  //     ( e_3  0  f_2 d_3 )
  for (int i=0; i<n; i++)
    printf ("%d = %f %f %f \n", i, gsl_vector_get(d,i)
                                 , gsl_vector_get(e,i)
                                 , gsl_vector_get(f,i));

  gsl_linalg_solve_cyc_tridiag (d, e, f, b, x);

  printf("\n cyclic \n\n");
  for (int i=0; i<n; i++)
    printf ("x[%d] = %f \n", i, gsl_vector_get(x,i));
}

void tridag(double a[],
            double b[],
            double c[],
            double r[],
            double u[],
            int n)
{
  double bet;
  double *gam;

  gam = (double *) malloc(n*sizeof(double));

  if (b[0] == 0.0) printf("Error 1 in tridag");

  u[0]=r[0]/(bet=b[0]);

  for (int j=1;j<n;j++)
  {
    gam[j]= c[j-1]/bet;
    bet=b[j]-a[j]*gam[j];
    if (bet == 0.0) printf("Error 2 in tridag");
    u[j]=(r[j]-a[j]*u[j-1])/bet;
  }

  for (int j=(n-2);j>=0;j--)
    u[j] -= gam[j+1]*u[j+1];

  free(gam);

}

void cyclic(double a[],
            double b[],
            double c[],
            double alpha,
            double beta,
            double r[],
            double x[], int n)

{
  double fact,gamma;
  double *bb;
  double *u;
  double *z;

  if (n <= 2) printf("n too small in cyclic");
  bb = (double *) malloc(n*sizeof(double));
  u  = (double *) malloc(n*sizeof(double));
  z  = (double *) malloc(n*sizeof(double));

  gamma   = -b[0];
  bb[0]   =  b[0]-gamma;
  bb[n-1] =  b[n-1]-alpha*beta/gamma;

  for (int i=1;i<n-1;i++)
    bb[i]=b[i];

  tridag(a,bb,c,r,x,n);
  u[0]  =gamma;
  u[n-1]=alpha;

  for (int i=1;i<n-1;i++)
    u[i]=0.0;

  tridag(a,bb,c,u,z,n);

  fact=(x[0]+beta*x[n-1]/gamma)/(1.0+z[0]+beta*z[n-1]/gamma);

  for (int i=0;i<n;i++)
    x[i] -= fact*z[i];

  free(z);
  free(u);
  free(bb);
}

void nr_solve( gsl_matrix *A,
               int n,
               double xmax,
               double xmin,
               double *B )
{

  printf("\n *** Solve with Numerical Recipes *** \n");

  double *a; a = (double *) malloc(n*sizeof(double));
  double *b; b = (double *) malloc(n*sizeof(double));
  double *c; c = (double *) malloc(n*sizeof(double));
  double *r; r = (double *) malloc(n*sizeof(double));
  double *u; u = (double *) malloc(n*sizeof(double));

  for (int i = 0; i < n; i++)
    {
      if (i>0) a[i] = gsl_matrix_get (A, i, i-1);
      b[i] = gsl_matrix_get (A, i, i);
      if (i<n-1) c[i] = gsl_matrix_get (A, i, i+1);
      r[i] = B[i];
    }

  for (int i=0; i<n; i++)
    printf ("%d = %f %f %f \n", i, a[i], b[i], c[i]);

  tridag ( a, b, c, r, u, n);
  for (int i=0; i<n; i++)
    printf ("u[%d] = %f \n", i, u[i]);

  double alpha = 0.5;
  double beta  = 0.5;

  cyclic( a, b, c, alpha, beta, r, u, n);

  printf("\n cyclic \n\n");
  for (int i=0; i<n; i++)
    printf ("u[%d] = %f \n", i, u[i]);
}

int
main (void)
{
  const size_t n = N;
  const size_t ncoeffs = NCOEFFS;
  const size_t nbreak = NBREAK;
  size_t i, j;
  gsl_bspline_workspace *bw;
  gsl_vector *B;
  double xmin=0.0, xmax=2.0;
  const size_t k = K;
  gsl_vector *tau;
  gsl_matrix *A;

  double *y;

  double dx = (xmax-xmin)/(n-1);
  double xc = (xmax+xmin)*0.5;
  y = (double *) malloc(ncoeffs*sizeof(double));

  tau = gsl_vector_alloc(n);

  /* allocate a cubic bspline workspace */
  bw = gsl_bspline_alloc(k, nbreak);
  B  = gsl_vector_alloc(ncoeffs);
  A  = gsl_matrix_alloc(n, ncoeffs);

  /* gsl_bspline_knots(breakpts, bw)
  Compute the knots from the given breakpoints:
   knots(1:k) = breakpts(1)
   knots(k+1:k+l-1) = breakpts(i), i = 2 .. l
   knots(n+1:n+k) = breakpts(l + 1) */

  /* use uniform breakpoints on [xmin, xmax] */
  gsl_bspline_knots_uniform(xmin, xmax, bw);
  //knots(1:k) = xmin
  //knots(k+1:k+l-1) = xmin + i*delta, i = 1 .. l - 1
  //knots(n+1:n+k) = xmax

  /* construct the fit matrix X */
  for (i = 0; i < n; ++i)
    {
      double taui = i * dx;
      double ti   = gsl_vector_get(bw->knots, i);

      gsl_vector_set(tau, i, taui);

      /* compute B_j(xi) for all j */
      gsl_bspline_eval(taui, B, bw);

      printf("%f %f\n", taui, ti);

      /* fill in row i of X */
      for (j = 0; j < ncoeffs; ++j)
      {
        double Bj = gsl_vector_get(B, j);
        gsl_matrix_set(A, i, j, Bj);
      }

    }

  for (int i = 0; i < n; i++)
    { for (int j = 0; j < n; j++)
        {
          printf ("%6.3f", gsl_matrix_get (A, i, j));
        }
      printf("\n");
    }

  for (int i=0; i<n;i++)
    y[i] = fabs(i*dx-xc);

  nr_solve( A, ncoeffs, xmax, xmin, y);

  gsl_solve( A, ncoeffs, xmax, xmin, y);

  gsl_bspline_free(bw);
  gsl_vector_free(B);
  gsl_matrix_free(A);
  return 0;

}

ext: c

max_stars_repo_name: utastudents/selalib

max_issues_repo_name: utastudents/selalib

max_forks_repo_name: utastudents/selalib


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "compearth.h"
#include "cmopad.h"
#ifdef COMPEARTH_USE_MKL
#ifdef __clang__
#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Wreserved-id-macro"
#pragma clang diagnostic ignored "-Wstrict-prototypes"
#endif
#include <mkl_cblas.h>
#ifdef __clang__
#pragma clang diagnostic pop
#endif
#else
#include <cblas.h>
#endif

struct mt_struct
{
    double mt[6];
    double tax[3];
    double bax[3];
    double pax[3];
    double sdr1[3]; //strike dip rake; fault plane 2
    double sdr2[3]; //strike dip rake; fault plane 1
    double mag;
    double exp;
};

/*
int compearth_standardDecomposition(const int nmt,
                                    const double *__restrict__ M,
                                    enum compearthCoordSystem_enum basis,
                                    double *__restrict__ M0,
                                    double *__restrict__ Mw,
                                    double *__restrict__ fp1,
                                    double *__restrict__ fp2,
                                    double *__restrict__ pAxis,
                                    double *__restrict__ bAxis,
                                    double *__restrict__ tAxis,
                                    double *__restrict__ isoPct,
                                    double *__restrict__ devPct,
                                    double *__restrict__ dcPct,
                                    double *__restrict__ clvdPct);
*/
int decompose(const int nmt, const double mtUSE[6]);

int decompose(const int nmt, const double mtUSE[6])
{
    //double K[3], N[3], S[3], lam[3], U[9], Uuse[9], b[3], p[3], t[3],
    //double gamma, delta, M0, kappa, theta, sigma, thetadc,
    //       k2, d2, s2;
    struct cmopad_struct src;
    const char *fcnm = "decompose\0";
    enum cmopad_basis_enum cin, cloc;
    double M[3][3], pax[3], bax[3], tax[3];
    int i, ierr;
    const int iverb = 0;
/*
    //------------------------------------------------------------------------//
    //                       Do the cMoPaD decomposition                      //
    //------------------------------------------------------------------------//
    // Copy moment tensor
    M[0][0] = mtUSE[0];           //Mrr
    M[1][1] = mtUSE[1];           //Mtt
    M[2][2] = mtUSE[2];           //Mpp
    M[0][1] = M[1][0] = mtUSE[3]; //Mrt
    M[0][2] = M[2][0] = mtUSE[4]; //Mrp
    M[1][2] = M[2][1] = mtUSE[5]; //Mtp
    // Mopad works in North, East, Down frame but inversion is Up, South, East
    cin = USE;
    cloc = NED;
    cmopad_basis_transformMatrixM33(M, cin, cloc); //USE -> NED
    // Compute the isotropic, CLVD, DC decomposition
    ierr = cmopad_standardDecomposition(M, &src);
    if (ierr != 0)
    {
        printf("%s: Error in decomposition!\n", fcnm);
        return EXIT_FAILURE;
    }
    // Compute the princple axes with corresponding strikes, dips, and rakes
    ierr = cmopad_MT2PrincipalAxisSystem(iverb, &src);
    if (ierr != 0)
    {
        printf("%s: Error computing principal axis\n",fcnm);
        return EXIT_FAILURE;
    }
    // Compute the pressure, null, and, tension principal axes as len,az,plunge
    cin = NED; //MoPaD is in north, east, down
    ierr = cmopad_Eigenvector2PrincipalAxis(cin, src.eig_pnt[0],
                                            src.p_axis,    pax);
    if (ierr != 0)
    {
        printf("%s: Error converting pax\n",fcnm);
        return EXIT_FAILURE;
    }
    ierr = cmopad_Eigenvector2PrincipalAxis(cin, src.eig_pnt[1],
                                            src.null_axis, bax);
    if (ierr != 0){
        printf("%s: Error converting bax\n",fcnm);
        return EXIT_FAILURE;
    }
    ierr = cmopad_Eigenvector2PrincipalAxis(cin, src.eig_pnt[2],
                                            src.t_axis,    tax);
    if (ierr != 0)
    {
        printf("%s: Error converting tax\n",fcnm);
        return EXIT_FAILURE;
    }
*/
    //printf("%f %f %f %f %f %f\n",
    //       mtUSE[0], mtUSE[1], mtUSE[2], mtUSE[3], mtUSE[4], mtUSE[5]);
    //ierr = compearth_CMT2TT(1, mtUSE, false,
    //                        &gamma, &delta, &M0,
    //                        &kappa, &theta, &sigma,
    //                        K, N, S, &thetadc,
    //                        lam, U);
    //if (ierr != 0)
    //{
    //    fprintf(stderr, "%s: CMT2TT failed\n", __func__);
    //    return EXIT_FAILURE;
    //}
    /*
    // Need to switch from SEU to USE.  This requires computing
    // R U inv(R) where:
    //  R = [0 0 1
    //       1 0 0
    //       0 1 0]
    // Carrying through the multiplication yields the permutation
    //  U_{use} = [U_{33} U_{31} U_{32}
    //             U_{13} U_{11} U_{12}
    //             U_{23} U_{21} U_{22}]
    //
    // Permute the column major matrix
    Uuse[0] = U[8];
    Uuse[1] = U[6];
    Uuse[2] = U[7];

    Uuse[3] = U[2];
    Uuse[4] = U[0];
    Uuse[5] = U[1];

    Uuse[6] = U[5];
    Uuse[7] = U[3];
    Uuse[8] = U[4];
printf("%f %f %f\n%f %f %f\n%f %f %f\n",
       U[0], U[3], U[6],
       U[1], U[4], U[7],
       U[2], U[5], U[8]);
*/
    //ierr = compearth_U2pa(1, U,
    //                      &p[2], &p[1], &b[2], &b[1], &t[2], &t[1]);
    //if (ierr != 0)
    //{
    //    fprintf(stderr, "%s: Error converting u to plunge/azimuth\n",
    //            __func__);
    //    return EXIT_FAILURE;
    //}
    //// Fill in the eigenvalues
    //p[0] = lam[0];
    //b[0] = lam[1];
    //t[0] = lam[2];
    //printf("p: %e %f %f\n", p[0], p[1], p[2]);
    //printf("b: %e %f %f\n", b[0], b[1], b[2]);
    //printf("t: %e %f %f\n", t[0] ,t[1], t[2]);
    // Compute the auxiliary fault plane
    //ierr = compearth_auxiliaryPlane(1, &kappa, &theta, &sigma,
    //                                &k2, &d2, &s2);
    // Just define some convention where smaller strike goes first
    //printf("%e %f %f %f\n", M0, kappa, theta, sigma);
    //printf("%f %f %f\n", k2, d2, s2);
    //------------------------------------------------------------------------//
    //                         do the standard decomposition                  //
    //------------------------------------------------------------------------//
    //double pAxis[3], bAxis[3], tAxis[3], fp1[3], fp2[3], isoPct, dcPct, devPct, clvdPct, Mw;
    double *pAxis, *bAxis, *tAxis, *fp1, *fp2, *isoPct, *dcPct, *devPct, *clvdPct, *Mw, *M0;
    M0 = (double *) calloc((size_t) nmt, sizeof(double));
    Mw = (double *) calloc((size_t) nmt, sizeof(double));
    fp1 = (double *) calloc((size_t) (3*nmt), sizeof(double));
    fp2 = (double *) calloc((size_t) (3*nmt), sizeof(double));
    pAxis = (double *) calloc((size_t) (3*nmt), sizeof(double));
    bAxis = (double *) calloc((size_t) (3*nmt), sizeof(double));
    tAxis = (double *) calloc((size_t) (3*nmt), sizeof(double));
    isoPct = (double *) calloc((size_t) nmt, sizeof(double));
    devPct = (double *) calloc((size_t) nmt, sizeof(double));
    dcPct = (double *) calloc((size_t) nmt, sizeof(double));
    clvdPct = (double *) calloc((size_t) nmt, sizeof(double));
    ierr = compearth_standardDecomposition(nmt, mtUSE, CE_USE,
                                           M0, Mw, fp1, fp2,
                                           pAxis, bAxis, tAxis,
                                           isoPct, devPct, dcPct, clvdPct);
    if (ierr != 0)
    {
        fprintf(stderr, "Failed to compute standard decomposition\n");
        return EXIT_FAILURE;
    }
    for (i=0; i<nmt; i++)
    {
        //--------------------------------------------------------------------//
        //                       Do the cMoPaD decomposition                  //
        //--------------------------------------------------------------------//
        // Copy moment tensor
        M[0][0] = mtUSE[6*i+0];           //Mrr
        M[1][1] = mtUSE[6*i+1];           //Mtt
        M[2][2] = mtUSE[6*i+2];           //Mpp
        M[0][1] = M[1][0] = mtUSE[6*i+3]; //Mrt
        M[0][2] = M[2][0] = mtUSE[6*i+4]; //Mrp
        M[1][2] = M[2][1] = mtUSE[6*i+5]; //Mtp
        // Mopad works in North, East, Down
        cin = USE;
        cloc = NED;
        cmopad_basis_transformMatrixM33(M, cin, cloc); //USE -> NED
        // Compute the isotropic, CLVD, DC decomposition
        ierr = cmopad_standardDecomposition(M, &src);
        if (ierr != 0)
        {
            printf("%s: Error in decomposition!\n", fcnm);
            return EXIT_FAILURE;
        }
        // Compute the princple axes and strikes, dips, and rakes
        ierr = cmopad_MT2PrincipalAxisSystem(iverb, &src);
        if (ierr != 0)
        {
            printf("%s: Error computing principal axis\n",fcnm);
            return EXIT_FAILURE;
        }
        // Compute the pressure, null, and, tension principal axes as
        // len,az,plunge
        cin = NED; //MoPaD is in north, east, down
        ierr = cmopad_Eigenvector2PrincipalAxis(cin, src.eig_pnt[0],
                                                src.p_axis,    pax);
        if (ierr != 0)
        {
            printf("%s: Error converting pax\n",fcnm);
            return EXIT_FAILURE;
        }
        ierr = cmopad_Eigenvector2PrincipalAxis(cin, src.eig_pnt[1],
                                                src.null_axis, bax);
        if (ierr != 0)
        {
            printf("%s: Error converting bax\n",fcnm);
            return EXIT_FAILURE;
        }
        ierr = cmopad_Eigenvector2PrincipalAxis(cin, src.eig_pnt[2],
                                                src.t_axis,    tax);
        if (ierr != 0)
        {
            printf("%s: Error converting tax\n",fcnm);
            return EXIT_FAILURE;
        }
        // Check it
        //printf("%d %d %e\n", i, i%9, Mw[i]);
        if (fabs(M0[i] - src.seismic_moment)/src.seismic_moment > 1.e-10)
        {
            fprintf(stderr, "failed M0 %f %f %e %d\n", M0[i], src.seismic_moment,
                    fabs(M0[i] - src.seismic_moment)/src.seismic_moment, i);
            return EXIT_FAILURE;
        }
        if (fabs(Mw[i] - src.moment_magnitude) > 1.e-10)
        {
            fprintf(stderr, "failed Mw %f %f %d\n", Mw[i], src.moment_magnitude, i);
            return EXIT_FAILURE;
        }
        if (fabs(isoPct[i] - src.ISO_percentage) > 1.e-10)
        {
            fprintf(stderr, "Failed iso pct %f %f\n",
                    isoPct[i], src.ISO_percentage);
            return EXIT_FAILURE;
        }
        if (fabs(dcPct[i] - src.DC_percentage) > 1.e-10)
        {
            fprintf(stderr, "Failed DC pct %f %f %d\n",
                    dcPct[i], src.DC_percentage, i);
            return EXIT_FAILURE;
        }
        if (fabs(devPct[i] - src.DEV_percentage) > 1.e-10)
        {
            fprintf(stderr, "Failed dev pct %f %f\n",
                    devPct[i], src.DEV_percentage);
            return EXIT_FAILURE;
        }
        if (fabs(clvdPct[i] - src.CLVD_percentage) > 1.e-10)
        {
            fprintf(stderr, "Failed clvd pct %f %f\n",
                    devPct[i], src.CLVD_percentage);
            return EXIT_FAILURE;
        }
        // Swap fault planes
        if (fabs(fp1[3*i+0] - src.fp1[0]) > 1.e-10)
        {
            if (fabs(fp1[3*i+0] - src.fp2[0]) > 1.e-10)
            {
                fprintf(stderr, "Failed strike 1 %f %f\n", fp1[3*i+0], src.fp2[0]);
                return EXIT_FAILURE;
            }
            if (fabs(fp1[3*i+1] - src.fp2[1]) > 1.e-10)
            {
                fprintf(stderr, "Failed dip 1 %f %f\n", fp1[3*i+1], src.fp2[1]);
                return EXIT_FAILURE;
            }
            if (fabs(fp1[3*i+2] - src.fp2[2]) > 1.e-10)
            {
               fprintf(stderr, "Failed rake 1 %f %f\n", fp1[3*i+2], src.fp2[2]);
               return EXIT_FAILURE;
            }
            if (fabs(fp2[3*i+0] - src.fp1[0]) > 1.e-10)
            {
                fprintf(stderr, "Failed strike 2 %f %f\n", fp2[3*i+0], src.fp1[0]);
                return EXIT_FAILURE;
            }
            if (fabs(fp2[3*i+1] - src.fp1[1]) > 1.e-10)
            {
                fprintf(stderr, "Failed dip 2 %f %f\n", fp2[3*i+1], src.fp1[1]);
                return EXIT_FAILURE;
            }
            if (fabs(fp2[3*i+2] - src.fp1[2]) > 1.e-10)
            {
                fprintf(stderr, "Failed rake 2 %f %f\n", fp2[3*i+2], src.fp1[2]);
                return EXIT_FAILURE;
            }
        }
        else
        {
            if (fabs(fp1[3*i+0] - src.fp1[0]) > 1.e-10)
            {
                fprintf(stderr, "Failed strike 1 %f %f\n",
                        fp1[3*i+0], src.fp1[0]);
                return EXIT_FAILURE;
            }
            if (fabs(fp1[3*i+1] - src.fp1[1]) > 1.e-10)
            {
                fprintf(stderr, "Failed dip 1 %f %f\n",
                        fp1[3*i+1], src.fp1[1]);
                return EXIT_FAILURE;
            }
            if (fabs(fp1[3*i+2] - src.fp1[2]) > 1.e-10)
            {
                fprintf(stderr, "Failed rake 1 %f %f\n",
                        fp1[3*i+2], src.fp1[2]);
                return EXIT_FAILURE;
            }
            if (fabs(fp2[3*i+0] - src.fp2[0]) > 1.e-10)
            {
                fprintf(stderr, "Failed strike 2 %f %f\n",
                       fp2[3*i+0], src.fp2[0]);
                return EXIT_FAILURE;
            }
            if (fabs(fp2[3*i+1] - src.fp2[1]) > 1.e-10)
            {
                fprintf(stderr, "Failed dip 2 %f %f\n",
                        fp2[3*i+1], src.fp2[1]);
                return EXIT_FAILURE;
            }
            if (fabs(fp2[3*i+2] - src.fp2[2]) > 1.e-10)
            {
                fprintf(stderr, "Failed rake 2 %f %f\n",
                        fp2[3*i+2], src.fp2[2]);
                return EXIT_FAILURE;
            }
        }
        if (fabs(pAxis[3*i+0] - pax[0]) > 1.e-10 ||
            fabs(pAxis[3*i+1] - pax[1]) > 1.e-10 ||
            fabs((pAxis[3*i+2] - pax[2])/pax[2]) > 1.e-10)
        {
            fprintf(stderr, "failed pressure axix: %f %f %f %f %f %f\n",
                    pAxis[0], pAxis[1], pAxis[2], pax[0], pax[1], pax[2]);
            return EXIT_FAILURE;
        }
        if (fabs(bAxis[3*i+0] - bax[0]) > 1.e-10 ||
            fabs(bAxis[3*i+1] - bax[1]) > 1.e-10 ||
            fabs((bAxis[3*i+2] - bax[2])/bax[2]) > 1.e-10)
        {
            fprintf(stderr, "Failed null axis: %f %f %f %f %f %f\n",
                    bAxis[0], bAxis[1], bAxis[2], bax[0], bax[1], bax[2]);
            return EXIT_FAILURE;
        }
        if (fabs(tAxis[3*i+0] - tax[0]) > 1.e-10 ||
            fabs(tAxis[3*i+1] - tax[1]) > 1.e-10 ||
            fabs((tAxis[3*i+2] - tax[2])/tax[2]) > 1.e-10)
        {
            fprintf(stderr, "Failed tension axis %f %f %f %f %f %f\n",
                    tAxis[0], tAxis[1], tAxis[2], tax[0], tax[1], tax[2]);
            return EXIT_FAILURE;
        }
        if (nmt == 1)
        {
            printf("Summary:\n");
            printf("Mw: %f\n", Mw[i]);
            printf("(strike,dip,rake): (%f,%f,%f)\n",
                   fp1[3*i+0], fp1[3*i+1], fp1[3*i+2]);
            printf("(strike,dip,rake): (%f,%f,%f)\n",
                   fp2[3*i+0], fp2[3*i+1], fp2[3*i+2]);
            printf("iso pct: %f\n", isoPct[i]);
            printf("dc pct: %f\n", dcPct[i]);
            printf("dev pct: %f\n", devPct[i]);
            printf("clvd pct: %f\n", clvdPct[i]);
            printf("\n");
        }
    }
    free(M0);
    free(Mw);
    free(fp1);
    free(fp2);
    free(pAxis);
    free(bAxis);
    free(tAxis);
    free(isoPct);
    free(devPct);
    free(dcPct);
    free(clvdPct);
    return EXIT_SUCCESS;
}

int main()
{
    struct mt_struct mt;
    const int nmt = 21*CE_CHUNKSIZE - 1;
    double *mtTest = (double *) calloc((size_t) (6*nmt), sizeof(double));
             //mtTest[2*6*CE_CHUNKSIZE];
    double xscal;
    int ierr, jmt;
    const int n6 = 6;
    const int incx = 1;
    //-----------------------------------1------------------------------------//
    mt.exp = 23.0 - 7.0;
    mt.mag = 5.1;
    mt.mt[0] =-4.360;
    mt.mt[1] = 3.200;
    mt.mt[2] = 1.170;
    mt.mt[3] =-0.794;
    mt.mt[4] = 1.650;
    mt.mt[5] = 2.060;
    mt.sdr1[0] = 77.0;
    mt.sdr1[1] = 47.0;
    mt.sdr1[2] =-62.0;
    mt.sdr2[0] = 219.0;
    mt.sdr2[1] = 50.0;
    mt.sdr2[2] =-117.0;
    mt.tax[0] = 4.49;
    mt.tax[1] = 1.0;
    mt.tax[2] = 328.0;
    mt.bax[0] = 0.56;
    mt.bax[1] = 20.0;
    mt.bax[2] = 237.0;
    mt.pax[0] =-5.04;
    mt.pax[1] = 70.0;
    mt.pax[2] = 61.0;
    xscal = pow(10.0, mt.exp);
    cblas_dscal(n6, xscal, mt.mt, incx);
    cblas_dcopy(n6, mt.mt, 1, &mtTest[0], 1);
    // Perform decomposition
    ierr = decompose(1, mt.mt);
    if (ierr != EXIT_SUCCESS)
    {
        fprintf(stderr, "Error decomposiing 1\n");
        return EXIT_FAILURE;
    }
    //--------------------------------------2----------------------------------//
    mt.exp = 23.0 - 7.0;
    mt.mag = 5.0;
    mt.mt[0] =-4.480;
    mt.mt[1] = 1.440;
    mt.mt[2] = 3.040;
    mt.mt[3] = 0.603;
    mt.mt[4] = 0.722;
    mt.mt[5] =-1.870;
    mt.sdr1[0] =316.0;
    mt.sdr1[1] = 44.0;
    mt.sdr1[2] =-105.0;
    mt.sdr2[0] = 157.0;
    mt.sdr2[1] = 48.0;
    mt.sdr2[2] =-76.0;
    mt.tax[0] = 4.28;
    mt.tax[1] = 2.0;
    mt.tax[2] = 237.0;
    mt.bax[0] = 0.37;
    mt.bax[1] = 10.0;
    mt.bax[2] = 327.0;
    mt.pax[0] =-4.66;
    mt.pax[1] = 79.0;
    mt.pax[2] =137.0;
    xscal = pow(10.0,mt.exp);
    cblas_dscal(n6, xscal, mt.mt, incx);
    cblas_dcopy(n6, mt.mt, 1, &mtTest[6], 1);
    // Perform decomposition
    ierr = decompose(1, mt.mt);
    if (ierr != EXIT_SUCCESS)
    {
        fprintf(stderr, "Error decomposiing 1\n");
        return EXIT_FAILURE;
    }
    //------------------------------------3-----------------------------------//
    mt.exp = 23.0 - 7.0;
    mt.mag = 4.9;
    mt.mt[0] =-2.460;
    mt.mt[1] = 0.207;
    mt.mt[2] = 2.250;
    mt.mt[3] = 0.793;
    mt.mt[4] = 0.267;
    mt.mt[5] =-0.363;
    mt.sdr1[0] =335.0;
    mt.sdr1[1] = 46.0;
    mt.sdr1[2] =-113.0;
    mt.sdr2[0] = 186.0;
    mt.sdr2[1] = 49.0;
    mt.sdr2[2] =-68.0;
    mt.tax[0] = 2.32;
    mt.tax[1] = 2.0;
    mt.tax[2] = 261.0;
    mt.bax[0] = 0.38;
    mt.bax[1] = 16.0;
    mt.bax[2] = 351.0;
    mt.pax[0] =-2.70;
    mt.pax[1] = 74.0;
    mt.pax[2] =165.0;
    xscal = pow(10.0,mt.exp);
    cblas_dscal(n6, xscal, mt.mt, incx);
    cblas_dcopy(n6, mt.mt, 1, &mtTest[12], 1);
    ierr = decompose(1, mt.mt);
    if (ierr != EXIT_SUCCESS)
    {
        fprintf(stderr, "Error decomposing 3\n");
        return EXIT_FAILURE;
    }
    //-----------------------------------4------------------------------------//
    mt.exp = 23.0 - 7.0;
    mt.mag = 4.9;
    mt.mt[0] =-2.270;
    mt.mt[1] = 0.058;
    mt.mt[2] = 2.210;
    mt.mt[3] = 0.079;
    mt.mt[4] =-0.737;
    mt.mt[5] = 0.246;
    mt.sdr1[0] =190.0;
    mt.sdr1[1] = 36.0;
    mt.sdr1[2] =-84.0;
    mt.sdr2[0] = 3.0;
    mt.sdr2[1] = 54.0;
    mt.sdr2[2] =-94.0;
    mt.tax[0] = 2.35;
    mt.tax[1] = 9.0;
    mt.tax[2] = 96.0;
    mt.bax[0] = 0.04;
    mt.bax[1] = 4.0;
    mt.bax[2] = 5.0;
    mt.pax[0] =-2.39;
    mt.pax[1] = 80.0;
    mt.pax[2] =253.0;
    xscal = pow(10.0,mt.exp);
    cblas_dscal(n6, xscal, mt.mt, incx);
    cblas_dcopy(n6, mt.mt, 1, &mtTest[18], 1);
    ierr = decompose(1, mt.mt);
    if (ierr != EXIT_SUCCESS)
    {
        fprintf(stderr, "Error dcomposing 4\n");
        return EXIT_FAILURE;
    }
    //-----------------------------------5------------------------------------//
    mt.exp = 26.0 - 7.0;
    mt.mag = 6.5;
    mt.mt[0] = 0.745;
    mt.mt[1] =-0.036;
    mt.mt[2] =-0.709;
    mt.mt[3] =-0.242;
    mt.mt[4] =-0.048;
    mt.mt[5] = 0.208;
    mt.sdr1[0] =181.0;
    mt.sdr1[1] = 47.0;
    mt.sdr1[2] =114.0;
    mt.sdr2[0] = 328.0;
    mt.sdr2[1] = 48.0;
    mt.sdr2[2] = 67.0;
    mt.tax[0] = 0.82;
    mt.tax[1] = 73.0;
    mt.tax[2] = 166.0;
    mt.bax[0] =-0.05;
    mt.bax[1] = 17.0;
    mt.bax[2] = 344.0;
    mt.pax[0] =-0.77;
    mt.pax[1] = 1.0;
    mt.pax[2] = 74.0;
    xscal = pow(10.0,mt.exp);
    cblas_dscal(n6, xscal, mt.mt, incx);
    cblas_dcopy(n6, mt.mt, 1, &mtTest[24], 1);
    ierr = decompose(1, mt.mt);
    if (ierr != EXIT_SUCCESS)
    {
        fprintf(stderr, "Error decomposing 5\n");
        return EXIT_FAILURE;
    }
    //-----------------------------------6------------------------------------//
    mt.exp = 24.0 - 7.0;
    mt.mag = 5.3;
    mt.mt[0] = 0.700;
    mt.mt[1] =-0.883;
    mt.mt[2] = 0.183;
    mt.mt[3] = 0.260;
    mt.mt[4] = 0.289;
    mt.mt[5] = 0.712;
    mt.sdr1[0] =329.0;
    mt.sdr1[1] = 54.0;
    mt.sdr1[2] =141.0;
    mt.sdr2[0] = 85.0;
    mt.sdr2[1] = 59.0;
    mt.sdr2[2] = 43.0;
    mt.tax[0] = 1.01;
    mt.tax[1] = 51.0;
    mt.tax[2] = 300.0;
    mt.bax[0] = 0.24;
    mt.bax[1] = 39.0;
    mt.bax[2] = 113.0;
    mt.pax[0] =-1.25;
    mt.pax[1] = 3.0;
    mt.pax[2] = 206.;
    xscal = pow(10.0,mt.exp);
    cblas_dscal(n6, xscal, mt.mt, incx);
    cblas_dcopy(n6, mt.mt, 1, &mtTest[30], 1);
    ierr = decompose(1, mt.mt);
    if (ierr != EXIT_SUCCESS)
    {
        fprintf(stderr, "Error decomposing 6\n");
        return EXIT_FAILURE;
    }
    //-----------------------------------7------------------------------------//
    mt.exp = 23.0 - 7.0;
    mt.mag = 5.0;
    mt.mt[0] = 3.150;
    mt.mt[1] =-2.470;
    mt.mt[2] =-0.676;
    mt.mt[3] = 1.650;
    mt.mt[4] = 1.880;
    mt.mt[5] =-1.470;
    mt.sdr1[0] =252.0;
    mt.sdr1[1] = 28.0;
    mt.sdr1[2] =111.0;
    mt.sdr2[0] = 49.0;
    mt.sdr2[1] = 64.0;
    mt.sdr2[2] = 79.0;
    mt.tax[0] = 4.08;
    mt.tax[1] = 69.0;
    mt.tax[2] = 297.0;
    mt.bax[0] = 0.01;
    mt.bax[1] = 10.0;
    mt.bax[2] = 54.0;
    mt.pax[0] =-4.08;
    mt.pax[1] = 18.0;
    mt.pax[2] = 147.0;
    xscal = pow(10.0,mt.exp);
    cblas_dscal(n6, xscal, mt.mt, incx);
    cblas_dcopy(n6, mt.mt, 1, &mtTest[36], 1);
    ierr = decompose(1, mt.mt);
    if (ierr != EXIT_SUCCESS)
    {
        fprintf(stderr, "Error decomposing 7\n");
        return EXIT_FAILURE;
    }
    //-----------------------------------8------------------------------------//
    mt.exp = 23.0 - 7.0;
    mt.mag = 4.8;
    mt.mt[0] =-1.870;
    mt.mt[1] = 0.488;
    mt.mt[2] = 1.380;
    mt.mt[3] =-0.057;
    mt.mt[4] = 0.600;
    mt.mt[5] = 0.664;
    mt.sdr1[0] = 36.0;
    mt.sdr1[1] = 38.0;
    mt.sdr1[2] =-76.0;
    mt.sdr2[0] = 199.0;
    mt.sdr2[1] = 53.0;
    mt.sdr2[2] =-101.0;
    mt.tax[0] = 1.80;
    mt.tax[1] = 8.0;
    mt.tax[2] = 296.0;
    mt.bax[0] = 0.18;
    mt.bax[1] = 9.0;
    mt.bax[2] = 205.0;
    mt.pax[0] =-1.99;
    mt.pax[1] = 78.0;
    mt.pax[2] = 69.0;
    xscal = pow(10.0,mt.exp);
    cblas_dscal(n6, xscal, mt.mt, incx);
    cblas_dcopy(n6, mt.mt, 1, &mtTest[42], 1);
    ierr = decompose(1, mt.mt);
    if (ierr != EXIT_SUCCESS)
    {
        fprintf(stderr, "Error decomposing 8\n");
        return EXIT_FAILURE;
    }
    //-----------------------------------9------------------------------------//
    mt.exp = 23.0 - 7.0;
    mt.mag = 5.0;
    mt.mt[0] =-3.540;
    mt.mt[1] = 1.040;
    mt.mt[2] = 2.510;
    mt.mt[3] =-0.474;
    mt.mt[4] = 1.640;
    mt.mt[5] = 1.530;
    mt.sdr1[0] = 47.0;
    mt.sdr1[1] = 38.0;
    mt.sdr1[2] =-63.0;
    mt.sdr2[0] = 195.0;
    mt.sdr2[1] = 56.0;
    mt.sdr2[2] =-109.0;
    mt.tax[0] = 3.66;
    mt.tax[1] = 10.0;
    mt.tax[2] = 299.0;
    mt.bax[0] = 0.45;
    mt.bax[1] = 16.0;
    mt.bax[2] =206.0;
    mt.pax[0] =-4.10;
    mt.pax[1] = 71.0;
    mt.pax[2] = 58.0;
    xscal = pow(10.0,mt.exp);
    cblas_dscal(n6, xscal, mt.mt, incx);
    cblas_dcopy(n6, mt.mt, 1, &mtTest[48], 1);
    // now copy this
    for (jmt=9; jmt<nmt; jmt++)
    {
        cblas_dcopy(n6, &mtTest[6*((jmt-9)%9)], 1, &mtTest[6*jmt], 1);
    }
    ierr = decompose(1, mt.mt);
    if (ierr != EXIT_SUCCESS)
    {
        fprintf(stderr, "Error decompsoing 9\n");
        return EXIT_FAILURE;
    }
    //-----------------------------Test 'em all-------------------------------//
    ierr = decompose(nmt, mtTest);
    if (ierr != EXIT_SUCCESS)
    {
        fprintf(stderr, "Error in bulk decomposition\n");
        return EXIT_FAILURE;
    }
    fprintf(stdout, "Success!\n");
    free(mtTest);
    return EXIT_SUCCESS;
}

max_stars_repo_name: OUCyf/mtbeach

max_issues_repo_name: carltape/mtbeach

max_forks_repo_name: carltape/mtbeach


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <pthread.h>

#include <gsl/gsl_sf.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>

#include "SIShift.h"

static int  verbose  = FALSE;

static int  nLines = 10;

static char *usage[] = {"SIShift - \n",
                        "Required parameters:\n",
			"   -a    abundance data\n",
			"   -in   filename      parameter file  \n",
                        "Optional:\n",
			"   -b    integer       length of sample file to ignore\n",
			"   -s    integer       sampling frequency\n",
			"   -seed long          seed random number generator\n",
			"   -v                  verbose\n"};

void updateExpectations(double* adExpect, int nMax, double dAlpha, double dBeta, double dGamma, int nS, t_Data *ptData)
{
  int i = 0;

  for(i = 1; i <= nMax; i++){
    int nA = i;
    double dLog = 0.0, dP = 0.0;

    dLog = logLikelihood(nA, dAlpha, dBeta, dGamma);

    dP = exp(dLog);

    adExpect[i - 1] += dP*nS;
  }
}

int main(int argc, char* argv[]){
  t_Params    tParams;
  t_SIParams* atSIParams;
  int         i = 0, nSamples = 0, nMax = 0;
  t_Data      tData;
  double*     adExpect = NULL;
  double      dShift = 0.0;

  gsl_set_error_handler_off();

  getCommandLineParams(&tParams, argc, argv);

  /*allocate memory for samples*/
  atSIParams = (t_SIParams *) malloc(MAX_SAMPLES*sizeof(t_SIParams));
  if(!atSIParams)
    goto memoryError;

  /*read in Monte-Carlo samples*/
  readSamples(&tParams, atSIParams, &nSamples);

  readAbundanceData(tParams.szAbundFile, &tData);

  nMax = 1000;//tData.aanAbund[tData.nNA - 1][0];
  nMax = floor(pow(2.0,ceil(log((double) nMax)/log(2.0)) + 2.0) + 1.0e-7);

  adExpect = (double *) malloc(sizeof(double)*nMax);

  for(i = 0; i < nMax; i++){
    adExpect[i] = 0.0;
  }

  dShift = 5.0e5/(double) tData.nJ;
  for(i = 0; i < nSamples; i++){
    updateExpectations(adExpect, nMax,
		       atSIParams[i].dAlpha*sqrt(dShift),
		       atSIParams[i].dBeta*dShift,
		       atSIParams[i].dGamma,
		       atSIParams[i].nS,
		       &tData);
  }

  for(i = 1; i <= nMax; i++){
    printf("%d %f\n",i,adExpect[i - 1]/((double) nSamples));
  }

  free(adExpect);
  exit(EXIT_SUCCESS);

 memoryError:
  fprintf(stderr, "Failed to allocate memory in main aborting ...\n");
  fflush(stderr);
  exit(EXIT_FAILURE);
}

int intCompare(const void *pvA, const void *pvB)
{
  int* pnA = (int *) pvA, *pnB = (int *) pvB;

  if(*pnA < *pnB)
    return +1;
  else if(*pnA == *pnB){
    return 0;
  }
  else{
    return -1;
  }
}

int doubleCompare(const void *pvA, const void *pvB)
{
  double* pnA = (double *) pvA, *pnB = (double *) pvB;

  if(*pnA < *pnB)
    return +1;
  else if(*pnA == *pnB){
    return 0;
  }
  else{
    return -1;
  }
}

void writeUsage(FILE* ofp)
{
  int i = 0;
  char *line;

  for(i = 0; i < nLines; i++){
    line = usage[i];
    fputs(line,ofp);
  }
}

char *extractParameter(int argc, char **argv, char *param,int when)
{
  int i = 0;

  while((i < argc) && (strcmp(param,argv[i]))){
    i++;
  }

  if(i < argc - 1){
    return(argv[i + 1]);
  }

  if((i == argc - 1) && (when == OPTION)){
    return "";
  }

  if(when == ALWAYS){
    fprintf(stdout,"Can't find asked option %s\n",param);
  }

  return (char *) NULL;
}

void getCommandLineParams(t_Params *ptParams,int argc,char *argv[])
{
  char *szTemp = NULL;
  char *cError = NULL;

  /*get parameter file name*/
  ptParams->szInputFile  = extractParameter(argc,argv, INPUT_FILE,ALWAYS);
  if(ptParams->szInputFile == NULL)
    goto error;

 /*get abundance file name*/
  ptParams->szAbundFile  = extractParameter(argc,argv, ABUND_FILE,ALWAYS);
  if(ptParams->szAbundFile == NULL)
    goto error;

  /*get long seed*/
  szTemp  = extractParameter(argc,argv,SEED,OPTION);
  if(szTemp != NULL){
    ptParams->lSeed = strtol(szTemp,&cError,10);
    if(*cError != '\0'){
      goto error;
    }
  }
  else{
    ptParams->lSeed = 0;
  }

  /*get burn*/
  szTemp  = extractParameter(argc, argv, BURN, OPTION);
  if(szTemp != NULL){
    ptParams->nBurn = strtol(szTemp,&cError,10);
    if(*cError != '\0'){
      goto error;
    }
  }
  else{
    ptParams->nBurn = DEF_BURN;
  }

 /*get long seed*/
  szTemp  = extractParameter(argc, argv, SAMPLE, OPTION);
  if(szTemp != NULL){
    ptParams->nSample = strtol(szTemp,&cError,10);
    if(*cError != '\0'){
      goto error;
    }
  }
  else{
    ptParams->nSample = DEF_SAMPLE;
  }

  /*verbosity*/
  szTemp = extractParameter(argc, argv, VERBOSE, OPTION);
  if(szTemp != NULL){
    verbose = TRUE;
  }

  return;

 error:
  writeUsage(stdout);
  exit(EXIT_FAILURE);
}

void readSamples(t_Params *ptParams, t_SIParams *atSIParams, int *pnSamples)
{
  int  nSamples = 0;
  char *szInputFile = ptParams->szInputFile;
  char szLine[MAX_LINE_LENGTH];
  FILE *ifp = NULL;

  ifp = fopen(szInputFile, "r");

  if(ifp){
    while(fgets(szLine, MAX_LINE_LENGTH, ifp)){
      char *szTok = NULL, *szBrk = NULL, *pcError = NULL;
      int  nTime  = 0;

      /*remove trailing new line*/
      szBrk = strpbrk(szLine, "\n"); (*szBrk) = '\0';

      szTok = strtok(szLine, DELIM);

      nTime = strtol(szTok, &pcError, 10);
      if(*pcError != '\0') goto fileFormatError;

      if(nTime > ptParams->nBurn && nTime % ptParams->nSample == 0){

	szTok = strtok(NULL,DELIM);
	atSIParams[nSamples].dAlpha = strtod(szTok, &pcError);
	if(*pcError != '\0') goto fileFormatError;

	szTok = strtok(NULL,DELIM);
	atSIParams[nSamples].dBeta = strtod(szTok, &pcError);
	if(*pcError != '\0') goto fileFormatError;

	szTok = strtok(NULL,DELIM);
	atSIParams[nSamples].dGamma = strtod(szTok, &pcError);
	if(*pcError != '\0') goto fileFormatError;

	szTok = strtok(NULL,DELIM);
	atSIParams[nSamples].nS = strtol(szTok, &pcError, 10);
	if(*pcError != '\0') goto fileFormatError;

	nSamples++;
      }
    }

    fclose(ifp);
  }
  else{
    fprintf(stderr, "Failed to open file %s aborting\n", szInputFile);
    fflush(stderr);
    exit(EXIT_FAILURE);
  }

  (*pnSamples) = nSamples;
  return;

 fileFormatError:
  fprintf(stderr, "Incorrectly formatted input file aborting\n");
  fflush(stderr);
  exit(EXIT_FAILURE);
}

void readAbundanceData(const char *szFile, t_Data *ptData)
{
  int **aanAbund = NULL;
  int  i = 0, nNA = 0, nA = 0, nC = 0;
  int  nL = 0, nJ = 0;
  char szLine[MAX_LINE_LENGTH];
  FILE* ifp = NULL;

  ifp = fopen(szFile, "r");

  if(ifp){
    char* szTok   = NULL;
    char* pcError = NULL;

    fgets(szLine, MAX_LINE_LENGTH, ifp);

    szTok = strtok(szLine, DELIM2);

    nNA = strtol(szTok,&pcError,10);
    if(*pcError != '\0'){
      goto formatError;
    }

    aanAbund = (int **) malloc(nNA*sizeof(int*));

    for(i = 0; i < nNA; i++){
      aanAbund[i] = (int *) malloc(sizeof(int)*2);

      fgets(szLine, MAX_LINE_LENGTH, ifp);

      szTok = strtok(szLine, DELIM2);

      nA = strtol(szTok,&pcError,10);
      if(*pcError != '\0'){
	goto formatError;
      }

      szTok = strtok(NULL, DELIM2);

      nC = strtol(szTok,&pcError,10);
      if(*pcError != '\0'){
	goto formatError;
      }

      nL += nC;
      nJ += nC*nA;

      aanAbund[i][0]  = nA;
      aanAbund[i][1]  = nC;
    }
  }
  else{
    fprintf(stderr, "Failed to open abundance data file %s aborting\n", szFile);
    fflush(stderr);
    exit(EXIT_FAILURE);
  }

  ptData->nJ          = nJ;
  ptData->nL          = nL;
  ptData->aanAbund    = aanAbund;
  ptData->nNA         = nNA;
  return;

 formatError:
  fprintf(stderr, "Incorrectly formatted abundance data file\n");
  fflush(stderr);
  exit(EXIT_FAILURE);
}

double fX(double x, double dA, double dB, double dNDash)
{
  double dTemp1 = (dA*(x - dB)*(x - dB))/x;

  return log(x) - (1.0/dNDash)*(x + dTemp1);
}

double f2X(double x, double dA, double dB, double dNDash)
{
  double dRet = 0.0, dTemp = 2.0*dA*dB*dB;

  dRet = (1.0/(x*x))*(1.0 + (1.0/dNDash)*(dTemp/x));

  return -dRet;
}

double sd(int n, double dAlpha, double dBeta, double dGamma)
{
 double dA = 0.5*(-1.0 + sqrt(1.0 + (dAlpha*dAlpha)/(dBeta*dBeta)));
  double dN = (double) n, dNDash = dN + dGamma - 1.0, dRN = 1.0/dN;
  double dTemp1 = (0.5*dN)/(1.0 + dA), dTemp2 = 4.0*dRN*dRN*(1.0 + dA)*dA*dBeta*dBeta;
  double dXStar = dTemp1*(1.0 + sqrt(1.0 + dTemp2));
  double dFX = fX(dXStar, dA, dBeta, dNDash);
  double d2FX = -dNDash*f2X(dXStar, dA, dBeta, dNDash);
  double dLogK = 0.0, dGamma1 = dGamma;

  if(dGamma1 < 0.0){
    dGamma1 *= -1.0;
  }

  dLogK = gsl_sf_bessel_lnKnu(dGamma1,2.0*dA*dBeta);

  return -2.0*dA*dBeta -log(2.0) -dLogK -dGamma*log(dBeta) + dNDash*dFX + 0.5*log(2.0*M_PI) - 0.5*log(d2FX);
}

int bessel(double* pdResult, int n, double dAlpha, double dBeta, double dGamma)
{
  double dResult = 0.0;
  double dOmega = 0.0, dGamma2 = 0.0;
  double dLogK1 = 0.0, dLogK2 = 0.0;
  double dN = (double) n, dNu = dGamma + dN;
  double dTemp1 = 0.0;

  if(dNu < 0.0){
    dNu = -dNu;
  }

  if(dGamma < 0.0){
    dGamma2 = -dGamma;
  }
  else{
    dGamma2 = dGamma;
  }

  dOmega = sqrt(dBeta*dBeta + dAlpha*dAlpha) - dBeta;

  dLogK2 = gsl_sf_bessel_lnKnu(dNu, dAlpha);

  if(!gsl_finite(dLogK2)){
    if(dAlpha < 0.1*sqrt(dNu + 1.0)){
      //printf("l ");
      dLogK2 = gsl_sf_lngamma(dNu) + (dNu - 1.0)*log(2.0) - dNu*log(dAlpha);
    }
    else{
      //printf("sd ");
      (*pdResult) = dResult;
      return FALSE;
    }
  }

  dLogK1 = dGamma*log(dOmega/dAlpha) -gsl_sf_bessel_lnKnu(dGamma2,dOmega);

  dTemp1 = log((dBeta*dOmega)/dAlpha);

  dResult = dN*dTemp1 + dLogK2 + dLogK1;
  (*pdResult) = dResult;
  return TRUE;
}

double logLikelihood(int n, double dAlpha, double dBeta, double dGamma)
{
  double dLogFacN = 0.0;
  int status      = 0;
  double dRet     = 0.0;

  if(n < 50){
    dLogFacN = gsl_sf_fact(n);
    dLogFacN = log(dLogFacN);
  }
  else{
    dLogFacN = gsl_sf_lngamma(((double) n) + 1.0);
  }

  status = bessel(&dRet,n, dAlpha,dBeta,dGamma);
  if(status == FALSE){
    dRet = sd(n, dAlpha,dBeta,dGamma);
  }

  return dRet - dLogFacN;
}

max_stars_repo_name: chrisquince/DiversityEstimates

max_issues_repo_name: chrisquince/DiversityEstimates

max_forks_repo_name: chrisquince/DiversityEstimates


text:
/* -*- Mode: C; tab-width: 4; c-basic-offset: 4; indent-tabs-mode: nil -*- */
#pragma once

#include <mcbp/protocol/status.h>
#include <cstdint>
#include <functional>
#include <gsl/gsl>

struct EngineIface;

/**
 * Callback for any function producing stats.
 *
 * @param key the stat's key
 * @param value the stat's value in an ascii form (e.g. text form of a number)
 * @param cookie magic callback cookie
 */
using AddStatFn = std::function<void(std::string_view key,
                                     std::string_view value,
                                     gsl::not_null<const void*> cookie)>;

/**
 * Callback for adding a response backet
 * @param key The key to put in the response
 * @param extras The data to put in the extended field in the response
 * @param body The data body
 * @param datatype This is currently not used and should be set to 0
 * @param status The status code of the return packet (see in protocol_binary
 *               for the legal values)
 * @param cas The cas to put in the return packet
 * @param cookie The cookie provided by the frontend
 * @return true if return message was successfully created, false if an
 *              error occured that prevented the message from being sent
 */
using AddResponseFn = std::function<bool(std::string_view key,
                                         std::string_view extras,
                                         std::string_view body,
                                         uint8_t datatype,
                                         cb::mcbp::Status status,
                                         uint64_t cas,
                                         const void* cookie)>;

ext: h


text:
//        Copyright The Authors 2018.
//    Distributed under the 3-Clause BSD License.
//    (See accompanying file LICENSE or copy at
//   https://opensource.org/licenses/BSD-3-Clause)

#pragma once

#include <cstdint>      // for std::uint8_t etc.
#include <type_traits>  // for std::underlying_type_t<>

#include <gsl/gsl>  // for gsl::span

#include <p2p/kademlia/buffer.h>
#include <p2p/kademlia/node.h>

namespace blocxxi {
namespace p2p {
namespace kademlia {

/*!
 * @brief Represents a message header block.
 *
 * In its serialized form, every message is composed of a header and a body. The
 * header contains the protocol version, the message body type, the ID (hash
 * 160) of the source node and a random token (hash 160) used to correlate
 * requests and responses.
 *
 * A serialized header uses the following layout:
 * ```
 *      <-----------------><-------...-------><-------...------->
 *            1 byte             20 bytes           20 bytes
 *        4 bits | 4 bits      source id          random token
 *        version  type
 *
 * ```
 */
struct Header final {
  /// Protocol version
  enum class Version : std::uint8_t {
    V1 = 1,
  } version_{Version::V1};

  /// Message body type
  enum class MessageType : std::uint8_t {
    PING_REQUEST,
    PING_RESPONSE,
    STORE_REQUEST,
    FIND_NODE_REQUEST,
    FIND_NODE_RESPONSE,
    FIND_VALUE_REQUEST,
    FIND_VALUE_RESPONSE,
  } type_{MessageType::PING_REQUEST};

  /// ID of the source node (Hash 160)
  Node::IdType source_id_;
  /// A random token (Hash 160)
  blocxxi::crypto::Hash160 random_token_;
};

/*!
 * @brief Generalized function that can take any scoped enumerator and return
 * a compile time constant of the underlying type.
 *
 * Used specifically for the serialization and deserialization of headers. For
 * more details on the rational see Scott Meyers, Effective Modern C++ Item 10.
 *
 * @tparam TEnum the scoped enum class type.
 * @param [in] e_value any value of the scoped enum class.
 * @return a compile time constant of the underlying type corresponding to
 * e_value.
 */
template <typename TEnum>
constexpr std::underlying_type_t<TEnum> ToUnderlying(TEnum e_value) noexcept {
  return static_cast<std::underlying_type_t<TEnum>>(e_value);
}

/*!
 * @brief A template trait class used to carry extra information about the
 * messages types in use.
 *
 * The traits class T<M> lets us record such extra information about a message
 * class M, without requiring any change at all to M. Such extra information is
 * used to generically serialize and deserialize messages.
 *
 * @tparam TMessageBody the type of the message body described by this traits
 * class.
 */
template <typename TMessageBody>
struct MessageTraits;

/*!
 * @brief Serialize the header into the given buffer, which is expected to be
 * able to accomodate the header binary size (41 bytes) or expand as needed.
 *
 * @param [in] header header object to be serialized.
 * @param [out] buffer destination buffer for the header serialized binary form.
 */
void Serialize(Header const &header, Buffer &buffer);

/*!
 * @brief Deserialize the contents of the buffer into a Header object.
 *
 * @param [in] buffer a span over a range of bytes (std::uint8_t) representing
 * the header binary data.
 * @param [out] header the deserialized resulting header object.
 * @return an error code indicating success or failure of the deserialization.
 */
std::size_t Deserialize(BufferReader const &buffer, Header &header);

/// FIND_NODE request message body.
struct FindNodeRequestBody final {
  /// The hash 160 id of the node to find.
  Node::IdType node_id_;
};

/// Message traits template specialization for the FIND_NODE request message.
template <>
struct MessageTraits<FindNodeRequestBody> {
  static constexpr Header::MessageType TYPE_ID =
      Header::MessageType::FIND_NODE_REQUEST;
};

/*!
 * @brief Serialize a FIND_NODE request body into the given buffer.
 *
 * @param [in] body the request body.
 * @param [out] buffer the destination buffer.
 */
void Serialize(FindNodeRequestBody const &body, Buffer &buffer);

/*!
 * @brief Deserialize a FIND_NODE request body from the given buffer.
 *
 * @param [in] buffer the input buffer.
 * @param [out] body the resulting request body.
 * @return the number of consumed bytes from the input buffer. Subsequent
 * deserialization from the buffer need to start after the consumed bytes.
 */
std::size_t Deserialize(BufferReader const &buffer, FindNodeRequestBody &body);

/// FIND_NODE response message body.
struct FindNodeResponseBody final {
  /// The list of nodes in the response.
  std::vector<Node> peers_;
};

/// Message traits template specialization for the FIND_NODE response message.
template <>
struct MessageTraits<FindNodeResponseBody> {
  static constexpr Header::MessageType TYPE_ID =
      Header::MessageType::FIND_NODE_RESPONSE;
};

/*!
 * @brief Serialize a FIND_NODE response body into the given buffer.
 *
 * @param [in] body the response body.
 * @param [out] buffer the destination buffer.
 */
void Serialize(FindNodeResponseBody const &body, Buffer &buffer);

/*!
 * @brief Deserialize a FIND_NODE response body from the given buffer.
 *
 * @param [in] buffer the input buffer.
 * @param [out] body the resulting response body.
 * @return the number of consumed bytes from the input buffer. Subsequent
 * deserialization from the buffer need to start after the consumed bytes.
 */
std::size_t Deserialize(BufferReader const &buffer, FindNodeResponseBody &body);

/// FIND_VALUE request message body.
struct FindValueRequestBody final {
  /// The hash 160 key of the value to find.
  blocxxi::crypto::Hash160 value_key_;
};

/// Message traits template specialization for the FIND_VALUE request message.
template <>
struct MessageTraits<FindValueRequestBody> {
  static constexpr Header::MessageType TYPE_ID =
      Header::MessageType::FIND_VALUE_REQUEST;
};

/*!
 * @brief Serialize a FIND_VALUE request body into the given buffer.
 *
 * @param [in] body the request body.
 * @param [out] buffer the destination buffer.
 */
void Serialize(FindValueRequestBody const &body, Buffer &buffer);

/*!
 * @brief Deserialize a FIND_VALUE request body from the given buffer.
 *
 * @param [in] buffer the input buffer.
 * @param [out] body the resulting request body.
 * @return the number of consumed bytes from the input buffer. Subsequent
 * deserialization from the buffer need to start after the consumed bytes.
 */
std::size_t Deserialize(BufferReader const &buffer, FindValueRequestBody &body);

/// FIND_VALUE response message body.
struct FindValueResponseBody final {
  /// The value data  as a vector of bytes.
  std::vector<std::uint8_t> data_;
};

/// Message traits template specialization for the FIND_VALUE response message.
template <>
struct MessageTraits<FindValueResponseBody> {
  static constexpr Header::MessageType TYPE_ID =
      Header::MessageType::FIND_VALUE_RESPONSE;
};

/*!
 * @brief Serialize a FIND_VALUE response body into the given buffer.
 *
 * @param [in] body the response body.
 * @param [out] buffer the destination buffer.
 */
void Serialize(FindValueResponseBody const &body, Buffer &buffer);

/*!
 * @brief Deserialize a FIND_VALUE response body from the given buffer.
 *
 * @param [in] buffer the input buffer.
 * @param [out] body the resulting response body.
 * @return the number of consumed bytes from the input buffer. Subsequent
 * deserialization from the buffer need to start after the consumed bytes.
 */
std::size_t Deserialize(BufferReader const &buffer,
                        FindValueResponseBody &body);

/// STORE_VALUE request message body.
struct StoreValueRequestBody final {
  /// The hash 160 of the data key.
  blocxxi::crypto::Hash160 data_key_;
  /// The data value as a vector of bytes.
  std::vector<std::uint8_t> data_value_;
};

/// Message traits template specialization for the STORE_VALUE request message.
template <>
struct MessageTraits<StoreValueRequestBody> {
  static constexpr Header::MessageType TYPE_ID =
      Header::MessageType::STORE_REQUEST;
};

/*!
 * @brief Serialize a STORE_VALUE request body into the given buffer.
 *
 * @param [in] body the request body.
 * @param [out] buffer the destination buffer.
 */
void Serialize(StoreValueRequestBody const &body, Buffer &buffer);

/*!
 * @brief Deserialize a STORE_VALUE request body from the given buffer.
 *
 * @param [in] buffer the input buffer.
 * @param [out] body the resulting request body.
 * @return the number of consumed bytes from the input buffer. Subsequent
 * deserialization from the buffer need to start after the consumed bytes.
 */
std::size_t Deserialize(BufferReader const &buffer,
                        StoreValueRequestBody &body);

}  // namespace kademlia
}  // namespace p2p
}  // namespace blocxxi


text:
// mathematical routines

#include "../include/bjutils.h"

double integrate_qag(double (*func)(double, void*),void *par,double low,double up,double *error)
{
  const int lim=2000;
  const int key=5;
  const double eabs=0.0;
  const double erel=1.e-4;
  double res,err;

  gsl_integration_workspace *work=gsl_integration_workspace_alloc(lim);

  gsl_function F;
  F.function=func;
  F.params=par;

  gsl_integration_qag(&F,low,up,eabs,erel,lim,key,work,&res,&err);

  gsl_integration_workspace_free(work);
  if (error!=NULL) *error=err;
  return(res);
}

double findroot(double (*func)(double, void*),void *par,double low,double hig,double tol,double *error)
{
  #include <gsl/gsl_errno.h>
  int status,iter=0;
  double root,rnext;
  const int max_iter=1000;    // maximum no. of iterations

  gsl_function F;
  F.function=func;
  F.params=par;

  const gsl_root_fsolver_type *T=gsl_root_fsolver_brent;
  gsl_root_fsolver *solv=gsl_root_fsolver_alloc(T);
  gsl_root_fsolver_set(solv,&F,low,hig);

  do {
    iter++;
    status=gsl_root_fsolver_iterate(solv);
    root=gsl_root_fsolver_root(solv);
    low=gsl_root_fsolver_x_lower(solv);
    hig=gsl_root_fsolver_x_upper(solv);
    status=gsl_root_test_interval(low,hig,0.0,tol);
    if (status==GSL_SUCCESS) printf ("Root finder converged after %i iterations\n",iter);
  } while (status==GSL_CONTINUE && iter<max_iter);

  gsl_root_fsolver_iterate(solv);
  rnext=gsl_root_fsolver_root(solv);
  *error=fabs(root-rnext);
  gsl_root_fsolver_free(solv);
  return(root);
}

// index search via bisection (log N process)
// array, dimension of array, value to be found, returned: index of array lower than x
void bisect_interpol(double *a,int dim,double x,int *erg)
{
  int l=0,u=dim-1,m;

  if ((x<a[l]) || (x>a[u])) {
    printf("Error: index out of range: %g not in %g-%g\n",x,a[l],a[u]);
    exit(-1);
  }
  while (u-l>1) {
    m=(u+l)/2;
    if (x>a[m]) l=m;
    else u=m;
  }
  *erg=l;
}

// 1D linear interpolation + extrapolation
double interpol_linear_extra_1D(double *ax,int NX,double *z,double x)
{
  if (NX==1) return(z[0]);  // assume z constant if only one data point
  if (x<=ax[0]) return((z[1]-z[0])/(ax[1]-ax[0])*(x-ax[0])+z[0]);
  else if (x>=ax[NX-1]) return((z[NX-1]-z[NX-2])/(ax[NX-1]-ax[NX-2])*(x-ax[NX-1])+z[NX-1]);
  else {
    int ix;
    bisect_interpol(ax,NX,x,&ix);
    double t=(x-ax[ix])/(ax[ix+1]-ax[ix]);
    return((1-t)*z[ix]+t*z[ix+1]);
  }
}

// 2D linear interpolation + extrapolation
double interpol_linear_extra_2D(double *ax,int NX,double *ay,int NY,double **z,double x,double y)
{
  if ((x<=ax[0])&&(y>ay[0])&&(y<ay[NY-1])) {
    double int1=interpol_linear_extra_1D(ay,NY,z[0],y);
    double int2=interpol_linear_extra_1D(ay,NY,z[1],y);
    return(int1+(int2-int1)/(ax[1]-ax[0])*(x-ax[0]));
  }
  else if ((x>=ax[NX-1])&&(y>ay[0])&&(y<ay[NY-1])) {
    double int1=interpol_linear_extra_1D(ay,NY,z[NX-2],y);
    double int2=interpol_linear_extra_1D(ay,NY,z[NX-1],y);
    return(int1+(int2-int1)/(ax[NX-1]-ax[NX-2])*(x-ax[NX-2]));
  }
  else if ((y<=ay[0])&&(x>ax[0])&&(x<ax[NX-1])) {
    int i;
    double tmp1[NX],tmp2[NX];
    for(i=0;i<NX;i++) {
      tmp1[i]=z[i][0];
      tmp2[i]=z[i][1];
    }
    double int1=interpol_linear_extra_1D(ax,NX,tmp1,x);
    double int2=interpol_linear_extra_1D(ax,NX,tmp2,x);
    return(int1+(int2-int1)/(ay[1]-ay[0])*(y-ay[0]));
  }
  else if ((y>=ay[NY-1])&&(x>ax[0])&&(x<ax[NX-1])) {
    int i;
    double tmp1[NX],tmp2[NX];
    for(i=0;i<NX;i++) {
      tmp1[i]=z[i][NY-2];
      tmp2[i]=z[i][NY-1];
    }
    double int1=interpol_linear_extra_1D(ax,NX,tmp1,x);
    double int2=interpol_linear_extra_1D(ax,NX,tmp2,x);
    return(int1+(int2-int1)/(ay[NY-1]-ay[NY-2])*(y-ay[NY-2]));
  }
  else {  // produces error via 'bisect' if both dimensions outside range
    int ix,iy;
    bisect_interpol(ax,NX,x,&ix);
    bisect_interpol(ay,NY,y,&iy);
    double t=(x-ax[ix])/(ax[ix+1]-ax[ix]);
    double u=(y-ay[iy])/(ay[iy+1]-ay[iy]);
    return((1-t)*(1-u)*z[ix][iy]+t*(1-u)*z[ix+1][iy]+(1-t)*u*z[ix][iy+1]+t*u*z[ix+1][iy+1]);
  }
}

// 3D linear interpolation + extrapolation
double interpol_linear_extra_3D(double *ax,int NX,double *ay,int NY,double *az,int NZ,double ***f,double x,double y,double z)
{
  if ((x<=ax[0])&&(y>ay[0])&&(y<ay[NY-1])&&(z>az[0])&&(z<az[NZ-1])) {   // x out of range, low
    double int1=interpol_linear_extra_2D(ay,NY,az,NZ,f[0],y,z);
    double int2=interpol_linear_extra_2D(ay,NY,az,NZ,f[1],y,z);
    return(int1+(int2-int1)/(ax[1]-ax[0])*(x-ax[0]));
  }
  else if ((x>=ax[NX-1])&&(y>ay[0])&&(y<ay[NY-1])&&(z>az[0])&&(z<az[NZ-1])) {   // x out of range, high
    double int1=interpol_linear_extra_2D(ay,NY,az,NZ,f[NX-2],y,z);
    double int2=interpol_linear_extra_2D(ay,NY,az,NZ,f[NX-1],y,z);
    return(int1+(int2-int1)/(ax[NX-1]-ax[NX-2])*(x-ax[NX-2]));
  }
  else if ((y<=ay[0])&&(x>ax[0])&&(x<ax[NX-1])&&(z>az[0])&&(z<az[NZ-1])) {   // y out of range, low
    int i,j;
    double **tmp1=malloc(NX*sizeof(double *));
    double **tmp2=malloc(NX*sizeof(double *));
    for(i=0;i<NX;i++) {
      tmp1[i]=malloc(NZ*sizeof(double));
      tmp2[i]=malloc(NZ*sizeof(double));
    }
    for(i=0;i<NX;i++) {
      for(j=0;j<NZ;j++) {
	tmp1[i][j]=f[i][0][j];
	tmp2[i][j]=f[i][1][j];
      }
    }
    double int1=interpol_linear_extra_2D(ax,NX,az,NZ,tmp1,x,z);
    double int2=interpol_linear_extra_2D(ax,NX,az,NZ,tmp2,x,z);
    for(i=0;i<NX;i++) {
      free(tmp1[i]);
      free(tmp2[i]);
    }
    free(tmp1);
    free(tmp2);
    return(int1+(int2-int1)/(ay[1]-ay[0])*(y-ay[0]));
  }
  else if ((y>=ay[NY-1])&&(x>ax[0])&&(x<ax[NX-1])&&(z>az[0])&&(z<az[NZ-1])) {   // y out of range, high
    int i,j;
    double **tmp1=malloc(NX*sizeof(double *));
    double **tmp2=malloc(NX*sizeof(double *));
    for(i=0;i<NX;i++) {
      tmp1[i]=malloc(NZ*sizeof(double));
      tmp2[i]=malloc(NZ*sizeof(double));
    }
    for(i=0;i<NX;i++) {
      for(j=0;j<NZ;j++) {
	tmp1[i][j]=f[i][NY-2][j];
	tmp2[i][j]=f[i][NY-1][j];
      }
    }
    double int1=interpol_linear_extra_2D(ax,NX,az,NZ,tmp1,x,z);
    double int2=interpol_linear_extra_2D(ax,NX,az,NZ,tmp2,x,z);
    for(i=0;i<NX;i++) {
      free(tmp1[i]);
      free(tmp2[i]);
    }
    free(tmp1);
    free(tmp2);
    return(int1+(int2-int1)/(ay[NY-1]-ay[NY-2])*(y-ay[NY-2]));
  }
  else if ((z<=az[0])&&(x>ax[0])&&(x<ax[NX-1])&&(y>ay[0])&&(y<ay[NY-1])) {   // z out of range, low
    int i,j;
    double **tmp1=malloc(NX*sizeof(double *));
    double **tmp2=malloc(NX*sizeof(double *));
    for(i=0;i<NX;i++) {
      tmp1[i]=malloc(NY*sizeof(double));
      tmp2[i]=malloc(NY*sizeof(double));
    }
    for(i=0;i<NX;i++) {
      for(j=0;j<NY;j++) {
	tmp1[i][j]=f[i][j][0];
	tmp2[i][j]=f[i][j][1];
      }
    }
    double int1=interpol_linear_extra_2D(ax,NX,ay,NY,tmp1,x,y);
    double int2=interpol_linear_extra_2D(ax,NX,ay,NY,tmp2,x,y);
    for(i=0;i<NX;i++) {
      free(tmp1[i]);
      free(tmp2[i]);
    }
    free(tmp1);
    free(tmp2);
    return(int1+(int2-int1)/(az[1]-az[0])*(z-az[0]));
  }
  else if ((z>=az[NZ-1])&&(x>ax[0])&&(x<ax[NX-1])&&(y>ay[0])&&(y<ay[NY-1])) {   // z out of range, high
    int i,j;
    double **tmp1=malloc(NX*sizeof(double *));
    double **tmp2=malloc(NX*sizeof(double *));
    for(i=0;i<NX;i++) {
      tmp1[i]=malloc(NY*sizeof(double));
      tmp2[i]=malloc(NY*sizeof(double));
    }
    for(i=0;i<NX;i++) {
      for(j=0;j<NY;j++) {
	tmp1[i][j]=f[i][j][NZ-2];
	tmp2[i][j]=f[i][j][NZ-1];
      }
    }
    double int1=interpol_linear_extra_2D(ax,NX,ay,NY,tmp1,x,y);
    double int2=interpol_linear_extra_2D(ax,NX,ay,NY,tmp2,x,y);
    for(i=0;i<NX;i++) {
      free(tmp1[i]);
      free(tmp2[i]);
    }
    free(tmp1);
    free(tmp2);
    return(int1+(int2-int1)/(az[NZ-1]-az[NZ-2])*(z-az[NZ-2]));
  }
  else {  // produces error via 'bisect' if more than 1 dimension outside range
    int ix,iy,iz;
    double f1a,f1b,f2a,f2b;

    bisect_interpol(ax,NX,x,&ix);
    bisect_interpol(ay,NY,y,&iy);
    bisect_interpol(az,NZ,z,&iz);

    double q1=(x-ax[ix])/(ax[ix+1]-ax[ix]);
    double q2=(y-ay[iy])/(ay[iy+1]-ay[iy]);
    double q3=(z-az[iz])/(az[iz+1]-az[iz]);
    double p1=1-q1;
    double p2=1-q2;

    f1a=p1*f[ix][iy][iz]+q1*f[ix+1][iy][iz];
    f1b=p1*f[ix][iy+1][iz]+q1*f[ix+1][iy+1][iz];
    f2a=p2*f1a+q2*f1b;

    f1a=p1*f[ix][iy][iz+1]+q1*f[ix+1][iy][iz+1];
    f1b=p1*f[ix][iy+1][iz+1]+q1*f[ix+1][iy+1][iz+1];
    f2b=p2*f1a+q2*f1b;

    return((1-q3)*f2a+q3*f2b);
  }
}

// 4D linear interpolation + extrapolation
double interpol_linear_extra_4D(double *ax,int NX,double *ay,int NY,double *az,int NZ,double *aw,int NW,double ****f,double x,double y,double z,double w)
{
  int ix;
  double f1,f2,t;

  if (((x<ax[0])||(x>ax[NX-1]))&&((y<ay[0])||(y>ay[NY-1])||(z<az[0])||(z>az[NZ-1])||(w<aw[0])||(w>aw[NW-1]))) {
    printf("Error in routine 'interpol_linear_extra_4D': attempt extrapolation in more than one dimension.\n");
    exit(-1);
  }

  if (x<=ax[0]) {
    f1=interpol_linear_extra_3D(ay,NY,az,NZ,aw,NW,f[0],y,z,w);
    f2=interpol_linear_extra_3D(ay,NY,az,NZ,aw,NW,f[1],y,z,w);
    return((f2-f1)/(ax[1]-ax[0])*(x-ax[0])+f1);
  }
  else if (x>=ax[NX-1]) {
    f1=interpol_linear_extra_3D(ay,NY,az,NZ,aw,NW,f[NX-2],y,z,w);
    f2=interpol_linear_extra_3D(ay,NY,az,NZ,aw,NW,f[NX-1],y,z,w);
    return((f2-f1)/(ax[NX-1]-ax[NX-2])*(x-ax[NX-1])+f2);
  }
  else {
    bisect_interpol(ax,NX,x,&ix);

    f1=interpol_linear_extra_3D(ay,NY,az,NZ,aw,NW,f[ix],y,z,w);
    f2=interpol_linear_extra_3D(ay,NY,az,NZ,aw,NW,f[ix+1],y,z,w);

    t=(x-ax[ix])/(ax[ix+1]-ax[ix]);
    return((1-t)*f1+t*f2);
  }
}

void mattimesmat_gen(gsl_matrix *a,gsl_matrix *b,gsl_matrix *erg,int dim1,int dim2,int dim3)
{
  int i,j,k;
  double sum;
  for (i=0;i<dim1;i++) {
    for (j=0;j<dim3;j++) {
      sum=0.0;
      for (k=0;k<dim2;k++) {
        sum+=gsl_matrix_get(a,i,k)*gsl_matrix_get(b,k,j);
      }
      gsl_matrix_set(erg,i,j,sum);
    }
  }
}

void mattimesmat(gsl_matrix *a,gsl_matrix *b,gsl_matrix *erg,int N)
{
  int i,j,k;
  double sum;
  for (i=0;i<N;i++) {
    for (j=0;j<N;j++) {
      sum=0.0;
      for (k=0;k<N;k++) {
        sum+=gsl_matrix_get(a,i,k)*gsl_matrix_get(b,k,j);
      }
      gsl_matrix_set(erg,i,j,sum);
    }
  }
}

void mattimesvec_gen(gsl_matrix *a,gsl_vector *b,gsl_vector *erg,int dim1,int dim2)
{
  int i,j;
  double sum;
  for (i=0;i<dim1;i++) {
    sum=0.0;
    for (j=0;j<dim2;j++) {
      sum+=gsl_matrix_get(a,i,j)*gsl_vector_get(b,j);
    }
    gsl_vector_set(erg,i,sum);
  }
}

void mattimesvec(gsl_matrix *a,gsl_vector *b,gsl_vector *erg,int N)
{
  int i,j;
  double sum;
  for (i=0;i<N;i++) {
    sum=0.0;
    for (j=0;j<N;j++) {
      sum+=gsl_matrix_get(a,i,j)*gsl_vector_get(b,j);
    }
    gsl_vector_set(erg,i,sum);
  }
}

int matrixinvers(gsl_matrix *m,gsl_matrix *inv,int DIM)
{
  int i,j,k,flag=0;
  double elem,sum;
  const int DOCHECK=1;

  gsl_matrix *v=gsl_matrix_alloc(DIM,DIM);
  gsl_matrix *id=gsl_matrix_alloc(DIM,DIM);
  gsl_matrix *u=gsl_matrix_alloc(DIM,DIM);
  gsl_vector *work=gsl_vector_alloc(DIM);
  gsl_vector *s=gsl_vector_alloc(DIM);

  gsl_matrix_memcpy(u,m);
  gsl_linalg_SV_decomp(u,v,s,work);

  for (i=0;i<DIM;i++) {
    elem=gsl_vector_get(s,i);
    if (fabs(elem)==0) {
      printf("Error: zero singular value!\n");
      exit(-1);
    }
    if (fabs(elem)<1.e-35) {
      printf("Possible singular value: %g\n",elem);
      flag=1;
    }
    //printf("eigenvalue: %g\n",elem);
  }

  for (i=0;i<DIM;i++) {
    for (j=0;j<DIM;j++) {
      sum=0.0;
      for (k=0;k<DIM;k++) {
	sum+=gsl_matrix_get(v,i,k)/gsl_vector_get(s,k)*gsl_matrix_get(u,j,k);
      }
      gsl_matrix_set(inv,i,j,sum);
    }
  }

  if (DOCHECK) {
    mattimesmat(m,inv,v,DIM);           //v reused
    gsl_matrix_set_identity(id);
    gsl_matrix_sub(v,id);                //result in v
    for (i=0;i<DIM;i++) {
      for (j=0;j<DIM;j++) {
	elem=gsl_matrix_get(v,i,j);
	if (fabs(elem)>1.e-2) {
	  printf("Problematic matrix inversion: %g\n",elem);
	  flag=2;
	  break;
	}
      }
    }
  }

  gsl_matrix_free(u);
  gsl_matrix_free(v);
  gsl_matrix_free(id);
  gsl_vector_free(work);
  gsl_vector_free(s);
  return(flag);
}

gsl_matrix *pseudoinvers(gsl_matrix *m,double thresh,double *eigen)
{
  int i,j,k;
  double sum;

  const int dim=(int)m->size1;
  if (m->size1!=m->size2) {
    printf("Error in 'pseudoinvers': non-square matrix - not implemented yet.\n");
    exit(-1);
  }

  gsl_matrix *res=gsl_matrix_alloc(dim,dim);
  gsl_matrix *v=gsl_matrix_alloc(dim,dim);
  gsl_matrix *cov=gsl_matrix_alloc(dim,dim);
  gsl_vector *work=gsl_vector_alloc(dim);
  gsl_vector *s=gsl_vector_alloc(dim);
  gsl_vector *invs=gsl_vector_alloc(dim);

  gsl_matrix_memcpy(cov,m);
  gsl_linalg_SV_decomp(cov,v,s,work);

  for (i=0;i<dim;i++) {
    eigen[i]=gsl_vector_get(s,i);
    if (eigen[i]<=thresh) gsl_vector_set(invs,i,0.0);
    else gsl_vector_set(invs,i,1./eigen[i]);
  }

  for (i=0;i<dim;i++) {
    for (j=0;j<dim;j++) {
      sum=0.0;
      for (k=0;k<dim;k++) {
	sum+=gsl_matrix_get(v,i,k)*gsl_vector_get(invs,k)*gsl_matrix_get(cov,j,k);
      }
      gsl_matrix_set(res,i,j,sum);
    }
  }

  gsl_matrix_free(cov);
  gsl_matrix_free(v);
  gsl_vector_free(work);
  gsl_vector_free(s);
  gsl_vector_free(invs);
  return(res);
}

int covinvers(gsl_matrix *m,gsl_matrix *inv,int DIM)
{
  #include<gsl/gsl_errno.h>
  gsl_set_error_handler_off();
  const int DOCHECK=1;

  int i,flag=0;
  gsl_matrix *check=gsl_matrix_alloc(DIM,DIM);
  gsl_vector *b=gsl_vector_alloc(DIM);

  gsl_matrix_memcpy(check,m);

  flag=gsl_linalg_cholesky_decomp(check);

  if (!flag) {
    for (i=0;i<DIM;i++) {
      gsl_vector_set_basis(b,i);
      gsl_linalg_cholesky_svx(check,b);
      gsl_matrix_set_col(inv,i,b);
    }

    if (DOCHECK) {
      gsl_matrix *id=gsl_matrix_alloc(DIM,DIM);
      gsl_matrix_set_identity(id);
      int j;
      double elem;

      mattimesmat(m,inv,check,DIM);           //check reused
      gsl_matrix_sub(check,id);               //result in check
      for (i=0;i<DIM;i++) {
	for (j=0;j<DIM;j++) {
	  elem=gsl_matrix_get(check,i,j);
	  if (fabs(elem)>1.e-6) {
	    printf("Problematic matrix inversion: %g\n",elem);
	    flag=2;
	    break;
	  }
	}
      }
      gsl_matrix_free(id);
    }
  }

  gsl_matrix_free(check);
  gsl_vector_free(b);
  return(flag);
}

double matrix_trace(gsl_matrix *a,int N)
{
  int i;
  double sum=0.0;
  for (i=0;i<N;i++) {
    sum+=gsl_matrix_get(a,i,i);
  }
  return(sum);
}

double matrix_det(gsl_matrix *a,int N)
{
  int sign;
  gsl_permutation *perm=gsl_permutation_alloc(N);
  gsl_matrix *LU=gsl_matrix_alloc(N,N);
  gsl_matrix_memcpy(LU,a);
  gsl_linalg_LU_decomp(LU,perm,&sign);
  gsl_permutation_free(perm);
  return(gsl_linalg_LU_det(LU,sign));
}

void matrix_transpose(gsl_matrix *a,gsl_matrix *t,int dim1,int dim2)
{
  int i,j;
  for (i=0;i<dim1;i++) {
    for (j=0;j<dim2;j++) {
      gsl_matrix_set(t,j,i,gsl_matrix_get(a,i,j));
    }
  }
  return;
}

void matrix_print(gsl_matrix *a,char *file)
{
  int i,j;
  FILE *dat;
  if ((dat=fopen(file,"w"))==NULL) {
    printf("Error in routine 'matrix_print': could not open file %s\n",file);
    exit(-1);
  }
  for (i=0;i<a->size1;i++) {
    for (j=0;j<a->size2;j++) {
      fprintf(dat,"%15.10g\t",gsl_matrix_get(a,i,j));
    }
    fprintf(dat,"\n");
  }
  fclose(dat);
  return;
}

void vector_print(gsl_vector *v,char *file)
{
  int i;
  FILE *dat;
  if ((dat=fopen(file,"w"))==NULL) {
    printf("Error in routine 'vector_print': could not open file %s\n",file);
    exit(-1);
  }
  for (i=0;i<v->size;i++) {
      fprintf(dat,"%15.10g\n",gsl_vector_get(v,i));
  }
  fclose(dat);
  return;
}

void solve_lse(gsl_matrix *m,gsl_vector *in,gsl_vector *res)
{
  if (m->size1!=res->size) {
    printf("Error in routine 'solve_lse': matrix row dimension does not match dimension of solution vector: %zu - %zu\n",m->size1,res->size);
    exit(-1);
  }
  if (m->size2!=in->size) {
    printf("Error in routine 'solve_lse': matrix column dimension does not match dimension of input vector: %zu - %zu\n",m->size2,in->size);
    exit(-1);
  }

  gsl_matrix *v=gsl_matrix_alloc(m->size2,m->size2);
  gsl_matrix *u=gsl_matrix_alloc(m->size1,m->size2);
  gsl_vector *work=gsl_vector_alloc(m->size2);
  gsl_vector *s=gsl_vector_alloc(m->size2);

  gsl_matrix_memcpy(u,m);
  gsl_linalg_SV_decomp(u,v,s,work);
  gsl_linalg_SV_solve(u,v,s,in,res);

  gsl_matrix_free(u);
  gsl_matrix_free(v);
  gsl_vector_free(work);
  gsl_vector_free(s);
  return;
}

// inverse of symmetric matrix via Cholesky decomposition
int inverse_symm(gsl_matrix *in,gsl_matrix *out)
{
  #include<gsl/gsl_errno.h>
  gsl_set_error_handler_off();

  int i,flag=0;
  const int dim=(int)in->size1;
  gsl_matrix *check=gsl_matrix_alloc(dim,dim);
  gsl_vector *b=gsl_vector_alloc(dim);
  gsl_matrix_memcpy(check,in);  // do not overwrite 'in'

  flag=gsl_linalg_cholesky_decomp(check);

  if (!flag) {
    for (i=0;i<dim;i++) {
      gsl_vector_set_basis(b,i);
      gsl_linalg_cholesky_svx(check,b);
      gsl_matrix_set_col(out,i,b);
    }
  }
  gsl_matrix_free(check);
  gsl_vector_free(b);
  return(flag);
}

int fisher_ellipse(gsl_matrix *fisher,int index1,int index2,double *semi_major,double *semi_minor,double *angle)
{
  int flag;

  // setup
  if (index1==index2) return(2);  // can only make ellipse for two different indices
  const int dim=fisher->size1;
  gsl_matrix *inv=gsl_matrix_alloc(dim,dim);
  gsl_matrix *sub=gsl_matrix_alloc(2,2);
  gsl_matrix *v=gsl_matrix_alloc(2,2);
  gsl_vector *work=gsl_vector_alloc(2);
  gsl_vector *s=gsl_vector_alloc(2);

  // invert Fisher matrix
  flag=inverse_symm(fisher,inv);

  // get ellipse plotting parameters
  gsl_matrix_set(sub,0,0,gsl_matrix_get(inv,index1,index1));
  gsl_matrix_set(sub,0,1,gsl_matrix_get(inv,index1,index2));
  gsl_matrix_set(sub,1,0,gsl_matrix_get(inv,index2,index1));
  gsl_matrix_set(sub,1,1,gsl_matrix_get(inv,index2,index2));

  gsl_linalg_SV_decomp(sub,v,s,work);
  *semi_major=sqrt(gsl_vector_get(s,0));
  *semi_minor=sqrt(gsl_vector_get(s,1));
  *angle=atan(gsl_matrix_get(v,1,0)/gsl_matrix_get(v,0,0));

  // clean up
  gsl_vector_free(work);
  gsl_matrix_free(v);
  gsl_matrix_free(inv);
  gsl_matrix_free(sub);
  gsl_vector_free(s);
  return(flag);
}

// provides KDE in 1D with Gaussians, using Silverman's rule if width<=0 on call
void kde_1D(double *dp,int nd,int nstep,double step_min,double step_max,double h,double *steps,double *func)
{
  int i,j;
  double width=1.;
  const double gauss_norm=1./sqrt(2.*PI);

  // set binning
  for (i=0;i<nstep;i++) {
    steps[i]=step_min+i*(step_max-step_min)/(1.*nstep-1.);
  }

  // get width of Gaussian
  if (h<=0.0) {  // apply Silverman's rule
    width=1.06*gsl_stats_sd(dp,1,nd)*pow(nd,-0.2);
  }
  else width=h;

  // perform KDE
  for (i=0;i<nstep;i++) {
    func[i]=0.0;
    for (j=0;j<nd;j++) {
      func[i]+=exp(-0.5*(steps[i]-dp[j])*(steps[i]-dp[j])/(width*width))/width;
    }
    func[i]*=gauss_norm;  // normalised to nd if integrated
  }

  return;
}


text:
/* -----------------------------------------------------------------------------
 * Copyright 2021 Jonathan Haigh
 * SPDX-License-Identifier: MIT
 * ---------------------------------------------------------------------------*/

#ifndef SQ_INCLUDE_GUARD_system_linux_SqPrimitiveTypeSchemaImpl_h_
#define SQ_INCLUDE_GUARD_system_linux_SqPrimitiveTypeSchemaImpl_h_

#include "core/typeutil.h"
#include "system/SqPrimitiveTypeSchema.gen.h"
#include "system/schema.h"

#include <gsl/gsl>

namespace sq::system::linux {

class SqPrimitiveTypeSchemaImpl
    : public SqPrimitiveTypeSchema<SqPrimitiveTypeSchemaImpl> {
public:
  explicit SqPrimitiveTypeSchemaImpl(
      const PrimitiveTypeSchema &primitive_type_schema);

  SQ_ND Result get_name() const;
  SQ_ND Result get_doc() const;

  SQ_ND Primitive to_primitive() const override;

private:
  gsl::not_null<const PrimitiveTypeSchema *> primitive_type_schema_;
};

} // namespace sq::system::linux

#endif // SQ_INCLUDE_GUARD_system_linux_SqPrimitiveTypeSchemaImpl_h_

max_stars_repo_name: jonathanhaigh/sq

max_issues_repo_name: jonathanhaigh/sq

max_forks_repo_name: jonathanhaigh/sq