text:
static const char help[] =
"Solves doubly-nonlinear obstacle problems in 2D.  Option prefix dnl_.\n"
"The PDE (interior condition) of such problems has solution u(x,y):\n"
"       - div (u^q |grad(u+b)|^{p-2} grad(u+b)) = f\n"
"subject to a obstacle constraint\n"
"       u >= psi\n"
"FIXME Includes the steady-state, nonlinear ice sheet problem in 2D in which u=H\n"
"is ice thickness,  b  is bed elevation, and  s = H + b  is surface elevation:\n"
"       - div (D grad H) - div(W H^{n+2}) = m\n"
"The diffusivity D and pseudo-velocity W (Bueler, 2016) are from the\n"
"nonsliding shallow ice approximation (SIA) flux:\n"
"       D = Gamma H^{n+2} |grad H + grad b|^{n-1}\n"
"       W = - Gamma |grad H + grad b|^{n-1} grad b\n"
"The climatic mass balance f = m(x,y,H) is from one of two models.\n"
"Constants are  n >= 1  and Gamma = 2 A (rho g)^n / (n+2)  where A is the ice\n"
"softness.  The domain is square  (0,L)^2  with zero Dirichlet boundary conditions.\n"
"The equation is discretized by a Q1 structured-grid FVE method (Bueler, 2016).\n"
"Requires SNESVI (-snes_type vinewton{rsls|ssls}) because of constraint;\n"
"defaults to SSLS.\n\n";

/*
1. shows basic success with SSLS but DIVERGES AT LEVEL 4:
   mpiexec -n 4 ./ice -ice_verif -snes_converged_reason -snes_grid_sequence LEV

2. consider making CMB model smooth

3. add CMB to dump and create plotting script (.py)

4. using exact init shows convergence depends strongly on eps for fine grids:
    for LEV in 1 2 3 4 5; do ./ice -ice_verif -ice_exact_init -snes_converged_reason -ksp_type gmres -pc_type gamg -da_refine $LEV -ice_eps EPS; done
result:
  (a) works at all levels if EPS=0.005; last KSP somewhat constant but SNES iters growing
  (b) fails on level 3 if EPS=0.003,0.002

5. convergent and nearly optimal GMG in flops *but cheating with exact init*, and *avoiding -snes_grid_sequence* and *significant eps=0.01 regularization*:
    for LEV in 1 2 3 4 5 6 7 8; do ./ice -ice_verif -ice_exact_init -snes_converged_reason -ksp_type gmres -pc_type mg -da_refine $LEV -snes_type vinewtonrsls -ice_eps 0.01; done

6. visualizing -snes_grid_sequence:
    ./ice -ice_verif -snes_grid_sequence 2 -ice_eps 0.005 -snes_converged_reason -snes_monitor_solution draw
(was -snes_grid_sequence bug with periodic BCs? see PETSc issue #300)

8. even seems to work in parallel:
    mpiexec -n 4 ./ice -ice_verif -snes_grid_sequence 5 -ice_eps 0.005 -snes_converged_reason -snes_monitor_solution draw

9. same outcome with -ice_exact_init and -da_refine 5
    mpiexec -n 4 ./ice -ice_verif -da_refine 5 -ice_eps 0.005 -snes_converged_reason -snes_monitor_solution draw -ice_exact_init

10. unpredictable response to changing -snes_linesearch_type bt|l2|basic  (cp seems rarely to work)
*/

/* see comments on runtime stuff in icet/icet.c, the time-dependent version */

#include <petsc.h>
#include "icecmb.h"

typedef struct {
    double    secpera,    // number of seconds in a year
              L,          // spatial domain is (0,L) x (0,L)
              g,          // acceleration of gravity
              rho_ice,    // ice density
              n_ice,      // Glen exponent for SIA flux term
              A_ice,      // ice softness
              Gamma,      // coefficient for SIA flux term
              D0,         // representative value of diffusivity (used in regularizing D)
              eps,        // regularization parameter for diffusivity D
              delta,      // dimensionless regularization for slope in SIA formulas
              lambda;     // amount of upwinding; lambda=0 is none and lambda=1 is "full"
    PetscBool verif,      // use dome formulas if true
              check_admissible; // check admissibility at start of FormFunctionLocal()
    CMBModel  *cmb;       // defined in cmbmodel.h
} AppCtx;

// compute radius from center of (0,L) x (0,L)
double radialcoord(double x, double y, AppCtx *user) {
  const double xc = x - user->L/2.0,
               yc = y - user->L/2.0;
  return PetscSqrtReal(xc * xc + yc * yc);
}

double DomeCMB(double x, double y, AppCtx *user) {
  const double  domeR  = 750.0e3,  // radius of exact ice sheet (m)
                domeH0 = 3600.0,   // center thickness of exact ice sheet (m)
                n  = user->n_ice,
                pp = 1.0 / n,
                CC = user->Gamma * PetscPowReal(domeH0,2.0*n+2.0)
                         / PetscPowReal(2.0 * domeR * (1.0-1.0/n),n);
  double        r, s, tmp1, tmp2;
  r = radialcoord(x, y, user);
  // avoid singularities at center and margin
  if (r < 0.01)
      r = 0.01;
  if (r > domeR - 0.01)
      r = domeR - 0.01;
  s = r / domeR;
  tmp1 = PetscPowReal(s,pp) + PetscPowReal(1.0-s,pp) - 1.0;
  tmp2 = 2.0 * PetscPowReal(s,pp) + PetscPowReal(1.0-s,pp-1.0) * (1.0 - 2.0*s) - 1.0;
  return (CC / r) * PetscPowReal(tmp1,n-1.0) * tmp2;
}

PetscErrorCode DomeThicknessLocal(DMDALocalInfo *info, double **aH, AppCtx *user) {
  const double   domeR  = 750.0e3,  // radius of exact ice sheet (m)
                 domeH0 = 3600.0,   // center thickness of exact ice sheet (m)
                 n  = user->n_ice,
                 mm = 1.0 + 1.0 / n,
                 qq = n / (2.0 * n + 2.0),
                 CC = domeH0 / PetscPowReal(1.0 - 1.0 / n,qq),
                 dx = user->L / (double)(info->mx-1),
                 dy = user->L / (double)(info->my-1);
  double         x, y, r, s, tmp;
  int            j, k;

  PetscFunctionBeginUser;
  for (k=info->ys; k<info->ys+info->ym; k++) {
      y = k * dy;
      for (j=info->xs; j<info->xs+info->xm; j++) {
          x = j * dx;
          r = radialcoord(x, y, user);
          // avoid singularities at margin and center
          if (r > domeR - 0.01)
              aH[k][j] = 0.0;
          else {
              if (r < 0.01)
                  r = 0.01;
              s = r / domeR;
              tmp = mm * s - (1.0/n) + PetscPowReal(1.0-s,mm) - PetscPowReal(s,mm);
              aH[k][j] = CC * PetscPowReal(tmp,qq);
          }
      }
  }
  PetscFunctionReturn(0);
}

extern PetscErrorCode SetFromOptionsAppCtx(AppCtx*);
extern PetscErrorCode FormBedLocal(DMDALocalInfo*, int, double**, AppCtx*);
extern PetscErrorCode FormBounds(SNES,Vec,Vec);
extern PetscErrorCode FormFunctionLocal(DMDALocalInfo*, double**, double**, AppCtx*);

