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

#include <config.h>
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>

int
main (void)
{
  double y, y_expected;
  int e, e_expected;

  gsl_ieee_env_setup ();

  /* Test for expm1 */

  y = gsl_expm1 (0.0);
  y_expected = 0.0;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_expm1(0.0)");

  y = gsl_expm1 (1e-10);
  y_expected = 1.000000000050000000002e-10;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_expm1(1e-10)");

  y = gsl_expm1 (-1e-10);
  y_expected = -9.999999999500000000017e-11;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_expm1(-1e-10)");

  y = gsl_expm1 (0.1);
  y_expected = 0.1051709180756476248117078264902;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_expm1(0.1)");

  y = gsl_expm1 (-0.1);
  y_expected = -0.09516258196404042683575094055356;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_expm1(-0.1)");

  y = gsl_expm1 (10.0);
  y_expected = 22025.465794806716516957900645284;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_expm1(10.0)");

  y = gsl_expm1 (-10.0);
  y_expected = -0.99995460007023751514846440848444;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_expm1(-10.0)");

  /* Test for log1p */

  y = gsl_log1p (0.0);
  y_expected = 0.0;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_log1p(0.0)");

  y = gsl_log1p (1e-10);
  y_expected = 9.9999999995000000000333333333308e-11;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_log1p(1e-10)");

  y = gsl_log1p (0.1);
  y_expected = 0.095310179804324860043952123280765;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_log1p(0.1)");

  y = gsl_log1p (10.0);
  y_expected = 2.3978952727983705440619435779651;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_log1p(10.0)");

  /* Test for gsl_hypot */

  y = gsl_hypot (0.0, 0.0);
  y_expected = 0.0;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_hypot(0.0, 0.0)");

  y = gsl_hypot (1e-10, 1e-10);
  y_expected = 1.414213562373095048801688e-10;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_hypot(1e-10, 1e-10)");

  y = gsl_hypot (1e-38, 1e-38);
  y_expected = 1.414213562373095048801688e-38;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_hypot(1e-38, 1e-38)");

  y = gsl_hypot (1e-10, -1.0);
  y_expected = 1.000000000000000000005;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_hypot(1e-10, -1)");

  y = gsl_hypot (-1.0, 1e-10);
  y_expected = 1.000000000000000000005;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_hypot(-1, 1e-10)");

  y = gsl_hypot (1e307, 1e301);
  y_expected = 1.000000000000499999999999e307;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_hypot(1e307, 1e301)");

  y = gsl_hypot (1e301, 1e307);
  y_expected = 1.000000000000499999999999e307;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_hypot(1e301, 1e307)");

  y = gsl_hypot (1e307, 1e307);
  y_expected = 1.414213562373095048801688e307;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_hypot(1e307, 1e307)");

  /* Test +-Inf, finite */

  y = gsl_hypot (GSL_POSINF, 1.2);
  y_expected = GSL_POSINF;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_hypot(GSL_POSINF, 1.2)");

  y = gsl_hypot (GSL_NEGINF, 1.2);
  y_expected = GSL_POSINF;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_hypot(GSL_NEGINF, 1.2)");

  y = gsl_hypot (1.2, GSL_POSINF);
  y_expected = GSL_POSINF;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_hypot(1.2, GSL_POSINF)");

  y = gsl_hypot (1.2, GSL_NEGINF);
  y_expected = GSL_POSINF;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_hypot(1.2, GSL_NEGINF)");

  /* Test NaN, finite */

  y = gsl_hypot (GSL_NAN, 1.2);
  y_expected = GSL_NAN;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_hypot(GSL_NAN, 1.2)");

  y = gsl_hypot (1.2, GSL_NAN);
  y_expected = GSL_NAN;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_hypot(1.2, GSL_NAN)");

  /* Test NaN, NaN */

  y = gsl_hypot (GSL_NAN, GSL_NAN);
  y_expected = GSL_NAN;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_hypot(GSL_NAN, GSL_NAN)");

  /* Test +Inf, NaN */

  y = gsl_hypot (GSL_POSINF, GSL_NAN);
  y_expected = GSL_POSINF;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_hypot(GSL_POSINF, GSL_NAN)");

  /* Test -Inf, NaN */

  y = gsl_hypot (GSL_NEGINF, GSL_NAN);
  y_expected = GSL_POSINF;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_hypot(GSL_NEGINF, GSL_NAN)");

  /* Test NaN, +Inf */

  y = gsl_hypot (GSL_NAN, GSL_POSINF);
  y_expected = GSL_POSINF;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_hypot(GSL_NAN, GSL_POSINF)");

  /* Test NaN, -Inf */

  y = gsl_hypot (GSL_NAN, GSL_NEGINF);
  y_expected = GSL_POSINF;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_hypot(GSL_NAN, GSL_NEGINF)");

  /* Test for gsl_hypot3 */

  y = gsl_hypot3 (0.0, 0.0, 0.0);
  y_expected = 0.0;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_hypot3(0.0, 0.0, 0.0)");

  y = gsl_hypot3 (1e-10, 1e-10, 1e-10);
  y_expected = 1.732050807568877293527446e-10;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_hypot3(1e-10, 1e-10, 1e-10)");

  y = gsl_hypot3 (1e-38, 1e-38, 1e-38);
  y_expected = 1.732050807568877293527446e-38;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_hypot3(1e-38, 1e-38, 1e-38)");

  y = gsl_hypot3 (1e-10, 1e-10, -1.0);
  y_expected = 1.000000000000000000099;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_hypot3(1e-10, 1e-10, -1)");

  y = gsl_hypot3 (1e-10, -1.0, 1e-10);
  y_expected = 1.000000000000000000099;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_hypot3(1e-10, -1, 1e-10)");

  y = gsl_hypot3 (-1.0, 1e-10, 1e-10);
  y_expected = 1.000000000000000000099;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_hypot3(-1, 1e-10, 1e-10)");

  y = gsl_hypot3 (1e307, 1e301, 1e301);
  y_expected = 1.0000000000009999999999995e307;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_hypot3(1e307, 1e301, 1e301)");

  y = gsl_hypot3 (1e307, 1e307, 1e307);
  y_expected = 1.732050807568877293527446e307;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_hypot3(1e307, 1e307, 1e307)");

  y = gsl_hypot3 (1e307, 1e-307, 1e-307);
  y_expected = 1.0e307;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_hypot3(1e307, 1e-307, 1e-307)");

  /* Test for acosh */

  y = gsl_acosh (1.0);
  y_expected = 0.0;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_acosh(1.0)");

  y = gsl_acosh (1.1);
  y_expected = 4.435682543851151891329110663525e-1;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_acosh(1.1)");

  y = gsl_acosh (10.0);
  y_expected = 2.9932228461263808979126677137742e0;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_acosh(10.0)");

  y = gsl_acosh (1e10);
  y_expected = 2.3718998110500402149594646668302e1;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_acosh(1e10)");

  /* Test for asinh */

  y = gsl_asinh (0.0);
  y_expected = 0.0;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_asinh(0.0)");

  y = gsl_asinh (1e-10);
  y_expected = 9.9999999999999999999833333333346e-11;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_asinh(1e-10)");

  y = gsl_asinh (-1e-10);
  y_expected = -9.9999999999999999999833333333346e-11;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_asinh(1e-10)");

  y = gsl_asinh (0.1);
  y_expected = 9.983407889920756332730312470477e-2;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_asinh(0.1)");

  y = gsl_asinh (-0.1);
  y_expected = -9.983407889920756332730312470477e-2;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_asinh(-0.1)");

  y = gsl_asinh (1.0);
  y_expected = 8.8137358701954302523260932497979e-1;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_asinh(1.0)");

  y = gsl_asinh (-1.0);
  y_expected = -8.8137358701954302523260932497979e-1;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_asinh(-1.0)");

  y = gsl_asinh (10.0);
  y_expected = 2.9982229502979697388465955375965e0;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_asinh(10)");

  y = gsl_asinh (-10.0);
  y_expected = -2.9982229502979697388465955375965e0;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_asinh(-10)");

  y = gsl_asinh (1e10);
  y_expected = 2.3718998110500402149599646668302e1;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_asinh(1e10)");

  y = gsl_asinh (-1e10);
  y_expected = -2.3718998110500402149599646668302e1;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_asinh(-1e10)");

  /* Test for atanh */

  y = gsl_atanh (0.0);
  y_expected = 0.0;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_atanh(0.0)");

  y = gsl_atanh (1e-20);
  y_expected = 1e-20;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_atanh(1e-20)");

  y = gsl_atanh (-1e-20);
  y_expected = -1e-20;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_atanh(-1e-20)");

  y = gsl_atanh (0.1);
  y_expected = 1.0033534773107558063572655206004e-1;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_atanh(0.1)");

  y = gsl_atanh (-0.1);
  y_expected = -1.0033534773107558063572655206004e-1;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_atanh(-0.1)");

  y = gsl_atanh (0.9);
  y_expected = 1.4722194895832202300045137159439e0;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_atanh(0.9)");

  y = gsl_atanh (-0.9);
  y_expected = -1.4722194895832202300045137159439e0;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_atanh(0.9)");

  /* Test for pow_int */

  y = gsl_pow_2 (-3.14);
  y_expected = pow (-3.14, 2.0);
  gsl_test_rel (y, y_expected, 1e-15, "gsl_pow_2(-3.14)");

  y = gsl_pow_3 (-3.14);
  y_expected = pow (-3.14, 3.0);
  gsl_test_rel (y, y_expected, 1e-15, "gsl_pow_3(-3.14)");

  y = gsl_pow_4 (-3.14);
  y_expected = pow (-3.14, 4.0);
  gsl_test_rel (y, y_expected, 1e-15, "gsl_pow_4(-3.14)");

  y = gsl_pow_5 (-3.14);
  y_expected = pow (-3.14, 5.0);
  gsl_test_rel (y, y_expected, 1e-15, "gsl_pow_5(-3.14)");

  y = gsl_pow_6 (-3.14);
  y_expected = pow (-3.14, 6.0);
  gsl_test_rel (y, y_expected, 1e-15, "gsl_pow_6(-3.14)");

  y = gsl_pow_7 (-3.14);
  y_expected = pow (-3.14, 7.0);
  gsl_test_rel (y, y_expected, 1e-15, "gsl_pow_7(-3.14)");

  y = gsl_pow_8 (-3.14);
  y_expected = pow (-3.14, 8.0);
  gsl_test_rel (y, y_expected, 1e-15, "gsl_pow_8(-3.14)");

  y = gsl_pow_9 (-3.14);
  y_expected = pow (-3.14, 9.0);
  gsl_test_rel (y, y_expected, 1e-15, "gsl_pow_9(-3.14)");

  {
    int n;
    for (n = -9; n < 10; n++)
      {
        y = gsl_pow_int (-3.14, n);
        y_expected = pow (-3.14, n);
        gsl_test_rel (y, y_expected, 1e-15, "gsl_pow_int(-3.14,%d)", n);
      }
  }

  {
    unsigned int n;
    for (n = 0; n < 10; n++)
      {
        y = gsl_pow_uint (-3.14, n);
        y_expected = pow (-3.14, n);
        gsl_test_rel (y, y_expected, 1e-15, "gsl_pow_uint(-3.14,%d)", n);
      }
  }

  /* Test case for n at INT_MAX, INT_MIN */

  {
    double u = 1.0000001;
    int n = INT_MAX;
    y = gsl_pow_int (u, n);
    y_expected = pow (u, n);
    gsl_test_rel (y, y_expected, 1e-6, "gsl_pow_int(%.7f,%d)", u, n);

    n = INT_MIN;
    y = gsl_pow_int (u, n);
    y_expected = pow (u, n);
    gsl_test_rel (y, y_expected, 1e-6, "gsl_pow_int(%.7f,%d)", u, n);
  }

  /* Test for ldexp */

  y = gsl_ldexp (M_PI, -2);
  y_expected = M_PI_4;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_ldexp(pi,-2)");

  y = gsl_ldexp (1.0, 2);
  y_expected = 4.000000;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_ldexp(1.0,2)");

  y = gsl_ldexp (0.0, 2);
  y_expected = 0.0;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_ldexp(0.0,2)");

  y = gsl_ldexp (9.999999999999998890e-01, 1024);
  y_expected = GSL_DBL_MAX;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_ldexp DBL_MAX");

  y = gsl_ldexp (1e308, -2000);
  y_expected = 8.7098098162172166755761e-295;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_ldexp(1e308,-2000)");

  y = gsl_ldexp (GSL_DBL_MIN, 2000);
  y_expected = 2.554675596204441378334779940e294;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_ldexp(DBL_MIN,2000)");

  /* Test subnormals */

  {
    int i = 0;
    volatile double x = GSL_DBL_MIN;
    y_expected = 2.554675596204441378334779940e294;

    x /= 2;
    while (x > 0)
      {
        i++ ;
        y = gsl_ldexp (x, 2000 + i);
        gsl_test_rel (y, y_expected, 1e-15, "gsl_ldexp(DBL_MIN/2**%d,%d)",i,2000+i);
        x /= 2;
      }
  }

  /* Test for frexp */

  y = gsl_frexp (0.0, &e);
  y_expected = 0;
  e_expected = 0;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_frexp(0) fraction");
  gsl_test_int (e, e_expected, "gsl_frexp(0) exponent");

  y = gsl_frexp (M_PI, &e);
  y_expected = M_PI_4;
  e_expected = 2;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_frexp(pi) fraction");
  gsl_test_int (e, e_expected, "gsl_frexp(pi) exponent");

  y = gsl_frexp (2.0, &e);
  y_expected = 0.5;
  e_expected = 2;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_frexp(2.0) fraction");
  gsl_test_int (e, e_expected, "gsl_frexp(2.0) exponent");

  y = gsl_frexp (1.0 / 4.0, &e);
  y_expected = 0.5;
  e_expected = -1;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_frexp(0.25) fraction");
  gsl_test_int (e, e_expected, "gsl_frexp(0.25) exponent");

  y = gsl_frexp (1.0 / 4.0 - 4.0 * GSL_DBL_EPSILON, &e);
  y_expected = 0.999999999999996447;
  e_expected = -2;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_frexp(0.25-eps) fraction");
  gsl_test_int (e, e_expected, "gsl_frexp(0.25-eps) exponent");

  y = gsl_frexp (GSL_DBL_MAX, &e);
  y_expected = 9.999999999999998890e-01;
  e_expected = 1024;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_frexp(DBL_MAX) fraction");
  gsl_test_int (e, e_expected, "gsl_frexp(DBL_MAX) exponent");

  y = gsl_frexp (-GSL_DBL_MAX, &e);
  y_expected = -9.999999999999998890e-01;
  e_expected = 1024;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_frexp(-DBL_MAX) fraction");
  gsl_test_int (e, e_expected, "gsl_frexp(-DBL_MAX) exponent");

  y = gsl_frexp (GSL_DBL_MIN, &e);
  y_expected = 0.5;
  e_expected = -1021;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_frexp(DBL_MIN) fraction");
  gsl_test_int (e, e_expected, "gsl_frexp(DBL_MIN) exponent");

  y = gsl_frexp (-GSL_DBL_MIN, &e);
  y_expected = -0.5;
  e_expected = -1021;
  gsl_test_rel (y, y_expected, 1e-15, "gsl_frexp(-DBL_MIN) fraction");
  gsl_test_int (e, e_expected, "gsl_frexp(-DBL_MIN) exponent");

  /* Test subnormals */

  {
    int i = 0;
    volatile double x = GSL_DBL_MIN;
    y_expected = 0.5;
    e_expected = -1021;

    x /= 2;

    while (x > 0)
      {
        e_expected--;
        i++ ;

        y = gsl_frexp (x, &e);
        gsl_test_rel (y, y_expected, 1e-15, "gsl_frexp(DBL_MIN/2**%d) fraction",i);
        gsl_test_int (e, e_expected, "gsl_frexp(DBL_MIN/2**%d) exponent", i);
        x /= 2;
      }
  }

  /* Test for approximate floating point comparison */
  {
    double x, y;
    int i;

    x = M_PI;
    y = 22.0 / 7.0;

    /* test the basic function */

    for (i = 0; i < 10; i++)
      {
        double tol = pow (10, -i);
        int res = gsl_fcmp (x, y, tol);
        gsl_test_int (res, -(i >= 4), "gsl_fcmp(%.5f,%.5f,%g)", x, y, tol);
      }

    for (i = 0; i < 10; i++)
      {
        double tol = pow (10, -i);
        int res = gsl_fcmp (y, x, tol);
        gsl_test_int (res, (i >= 4), "gsl_fcmp(%.5f,%.5f,%g)", y, x, tol);
      }
  }

#if HAVE_IEEE_COMPARISONS
  /* Test for isinf, isnan, finite */

  {
    double zero, one, inf, nan;
    int s;

    zero = 0.0;
    one = 1.0;
    inf = exp (1.0e10);
    nan = inf / inf;

    s = gsl_isinf (zero);
    gsl_test_int (s, 0, "gsl_isinf(0)");

    s = gsl_isinf (one);
    gsl_test_int (s, 0, "gsl_isinf(1)");

    s = gsl_isinf (inf);
    gsl_test_int (s, 1, "gsl_isinf(inf)");

    s = gsl_isinf (-inf);
    gsl_test_int (s, -1, "gsl_isinf(-inf)");

    s = gsl_isinf (nan);
    gsl_test_int (s, 0, "gsl_isinf(nan)");

    s = gsl_isnan (zero);
    gsl_test_int (s, 0, "gsl_isnan(0)");

    s = gsl_isnan (one);
    gsl_test_int (s, 0, "gsl_isnan(1)");

    s = gsl_isnan (inf);
    gsl_test_int (s, 0, "gsl_isnan(inf)");

    s = gsl_isnan (-inf);
    gsl_test_int (s, 0, "gsl_isnan(-inf)");

    s = gsl_isnan (nan);
    gsl_test_int (s, 1, "gsl_isnan(nan)");

    s = gsl_finite (zero);
    gsl_test_int (s, 1, "gsl_finite(0)");

    s = gsl_finite (one);
    gsl_test_int (s, 1, "gsl_finite(1)");

    s = gsl_finite (inf);
    gsl_test_int (s, 0, "gsl_finite(inf)");

    s = gsl_finite (-inf);
    gsl_test_int (s, 0, "gsl_finite(-inf)");

    s = gsl_finite (nan);
    gsl_test_int (s, 0, "gsl_finite(nan)");
  }
#endif

  {
    double x = gsl_fdiv (2.0, 3.0);
    gsl_test_rel (x, 2.0 / 3.0, 4 * GSL_DBL_EPSILON, "gsl_fdiv(2,3)");
  }

  /* Test constants in gsl_math.h */

  {
    double x = log(M_E);
    gsl_test_rel (x, 1.0, 4 * GSL_DBL_EPSILON, "ln(M_E)");
  }

  {
    double x=pow(2.0,M_LOG2E);
    gsl_test_rel (x, exp(1.0), 4 * GSL_DBL_EPSILON, "2^M_LOG2E");
  }

  {
    double x=pow(10.0,M_LOG10E);
    gsl_test_rel (x, exp(1.0), 4 * GSL_DBL_EPSILON, "10^M_LOG10E");
  }

  {
    double x=pow(M_SQRT2, 2.0);
    gsl_test_rel (x, 2.0, 4 * GSL_DBL_EPSILON, "M_SQRT2^2");
  }

  {
    double x=pow(M_SQRT1_2, 2.0);
    gsl_test_rel (x, 1.0/2.0, 4 * GSL_DBL_EPSILON, "M_SQRT1_2");
  }

  {
    double x=pow(M_SQRT3, 2.0);
    gsl_test_rel (x, 3.0, 4 * GSL_DBL_EPSILON, "M_SQRT3^2");
  }

  {
    double x = M_PI;
    gsl_test_rel (x, 3.1415926535897932384626433832795, 4 * GSL_DBL_EPSILON, "M_PI");
  }

  {
    double x = 2 * M_PI_2;
    gsl_test_rel (x, M_PI, 4 * GSL_DBL_EPSILON, "2*M_PI_2");
  }

  {
    double x = 4 * M_PI_4;
    gsl_test_rel (x, M_PI, 4 * GSL_DBL_EPSILON, "4*M_PI_4");
  }

  {
    double x = pow(M_SQRTPI, 2.0);
    gsl_test_rel (x, M_PI, 4 * GSL_DBL_EPSILON, "M_SQRTPI^2");
  }

  {
    double x = pow(M_2_SQRTPI, 2.0);
    gsl_test_rel (x, 4/M_PI, 4 * GSL_DBL_EPSILON, "M_SQRTPI^2");
  }

  {
    double x = M_1_PI;
    gsl_test_rel (x, 1/M_PI, 4 * GSL_DBL_EPSILON, "M_1_SQRTPI");
  }

  {
    double x = M_2_PI;
    gsl_test_rel (x, 2.0/M_PI, 4 * GSL_DBL_EPSILON, "M_2_PI");
  }

  {
    double x = exp(M_LN10);
    gsl_test_rel (x, 10, 4 * GSL_DBL_EPSILON, "exp(M_LN10)");
  }

  {
    double x = exp(M_LN2);
    gsl_test_rel (x, 2, 4 * GSL_DBL_EPSILON, "exp(M_LN2)");
  }

  {
    double x = exp(M_LNPI);
    gsl_test_rel (x, M_PI, 4 * GSL_DBL_EPSILON, "exp(M_LNPI)");
  }

  {
    double x = M_EULER;
    gsl_test_rel (x, 0.5772156649015328606065120900824, 4 * GSL_DBL_EPSILON, "M_EULER");
  }

  exit (gsl_test_summary ());
}

