text:
/*System includes*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <sys/stat.h>
#include <float.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_cblas.h>
#include <pthread.h>

#include <Python.h>

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

static PyObject *vbgmm_fit(PyObject *self, PyObject *args)
{
    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;
    sts = driver(szFileStub,nKStart,nLMin,lSeed,nMaxIter,dEpsilon,bCout);
    return Py_BuildValue("i", sts);
}

static PyMethodDef VbgmmMethods[] = {
    {"fit",  vbgmm_fit, METH_VARARGS,
     "Fit a variational Bayesian Gaussian mixture."},
    {NULL, NULL, 0, NULL}        /* Sentinel */
};

/* Python 2.x code */

PyMODINIT_FUNC
initvbgmm(void)
{
    (void) Py_InitModule("vbgmm", VbgmmMethods);
}

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 i = 0, k = 0, nD = 0, nN = 0;
  char szOFile[MAX_LINE_LENGTH];
  FILE *ofp = NULL;
  t_VBParams tVBParams;
  t_Cluster  *ptBestCluster = NULL;
  gsl_matrix *ptTemp = NULL;
  gsl_matrix *ptTVar = NULL;
  /*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 = szFileStub;
  }
  else{
	ptBestCluster->szCOutFile = NULL;
  }
  runRThreads((void *) &ptBestCluster);

  compressCluster(ptBestCluster);

  calcCovarMatrices(ptBestCluster,&tData);

  sprintf(szOFile,"%sclustering_gt%d.csv",tParams.szOutFileStub,tParams.nLMin);
  writeClusters(szOFile,ptBestCluster,&tData);

  sprintf(szOFile,"%spca_means_gt%d.csv",tParams.szOutFileStub,tParams.nLMin);
  writeMeans(szOFile,ptBestCluster);

  sprintf(szOFile,"%smeans_gt%d.csv",tParams.szOutFileStub,tParams.nLMin);
  writeTMeans(szOFile,ptBestCluster,&tData);

  for(k = 0; k < ptBestCluster->nK; k++){
    sprintf(szOFile,"%spca_variances_gt%d_dim%d.csv",tParams.szOutFileStub,tParams.nLMin,k);

    writeSquareMatrix(szOFile, ptBestCluster->aptSigma[k], nD);

    /*not entirely sure this is correct?*/
    gsl_blas_dgemm (CblasNoTrans,CblasNoTrans,1.0,tData.ptTMatrix,ptBestCluster->aptSigma[k],0.0,ptTemp);

    gsl_blas_dgemm (CblasNoTrans,CblasTrans,1.0,ptTemp,tData.ptTMatrix,0.0,ptTVar);

    sprintf(szOFile,"%svariances_gt%d_dim%d.csv",tParams.szOutFileStub,tParams.nLMin,k);

    writeSquareMatrix(szOFile, ptTVar, nD);
  }

  sprintf(szOFile,"%sresponsibilities.csv",tParams.szOutFileStub);

  ofp = fopen(szOFile,"w");
  if(ofp){

    for(i = 0; i < nN; i++){
      for(k = 0; k < ptBestCluster->nK - 1; k++){
	fprintf(ofp,"%f,",ptBestCluster->aadZ[i][k]);
      }
      fprintf(ofp,"%f\n",ptBestCluster->aadZ[i][ptBestCluster->nK - 1]);
    }

    fclose(ofp);
  }
  else{
    fprintf(stderr,"Failed openining %s in main\n", szOFile);
    fflush(stderr);
  }

  sprintf(szOFile,"%svbl.csv",tParams.szOutFileStub);

  ofp = fopen(szOFile,"w");
  if(ofp){
    fprintf(ofp,"%d,%f,%d\n",ptBestCluster->nK,ptBestCluster->dVBL,ptBestCluster->nThread);

    fclose(ofp);
  }
  else{
    fprintf(stderr,"Failed openining %s in main\n", szOFile);
    fflush(stderr);
  }

  /*free up memory in data object*/
  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);

  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->szInputFile)
    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;
	}
      }
    }
  }
  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);
      }
    }
  }
  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;
}

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;
  FILE   *ofp = NULL;

  if(ptCluster->szCOutFile){
	ofp = fopen(ptCluster->szCOutFile,"w");
	if(!ofp){
		fprintf(stderr, "Failed to open file %s in gmmTrainVB\n",ptCluster->szCOutFile);
		fflush(stderr);
	}
  }

  /*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);

    if(ofp){
    	fprintf(ofp,"%d,%f,%f,",nIter, ptCluster->dVBL, dDelta);
    	for(k = 0; k < nK-1; k++){
    		fprintf(ofp,"%f,",ptCluster->adPi[k]);
    	}
    	fprintf(ofp,"%f\n",ptCluster->adPi[nK - 1]);
	fflush(ofp);
    }
    nIter++;
  }

  if(ofp){
	fclose(ofp);
  }

  /*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]);
    }
  }
  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]);
    }
  }
  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));
    }
  }
  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;

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

  gsl_rng_set (ptGSLRNG, ptCluster->lSeed);

  initKMeans(ptGSLRNG, ptCluster, ptCluster->ptData);

  gmmTrainVB(ptCluster, ptCluster->ptData);

  gsl_rng_free(ptGSLRNG);

  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         iret[N_RTHREADS];
  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;
    iret[r] = pthread_create(&atRestarts[r], NULL, fitEM, (void*) aptCluster[r]);
  }

  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;
    }
  }

  *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: binnisb/CONCOCT

max_issues_repo_name: binnisb/CONCOCT

max_forks_repo_name: binnisb/CONCOCT


text:
// The top most header - contains switches and constants
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <gsl/gsl_rng.h>

#define Pi   acos(-1.0)

#define real double

// Cell shape and size
//#define circle
#define square
#define Dia 40  // Cell radius or square side
#define thickness 40

// clamped or hinged MTs

int iprint, iter, t;
char datadir[64] ;

ext: h

max_stars_repo_name: misragaurav/SingleMicrotubuleBuckling

max_issues_repo_name: misragaurav/SingleMicrotubuleBuckling

max_forks_repo_name: misragaurav/SingleMicrotubuleBuckling


text:
/* ode-initval/gear2.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 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 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Gear 2 */

/* Author:  G. Jungman
 */
#include <config.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include "odeiv_util.h"
#include "gsl_odeiv.h"

/* gear2 state object */
typedef struct {
  int primed;                /* flag indicating that yim1 is ready */
  double t_primed;           /* system was primed for this value of t */
  double last_h;             /* last step size */
  gsl_odeiv_step * primer;   /* stepper to use for priming */
  double * yim1;             /* y_{i-1}    */
  double * k;                /* work space */
  double * y0;               /* work space */
  int stutter;
}
gear2_state_t;

static void *
gear2_alloc (size_t dim)
{
  gear2_state_t *state = (gear2_state_t *) malloc (sizeof (gear2_state_t));

  if (state == 0)
    {
      GSL_ERROR_NULL ("failed to allocate space for gear2_state", GSL_ENOMEM);
    }

  state->yim1 = (double *) malloc (dim * sizeof (double));

  if (state->yim1 == 0)
    {
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for yim1", GSL_ENOMEM);
    }

  state->k = (double *) malloc (dim * sizeof (double));

  if (state->k == 0)
    {
      free (state->yim1);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for k", GSL_ENOMEM);
    }

  state->y0 = (double *) malloc (dim * sizeof (double));

  if (state->y0 == 0)
    {
      free (state->k);
      free (state->yim1);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for y0", GSL_ENOMEM);
    }

  state->primed = 0;
  state->primer = gsl_odeiv_step_alloc (gsl_odeiv_step_rk4imp, dim);
  state->last_h = 0.0;

  return state;
}