int main(int argc,char **argv) {
  PetscErrorCode      ierr;
  DM                  da;
  SNES                snes;
  KSP                 ksp;
  Vec                 H;
  AppCtx              user;
  CMBModel            cmb;
  PetscBool           exact_init = PETSC_FALSE, // initialize using dome exact solution
                      dump = PETSC_FALSE;       // dump state (H,b) in binary file ice_MXxMY.dat after solve
  DMDALocalInfo       info;
  double              **aH;
  SNESConvergedReason reason;
  int                 snesit,kspit;

  PetscInitialize(&argc,&argv,(char*)0,help);

  user.secpera    = 31556926.0;  // number of seconds in a year
  user.L          = 1800.0e3;    // m; compare domeR=750.0e3 radius
  user.g          = 9.81;        // m/s^2
  user.rho_ice    = 910.0;       // kg/m^3
  user.n_ice      = 3.0;         // Glen exponent
  user.A_ice      = 3.1689e-24;  // 1/(Pa^3 s); EISMINT I value
  user.D0         = 1.0;         // m^2 / s
  user.eps        = 0.001;
  user.delta      = 1.0e-4;
  user.lambda     = 0.25;
  user.verif      = PETSC_FALSE;
  user.check_admissible = PETSC_FALSE;
  user.cmb        = NULL;

  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"ice_","options to ice","");CHKERRQ(ierr);
  ierr = PetscOptionsReal(
      "-A", "set value of ice softness A in units Pa-3 s-1",
      "ice.c",user.A_ice,&user.A_ice,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsBool(
      "-check_admissible", "check admissibility of iterate at start of residual evaluation FormFunctionLocal()",
      "ice.c",user.check_admissible,&user.check_admissible,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsReal(
      "-D0", "representative value of diffusivity (used in regularizing D) in units m2 s-1",
      "ice.c",user.D0,&user.D0,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsReal(
      "-delta", "dimensionless regularization for slope in SIA formulas",
      "ice.c",user.delta,&user.delta,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsBool(
      "-dump", "save final state (H, b)",
      "ice.c",dump,&dump,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsReal(
      "-eps", "dimensionless regularization for diffusivity D",
      "ice.c",user.eps,&user.eps,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsBool(
      "-exact_init", "initialize with dome exact solution",
      "ice.c",exact_init,&exact_init,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsReal(
      "-L", "side length of domain in meters",
      "ice.c",user.L,&user.L,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsReal(
      "-lambda", "amount of upwinding; lambda=0 is none and lambda=1 is full",
      "ice.c",user.lambda,&user.lambda,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsReal(
      "-n", "value of Glen exponent n",
      "ice.c",user.n_ice,&user.n_ice,NULL);CHKERRQ(ierr);
  if (user.n_ice <= 1.0) {
      SETERRQ1(PETSC_COMM_WORLD,1,
          "ERROR: n = %f not allowed ... n > 1.0 is required\n",user.n_ice);
  }
  ierr = PetscOptionsReal(
      "-rho", "ice density in units kg m3",
      "ice.c",user.rho_ice,&user.rho_ice,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsBool(
      "-verif", "use dome exact solution for verification",
      "ice.c",user.verif,&user.verif,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsEnd();CHKERRQ(ierr);

  // derived constant computed after other ice properties are set
  user.Gamma = 2.0 * PetscPowReal(user.rho_ice*user.g,user.n_ice)
                   * user.A_ice / (user.n_ice+2.0);

  ierr = SetFromOptions_CMBModel(&cmb,user.secpera);
  user.cmb = &cmb;

  // DMDA for the cell-centered grid
  ierr = DMDACreate2d(PETSC_COMM_WORLD,
                      DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,
                      DMDA_STENCIL_BOX,
                      5,5, PETSC_DECIDE,PETSC_DECIDE,
                      1, 1,        // dof=1, stencilwidth=1
                      NULL,NULL,&da);
  ierr = DMSetFromOptions(da); CHKERRQ(ierr);
  ierr = DMSetUp(da); CHKERRQ(ierr);  // this must be called BEFORE SetUniformCoordinates
  ierr = DMDASetUniformCoordinates(da,0.0,user.L,0.0,user.L,-1.0,-1.0);CHKERRQ(ierr);
  ierr = DMSetApplicationContext(da, &user);CHKERRQ(ierr);

  // create and configure the SNES to solve a NCP/VI at each step
  ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr);
  ierr = SNESSetDM(snes,da);CHKERRQ(ierr);
  ierr = SNESSetApplicationContext(snes,&user);CHKERRQ(ierr);
  ierr = DMDASNESSetFunctionLocal(da,INSERT_VALUES,
               (DMDASNESFunction)FormFunctionLocal,&user); CHKERRQ(ierr);
  ierr = SNESSetType(snes,SNESVINEWTONSSLS); CHKERRQ(ierr);
  ierr = SNESVISetComputeVariableBounds(snes,&FormBounds); CHKERRQ(ierr);
  ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);

  // set up initial iterate
  ierr = DMCreateGlobalVector(da,&H);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject)H,"H"); CHKERRQ(ierr);
  if (exact_init) {
      ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);
      ierr = DMDAVecGetArray(da,H,&aH); CHKERRQ(ierr);
      ierr = DomeThicknessLocal(&info,aH,&user); CHKERRQ(ierr);
      ierr = DMDAVecRestoreArray(da,H,&aH); CHKERRQ(ierr);
  } else {
      ierr = VecSet(H,0.0); CHKERRQ(ierr);
  }

  // solve
  ierr = SNESSolve(snes,NULL,H); CHKERRQ(ierr);
  ierr = SNESGetConvergedReason(snes,&reason); CHKERRQ(ierr);
  if (reason <= 0) {
      ierr = PetscPrintf(PETSC_COMM_WORLD,
          "WARNING: SNES not converged ... use -snes_converged_reason to check\n"); CHKERRQ(ierr);
  }

  // get solution & DM on fine grid (which may have changed) after solve
  ierr = VecDestroy(&H); CHKERRQ(ierr);
  ierr = DMDestroy(&da); CHKERRQ(ierr);
  ierr = SNESGetDM(snes,&da); CHKERRQ(ierr); /* do not destroy da */
  ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);
  ierr = SNESGetSolution(snes,&H); CHKERRQ(ierr); /* do not destroy H */
  ierr = PetscObjectSetName((PetscObject)H,"H"); CHKERRQ(ierr);

  // compute performance measures; note utility of reporting last grid,
  //   last snesit/kspit when doing -snes_grid_sequence
  ierr = SNESGetIterationNumber(snes,&snesit); CHKERRQ(ierr);  //
  ierr = SNESGetKSP(snes,&ksp); CHKERRQ(ierr);
  ierr = KSPGetIterationNumber(ksp,&kspit); CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,
      "done on %d x %d grid ... SNES iters = %d, last KSP iters = %d\n",
      info.mx,info.my,snesit,kspit); CHKERRQ(ierr);

  // dump state (H,b) if requested
  if (dump) {
      char           filename[1024];
      PetscViewer    viewer;
      Vec            b;
      double         **ab;
      ierr = VecDuplicate(H,&b); CHKERRQ(ierr);
      ierr = PetscObjectSetName((PetscObject)b,"b"); CHKERRQ(ierr);
      if (user.verif) {
          ierr = VecSet(b,0.0); CHKERRQ(ierr);
      } else {
          ierr = DMDAVecGetArray(da,b,&ab); CHKERRQ(ierr);
          ierr = FormBedLocal(&info,0,ab,&user); CHKERRQ(ierr);
          ierr = DMDAVecRestoreArray(da,b,&ab); CHKERRQ(ierr);
      }
      ierr = sprintf(filename,"ice_%dx%d.dat",info.mx,info.my);
      ierr = PetscPrintf(PETSC_COMM_WORLD,"writing PETSC binary file %s ...\n",filename); CHKERRQ(ierr);
      ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_WRITE,&viewer); CHKERRQ(ierr);
      ierr = VecView(b,viewer); CHKERRQ(ierr);
      ierr = VecView(H,viewer); CHKERRQ(ierr);
      ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
      VecDestroy(&b);
  }

  // compute error in verification case
  if (user.verif) {
      Vec Hexact;
      double infnorm, onenorm;
      ierr = VecDuplicate(H,&Hexact); CHKERRQ(ierr);
      ierr = DMDAVecGetArray(da,Hexact,&aH); CHKERRQ(ierr);
      ierr = DomeThicknessLocal(&info,aH,&user); CHKERRQ(ierr);
      ierr = DMDAVecRestoreArray(da,Hexact,&aH); CHKERRQ(ierr);
      ierr = VecAXPY(H,-1.0,Hexact); CHKERRQ(ierr);    // H <- H + (-1.0) Hexact
      VecDestroy(&Hexact);
      ierr = VecNorm(H,NORM_INFINITY,&infnorm); CHKERRQ(ierr);
      ierr = VecNorm(H,NORM_1,&onenorm); CHKERRQ(ierr);
      ierr = PetscPrintf(PETSC_COMM_WORLD,
          "numerical errors: |H-Hexact|_inf = %.3f, |H-Hexact|_average = %.3f\n",
          infnorm,onenorm/(double)(info.mx*info.my)); CHKERRQ(ierr);
  }

  SNESDestroy(&snes);
  return PetscFinalize();
}

PetscErrorCode FormBedLocal(DMDALocalInfo *info, int stencilwidth, double **ab, AppCtx *user) {
  int          j,k,r,s;
  const double dx = user->L / (double)(info->mx-1),
               dy = user->L / (double)(info->my-1),
               Z = PETSC_PI / user->L;
  double       x, y, b;
  // vaguely-random frequencies and coeffs generated by fiddling; see randbed.py
  const int    nc = 4,
               jc[4] = {1, 3, 6, 8},
               kc[4] = {1, 3, 4, 7};
  const double scalec = 750.0,
               C[4][4] = { { 2.00000000,  0.33000000, -0.55020034,  0.54495520},
                           { 0.50000000,  0.45014486,  0.60551833, -0.52250644},
                           { 0.93812068,  0.32638429, -0.24654812,  0.33887052},
                           { 0.17592361, -0.35496741,  0.22694547, -0.05280704} };
  PetscFunctionBeginUser;
  // go through owned portion of grid and compute  b(x,y)
  for (k = info->ys-stencilwidth; k < info->ys + info->ym+stencilwidth; k++) {
      y = k * dy;
      for (j = info->xs-stencilwidth; j < info->xs + info->xm+stencilwidth; j++) {
          if (j < 0 || j >= info->mx-1 || k < 0 || k >= info->my-1)
              continue;
          x = j * dx;
          // b(x,y) is sum of a few sines
          b = 0.0;
          for (r = 0; r < nc; r++) {
              for (s = 0; s < nc; s++) {
                  b += C[r][s] * sin(jc[r] * Z * x) * sin(kc[s] * Z * y);
              }
          }
          ab[k][j] = scalec * b;
      }
  }
  PetscFunctionReturn(0);
}

//  for call-back: tell SNESVI (variational inequality) that we want
//    0.0 <= H < +infinity
PetscErrorCode FormBounds(SNES snes, Vec Xl, Vec Xu) {
  PetscErrorCode ierr;
  PetscFunctionBeginUser;
  ierr = VecSet(Xl,0.0); CHKERRQ(ierr);
  ierr = VecSet(Xu,PETSC_INFINITY); CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

// value of gradient at a point
typedef struct {
    double x,y;
} Grad;

/* We factor the SIA flux as
    q = - H^{n+2} sigma(|grad s|) grad s
where sigma is the slope-dependent part
    sigma(z) = Gamma z^{n-1}.
Also
    D = H^{n+2} sigma(|grad s|)
so that q = - D grad s.  */
static double sigma(Grad gH, Grad gb, const AppCtx *user) {
    const double sx = gH.x + gb.x,
                 sy = gH.y + gb.y,
                 slopesqr = sx * sx + sy * sy + user->delta * user->delta;
    return user->Gamma * PetscPowReal(slopesqr,(user->n_ice-1.0)/2);
}

/* Pseudo-velocity from bed slope:  W = - sigma * grad b. */
static Grad W(double sigma, Grad gb) {
    Grad W;
    W.x = - sigma * gb.x;
    W.y = - sigma * gb.y;
    return W;
}

/* DCS = diffusivity from the continuation scheme:
     D(eps) = (1-eps) sigma H^{n+2} + eps D_0
so D(1)=D_0 and D(0)=sigma H^{n+2}. */
static double DCS(double sigma, double H, const AppCtx *user) {
  return (1.0 - user->eps) * sigma * PetscPowReal(PetscAbsReal(H),user->n_ice+2.0)
         + user->eps * user->D0;
}

/* Flux component from the non-sliding SIA on a general bed. */
PetscErrorCode SIAflux(Grad gH, Grad gb, double H, double Hup, PetscBool xdir,
                       double *D, double *q, const AppCtx *user) {
  const double mysig = sigma(gH,gb,user),
               myD   = DCS(mysig,H,user);
  const Grad   myW   = W(mysig,gb);
  PetscFunctionBeginUser;
  if (D) {
      *D = myD;
  }
  if (xdir && q) {
      *q = - myD * gH.x + myW.x * PetscPowReal(PetscAbsReal(Hup),user->n_ice+2.0);
  } else {
      *q = - myD * gH.y + myW.y * PetscPowReal(PetscAbsReal(Hup),user->n_ice+2.0);
  }
  PetscFunctionReturn(0);
}

// gradients of weights for Q^1 interpolant
static const double gx[4] = {-1.0,  1.0, 1.0, -1.0},
                    gy[4] = {-1.0, -1.0, 1.0,  1.0};

static double fieldatpt(double xi, double eta, double f[4]) {
  // weights for Q^1 interpolant
  double x[4] = { 1.0-xi,      xi,  xi, 1.0-xi},
         y[4] = {1.0-eta, 1.0-eta, eta,    eta};
  return   x[0] * y[0] * f[0] + x[1] * y[1] * f[1]
         + x[2] * y[2] * f[2] + x[3] * y[3] * f[3];
}

static double fieldatptArray(int u, int v, double xi, double eta, double **f) {
  double ff[4] = {f[v][u], f[v][u+1], f[v+1][u+1], f[v+1][u]};
  return fieldatpt(xi,eta,ff);
}

static Grad gradfatpt(double xi, double eta, double dx, double dy, double f[4]) {
  Grad gradf;
  double x[4] = { 1.0-xi,      xi,  xi, 1.0-xi},
         y[4] = {1.0-eta, 1.0-eta, eta,    eta};
  gradf.x =   gx[0] * y[0] * f[0] + gx[1] * y[1] * f[1]
            + gx[2] * y[2] * f[2] + gx[3] * y[3] * f[3];
  gradf.y =    x[0] *gy[0] * f[0] +  x[1] *gy[1] * f[1]
            +  x[2] *gy[2] * f[2] +  x[3] *gy[3] * f[3];
  gradf.x /= dx;
  gradf.y /= dy;
  return gradf;
}

static Grad gradfatptArray(int u, int v, double xi, double eta, double dx, double dy,
                           double **f) {
  double ff[4] = {f[v][u], f[v][u+1], f[v+1][u+1], f[v+1][u]};
  return gradfatpt(xi,eta,dx,dy,ff);
}

// indexing of the 8 quadrature points along the boundary of the control volume in M*
// point s=0,...,7 is in element (j,k) = (j+je[s],k+ke[s])
static const int  je[8] = {0,  0, -1, -1, -1, -1,  0,  0},
                  ke[8] = {0,  0,  0,  0, -1, -1, -1, -1},
                  ce[8] = {0,  3,  1,  0,  2,  1,  3,  2};

// direction of flux at 4 points in each element
static const PetscBool xdire[4] = {PETSC_TRUE, PETSC_FALSE, PETSC_TRUE, PETSC_FALSE};

// local (element-wise) coords of quadrature points for M*
static const double locx[4] = {  0.5, 0.75,  0.5, 0.25},
                    locy[4] = { 0.25,  0.5, 0.75,  0.5};

/* FormFunctionLocal  =  call-back by SNES using DMDA info.

Evaluates residual FF on local process patch:
   FF_{j,k} = \int_{\partial V_{j,k}} \mathbf{q} \cdot \mathbf{n}
              - m_{j,k} \Delta x \Delta y
where V_{j,k} is the control volume centered at (x_j,y_k).

Regarding indexing locations along the boundary of the control volume where
flux is evaluated, this figure shows the control volume centered at (x_j,y_k)
and the four elements it meets.  Quadrature uses 8 points on the boundary of
the control volume, numbered s=0,...,7:

     -------------------
    |         |         |
    |    ..2..|..1..    |
    |   3:    |    :0   |
  k |--------- ---------|
    |   4:    |    :7   |
    |    ..5..|..6..    |
    |         |         |
     -------------------
              j

Regarding flux-component indexing on the element indexed by (j,k) node,
the value  (aqquad[c])[k][j] for c=0,1,2,3 is an x-component at "*" and
a y-component at "%"; note (x_j,y_k) is lower-left corner:

     -------------------
    |         :         |
    |         *2        |
    |    3    :    1    |
    |....%.... ....%....|
    |         :         |
    |         *0        |
    |         :         |
    @-------------------
  (j,k)

*/
PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, double **aHin,
                                        double **FF, AppCtx *user) {
  PetscErrorCode  ierr;
  const double    dx = user->L / (double)(info->mx-1),
                  dy = user->L / (double)(info->my-1);
  // coefficients of quadrature evaluations along the boundary of the control volume in M*
  const double    coeff[8] = {dy/2, dx/2, dx/2, -dy/2, -dy/2, -dx/2, -dx/2, dy/2};
  const PetscBool upwind = (user->lambda > 0.0);
  const double    upmin = (1.0 - user->lambda) * 0.5,
                  upmax = (1.0 + user->lambda) * 0.5;
  int             c, j, k, s;
  double          H, Hup, lxup, lyup, **aqquad[4], **ab, **aH, DSIA_ckj, qSIA_ckj,
                  M, x, y;
  Grad            gH, gb;
  Vec             qquad[4], Hcopy, b;

  PetscFunctionBeginUser;

  // copy and set boundary conditions to zero
  ierr = DMGetLocalVector(info->da, &Hcopy); CHKERRQ(ierr);
  ierr = DMDAVecGetArray(info->da,Hcopy,&aH); CHKERRQ(ierr);
  for (k = info->ys-1; k <= info->ys + info->ym; k++) {
      for (j = info->xs-1; j <= info->xs + info->xm; j++) {
          if (j < 0 || j > info->mx-1 || k < 0 || k > info->my-1)
              continue;
          if (user->check_admissible && aHin[k][j] < 0.0) {
              SETERRQ3(PETSC_COMM_WORLD,1,
                       "ERROR: non-admissible value H[k][j] = %.3e < 0.0 at j,k = %d,%d\n",
                       aHin[k][j],j,k);
          }
          if (j == 0 || j == info->mx-1 || k == 0 || k == info->my-1) {
              if (j >= info->xs && j < info->xs+info->xm && k >= info->ys && k < info->ys+info->ym)
                  FF[k][j] = aHin[k][j];   // FIXME scaling?
              aH[k][j] = 0.0;
          } else
              aH[k][j] = aHin[k][j];
      }
  }

  // get bed elevation b(x,y) on this grid
  ierr = DMGetLocalVector(info->da, &b); CHKERRQ(ierr);
  ierr = DMDAVecGetArray(info->da,b,&ab); CHKERRQ(ierr);
  if (user->verif) {
      ierr = VecSet(b,0.0); CHKERRQ(ierr);
  } else {
      ierr = FormBedLocal(info,1,ab,user); CHKERRQ(ierr);  // get stencil width
  }

  // working space for fluxes; see text for face location of flux evaluation
  for (c = 0; c < 4; c++) {
      ierr = DMGetLocalVector(info->da, &(qquad[c])); CHKERRQ(ierr);
      ierr = DMDAVecGetArray(info->da,qquad[c],&(aqquad[c])); CHKERRQ(ierr);
  }

  // loop over locally-owned elements, including ghosts, to get fluxes q at
  // c = 0,1,2,3 points in element;  note start at (xs-1,ys-1)
  for (k = info->ys-1; k < info->ys + info->ym; k++) {
      for (j = info->xs-1; j < info->xs + info->xm; j++) {
          if (j < 0 || j >= info->mx-1 || k < 0 || k >= info->my-1)
              continue;
          for (c=0; c<4; c++) {
              H  = fieldatptArray(j,k,locx[c],locy[c],aH);
              gH = gradfatptArray(j,k,locx[c],locy[c],dx,dy,aH);
              gb = gradfatptArray(j,k,locx[c],locy[c],dx,dy,ab);
              if (upwind) {
                  if (xdire[c] == PETSC_TRUE) {
                      lxup = (gb.x <= 0.0) ? upmin : upmax;
                      lyup = locy[c];
                  } else {
                      lxup = locx[c];
                      lyup = (gb.y <= 0.0) ? upmin : upmax;
                  }
                  Hup = fieldatptArray(j,k,lxup,lyup,aH);
              } else
                  Hup = H;
              ierr = SIAflux(gH,gb,H,Hup,xdire[c],
                             &DSIA_ckj,&qSIA_ckj,user); CHKERRQ(ierr);
              aqquad[c][k][j] = qSIA_ckj;
          }
      }
  }

  // loop over nodes, not including ghosts, to get function F(H) from quadature over
  // s = 0,1,...,7 points on boundary of control volume (rectangle) around node
  for (k=info->ys; k<info->ys+info->ym; k++) {
      for (j=info->xs; j<info->xs+info->xm; j++) {
          if (j == 0 || j == info->mx-1 || k == 0 || k == info->my-1)
              continue;
          // climatic mass balance
          if (user->verif) {
              x = j * dx;
              y = k * dy;
              M = DomeCMB(x,y,user);
          } else {
              M = M_CMBModel(user->cmb,ab[k][j] + aH[k][j]);  // s=b+H is surface elevation
          }
          FF[k][j] = - M * dx * dy;
          // now add integral over control volume boundary using two
          // quadrature points on each side
          for (s=0; s<8; s++)
              FF[k][j] += coeff[s] * aqquad[ce[s]][k+ke[s]][j+je[s]];
      }
  }

  // restore working space and bed
  for (c = 0; c < 4; c++) {
      ierr = DMDAVecRestoreArray(info->da,qquad[c],&(aqquad[c])); CHKERRQ(ierr);
      ierr = DMRestoreLocalVector(info->da, &(qquad[c])); CHKERRQ(ierr);
  }
  ierr = DMDAVecRestoreArray(info->da,Hcopy,&aH); CHKERRQ(ierr);
  ierr = DMRestoreLocalVector(info->da, &Hcopy); CHKERRQ(ierr);
  ierr = DMDAVecRestoreArray(info->da,b,&ab); CHKERRQ(ierr);
  ierr = DMRestoreLocalVector(info->da, &b); CHKERRQ(ierr);
  PetscFunctionReturn(0);
}


text:
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_rng.h>
#include <math.h>
#include <time.h>
#include <stdbool.h>

#define RND gsl_rng_uniform(r)

double gsl_ran_gaussian(const gsl_rng *r, double sigma);
double gsl_ran_gaussian_pdf(double x, double sigma);
double gsl_ran_beta_pdf(double x, double a, double b);
double gsl_ran_flat(const gsl_rng *r, double a, double b);
double gsl_ran_flat_pdf(double x, double a, double b);
unsigned int gsl_ran_binomial(const gsl_rng * r, double p, unsigned int n);

double mean_int(int array[], int array_size) {
	double avg = 0;
	int h;
	for (h=0; h<array_size; h++) {
		avg += array[h];
	}
	avg /= array_size;
	return avg;
}

// double ABC_distance_BLR(double y[], double y_s[], int size_data) {
// 	double dist = 0;
// 	int k;
// 	for (k=0; k<size_data; k++) {
// 		dist += (y[k] - y_s[k]) * (y[k] - y_s[k]);
// 	}
// 	return dist;
// }

// double prior(int l, double theta_star, double a, double b) {
// 	double res;
// 	if (l==0) {
// 		res = gsl_ran_flat_pdf(theta_star, a, b);
// 		return res;
// 	}
// 	if (l==1) {
// 		res = gsl_ran_flat_pdf(theta_star, a, b);
// 		return res;
// 	}
// 	else {return -1;}
// }

double prior_BB(double theta_star, double a, double b) {

	double res;
	res = gsl_ran_beta_pdf(theta_star, a, b);
	return res;

}

// double sort(double number[], int n) {

// 	int i, j;
// 	double a;
//     for (i = 0; i < n; i++) {
//         for (j = i + 1; j < n; j++) {
//             if (number[i] > number[j]) {
//                 a =  number[i];
//                 number[i] = number[j];
//                 number[j] = a;
//             }
//         }
//     }
//     return 0;
// }

int select_particle(double *weights, double *theta, int time, int n_particles, const gsl_rng *r) {

	FILE * w_theta;
	char weights_theta[500];
	double r2;
	double aux = 0.;
	int i = 0;
	int indicator = -1;

	double* theta_vec;
	theta_vec = malloc((n_particles) * sizeof(double*));

	sprintf(weights_theta, "w_theta_time_%d_part_%d.txt", time, n_particles);
	w_theta = fopen(weights_theta, "r");
	if (w_theta == NULL) {
		printf("File %s not found!!!\n", weights_theta);
		return 1;
	}
	fscanf(w_theta, "%*[^\n]");

	r2 = RND;
	for (i=0; i<n_particles; i++) {
		fscanf(w_theta, "%lf %lf", &weights[i], &theta_vec[i]);
		aux += weights[i];
			if (r2 > aux && indicator < 0) {
				*theta = theta_vec[i];
				indicator = 1;
		}
	}

	// for (i=0; i<n_particles; i++) {
	// 	fscanf(w_theta, "%lf %lf", &weights[i], &theta_vec[i]);
	// 	aux += weights[i];
	// 	if (r2<aux) {
	// 		if (i_th<0) {
	// 			i_th = i;
	// 		}
	// 	}
	// }

	fclose(w_theta);
	free(theta_vec);

	return 0;

}

void calculate_weight_BB(double *weight, double *weights, double theta, double theta_star, int n_particles, double a0, double b0, double std_ker) {
	// Calculates the new weights for the successive ABC round
	int i;
	double tot_prior_pdf = 1.;
	double kernel_value_pdf[n_particles];
	for (i=0; i<n_particles; i++) {kernel_value_pdf[i] = 1.;}

	tot_prior_pdf *= prior_BB(theta_star, a0, b0);

	double denominator = 0.;
	for (i=0; i<n_particles; i++) {
		kernel_value_pdf[i] *= gsl_ran_gaussian_pdf(theta_star - theta, std_ker);
		denominator += weights[i] * kernel_value_pdf[i];
	}
	*weight = tot_prior_pdf/(float)(denominator);

}

int ker_perturb(double *theta_star, const gsl_rng *r, double theta, double a0, double b0, double std_ker) {
	//Perturb the i_th parameter theta via a gaussian kernel with std = std_ker
	*theta_star = gsl_ran_gaussian(r, std_ker) + theta;
	//If the prior on theta_star[l] is zero when calculated in theta_star[l],
	//the new particle is rejected and we start over
	if (prior_BB(*theta_star, a0, b0) == 0.) {
		return 1;
	}

	return 0;

}

int read_threshold(double* epsilon_t, int time) {

	FILE * epsilon;
	char epsilon_file[500];
	sprintf(epsilon_file, "epsilon_%d.txt", time);
	epsilon = fopen(epsilon_file, "r");
	if (epsilon == NULL) {
		printf("File %s not found!!!\n", epsilon_file);
		return 1;
	}
	fscanf(epsilon, "%lf", epsilon_t);
	fclose(epsilon);

	return 0;

}

// Distance function equal to zero for sample to be accepted
int verify_par_BB(double theta_star, double mean_y, const gsl_rng *r, int size_data, int n, double epsilon_t, int time, int seed) {

	int k;
	int g[size_data];
	double mean_g;
	for (k=0; k<size_data; k++) {g[k] = gsl_ran_binomial(r, theta_star, n);}
	mean_g = mean_int(g, size_data);
	printf("mean g: %f\nmean_y: %f\nepsilon_t: %f\n\n", mean_g, mean_y, epsilon_t);
	if ( fabs(mean_g - mean_y) <= epsilon_t ) {
		FILE* epsilon_file;
		char epsilon_char[500];
		sprintf(epsilon_char, "epsilon_%d_%d.txt", time+1, seed);
		epsilon_file = fopen(epsilon_char, "w");
		fprintf(epsilon_file, "%f", fabs(mean_g - mean_y));
		fclose(epsilon_file);
		return 0;
	}

	return 1;

}

void write_particles_BB(double theta_star, int time, int seed, int striding_size, double weight) {
	// Writes ABC particles on file
	FILE *output_file;
	char output_char[500];

	sprintf(output_char, "w_theta_time_%d_%d.txt", time+1, 1+(seed-1)/striding_size);
	output_file = fopen(output_char, "a");
	fprintf(output_file, "%.10f %.5f\n", weight, theta_star);
	fclose(output_file);

}

int import_data_BB(int *y, int size_data) {
	// Import data for beta-binomial model
	int d;

	FILE *file_data;
	char file_name[500];

	sprintf(file_name, "data_bb.txt");
	file_data = fopen("data_bb.txt", "r");
	if (file_data == NULL) {
		return 1;
	}
	for (d=0; d<size_data; d++) {
		fscanf(file_data, "%d", &y[d]);
	}
	fclose(file_data);

	return 0;

}

void write_epsilon(double distance, int time, int seed, int striding_size) {
	// Write epsilon for successive ABC round to file
	FILE *epsilon_;
	char eps_char[500];
	sprintf(eps_char, "epsilon_%d_%d.txt", time+1, 1+(seed-1)/striding_size);
	epsilon_ = fopen(eps_char, "a");
	fprintf(epsilon_, "%f\n", distance);
	fclose(epsilon_);

}

int read_epsilon_3d(double *eps_t, int time) {
	// Read 3d epsilon
	FILE *epsilon_;
	char eps_char[500];

	sprintf(eps_char, "epsilon_%d.txt", time);
	epsilon_ = fopen(eps_char, "r");
	if (epsilon_ == NULL) {
		return 1;
	}

	fscanf(epsilon_, "%lf %lf %lf", &eps_t[0], &eps_t[1], &eps_t[2]);
	fclose(epsilon_);

	return 0;

}

void write_epsilon_3d(double *distance, int time, int seed, int striding_size) {

	FILE *epsilon_;
	char eps_char[500];
	sprintf(eps_char, "epsilon_%d_%d.txt", time+1, 1+(seed-1)/striding_size);
	epsilon_ = fopen(eps_char, "a");
	fprintf(epsilon_, "%f %f %f\n", distance[0], distance[1], distance[2]);
	fclose(epsilon_);

}

ext: h


text:
// Copyright (c) Microsoft Corporation.
// Licensed under the MIT License.

#pragma once

#include "BaseFeatureMatcher.h"
#include "OnlineBow.h"

#include <gsl\gsl>
#include <unordered_map>

namespace mage
{
    class Keyframe;

    class OnlineBowFeatureMatcher : public BaseFeatureMatcher
    {
    public:
        OnlineBowFeatureMatcher(const OnlineBow& onlineBow, const Id<Keyframe>& id, gsl::span<const ORBDescriptor> features);

        virtual size_t QueryFeatures(const ORBDescriptor& descriptor, std::vector<ptrdiff_t>& matches) const override;

        virtual const Id<Keyframe>& GetId() const override;

    private:
        // a Map of NodeId <==> the indexes of the features which are assigned to the Node.
        std::unordered_map<ptrdiff_t, std::vector<ptrdiff_t>> m_featureMap;
        const OnlineBow& m_onlineBow;
        const Id<Keyframe> m_id;
    };
}

max_stars_repo_name: syntheticmagus/mageslam

max_issues_repo_name: syntheticmagus/mageslam

max_forks_repo_name: syntheticmagus/mageslam


text:
#ifndef __SPATIAL_STATISTICS_H__
#define __SPATIAL_STATISTICS_H__

#include <cmath>
#include <gsl/gsl_sf_bessel.h>

extern "C" double gsl_sf_gamma(const double);
extern "C" double gsl_sf_hyperg_1F1(double, double, double);

template <typename T> struct Spatial_Statistics
{
  private:
    T phi, nu;
    T con;
    int d;

  public:
    Spatial_Statistics(T phi, T nu, int dim)
    {
        this->phi = phi;
        this->nu = nu;
        // this->alpha = alpha;
        this->d = dim;

        con = pow(2, (nu - 1)) * tgamma(nu);
        con = 1.0 / con;

        // c1 = pow(h / sigma, d) / pow(sigma, alpha);
        // c2 = -pow(2.0, alpha) * gsl_sf_gamma((alpha + (T)d) / 2) / pow(M_PI, (T)d / 2) / gsl_sf_gamma((T)d / 2);
    }

    // T G(T x) const
    //{
    //    return c2 * gsl_sf_hyperg_1F1((alpha + (T)d) / 2, (T)d / 2, -x * x);
    //}

    T operator()(T *pt_x, T *pt_y) const
    {
        T diff = 0.0;
        for (int i = 0; i < d; i++)
            diff += (pt_x[i] - pt_y[i]) * (pt_x[i] - pt_y[i]);
        T dist = sqrt(diff);

        dist = 4.0 * sqrt(2.0 * nu) * (dist / phi);
        if (dist == 0.0)
        {
            return 1.0;
        }
        else
        {
            return con * pow(dist, nu) * gsl_sf_bessel_Knu(nu, dist);
        }
    }
};

template <typename T> struct MaternKernel
{
  private:
    T h, sigma, alpha;
    T c1, c2;
    int d;

  public:
    MaternKernel(T h, T sigma, T alpha, int dim)
    {
        this->h = h;
        this->sigma = sigma;
        this->alpha = alpha;
        this->d = dim;

        c1 = pow(h / sigma, d) / pow(sigma, alpha);
        c2 = -pow(2.0, alpha) * gsl_sf_gamma((alpha + (T)d) / 2.0) / pow(M_PI, (T)d / 2.0) / gsl_sf_gamma((T)d / 2.0);
    }

    T G(T x) const
    {
        return c2 * gsl_sf_hyperg_1F1((alpha + (T)d) / 2.0, (T)d / 2.0, -x * x);
    }

    T operator()(T *pt_x, T *pt_y) const
    {
        T diff = 0.0;
        for (int i = 0; i < d; i++)
            diff += (pt_x[i] - pt_y[i]) * (pt_x[i] - pt_y[i]);
        T dist = sqrt(diff);

        return c1 * G(dist / sigma);
    }
};

#endif


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h> // high level interface to blas
#include <gsl/gsl_linalg.h> // lapack

void print_matrix(gsl_matrix* m, int nx, int ny);

int main (void) {
	int i, j;
	int signum = 0;
	double determinant_m = 0.0;
	gsl_permutation  *m_LU_permutation = gsl_permutation_alloc(8);
	gsl_matrix *m_test = gsl_matrix_alloc(8, 8);
	gsl_matrix *m_original = gsl_matrix_alloc(8,8);
	gsl_matrix *m_inverted = gsl_matrix_alloc(8,8);

	/*for (i = 0; i < 8; i++){
		for(j = 0; j < 8; j++){
			gsl_matrix_set(m_test, i, j, 1);
		}
	}*/

	gsl_matrix_set(m_test, 0, 7, 1);
	gsl_matrix_set(m_test, 7,0, 1);
	gsl_matrix_set(m_test, 3,2, 5);
	gsl_matrix_set(m_test, 5,2, 5);

	for( i = 0; i < 8; i++){
		gsl_matrix_set(m_test, i, i, 3);
	}

	// set the oritinal one too
	gsl_matrix_memcpy(m_original, m_test);
	// now print the matrix out
	//print_matrix(m_test, 8, 8);

	printf("\n");
	print_matrix(m_original, 8,8);

	// simple matrix operations
	//gsl_matrix_add(m_test, m_test);
	//gsl_matrix_add_constant(m_test, 1.05);

	//print_matrix(m_test, 8,8);

	// creates the lu decomp of the matrix
	// need to allocate the permutation
	gsl_linalg_LU_decomp(m_test, m_LU_permutation, &signum);

	// calculate the determinant of m_test from the lu
	determinant_m = gsl_linalg_LU_det(m_test, signum);
	printf("det_m = %g\n", determinant_m);
	// invert the matrix
	gsl_linalg_LU_invert(m_test, m_LU_permutation, m_inverted);

	// multiply the inverse and the matrix together  (can't use the LU one bc it's fucked up)
	// m_test should be I
	gsl_blas_dgemm( CblasNoTrans,CblasNoTrans , 1.0, m_original, m_inverted, 0.0, m_test);

	print_matrix(m_test, 8,8);

	gsl_permutation_free(m_LU_permutation);
	gsl_matrix_free(m_inverted);
	gsl_matrix_free(m_test);
	return(0);
}

void print_matrix(gsl_matrix* m, int nx, int ny){
	int i,j;
	for(i = 0; i < nx; i++){
		for(j= 0; j < ny; j++){
			printf("%g ", gsl_matrix_get(m, i, j));
		}
		printf("\n");
 	}
}

ext: c

max_stars_repo_name: jackdawjackdaw/emulator

max_issues_repo_name: MADAI/MADAIEmulator

max_forks_repo_name: MADAI/MADAIEmulator


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>

// convex hull is obtained as cvh(n_bin[0],...,n_bin[end])
// inputs are y = Ef(x). x:n_data by 4, y:n_data by 1
void obtain_convex_h_Nd(int n_data, double *ef2cvh, int *x, double *y, int Ndim, int *nCat_min, int *nCat_max, double *E0) {

    int i, j, k, l, nid, id1, id2, ctr;
    double tmp;
    int convex_h_ids[n_data];   //preallocated, not all filled
    double convex_h_x[n_data][Ndim];   //preallocated, not all filled
    double convex_h_y[n_data];  //preallocated, not all filled

    printf("Convex hull is created.\n");
    printf("Input data\n");
    printf("---------------------------\n");
    printf("  nL   nN   nM   nC   nV    Ef\n");
    for (i=0;i<n_data;i++) {
        for (j=0;j<Ndim;j++)
            printf("%4d ",x[Ndim*i+j]);
        printf("  % .4f\n",y[i]);
    }
    printf("---------------------------\n");

    // calculate the gradient and find Ndim base convex hull points
    double dEdx[n_data][Ndim];
    double dx[n_data][Ndim];
    for (i=0;i<n_data;i++) {
        printf("%d-data:%d,%d,%d,%d,%d, dE/dx=[",i,x[i*Ndim],x[i*Ndim+1],x[i*Ndim+2],x[i*Ndim+3],x[i*Ndim+4]);
        for (j=0;j<Ndim;j++) {
            dx[i][j] = (double) (nCat_max[j]-x[i*Ndim+j])*1.0/(nCat_max[j]-nCat_min[j]);
            if (dx[i][j] == 0.0)  // case of convex hull
                dEdx[i][j] = 0.0;
            else
                dEdx[i][j] = (E0[j]-y[i]) / dx[i][j];
            printf("% .2e ",dEdx[i][j]);
        }
        printf("]\n");
    }
    // check if the point is the convex hull
    int ctr_set[Ndim];
    nid = 0;
    for (i=0;i<n_data;i++) {
        if (dx[i][j] == 0.0) {
            convex_h_ids[nid] = i;
            for (j=0;j<Ndim;j++)
                convex_h_x[nid][j] = (double) x[i*Ndim+j]*1.0;
            convex_h_y[nid] = y[i];
            nid++;
            continue;
        }
        for (k=0;k<Ndim;k++) {
            ctr_set[k] = 1;
            for (j=0;j<n_data;j++) {
                if (j==i)
                    continue;
                if (dEdx[j][k] < dEdx[i][k]) {
                    ctr_set[k] = 0;
                    break;
                }
            }
        }
        ctr = 0;
        for (k=0;k<Ndim;k++)
            ctr += ctr_set[k];
        if (ctr > 0) {
            convex_h_ids[nid] = i;
            for (j=0;j<Ndim;j++)
                convex_h_x[nid][j] = (double) x[i*Ndim+j]*1.0;
            convex_h_y[nid] = y[i];
            nid++;
        }
    }
    printf("%d cvh points were identified.\n",nid);
    for (i=0;i<nid;i++) {
        printf("%d ",convex_h_ids[i]);
        for (j=0;j<Ndim;j++)
            printf("%f ",convex_h_x[i][j]);
        printf("%f\n",convex_h_y[i]);
    }

    // A*c = b
    gsl_matrix * A = gsl_matrix_alloc (Ndim,Ndim);
    gsl_matrix * V = gsl_matrix_alloc (Ndim,Ndim);
    gsl_vector * S = gsl_vector_alloc (Ndim);
    gsl_vector * work = gsl_vector_alloc (Ndim);
    gsl_vector * b = gsl_vector_alloc (Ndim);
    gsl_vector * c = gsl_vector_alloc (Ndim);

    double dij2[nid];
    int ids[nid];
    double ef_cvh;
    for (i=0;i<n_data;i++) {
        for (j=0;j<nid;j++) {
            dij2[j] = 0.0;
            for (k=0;k<Ndim;k++) {
                tmp = (double) x[i*Ndim+k]*1.0;
                dij2[j] += pow(tmp - convex_h_x[j][k],2);
            }
            ids[j] = j;
        }
        sort_array_ids(dij2,ids,nid);
        if (dij2[0]==0.0)
            ef2cvh[i] = 0.0;
        else {
            for (j=0;j<Ndim;j++) {
                for (k=0;k<Ndim;k++)
                    gsl_matrix_set (A,k,j,convex_h_x[ids[j]][k]);
                gsl_vector_set(b,j,x[i*Ndim+j]*1.0);
            }
            gsl_linalg_SV_decomp(A, V, S, work);
            gsl_linalg_SV_solve (A, V, S, b, c);
            ef_cvh = 0.0;
            for (j=0;j<Ndim;j++)
                ef_cvh += gsl_vector_get(c,j)*convex_h_y[ids[j]];
            ef2cvh[i] = y[i]-ef_cvh;
        }
    }

    gsl_matrix_free(A);
    gsl_matrix_free(V);
    gsl_vector_free(S);
    gsl_vector_free(work);
    gsl_vector_free(b);
    gsl_vector_free(c);

}


text:
/* err/gsl_errno.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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

#ifndef __GSL_ERRNO_H__
#define __GSL_ERRNO_H__

#include <stdio.h>
#include <errno.h>
#include <gsl/gsl_types.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

enum {
  GSL_SUCCESS  = 0,
  GSL_FAILURE  = -1,
  GSL_CONTINUE = -2,  /* iteration has not converged */
  GSL_EDOM     = 1,   /* input domain error, e.g sqrt(-1) */
  GSL_ERANGE   = 2,   /* output range error, e.g. exp(1e100) */
  GSL_EFAULT   = 3,   /* invalid pointer */
  GSL_EINVAL   = 4,   /* invalid argument supplied by user */
  GSL_EFAILED  = 5,   /* generic failure */
  GSL_EFACTOR  = 6,   /* factorization failed */
  GSL_ESANITY  = 7,   /* sanity check failed - shouldn't happen */
  GSL_ENOMEM   = 8,   /* malloc failed */
  GSL_EBADFUNC = 9,   /* problem with user-supplied function */
  GSL_ERUNAWAY = 10,  /* iterative process is out of control */
  GSL_EMAXITER = 11,  /* exceeded max number of iterations */
  GSL_EZERODIV = 12,  /* tried to divide by zero */
  GSL_EBADTOL  = 13,  /* user specified an invalid tolerance */
  GSL_ETOL     = 14,  /* failed to reach the specified tolerance */
  GSL_EUNDRFLW = 15,  /* underflow */
  GSL_EOVRFLW  = 16,  /* overflow  */
  GSL_ELOSS    = 17,  /* loss of accuracy */
  GSL_EROUND   = 18,  /* failed because of roundoff error */
  GSL_EBADLEN  = 19,  /* matrix, vector lengths are not conformant */
  GSL_ENOTSQR  = 20,  /* matrix not square */
  GSL_ESING    = 21,  /* apparent singularity detected */
  GSL_EDIVERGE = 22,  /* integral or series is divergent */
  GSL_EUNSUP   = 23,  /* requested feature is not supported by the hardware */
  GSL_EUNIMPL  = 24,  /* requested feature not (yet) implemented */
  GSL_ECACHE   = 25,  /* cache limit exceeded */
  GSL_ETABLE   = 26,  /* table limit exceeded */
  GSL_ENOPROG  = 27,  /* iteration is not making progress towards solution */
  GSL_ENOPROGJ = 28,  /* jacobian evaluations are not improving the solution */
  GSL_ETOLF    = 29,  /* cannot reach the specified tolerance in F */
  GSL_ETOLX    = 30,  /* cannot reach the specified tolerance in X */
  GSL_ETOLG    = 31,  /* cannot reach the specified tolerance in gradient */
  GSL_EOF      = 32   /* end of file */
} ;

void gsl_error (const char * reason, const char * file, int line,
                int gsl_errno);

void gsl_stream_printf (const char *label, const char *file,
                        int line, const char *reason);

const char * gsl_strerror (const int gsl_errno);

typedef void gsl_error_handler_t (const char * reason, const char * file,
                                  int line, int gsl_errno);

typedef void gsl_stream_handler_t (const char * label, const char * file,
                                   int line, const char * reason);

gsl_error_handler_t *
gsl_set_error_handler (gsl_error_handler_t * new_handler);

gsl_error_handler_t *
gsl_set_error_handler_off (void);

gsl_stream_handler_t *
gsl_set_stream_handler (gsl_stream_handler_t * new_handler);

FILE * gsl_set_stream (FILE * new_stream);

/* GSL_ERROR: call the error handler, and return the error code */

#define GSL_ERROR(reason, gsl_errno) \
       do { \
       gsl_error (reason, __FILE__, __LINE__, gsl_errno) ; \
       return gsl_errno ; \
       } while (0)

/* GSL_ERROR_VAL: call the error handler, and return the given value */

#define GSL_ERROR_VAL(reason, gsl_errno, value) \
       do { \
       gsl_error (reason, __FILE__, __LINE__, gsl_errno) ; \
       return value ; \
       } while (0)

/* GSL_ERROR_VOID: call the error handler, and then return
   (for void functions which still need to generate an error) */

#define GSL_ERROR_VOID(reason, gsl_errno) \
       do { \
       gsl_error (reason, __FILE__, __LINE__, gsl_errno) ; \
       return ; \
       } while (0)

/* GSL_ERROR_NULL suitable for out-of-memory conditions */

#define GSL_ERROR_NULL(reason, gsl_errno) GSL_ERROR_VAL(reason, gsl_errno, 0)

/* Sometimes you have several status results returned from
 * function calls and you want to combine them in some sensible
 * way. You cannot produce a "total" status condition, but you can
 * pick one from a set of conditions based on an implied hierarchy.
 *
 * In other words:
 *    you have: status_a, status_b, ...
 *    you want: status = (status_a if it is bad, or status_b if it is bad,...)
 *
 * In this example you consider status_a to be more important and
 * it is checked first, followed by the others in the order specified.
 *
 * Here are some dumb macros to do this.
 */
#define GSL_ERROR_SELECT_2(a,b)       ((a) != GSL_SUCCESS ? (a) : ((b) != GSL_SUCCESS ? (b) : GSL_SUCCESS))
#define GSL_ERROR_SELECT_3(a,b,c)     ((a) != GSL_SUCCESS ? (a) : GSL_ERROR_SELECT_2(b,c))
#define GSL_ERROR_SELECT_4(a,b,c,d)   ((a) != GSL_SUCCESS ? (a) : GSL_ERROR_SELECT_3(b,c,d))
#define GSL_ERROR_SELECT_5(a,b,c,d,e) ((a) != GSL_SUCCESS ? (a) : GSL_ERROR_SELECT_4(b,c,d,e))

#define GSL_STATUS_UPDATE(sp, s) do { if ((s) != GSL_SUCCESS) *(sp) = (s);} while(0)

__END_DECLS

#endif /* __GSL_ERRNO_H__ */


text:
/*
 *  vbHmmGauss.c
 *  Model-specific core functions for VB-HMM-GAUSS.
 *
 *  Created by OKAMOTO Kenji, SAKO Yasushi and RIKEN
 *  Copyright 2011-2015
 *  Cellular Informatics Laboratory, Advance Science Institute, RIKEN, Japan.
 *  All rights reserved.
 *
 *  Ver. 1.1.0
 *  Last modified on 2016.11.04
 */

#include "vbHmmGauss.h"
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_psi.h>
#include <string.h>
#include "rand.h"

#ifdef _OPENMP
#include "omp.h"
#endif

#define  MAX(a,b)  ((a)>(b)?(a):(b))
#define  MIN(a,b)  ((a)<(b)?(a):(b))

static int isGlobalAnalysis = 0;

void setFunctions_gauss(){
    commonFunctions funcs;
    funcs.newModelParameters    = newModelParameters_gauss;
    funcs.freeModelParameters   = freeModelParameters_gauss;
    funcs.newModelStats         = newModelStats_gauss;
    funcs.freeModelStats        = freeModelStats_gauss;
    funcs.initializeVbHmm       = initializeVbHmm_gauss;
    funcs.pTilde_z1             = pTilde_z1_gauss;
    funcs.pTilde_zn_zn1         = pTilde_zn_zn1_gauss;
    funcs.pTilde_xn_zn          = pTilde_xn_zn_gauss;
    funcs.calcStatsVars         = calcStatsVars_gauss;
    funcs.maximization          = maximization_gauss;
    funcs.varLowerBound         = varLowerBound_gauss;
    funcs.reorderParameters     = reorderParameters_gauss;
    funcs.outputResults         = outputResults_gauss;
    setFunctions( funcs );
}

void setGFunctions_gauss(){
    gCommonFunctions funcs;
    funcs.newModelParameters    = newModelParameters_gauss;
    funcs.freeModelParameters   = freeModelParameters_gauss;
    funcs.newModelStats         = newModelStats_gauss;
    funcs.freeModelStats        = freeModelStats_gauss;
    funcs.newModelStatsG        = newModelStatsG_gauss;
    funcs.freeModelStatsG       = freeModelStatsG_gauss;
    funcs.initializeVbHmmG      = initializeVbHmmG_gauss;
    funcs.pTilde_z1             = pTilde_z1_gauss;
    funcs.pTilde_zn_zn1         = pTilde_zn_zn1_gauss;
    funcs.pTilde_xn_zn          = pTilde_xn_zn_gauss;
    funcs.calcStatsVarsG        = calcStatsVarsG_gauss;
    funcs.maximizationG         = maximizationG_gauss;
    funcs.varLowerBoundG        = varLowerBoundG_gauss;
    funcs.reorderParametersG    = reorderParametersG_gauss;
    funcs.outputResultsG        = outputResultsG_gauss;
    setGFunctions( funcs );
    isGlobalAnalysis = 1;
}

void outputResults_gauss( xn, gv, iv, logFP )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
FILE *logFP;
{
    outputGaussResults( xn, gv, iv, logFP );
}

void outputResultsG_gauss( xns, gv, ivs, logFP )
xnDataBundle *xns;
globalVars *gv;
indVarBundle *ivs;
FILE *logFP;
{
    outputGaussResultsG( xns, gv, ivs, logFP );
}

void *newModelParameters_gauss( xn, sNo )
xnDataSet *xn;
int sNo;
{
    int i;
    gaussParameters *p = (void*)malloc( sizeof(gaussParameters) );

    p->uPiArr = (double*)malloc( sNo * sizeof(double) );
    p->sumUPi = 0.0;
    p->uAMat = (double**)malloc( sNo * sizeof(double*) );
    for( i = 0 ; i < sNo ; i++ ){
        p->uAMat[i] = (double*)malloc( sNo * sizeof(double) );
    }
    p->sumUAArr = (double*)malloc( sNo * sizeof(double) );

    p->avgPi = (double *)malloc( sNo * sizeof(double) );
    p->avgLnPi = (double *)malloc( sNo * sizeof(double) );
    p->avgA = (double **)malloc( sNo * sizeof(double*) );
    p->avgLnA = (double **)malloc( sNo * sizeof(double*) );
    for( i = 0 ; i < sNo ; i++ ){
        p->avgA[i] = (double *)malloc( sNo * sizeof(double) );
        p->avgLnA[i] = (double *)malloc( sNo * sizeof(double) );
    }
    p->avgMu = (double *)malloc( sNo * sizeof(double) );
    p->avgLm = (double *)malloc( sNo * sizeof(double) );
    p->avgLnLm = (double *)malloc( sNo * sizeof(double) );

    p->uBtArr = (double *)malloc( sNo * sizeof(double) );
    p->uMuArr = (double *)malloc( sNo * sizeof(double) );
    p->uAArr = (double *)malloc( sNo * sizeof(double) );
    p->uBArr = (double *)malloc( sNo * sizeof(double) );
    p->btMu = (double *)malloc( sNo * sizeof(double) );
    p->aLm = (double *)malloc( sNo * sizeof(double) );
    p->bLm = (double *)malloc( sNo * sizeof(double) );
    p->mu0 = (double *)malloc( sNo * sizeof(double) );

    return p;
}

void freeModelParameters_gauss( p, xn, sNo )
void **p;
xnDataSet *xn;
int sNo;
{
    gaussParameters *gp = *p;
    int i;

    free( gp->uPiArr );
    for( i = 0 ; i < sNo ; i++ ){
        free( gp->uAMat[i] );
    }
    free( gp->uAMat );
    free( gp->sumUAArr );

    free( gp->avgPi );
    free( gp->avgLnPi );
    for( i = 0 ; i < sNo ; i++ ){
        free( gp->avgA[i] );
        free( gp->avgLnA[i] );
    }
    free( gp->avgA );
    free( gp->avgLnA );
    free( gp->avgMu );
    free( gp->avgLm );
    free( gp->avgLnLm );

    free( gp->uBtArr );
    free( gp->uMuArr );
    free( gp->uAArr );
    free( gp->uBArr );
    free( gp->btMu );
    free( gp->aLm );
    free( gp->bLm );
    free( gp->mu0 );

    free( *p );
    *p = NULL;
}

void *newModelStats_gauss( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    if( isGlobalAnalysis == 0 ){
        int sNo = gv->sNo;
        gaussStats *s = (gaussStats*)malloc( sizeof(gaussStats) );

        int i;
        s->Ni = (double *)malloc( sNo * sizeof(double) );
        s->Nij = (double **)malloc( sNo * sizeof(double*) );
        for( i = 0 ; i < sNo ; i++ )
        {   s->Nij[i] = (double *)malloc( sNo * sizeof(double) );   }
        s->Nii = (double *)malloc( sNo * sizeof(double) );
        s->barX = (double *)malloc( sNo * sizeof(double) );
        s->NiSi = (double *)malloc( sNo * sizeof(double) );

        return s;

    } else {

        return NULL;

    }
}

void freeModelStats_gauss( s, xn, gv, iv )
void **s;
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    if( isGlobalAnalysis == 0 ){
        int sNo = gv->sNo;
        gaussStats *gs = *s;
        int i;
        free( gs->Ni );
        for( i = 0 ; i < sNo ; i++ )
        {   free( gs->Nij[i] );   }
        free( gs->Nij );
        free( gs->Nii );
        free( gs->barX );
        free( gs->NiSi );

        free( gs );
        *s = NULL;
    }
}