max_stars_repo_name: ruslankuzmin/julia

max_issues_repo_name: ruslankuzmin/julia

max_forks_repo_name: ruslankuzmin/julia


text:
/* coupled_oscillators.c Kuramoto model of coupled phase oscillators */
/*
 * Dynamique Populations Cellulaires
 * Université Lyon 1, 2014
 *
 * Fichiers:
 *  coupled_oscillators.c
 *  makefile
 *
 * Ce code utilise la librairie gsl (GNU Scientific Library)
 *
 * Pour compiler sur mac ou linux:
 *
 * ouvrir un terminal
 * >> make
 * >> ./coupled_oscillators
 *
 * Deux fichiers sont générés en sortie:
 *  osc.txt : series temporelles
 *  order.txt: paramètre d'ordre
 *
 * Le fichier osc.txt contient est structuré en colonnes::
 * 0   w1  w2 ...      wNOSC
 * t1  1   2   3   ...  NOSC
 *
 * La premiere ligne est composee des frequences des oscillateurs. Les
 * oscillateurs sont places en ordre croissant de frequence.
 * t est le temps, r le module du paramètre d'ordre, psi l'argument du paramètre
 * d'ordre, 1...NOSC les phases des oscillateur.
 * Les phases des oscillateurs sont dans R (pas dans [0, 2*pi]).
 *
 * Pour visualiser les résultats, on peut utiliser gnuplot ou Matlab
 *
 * GNUPLOT:
 * plot 'osc.txt' u 1:2:(sin($3)) nonuniform matrix with image
 *
 * Pour tracer le paramètre d'ordre
 *
 * pour r:
 * plot 'order.txt' using 1:2 w lines
 *
 * pour psi:
 * plot 'order.txt' using 1:3 w lines
 *
 * Pour Matlab, utiliser le script coupled_oscillators.m
 *
 */

/* ================================================================= */
/*                              Libraries                            */
/* ================================================================= */
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_odeiv2.h>

struct par {
    uint32_t NOSC;
    double K;
    double *w;
};

int ode_rhs(double t, const double y[], double f[], void *params);
int compare(const void *, const void *);

int main ( int argc, char *argv[] )
{
    uint32_t NOSC=50; /* nbr oscillators */
    double K=.5;  /* coupling stength */
    double sigma=0.5; /* std w distribution */
    double *y;     /* phase of oscillators */
    struct par mu; /* parameters */
    FILE *file_order,*file_out;
    uint32_t i;
    double t,tfinal;
    double h,hmin;
    double w_mean;
    double r,psi, rx, ry;

	/* get command line input parameters */
	if (argc == 4)
    {
        NOSC = strtoul(argv[1],NULL,10);
    	K = strtod(argv[2],NULL);
        tfinal = strtod(argv[3],NULL);
    }
	else
    {
    	fprintf(stderr, "Error: call to coupled_oscillators must have the following three mandatory arguments: unsigned long NOSC, double K, and double tfinal\n");
    	exit(EXIT_FAILURE);
    }

	/* init random number generator */
	const gsl_rng_type * rngT;
	gsl_rng * rand_nbr_gen;
	gsl_rng_env_setup();
	rngT = gsl_rng_default;
	rand_nbr_gen = gsl_rng_alloc(rngT);

    /* gsl ode */
    const gsl_odeiv2_step_type * odeT;
    gsl_odeiv2_step * s;
    gsl_odeiv2_control * c;
    gsl_odeiv2_evolve * e;
    gsl_odeiv2_system sys;
    int status;

    mu.NOSC = NOSC;
    mu.K = K;
    mu.w = malloc(NOSC*sizeof(double));
    y = malloc(NOSC*sizeof(double));
    w_mean = 0.0;
    for (i = 0; i < NOSC; i++)
    {
        y[i] = 2*M_PI*gsl_rng_uniform(rand_nbr_gen)-M_PI;
        mu.w[i] = gsl_ran_gaussian(rand_nbr_gen,sigma*sigma);
        w_mean += mu.w[i];
    }
    w_mean /= NOSC;

    /* sort the frequencies from slowest to fastest */
    qsort(mu.w, NOSC, sizeof(double), compare);

    /* remove the mean to mu.w */
    file_out = fopen("osc.txt","w");
    fprintf(file_out,"%.5e ",0.0);
    for ( i = 0; i < NOSC; i++)
    {
        mu.w[i] -= w_mean;
        fprintf(file_out,"%.5e ",mu.w[i]);
    }
    fprintf(file_out,"\n");
    /* DO NOT TOUCH mu.w BELOW */

    /* initialise the ode solver */
    /* the solver used is a Runge-Kutta-Fehlberg (4, 5) method */
    odeT = gsl_odeiv2_step_rkf45;
    s = gsl_odeiv2_step_alloc(odeT,NOSC);
    c = gsl_odeiv2_control_y_new(1e-6,0.0);
    e = gsl_odeiv2_evolve_alloc(NOSC);

    h = 1e-2;
    hmin = 1e-5;
    t = 0.0;

    file_order = fopen("order.txt","w");
    while (t < tfinal )
    {

        /* apply one step of the ode solver */
        sys = (gsl_odeiv2_system) {ode_rhs, NULL, NOSC, &mu};
        status = gsl_odeiv2_evolve_apply(e,c,s,&sys,&t,tfinal,&h,y);
        if (status != GSL_SUCCESS)
            break;

        rx = 0;
        ry = 0;
        for ( i = 0; i < NOSC; i++)
        {
            rx += cos(y[i]);
            ry += sin(y[i]);

        }
        rx /= NOSC;
        ry /= NOSC;
        r = hypot(rx,ry);
        psi = atan2(ry,rx);

        if (h < hmin) h = hmin;
        printf(".\n");

        fprintf(file_order,"%.5e %.5e %.5e\n",t,r,psi);
        fprintf(file_out,"%.5e ",t);

        for (i = 0; i < NOSC; i++)
        {
            fprintf (file_out,"%.5e ",y[i]);
        }
        fprintf(file_out,"\n");

    }

    gsl_odeiv2_evolve_free(e);
    gsl_odeiv2_control_free(c);
    gsl_odeiv2_step_free(s);

    fclose(file_out);
    fclose(file_order);

    free(y);
    free(mu.w);

	return 0;
}

int ode_rhs(double t, const double y[], double f[], void *params)
{

    struct par mu = *(struct par*)params;
    uint32_t NOSC = mu.NOSC;
    double K  = mu.K;
    double *w = mu.w;
    uint32_t i,j;
    double *c;

    /* printf("call to ode_rhs, NOSC=%d\n",NOSC); */
    c = malloc(NOSC*sizeof(double));
    /* compute coupling */
    for ( i = 0; i < NOSC; i++)
    {
        c[i] = 0.0;
        for ( j = 0; j < NOSC; j++)
        {
            c[i] += sin(y[j]-y[i]);
        }

        f[i] = w[i] + K/NOSC*c[i];
        /* printf("f[%d]=%f\n",i,f[i]); */

    }

    free(c);

    return GSL_SUCCESS;

}

int compare(const void *a, const void *b)
{
    const double *da = (const double *) a;
    const double *db = (const double *) b;

    return (*da > *db) - (*da < *db);
}

ext: c

max_stars_repo_name: samubernard/popdyn

max_issues_repo_name: samubernard/popdyn

max_forks_repo_name: samubernard/popdyn


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

#ifndef __GSL_FFT_COMPLEX_FLOAT_H__
#define __GSL_FFT_COMPLEX_FLOAT_H__

#include <stddef.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_fft.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

/*  Power of 2 routines  */

int gsl_fft_complex_float_radix2_forward (gsl_complex_packed_array_float data,
                                          const size_t stride,
                                          const size_t n);

int gsl_fft_complex_float_radix2_backward (gsl_complex_packed_array_float data,
                                           const size_t stride,
                                           const size_t n);

int gsl_fft_complex_float_radix2_inverse (gsl_complex_packed_array_float data,
                                          const size_t stride,
                                          const size_t n);

int gsl_fft_complex_float_radix2_transform (gsl_complex_packed_array_float data,
                                            const size_t stride,
                                            const size_t n,
                                            const gsl_fft_direction sign);

int gsl_fft_complex_float_radix2_dif_forward (gsl_complex_packed_array_float data,
                                              const size_t stride,
                                              const size_t n);

int gsl_fft_complex_float_radix2_dif_backward (gsl_complex_packed_array_float data,
                                               const size_t stride,
                                               const size_t n);

int gsl_fft_complex_float_radix2_dif_inverse (gsl_complex_packed_array_float data,
                                              const size_t stride,
                                              const size_t n);

int gsl_fft_complex_float_radix2_dif_transform (gsl_complex_packed_array_float data,
                                                const size_t stride,
                                                const size_t n,
                                                const gsl_fft_direction sign);

/*  Mixed Radix general-N routines  */

typedef struct
  {
    size_t n;
    size_t nf;
    size_t factor[64];
    gsl_complex_float *twiddle[64];
    gsl_complex_float *trig;
  }
gsl_fft_complex_wavetable_float;

typedef struct
{
    size_t n;
    float *scratch;
}
gsl_fft_complex_workspace_float;