static int
gear2_apply(void * vstate,
            size_t dim,
            double t,
            double h,
            double y[],
            double yerr[],
            const double dydt_in[],
            double dydt_out[],
            const gsl_odeiv_system * sys)
{
  gear2_state_t *state = (gear2_state_t *) vstate;

  state->stutter = 0;

  if(state->primed == 0 || t == state->t_primed || h != state->last_h) {

    /* Execute a single-step method to prime the process.  Note that
     * we do this if the step size changes, so frequent step size
     * changes will cause the method to stutter.
     *
     * Note that we reuse this method if the time has not changed,
     * which can occur when the adaptive driver is attempting to find
     * an appropriate step-size on its first iteration */

    int status;
    DBL_MEMCPY(state->yim1, y, dim);

    status = gsl_odeiv_step_apply(state->primer, t, h, y, yerr, dydt_in, dydt_out, sys);

    /* Make note of step size and indicate readiness for a Gear step. */

    state->primed = 1;
    state->t_primed = t;
    state->last_h = h;
    state->stutter = 1;

    return status;
  } else {
    /* We have a previous y value in the buffer, and the step
     * sizes match, so we go ahead with the Gear step.
     */

    double * const k  = state->k;
    double * const y0 = state->y0;
    double * const yim1 = state->yim1;

    const int iter_steps = 3;
    int status = 0;
    int nu;
    size_t i;

    DBL_MEMCPY(y0, y, dim);

    /* iterative solution */

    if(dydt_out != NULL) {
      DBL_MEMCPY(k, dydt_out, dim);
    }

    for(nu=0; nu<iter_steps; nu++) {
      int s = GSL_ODEIV_FN_EVAL(sys, t + h, y, k);
      GSL_STATUS_UPDATE(&status, s);
      for(i=0; i<dim; i++) {
        y[i] = ((4.0 * y0[i] - yim1[i]) + 2.0 * h * k[i]) / 3.0;
      }
    }

    /* Estimate error and update the state buffer. */
    for(i=0; i<dim; i++) {
      yerr[i] = h * h * (y[i] - y0[i]); /* error is third order */
      yim1[i] = y0[i];
    }

    /* Make note of step size. */
    state->last_h = h;

    return status;
  }
}

static int
gear2_reset(void * vstate, size_t dim)
{
  gear2_state_t *state = (gear2_state_t *) vstate;

  DBL_ZERO_MEMSET (state->yim1, dim);
  DBL_ZERO_MEMSET (state->k, dim);
  DBL_ZERO_MEMSET (state->y0, dim);

  state->primed = 0;
  state->last_h = 0.0;
  return GSL_SUCCESS;
}

static unsigned int
gear2_order (void *vstate)
{
  gear2_state_t *state = (gear2_state_t *) vstate;
  state = 0; /* prevent warnings about unused parameters */
  return 3;
}

static void
gear2_free(void *vstate)
{
  gear2_state_t *state = (gear2_state_t *) vstate;

  free(state->yim1);
  free(state->k);
  free(state->y0);
  gsl_odeiv_step_free(state->primer);

  free(state);
}

static const gsl_odeiv_step_type gear2_type = { "gear2",	/* name */
  1,				/* can use dydt_in */
  0,				/* gives exact dydt_out */
  &gear2_alloc,
  &gear2_apply,
  &gear2_reset,
  &gear2_order,
  &gear2_free
};

const gsl_odeiv_step_type *gsl_odeiv_step_gear2 = &gear2_type;

ext: c


text:
// Authors: David Blei (blei@cs.princeton.edu)
//          Sean Gerrish (sgerrish@cs.princeton.edu)
//
// Copyright 2011 Sean Gerrish and David Blei
// All Rights Reserved.
//
// See the README for this package for details about modifying or
// distributing this software.

/*
 * state space language model variational inference
 *
 */

#ifndef SSLM_H
#define SSLM_H

#include "gsl-wrappers.h"
#include "params.h"
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <assert.h>
#include <math.h>

#include "data.h"

#define SSLM_MAX_ITER 2 // maximum number of optimization iters
#define SSLM_FIT_THRESHOLD 1e-6 // convergence criterion for fitting sslm
#define INIT_MULT 1000 // multiplier to variance for first obs
// #define OBS_NORM_CUTOFF 10 // norm cutoff after which we use all 0 obs
//#define OBS_NORM_CUTOFF 8 // norm cutoff after which we use all 0 obs
#define OBS_NORM_CUTOFF 2 // norm cutoff after which we use all 0 obs

/*
 * functions for variational inference
 *
 */

// allocate new state space language model variational posterior
sslm_var* sslm_var_alloc(int W, int T);

// allocate extra parameters for inference
void sslm_inference_alloc(sslm_var* var);

// free extra parameters for inference
void sslm_inference_free(sslm_var* var);

// initialize with zero observations
void sslm_zero_init(sslm_var* var,
                    double obs_variance,
                    double chain_variance);

// initialize with counts
void sslm_counts_init(sslm_var* var,
                      double obs_variance,
                      double chain_variance,
                      const gsl_vector* counts);

// initialize from variational observations
void sslm_obs_file_init(sslm_var* var,
                        double obs_variance,
                        double chain_variance,
                        const char* filename);

// compute E[\beta_{w,t}] for t = 1:T
void compute_post_mean(int w, sslm_var* var, double chain_variance);

// compute Var[\beta_{w,t}] for t = 1:T
void compute_post_variance(int w, sslm_var* var, double chain_variance);

// optimize \hat{beta}
void optimize_var_obs(sslm_var* var);

// compute dE[\beta_{w,t}]/d\obs_{w,s} for t = 1:T
void compute_mean_deriv(int word, int time, sslm_var* var,
                        gsl_vector* deriv);

// compute d bound/d obs_{w, t} for t=1:T.
void compute_obs_deriv(int word, gsl_vector* word_counts,
                       gsl_vector* total_counts, sslm_var* var,
                       gsl_matrix* mean_deriv_mtx, gsl_vector* deriv);

// update observations
void update_obs(gsl_matrix* word_counts, gsl_vector* totals,
                sslm_var* var);

// log probability bound
double compute_bound(gsl_matrix* word_counts, gsl_vector* totals,
                     sslm_var* var);

// fit variational distribution
double fit_sslm(sslm_var* var, gsl_matrix* word_counts);

// read and write variational distribution
void write_sslm_var(sslm_var* var, char* out);
sslm_var* read_sslm_var(char* in);

void compute_expected_log_prob(sslm_var* var);
// !!! old function (from doc mixture...)
double expected_log_prob(int w, int t, sslm_var* var);

// update zeta
void update_zeta(sslm_var* var);

#endif

max_stars_repo_name: iwangjian/dtm-lab

max_issues_repo_name: iwangjian/topic-extractor

max_forks_repo_name: iwangjian/topic-extractor


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

/***************************************************************************//**
 *
 * @ingroup CORE_double
 *
 *  CORE_PLASMA_dlag2s converts a double matrix, A, to a
 *  float matrix, B.
 *
 *******************************************************************************
 *
 * @param[in] m
 *          The number of rows of the matrices A and B. m >= 0.
 *
 * @param[in] n
 *          The number of columns of the matrices A and B. n >= 0.
 *
 * @param[in] A
 *          The double m-by-n matrix to convert.
 *
 * @param[in] lda
 *          The leading dimension of the array A. lda >= max(1,m).
 *
 * @param[out] B
 *          The float m-by-n matrix to convert.
 *
 * @param[in] ldb
 *          The leading dimension of the array B. ldb >= max(1,m).
 *
 * @param[out] info
 *         - 0 on successful exit.
 *         - 1 if an entry of the matrix A is greater than the SINGLE
 *         PRECISION overflow threshold, in this case, the content
 *         of B in exit is unspecified.
 *
 ******************************************************************************/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_dlag2s = PCORE_dlag2s