void *newModelStatsG_gauss( xns, gv, ivs)
xnDataBundle *xns;
globalVars *gv;
indVarBundle *ivs;
{
    int sNo = gv->sNo;
    gaussGlobalStats *gs = (gaussGlobalStats*)malloc( sizeof(gaussGlobalStats) );

    int i;
    gs->NiR = (double *)malloc( sNo * sizeof(double) );
    gs->NijR = (double **)malloc( sNo * sizeof(double*) );
    for( i = 0 ; i < sNo ; i++ )
    {   gs->NijR[i] = (double *)malloc( sNo * sizeof(double) );   }
    gs->NiiR = (double *)malloc( sNo * sizeof(double) );
    gs->barXR = (double *)malloc( sNo * sizeof(double) );
    gs->NiSiR = (double *)malloc( sNo * sizeof(double) );
    gs->z1iR = (double *)malloc( sNo * sizeof(double) );

    return gs;
}

void freeModelStatsG_gauss( gs, xns, gv, ivs )
void **gs;
xnDataBundle *xns;
globalVars *gv;
indVarBundle *ivs;
{
    int sNo = gv->sNo;
    gaussGlobalStats *ggs = *gs;
    int i;
    free( ggs->NiR );
    for( i = 0 ; i < sNo ; i++ )
    {   free( ggs->NijR[i] );   }
    free( ggs->NijR );
    free( ggs->NiiR );
    free( ggs->barXR );
    free( ggs->NiSiR );
    free( ggs->z1iR );

    free( *gs );
    *gs = NULL;
}