gsl_fft_complex_wavetable_float *gsl_fft_complex_wavetable_float_alloc (size_t n);

void gsl_fft_complex_wavetable_float_free (gsl_fft_complex_wavetable_float * wavetable);

gsl_fft_complex_workspace_float *gsl_fft_complex_workspace_float_alloc (size_t n);

void gsl_fft_complex_workspace_float_free (gsl_fft_complex_workspace_float * workspace);

int gsl_fft_complex_float_memcpy (gsl_fft_complex_wavetable_float * dest,
                               gsl_fft_complex_wavetable_float * src);

int gsl_fft_complex_float_forward (gsl_complex_packed_array_float data,
                                   const size_t stride,
                                   const size_t n,
                                   const gsl_fft_complex_wavetable_float * wavetable,
                                   gsl_fft_complex_workspace_float * work);

int gsl_fft_complex_float_backward (gsl_complex_packed_array_float data,
                                    const size_t stride,
                                    const size_t n,
                                    const gsl_fft_complex_wavetable_float * wavetable,
                                    gsl_fft_complex_workspace_float * work);

int gsl_fft_complex_float_inverse (gsl_complex_packed_array_float data,
                                   const size_t stride,
                                   const size_t n,
                                   const gsl_fft_complex_wavetable_float * wavetable,
                                   gsl_fft_complex_workspace_float * work);

int gsl_fft_complex_float_transform (gsl_complex_packed_array_float data,
                                     const size_t stride, const size_t n,
                                     const gsl_fft_complex_wavetable_float * wavetable,
                                     gsl_fft_complex_workspace_float * work,
                                     const gsl_fft_direction sign);

__END_DECLS

#endif /* __GSL_FFT_COMPLEX_FLOAT_H__ */

ext: h

max_stars_repo_name: shi-bash-cmd/qtTest

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#pragma once
#include <cstddef>
#include <gsl/gsl.h>

namespace util{

template<typename T>
using span_dyn = gsl::span<T>;
template<typename T, int64_t M>
using span_1d = gsl::span<T,M>;
template<typename T, int64_t M, int64_t N>
using span_2d = gsl::span<T,M,N>;
template<typename T, int64_t L, int64_t M, int64_t N>
using span_3d = gsl::span<T,L,M,N>;

//Not works and param.cpp directly uses gsl.h
template<typename... Args>
auto as_span(Args&&... args){
    return gsl::as_span(std::forward<Args>(args)...);
}

template<int64_t T>
using dim = gsl::dim<T>;

template<std::ptrdiff_t Extent = gsl::dynamic_range>
using cstring_span = gsl::basic_string_span<const char, Extent>;

}//namespace util

max_stars_repo_name: uphere-co/nlp-prototype

max_issues_repo_name: uphere-co/nlp-prototype

max_forks_repo_name: uphere-co/nlp-prototype


text:
#pragma once

#include "Cesium3DTiles/Library.h"
#include "CesiumGltf/Model.h"
#include <functional>
#include <glm/mat4x4.hpp>
#include <gsl/span>
#include <optional>

namespace Cesium3DTiles {

/**
 * @brief Functions for loading and processing glTF models.
 *
 * This class offers basic functions for loading glTF models as
 * `CesiumGltf::Model` instances, and for processing the mesh primitives that
 * appear in the resulting models.
 */
class CESIUM3DTILES_API Gltf final {

public:
  /** @brief This class cannot be instantiated */
  Gltf() = delete;

  /**
   * @brief A callback function for {@link Gltf::forEachPrimitiveInScene}.
   */
  typedef void ForEachPrimitiveInSceneCallback(
      CesiumGltf::Model& gltf,
      CesiumGltf::Node& node,
      CesiumGltf::Mesh& mesh,
      CesiumGltf::MeshPrimitive& primitive,
      const glm::dmat4& transform);

  /**
   * @brief Apply the given callback to all relevant primitives.
   *
   * If the given `sceneID` is non-negative and exists in the given glTF,
   * then the given callback will be applied to all meshes of this scene.
   *
   * If the given `sceneId` is negative, then the meshes that the callback
   * will be applied to depends on the structure of the glTF model:
   *
   * * If the glTF model has a default scene, then it will
   *   be applied to all meshes of the default scene.
   * * Otherwise, it will be applied to all meshes of the the first scene.
   * * Otherwise (if the glTF model does not contain any scenes), it will
   *   be applied to all meshes that can be found by starting a traversal
   *   at the root node.
   * * Otherwise (if there are no scenes and no nodes), then all meshes
   *   will be traversed.
   *
   * @param gltf The glTF model.
   * @param sceneID The scene ID (index)
   * @param callback The callback to apply
   */
  static void forEachPrimitiveInScene(
      CesiumGltf::Model& gltf,
      int sceneID,
      std::function<ForEachPrimitiveInSceneCallback>&& callback);

  /**
   * @brief A callback function for {@link Gltf::forEachPrimitiveInScene}.
   */
  typedef void ForEachPrimitiveInSceneConstCallback(
      const CesiumGltf::Model& gltf,
      const CesiumGltf::Node& node,
      const CesiumGltf::Mesh& mesh,
      const CesiumGltf::MeshPrimitive& primitive,
      const glm::dmat4& transform);

  /** @copydoc Gltf::forEachPrimitiveInScene() */
  static void forEachPrimitiveInScene(
      const CesiumGltf::Model& gltf,
      int sceneID,
      std::function<ForEachPrimitiveInSceneConstCallback>&& callback);
};

} // namespace Cesium3DTiles


text:
#ifndef PDCD_SMSVM_H
#define PDCD_SMSVM_H

#include <string>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <stdio.h>      /* printf */
#include <time.h>
#include <fstream>
#include <algorithm>
#include <iomanip>
#include <ctime>
#include <math.h>

//This class implements the method SMART_CD.

/*
The optimization problem to solve is of the form f(x)+g(x)+h(Mx) ;
 where f(x)=sum_{j=1}^m lambda_f[j] phi_j(<A_j,x>) and g(x)=sum_{i=1}^n g_i(x_i), h(x)= |x|_1. We all assume that each phi_j is 1-smooth.
*/

template<typename L, typename D>
class PDCD_SMSVM
{
private:

  std::vector<D> Au;
  std::vector<D> Az;

  std::vector<D> Mu;
  std::vector<D> Mz;

  std::vector<D> Mx;
  std::vector<D> Ax;
  std::vector<D> lambda_f;

protected:

  // parameters

  D mu_f;

  D mu_psi;

  D mu_g;

  L n;   // x\in \R^n

  L tau;  //number of threads on each node/computer

  D sumofLi;

  // variables
  std::vector<D> u;

  std::vector<D> z;

  std::vector<D> x;

  std::vector<D> v;

  std::vector<D> L_M;

  std::vector<D> L_f;

  std::vector<D> t;

  D gamma;

  D theta;

  D theta0;

  L m_1;

  L m_2;

  // sampling variables

  std::vector<D> proba_vector;

  std::vector<D> S;

  std::vector<D> all_n;

  std::vector<D> sampled;

  D max_p;

  D min_p;

  // auxiliary variables

  L max_nb_loops;

  L evaluation;

  ofstream samp;

  L nb_iters;

  L nb_of_iters_per_loop;

  L print_every_N;

  D running_time;

  Matrix<L,D> data_A;

  Matrix<L, D> data_M;

  std::vector<D> lambda;

  std::vector<D> M_tlambda;

  D beta_s;

  D function_value;

  D residual1;

  D residual2;

public:

  gsl_rng * rng;

  virtual inline D gradient_of_phi_j(D, L){return D(NULL);}

  virtual inline D value_of_g_i(D, L){return D(NULL);}
  virtual inline D value_of_phi_j(D, L){return D(NULL);}

  virtual inline D prox_of_g_i(D,D,D, L){return D(NULL);}
  virtual inline D value_of_phistar_i(D,L) {return D(NULL);}

  virtual inline void set_matrix_M(){}

  virtual inline void set_matrix_A(){}
  PDCD_SMSVM()
  {

  }

  void set_rng()
  {
    gsl_rng_env_setup();
    const gsl_rng_type * T;
    T = gsl_rng_default;
    rng = gsl_rng_alloc(T);
    gsl_rng_set(rng,time(NULL));
    //gsl_rng_set(rng, 27432042);

  }

// sample i with probability pi=proba_vector[i]
  L sampling()
  {
    //L i=(floor)(gsl_rng_uniform(rng)*n);
    L i=gsl_rng_uniform_int(rng, n);
    if(tau==1)
    {
      D y=gsl_rng_uniform(rng);
      while(y*max_p>proba_vector[i])
      {
        i=(floor)(gsl_rng_uniform(rng)*n);
        y=gsl_rng_uniform(rng);
      }
    }
    return i;
  }

// sample S
  void batch_sampling()
  {
    if(tau<n)
    {
      L i=sampling();
      for(L k=0;k<tau;k++)
      {
        while(sampled[i]==1)
        {
          i=sampling();
        }
        sampled[i]=1;
        S[k]=i;
      }
      for(L k=0;k<tau;k++)
      {
        sampled[S[k]]=0;
      }
    }
    else {
      S=all_n;
      //cout<<"s=all_n"<<endl;
    }
  }

  void compute_x()
  {
    for(L i=0;i<n;i++){
    x[i]=gamma*u[i]+z[i];
	}
	for(L j=0;j<m_2;j++){
        Mx[j]=gamma*Mu[j]+Mz[j];
    }
    for(L j=0;j<m_1;j++){
        Ax[j]=gamma*Au[j]+Az[j];
    }
  }

  inline void compute_function_value() {
     D res=0;
     D res2= 0;
     D tmp1= 0;
     for(L i=0;i<this->n;i++){
       res+=value_of_g_i(x[i],i);
     }
     for(L j=0;j<m_1;j++){
       res+=lambda_f[j]*value_of_phi_j(Ax[j],j);
     }
     for(L j=0;j<m_2;j++){
     	if (1- Mx[j]> tmp1){
     		tmp1= 1- Mx[j];
		 }
        res2+=max(0.0,1- Mx[j]);
     }
     residual1= tmp1;
     residual2= res2/m_2;
     function_value=res+ m_2*residual2;;
   }

  void initialize(D beta_0, vector<D> & x0, vector<D> & y0, L val_tau, D p_N, D val_lambda_f)
  {
  	set_matrix_M();
    set_matrix_A();

    this->tau=val_tau;
    m_1=data_A.get_n();
    m_2=data_M.get_n();
    cout<<"m_1="<<m_1<<endl;
    cout<<"m_2="<<m_2<<endl;
    lambda_f.resize(m_1,val_lambda_f);
    n=data_A.nfeatures;
    tau=val_tau;
    all_n.resize(n,0);
    for(L i=0;i<n;i++)
    all_n[i]=i;
    print_every_N= p_N;
    nb_of_iters_per_loop=floor(max(1.,n/(tau+0.0)));
    lambda.resize(m_2,0);
    for(L j=0;j<m_2;j++)
           lambda[j]=y0[j];

    u.clear();
    u.resize(n,0);
    z.clear();
    z.resize(n,0);
    x.clear();
    x.resize(n,0);
    for(L i=0;i<n;i++)
    {
      z[i]=x0[i];
      x[i]=x0[i];
    }
    sampled.clear();
    sampled.resize(n,0);
    S.clear();
    S.resize(tau,0);
    t.clear();
    t.resize(n,0);
    M_tlambda.resize(n,0);
    compute_Mty();

    Au.clear();
    Au.resize(m_1,0);
    Az.clear();
    Az.resize(m_1,0);
    Ax.clear();
    Ax.resize(m_1,0);
    Mu.clear();
    Mu.resize(m_2,0);
    Mz.clear();
    Mz.resize(m_2,0);
    Mx.clear();
    Mx.resize(m_2,0);
    compute_Az(x0);
    compute_Mz(x0);
    compute_Ax(x0);
    compute_Mx(x0);
    beta_s= beta_0;
    set_v();
    set_p();
    set_rng();

    theta0= min_p;
    theta= theta0;
    cout<< "theta0= "<< theta0<< endl;
    gamma= 1- theta0;

    cout<<"finished SMART_CD initializing"<<endl;
  }

  inline void compute_and_record_res(){
  	    if (nb_iters%print_every_N== 0){
  	    cout<< "gamma= "<< gamma<< " theta= "<< theta<< " beta_s= "<< beta_s<< endl;
        //compute_KKT_residual();
        compute_function_value();
        cout<<setprecision(9)<<"Iteration: "<<nb_iters<<"; time="<<running_time<<"; infeasibility= "<< residual1<<" "<< residual2<<"; function value="<<function_value<< endl;
        samp<<setprecision(9)<< nb_iters<<" "<<running_time<<" "<<residual1<< " "<< residual2<< " "<< function_value<<" "<<endl;
        }
   }

  inline D partial_i_of_f(L i)
  {
  	D res=0;
      for (L k = data_A.ptr_t[i]; k < data_A.ptr_t[i + 1];k++)
      {
        L j=data_A.col_idx[k];
        D tmp=lambda_f[j]*gradient_of_phi_j(gamma*Au[j]+Az[j], j);
        res+=data_A.A_t[k]*tmp;
      }
      for (L k = data_M.ptr_t[i]; k < data_M.ptr_t[i + 1];k++)
      {
        L j=data_M.col_idx[k];
        D tmp=(gamma*Mu[j]+Mz[j]- 1)/beta_s+lambda[j];
        if (tmp> 0){
        	tmp= 0;
		}
		else if( tmp< -1){
			tmp= -1;
		}
        res+=data_M.A_t[k]*tmp;
      }
      return res;
  }

  inline void set_v()
   {
     v.resize(n,0);
     L_f.resize(n,0);
     L_M.resize(n,0);
     D maxv=0;
     D minv=std::numeric_limits<double>::max();
     D sumv=0;
     D sumvi1=0;
     L sumw=0;
     L maxw=0;
     L minw=n;
     sumofLi=0;
     for(L j=0;j<m_1;j++)
     {
       sumw+=data_A.w_t[j];
       maxw=max(maxw,data_A.w_t[j]);
       minw=min(minw,data_A.w_t[j]);
     }
     cout<<"sumw="<<sumw<<";  maxw="<<maxw<<"; minw="<<minw<<endl;
     for(L i=0;i<n;i++)
     {
       D lfi=0;
       D lmi=0;
       D lfi1= 0;
       D lmi1= 0;
       for (L k = data_A.ptr_t[i]; k < data_A.ptr_t[i + 1];k++)
       {
         L j=data_A.col_idx[k];
         lfi+=(1.+(data_A.w_t[j]-1.)*(tau-1.)/max(n-1.,1.))*data_A.A_t[k]*data_A.A_t[k]*lambda_f[j];
         lfi1+=data_A.A_t[k]*data_A.A_t[k]*lambda_f[j];
       }
       for (L k = data_M.ptr_t[i]; k < data_M.ptr_t[i + 1];k++)
       {
         L j=data_M.col_idx[k];
         lmi+=(1.+(data_M.w_t[j]-1.)*(tau-1.)/max(n-1.,1.))*data_M.A_t[k]*data_M.A_t[k];
         lmi1+=data_M.A_t[k]*data_M.A_t[k];
       }
       L_f[i]=lfi;
       L_M[i]= lmi;
       sumv+=lfi;
       sumvi1+=lfi1;
       if(maxv<lfi) maxv=lfi;
       if(minv>lfi) minv=lfi;
     }
     if(tau==n){
          for(L i=0;i<n;i++)
              L_f[i]=sumvi1;
     }
     sumofLi=sumvi1;
     cout<<"  max of v: "<<maxv<<" ;  min of v: "<<minv<<" ;  sumofv: "<<sumv<<" sumofLi="<<sumofLi<<endl;
   }