#define CORE_dlag2s PCORE_dlag2s
#endif
void CORE_dlag2s(int m, int n,
                 const double *A, int lda,
                 float *B, int ldb, int *info)
{
    *info = LAPACKE_dlag2s_work(LAPACK_COL_MAJOR, m, n, A, lda, B, ldb);
}

/***************************************************************************//**
 *
 * @ingroup CORE_double
 *
 *  CORE_PLASMA_slag2d converts a float matrix, A, to a
 *  double matrix, B.
 *
 *******************************************************************************
 *
 * @param[in] m
 *          The number of rows of the matrices A and B. m >= 0.
 *
 * @param[in] n
 *          The number of columns of the matrices A and B. n >= 0.
 *
 * @param[in] A
 *          The float m-by-n matrix to convert.
 *
 * @param[in] lda
 *          The leading dimension of the array A. lda >= max(1,m).
 *
 * @param[out] B
 *          The double m-by-n matrix to convert.
 *
 * @param[in] ldb
 *          The leading dimension of the array B. ldb >= max(1,m).
 *
 ******************************************************************************/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_slag2d = PCORE_slag2d
#define CORE_slag2d PCORE_slag2d
#endif
void CORE_slag2d(int m, int n,
                 const float *A, int lda,
                 double *B, int ldb)
{
    LAPACKE_slag2d_work(LAPACK_COL_MAJOR, m, n, A, lda, B, ldb);
}


text:
/* min/gsl_min.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007, 2009 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_MIN_H__
#define __GSL_MIN_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_math.h>

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

__BEGIN_DECLS

typedef struct
  {
    const char *name;
    size_t size;
    int (*set) (void *state, gsl_function * f, double x_minimum, double f_minimum, double x_lower, double f_lower, double x_upper, double f_upper);
    int (*iterate) (void *state, gsl_function * f, double * x_minimum, double * f_minimum, double * x_lower, double * f_lower, double * x_upper, double * f_upper);
  }
gsl_min_fminimizer_type;

typedef struct
  {
    const gsl_min_fminimizer_type * type;
    gsl_function * function ;
    double x_minimum ;
    double x_lower ;
    double x_upper ;
    double f_minimum, f_lower, f_upper;
    void *state;
  }
gsl_min_fminimizer;

GSL_FUN gsl_min_fminimizer *
gsl_min_fminimizer_alloc (const gsl_min_fminimizer_type * T) ;

GSL_FUN void gsl_min_fminimizer_free (gsl_min_fminimizer * s);

GSL_FUN int gsl_min_fminimizer_set (gsl_min_fminimizer * s,
                            gsl_function * f, double x_minimum,
                            double x_lower, double x_upper);

GSL_FUN int gsl_min_fminimizer_set_with_values (gsl_min_fminimizer * s,
                                        gsl_function * f,
                                        double x_minimum, double f_minimum,
                                        double x_lower, double f_lower,
                                        double x_upper, double f_upper);

GSL_FUN int gsl_min_fminimizer_iterate (gsl_min_fminimizer * s);

GSL_FUN const char * gsl_min_fminimizer_name (const gsl_min_fminimizer * s);

GSL_FUN double gsl_min_fminimizer_x_minimum (const gsl_min_fminimizer * s);
GSL_FUN double gsl_min_fminimizer_x_lower (const gsl_min_fminimizer * s);
GSL_FUN double gsl_min_fminimizer_x_upper (const gsl_min_fminimizer * s);
GSL_FUN double gsl_min_fminimizer_f_minimum (const gsl_min_fminimizer * s);
GSL_FUN double gsl_min_fminimizer_f_lower (const gsl_min_fminimizer * s);
GSL_FUN double gsl_min_fminimizer_f_upper (const gsl_min_fminimizer * s);

/* Deprecated, use x_minimum instead */
GSL_FUN double gsl_min_fminimizer_minimum (const gsl_min_fminimizer * s);

GSL_FUN int
gsl_min_test_interval (double x_lower, double x_upper, double epsabs, double epsrel);

GSL_VAR const gsl_min_fminimizer_type  * gsl_min_fminimizer_goldensection;
GSL_VAR const gsl_min_fminimizer_type  * gsl_min_fminimizer_brent;
GSL_VAR const gsl_min_fminimizer_type  * gsl_min_fminimizer_quad_golden;

typedef
int (*gsl_min_bracketing_function)(gsl_function *f,
                                   double *x_minimum,double * f_minimum,
                                   double *x_lower, double * f_lower,
                                   double *x_upper, double * f_upper,
                                   size_t eval_max);

GSL_FUN int
gsl_min_find_bracket(gsl_function *f,double *x_minimum,double * f_minimum,
                     double *x_lower, double * f_lower,
                     double *x_upper, double * f_upper,
                     size_t eval_max);

__END_DECLS

#endif /* __GSL_MIN_H__ */

max_stars_repo_name: mgreter/astrometrylib

max_issues_repo_name: zzpwahaha/VimbaCamJILA

max_forks_repo_name: zzpwahaha/VimbaCamJILA


text:
#pragma once
/*
 * (C) Copyright 2021 UCAR
 *
 * This software is licensed under the terms of the Apache Licence Version 2.0
 * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
 */
/*! \addtogroup ioda_internals_engines_hh
 *
 * @{
 * \file HH-util.h
 * \brief Utility functions for HDF5.
 */

#include <string>
#include <utility>
#include <vector>

#include <hdf5.h>
#include <gsl/gsl-lite.hpp>

#include "./Handles.h"
#include "ioda/defs.h"
#include "ioda/Exception.h"
#include "ioda/Types/Marshalling.h"