void initializeVbHmm_gauss( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    gaussData *d = xn->data;
    size_t dLen = xn->N;
    int sNo = gv->sNo;
    gaussParameters *p = gv->params;

    int i, j;
    double totalX = 0.0, precX, varX = 0.0;
    for( i = 0 ; i < dLen ; i++ ){
        totalX += d->v[i];
    }
    double meanX = totalX / (double)dLen;
    for( i = 0 ; i < dLen ; i++ ){
        varX += pow(d->v[i] - meanX, 2.0);
    }
    varX /= (double)(dLen - 1);
    precX = 1.0 / varX;

    // hyper parameter for p( pi(i) )
    p->sumUPi = 0.0;
    for( i = 0 ; i < sNo ; i++ ){
        p->uPiArr[i] = 1.0;
        p->sumUPi += p->uPiArr[i];
    }

    // hyper parameter for p( A(i,j) )
    for( i = 0 ; i < sNo ; i++ ){
        p->sumUAArr[i] = 0.0;
        for( j = 0 ; j < sNo ; j++ ){
            if( j == i ){
                p->uAMat[i][j] = 5.0;
            } else {
                p->uAMat[i][j] = 1.0;
            }
            p->sumUAArr[i] += p->uAMat[i][j];
        }
    }

    // hyper parameter for p( mu(k), lm(k) )
    for( i = 0 ; i < sNo ; i++ ){
        p->uBtArr[i] = 0.25;
        p->uMuArr[i] = meanX;
        p->uAArr[i] = 0.01 * varX;
        p->uBArr[i] = 0.01;
    }

    initialize_indVars_gauss( xn, gv, iv );

    calcStatsVars_gauss( xn, gv, iv );
    maximization_gauss( xn, gv, iv );
}