   inline void set_p(){
     proba_vector.resize(n,0.0);
     D res= 0;
     D tmp= 1;
     D tmp2= 0;
     for (L i=0; i< n; i++){
     	res+= sqrt(L_f[i]+ L_M[i]/beta_s);
	 }
	 for (L i= 0; i< n; i++){
	 	proba_vector[i]= sqrt(L_f[i]+ L_M[i]/beta_s)/res;
	 	if (proba_vector[i]< tmp){
	 		tmp= proba_vector[i];
		 }
		if (proba_vector[i]> tmp2){
	 		tmp2= proba_vector[i];
		 }
	 }
	 if (tmp== 0){
	 	for (L i= 0; i< n; i++){
	 	proba_vector[i]= 1.0/n;
	    }
	    min_p= 1.0/n;
	    max_p= 1.0/n;
	 }
	 else{
	    min_p= tmp;
	    max_p= tmp2;
     }
	 cout<< "max_p= "<< max_p<< "; min_p= "<< min_p<< endl;
   }

  inline void update_z_coordinate( L i, D dz){
     z[i]+=dz;
     for (L k = data_A.ptr_t[i]; k < data_A.ptr_t[i + 1];k++)
     {
       L j=data_A.col_idx[k];
       Az[j]+=dz*data_A.A_t[k];
     }
     for (L k = data_M.ptr_t[i]; k < data_M.ptr_t[i + 1];k++)
     {
       L j=data_M.col_idx[k];
       Mz[j]+=dz*data_M.A_t[k];
       //if(j==2869) cout<<"gamma="<<gamma<<"; "<<dz<<"; "<<Mz[j]<<endl;
     }
   }

    inline void update_x_coordinate( L i, D dx){
     x[i]+=dx;
     x[i]=x[i];
     L j;
     for (L k = data_M.ptr_t[i]; k < data_M.ptr_t[i + 1];k++)
     {
       j=data_M.col_idx[k];
       Mx[j]+=dx*data_M.A_t[k];
     }
      for (L k = data_A.ptr_t[i]; k < data_A.ptr_t[i + 1];k++)
     {
       j=data_A.col_idx[k];
       Ax[j]+=dx*data_A.A_t[k];
     }
   }

   inline void update_u_coordinate( L i, D du){
     u[i]+=du;
     for (L k = data_A.ptr_t[i]; k < data_A.ptr_t[i + 1];k++)
     {
       L j=data_A.col_idx[k];
       Au[j]+=du*data_A.A_t[k];
     }
     for (L k = data_M.ptr_t[i]; k < data_M.ptr_t[i + 1];k++)
     {
       L j=data_M.col_idx[k];
       Mu[j]+=du*data_M.A_t[k];
     }
   }

   void compute_Au(){
     for(L j=0;j<m_1;j++)
       for(L k = data_A.ptr[j]; k < data_A.ptr[j + 1];k++){
         L kj=data_A.row_idx[k];
         Au[j]+=u[kj]*data_A.A[k];
       }
   }

   void compute_Mu(){
     for(L j=0;j<m_2;j++)
       for(L k = data_M.ptr[j]; k < data_M.ptr[j + 1];k++){
         L kj=data_M.row_idx[k];
         Mu[j]+=u[kj]*data_M.A[k];
       }
   }

   void compute_Az(){
     for(L j=0;j<m_1;j++)
       for(L k = data_A.ptr[j]; k < data_A.ptr[j + 1];k++){
         L kj=data_A.row_idx[k];
         Az[j]+=z[kj]*data_A.A[k];
       }
   }
   void compute_Mz(){
     for(L j=0;j<m_2;j++)
       for(L k = data_M.ptr[j]; k < data_M.ptr[j + 1];k++){
         L kj=data_M.row_idx[k];
         Mz[j]+=z[kj]*data_M.A[k];
       }
   }

   void compute_Au(vector<D> & x0){
     for(L j=0;j<m_1;j++)
       for(L k = data_A.ptr[j]; k < data_A.ptr[j + 1];k++){
         L kj=data_A.row_idx[k];
         Au[j]+=x0[kj]*data_A.A[k];
       }

   }

   void compute_Az(vector<D> & x0){
     for(L j=0;j<m_1;j++)
       for(L k = data_A.ptr[j]; k < data_A.ptr[j + 1];k++){
         L kj=data_A.row_idx[k];
         Az[j]+=x0[kj]*data_A.A[k];
       }
   }
   void compute_Mu(vector<D> & x0){
     for(L j=0;j<m_2;j++)
       for(L k = data_M.ptr[j]; k < data_M.ptr[j + 1];k++){
         L kj=data_M.row_idx[k];
         Mu[j]+=x0[kj]*data_M.A[k];
       }

   }

   void compute_Mz(vector<D> & x0){
     for(L j=0;j<m_2;j++)
       for(L k = data_M.ptr[j]; k < data_M.ptr[j + 1];k++){
         L kj=data_M.row_idx[k];
         Mz[j]+=x0[kj]*data_M.A[k];
       }
   }

     void compute_Mx(vector<D> & x0){
     for(L j=0;j<m_2;j++)
       for(L k = data_M.ptr[j]; k < data_M.ptr[j + 1];k++){
         L kj=data_M.row_idx[k];
         Mx[j]+=x0[kj]*data_M.A[k];
       }
   }

   void compute_Ax(vector<D> & x0){
     for(L j=0;j<m_1;j++)
       for(L k = data_A.ptr[j]; k < data_A.ptr[j + 1];k++){
         L kj=data_A.row_idx[k];
         Ax[j]+=x0[kj]*data_A.A[k];
       }
   }

  void compute_Mty(){
        M_tlambda.clear();
        M_tlambda.resize(this->n,0);
        for(L i=0;i<this->n;i++)
        {
          for(L k = data_M.ptr_t[i]; k < data_M.ptr_t[i + 1];k++){
            L j=data_M.col_idx[k];
            M_tlambda[i]+=lambda[j]*data_M.A_t[k];
           }
        }

    }

    void update_theta()
  {

	 theta= (sqrt(theta*theta*theta*theta+ 4*theta*theta)- theta*theta)/2;

  }

  void reset(){
  	beta_s= beta_s/(1+ theta);
  }

  void PDCD_SMSVM_solver(D beta_0, vector<D> & x0,vector<D> & y0, L val_tau, L max_nb_epoch, L p_N, D val_lambda_f,string filename1, D time2)
  {
    initialize(beta_0, x0, y0, val_tau,  p_N, val_lambda_f);
    cout<<"running SMART_CD"<<" ; "<<filename1<<" max_nb_epoch "<<max_nb_epoch<<endl;
    nb_iters=0;
    filename1="results/PDCD_SMSVM_"+filename1;
    samp.open(filename1.c_str());
    //srand48(27432042);
    srand(time(NULL));
    running_time= 0;
    print_every_N=p_N;
    compute_and_record_res();
    D start;
    start = std::clock();
    while(nb_iters<max_nb_epoch)
    {
      start = std::clock();
      D ti=0;
      //D si=0;
      D gi=0;
      D Li=0;
      L i=0;
      for(L it=0;it<nb_of_iters_per_loop;it++)
      {
        gamma*=(1-theta);
        if(theta==1) gamma=1;
        batch_sampling();
        for(L it_S=0;it_S<tau;it_S++)
        {
          i=S[it_S];
          gi=partial_i_of_f(i);
          Li=(L_f[i]+ L_M[i]/beta_s)*theta/theta0;
          t[i]=prox_of_g_i(gi, Li, z[i],i);
        }
        for(L it_S=0;it_S<tau;it_S++)
        {
          i=S[it_S];
          ti=t[i];
          update_z_coordinate(i, ti);
          update_u_coordinate(i, -(1- theta/theta0)/gamma*ti);
        }
        update_theta();
        reset();
      }
      running_time+=( std::clock() - start ) / (double) CLOCKS_PER_SEC;
      nb_iters++;
      compute_x();
      compute_and_record_res();
      if (running_time> time2){
      	break;
	  }
    }

  }

};

#endif


text:
/*Python3 includes*/
#include <Python.h>

/*System includes*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <sys/stat.h>
#include <float.h>
#include <sys/time.h>

/*GSL includes*/
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_roots.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_statistics_double.h>
#include <gsl/gsl_fft_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_deriv.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_cblas.h>
#include <pthread.h>

/*User includes*/
#include "vbgmmmodule.h"

char szMLInputFile[BUFSIZ]; // mustlink file path

// Module method definition
static PyObject *get_wrapper_version(PyObject *self, PyObject *args)
{
	(void)self;
	(void)args;
	return Py_BuildValue("s", "0.3");
}

static PyObject *vbgmm_get_n_jobs(PyObject *self, PyObject *args)
{
	(void)self;
	(void)args;
	return PyLong_FromLong(N_RTHREADS);
}

static PyObject *vbgmm_testfit(PyObject *self, PyObject *args)
{
	(void)self;
	const char *szFileStub = "/home/kgravouil/THESE/data/FENNEC/fennec/c-concoct/tests";
	int nKStart = 60, nLMin = 1000, nMaxIter = 30;
	unsigned long lSeed = 666;
	double dEpsilon = 1e-1;
	int bCout = 0;
	int sts;

	sts = driver(szFileStub, nKStart, nLMin, lSeed, nMaxIter, dEpsilon, bCout);
	return Py_BuildValue("i", sts);
}

static PyObject *vbgmm_fit(PyObject *self, PyObject *args)
{
	(void)self;
	const char *szFileStub = NULL;
	int nKStart = 0, nLMin = 0, nMaxIter;
	unsigned long lSeed;
	double dEpsilon;
	int bCout = 0;
	int sts;

	if (!PyArg_ParseTuple(args, "siikidi", &szFileStub, &nKStart, &nLMin, &lSeed, &nMaxIter,
						  &dEpsilon, &bCout))
		return NULL;
	//TODO: return the np.array of clustering result
	sts = driver(szFileStub, nKStart, nLMin, lSeed, nMaxIter, dEpsilon, bCout);
	return Py_BuildValue("i", sts);
}

//Method definition object for this extension, these argumens mean:
//ml_name: The name of the method
//ml_meth: Function pointer to the method implementation
//ml_flags: Flags indicating special features of this method, such as
//          accepting arguments, accepting keyword arguments, being a
//          class method, or being a static method of a class.
//ml_doc:  Contents of this method's docstring
static PyMethodDef vbgmm_methods[] = {
	{"get_wrapper_version",
	 get_wrapper_version,
	 METH_NOARGS,
	 "Return the python wrapper version"},
	{"get_n_jobs",
	 vbgmm_get_n_jobs,
	 METH_NOARGS,
	 "Return number of jobs to run."},
	{"fit",
	 vbgmm_fit,
	 METH_VARARGS,
	 "Fit a variational Bayesian Gaussian mixture."},
	{"test_fit",
	 vbgmm_testfit,
	 METH_VARARGS,
	 "Quick test of fit function."},
	{NULL, NULL, 0, NULL} /* Sentinel */
};

//Module definition
//The arguments of this structure tell Python what to call your extension,
//what its methods are and where to look for its method definitions
static struct PyModuleDef vbgmm_definition = {
	PyModuleDef_HEAD_INIT,
	"vbgmm",
	"A Python module that compute VBGMM",
	-1,
	vbgmm_methods};

//Module initialization
//Python calls this function when importing your extension. It is important
//that this function is named PyInit_[[your_module_name]] exactly, and matches
//the name keyword argument in setup.py's setup() call.
PyMODINIT_FUNC PyInit_vbgmm(void)
{
	Py_Initialize();
	return PyModule_Create(&vbgmm_definition);
}

/*** core functions ***/
int driver(const char *szFileStub, int nKStart, int nLMin, unsigned long lSeed,
		   int nMaxIter, double dEpsilon, int bCOut)
{
	t_Params tParams;
	t_Data tData;
	gsl_rng *ptGSLRNG = NULL;
	const gsl_rng_type *ptGSLRNGType = NULL;
	int nD = 0, nN = 0;
	char szOFile[MAX_LINE_LENGTH];
	t_VBParams tVBParams;
	t_Cluster *ptBestCluster = NULL;
	gsl_matrix *ptTemp = NULL;
	gsl_matrix *ptTVar = NULL;

	sprintf(szMLInputFile, "%s%s%d.dat", szFileStub, MLINPUT_FILE, nLMin);

	/*initialise GSL RNG */
	gsl_rng_env_setup();

	gsl_set_error_handler_off();

	ptGSLRNGType = gsl_rng_default;
	ptGSLRNG = gsl_rng_alloc(ptGSLRNGType);

	/*get command line params */
	tParams.nKStart = nKStart;
	tParams.nLMin = nLMin;
	tParams.nMaxIter = nMaxIter;
	tParams.dEpsilon = dEpsilon;
	tParams.lSeed = lSeed;

	setParams(&tParams, szFileStub);

	/*read in input data */
	readInputData(tParams.szInputFile, &tData);
	readPInputData(tParams.szPInputFile, &tData);

	nD = tData.nD;
	nN = tData.nN;

	ptTemp = gsl_matrix_alloc(tData.nT, nD);
	ptTVar = gsl_matrix_alloc(tData.nT, tData.nT);

	setVBParams(&tVBParams, &tData);

	ptBestCluster = (t_Cluster *)malloc(sizeof(t_Cluster));

	ptBestCluster->nN = nN;
	ptBestCluster->nK = tParams.nKStart;
	ptBestCluster->nD = nD;
	ptBestCluster->ptData = &tData;
	ptBestCluster->ptVBParams = &tVBParams;
	ptBestCluster->lSeed = tParams.lSeed;
	ptBestCluster->nMaxIter = tParams.nMaxIter;
	ptBestCluster->dEpsilon = tParams.dEpsilon;

	if (bCOut > 0)
		ptBestCluster->szCOutFile = (char *)szFileStub;
	else
		ptBestCluster->szCOutFile = NULL;

	runRThreads((void *)&ptBestCluster);

	compressCluster(ptBestCluster);

	calcCovarMatrices(ptBestCluster, &tData);

	// printing clustering results
	printf("[DEBUG] Printing clustering\n");
	sprintf(szOFile, "%sclustering_gt%d.csv", tParams.szOutFileStub,
			tParams.nLMin);
	writeClusters(szOFile, ptBestCluster, &tData);

	/*free up memory in data object */
	printf("[DEBUG] Cleaning memory\n");
	destroyData(&tData);

	/*free up best BIC clusters */
	destroyCluster(ptBestCluster);
	free(ptBestCluster);

	destroyParams(&tParams);
	gsl_rng_free(ptGSLRNG);
	gsl_matrix_free(tVBParams.ptInvW0);
	gsl_matrix_free(ptTemp);
	gsl_matrix_free(ptTVar);

	printf("[INFO] done\n");
	return EXIT_SUCCESS;
}