namespace ioda {
namespace detail {
namespace Engines {
namespace HH {
class HH_Attribute;
class HH_Variable;

//#if H5_VERSION_GE(1, 12, 0)
//typedef H5R_ref_t ref_t;
//#else
typedef hobj_ref_t ref_t;
//#endif

// brief Conveys information that a variable (or scale) is attached along a specified axis.
//typedef std::pair<ref_t, unsigned> ref_axis_t;

/// @brief Duplicate the HDF5 dataset list structure for REFERENCE_LISTs.
struct ds_list_t {
  hobj_ref_t ref;       /* object reference  */
  unsigned int dim_idx; /* dimension index of the dataset */
};

/// Data to pass to/from iterator classes.
struct IODA_HIDDEN Iterator_find_attr_data_t {
  std::string search_for;
  hsize_t idx  = 0;
  bool success = false;
};

/// Callback function for H5Aiterate / H5Aiterate2 / H5Aiterate1.
#if H5_VERSION_GE(1, 8, 0)
IODA_HIDDEN herr_t iterate_find_attr(hid_t loc_id, const char* name, const H5A_info_t* info,
                                     void* op_data);
#else
IODA_HIDDEN herr_t iterate_find_attr(hid_t loc_id, const char* name, void* op_data);
#endif

/*! @brief Determine attribute creation order for a dataset.
 *  @param obj is the dataset being queried
 *  @param objType is the type of object (Dataset or Group).
 *  @return H5_INDEX_CRT_ORDER if creation order is tracked.
 *  @return H5_INDEX_NAME if creation order is not tracked.
 *  @details Check if the variable has attribute creation order stored and/or indexed.
 * This is not the default, but it can speed up object list accesses considerably.
 */
IODA_HIDDEN H5_index_t getAttrCreationOrder(hid_t obj, H5O_type_t objType);

/// @brief Linear search to find an attribute.
/// @param baseObject is the object that could contain the attribute.
/// @param attname is the name of the attribute.
/// @param iteration_type is the type of iteration for the search. See getAttrCreationOrder.
/// @return A pair of (success_flag, index).
IODA_HIDDEN std::pair<bool, hsize_t> iterativeAttributeSearch(hid_t baseObject,
                                                              const char* attname,
                                                              H5_index_t iteration_type);

/// @brief Linear search to find and open an attribute, if it exists.
/// @param baseObject is the object that could contain the attribute.
/// @param objType is the type of object (Dataset or Group).
/// @param attname is the name of the attribute.
/// @return An open handle to the attribute, if it exists.
/// @return An invalid handle if the attribute does not exist or upon general failure.
/// @details This function is useful because it is faster than the regular attribute
///   open by name routine, which does not take advantage of attribute creation order
///   indexing. Performance is particularly good when there are few attributes attached
///   to the base object.
IODA_HIDDEN HH_Attribute iterativeAttributeSearchAndOpen(hid_t baseObject, H5O_type_t objType,
                                                         const char* attname);

/*! @brief Attribute DIMENSION_LIST update function
*
* @details This function exists to update DIMENSION_LISTs without updating the
* mirrored REFERENCE_LIST entry in the variable's scales. This is done for
* performance reasons, as attaching dimension scales for hundreds of
* variables sequentially is very slow.
*
* NOTE: This code does not use the regular atts.open(...) call
* for performance reasons when we have to repeat this call for hundreds or
* thousands of variables. We instead do a creation-order-preferred search.
*
* @param var is the variable of interest.
* @param new_dim_list is the mapping of dimensions that should be added to the variable.
*/
IODA_HIDDEN void attr_update_dimension_list(HH_Variable* var,
                                            const std::vector<std::vector<ref_t>>& new_dim_list);

/*! @brief Attribute REFERENCE_LIST update function
*
* @details This function exists to update REFERENCE_LISTs without updating the
* mirrored DIMENSION_LIST entry. This is done for
* performance reasons, as attaching dimension scales for hundreds of
* variables sequentially is very slow.
*
* NOTE: This code does not use the regular atts.open(...) call
* for performance reasons when we have to repeat this call for hundreds or
* thousands of variables. We instead do a creation-order-preferred search.
*
* @param scale is the scale of interest.
* @param ref_var_axis_list is the mapping of variables-dimension numbers that
*   should be added to the scale's REFERENCE_LIST attribute.
*/
IODA_HIDDEN void attr_update_reference_list(HH_Variable* scale,
                                            const std::vector<ds_list_t>& ref_var_axis);

/// @brief A "view" of hvl_t objects. Adds C++ conveniences to an otherwise troublesome class.
/// @tparam Inner is the datatype that we are really manipulating.
template <class Inner>
struct IODA_HIDDEN View_hvl_t {
  hvl_t &obj;
  View_hvl_t(hvl_t& obj) : obj{obj} {}
  size_t size() const { return obj.len; }
  void resize(size_t newlen) {
    if (newlen) {
      obj.p = (obj.p) ? H5resize_memory(obj.p, newlen * sizeof(Inner))
        : H5allocate_memory(newlen * sizeof(Inner), false);
      if (!obj.p) throw Exception("Failed to allocate memory", ioda_Here());
    }
    else {
      if (obj.p)
        if (H5free_memory(obj.p) < 0) throw Exception("Failed to free memory", ioda_Here());
      obj.p = nullptr;
    }
    obj.len = newlen;
  }
  void clear() { resize(0); }
  Inner* at(size_t i) {
    if (i >= obj.len) throw Exception("Out-of-bounds access", ioda_Here())
      .add("i", i).add("obj.len", obj.len);
    return operator[](i);
  }
  Inner* operator[](size_t i) {
    return &(static_cast<Inner*>(obj.p)[i]);
  }
};

/*! Internal structure to encapsulate resources and prevent leaks.
*
* When reading dimension scales, calls to H5Aread return variable-length arrays
* that must be reclaimed. We use a custom object to encapsulate this.
*/
struct IODA_HIDDEN Vlen_data {
  std::unique_ptr<hvl_t[]> buf;  // NOLINT: C array and visibility warnings.
  HH_hid_t typ, space;           // NOLINT: C visibility warnings.
  size_t sz;
  Vlen_data(size_t sz, HH_hid_t typ, HH_hid_t space)
      : buf(new hvl_t[sz]), typ{typ}, space{space}, sz{sz} {
    if (!buf) throw Exception("Failed to allocate buf", ioda_Here());
    for (size_t i = 0; i < sz; i++) {
      buf[i].len = 0;
      buf[i].p   = nullptr;
    }
  }
  ~Vlen_data() {
    if (buf) {
      /*
      for (size_t i = 0; i < sz; i++) {
        if (buf[i].p) delete[] buf[i].p;
        buf[i].len = 0;
        buf[i].p   = nullptr;
      }
      */

      H5Dvlen_reclaim(typ.get(), space.get(), H5P_DEFAULT, reinterpret_cast<void*>(buf.get()));
    }
  }
  Vlen_data(const Vlen_data&) = delete;
  Vlen_data(Vlen_data&&)      = delete;
  Vlen_data operator=(const Vlen_data&) = delete;
  Vlen_data operator=(Vlen_data&&) = delete;

  hvl_t& operator[](size_t idx) { return (buf).get()[idx]; }
};

/// @brief Gets a variable / group / link name from an id. Useful for debugging.
/// @param obj_id is the object.
/// @return One of the possible object names.
/// @throws ioda::Exception if obj_id is invalid.
IODA_HIDDEN std::string getNameFromIdentifier(hid_t obj_id);

/// @brief Convert from variable-length data to fixed-length data.
/// @param in_buf is the input buffer. Buffer has a sequence of pointers, serialized as chars.
/// @param unitLength is the length of each fixed-length element.
/// @param lengthsAreExact signifies that padding is not needed between elements. Each
///   variable-length string is already of length unitLength. If false, then strlen
///   is calculated for each element (which finds the first null byte).
/// @returns the output buffer.
IODA_HIDDEN std::vector<char> convertVariableLengthToFixedLength(
  gsl::span<const char> in_buf, size_t unitLength, bool lengthsAreExact);

/// @brief Convert from fixed-length data to variable-length data.
/// @param in_buf is the input buffer. Buffer is a sequence of fixed-length elements (*not* pointers).
/// @param unitLength is the length of each fixed-length element.
/// @returns the converted buffer.
IODA_HIDDEN Marshalled_Data<char*, char*, true> convertFixedLengthToVariableLength(
  gsl::span<const char> in_buf, size_t unitLength);

}  // namespace HH
}  // namespace Engines
}  // namespace detail
}  // namespace ioda

max_stars_repo_name: NOAA-EMC/ioda

max_issues_repo_name: NOAA-EMC/ioda

max_forks_repo_name: NOAA-EMC/ioda


text:
#ifndef VKST_WSI_INPUT_H
#define VKST_WSI_INPUT_H

#include <gsl.h>
#include <wsi/types.h>
#include <bitset>

namespace wsi {

enum class keys {
  eUnknown = 0,