void initializeVbHmmG_gauss( xns, gv, ivs )
xnDataBundle *xns;
globalVars *gv;
indVarBundle *ivs;
{
    gaussData *d;
    size_t dLen, totalN;
    int sNo = gv->sNo, rNo = xns->R;
    gaussParameters *p = gv->params;

    int i, j, r;
    double totalX, precX, varX, meanX;

    totalX = 0.0;
    totalN = 0;
    for( r = 0 ; r < rNo ; r++ ){
        d = xns->xn[r]->data;
        dLen = xns->xn[r]->N;
        for( i = 0 ; i < dLen ; i++ ){
            totalX += d->v[i];
        }
        totalN += dLen;
    }
    meanX = totalX / (double)totalN;
    varX = 0.0;
    for( r = 0 ; r < rNo ; r++ ){
        d = xns->xn[r]->data;
        dLen = xns->xn[r]->N;
        for( i = 0 ; i < dLen ; i++ ){
            varX += pow(d->v[i] - meanX, 2.0);
        }
    }
    varX /= (double)(totalN - 1);
    precX = 1.0 / varX;

    p->sumUPi = 0.0;
    for( i = 0 ; i < sNo ; i++ ){
        p->uPiArr[i] = 1.0;
        p->sumUPi += p->uPiArr[i];
    }

    for( i = 0 ; i < sNo ; i++ ){
        p->sumUAArr[i] = 0.0;
        for( j = 0 ; j < sNo ; j++ ){
            if( j == i ){
                p->uAMat[i][j] = 5.0;
            } else {
                p->uAMat[i][j] = 1.0;
            }
            p->sumUAArr[i] += p->uAMat[i][j];
        }
    }

    // hyper parameter for p( mu(k), lm(k) )
    for( i = 0 ; i < sNo ; i++ ){
        p->uBtArr[i] = 0.25;
        p->uMuArr[i] = meanX;
        p->uAArr[i] = 0.01 * varX;
        p->uBArr[i] = 0.01;
    }

    for( r = 0 ; r < rNo ; r++ ){
        initialize_indVars_gauss( xns->xn[r], gv, ivs->indVars[r] );
    }

    calcStatsVarsG_gauss( xns, gv, ivs );
    maximizationG_gauss( xns, gv, ivs );
}

void initialize_indVars_gauss( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    size_t dLen = xn->N;
    int sNo = gv->sNo;
    double **gmMat = iv->gmMat;

    int i;
    size_t n;
    double sumPar;
    for( n = 0 ; n < dLen ; n++ ){
        sumPar = 0.0;
        for( i = 0 ; i < sNo ; i++ ){
            gmMat[n][i] = enoise(1.0) + 1.0;
            sumPar += gmMat[n][i];
        }
        for( i = 0 ; i < sNo ; i++ ){
            gmMat[n][i] /= sumPar;
        }
    }
}

xnDataSet *newXnDataSet_gauss( filename )
const char *filename;
{
    xnDataSet *xn = (xnDataSet*)malloc( sizeof(xnDataSet) );
    xn->name = (char*)malloc( strlen(filename) + 2 );
    strncpy( xn->name, filename, strlen(filename)+1 );
    xn->data = (gaussData*)malloc( sizeof(gaussData) );
    gaussData *d = (gaussData*)xn->data;
    d->v = NULL;
    return xn;
}

void freeXnDataSet_gauss( xn )
xnDataSet **xn;
{
    gaussData *d = (gaussData*)(*xn)->data;
    free( d->v );
    free( (*xn)->data );
    free( (*xn)->name );
    free( *xn );
    *xn = NULL;
}

double pTilde_z1_gauss( i, params )
int i;
void *params;
{
    gaussParameters *p = (gaussParameters*)params;
    return exp( p->avgLnPi[i] );
}

double pTilde_zn_zn1_gauss( i, j, params )
int i, j;
void *params;
{
    gaussParameters *p = (gaussParameters*)params;
    return exp( p->avgLnA[i][j] );
}

double pTilde_xn_zn_gauss( xn, n, i, params )
xnDataSet *xn;
size_t n;
int i;
void *params;
{
    gaussParameters *p = (gaussParameters*)params;
    gaussData *d = (gaussData*)xn->data;
    double val;
    val  = p->avgLnLm[i] - log(2.0 * M_PI);
    val -= 1.0/p->btMu[i] + p->aLm[i] / p->bLm[i] * pow( d->v[n] - p->avgMu[i], 2.0);
    return exp(val / 2.0);
}