void setParams(t_Params *ptParams, const char *szFileStub)
{
	ptParams->szInputFile = (char *)malloc(MAX_LINE_LENGTH * sizeof(char));
	if (!ptParams->szInputFile)
		goto memoryError;

	ptParams->szPInputFile = (char *)malloc(MAX_LINE_LENGTH * sizeof(char));
	if (!ptParams->szPInputFile)
		goto memoryError;

	sprintf(ptParams->szInputFile, "%s%s%d.csv", szFileStub, INPUT_FILE,
			ptParams->nLMin);

	sprintf(ptParams->szPInputFile, "%s%s%d.csv", szFileStub, PINPUT_FILE,
			ptParams->nLMin);

	ptParams->szOutFileStub = szFileStub;

	return;

memoryError:
	fprintf(stderr, "Failed allocating memory in setParams\n");
	fflush(stderr);
	exit(EXIT_FAILURE);
}

void destroyParams(t_Params *ptParams)
{
	free(ptParams->szInputFile);
	free(ptParams->szPInputFile);
}

void setVBParams(t_VBParams *ptVBParams, t_Data *ptData)
{
	int i = 0, nD = ptData->nD;
	double adVar[nD], adMu[nD];

	ptVBParams->dBeta0 = DEF_BETA0;
	ptVBParams->dNu0 = (double)nD;
	ptVBParams->ptInvW0 = gsl_matrix_alloc(nD, nD);

	calcSampleVar(ptData, adVar, adMu);
	gsl_matrix_set_zero(ptVBParams->ptInvW0);

	for (i = 0; i < nD; i++)
	{
		double dRD = adVar[i] * ((double)nD);

		gsl_matrix_set(ptVBParams->ptInvW0, i, i, dRD);
	}

	ptVBParams->dLogWishartB =
		dLogWishartB(ptVBParams->ptInvW0, nD, ptVBParams->dNu0, TRUE);
}

void readInputData(const char *szFile, t_Data *ptData)
{
	double **aadX = NULL;
	int i = 0, j = 0, nD = 0, nN = 0;
	char szLine[MAX_LINE_LENGTH];
	FILE *ifp = NULL;

	ifp = fopen(szFile, "r");

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

		if (fgets(szLine, MAX_LINE_LENGTH, ifp) == NULL)
			goto formatError;

		szTok = strtok(szLine, DELIM);
		/*count dimensions */
		while (strtok(NULL, DELIM) != NULL)
		{
			nD++;
		}
		/*count data points */
		while (fgets(szLine, MAX_LINE_LENGTH, ifp) != NULL)
		{
			nN++;
		}
		fclose(ifp);

		/*reopen input file */
		ifp = fopen(szFile, "r");

		if (fgets(szLine, MAX_LINE_LENGTH, ifp) == NULL)
			goto formatError;

		/*allocate memory for dimension names */
		ptData->aszDimNames = (char **)malloc(nD * sizeof(char *));
		if (!ptData->aszDimNames)
			goto memoryError;

		szTok = strtok(szLine, DELIM);
		/*read in dim names */
		for (i = 0; i < nD; i++)
		{
			szTok = strtok(NULL, DELIM);
			ptData->aszDimNames[i] = strdup(szTok);
		}

		/*allocate memory for data matrix */
		aadX = (double **)malloc(nN * sizeof(double *));
		if (!aadX)
			goto memoryError;
		for (i = 0; i < nN; i++)
		{
			aadX[i] = (double *)malloc(nD * sizeof(double));
			if (!aadX[i])
				goto memoryError;
		}

		/*read in input data */
		ptData->aszSampleNames = (char **)malloc(nN * sizeof(char *));
		if (!ptData->aszSampleNames)
			goto memoryError;

		for (i = 0; i < nN; i++)
		{

			if (fgets(szLine, MAX_LINE_LENGTH, ifp) == NULL)
				goto formatError;

			szTok = strtok(szLine, DELIM);
			ptData->aszSampleNames[i] = strdup(szTok);
			for (j = 0; j < nD; j++)
			{
				szTok = strtok(NULL, DELIM);

				aadX[i][j] = strtod(szTok, &pcError);

				if (*pcError != '\0')
				{
					goto formatError;
				}
			}
		}
		fclose(ifp);
	}
	else
	{
		fprintf(stderr,
				"Failed to open abundance data file %s aborting\n",
				szFile);
		fflush(stderr);
		exit(EXIT_FAILURE);
	}

	ptData->nD = nD;
	ptData->nN = nN;
	ptData->aadX = aadX;
	return;

memoryError:
	fprintf(stderr, "Failed allocating memory in readInputData\n");
	fflush(stderr);
	exit(EXIT_FAILURE);

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

void readPInputData(const char *szFile, t_Data *ptData)
{
	int i = 0, j = 0, nD = ptData->nD, nT = 0;
	char szLine[MAX_LINE_LENGTH];
	FILE *ifp = NULL;

	ifp = fopen(szFile, "r");

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

		if (fgets(szLine, MAX_LINE_LENGTH, ifp) == NULL)
			goto memoryError;

		szTok = strtok(szLine, DELIM);
		/*count dimensions */
		while (strtok(NULL, DELIM) != NULL)
		{
			nT++;
		}

		ptData->ptTMatrix = gsl_matrix_alloc(nT, nD);

		fclose(ifp);

		/*reopen input file */
		ifp = fopen(szFile, "r");

		for (i = 0; i < nD; i++)
		{
			double dTemp = 0.0;

			if (fgets(szLine, MAX_LINE_LENGTH, ifp) == NULL)
				goto formatError;

			szTok = strtok(szLine, DELIM);

			dTemp = strtod(szTok, &pcError);
			if (*pcError != '\0')
			{
				goto formatError;
			}

			gsl_matrix_set(ptData->ptTMatrix, 0, i, dTemp);

			for (j = 1; j < nT; j++)
			{
				szTok = strtok(NULL, DELIM);

				dTemp = strtod(szTok, &pcError);
				if (*pcError != '\0')
				{
					goto formatError;
				}

				gsl_matrix_set(ptData->ptTMatrix, j, i, dTemp);
			}
		}
		fclose(ifp);
	}
	else
	{
		fprintf(stderr,
				"Failed to open abundance data file %s aborting\n",
				szFile);
		fflush(stderr);
		exit(EXIT_FAILURE);
	}

	ptData->nT = nT;
	return;

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

memoryError:
	fprintf(stderr, "Failed allocating memory in readPInputFile\n");
	fflush(stderr);
	exit(EXIT_FAILURE);
}

void destroyData(t_Data *ptData)
{
	int nN = ptData->nN, nD = ptData->nD;
	int i = 0;

	gsl_matrix_free(ptData->ptTMatrix);

	for (i = 0; i < nD; i++)
	{
		free(ptData->aszDimNames[i]);
	}
	free(ptData->aszDimNames);

	for (i = 0; i < nN; i++)
	{
		free(ptData->aadX[i]);
	}
	free(ptData->aadX);

	for (i = 0; i < nN; i++)
	{
		free(ptData->aszSampleNames[i]);
	}
	free(ptData->aszSampleNames);

	return;
}

void destroyCluster(t_Cluster *ptCluster)
{
	int i = 0, nN = ptCluster->nN, nKSize = ptCluster->nKSize;

	if (ptCluster->szCOutFile != NULL)
	{
		free(ptCluster->szCOutFile);
	}

	free(ptCluster->anMaxZ);

	free(ptCluster->anW);

	for (i = 0; i < nN; i++)
	{
		free(ptCluster->aadZ[i]);
	}
	free(ptCluster->aadZ);

	free(ptCluster->adLDet);
	free(ptCluster->adPi);
	free(ptCluster->adBeta);
	free(ptCluster->adNu);

	for (i = 0; i < nKSize; i++)
	{
		free(ptCluster->aadMu[i]);
		free(ptCluster->aadM[i]);
	}

	free(ptCluster->aadMu);
	free(ptCluster->aadM);

	for (i = 0; i < nKSize; i++)
	{
		gsl_matrix_free(ptCluster->aptSigma[i]);
		gsl_matrix_free(ptCluster->aptCovar[i]);
	}
	free(ptCluster->aptSigma);
	free(ptCluster->aptCovar);
	return;
}

void allocateCluster(t_Cluster *ptCluster, int nN, int nK, int nD,
					 t_Data *ptData, long lSeed, int nMaxIter, double dEpsilon,
					 char *szCOutFile)
{
	int i = 0, j = 0, k = 0;

	ptCluster->szCOutFile = szCOutFile;
	ptCluster->ptVBParams = NULL;
	ptCluster->lSeed = lSeed;
	ptCluster->nMaxIter = nMaxIter;
	ptCluster->dEpsilon = dEpsilon;
	ptCluster->ptData = ptData;

	ptCluster->nN = nN;
	ptCluster->nK = nK;
	ptCluster->nKSize = nK;
	ptCluster->nD = nD;

	ptCluster->dVBL = 0.0;

	ptCluster->anMaxZ = (int *)malloc(nN * sizeof(int)); /*destroyed */
	if (!ptCluster->anMaxZ)
		goto memoryError;

	ptCluster->anW = (int *)malloc(nK * sizeof(int)); /*destroyed */
	if (!ptCluster->anW)
		goto memoryError;

	for (i = 0; i < nN; i++)
	{
		ptCluster->anMaxZ[i] = NOT_SET;
	}

	for (i = 0; i < nK; i++)
	{
		ptCluster->anW[i] = 0;
	}

	ptCluster->aadZ = (double **)malloc(nN * sizeof(double *)); /*destroyed */
	if (!ptCluster->aadZ)
		goto memoryError;

	for (i = 0; i < nN; i++)
	{
		ptCluster->aadZ[i] = (double *)malloc(nK * sizeof(double)); /*destroyed */
		if (!ptCluster->aadZ[i])
			goto memoryError;

		for (j = 0; j < nK; j++)
		{
			ptCluster->aadZ[i][j] = 0.0;
		}
	}

	ptCluster->adLDet = (double *)malloc(nK * sizeof(double)); /*all */
	ptCluster->adPi = (double *)malloc(nK * sizeof(double));
	ptCluster->adBeta = (double *)malloc(nK * sizeof(double));
	ptCluster->adNu = (double *)malloc(nK * sizeof(double)); /*destroyed */

	if (!ptCluster->adLDet || !ptCluster->adPi)
		goto memoryError;

	if (!ptCluster->adBeta || !ptCluster->adNu)
		goto memoryError;

	for (k = 0; k < nK; k++)
	{
		ptCluster->adLDet[k] = 0.0;
		ptCluster->adPi[k] = 0.0;
		ptCluster->adBeta[k] = 0.0;
		ptCluster->adNu[k] = 0.0;
	}

	ptCluster->aadMu = (double **)malloc(nK * sizeof(double *));
	if (!ptCluster->aadMu)
		goto memoryError;

	ptCluster->aadM = (double **)malloc(nK * sizeof(double *));
	if (!ptCluster->aadM)
		goto memoryError;

	for (i = 0; i < nK; i++)
	{
		ptCluster->aadM[i] = (double *)malloc(nD * sizeof(double));
		if (!ptCluster->aadM[i])
			goto memoryError;

		ptCluster->aadMu[i] = (double *)malloc(nD * sizeof(double));
		if (!ptCluster->aadMu[i])
			goto memoryError;
	}

	ptCluster->aptSigma = (gsl_matrix **)malloc(nK * sizeof(gsl_matrix *));
	if (!ptCluster->aptSigma)
		goto memoryError;

	for (i = 0; i < nK; i++)
	{
		ptCluster->aptSigma[i] =
			(gsl_matrix *)gsl_matrix_alloc(nD, nD);
	}

	ptCluster->aptCovar = (gsl_matrix **)malloc(nK * sizeof(gsl_matrix *));
	if (!ptCluster->aptCovar)
		goto memoryError;

	for (i = 0; i < nK; i++)
	{
		ptCluster->aptCovar[i] =
			(gsl_matrix *)gsl_matrix_alloc(nD, nD);
	}

	return;

memoryError:
	fprintf(stderr, "Failed allocating memory in allocateCluster\n");
	fflush(stderr);
	exit(EXIT_FAILURE);
}

void writeSquareMatrix(char *szFile, gsl_matrix *ptMatrix, int nD)
{
	int i = 0, j = 0;
	FILE *ofp = fopen(szFile, "w");

	if (ofp)
	{
		for (i = 0; i < nD; i++)
		{
			for (j = 0; j < nD - 1; j++)
			{
				fprintf(ofp, "%f,",
						gsl_matrix_get(ptMatrix, i, j));
			}
			fprintf(ofp, "%f\n", gsl_matrix_get(ptMatrix, i, j));
		}
	}
	else
	{
		fprintf(stderr,
				"Failed to open %s for writing in writeSquareMatrix\n",
				szFile);
		fflush(stderr);
	}
}

double decomposeMatrix(gsl_matrix *ptSigmaMatrix, int nD)
{
	double dDet = 0.0;
	int status;
	int l = 0;

	status = gsl_linalg_cholesky_decomp(ptSigmaMatrix);

	if (status == GSL_EDOM)
	{
		fprintf(stderr,
				"Failed Cholesky decomposition in decomposeMatrix\n");
		fflush(stderr);
		exit(EXIT_FAILURE);
	}
	else
	{
		for (l = 0; l < nD; l++)
		{
			double dT = gsl_matrix_get(ptSigmaMatrix, l, l);
			dDet += 2.0 * log(dT);
		}
		gsl_linalg_cholesky_invert(ptSigmaMatrix);
		return dDet;
	}
}