  eSpace = 32,
  eApostrophe = 39,
  eComma = 44,
  eMinus = 45,
  ePeriod = 46,
  eSlash = 47,
  e0 = 48,
  e1 = 49,
  e2 = 50,
  e3 = 51,
  e4 = 52,
  e5 = 53,
  e6 = 54,
  e7 = 55,
  e8 = 56,
  e9 = 57,
  eSemicolon = 59,
  eEqual = 61,
  eA = 65,
  eB = 66,
  eC = 67,
  eD = 68,
  eE = 69,
  eF = 70,
  eG = 71,
  eH = 72,
  eI = 73,
  eJ = 74,
  eK = 75,
  eL = 76,
  eM = 77,
  eN = 78,
  eO = 79,
  eP = 80,
  eQ = 81,
  eR = 82,
  eS = 83,
  eT = 84,
  eU = 85,
  eV = 86,
  eW = 87,
  eX = 88,
  eY = 89,
  eZ = 90,
  eLeftBracket = 91,
  eBackslash = 92,
  eRightBracket = 93,
  eGraveAccent = 96,

  eEscape = 156,
  eEnter = 157,
  eTab = 158,
  eBackspace = 159,
  eInsert = 160,
  eDelete = 161,
  eRight = 162,
  eLeft = 163,
  eDown = 164,
  eUp = 165,
  ePageUp = 166,
  ePageDown = 167,
  eHome = 168,
  eEnd = 169,

  eCapsLock = 180,
  eScrollLock = 181,
  eNumLock = 182,
  ePrintScreen = 183,
  ePause = 184,

  eF1 = 190,
  eF2 = 191,
  eF3 = 192,
  eF4 = 193,
  eF5 = 194,
  eF6 = 195,
  eF7 = 196,
  eF8 = 197,
  eF9 = 198,
  eF10 = 199,
  eF11 = 200,
  eF12 = 201,
  eF13 = 202,
  eF14 = 203,
  eF15 = 204,
  eF16 = 205,
  eF17 = 206,
  eF18 = 207,
  eF19 = 208,
  eF20 = 209,
  eF21 = 210,
  eF22 = 211,
  eF23 = 212,
  eF24 = 213,
  eF25 = 214,

  eKeypad0 = 220,
  eKeypad1 = 221,
  eKeypad2 = 222,
  eKeypad3 = 223,
  eKeypad4 = 224,
  eKeypad5 = 225,
  eKeypad6 = 226,
  eKeypad7 = 227,
  eKeypad8 = 228,
  eKeypad9 = 229,
  eKeypadDecimal = 230,
  eKeypadDivide = 231,
  eKeypadMultiply = 232,
  eKeypadSubtract = 233,
  eKeypadAdd = 234,
  eKeypadEnter = 235,
  eKeypadEqual = 236,

  eLeftShift = 240,
  eLeftControl = 241,
  eLeftAlt = 242,
  eLeftSuper = 243,
  eRightShift = 244,
  eRightControl = 245,
  eRightAlt = 246,
  eRightSuper = 247,
  eMenu = 248,
}; // enum class keys

class keyset final {
  std::bitset<256> _keys;

public:
  using reference = decltype(_keys)::reference;

  bool all() const { return _keys.all(); }
  bool any() const { return _keys.any(); }
  bool none() const { return _keys.none(); }

  constexpr bool operator[](keys key) const {
    return _keys[static_cast<std::size_t>(key)];
  }

  reference operator[](keys key) {
    return _keys[static_cast<std::size_t>(key)];
  }

  keyset& set(keys key, bool value = true) noexcept {
    _keys[static_cast<std::size_t>(key)] = value;
    return *this;
  }

  keyset& reset(keys key) noexcept {
    _keys[static_cast<std::size_t>(key)] = false;
    return *this;
  }
}; // class keyset

enum class buttons {
  e1 = 1,
  e2,
  e3,
  e4,
  e5,
  e6,
  e7,
  e8,
  e9,
  e10,
}; // enum class buttons

class buttonset final {
  std::bitset<32> _buttons;

public:
  using reference = decltype(_buttons)::reference;

  bool all() const { return _buttons.all(); }
  bool any() const { return _buttons.any(); }
  bool none() const { return _buttons.none(); }

  constexpr bool operator[](buttons button) const {
    return _buttons[static_cast<std::size_t>(button)];
  }

  reference operator[](buttons button) {
    return _buttons[static_cast<std::size_t>(button)];
  }

  buttonset& set(buttons button, bool value = true) noexcept {
    _buttons[static_cast<std::size_t>(button)] = value;
    return *this;
  }

  buttonset& reset(buttons button) noexcept {
    _buttons[static_cast<std::size_t>(button)] = false;
    return *this;
  }
}; // class buttonset

class window;

class input final {
public:
  input(gsl::not_null<window*> win) noexcept;

  bool key_pressed(keys key) const noexcept {
    return !_prev_keys[key] && _curr_keys[key];
  }

  bool key_released(keys key) const noexcept {
    return _prev_keys[key] && !_curr_keys[key];
  }

  bool key_down(keys key) const noexcept {
    return _prev_keys[key] && _curr_keys[key];
  }

  bool button_pressed(buttons button) const noexcept {
    return !_prev_buttons[button] && _curr_buttons[button];
  }

  bool button_released(buttons button) const noexcept {
    return _prev_buttons[button] && !_curr_buttons[button];
  }

  bool button_down(buttons button) const noexcept {
    return _prev_buttons[button] && _curr_buttons[button];
  }

  int prev_scroll() const noexcept { return _prev_scroll; }
  int curr_scroll() const noexcept { return _curr_scroll; }

  int scroll_delta() const noexcept { return _curr_scroll - _prev_scroll; }

  offset2d cursor_pos() const noexcept;

  offset2d prev_cursor_pos() const noexcept { return _prev_pos; }
  offset2d curr_cursor_pos() const noexcept { return _curr_pos; }

  offset2d cursor_delta() const noexcept { return _curr_pos - _prev_pos; }

  void tick() noexcept;

private:
  window* _win;
  keyset _prev_keys, _curr_keys;
  buttonset _prev_buttons, _curr_buttons;
  int _prev_scroll, _curr_scroll;
  offset2d _prev_pos, _curr_pos;
}; // class input

} // namespace wsi

#endif // VKST_WSI_INPUT_H

max_stars_repo_name: wesleygriffin/vkst

max_issues_repo_name: wesleygriffin/vkst

max_forks_repo_name: wesleygriffin/vkst


text:
#ifdef DFFTW
#include <dfftw.h>
#endif

#ifdef SFFTW
#include <sfftw.h>
#endif

#ifdef FFTW
#include <fftw.h>
#endif

int
main()
{
#ifdef DOUBLE
  /* Cause a compile-time error if fftw_real is not double */
  int _array_ [1 - 2 * !((sizeof(fftw_real)) == sizeof(double))];
#else
  /* Cause a compile-time error if fftw_real is not float */
  int _array_ [1 - 2 * !((sizeof(fftw_real)) == sizeof(float))];
#endif
  return 0;
}


text:
#pragma once
#ifndef _NOGSL
#include <gsl/gsl_vector.h>
#else
#include "CustomSolver.h"
#endif
#include "ValueGrad.h"
#include "BooleanNodes.h"
#include "BooleanDAG.h"
#include "NodeVisitor.h"
#include "VarStore.h"
#include <map>
#include "FloatSupport.h"
#include "SymbolicEvaluator.h"

#include <iostream>

using namespace std;