void calcStatsVars_gauss( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    gaussData *d = (gaussData*)xn->data;
    gaussStats *s = (gaussStats*)iv->stats;
    size_t dLen = xn->N;
    int sNo = gv->sNo;
    double **gmMat = iv->gmMat, ***xiMat = iv->xiMat;
    double *Nii = s->Nii, **Nij = s->Nij, *Ni = s->Ni;
    double *barX = s->barX, *NiSi = s->NiSi;
    size_t n;
    int i, j;

    for( i = 0 ; i < sNo ; i++ ){
        Ni[i]   = 1e-10;
        Nii[i]  = 1e-10;
        barX[i] = 1e-10;
        NiSi[i] = 1e-10;
        for( j = 0 ; j < sNo ; j++ ){
            Nij[i][j] = 1e-10;
        }

        for( n = 0 ; n < dLen ; n++ ){
            Ni[i]   += gmMat[n][i];
            barX[i] += gmMat[n][i] * d->v[n];
            for( j = 0 ; j < sNo ; j++ ){
                Nii[i]    += xiMat[n][i][j];
                Nij[i][j] += xiMat[n][i][j];
            }
        }
        barX[i] /= Ni[i];
        for( n = 0 ; n < dLen ; n++ ){
            NiSi[i] += gmMat[n][i] * pow( d->v[n] - barX[i], 2.0);
        }
    }
}

void calcStatsVarsG_gauss( xns, gv, ivs )
xnDataBundle *xns;
globalVars *gv;
indVarBundle *ivs;
{
    gaussData *d;
    gaussGlobalStats *gs = (gaussGlobalStats*)ivs->stats;
    int sNo = gv->sNo, rNo = xns->R;
    double **gmMat, ***xiMat;
    double *NiiR = gs->NiiR, **NijR = gs->NijR, *NiR = gs->NiR;
    double *barXR = gs->barXR, *NiSiR = gs->NiSiR, *z1iR = gs->z1iR;
    size_t dLen, n;
    int i, j, r;

    for( i = 0 ; i < sNo ; i++ ){
        NiR[i]   = 1e-10;
        NiiR[i]  = 1e-10;
        barXR[i] = 1e-10;
        NiSiR[i] = 1e-10;
        for( j = 0 ; j < sNo ; j++ ){
            NijR[i][j] = 1e-10;
        }
        z1iR[i] = 1e-10;
    }
    for( r = 0 ; r < rNo ; r++ ){
        d = (gaussData*)xns->xn[r]->data;
        dLen = xns->xn[r]->N;
        gmMat = ivs->indVars[r]->gmMat;
        xiMat = ivs->indVars[r]->xiMat;
        for( i = 0 ; i < sNo ; i++ ){
            z1iR[i] += gmMat[0][i];
            for( n = 0 ; n < dLen ; n++ ){
                NiR[i]   += gmMat[n][i];
                barXR[i] += gmMat[n][i] * d->v[n];
                for( j = 0 ; j < sNo ; j++ ){
                    NiiR[i]  += xiMat[n][i][j];
                    NijR[i][j] += xiMat[n][i][j];
                }
            }
        }
    }
    for( i = 0 ; i < sNo ; i++ ){
        barXR[i] /= NiR[i];
        for( r = 0 ; r < rNo ; r++ ){
            d = (gaussData*)xns->xn[r]->data;
            dLen = xns->xn[r]->N;
            gmMat = ivs->indVars[r]->gmMat;
            for( n = 0 ; n < dLen ; n++ ){
                NiSiR[i] += gmMat[n][i] * pow( d->v[n] - barXR[i], 2.0);
            }
        }
    }
}

void maximization_gauss( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    gaussParameters *p = (gaussParameters*)gv->params;
    gaussStats *s = (gaussStats*)iv->stats;
    int sNo = gv->sNo;
    double **gmMat = iv->gmMat;
    double *uPiArr = p->uPiArr, sumUPi = p->sumUPi;
    double **uAMat = p->uAMat, *sumUAArr = p->sumUAArr;
    double *uBtArr = p->uBtArr, *uMuArr = p->uMuArr, *uAArr = p->uAArr, *uBArr = p->uBArr;
    double *avgPi = p->avgPi, *avgLnPi = p->avgLnPi, **avgA = p->avgA, **avgLnA = p->avgLnA;
    double *avgMu = p->avgMu, *avgLm = p->avgLm, *avgLnLm = p->avgLnLm;
    double *mu0 = p->mu0, *btMu = p->btMu, *aLm = p->aLm, *bLm = p->bLm;
    double *Ni = s->Ni, *Nii = s->Nii, **Nij = s->Nij, *barX = s->barX, *NiSi = s->NiSi;
    int i, j;

    for( i = 0 ; i < sNo ; i++ ){
        avgPi[i] = ( uPiArr[i] + gmMat[0][i] ) / ( sumUPi + 1.0 );
        avgLnPi[i] = gsl_sf_psi( uPiArr[i] + gmMat[0][i] ) - gsl_sf_psi( sumUPi + 1.0 );

        for( j = 0 ; j < sNo ; j++ ){
            avgA[i][j] = ( uAMat[i][j] + Nij[i][j] ) / ( sumUAArr[i] + Nii[i] );
            avgLnA[i][j] = gsl_sf_psi( uAMat[i][j] + Nij[i][j] ) - gsl_sf_psi( sumUAArr[i] + Nii[i] );
        }

        btMu[i] = uBtArr[i] + Ni[i];
        mu0[i]  = (uBtArr[i] * uMuArr[i] + Ni[i] * barX[i]) / btMu[i];
        aLm[i]  = uAArr[i] + Ni[i] / 2.0;
        bLm[i]  = uBArr[i] + (NiSi[i] / 2.0);
        bLm[i] += uBtArr[i] * Ni[i] * pow( barX[i] - uMuArr[i], 2.0) / 2.0 / (uBtArr[i] + Ni[i]);

        avgMu[i]    = mu0[i];
        avgLm[i]    = aLm[i] / bLm[i];
        avgLnLm[i]  = gsl_sf_psi( aLm[i] ) - log( bLm[i] );
    }
}

void maximizationG_gauss( xns, gv, ivs )
xnDataBundle *xns;
globalVars *gv;
indVarBundle *ivs;
{
    gaussParameters *p = (gaussParameters*)gv->params;
    double *uPiArr = p->uPiArr, sumUPi = p->sumUPi;
    double **uAMat = p->uAMat, *sumUAArr = p->sumUAArr;
    double *uBtArr = p->uBtArr, *uMuArr = p->uMuArr, *uAArr = p->uAArr, *uBArr = p->uBArr;
    double *avgPi = p->avgPi, *avgLnPi = p->avgLnPi, **avgA = p->avgA, **avgLnA = p->avgLnA;
    double *avgMu = p->avgMu, *avgLm = p->avgLm, *avgLnLm = p->avgLnLm;
    double *mu0 = p->mu0, *btMu = p->btMu, *aLm = p->aLm, *bLm = p->bLm;
    gaussGlobalStats *gs = (gaussGlobalStats*)ivs->stats;
    double *NiR = gs->NiR, *NiiR = gs->NiiR, **NijR = gs->NijR, *barXR = gs->barXR, *NiSiR = gs->NiSiR;
    double *z1iR = gs->z1iR, dR = (double)(xns->R);
    int sNo = gv->sNo;
    int i, j;

    for( i = 0 ; i < sNo ; i++ ){
        avgPi[i] = ( uPiArr[i] + z1iR[i] ) / ( sumUPi + dR );
        avgLnPi[i] = gsl_sf_psi( uPiArr[i] + z1iR[i] ) - gsl_sf_psi( sumUPi + dR );

        for( j = 0 ; j < sNo ; j++ ){
            avgA[i][j] = ( uAMat[i][j] + NijR[i][j] ) / ( sumUAArr[i] + NiiR[i] );
            avgLnA[i][j] = gsl_sf_psi( uAMat[i][j] + NijR[i][j] ) - gsl_sf_psi( sumUAArr[i] + NiiR[i] );
        }

        btMu[i] = uBtArr[i] + NiR[i];
        mu0[i]  = (uBtArr[i] * uMuArr[i] + NiR[i] * barXR[i]) / btMu[i];
        aLm[i]  = uAArr[i] + NiR[i] / 2.0;
        bLm[i]  = uBArr[i] + (NiSiR[i] / 2.0);
        bLm[i] += uBtArr[i] * NiR[i] * pow( barXR[i] - uMuArr[i], 2.0) / 2.0 / (uBtArr[i] + NiR[i]);

        avgMu[i]    = mu0[i];
        avgLm[i]    = aLm[i] / bLm[i];
        avgLnLm[i]  = gsl_sf_psi( aLm[i] ) - log( bLm[i] );
    }
}

double varLowerBound_gauss( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    gaussParameters *p = (gaussParameters*)gv->params;
    gaussStats *s = (gaussStats*)iv->stats;
    size_t dLen = xn->N;
    int sNo = gv->sNo;
    double **gmMat = iv->gmMat, *cn = iv->cn;
    double *uPiArr = p->uPiArr, sumUPi = p->sumUPi;
    double **uAMat = p->uAMat, *sumUAArr = p->sumUAArr;
    double *uBtArr = p->uBtArr, *uMuArr = p->uMuArr, *uAArr = p->uAArr, *uBArr = p->uBArr;
    double *avgLnPi = p->avgLnPi, **avgLnA = p->avgLnA;
    double *avgLm = p->avgLm, *avgLnLm = p->avgLnLm;
    double *mu0 = p->mu0, *btMu = p->btMu, *aLm = p->aLm, *bLm = p->bLm;
    double *Nii = s->Nii, **Nij = s->Nij;
    size_t n;
    int i, j;

    double lnpPi = gsl_sf_lngamma(sumUPi);
    double lnpA = 0.0;
    double lnpMuLm = 0.0;
    double lnqPi = gsl_sf_lngamma(sumUPi + 1.0);
    double lnqA = 0.0;
    double lnqMuLm = - sNo / 2.0;
    for( i = 0 ; i < sNo ; i++ ){
        lnpPi += (uPiArr[i]-1.0) * avgLnPi[i] - gsl_sf_lngamma(uPiArr[i]);

        lnpMuLm += log(uBtArr[i]) / 2.0;
        lnpMuLm -= uBtArr[i] / 2.0 *( 1.0 / btMu[i] + aLm[i] / bLm[i] * pow(mu0[i] - uMuArr[i], 2.0) );
        lnpMuLm += - gsl_sf_lngamma(uAArr[i]) + uAArr[i] * log(uBArr[i]);
        lnpMuLm += (uAArr[i] - 0.5) * avgLnLm[i] - uBArr[i] * avgLm[i];

        lnqPi += (uPiArr[i]+gmMat[0][i]-1.0) * (gsl_sf_psi(uPiArr[i]+gmMat[0][i]) - gsl_sf_psi(sumUPi+1.0));
        lnqPi -= gsl_sf_lngamma(uPiArr[i] + gmMat[0][i]);

        lnpA += gsl_sf_lngamma(sumUAArr[i]);
        lnqA += gsl_sf_lngamma(sumUAArr[i] + Nii[i]);
        for( j = 0 ; j < sNo ; j++ ){
            lnpA += (uAMat[i][j]-1.0)*avgLnA[i][j] - gsl_sf_lngamma(uAMat[i][j]);

            lnqA += (uAMat[i][j] + Nij[i][j] - 1.0) * (gsl_sf_psi(uAMat[i][j]+Nij[i][j]) - gsl_sf_psi(sumUAArr[i]+Nii[i]));
            lnqA -= gsl_sf_lngamma( uAMat[i][j] + Nij[i][j] );
        }

        lnqMuLm += log(btMu[i]) / 2.0 - gsl_sf_lngamma(aLm[i]) + aLm[i] * log(bLm[i]);
        lnqMuLm += (aLm[i] - 0.5) * avgLnLm[i] - aLm[i];
    }

    double lnpX = 0.0;
    for( n = 0 ; n < dLen ; n++ ){
        lnpX += log( cn[n] );
    }

    double val;
    val  = lnpPi + lnpA + lnpMuLm;
    val -= lnqPi + lnqA + lnqMuLm;
    val += lnpX;
    val += log(gsl_sf_fact(sNo));

    return val;
}

double varLowerBoundG_gauss( xns, gv, ivs )
xnDataBundle *xns;
globalVars *gv;
indVarBundle *ivs;
{
    int sNo = gv->sNo, rNo = xns->R;
    gaussParameters *p = (gaussParameters*)gv->params;
    double *uPiArr = p->uPiArr, sumUPi = p->sumUPi;
    double **uAMat = p->uAMat, *sumUAArr = p->sumUAArr;
    double *uBtArr = p->uBtArr, *uMuArr = p->uMuArr, *uAArr = p->uAArr, *uBArr = p->uBArr;
    double *avgLnPi = p->avgLnPi, **avgLnA = p->avgLnA;
    double *avgLm = p->avgLm, *avgLnLm = p->avgLnLm;
    double *mu0 = p->mu0, *btMu = p->btMu, *aLm = p->aLm, *bLm = p->bLm;
    gaussGlobalStats *gs = (gaussGlobalStats*)ivs->stats;
    double *NiiR = gs->NiiR, **NijR = gs->NijR, *z1iR = gs->z1iR, dR = (double)xns->R;
    size_t n;
    int i, j, r;

    double lnpPi = gsl_sf_lngamma(sumUPi);
    double lnpA = 0.0;
    double lnpMuLm = 0.0;
    double lnqPi = gsl_sf_lngamma(sumUPi + dR);
    double lnqA = 0.0;
    double lnqMuLm = - sNo / 2.0;
    for( i = 0 ; i < sNo ; i++ ){
        lnpPi += (uPiArr[i]-1.0) * avgLnPi[i] - gsl_sf_lngamma(uPiArr[i]);

        lnpMuLm += log(uBtArr[i]) / 2.0;
        lnpMuLm -= uBtArr[i] / 2.0 *( 1.0 / btMu[i] + aLm[i] / bLm[i] * pow(mu0[i] - uMuArr[i], 2.0) );
        lnpMuLm += - gsl_sf_lngamma(uAArr[i]) + uAArr[i] * log(uBArr[i]);
        lnpMuLm += (uAArr[i] - 0.5) * avgLnLm[i] - uBArr[i] * avgLm[i];

        lnqPi += (uPiArr[i]+z1iR[i]-1.0) * (gsl_sf_psi(uPiArr[i]+z1iR[i]) - gsl_sf_psi(sumUPi+dR));
        lnqPi -= gsl_sf_lngamma(uPiArr[i] + z1iR[i]);

        lnpA += gsl_sf_lngamma(sumUAArr[i]);
        lnqA += gsl_sf_lngamma(sumUAArr[i] + NiiR[i]);
        for( j = 0 ; j < sNo ; j++ ){
            lnpA += (uAMat[i][j]-1.0)*avgLnA[i][j] - gsl_sf_lngamma(uAMat[i][j]);

            lnqA += (uAMat[i][j] + NijR[i][j] - 1.0) * (gsl_sf_psi(uAMat[i][j]+NijR[i][j]) - gsl_sf_psi(sumUAArr[i]+NiiR[i]));
            lnqA -= gsl_sf_lngamma( uAMat[i][j] + NijR[i][j] );
        }

        lnqMuLm += log(btMu[i]) / 2.0 - gsl_sf_lngamma(aLm[i]) + aLm[i] * log(bLm[i]);
        lnqMuLm += (aLm[i] - 0.5) * avgLnLm[i] - aLm[i];
    }

    double lnpX = 0.0;
    for( r = 0 ; r < rNo ; r++ ){
        size_t dLen = xns->xn[r]->N;
        for( n = 0 ; n < dLen ; n++ ){
            lnpX += log( ivs->indVars[r]->cn[n] );
        }
    }

    double val;
    val  = lnpPi + lnpA + lnpMuLm;
    val -= lnqPi + lnqA + lnqMuLm;
    val += lnpX;
    val += log(gsl_sf_fact(sNo));

    return val;
}