void performMStep(t_Cluster *ptCluster, t_Data *ptData)
{
	int i = 0, j = 0, k = 0, l = 0, m = 0;
	int nN = ptData->nN, nK = ptCluster->nK, nD = ptData->nD;
	double **aadZ = ptCluster->aadZ, **aadX = ptData->aadX;
	double *adLDet = ptCluster->adLDet, *adPi = ptCluster->adPi;
	double **aadCovar = NULL, **aadInvWK = NULL;
	t_VBParams *ptVBParams = ptCluster->ptVBParams;

	// gsl_matrix *ptSigmaMatrix = gsl_matrix_alloc(nD, nD);

	aadCovar = (double **)malloc(nD * sizeof(double *));
	if (!aadCovar)
		goto memoryError;

	for (i = 0; i < nD; i++)
	{
		aadCovar[i] = (double *)malloc(nD * sizeof(double));
		if (!aadCovar[i])
			goto memoryError;
	}

	aadInvWK = (double **)malloc(nD * sizeof(double *));
	if (!aadInvWK)
		goto memoryError;

	for (i = 0; i < nD; i++)
	{
		aadInvWK[i] = (double *)malloc(nD * sizeof(double));
		if (!aadInvWK[i])
			goto memoryError;
	}

	/*perform M step */
	for (k = 0; k < nK; k++)
	{ /*loop components */
		double *adMu = ptCluster->aadMu[k];
		gsl_matrix *ptSigmaMatrix = ptCluster->aptSigma[k];
		double dF = 0.0;
		/*recompute mixture weights and means */
		for (j = 0; j < nD; j++)
		{
			adMu[j] = 0.0;
			for (l = 0; l < nD; l++)
			{
				aadCovar[j][l] = 0.0;
			}
		}

		/* compute weight associated with component k */
		adPi[k] = 0.0;
		for (i = 0; i < nN; i++)
		{
			if (aadZ[i][k] > MIN_Z)
			{
				adPi[k] += aadZ[i][k];
				for (j = 0; j < nD; j++)
				{
					adMu[j] += aadZ[i][k] * aadX[i][j];
				}
			}
		}

		/*normalise means */
		if (adPi[k] > MIN_PI)
		{
			/*Equation 10.60 */
			ptCluster->adBeta[k] = ptVBParams->dBeta0 + adPi[k];

			for (j = 0; j < nD; j++)
			{
				/*Equation 10.61 */
				ptCluster->aadM[k][j] =
					adMu[j] / ptCluster->adBeta[k];
				adMu[j] /= adPi[k];
			}

			ptCluster->adNu[k] = ptVBParams->dNu0 + adPi[k];

			/*calculate covariance matrices */
			for (i = 0; i < nN; i++)
			{
				if (aadZ[i][k] > MIN_Z)
				{
					double adDiff[nD];

					for (j = 0; j < nD; j++)
					{
						adDiff[j] =
							aadX[i][j] - adMu[j];
					}

					for (l = 0; l < nD; l++)
					{
						for (m = 0; m <= l; m++)
						{
							aadCovar[l][m] +=
								aadZ[i][k] *
								adDiff[l] *
								adDiff[m];
						}
					}
				}
			}

			for (l = 0; l < nD; l++)
			{
				for (m = l + 1; m < nD; m++)
				{
					aadCovar[l][m] = aadCovar[m][l];
				}
			}

			/*save sample covariances for later use */
			for (l = 0; l < nD; l++)
			{
				for (m = 0; m < nD; m++)
				{
					double dC = aadCovar[l][m] / adPi[k];
					gsl_matrix_set(ptCluster->aptCovar[k],
								   l, m, dC);
				}
			}

			/*Now perform equation 10.62 */
			dF = (ptVBParams->dBeta0 * adPi[k]) /
				 ptCluster->adBeta[k];
			for (l = 0; l < nD; l++)
			{
				for (m = 0; m <= l; m++)
				{
					aadInvWK[l][m] =
						gsl_matrix_get(ptVBParams->ptInvW0,
									   l,
									   m) +
						aadCovar[l][m] +
						dF * adMu[l] * adMu[m];
				}
			}

			for (l = 0; l < nD; l++)
			{
				for (m = 0; m <= l; m++)
				{
					aadCovar[l][m] /= adPi[k];
					gsl_matrix_set(ptSigmaMatrix, l, m,
								   aadInvWK[l][m]);
					gsl_matrix_set(ptSigmaMatrix, m, l,
								   aadInvWK[l][m]);
				}
			}

			/*Implement Equation 10.65 */
			adLDet[k] = ((double)nD) * log(2.0);

			for (l = 0; l < nD; l++)
			{
				double dX =
					0.5 * (ptCluster->adNu[k] - (double)l);
				adLDet[k] += gsl_sf_psi(dX);
			}

			adLDet[k] -= decomposeMatrix(ptSigmaMatrix, nD);
		}
		else
		{
			/*Equation 10.60 */
			adPi[k] = 0.0;

			ptCluster->adBeta[k] = ptVBParams->dBeta0;

			for (j = 0; j < nD; j++)
			{
				/*Equation 10.61 */
				ptCluster->aadM[k][j] = 0.0;
				adMu[j] = 0.0;
			}

			ptCluster->adNu[k] = ptVBParams->dNu0;

			for (l = 0; l < nD; l++)
			{
				for (m = 0; m <= l; m++)
				{
					aadInvWK[l][m] =
						gsl_matrix_get(ptVBParams->ptInvW0,
									   l, m);
				}
			}

			for (l = 0; l < nD; l++)
			{
				for (m = 0; m <= l; m++)
				{
					aadInvWK[l][m] =
						gsl_matrix_get(ptVBParams->ptInvW0,
									   l, m);
				}
			}

			for (l = 0; l < nD; l++)
			{
				for (m = 0; m <= l; m++)
				{
					gsl_matrix_set(ptSigmaMatrix, l, m,
								   aadInvWK[l][m]);
					gsl_matrix_set(ptSigmaMatrix, m, l,
								   aadInvWK[l][m]);
				}
			}

			/*Implement Equation 10.65 */
			adLDet[k] = ((double)nD) * log(2.0);

			for (l = 0; l < nD; l++)
			{
				double dX =
					0.5 * (ptCluster->adNu[k] - (double)l);
				adLDet[k] += gsl_sf_psi(dX);
			}

			adLDet[k] -= decomposeMatrix(ptSigmaMatrix, nD);
		}
	}

	/*Normalise pi */
	if (1)
	{
		double dNP = 0.0;

		for (k = 0; k < nK; k++)
		{
			dNP += adPi[k];
		}

		for (k = 0; k < nK; k++)
		{
			adPi[k] /= dNP;
		}
	}

	/*free up memory */
	for (i = 0; i < nD; i++)
	{
		free(aadCovar[i]);
		free(aadInvWK[i]);
	}

	free(aadCovar);
	free(aadInvWK);

	return;

memoryError:
	fprintf(stderr, "Failed allocating memory in performMStep\n");
	fflush(stderr);
	exit(EXIT_FAILURE);
}

void calcCovarMatrices(t_Cluster *ptCluster, t_Data *ptData)
{
	int i = 0, j = 0, k = 0, l = 0, m = 0;
	int nN = ptData->nN, nK = ptCluster->nK, nD = ptData->nD;
	double **aadZ = ptCluster->aadZ, **aadX = ptData->aadX;
	double *adPi = ptCluster->adPi, **aadCovar = NULL;
	double dN = (double)nN;

	aadCovar = (double **)malloc(nD * sizeof(double *));
	if (!aadCovar)
		goto memoryError;

	for (i = 0; i < nD; i++)
	{
		aadCovar[i] = (double *)malloc(nD * sizeof(double));
		if (!aadCovar[i])
			goto memoryError;
	}

	for (k = 0; k < nK; k++)
	{ /*loop components */
		double *adMu = ptCluster->aadMu[k];
		gsl_matrix *ptSigmaMatrix = ptCluster->aptSigma[k];
		/*recompute mixture weights and means */
		for (j = 0; j < nD; j++)
		{
			adMu[j] = 0.0;
			for (l = 0; l < nD; l++)
			{
				aadCovar[j][l] = 0.0;
			}
			/*prevents singularities */
			aadCovar[j][j] = MIN_COVAR;
		}

		/* compute weight associated with component k */
		adPi[k] = 0.0;
		for (i = 0; i < nN; i++)
		{
			adPi[k] += aadZ[i][k];
			for (j = 0; j < nD; j++)
			{
				adMu[j] += aadZ[i][k] * aadX[i][j];
			}
		}
		/*normalise means */
		for (j = 0; j < nD; j++)
		{
			adMu[j] /= adPi[k];
		}

		/*calculate covariance matrices */
		for (i = 0; i < nN; i++)
		{
			double adDiff[nD];

			for (j = 0; j < nD; j++)
			{
				adDiff[j] = aadX[i][j] - adMu[j];
			}

			for (l = 0; l < nD; l++)
			{
				for (m = 0; m <= l; m++)
				{
					aadCovar[l][m] +=
						aadZ[i][k] * adDiff[l] * adDiff[m];
				}
			}
		}

		for (l = 0; l < nD; l++)
		{
			for (m = l + 1; m < nD; m++)
			{
				aadCovar[l][m] = aadCovar[m][l];
			}
		}

		for (l = 0; l < nD; l++)
		{
			for (m = 0; m < nD; m++)
			{
				aadCovar[l][m] /= adPi[k];
				gsl_matrix_set(ptSigmaMatrix, l, m,
							   aadCovar[l][m]);
			}
		}

		adPi[k] /= dN; /*normalise weights */
	}
	/*free up memory */
	for (i = 0; i < nD; i++)
	{
		free(aadCovar[i]);
	}

	//gsl_matrix_free(ptSigmaMatrix);
	free(aadCovar);

	return;
memoryError:
	fprintf(stderr, "Failed allocating memory in performMStep\n");
	fflush(stderr);
	exit(EXIT_FAILURE);
}

void calcCovarMatricesVB(t_Cluster *ptCluster, t_Data *ptData)
{
	int i = 0, j = 0, k = 0, l = 0, m = 0;
	int nN = ptData->nN, nK = ptCluster->nK, nD = ptData->nD;
	double **aadZ = ptCluster->aadZ, **aadX = ptData->aadX;
	double *adLDet = ptCluster->adLDet, *adPi = ptCluster->adPi;
	double **aadCovar = NULL, **aadInvWK = NULL;
	t_VBParams *ptVBParams = ptCluster->ptVBParams;

	aadCovar = (double **)malloc(nD * sizeof(double *));
	if (!aadCovar)
		goto memoryError;

	for (i = 0; i < nD; i++)
	{
		aadCovar[i] = (double *)malloc(nD * sizeof(double));
		if (!aadCovar[i])
			goto memoryError;
	}

	aadInvWK = (double **)malloc(nD * sizeof(double *));
	if (!aadInvWK)
		goto memoryError;

	for (i = 0; i < nD; i++)
	{
		aadInvWK[i] = (double *)malloc(nD * sizeof(double));
		if (!aadInvWK[i])
			goto memoryError;
	}

	/*perform M step */
	for (k = 0; k < nK; k++)
	{ /*loop components */
		double *adMu = ptCluster->aadMu[k];
		gsl_matrix *ptSigmaMatrix = ptCluster->aptSigma[k];
		double dF = 0.0;

		/*recompute mixture weights and means */
		for (j = 0; j < nD; j++)
		{
			adMu[j] = 0.0;
			for (l = 0; l < nD; l++)
			{
				aadCovar[j][l] = 0.0;
			}
		}

		/* compute weight associated with component k */
		adPi[k] = 0.0;
		for (i = 0; i < nN; i++)
		{
			adPi[k] += aadZ[i][k];
			for (j = 0; j < nD; j++)
			{
				adMu[j] += aadZ[i][k] * aadX[i][j];
			}
		}
		/*normalise means */
		if (adPi[k] > MIN_PI)
		{

			/*Equation 10.60 */
			ptCluster->adBeta[k] = ptVBParams->dBeta0 + adPi[k];

			for (j = 0; j < nD; j++)
			{
				/*Equation 10.61 */
				ptCluster->aadM[k][j] =
					adMu[j] / ptCluster->adBeta[k];
				adMu[j] /= adPi[k];
			}

			ptCluster->adNu[k] = ptVBParams->dNu0 + adPi[k];

			/*calculate covariance matrices */
			for (i = 0; i < nN; i++)
			{
				double adDiff[nD];

				for (j = 0; j < nD; j++)
				{
					adDiff[j] = aadX[i][j] - adMu[j];
				}

				for (l = 0; l < nD; l++)
				{
					for (m = 0; m <= l; m++)
					{
						aadCovar[l][m] +=
							aadZ[i][k] * adDiff[l] *
							adDiff[m];
					}
				}
			}

			for (l = 0; l < nD; l++)
			{
				for (m = l + 1; m < nD; m++)
				{
					aadCovar[l][m] = aadCovar[m][l];
				}
			}

			/*save sample covariances for later use */
			for (l = 0; l < nD; l++)
			{
				for (m = 0; m < nD; m++)
				{
					double dC = aadCovar[l][m] / adPi[k];
					gsl_matrix_set(ptCluster->aptCovar[k],
								   l, m, dC);
				}
			}

			/*Now perform equation 10.62 */
			dF = (ptVBParams->dBeta0 * adPi[k]) /
				 ptCluster->adBeta[k];
			for (l = 0; l < nD; l++)
			{
				for (m = 0; m <= l; m++)
				{
					aadInvWK[l][m] =
						gsl_matrix_get(ptVBParams->ptInvW0,
									   l,
									   m) +
						aadCovar[l][m] +
						dF * adMu[l] * adMu[m];
				}
			}

			for (l = 0; l < nD; l++)
			{
				for (m = 0; m <= l; m++)
				{
					//aadCovar[l][m] /= adPi[k];
					gsl_matrix_set(ptSigmaMatrix, l, m,
								   aadInvWK[l][m]);
					gsl_matrix_set(ptSigmaMatrix, m, l,
								   aadInvWK[l][m]);
				}
			}
		}
		else
		{
			/*Equation 10.60 */
			adPi[k] = 0.0;

			ptCluster->adBeta[k] = ptVBParams->dBeta0;

			for (j = 0; j < nD; j++)
			{
				/*Equation 10.61 */
				ptCluster->aadM[k][j] = 0.0;
				adMu[j] = 0.0;
			}

			ptCluster->adNu[k] = ptVBParams->dNu0;

			for (l = 0; l < nD; l++)
			{
				for (m = 0; m <= l; m++)
				{
					aadInvWK[l][m] =
						gsl_matrix_get(ptVBParams->ptInvW0,
									   l, m);
				}
			}

			for (l = 0; l < nD; l++)
			{
				for (m = 0; m <= l; m++)
				{
					//aadCovar[l][m] /= adPi[k];
					gsl_matrix_set(ptSigmaMatrix, l, m,
								   aadInvWK[l][m]);
					gsl_matrix_set(ptSigmaMatrix, m, l,
								   aadInvWK[l][m]);
				}
			}

			/*Implement Equation 10.65 */
			adLDet[k] = ((double)nD) * log(2.0);

			for (l = 0; l < nD; l++)
			{
				double dX =
					0.5 * (ptCluster->adNu[k] - (double)l);
				adLDet[k] += gsl_sf_psi(dX);
			}

			adLDet[k] -= decomposeMatrix(ptSigmaMatrix, nD);
		}
	}

	/*Normalise pi */
	if (1)
	{
		double dNP = 0.0;

		for (k = 0; k < nK; k++)
		{
			dNP += adPi[k];
		}

		for (k = 0; k < nK; k++)
		{
			adPi[k] /= dNP;
		}
	}

	/*free up memory */
	for (i = 0; i < nD; i++)
	{
		free(aadCovar[i]);
		free(aadInvWK[i]);
	}

	//gsl_matrix_free(ptSigmaMatrix);
	free(aadCovar);
	free(aadInvWK);
	return;
memoryError:
	fprintf(stderr, "Failed allocating memory in performMStep\n");
	fflush(stderr);
	exit(EXIT_FAILURE);
}

void updateMeans(t_Cluster *ptCluster, t_Data *ptData)
{
	int i = 0, j = 0, k = 0;
	int nN = ptData->nN, nK = ptCluster->nK, nD = ptData->nD;
	int *anMaxZ = ptCluster->anMaxZ;
	int *anW = ptCluster->anW;
	double **aadX = ptData->aadX, **aadMu = ptCluster->aadMu;

	for (k = 0; k < nK; k++)
	{

		for (j = 0; j < nD; j++)
		{
			aadMu[k][j] = 0.0;
		}
	}

	for (i = 0; i < nN; i++)
	{
		int nZ = anMaxZ[i];

		for (j = 0; j < nD; j++)
		{
			aadMu[nZ][j] += aadX[i][j];
		}
	}

	for (k = 0; k < nK; k++)
	{ /*loop components */

		/*normalise means */
		if (anW[k] > 0)
		{
			for (j = 0; j < nD; j++)
			{
				aadMu[k][j] /= (double)anW[k];
			}
		}
		else
		{
			for (j = 0; j < nD; j++)
			{
				aadMu[k][j] = 0.0;
			}
		}
	}

	return;
}