// AutoDiff through only the numerical structure - ignores all boolean structure
class IteAutoDiff: public NodeVisitor, public SymbolicEvaluator
{
	FloatManager& floats;
	BooleanDAG& bdag;
	map<string, int> floatCtrls; // Maps float ctrl names to indices within grad vector
	int nctrls; // number of float ctrls
	gsl_vector* ctrls; // ctrl values
	vector<ValueGrad*> values; // Keeps track of values along with gradients for each node
	map<int, int> inputValues; // Maps node id to values set by the SAT solver
    map<int, int> iteCtrls;
	double error = 0.0;
	gsl_vector* errorGrad;

	int DEFAULT_INP = -1;

public:
	IteAutoDiff(BooleanDAG& bdag_p, FloatManager& _floats, const map<string, int>& floatCtrls_p, const map<int, int>& iteCtrls_p);
	~IteAutoDiff(void);

	virtual void visit( SRC_node& node );
	virtual void visit( DST_node& node );
	virtual void visit( CTRL_node& node );
	virtual void visit( PLUS_node& node );
	virtual void visit( TIMES_node& node );
	virtual void visit( ARRACC_node& node );
	virtual void visit( DIV_node& node );
	virtual void visit( MOD_node& node );
	virtual void visit( NEG_node& node );
	virtual void visit( CONST_node& node );
	virtual void visit( LT_node& node );
	virtual void visit( EQ_node& node );
	virtual void visit( AND_node& node );
	virtual void visit( OR_node& node );
	virtual void visit( NOT_node& node );
	virtual void visit( ARRASS_node& node );
	virtual void visit( UFUN_node& node );
	virtual void visit( TUPLE_R_node& node );
	virtual void visit( ASSERT_node& node );

	virtual void run(const gsl_vector* ctrls_p, const map<int, int>& inputValues_p);
    virtual bool checkAll(const gsl_vector* ctrls_p, const map<int, int>& inputValues_p, bool onlyBool = false){
        Assert(false, "NYI");
        return true;
    }
	virtual bool check(bool_node* n, int expected);
	virtual double computeError(bool_node* n, int expected, gsl_vector* errorGrad);
	virtual double computeDist(bool_node*, gsl_vector* distgrad);
	virtual bool hasDist(bool_node* n);
    virtual double computeVal(bool_node*, gsl_vector* distgrad);
    virtual bool hasVal(bool_node* n);
    virtual bool hasSqrtDist(bool_node* n);
    virtual double computeSqrtError(bool_node* n, gsl_vector* errorGrad);
    virtual double computeSqrtDist(bool_node* n, gsl_vector* errorGrad);
    virtual double getVal(bool_node* n);
    virtual gsl_vector* getGrad(bool_node* n);
    //virtual set<int> getConflicts(int nid);

	void setvalue(bool_node& bn, ValueGrad* v) {
		values[bn.id] = v;
	}

	ValueGrad* v(bool_node& bn) {
		ValueGrad* val = values[bn.id];
		if (val == NULL) {
			gsl_vector* g = gsl_vector_alloc(nctrls);
			val = new ValueGrad(0, g);
			setvalue(bn, val);
		}
		return val;
	}

	ValueGrad* v(bool_node* bn) {
		return v(*bn);
	}

	virtual void print() {
		for (int i = 0; i < bdag.size(); i++) {
			cout << bdag[i]->lprint() << " ";
			ValueGrad* val = v(bdag[i]);
			if (val->set) {
				cout << val->print() << endl;
			} else {
				cout << "UNSET" << endl;
			}
		}
	}
	virtual void printFull() {
		for (int i = 0; i < bdag.size(); i++) {
			cout << bdag[i]->lprint() << endl;
			ValueGrad* val = v(bdag[i]);
			if (val->set) {
				cout << val->printFull() << endl;
			} else {
				cout << "UNSET" << endl;
			}
		}
	}

	bool isFloat(bool_node& bn) {
		return (bn.getOtype() == OutType::FLOAT);
	}

	bool isFloat(bool_node* bn) {
		return (bn->getOtype() == OutType::FLOAT);
	}

	int getInputValue(bool_node& bn) {
		if (inputValues.find(bn.id) != inputValues.end()) {
			int val = inputValues[bn.id];
			//Assert(val == 0 || val == 1, "NYI: Integer values");
			return val;
		} else {
			return DEFAULT_INP;
		}
	}

	int getInputValue(bool_node* bn) {
		return getInputValue(*bn);
	}

};

max_stars_repo_name: natebragg/sketch-backend

max_issues_repo_name: natebragg/sketch-backend

max_forks_repo_name: natebragg/sketch-backend


text:
#include <gsl/gsl_math.h>
#include "gsl_cblas.h"
#include "cblas.h"

void
cblas_dspmv (const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
	     const int N, const double alpha, const double *Ap,
	     const double *X, const int incX, const double beta, double *Y,
	     const int incY)
{
#define BASE double
#include "source_spmv.h"
#undef BASE
}


text:
//===--- Sudoku/Solver_remove_option.h                                  ---===//
//
// Helper functions, to remove options
//===----------------------------------------------------------------------===//
#pragma once

#include "Sudoku/Board.h"
#include "Sudoku/Location.h"
#include "Sudoku/Location_Utilities.h"
#include "Sudoku/Options.h"
#include "Sudoku/Size.h"
#include "Sudoku/Value.h"
#include "Sudoku/exceptions.h"
#include "Sudoku/traits.h"
#include <gsl/gsl>

#include <vector>

#include <algorithm>   // none_of
#include <type_traits> // is_base_of

#include "Solver.fwd.h" // Forward declaration - single_option

#include <cassert>

/* experiment flags */
//! DO #undef at end of file!
// activate algorithms on removing option
#define DUAL_ON_REMOVE false     // NOLINT(cppcoreguidelines-macro-usage)
#define MULTIPLE_ON_REMOVE false // NOLINT(cppcoreguidelines-macro-usage)