void reorderParameters_gauss( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    gaussParameters *p = (gaussParameters*)gv->params;
    gaussStats *s = (gaussStats*)iv->stats;
    size_t dLen = xn->N;
    int sNo = gv->sNo;
    double **gmMat = iv->gmMat, ***xiMat = iv->xiMat;
    double *avgPi = p->avgPi, *avgLnPi = p->avgLnPi, **avgA = p->avgA, **avgLnA = p->avgLnA;
    double *avgMu = p->avgMu, *avgLm = p->avgLm, *avgLnLm = p->avgLnLm;
    double *btMu = p->btMu, *aLm = p->aLm, *bLm = p->bLm;
    double *Ni = s->Ni;
    size_t n;
    int i, j;

    int *index = (int*)malloc( sNo * sizeof(int) );
    double *store = (double*)malloc( sNo * sizeof(double) );
    double **s2D = (double**)malloc( sNo * sizeof(double*) );
    for( i = 0 ; i < sNo ; i++ )
    {   s2D[i] = (double*)malloc( MAX(sNo,2) * sizeof(double) );   }

    // index indicates order of avgMu values (0=biggest avgMu -- sNo=smallest avgMu).
    for( i = 0 ; i < sNo ; i++ ){
        index[i] = sNo - 1;
        for( j = 0 ; j < sNo ; j++ ){
            if( j != i ){
                if( avgMu[i] < avgMu[j] ){
                    index[i]--;
                } else if( avgMu[i] == avgMu[j] ){
                    if( j > i )
                    {   index[i]--;   }
                }
            }
        }
    }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = avgPi[i];   }
    for( i = 0 ; i < sNo ; i++ ){   avgPi[i] = store[i];   }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = avgLnPi[i];   }
    for( i = 0 ; i < sNo ; i++ ){   avgLnPi[i] = store[i];   }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = avgMu[i];   }
    for( i = 0 ; i < sNo ; i++ ){   avgMu[i] = store[i];   }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = avgLm[i];   }
    for( i = 0 ; i < sNo ; i++ ){   avgLm[i] = store[i];   }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = avgLnLm[i];   }
    for( i = 0 ; i < sNo ; i++ ){   avgLnLm[i] = store[i];   }

    //
    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = btMu[i];   }
    for( i = 0 ; i < sNo ; i++ ){   btMu[i] = store[i];   }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = aLm[i];   }
    for( i = 0 ; i < sNo ; i++ ){   aLm[i] = store[i];   }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = bLm[i];   }
    for( i = 0 ; i < sNo ; i++ ){   bLm[i] = store[i];   }
    //

    for( j = 0 ; j < sNo ; j++ ){
        for( i = 0 ; i < sNo ; i++ ){   s2D[index[i]][index[j]] = avgA[i][j];   }
    }
    for( j = 0 ; j < sNo ; j++ ){
        for( i = 0 ; i < sNo ; i++ ){   avgA[i][j] = s2D[i][j];   }
    }

    for( j = 0 ; j < sNo ; j++ ){
        for( i = 0 ; i < sNo ; i++ ){   s2D[index[i]][index[j]] = avgLnA[i][j];   }
    }
    for( j = 0 ; j < sNo ; j++ ){
        for( i = 0 ; i < sNo ; i++ ){   avgLnA[i][j] = s2D[i][j];   }
    }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = Ni[i];   }
    for( i = 0 ; i < sNo ; i++ ){   Ni[i] = store[i];   }

    for( n = 0 ; n < dLen ; n++ ){
        for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = gmMat[n][i];   }
        for( i = 0 ; i < sNo ; i++ ){   gmMat[n][i] = store[i];   }
    }

    for( n = 0 ; n < dLen ; n++ ){
        for( j = 0 ; j < sNo ; j++ ){
            for( i = 0 ; i < sNo ; i++ ){   s2D[index[i]][index[j]] = xiMat[n][i][j];   }
        }
        for( j = 0 ; j < sNo ; j++ ){
            for( i = 0 ; i < sNo ; i++ ){   xiMat[n][i][j] = s2D[i][j];   }
        }
    }

    for( i = 0 ; i < sNo ; i++ ){   free( s2D[i] );   }
    free( s2D );
    free( store );
    free( index );
}

void reorderParametersG_gauss( xns, gv, ivs )
xnDataBundle *xns;
globalVars *gv;
indVarBundle *ivs;
{
    size_t dLen;
    int sNo = gv->sNo, rNo = xns->R;
    gaussParameters *p = (gaussParameters*)gv->params;
    double *avgPi = p->avgPi, *avgLnPi = p->avgLnPi, **avgA = p->avgA, **avgLnA = p->avgLnA;
    double *avgMu = p->avgMu, *avgLm = p->avgLm, *avgLnLm = p->avgLnLm;
    double *btMu = p->btMu, *aLm = p->aLm, *bLm = p->bLm;
    size_t n;
    int i, j, r;

    int *index = (int*)malloc( sNo * sizeof(int) );
    double *store = (double*)malloc( sNo * sizeof(double) );
    double **s2D = (double**)malloc( sNo * sizeof(double*) );
    for( i = 0 ; i < sNo ; i++ )
    {   s2D[i] = (double*)malloc( MAX(sNo,2) * sizeof(double) );   }

    // index indicates order of avgMu values (0=biggest avgMu -- sNo=smallest avgMu).
    for( i = 0 ; i < sNo ; i++ ){
        index[i] = sNo - 1;
        for( j = 0 ; j < sNo ; j++ ){
            if( j != i ){
                if( avgMu[i] < avgMu[j] ){
                    index[i]--;
                } else if( avgMu[i] == avgMu[j] ){
                    if( j > i )
                    {   index[i]--;   }
                }
            }
        }
    }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = avgPi[i];   }
    for( i = 0 ; i < sNo ; i++ ){   avgPi[i] = store[i];   }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = avgLnPi[i];   }
    for( i = 0 ; i < sNo ; i++ ){   avgLnPi[i] = store[i];   }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = avgMu[i];   }
    for( i = 0 ; i < sNo ; i++ ){   avgMu[i] = store[i];   }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = avgLm[i];   }
    for( i = 0 ; i < sNo ; i++ ){   avgLm[i] = store[i];   }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = avgLnLm[i];   }
    for( i = 0 ; i < sNo ; i++ ){   avgLnLm[i] = store[i];   }

    //
    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = btMu[i];   }
    for( i = 0 ; i < sNo ; i++ ){   btMu[i] = store[i];   }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = aLm[i];   }
    for( i = 0 ; i < sNo ; i++ ){   aLm[i] = store[i];   }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = bLm[i];   }
    for( i = 0 ; i < sNo ; i++ ){   bLm[i] = store[i];   }
    //

    for( j = 0 ; j < sNo ; j++ ){
        for( i = 0 ; i < sNo ; i++ ){   s2D[index[i]][index[j]] = avgA[i][j];   }
    }
    for( j = 0 ; j < sNo ; j++ ){
        for( i = 0 ; i < sNo ; i++ ){   avgA[i][j] = s2D[i][j];   }
    }

    for( j = 0 ; j < sNo ; j++ ){
        for( i = 0 ; i < sNo ; i++ ){   s2D[index[i]][index[j]] = avgLnA[i][j];   }
    }
    for( j = 0 ; j < sNo ; j++ ){
        for( i = 0 ; i < sNo ; i++ ){   avgLnA[i][j] = s2D[i][j];   }
    }

    double *NiR = ((gaussGlobalStats*)ivs->stats)->NiR;
    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = NiR[i];   }
    for( i = 0 ; i < sNo ; i++ ){   NiR[i] = store[i];   }

    for( r = 0 ; r < rNo ; r++ ){
        double **gmMat = ivs->indVars[r]->gmMat;
        dLen = xns->xn[r]->N;
        for( n = 0 ; n < dLen ; n++ ){
            for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = gmMat[n][i];   }
            for( i = 0 ; i < sNo ; i++ ){   gmMat[n][i] = store[i];   }
        }
    }

    for( r = 0 ; r < rNo ; r++ ){
        double ***xiMat = ivs->indVars[r]->xiMat;
        dLen = xns->xn[r]->N;
        for( n = 0 ; n < dLen ; n++ ){
            for( j = 0 ; j < sNo ; j++ ){
                for( i = 0 ; i < sNo ; i++ ){   s2D[index[i]][index[j]] = xiMat[n][i][j];   }
            }
            for( j = 0 ; j < sNo ; j++ ){
                for( i = 0 ; i < sNo ; i++ ){   xiMat[n][i][j] = s2D[i][j];   }
            }
        }
    }

    for( i = 0 ; i < sNo ; i++ ){   free( s2D[i] );   }
    free( s2D );
    free( store );
    free( index );
}

void outputGaussResults( xn, gv, iv, logFP )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
FILE *logFP;
{
    gaussParameters *p = (gaussParameters*)gv->params;
    int sNo = gv->sNo;

    int i, j;
    fprintf(logFP, "  results: K = %d \n", sNo);

    fprintf(logFP, "   means: ( %g", p->avgMu[0]);
    for( i = 1 ; i < sNo ; i++ )
    {   fprintf(logFP, ", %g", p->avgMu[i]);   }
    fprintf(logFP, " ) \n");

    fprintf(logFP, "   lambda: ( %g", p->avgLm[0]);
    for( i = 1 ; i < sNo ; i++ ){
        fprintf(logFP, ", %g", p->avgLm[i]);
    }
    fprintf(logFP, " ) \n");

    fprintf(logFP, "   pi: ( %g", p->avgPi[0]);
    for( i = 1 ; i < sNo ; i++ ){
        fprintf(logFP, ", %g", p->avgPi[i]);
    }
    fprintf(logFP, " ) \n");

    fprintf(logFP, "   A_matrix: [");
    for( i = 0 ; i < sNo ; i++ ){
        fprintf(logFP, " ( %g", p->avgA[i][0]);
        for( j = 1 ; j < sNo ; j++ )
        {   fprintf(logFP, ", %g", p->avgA[i][j]);   }
        fprintf(logFP, ")");
    }
    fprintf(logFP, " ] \n\n");

    char fn[256];
    FILE *fp;
    size_t n;

    sprintf( fn, "%s.param%03d", xn->name, sNo );
    if( (fp = fopen( fn, "w")) != NULL ){
        fprintf(fp, "mu, lambda, pi");
        for( i = 0 ; i < sNo ; i++ )
        {   fprintf(fp, ", A%dx", i);   }
        fprintf(fp, "\n");

        for( i = 0 ; i < sNo ; i++ ){
            fprintf(fp, "%g, %g, %g", p->avgMu[i], p->avgLm[i], p->avgPi[i]);
            for( j = 0 ; j < sNo ; j++ )
            {   fprintf(fp, ", %g", p->avgA[j][i]);   }
            fprintf(fp, "\n");
        }
        fclose(fp);
    }

    sprintf( fn, "%s.Lq%03d", xn->name, sNo );
    if( (fp = fopen( fn, "w")) != NULL ){
        for( n = 0 ; n < gv->iteration ; n++ ){
            fprintf( fp, "%24.20e\n", gv->LqArr[n] );
        }
        fclose(fp);
    }

#ifdef OUTPUT_MAX_GAMMA
    sprintf( fn, "%s.maxG%03d", xn->name, sNo );
    if( (fp = fopen( fn, "w")) != NULL ){
        for( n = 0 ; n < xn->N ; n++ ){
            fprintf( fp, "%d\n", iv->gammaTraj[n] );
        }
        fclose(fp);
    }
#endif

    sprintf( fn, "%s.maxS%03d", xn->name, sNo );
    if( (fp = fopen( fn, "w")) != NULL ){
        for( n = 0 ; n < xn->N ; n++ ){
            fprintf( fp, "%d\n", iv->stateTraj[n] );
        }
        fclose(fp);
    }

}