void initRandom(gsl_rng *ptGSLRNG, t_Cluster *ptCluster, t_Data *ptData)
{
	/*very simple initialisation assign each data point to random cluster */
	int i = 0, k = 0;

	for (i = 0; i < ptData->nN; i++)
	{
		int nIK = -1;

		for (k = 0; k < ptCluster->nK; k++)
		{
			ptCluster->aadZ[i][k] = 0.0;
		}

		nIK = gsl_rng_uniform_int(ptGSLRNG, ptCluster->nK);

		ptCluster->aadZ[i][nIK] = 1.0;
	}

	performMStep(ptCluster, ptData);

	return;
}

void initMustlink(gsl_rng *ptGSLRNG, t_Cluster *ptCluster, t_Data *ptData)
{
	/* initialiation which assign set of individuals which must link to random cluster */
	int i = 0, k = 0;
	char buffer[MAX_LINE_LENGTH];
	FILE *f;

	for (i = 0; i < ptData->nN; i++)
		for (k = 0; k < ptCluster->nK; k++)
			ptCluster->aadZ[i][k] = 0.0;

	f = fopen(szMLInputFile, "r");
	if (f == NULL)
		return;

	while (fgets(buffer, MAX_LINE_LENGTH, f))
	{
		int nIK = gsl_rng_uniform_int(ptGSLRNG, ptCluster->nK);
		char *tok = buffer, *pos;
		while ((pos = strchr(tok, ',')))
		{
			*pos = 0;
			i = atoi(tok);
			tok = pos + 1;

			ptCluster->aadZ[i][nIK] = 1.0;
		}
		i = atoi(tok);
		ptCluster->aadZ[i][nIK] = 1.0;
	}

	fclose(f);

	performMStep(ptCluster, ptData);
}

double calcDist(double *adX, double *adMu, int nD)
{
	double dDist = 0.0;
	int i = 0;

	for (i = 0; i < nD; i++)
	{
		double dV = adX[i] - adMu[i];
		dDist += dV * dV;
	}

	return sqrt(dDist);
}

void initKMeans(gsl_rng *ptGSLRNG, t_Cluster *ptCluster, t_Data *ptData)
{
	/*very simple initialisation assign each data point to random cluster */
	int i = 0, k = 0, nN = ptData->nN, nK = ptCluster->nK, nD = ptData->nD;
	double **aadMu = ptCluster->aadMu, **aadX = ptData->aadX;
	int *anMaxZ = ptCluster->anMaxZ, *anW = ptCluster->anW, nChange = nN;
	int nIter = 0, nMaxIter = ptCluster->nMaxIter;
	for (i = 0; i < nN; i++)
	{
		int nIK = gsl_rng_uniform_int(ptGSLRNG, nK);

		ptCluster->anMaxZ[i] = nIK;
		anW[nIK]++;
	}

	updateMeans(ptCluster, ptData);

	while (nChange > 0 && nIter < nMaxIter)
	{
		nChange = 0;
		/*reassign vectors */
		for (i = 0; i < nN; i++)
		{
			double dMinDist = DBL_MAX;
			int nMinK = NOT_SET;

			for (k = 0; k < nK; k++)
			{
				double dDist = calcDist(aadX[i], aadMu[k], nD);

				if (dDist < dMinDist)
				{
					nMinK = k;
					dMinDist = dDist;
				}
			}

			if (nMinK != anMaxZ[i])
			{
				int nCurr = anMaxZ[i];
				nChange++;
				anW[nCurr]--;
				anW[nMinK]++;
				anMaxZ[i] = nMinK;

				/*check for empty clusters */
				if (anW[nCurr] == 0)
				{
					int nRandI =
						gsl_rng_uniform_int(ptGSLRNG, nN);
					int nKI = 0;
					/*select at random from non empty clusters */

					while (anW[anMaxZ[nRandI]] == 1)
					{
						nRandI =
							gsl_rng_uniform_int(ptGSLRNG, nN);
					}

					nKI = anMaxZ[nRandI];
					anW[nKI]--;
					anW[nCurr] = 1;
					anMaxZ[nRandI] = nCurr;
				}
			}
		}
		//printf("%d %d\n",nIter,nChange);
		nIter++;
		updateMeans(ptCluster, ptData);
	}

	for (i = 0; i < nN; i++)
	{
		for (k = 0; k < nK; k++)
		{
			ptCluster->aadZ[i][k] = 0.0;
		}
		ptCluster->aadZ[i][anMaxZ[i]] = 1.0;
	}

	performMStep(ptCluster, ptData);
	return;
}

/*note assuming you are using inverse W matrix*/
double dLogWishartB(gsl_matrix *ptInvW, int nD, double dNu, int bInv)
{
	int i = 0;
	double dRet = 0.0, dT = 0.0;
	double dLogDet = 0.0, dD = (double)nD;
	gsl_matrix *ptTemp = gsl_matrix_alloc(nD, nD);

	gsl_matrix_memcpy(ptTemp, ptInvW);

	dLogDet = decomposeMatrix(ptTemp, nD);

	if (bInv == TRUE)
	{
		dRet = 0.5 * dNu * dLogDet;
	}
	else
	{
		dRet = -0.5 * dNu * dLogDet;
	}

	dT = 0.5 * dNu * dD * log(2.0);

	dT += 0.25 * dD * (dD - 1.0) * log(M_PI);

	for (i = 0; i < nD; i++)
	{
		dT += gsl_sf_lngamma(0.5 * (dNu - (double)i));
	}

	gsl_matrix_free(ptTemp);

	return dRet - dT;
}

double dWishartExpectLogDet(gsl_matrix *ptW, double dNu, int nD)
{
	int i = 0;
	double dRet = 0.0, dLogDet = 0.0, dD = (double)nD;
	gsl_matrix *ptTemp = gsl_matrix_alloc(nD, nD);

	gsl_matrix_memcpy(ptTemp, ptW);

	dLogDet = decomposeMatrix(ptW, nD);

	dRet = dD * log(2.0) + dLogDet;

	for (i = 0; i < nD; i++)
	{
		dRet += gsl_sf_psi(0.5 * (dNu - (double)i));
	}

	gsl_matrix_free(ptTemp);

	return dRet;
}

double dWishartEntropy(gsl_matrix *ptW, double dNu, int nD)
{
	double dRet = -dLogWishartB(ptW, nD, dNu, FALSE);
	double dExpLogDet = dWishartExpectLogDet(ptW, dNu, nD);
	double dD = (double)nD;

	dRet -= 0.5 * (dNu - dD - 1.0) * dExpLogDet;

	dRet += 0.5 * dNu * dD;

	return dRet;
}

double calcVBL(t_Cluster *ptCluster)
{
	int i = 0, k = 0, l = 0, nN = ptCluster->nN;
	int nK = ptCluster->nK, nD = ptCluster->nD;
	double dBishop1 = 0.0, dBishop2 = 0.0, dBishop3 = 0.0, dBishop4 = 0.0, dBishop5 = 0.0; /*Bishop equations 10.71... */
	gsl_matrix *ptRes = gsl_matrix_alloc(nD, nD);
	gsl_vector *ptDiff = gsl_vector_alloc(nD);
	gsl_vector *ptR = gsl_vector_alloc(nD);
	double dD = (double)nD;
	double **aadMu = ptCluster->aadMu, **aadM = ptCluster->aadM, **aadZ = ptCluster->aadZ;
	double *adBeta = ptCluster->adBeta, *adNu = ptCluster->adNu, *adLDet = ptCluster->adLDet, *adPi = ptCluster->adPi;
	double adNK[nK];
	double d2Pi = 2.0 * M_PI, dBeta0 = ptCluster->ptVBParams->dBeta0, dNu0 = ptCluster->ptVBParams->dNu0, dRet = 0.0;
	double dK = 0.0;

	for (k = 0; k < nK; k++)
	{
		adNK[k] = 0.0;
	}

	/*Equation 10.72 */
	for (i = 0; i < nN; i++)
	{
		for (k = 0; k < nK; k++)
		{
			adNK[k] += aadZ[i][k];
			if (adPi[k] > 0.0)
			{
				dBishop2 += aadZ[i][k] * log(adPi[k]);
			}
		}
	}

	for (k = 0; k < nK; k++)
	{
		if (adNK[k] > 0.0)
		{
			dK++;
		}
	}

	/*Equation 10.71 */
	for (k = 0; k < nK; k++)
	{
		if (adNK[k] > 0.0)
		{
			double dT1 = 0.0, dT2 = 0.0, dF = 0.0;

			gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0,
						   ptCluster->aptCovar[k],
						   ptCluster->aptSigma[k], 0.0, ptRes);

			for (l = 0; l < nD; l++)
			{
				dT1 += gsl_matrix_get(ptRes, l, l);
			}

			for (l = 0; l < nD; l++)
			{
				gsl_vector_set(ptDiff, l,
							   aadMu[k][l] - aadM[k][l]);
			}

			gsl_blas_dsymv(CblasLower, 1.0, ptCluster->aptSigma[k],
						   ptDiff, 0.0, ptR);

			gsl_blas_ddot(ptDiff, ptR, &dT2);

			dF = adLDet[k] - adNu[k] * (dT1 + dT2) -
				 dD * (log(d2Pi) + (1.0 / adBeta[k]));

			dBishop1 += 0.5 * adNK[k] * dF;
		}
	}

	/*Equation 10.74 */
	for (k = 0; k < nK; k++)
	{
		if (adNK[k] > 0.0)
		{
			double dT1 = 0.0, dT2 = 0.0, dF = 0.0;

			gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0,
						   ptCluster->ptVBParams->ptInvW0,
						   ptCluster->aptSigma[k], 0.0, ptRes);

			for (l = 0; l < nD; l++)
			{
				dT1 += gsl_matrix_get(ptRes, l, l);
			}

			for (l = 0; l < nD; l++)
			{
				gsl_vector_set(ptDiff, l, aadM[k][l]);
			}

			gsl_blas_dsymv(CblasLower, 1.0, ptCluster->aptSigma[k],
						   ptDiff, 0.0, ptR);

			gsl_blas_ddot(ptDiff, ptR, &dT2);

			dF = dD * log(dBeta0 / d2Pi) + adLDet[k] -
				 ((dD * dBeta0) / adBeta[k]) -
				 dBeta0 * adNu[k] * dT2 - adNu[k] * dT1;

			dBishop3 += 0.5 * (dF + (dNu0 - dD - 1.0) * adLDet[k]);
		}
	}

	dBishop3 += dK * ptCluster->ptVBParams->dLogWishartB;

	/*Equation 10.75 */
	for (i = 0; i < nN; i++)
	{
		for (k = 0; k < nK; k++)
		{
			if (aadZ[i][k] > 0.0)
			{
				dBishop4 += aadZ[i][k] * log(aadZ[i][k]);
			}
		}
	}

	/*Equation 10.77 */
	for (k = 0; k < nK; k++)
	{
		if (adNK[k] > 0.0)
		{
			dBishop5 +=
				0.5 * adLDet[k] + 0.5 * dD * log(adBeta[k] / d2Pi) -
				0.5 * dD -
				dWishartExpectLogDet(ptCluster->aptSigma[k],
									 adNu[k], nD);
		}
	}

	gsl_matrix_free(ptRes);
	gsl_vector_free(ptDiff);
	gsl_vector_free(ptR);

	dRet = dBishop1 + dBishop2 + dBishop3 - dBishop4 - dBishop5;

	return dRet;
}

void calcZ(t_Cluster *ptCluster, t_Data *ptData)
{
	double **aadX = ptData->aadX, **aadZ = ptCluster->aadZ;
	int i = 0, k = 0, l = 0;
	int nK = ptCluster->nK, nD = ptCluster->nD, nN = ptData->nN;
	gsl_vector *ptDiff = gsl_vector_alloc(nD);
	gsl_vector *ptRes = gsl_vector_alloc(nD);
	double adDist[nK], dD = (double)nD;
	double **aadM = ptCluster->aadM, *adPi = ptCluster->adPi;

	for (i = 0; i < nN; i++)
	{
		double dMinDist = DBL_MAX;
		double dTotalZ = 0.0;
		double dNTotalZ = 0.0;

		for (k = 0; k < nK; k++)
		{
			if (adPi[k] > 0.)
			{
				/*set vector to data point */
				for (l = 0; l < nD; l++)
				{
					gsl_vector_set(ptDiff, l,
								   aadX[i][l] - aadM[k][l]);
				}

				gsl_blas_dsymv(CblasLower, 1.0,
							   ptCluster->aptSigma[k], ptDiff,
							   0.0, ptRes);

				gsl_blas_ddot(ptDiff, ptRes, &adDist[k]);

				adDist[k] *= ptCluster->adNu[k];

				adDist[k] -= ptCluster->adLDet[k];

				adDist[k] += dD / ptCluster->adBeta[k];

				if (adDist[k] < dMinDist)
				{
					dMinDist = adDist[k];
				}
			}
		}

		for (k = 0; k < nK; k++)
		{
			if (adPi[k] > 0.)
			{
				aadZ[i][k] =
					adPi[k] * exp(-0.5 *
								  (adDist[k] - dMinDist));
				dTotalZ += aadZ[i][k];
			}
			else
			{
				aadZ[i][k] = 0.0;
			}
		}

		for (k = 0; k < nK; k++)
		{
			double dF = aadZ[i][k] / dTotalZ;
			if (dF < MIN_Z)
			{
				aadZ[i][k] = 0.0;
			}
			dNTotalZ += aadZ[i][k];
		}
		if (dNTotalZ > 0.)
		{
			for (k = 0; k < nK; k++)
			{
				aadZ[i][k] /= dNTotalZ;
			}
		}
	}

	gsl_vector_free(ptRes);
	gsl_vector_free(ptDiff);
	return;
}

void gmmTrainVB(t_Cluster *ptCluster, t_Data *ptData)
{
	int i = 0, k = 0, nIter = 0;
	int nN = ptData->nN, nK = ptCluster->nK;
	/*change in log-likelihood */
	double dLastVBL = 0.0, dDelta = DBL_MAX;
	double **aadZ = ptCluster->aadZ;
	int nMaxIter = ptCluster->nMaxIter;
	double dEpsilon = ptCluster->dEpsilon;

	/*calculate data likelihood */
	calcZ(ptCluster, ptData);
	ptCluster->dVBL = calcVBL(ptCluster);

	while (nIter < nMaxIter && dDelta > dEpsilon)
	{

		/*update parameter estimates */
		performMStep(ptCluster, ptData);

		/*calculate responsibilities */
		calcZ(ptCluster, ptData);

		dLastVBL = ptCluster->dVBL;
		ptCluster->dVBL = calcVBL(ptCluster);
		dDelta = fabs(ptCluster->dVBL - dLastVBL);

		printf("[DEBUG-%d] nIter=%d ; maxIter=%d ; dVBL=%.4f ; delta=%.4e ; epsilon=%.4e\n",
			   ptCluster->nThread, nIter, nMaxIter, ptCluster->dVBL, dDelta, dEpsilon);
		nIter++;
	}

	/*assign to best clusters */
	for (i = 0; i < nN; i++)
	{
		double dMaxZ = aadZ[i][0];
		int nMaxK = 0;
		for (k = 1; k < nK; k++)
		{
			if (aadZ[i][k] > dMaxZ)
			{
				nMaxK = k;
				dMaxZ = aadZ[i][k];
			}
		}
		ptCluster->anMaxZ[i] = nMaxK;
	}

	return;
}