namespace Sudoku
{
//===----------------------------------------------------------------------===//
template<int N, typename Options = Options<elem_size<N>>>
int remove_option(Board<Options, N>&, Location<N>, Value);
template<int N, typename Options = Options<elem_size<N>>>
int remove_option(Board<Options, N>&, Location<N>, Options mask);

template<int N, typename Options = Options<elem_size<N>>, typename SectionT>
int remove_option_section(
	Board<Options, N>&, SectionT, Location<N> ignore, Value);
template<int N, typename Options = Options<elem_size<N>>, typename SectionT>
int remove_option_section(
	Board<Options, N>&,
	SectionT,
	const std::vector<Location<N>>& ignore,
	Value);
template<int N, typename Options = Options<elem_size<N>>, typename SectionT>
int remove_option_section(
	Board<Options, N>&,
	SectionT,
	const std::vector<Location<N>>& ignore,
	const Options& values);

template<int N, typename Options = Options<elem_size<N>>, typename SectionT>
int remove_option_outside_block(
	Board<Options, N>&, SectionT, Location<N> block_loc, Value);

//===----------------------------------------------------------------------===//

//	remove option from element, make answer if last option
template<int N, typename Options>
int remove_option(
	Board<Options, N>& board, // NOLINT(runtime/references)
	const Location<N> loc,
	const Value value)
{
	assert(is_valid<N>(value));

	int changes{};
	auto& item{board.at(loc)};

	if (item.test(value))
	{ // repeated answer in section (faster here)
		if (!is_option(item, value))
			throw error::invalid_Board();

		++changes;
		const auto count = item.remove_option(value).count();
		assert(count > 0); // never trigger, removed last option

#if MULTIPLE_ON_REMOVE == false
		if (count == 1)
		{
			changes += single_option(board, loc);
		}
#else
		// if (count < 4)
		{
			changes += multi_option(board, loc, count);
		}
#endif // multiple
#if DUAL_ON_REMOVE == true
		if (count == 2)
		{
			changes += dual_option(board, loc);
		}
#endif // dual
	}
	return changes;
}

// remove all options given by the mask
template<int N, typename Options>
int remove_option(
	Board<Options, N>& board, // NOLINT(runtime/references)
	const Location<N> loc,
	const Options mask)
{
	assert(is_answer_fast(mask)); // don't remove answer-bit
	assert(!mask.is_empty());     // useless

	auto& item{board.at(loc)};
	auto changes = gsl::narrow_cast<int, size_t>(item.count_all());
	if (!changes)
	{ // ignore empty elements
		return 0;
	}

	// remove options
	item -= mask;

	const auto count = gsl::narrow_cast<int, size_t>(item.count_all());
	if (count == 0)
	{ // removed last option
		throw error::invalid_Board();
	}
	changes -= count;

	if (changes && count == 1)
	{
		changes += single_option(board, loc);
	}
	return changes;
}

//	remove [value] in [ignore] from other elements in set
template<int N, typename Options, typename SectionT>
int remove_option_section(
	Board<Options, N>& board, // NOLINT(runtime/references)
	const SectionT section,
	const Location<N> ignore,
	const Value value)
{
	{
		static_assert(Board_Section::traits::is_Section_v<SectionT>);
		static_assert(std::is_same_v<Options, typename SectionT::value_type>);
		static_assert(traits::is_input<typename SectionT::iterator>);
		assert(is_valid(ignore));
		assert(is_valid<N>(value));
		assert(is_same_section(section, ignore));
		assert(is_answer(board.at(ignore), value)); // first set as answer!
	}
	int changes{0};
	const auto end = section.cend();

	for (auto itr = section.cbegin(); itr != end; ++itr)
	{
		if (itr.location() != ignore)
		{
			changes += remove_option(board, itr.location(), value);
		}
	}
	return changes;
}

//	remove [value] from set if not part of [loc]s
template<int N, typename Options, typename SectionT>
inline int remove_option_section(
	Board<Options, N>& board, // NOLINT(runtime/references)
	const SectionT section,
	const std::vector<Location<N>>& ignore,
	const Value value)
{
	{
		static_assert(Board_Section::traits::is_Section_v<SectionT>);
		static_assert(std::is_same_v<Options, typename SectionT::value_type>);
		static_assert(traits::is_input<typename SectionT::iterator>);
		assert(is_valid(ignore));
		assert(is_valid<N>(value));
		assert(is_same_section(section, ignore));
	}
	int changes{0};
	const auto begin = ignore.cbegin();
	const auto end   = ignore.cend();

	// TODO maybe faster to run get_same_row/col/block first?
	// and than not to check is_option(), since it is already in remove_option

	for (auto itr = section.cbegin(); itr != section.cend(); ++itr)
	{
		// TODO is the is_option check really faster?
		if (std::none_of(
				begin,
				end,
				[L1 = itr.location()](Location<N> L2) { return L1 == L2; }) &&
			is_option(*itr, value)) // <algorithm>
		{
			changes += remove_option(board, itr.location(), value);
		}
	}
	return changes;
}

// remove [value]s from [section] if not part of [loc]s
// the ignore Location vector must be sorted
template<int N, typename Options, typename SectionT>
int remove_option_section(
	Board<Options, N>& board, // NOLINT(runtime/references)
	const SectionT section,
	const std::vector<Location<N>>& ignore,
	const Options& values)
{
	{
		static_assert(Board_Section::traits::is_Section_v<SectionT>);
		static_assert(std::is_same_v<Options, typename SectionT::value_type>);
		static_assert(traits::is_input<typename SectionT::iterator>);
		assert(is_valid(ignore));
		assert(!values.all());
		assert(!values.is_empty());
		assert(is_same_section(section, ignore));
	}
	int changes{0};
	const auto begin = ignore.cbegin();
	const auto end   = ignore.cend();

	for (auto itr = section.cbegin(); itr != section.cend(); ++itr)
	{
		if (std::none_of(
				begin,
				end,
				[L1 = itr.location()](Location<N> L2) { return L1 == L2; }) &&
			not(itr->is_answer())) // <algorithm>
		{
			changes += remove_option(board, itr.location(), values);
		}
	}
	return changes;
}

//	remove [value] from elements not in same block as [block_loc]
template<int N, typename Options, typename SectionT>
int remove_option_outside_block(
	Board<Options, N>& board, // NOLINT(runtime/references)
	const SectionT section,
	const Location<N> block_loc,
	const Value value)
{
	{
		static_assert(Board_Section::traits::is_Section_v<SectionT>);
		static_assert(
			!Board_Section::traits::is_Block_v<SectionT>,
			"remove_option_outside_block is useless on a Block");
		static_assert(std::is_same_v<Options, typename SectionT::value_type>);
		static_assert(traits::is_input<typename SectionT::iterator>);
		assert(is_valid(block_loc));
		assert(is_valid<N>(value));
		assert(intersect_block(section, block_loc));
	}
	int changes{0};
	const auto end = section.cend();

	for (auto itr = section.cbegin(); itr != end; ++itr)
	{
		if (not is_same_block(itr.location(), block_loc))
		{
			changes += remove_option(board, itr.location(), value);
		}
	}
	return changes;
}

/*	Remove local macros */
#undef DUAL_ON_REMOVE
#undef MULTIPLE_ON_REMOVE

} // namespace Sudoku

max_stars_repo_name: Farwaykorse/fwkSudoku

max_issues_repo_name: Farwaykorse/fwkSudoku

max_forks_repo_name: Farwaykorse/fwkSudoku


text:
#include <stdio.h>
#include <stdarg.h>
#include <string.h>
#include <math.h>
#include <gbpLib.h>
#include <gbpRNG.h>
#include <gbpMCMC.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_fit.h>
#include <gsl/gsl_interp.h>

void compute_MCMC_ln_likelihood_default(MCMC_info *MCMC,
                                        double **  M,
                                        double *   P,
                                        double *   ln_likelihood_DS,
                                        int *      n_DoF_DS,
                                        double *   ln_likelihood_all,
                                        int *      n_DoF_all) {
    int           i_DS;
    int           i_M;
    int           n_M;
    double *      M_target;
    double *      dM_target;
    MCMC_DS_info *current_DS;
    MCMC_DS_info *next_DS;
    (*ln_likelihood_all) = 0.;
    (*n_DoF_all)         = 0;
    i_DS                 = 0;
    current_DS           = MCMC->DS;
    while(current_DS != NULL) {
        next_DS                = current_DS->next;
        n_M                    = current_DS->n_M;
        M_target               = current_DS->M_target;
        dM_target              = current_DS->dM_target;
        ln_likelihood_DS[i_DS] = 0.;
        n_DoF_DS[i_DS]         = 0;
        for(i_M = 0; i_M < n_M; i_M++) {
            ln_likelihood_DS[i_DS] += pow((M_target[i_M] - M[i_DS][i_M]) / dM_target[i_M], 2.);
            n_DoF_DS[i_DS]++;
        }
        (*ln_likelihood_all) += ln_likelihood_DS[i_DS];
        ln_likelihood_DS[i_DS] *= -0.5;
        (*n_DoF_all) += n_DoF_DS[i_DS];
        n_DoF_DS[i_DS] -= MCMC->n_P;
        i_DS++;
        current_DS = next_DS;
    }
    (*n_DoF_all) -= MCMC->n_P;
    (*ln_likelihood_all) *= -0.5;
}