void outputGaussResultsG( xns, gv, ivs, logFP )
xnDataBundle *xns;
globalVars *gv;
indVarBundle *ivs;
FILE *logFP;
{
    int sNo = gv->sNo, rNo = xns->R;
    gaussParameters *p = (gaussParameters*)gv->params;
    int i, j, r;

    fprintf(logFP, "  results: K = %d \n", sNo);

    fprintf(logFP, "   means: ( %g", p->avgMu[0]);
    for( i = 1 ; i < sNo ; i++ )
    {   fprintf(logFP, ", %g", p->avgMu[i]);   }
    fprintf(logFP, " ) \n");

    fprintf(logFP, "   lambda: ( %g", p->avgLm[0]);
    for( i = 1 ; i < sNo ; i++ ){
        fprintf(logFP, ", %g", p->avgLm[i]);
    }
    fprintf(logFP, " ) \n");

    fprintf(logFP, "   pi: ( %g", p->avgPi[0]);
    for( i = 1 ; i < sNo ; i++ ){
        fprintf(logFP, ", %g", p->avgPi[i]);
    }
    fprintf(logFP, " ) \n");

    fprintf(logFP, "   A_matrix: [");
    for( i = 0 ; i < sNo ; i++ ){
        fprintf(logFP, " ( %g", p->avgA[i][0]);
        for( j = 1 ; j < sNo ; j++ )
        {   fprintf(logFP, ", %g", p->avgA[i][j]);   }
        fprintf(logFP, ")");
    }
    fprintf(logFP, " ] \n\n");

    char fn[256];
    FILE *fp;
    size_t n;

    sprintf( fn, "%s.g.param%03d", xns->xn[0]->name, sNo );
    if( (fp = fopen( fn, "w")) != NULL ){
        fprintf(fp, "mu, lambda, pi");
        for( i = 0 ; i < sNo ; i++ )
        {   fprintf(fp, ", A%dx", i);   }
        fprintf(fp, "\n");

        for( i = 0 ; i < sNo ; i++ ){
            fprintf(fp, "%g, %g, %g", p->avgMu[i], p->avgLm[i], p->avgPi[i]);
            for( j = 0 ; j < sNo ; j++ )
            {   fprintf(fp, ", %g", p->avgA[j][i]);   }
            fprintf(fp, "\n");
        }
        fclose(fp);
    }

    sprintf( fn, "%s.g.Lq%03d", xns->xn[0]->name, sNo );
    if( (fp = fopen( fn, "w")) != NULL ){
        for( n = 0 ; n < gv->iteration ; n++ ){
            fprintf( fp, "%24.20e\n", gv->LqArr[n] );
        }
        fclose(fp);
    }

    int flag;
#ifdef OUTPUT_MAX_GAMMA
    sprintf( fn, "%s.g.maxG%03d", xns->xn[0]->name, sNo );
    flag = 0;
    if( (fp = fopen( fn, "w")) != NULL ){
        n = 0;
        do{
            flag = 1;
            for( r = 0 ; r < rNo ; r++ ){
                xnDataSet *xn = xns->xn[r];
                indVars *iv = ivs->indVars[r];

                if( r > 0 ){
                    fprintf( fp, "," );
                }
                if( n < xn->N ){
                    fprintf( fp, "%d", iv->gammaTraj[n] );
                }
                flag &= (n >= (xn->N - 1));
            }
            fprintf( fp, "\n" );
            n++;
        }while( !flag );
        fclose(fp);
    }
#endif

    sprintf( fn, "%s.g.maxS%03d", xns->xn[0]->name, sNo );
    flag = 0;
    if( (fp = fopen( fn, "w")) != NULL ){
        n = 0;
        do{
            flag = 1;
            for( r = 0 ; r < rNo ; r++ ){
                xnDataSet *xn = xns->xn[r];
                indVars *iv = ivs->indVars[r];

                if( r > 0 ){
                    fprintf( fp, "," );
                }
                if( n < xn->N ){
                    fprintf( fp, "%d", iv->stateTraj[n] );
                }
                flag &= (n >= (xn->N - 1));
            }
            fprintf( fp, "\n" );
            n++;
        }while( !flag );
        fclose(fp);
    }
}

//

max_stars_repo_name: okamoto-kenji/varBayes-HMM

max_issues_repo_name: okamoto-kenji/varBayes-HMM

max_forks_repo_name: okamoto-kenji/varBayes-HMM


text:
/**
 * File: trfit_utils.h
 *
 */
#ifndef _TRFIT_UTILS_H
#define _TRFIT_UTILS_H
#include "aper_conf.h"
#include <gsl/gsl_multifit_nlin.h>

#define N_COLUMNS_FIT 6
#define N_KAPPASIG_ITER 3
#define N_KAPPASIG_SIG  3.0
#define N_ITER_IPC  5
#define KAPPA_IPC  3.0

// the ipc fit to the
// data must be lower than this number
// to initialize a correction of the spectrum
#define MIN_IPC_LENGTH 60.0

typedef struct
{
  int n_data;

  double *x_values;
  double *y_values;
  double *e_values;
}
  fit_data;

struct function_data {
  int      n;
  double * x;
  double * y;
  double * sig;
};

extern gsl_vector *
fit_beamtrace(const aperture_conf  *conf, observation *obs,
	      const int beamID, beam act_beam);

double
gagauss(const double x, const double *params);

extern gsl_vector *
fit_wuttke(const fit_data *f_data);

extern gsl_vector *
fit_wuttke_talk(const fit_data *f_data);

extern d_point
find_grav_center(const fit_data *f_data);

extern px_point
get_fit_xrange(const aperture_conf  *conf, const observation *obs,
	       beam act_beam);

extern gsl_vector *
fit_function(const fit_data *f_data);

extern fit_data *
get_fitdata(observation *obs, int np, px_point tr_point);

extern fit_data *
alloc_fit_data(const int n_data);

void
print_fit_data(fit_data *f_data);

extern void
free_fit_data(fit_data *f_data);

int
gauss_f(const gsl_vector *x, void *params, gsl_vector *f);

int
gauss_df(const gsl_vector *x, void *params, gsl_matrix *J);

int
gauss_fdf (const gsl_vector *x, void *params, gsl_vector *f, gsl_matrix *J);

extern void
print_fit_state (size_t iter, gsl_multifit_fdfsolver * s);

extern double
comp_intrel_max(const fit_data *f_data);

extern gsl_matrix *
bcksub_observation(observation *obs, observation *bck);

extern double *
evaluate_ipc_fit(fit_data *f_data, double *gpars, double *fpars);

extern void
write_ipc_fit(fit_data *f_data, double *y_values, const char *ipc_file, const double qual, double *yfract);

extern double
fit_ipc_data(fit_data *fdata, double *gpars, int n_gpars,
	     double *fpars, int n_fpars);

extern fit_data *
get_ipc_fdata(const aperture_conf  *conf, observation *obs,
	      gsl_matrix *data_vals, beam act_beam);

extern double *
get_wf_special(const aperture_conf  *conf,  observation *obs,
	       gsl_matrix *data_vals, beam act_beam, double yoffs);

extern gsl_vector *
get_ipc_coldata(gsl_matrix *img_data, int np, px_point tr_point);

extern fit_data *
strip_fitdata(fit_data *old_fdata, int i_max);

extern double *
strip_wf_special(fit_data *old_fdata, double *y_fract, int i_max);

extern d_point
kappa_sigma_klipp_ipc(const aperture_conf  *conf, observation *obs,
		      gsl_matrix *data_matrix, beam act_beam,
		      const char *ipc_file);

extern void
reject_kappasigma(fit_data *f_data, double *y_values, double max_diff);

extern double
calculate_sigma(fit_data *f_data, double *y_values);
#endif

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe


text:
/*  Copyright 2017 International Business Machines Corporation

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.  */

#pragma once

#include <kernelpp/types.h>
#include <kernelpp/kernel.h>

#include <gsl.h>

namespace mwe {
namespace kernels
{
    using kernelpp::compute_mode;
    using kernelpp::error_code;

    /*  Declare a kernel 'add' which is overloaded to operate
     *  on single or array-like inputs.
     *
     *  Also, declare the compute modes (CPU,GPU) which
     *  this kernel will support.
     */
    KERNEL_DECL(add,
        compute_mode::CPU, compute_mode::CUDA)
    {
        template <compute_mode> static kernelpp::variant<int, error_code> op(
            int a, int b
            );

        template <compute_mode> static error_code op(
            const gsl::span<int> a, const gsl::span<int> b, gsl::span<int> result
            );
    };
}
}

max_stars_repo_name: rayglover-ibm/cuda-bindings

max_issues_repo_name: rayglover-ibm/cuda-bindings

max_forks_repo_name: rayglover-ibm/cuda-bindings


text:
/*
 * harvest.c -- compute solutions of harvest equation
 *
 * (c) 2016 Prof Dr Andreas Mueller,  Hochschule Rapperswil
 */
#include <stdlib.h>
#include <stdio.h>
#define __USE_XOPEN
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv2.h>

typedef struct param_s {
	double	a;
	double	h;
} param_t;

param_t	p0 = { .a = 5, .h = 0.8 };

int	f(double x, const double y[], double f[], void *params) {
	param_t	*p = (param_t *)params;
	f[0] = p->a * y[0] * (1 - y[0]) - p->h * (1 + sin(2 * M_PI * x));
	return GSL_SUCCESS;
}

int	singlecurve(FILE *out, const char *name, double x0, double y0,
		int direction) {
	printf("%s: %f\n", name, y0);
	gsl_odeiv2_system	system = { f, NULL, 1, &p0 };
	gsl_odeiv2_driver	*driver
		= gsl_odeiv2_driver_alloc_y_new(&system, gsl_odeiv2_step_rk8pd,
			1e-6 * direction, 1e-6, 0.0);
	fprintf(out, "path %s;\n", name);
	double	x = x0;
	double	y[1] = { y0 };
	fprintf(out, "%s = (%.3f,%.3f)", name, x, y[0]);
	for (int i = 1; i < 500; i++) {
		double	xnext = x + 0.01 * direction;
		int	status = gsl_odeiv2_driver_apply(driver, &x, xnext, y);
		if (status != GSL_SUCCESS) {
			fprintf(stderr, "error: rc = %d\n", status);
			fprintf(out, ";\n");
			return -1;
		}
		if (fabs(y[0]) > 100) {
			goto end;
		}
		fprintf(out, "\n    --(%.3f,%.3f)", x, y[0]);
	}
end:
	fprintf(out, ";\n");
	gsl_odeiv2_driver_free(driver);
	return 1;
}

int	centercurve(FILE *out, const char *name, double y0) {
	printf("%s: %f\n", name, y0);
	gsl_odeiv2_system	system = { f, NULL, 1, &p0 };
	gsl_odeiv2_driver	*driver
		= gsl_odeiv2_driver_alloc_y_new(&system, gsl_odeiv2_step_rk8pd,
			-1e-6, 1e-6, 0.0);
	fprintf(out, "path %s;\n", name);
	double	x = 2.5;
	double	y[1] = { y0 };
	// integrate backwards to the initial value
	gsl_odeiv2_driver_apply(driver, &x, 0, y);
	gsl_odeiv2_driver_free(driver);
	driver = gsl_odeiv2_driver_alloc_y_new(&system, gsl_odeiv2_step_rk8pd,
			1e-6, 1e-6, 0.0);
	if (fabs(y[0]) > 100) {
		goto end;
	}
	fprintf(out, "%s = (%.3f,%.3f)", name, x, y[0]);
	for (int i = 1; i < 501; i++) {
		double	xnext = x + 0.01;
		int	status = gsl_odeiv2_driver_apply(driver, &x, xnext, y);
		if (status != GSL_SUCCESS) {
			fprintf(stderr, "error: rc = %d\n", status);
			fprintf(out, ";\n");
			return -1;
		}
		if (fabs(y[0]) > 100) {
			goto end;
		}
		fprintf(out, "\n    --(%.3f,%.3f)", x, y[0]);
	}
end:
	fprintf(out, ";\n");
	gsl_odeiv2_driver_free(driver);
	return 1;
}

int	limitcurve(FILE *out, const char *name, int direction) {
	gsl_odeiv2_system	system = { f, NULL, 1, &p0 };
	gsl_odeiv2_driver	*driver
		= gsl_odeiv2_driver_alloc_y_new(&system, gsl_odeiv2_step_rk8pd,
			direction * 1e-6, 1e-6, 0.0);
	double	x = 0;
	double	y[1] = { 0.5 };
	double	xnext = 1000 * direction;
	int	status = gsl_odeiv2_driver_apply(driver, &x, xnext, y);
	printf("integration from %f to %f: %f (%d)\n", x, xnext, y[0], status);
	gsl_odeiv2_driver_free(driver);
	// now use this as the initial value
	return singlecurve(out, name, 2.5 * (1 - direction), y[0], direction);
}

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

	FILE	*out = fopen("harvest.mp", "w");
	char	suffix = 'a';
	char	curvename[10];
	for (int i = 0; i < 5; i++) {
		snprintf(curvename, sizeof(curvename), "path%c", suffix++);
		double	y0 = 0.7 + 0.2 * i;
		singlecurve(out, curvename, 0, y0, 1);
	}
	snprintf(curvename, sizeof(curvename), "path%c", suffix++);
	singlecurve(out, curvename, 0, 0.318, 1);
	snprintf(curvename, sizeof(curvename), "path%c", suffix++);
	singlecurve(out, curvename, 0, 0.315, 1);
	snprintf(curvename, sizeof(curvename), "path%c", suffix++);
	singlecurve(out, curvename, 0, 0.300, 1);
	for (int i = 0; i < 5; i++) {
		snprintf(curvename, sizeof(curvename), "path%c", suffix++);
		double	y0 = 0.3 - 0.2 * i;
		singlecurve(out, curvename, 5, y0, -1);
	}
	snprintf(curvename, sizeof(curvename), "path%c", suffix++);
	singlecurve(out, curvename, 5, 0.892, -1);
	snprintf(curvename, sizeof(curvename), "path%c", suffix++);
	singlecurve(out, curvename, 5, 0.895, -1);
	snprintf(curvename, sizeof(curvename), "path%c", suffix++);
	singlecurve(out, curvename, 5, 0.91, -1);
	for (int i = -3; i <= 3; i++) {
		snprintf(curvename, sizeof(curvename), "path%c", suffix++);
		double	y0 = 0.395 + 0.0945 * i;
		centercurve(out, curvename, y0);
	}

	snprintf(curvename, sizeof(curvename), "path%c", suffix++);
	limitcurve(out, curvename, -1);
	snprintf(curvename, sizeof(curvename), "path%c", suffix++);
	limitcurve(out, curvename, 1);

	fclose(out);

	return EXIT_SUCCESS;
}

max_stars_repo_name: MatthiasRubin/SeminarDGL

max_issues_repo_name: MatthiasRubin/SeminarDGL

max_forks_repo_name: MatthiasRubin/SeminarDGL