void writeClusters(char *szOutFile, t_Cluster *ptCluster, t_Data *ptData)
{
	int nN = ptCluster->nN, i = 0;
	FILE *ofp = fopen(szOutFile, "w");

	if (ofp)
	{
		for (i = 0; i < nN; i++)
		{
			fprintf(ofp, "%s,%d\n", ptData->aszSampleNames[i],
					ptCluster->anMaxZ[i]);
		}
		fclose(ofp);
		printf("[DEBUG] file '%s' closed\n", szOutFile);
	}
	else
	{
		fprintf(stderr,
				"Failed to open %s for writing in writeClusters\n",
				szOutFile);
		fflush(stderr);
	}
}

void writeMeans(char *szOutFile, t_Cluster *ptCluster)
{
	int nK = ptCluster->nK, nD = ptCluster->nD, i = 0, j = 0;
	FILE *ofp = fopen(szOutFile, "w");

	if (ofp)
	{
		for (i = 0; i < nK; i++)
		{
			for (j = 0; j < nD - 1; j++)
			{
				fprintf(ofp, "%f,", ptCluster->aadMu[i][j]);
			}
			fprintf(ofp, "%f\n", ptCluster->aadMu[i][nD - 1]);
		}
		fclose(ofp);
	}
	else
	{
		fprintf(stderr,
				"Failed to open %s for writing in writeMeanss\n",
				szOutFile);
		fflush(stderr);
	}
}

void writeTMeans(char *szOutFile, t_Cluster *ptCluster, t_Data *ptData)
{
	int nK = ptCluster->nK, nD = ptCluster->nD, nT = ptData->nT, i = 0, j = 0;
	FILE *ofp = fopen(szOutFile, "w");
	gsl_vector *ptVector = gsl_vector_alloc(nD);
	gsl_vector *ptTVector = gsl_vector_alloc(nT);

	if (ofp)
	{
		for (i = 0; i < nK; i++)
		{
			for (j = 0; j < nD; j++)
			{
				gsl_vector_set(ptVector, j,
							   ptCluster->aadMu[i][j]);
			}

			gsl_blas_dgemv(CblasNoTrans, 1.0, ptData->ptTMatrix,
						   ptVector, 0.0, ptTVector);

			for (j = 0; j < nT - 1; j++)
			{
				fprintf(ofp, "%f,",
						gsl_vector_get(ptTVector, j));
			}
			fprintf(ofp, "%f\n", gsl_vector_get(ptTVector, nD - 1));
		}
		fclose(ofp);
	}
	else
	{
		fprintf(stderr,
				"Failed to open %s for writing in writeMeanss\n",
				szOutFile);
		fflush(stderr);
	}

	gsl_vector_free(ptVector);
	gsl_vector_free(ptTVector);
}

void *fitEM(void *pvCluster)
{
	t_Cluster *ptCluster = (t_Cluster *)pvCluster;
	gsl_rng *ptGSLRNG = NULL;
	const gsl_rng_type *ptGSLRNGType = NULL;
	time_t start_t, end_t;

	time(&start_t);
	/*initialise GSL RNG*/
	ptGSLRNGType = gsl_rng_default;
	ptGSLRNG = gsl_rng_alloc(ptGSLRNGType);

	gsl_rng_set(ptGSLRNG, ptCluster->lSeed);

	if (ptCluster->szCOutFile == NULL)
	{
		printf("[DEBUG-%d] init KMeans...\n", ptCluster->nThread);
		initKMeans(ptGSLRNG, ptCluster, ptCluster->ptData);
	}
	else
	{
		printf("[DEBUG-%d] init using must-link...\n", ptCluster->nThread);
		initMustlink(ptGSLRNG, ptCluster, ptCluster->ptData);
	}

	printf("[DEBUG-%d] train VBGMM...\n", ptCluster->nThread);
	gmmTrainVB(ptCluster, ptCluster->ptData);

	gsl_rng_free(ptGSLRNG);
	time(&end_t);
	printf("[DEBUG-%d] Lower bound=%.2f ; Execution time fitEM=%.2f sec\n",
		   ptCluster->nThread, ptCluster->dVBL, difftime(end_t, start_t));

	return NULL;
}

void *runRThreads(void *pvpDCluster)
{
	t_Cluster **pptDCluster = (t_Cluster **)pvpDCluster;
	t_Cluster *ptDCluster = (t_Cluster *)*pptDCluster;
	double dBestVBL = -DBL_MAX;
	t_Cluster **aptCluster = NULL;
	pthread_t atRestarts[N_RTHREADS]; /*run each restart on a separate thread */
	int r = 0, nBestR = -1;
	char *szCOutFile = NULL;
	aptCluster = (t_Cluster **)malloc(N_RTHREADS * sizeof(t_Cluster *));
	if (!aptCluster)
		goto memoryError;

	for (r = 0; r < N_RTHREADS; r++)
	{
		if (ptDCluster->szCOutFile != NULL)
		{
			szCOutFile =
				(char *)malloc(sizeof(char) * MAX_LINE_LENGTH);
			sprintf(szCOutFile, "%sr%d.csv", ptDCluster->szCOutFile,
					r);
		}

		aptCluster[r] = (t_Cluster *)malloc(sizeof(t_Cluster));

		allocateCluster(aptCluster[r], ptDCluster->nN, ptDCluster->nK,
						ptDCluster->nD, ptDCluster->ptData,
						ptDCluster->lSeed + r * R_PRIME,
						ptDCluster->nMaxIter, ptDCluster->dEpsilon,
						szCOutFile);
		aptCluster[r]->ptVBParams = ptDCluster->ptVBParams;
		aptCluster[r]->nThread = r;
		pthread_create(&atRestarts[r], NULL, fitEM, (void *)aptCluster[r]);
	}
	//TODO: use pthread_tryjoin_np to join pthreads or display some info
	for (r = 0; r < N_RTHREADS; r++)
	{
		pthread_join(atRestarts[r], NULL);
	}

	/*free up memory associated with input cluster */
	free(ptDCluster);

	for (r = 0; r < N_RTHREADS; r++)
	{
		if (aptCluster[r]->dVBL > dBestVBL)
		{
			nBestR = r;
			dBestVBL = aptCluster[r]->dVBL;
		}
	}
	printf("[DEBUG] best clustering result dVBL = %.4f\n", dBestVBL);

	*pptDCluster = aptCluster[nBestR];
	for (r = 0; r < N_RTHREADS; r++)
	{
		if (r != nBestR)
		{
			destroyCluster(aptCluster[r]);
			free(aptCluster[r]);
		}
	}
	free(aptCluster);

	return NULL;

memoryError:
	fprintf(stderr, "Failed allocating memory in runRThreads\n");
	fflush(stderr);
	exit(EXIT_FAILURE);
}

void calcSampleVar(t_Data *ptData, double *adVar, double *adMu)
{
	double **aadX = ptData->aadX;
	int i = 0, n = 0;
	int nD = ptData->nD, nN = ptData->nN;
	/*sample means */
	double dN = (double)nN;

	for (i = 0; i < nD; i++)
	{
		adMu[i] = 0.0;
		adVar[i] = 0.0;
	}

	for (i = 0; i < nD; i++)
	{
		for (n = 0; n < nN; n++)
		{
			adMu[i] += aadX[n][i];
			adVar[i] += aadX[n][i] * aadX[n][i];
		}

		adMu[i] /= dN;

		adVar[i] = (adVar[i] - dN * adMu[i] * adMu[i]) / (dN - 1.0);
	}

	return;
}

void compressCluster(t_Cluster *ptCluster)
{
	int i = 0, k = 0, nNewK = 0, nN = ptCluster->nN;
	double **aadNewZ = NULL, dN = (double)nN;

	for (i = 0; i < ptCluster->nK; i++)
	{
		if (ptCluster->adPi[i] > 0.0)
		{
			nNewK++;
		}
	}

	aadNewZ = (double **)malloc(nN * sizeof(double *));
	if (!aadNewZ)
		goto memoryError;

	for (i = 0; i < nN; i++)
	{
		aadNewZ[i] = (double *)malloc(nNewK * sizeof(double));
		if (!aadNewZ[i])
			goto memoryError;
	}

	for (i = 0; i < nN; i++)
	{
		int nC = 0;
		for (k = 0; k < ptCluster->nK; k++)
		{
			if (ptCluster->adPi[k] > 0.0)
			{
				aadNewZ[i][nC] = ptCluster->aadZ[i][k];
				nC++;
			}
		}
	}

	for (i = 0; i < nN; i++)
	{
		free(ptCluster->aadZ[i]);
	}
	free(ptCluster->aadZ);

	/*reset Z and K */
	ptCluster->aadZ = aadNewZ;
	ptCluster->nK = nNewK;

	/*recalculate Pi */
	for (k = 0; k < ptCluster->nK; k++)
	{
		ptCluster->adPi[k] = 0.0;
		for (i = 0; i < nN; i++)
		{
			ptCluster->adPi[k] += ptCluster->aadZ[i][k];
		}
		ptCluster->adPi[k] /= dN;
	}

	/*assign to best clusters */
	for (i = 0; i < nN; i++)
	{
		double dMaxZ = ptCluster->aadZ[i][0];
		int nMaxK = 0;
		for (k = 1; k < ptCluster->nK; k++)
		{
			if (ptCluster->aadZ[i][k] > dMaxZ)
			{
				nMaxK = k;
				dMaxZ = ptCluster->aadZ[i][k];
			}
		}
		ptCluster->anMaxZ[i] = nMaxK;
	}

	for (k = 0; k < ptCluster->nK; k++)
	{
		ptCluster->anW[k] = 0;
	}

	for (i = 0; i < nN; i++)
	{
		ptCluster->anW[ptCluster->anMaxZ[i]]++;
	}

	return;

memoryError:
	fprintf(stderr, "Failed allocating memory in main\n");
	fflush(stderr);
	exit(EXIT_FAILURE);
}

max_stars_repo_name: keuv-grvl/fennec

max_issues_repo_name: keuv-grvl/fennec

max_forks_repo_name: keuv-grvl/fennec


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

/* Author:  G. Jungman */

#ifndef __GSL_SF_LAMBERT_H__
#define __GSL_SF_LAMBERT_H__

#include <gsl/gsl_sf_result.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

/* Lambert's Function W_0(x)
 *
 * W_0(x) is the principal branch of the
 * implicit function defined by W e^W = x.
 *
 * -1/E < x < \infty
 *
 * exceptions: GSL_EMAXITER;
 */
int     gsl_sf_lambert_W0_e(double x, gsl_sf_result * result);
double  gsl_sf_lambert_W0(double x);

/* Lambert's Function W_{-1}(x)
 *
 * W_{-1}(x) is the second real branch of the
 * implicit function defined by W e^W = x.
 * It agrees with W_0(x) when x >= 0.
 *
 * -1/E < x < \infty
 *
 * exceptions: GSL_MAXITER;
 */
int     gsl_sf_lambert_Wm1_e(double x, gsl_sf_result * result);
double  gsl_sf_lambert_Wm1(double x);

__END_DECLS

#endif /* __GSL_SF_LAMBERT_H__ */


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

#include <config.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spmatrix.h>

/*
gsl_spmatrix_ccs()
  Create a sparse matrix in compressed column format

Inputs: T - sparse matrix in triplet format

Return: pointer to new matrix (should be freed when finished with it)
*/

gsl_spmatrix *
gsl_spmatrix_ccs(const gsl_spmatrix *T)
{
  if (!GSL_SPMATRIX_ISTRIPLET(T))
    {
      GSL_ERROR_NULL("matrix must be in triplet format", GSL_EINVAL);
    }
  else
    {
      const size_t *Tj; /* column indices of triplet matrix */
      size_t *Cp;       /* column pointers of compressed column matrix */
      size_t *w;        /* copy of column pointers */
      gsl_spmatrix *m;
      size_t n;

      m = gsl_spmatrix_alloc_nzmax(T->size1, T->size2, T->nz,
                                   GSL_SPMATRIX_CCS);
      if (!m)
        return NULL;

      Tj = T->p;
      Cp = m->p;

      /* initialize column pointers to 0 */
      for (n = 0; n < m->size2 + 1; ++n)
        Cp[n] = 0;

      /*
       * compute the number of elements in each column:
       * Cp[j] = # non-zero elements in column j
       */
      for (n = 0; n < T->nz; ++n)
        Cp[Tj[n]]++;

      /* compute column pointers: p[j] = p[j-1] + nnz[j-1] */
      gsl_spmatrix_cumsum(m->size2, Cp);

      /* make a copy of the column pointers */
      w = (size_t *) m->work;
      for (n = 0; n < m->size2; ++n)
        w[n] = Cp[n];

      /* transfer data from triplet format to CCS */
      for (n = 0; n < T->nz; ++n)
        {
          size_t k = w[Tj[n]]++;
          m->i[k] = T->i[n];
          m->data[k] = T->data[n];
        }

      m->nz = T->nz;

      return m;
    }
}

gsl_spmatrix *
gsl_spmatrix_compcol(const gsl_spmatrix *T)
{
  return gsl_spmatrix_ccs(T);
}

/*
gsl_spmatrix_crs()
  Create a sparse matrix in compressed row format

Inputs: T - sparse matrix in triplet format

Return: pointer to new matrix (should be freed when finished with it)
*/

gsl_spmatrix *
gsl_spmatrix_crs(const gsl_spmatrix *T)
{
  if (!GSL_SPMATRIX_ISTRIPLET(T))
    {
      GSL_ERROR_NULL("matrix must be in triplet format", GSL_EINVAL);
    }
  else
    {
      const size_t *Ti; /* row indices of triplet matrix */
      size_t *Cp;       /* row pointers of compressed row matrix */
      size_t *w;        /* copy of column pointers */
      gsl_spmatrix *m;
      size_t n;

      m = gsl_spmatrix_alloc_nzmax(T->size1, T->size2, T->nz,
                                   GSL_SPMATRIX_CRS);
      if (!m)
        return NULL;

      Ti = T->i;
      Cp = m->p;

      /* initialize row pointers to 0 */
      for (n = 0; n < m->size1 + 1; ++n)
        Cp[n] = 0;

      /*
       * compute the number of elements in each row:
       * Cp[i] = # non-zero elements in row i
       */
      for (n = 0; n < T->nz; ++n)
        Cp[Ti[n]]++;

      /* compute row pointers: p[i] = p[i-1] + nnz[i-1] */
      gsl_spmatrix_cumsum(m->size1, Cp);

      /* make a copy of the row pointers */
      w = (size_t *) m->work;
      for (n = 0; n < m->size1; ++n)
        w[n] = Cp[n];

      /* transfer data from triplet format to CRS */
      for (n = 0; n < T->nz; ++n)
        {
          size_t k = w[Ti[n]]++;
          m->i[k] = T->p[n];
          m->data[k] = T->data[n];
        }

      m->nz = T->nz;

      return m;
    }
}

/*
gsl_spmatrix_cumsum()

Compute the cumulative sum:

p[j] = Sum_{k=0...j-1} c[k]

0 <= j < n + 1

Alternatively,
p[0] = 0
p[j] = p[j - 1] + c[j - 1]

Inputs: n - length of input array
        c - (input/output) array of size n + 1
            on input, contains the n values c[k]
            on output, contains the n + 1 values p[j]

Return: success or error
*/

void
gsl_spmatrix_cumsum(const size_t n, size_t *c)
{
  size_t sum = 0;
  size_t k;

  for (k = 0; k < n; ++k)
    {
      size_t ck = c[k];
      c[k] = sum;
      sum += ck;
    }

  c[n] = sum;
} /* gsl_spmatrix_cumsum() */