max_stars_repo_name: gbpoole/gbpCode

max_issues_repo_name: gbpoole/gbpCode

max_forks_repo_name: gbpoole/gbpCode


text:
/* randist/sphere.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004, 2007 James Theiler, 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 <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

void
gsl_ran_dir_2d (const gsl_rng * r, double *x, double *y)
{
  /* This method avoids trig, but it does take an average of 8/pi =
   * 2.55 calls to the RNG, instead of one for the direct
   * trigonometric method.  */

  double u, v, s;
  do
    {
      u = -1 + 2 * gsl_rng_uniform (r);
      v = -1 + 2 * gsl_rng_uniform (r);
      s = u * u + v * v;
    }
  while (s > 1.0 || s == 0.0);

  /* This is the Von Neumann trick. See Knuth, v2, 3rd ed, p140
   * (exercise 23).  Note, no sin, cos, or sqrt !  */

  *x = (u * u - v * v) / s;
  *y = 2 * u * v / s;

  /* Here is the more straightforward approach,
   *     s = sqrt (s);  *x = u / s;  *y = v / s;
   * It has fewer total operations, but one of them is a sqrt */
}

void
gsl_ran_dir_2d_trig_method (const gsl_rng * r, double *x, double *y)
{
  /* This is the obvious solution... */
  /* It ain't clever, but since sin/cos are often hardware accelerated,
   * it can be faster -- it is on my home Pentium -- than von Neumann's
   * solution, or slower -- as it is on my Sun Sparc 20 at work
   */
  double t = 6.2831853071795864 * gsl_rng_uniform (r);          /* 2*PI */
  *x = cos (t);
  *y = sin (t);
}

void
gsl_ran_dir_3d (const gsl_rng * r, double *x, double *y, double *z)
{
  double s, a;

  /* This is a variant of the algorithm for computing a random point
   * on the unit sphere; the algorithm is suggested in Knuth, v2,
   * 3rd ed, p136; and attributed to Robert E Knop, CACM, 13 (1970),
   * 326.
   */

  /* Begin with the polar method for getting x,y inside a unit circle
   */
  do
    {
      *x = -1 + 2 * gsl_rng_uniform (r);
      *y = -1 + 2 * gsl_rng_uniform (r);
      s = (*x) * (*x) + (*y) * (*y);
    }
  while (s > 1.0);

  *z = -1 + 2 * s;              /* z uniformly distributed from -1 to 1 */
  a = 2 * sqrt (1 - s);         /* factor to adjust x,y so that x^2+y^2
                                 * is equal to 1-z^2 */
  *x *= a;
  *y *= a;
}

void
gsl_ran_dir_nd (const gsl_rng * r, size_t n, double *x)
{
  double d;
  size_t i;
  /* See Knuth, v2, 3rd ed, p135-136.  The method is attributed to
   * G. W. Brown, in Modern Mathematics for the Engineer (1956).
   * The idea is that gaussians G(x) have the property that
   * G(x)G(y)G(z)G(...) is radially symmetric, a function only
   * r = sqrt(x^2+y^2+...)
   */
  d = 0;
  do
    {
      for (i = 0; i < n; ++i)
        {
          x[i] = gsl_ran_gaussian (r, 1.0);
          d += x[i] * x[i];
        }
    }
  while (d == 0);
  d = sqrt (d);
  for (i = 0; i < n; ++i)
    {
      x[i] /= d;
    }
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


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

#ifndef SQ_INCLUDE_GUARD_system_linux_udev_inl_h_
#define SQ_INCLUDE_GUARD_system_linux_udev_inl_h_

#include "core/errors.h"
#include "core/typeutil.h"

#include <fmt/format.h>
#include <gsl/gsl>

namespace sq::system::linux {

template <typename T> UdevPtr<T>::UdevPtr(SQ_MU std::nullptr_t np) noexcept {}

template <typename T>
UdevPtr<T>::UdevPtr(const UdevPtr &other) noexcept : t_{other.t_.p_} {
  if (t_.p_ != nullptr) {
    t_.add_ref();
  }
}

template <typename T>
UdevPtr<T>::UdevPtr(UdevPtr &&other) noexcept : t_{other.t_.p_} {
  other.t_.p_ = nullptr;
}

// When a raw pointer is first obtained from libudev it already has a reference
// count of 1 so we don't need to adjust it here.
template <typename T> UdevPtr<T>::UdevPtr(RawPtr p) noexcept : t_{p} {
  Expects(p != nullptr);
  Ensures(t_.p_ != nullptr);
}

template <typename T>
UdevPtr<T> &UdevPtr<T>::operator=(const UdevPtr &other) noexcept {
  // Note: must call add_ref before remove_ref to deal with the case where
  // *this == other
  if (other.t_.p_ != nullptr) {
    other.t_.add_ref();
  }
  if (t_.p_ != nullptr) {
    t_.remove_ref();
  }
  t_.p_ = other.t_.p_;
  return *this;
}

template <typename T> UdevPtr<T>::~UdevPtr() noexcept {
  if (t_.p_ != nullptr) {
    t_.remove_ref();
  }
}

template <typename T>
UdevPtr<T> &UdevPtr<T>::operator=(UdevPtr &&other) noexcept {
  std::swap(t_.p_, other.t_.p_);
  return *this;
}

template <typename T> void UdevPtr<T>::reset() noexcept {
  UdevPtr<T>{}.swap(*this);
}

template <typename T> void UdevPtr<T>::swap(UdevPtr &other) noexcept {
  std::swap(t_.p_, other.t_.p_);
}

template <typename T> T *UdevPtr<T>::get() const noexcept {
  if (t_.p_ == nullptr) {
    return nullptr;
  }
  return &t_;
}

template <typename T> T &UdevPtr<T>::operator*() const noexcept {
  Expects(t_.p_ != nullptr);
  return t_;
}

template <typename T> T *UdevPtr<T>::operator->() const noexcept {
  return get();
}

template <typename T> UdevPtr<T>::operator bool() const noexcept {
  return t_.p_ != nullptr;
}

template <typename T>
auto UdevPtr<T>::operator<=>(const UdevPtr<T> &other) const noexcept {
  return std::compare_three_way{}(t_.p_, other.t_.p_);
}

template <typename T>
bool UdevPtr<T>::operator==(const UdevPtr<T> &other) const noexcept {
  return t_.p_ == other.t_.p_;
}

template <typename T>
auto UdevPtr<T>::operator<=>(std::nullptr_t np) const noexcept {
  return std::compare_three_way{}(t_.p_, static_cast<decltype(t_.p_)>(np));
}

template <typename T>
bool UdevPtr<T>::operator==(std::nullptr_t np) const noexcept {
  return t_.p_ == static_cast<decltype(t_.p_)>(np);
}

template <typename U, typename... Args> UdevPtr<U> make_udev(Args &&...args) {
  auto *const p = U::create_new(SQ_FWD(args)...);
  if (p == nullptr) {
    throw UdevError{fmt::format("Failed to create new {}", base_type_name(p))};
  }
  return UdevPtr<U>{p};
}

} // namespace sq::system::linux

#endif // SQ_INCLUDE_GUARD_system_linux_udev_inl_h_

max_stars_repo_name: jonathanhaigh/sq

max_issues_repo_name: jonathanhaigh/sq

max_forks_repo_name: jonathanhaigh/sq