text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdarg.h>
#include <math.h>
#include <unistd.h>
#include <limits.h>
#include <sys/time.h>
#include <zlib.h>

#ifdef INTEL_COMPILER
#include "mkl.h"
#else
#include <cblas.h>
#endif

//#define FIBS_UNIT 1000
#define DEFAULT_NDIGITS 10
#define SZBUF 1024
#define SZBYTE 256
#define N_GENOTYPES 4
#define NA_GENO_CHAR (N_GENOTYPES-1)
#define DEFAULT_ROW_SIZE 100000
#define DEFAULT_SIZE_MATRIX 1000000
#define DEFAULT_SIZE_HEADER 100000
#define DEFAULT_DELIMS " \t\r\n"
#define SZ_LONG_BUF 1000000
#define DEFAULT_TPED_NUM_HEADER_COLS 4
#define DEFAULT_TFAM_NUM_HEADER_COLS 6
#define DEFAULT_TPED_SNPID_INDEX 1
#define DEFAULT_PHENO_NUM_HEADER_COLS 2

struct HFILE {
  int gzflag;       // 1 if gz if used
  int wflag;        // r(0)/w(1) for plain, rb(0)/wb(1) for gz
  int nheadercols;  // # of header columns (0 if nrows=0)
  int nvaluecols;   // # of value cols (0 if nrows=0)
  int nrows;        // # of rows
  FILE* fp;         // plain file handle
  gzFile gzfp;      // gzip file handle
};

// Input routines
void close_file (struct HFILE* fhp);
struct HFILE open_file(char* filename, int gzflag, int wflag);
struct HFILE open_file_with_suffix(char* prefix, char* suffix, int gzflag, int wflag);
//void read_matrix_with_col_headers( struct HFILE* fhp, int nheadercols, char* delims, int symmetric, int* p_nmiss, unsigned char** matrix, char*** headers);
void read_matrix_with_col_headers( struct HFILE* fhp, int nheadercols, char* delims, int* p_nmiss, unsigned char** matrix, char*** headers);
unsigned char* tokenize_tped_line_with_col_headers( struct HFILE* fhp, int nheadercols, char* delims, char* lbuf, unsigned char* values, char** headers, int* p_nvalues, int* p_nmiss );

void emmax_error( const char* format, ... );
void print_help(void);
FILE* readfile(char* filename);

void print_help(void) {
  fprintf(stderr,"Usage: emmax_kin [tpedf]\n");
  fprintf(stderr,"Required parameters\n");
  fprintf(stderr,"\t[tpedf]     : tped file\n");
  fprintf(stderr,"Optional parameters\n");
  fprintf(stderr,"\t-d [# digits]  : precision of the kinship values (default : 10)\n");
  fprintf(stderr,"\t-M [float] : maximum memory in GB (default: 4.0)\n");
  fprintf(stderr,"\t-s : compute IBS kinship matrix (default is Balding-Nicholas)\n");
  fprintf(stderr,"\t-v : turn on verbose mode\n");
  fprintf(stderr,"\t-r : randomly fill missing genotypes (default is imputation by average)\n");
  fprintf(stderr,"\t-x : include non-autosomal chromosomes in computing kinship matrices\n");
  fprintf(stderr,"\t-S [int] : set random seed\n");
  fprintf(stderr,"\t-m [float] : MAF threshold (default is 0)\n");
  fprintf(stderr,"\t-c [float] : Call rate threshold (default is 0)\n");
}

//void read_matrix_with_col_headers( struct HFILE* fhp, int nheadercols, char* delims, int symmetric, int* p_nmiss, unsigned char** matrix, char*** headers) {
void read_matrix_with_col_headers( struct HFILE* fhp, int nheadercols, char* delims, int* p_nmiss, unsigned char** matrix, char*** headers) {
  char* lbuf = (char*) malloc(sizeof(char*) * SZ_LONG_BUF);
  int szmat = DEFAULT_SIZE_MATRIX;
  int szheader = DEFAULT_SIZE_HEADER;
  unsigned char* cmat = (unsigned char*) malloc(sizeof(unsigned char) * szmat );
  char** cheaders = (char**) malloc(sizeof(char*) * szheader );
  int nvalues, i, j, nmiss;

  fhp->nheadercols = nheadercols;
  nmiss = 0;

  while( tokenize_tped_line_with_col_headers(fhp, nheadercols, delims, lbuf, &cmat[fhp->nrows*fhp->nvaluecols], &cheaders[fhp->nrows*fhp->nheadercols], &nvalues, &nmiss) != NULL ) {
    if ( fhp->nrows == 1 ) {
      fhp->nvaluecols = nvalues;
    }
    else if ( fhp->nvaluecols != nvalues ) {
      emmax_error("The column size %d do not match to %d at line %d\n",nvalues,fhp->nvaluecols,fhp->nrows);
    }

    if ( (fhp->nrows+1)*(fhp->nvaluecols) > szheader ) {
      szheader *= 2;
      fprintf(stderr,"Header size is doubled to %d\n",szheader);
      cheaders = (char**) realloc( cheaders, sizeof(char*) * szheader );
    }

    if ( (fhp->nrows+1)*(fhp->nvaluecols) > szheader ) {
      szmat *= 2;
      fprintf(stderr,"Matrix size is doubled to %d\n",szmat);
      cmat = (unsigned char*) realloc( cmat, sizeof(unsigned char) * szmat );
    }
  }
  free(lbuf);

  *p_nmiss = nmiss;

  unsigned char* fmat = (unsigned char*) malloc(sizeof(unsigned char)*fhp->nrows*fhp->nvaluecols);
  char** fheaders = (char**) malloc(sizeof(char*)*fhp->nrows*fhp->nheadercols);
  for(i=0; i < fhp->nrows; ++i) {
    for(j=0; j < fhp->nvaluecols; ++j) {
      fmat[i+j*fhp->nrows] = cmat[i*fhp->nvaluecols+j];
    }
    for(j=0; j < fhp->nheadercols; ++j) {
      fheaders[i+j*fhp->nrows] = cheaders[i*fhp->nheadercols+j];
    }
  }
  free(cmat);
  free(cheaders);

  if ( matrix != NULL ) {
    if ( *matrix != NULL ) {
      free(*matrix);
    }
    *matrix = fmat;
  }

  if ( headers != NULL ) {
    if ( *headers != NULL ) {
      free(*headers);
    }
    *headers = fheaders;
  }
}

unsigned char* tokenize_tped_line_with_col_headers( struct HFILE* fhp, int nheadercols, char* delims, char* lbuf, unsigned char* values, char** headers, int* p_nvalues, int* p_nmiss ) {
  int j;
  char *token;
  unsigned char ctoken;

  char *ret = (fhp->gzflag == 1) ? gzgets(fhp->gzfp, lbuf, SZ_LONG_BUF) : fgets( lbuf, SZ_LONG_BUF, fhp->fp );
  int nmiss = 0;

  if ( ret == NULL ) {
    return NULL;
  }

  if ( fhp->nheadercols != nheadercols ) {
    emmax_error("# of header columns mismatch (%d vs %d) at line %d",fhp->nheadercols,nheadercols,fhp->nrows);
  }

  //fprintf(stderr,"tokenize-line called %s\n",lbuf);

  token = strtok(lbuf, delims);
  for( j=0; token != NULL; ++j ) {
    if ( j < nheadercols ) {
      headers[j] = strdup(token);
    }
    // if zero_miss_flag is set, assume the genotypes are encoded 0,1,2
    // Additively encodes the two genotypes in the following way
    // when (j-nheadercols) is even, 0->MISSING, add 1->0, 2->1
    // when (j-nheadercols) is odd, check 0-0 consistency, and add 1->0, 2->1
    else {
      ctoken = (unsigned char)(token[0]-'0');

      if ( ctoken > 2 ) {
	fprintf(stderr,"Unrecognized token %s\n",token);
	abort();
      }

      if ( (j-nheadercols) % 2 == 0 ) {
	values[(j-nheadercols)/2] = ctoken;
      }
      else {
	if ( ( ctoken > 0 ) && ( values[(j-nheadercols)/2] == 0 ) ) {
	  fprintf(stderr,"Unmatched token pair 0 %s\n",token);
	  abort();
	}
	else if ( ( ctoken == 0 ) && ( values[(j-nheadercols)/2] > 0 ) ) {
	  fprintf(stderr,"Unmatched token pair - %d 0\n",(int)values[(j-nheadercols)/2]);
	  abort();
	}
	values[(j-nheadercols)/2] = (unsigned char)(values[(j-nheadercols)/2]+ctoken);
      }
    }
    token = strtok(NULL, delims);
  }
  //fprintf(stderr,"tokenize-line ended %d %d\n",j,nheadercols);
  if ( (j-nheadercols) % 2 != 0 ) {
    fprintf(stderr,"Number of value tokens are not even %d\n",j-nheadercols);
    abort();
  }

  *p_nvalues = (j-nheadercols)/2;
  *p_nmiss += nmiss;
  ++(fhp->nrows);

  if ( j < nheadercols ) {
    fprintf(stderr,"Number of header columns are %d, but only %d columns were observed\n", nheadercols, j);
    abort();
  }

  return values;
}

// open_file_with_suffix()
// - [prefix].[suffix] : file name to open
// - gzflag : gzip flag (use gzfp if gzflag=1, otherwise use fp)
// - wflag : write flag (1 if write mode otherwise read mode
struct HFILE open_file_with_suffix(char* prefix, char* suffix, int gzflag, int wflag) {
  char filename[SZBUF];
  sprintf(filename,"%s.%s",prefix,suffix);
  return open_file(filename,gzflag,wflag);
}

// open_file()
// - filename : file name to open
// - gzflag : gzip flag (use gzfp if gzflag=1, otherwise use fp)
// - wflag : write flag (1 if write mode otherwise read mode)
struct HFILE open_file(char* filename, int gzflag, int wflag) {
  struct HFILE fh;
  fh.gzflag = gzflag;
  fh.wflag = wflag;
  fh.nheadercols = 0;
  fh.nvaluecols = 0;
  fh.nrows = 0;
  if ( gzflag == 1 ) {
    char* mode = (wflag == 1) ? "wb" : "rb";
    fh.gzfp = gzopen(filename,mode);
    fh.fp = NULL;

    if ( fh.gzfp == NULL ) {
      emmax_error("Cannot open file %s for reading",filename);
    }
  }
  else {
    char* mode = (wflag == 1) ? "w" : "r";
    fh.gzfp = (gzFile) NULL;
    fh.fp = fopen(filename,mode);

    if ( fh.fp == NULL ) {
      emmax_error("Cannot open file %s for writing",filename);
    }
  }
  return fh;
}

void emmax_error( const char* format, ... ) {
  va_list args;
  fprintf(stderr, "ERROR: ");
  va_start (args, format);
  vfprintf(stderr, format, args);
  va_end (args);
  fprintf(stderr,"\n");
  abort();
}

void close_file(struct HFILE* fhp) {
  if ( fhp->gzflag == 1 ) {
    gzclose(fhp->gzfp);
    fhp->gzfp = NULL;
  }
  else {
    fclose(fhp->fp);
    fhp->fp = NULL;
  }
}

int main(int argc, char** argv) {
  int i, j, n, c, ac0, ac1, ac2, nmiss, nelems, nex, n_sum_nin, nin, nex_maf, nex_call_rate, nex_autosomal;
  int verbose, ndigits, tped_nheadercols, tfam_nheadercols, flag_autosomal, rand_fill_flag, ibs_flag, n_unit_lines;
  unsigned char *snprow;
  char *suffix, buf[SZBUF];
  double f, max_memory_GB, maf_thres, call_rate_thres, aaf, call_rate;
  //long *fibs_sums, *scores, mean_score;
  double *kin, *snpunit;
  char *tpedf, *delims, *lbuf;
  char **tfam_headers, **tped_headers;
  struct HFILE tpedh, tfamh, kinsh;
  struct timeval tv;

  // set default params
  gettimeofday(&tv, NULL);
  srand((unsigned int)tv.tv_usec);
  delims = DEFAULT_DELIMS;
  tped_nheadercols = DEFAULT_TPED_NUM_HEADER_COLS;
  tfam_nheadercols = DEFAULT_TFAM_NUM_HEADER_COLS;
  tped_headers = tfam_headers = NULL;
  tpedf = lbuf = 0;
  flag_autosomal = 1;
  rand_fill_flag = 0;
  ibs_flag = 0;
  verbose = 0;
  ndigits = DEFAULT_NDIGITS;
  max_memory_GB = 4.0; // 4.0GB
  maf_thres = 0.0;
  call_rate_thres = 0.0;

  // read arguments and update params
  while ((c = getopt(argc, argv, "d:rsc:vxS:M:m:c:")) != -1 ) {
    switch(c) {
    case 'd': // precision of digits
      ndigits = atoi(optarg);
      break;
    case 'r':
      rand_fill_flag = 1;
      break;
    case 's':
      ibs_flag = 1;
      break;
    case 'v':
      verbose = 1;
      break;
    case 'x':
      flag_autosomal = 0;
      break;
    case 'S':
      srand(atoi(optarg));
      break;
    case 'M':
      max_memory_GB = atof(optarg);
      break;
    case 'm':
      maf_thres = atof(optarg);
      break;
    case 'c':
      call_rate_thres = atof(optarg);
      break;
    default:
      fprintf(stderr,"Error : Unknown option unsigned character %c\n",c);
      abort();
    }
  }

  // Sanity check for the number of required parameters
  if ( argc != optind + 1 ) {
    print_help();
    abort();
  }

  // Read required parameters
  tpedf = argv[optind++];

  if ( verbose) fprintf(stderr,"\nReading TFAM file %s.tfam ....\n",tpedf);

  tfamh = open_file_with_suffix(tpedf, "tfam", 0, 0);
  //read_matrix_with_col_headers( &tfamh, tfam_nheadercols, delims, 0, &nmiss, NULL, &tfam_headers);
  read_matrix_with_col_headers( &tfamh, tfam_nheadercols, delims, &nmiss, NULL, &tfam_headers);
  n = tfamh.nrows; // n is # of individuals
  if ( verbose ) fprintf(stderr,"Identified %d individuals from TFAM file\n",n);

  // compute the # of lines to read together
  // we would need two n*n matrix, and n*m matrix
  // (n*n + n*m)*sizeof(double) < M*1e9
  // m < M*1e9/sizeof(double)/n - n
  n_unit_lines = (int)floor((max_memory_GB * 1.0e9 / sizeof(double) / n - n)/2)*2;
  n_sum_nin = 0;
#ifdef INTEL_COMPILER
  char cn = 'N', ct = 'T';
  double one = 1.;
#endif

  if ( verbose ) fprintf(stderr,"Setting # unit lines = %d to fit the memory requirement\n",n_unit_lines);

  snprow = (unsigned char*)malloc(sizeof(unsigned char)*n);
  tped_headers = (char**)malloc(sizeof(char*)*n);
  lbuf = (char*) malloc(sizeof(char*) * SZ_LONG_BUF);

  kin = (double*)calloc(n*n, sizeof(double));
  snpunit = (double*)malloc(n*n_unit_lines * sizeof(double));

  if ( verbose) fprintf(stderr,"Reading TPED file %s.tped ....\n",tpedf);

  tpedh = open_file_with_suffix(tpedf, "tped", 0, 0);
  tpedh.nheadercols = tped_nheadercols;

  nex_autosomal = nex_maf = nex_call_rate = 0;

  for ( i=0, nin=0, nex = 0; tokenize_tped_line_with_col_headers( &tpedh, tped_nheadercols, delims, lbuf, snprow, tped_headers, &nelems, &nmiss) != NULL; ++i) {
    if ( ( verbose ) && ( i % 10000 ) == 0 ) fprintf(stderr,"Reading %d SNPs\n",i);

    if ( ( flag_autosomal == 1 ) && ( ( atoi(tped_headers[0]) == 0 ) || ( atoi(tped_headers[0]) > 22 ) ) ) // if SNP is not in autosomal chromosomes
      //if ( ( flag_autosomal == 1 ) && ( atoi(tped_headers[0]) > 22 ) ) // if SNP is not in autosomal chromosomes
    {
      ++nex; // # excluded snps from the last unit
      ++nex_autosomal;
      continue;
    }

    if ( nelems != n ) {
      emmax_error("Number of values %d in line %d do not match to %d, the number of columns\n", nelems, tpedh.nvaluecols, n);
    }

  /*
    Perform rapid kinship generate IBS, BN, NCOR matrix
--------------------
IBS pairwise matrix
--------------------
i/j  0  1  2  3
0    NA NA NA NA
1    NA 2  1  0
2    NA 1  2  1
3    NA 0  1  2

In fact is what it does is
X = (m*n) genotype matrix (1,2,3 coded)
Xn = X-2
K = (t(Xn) %*% Xn)/(2*m)
------------------
* NA column may be just averaged or predicted based on r2 with previous SNP
  with a certain window size

---------------------
pairwise BN matrix
--------------------
X : (m*n) genotype matrix
Xn : (m*n) matrix each row standardized, missing assigned to 0
K = t(Xn) %*% Xn / L
--------------------
*/

    ac0 = ac1 = ac2 = 0;
    for(j=0; j < n; ++j) {
      if ( snprow[j] > 0 ) {
	if ( snprow[j] < 2 ) {
	  fprintf(stderr,"Invalid snprow[%d] value %d at line %d, individual %d\n",j,(int)snprow[j],i,j);
	}
	snprow[j] -= 2; // from 2,3,4 to 0,1,2 coding

	switch(snprow[j]) {
	case 0:
	  ++ac0;
	  break;
	case 1:
	  ++ac1;
	  break;
	case 2:
	  ++ac2;
	  break;
	default:
	  emmax_error("Unknown allele %s, converted to %d\n",buf,(int)snprow[j]);
	  break;
	}
      }
      else {
	snprow[j] = (unsigned char)NA_GENO_CHAR;
      }
    }

    call_rate = (double)(ac0+ac1+ac2)/(double)n;
    //fprintf(stderr,"CallRate = %lf\n", call_rate);
    if ( call_rate <= call_rate_thres ) {
      ++nex;
      ++nex_call_rate;
      continue;
    }
    aaf = (double)(ac1+2*ac2)/(double)(2*(ac0+ac1+ac2));
    if ( ( aaf <= maf_thres ) || ( 1.-aaf <= maf_thres ) ) {
      ++nex;
      ++nex_maf;
      continue;
    }

    if ( rand_fill_flag == 1 ) {
      for(j=0; j < n; ++j) {
	if ( snprow[j] == (unsigned char)NA_GENO_CHAR ) {
	  if ( (rand() / (double) RAND_MAX) > aaf ) {
	    if ( (rand() / (double) RAND_MAX) > aaf ) {
	      snprow[j] = (unsigned char)2;
	      //ac1 += 2;
	      ++ac2;
	    }
	    else {
	      snprow[j] = (unsigned char)1;
	      ++ac1;
	      //++ac0;
	      //++ac1;
	    }
	  }
	  else {
	    if ( (rand() / (double) RAND_MAX) > aaf ) {
	      snprow[j] = (unsigned char)1;
	      ++ac1;
	      //++ac0;
	      //++ac1;
	    }
	    else {
	      snprow[j] = (unsigned char)0;
	      ++ac0;
	      //ac0 += 2;
	    }
	  }
	}
      }
      aaf = (double)(ac1+2*ac2)/(double)(2*(ac0+ac1+ac2));
      //aaf = (double)ac1/(double)(ac0+ac1);
    }

    // copy current values to arrays
    // Xn = [-1,1] - two rows per SNP : IBS matrix
    // Xn ~ N(0,1) : BN matrix
    if ( ibs_flag == 1 ) {
      for(j=0; j < n; ++j) {
	if ( snprow[j] == (unsigned char)NA_GENO_CHAR ) {
	  snpunit[nin+  j*n_unit_lines] = 2.*aaf-1.;
	  snpunit[nin+1+j*n_unit_lines] = 2.*aaf-1.;
	}
	else {
	  if ( snprow[j] == 0 ) {
	    snpunit[nin+  j*n_unit_lines] = -1.;
	    snpunit[nin+1+j*n_unit_lines] = -1.;
	  }
	  else if ( snprow[j] == 1 ) {
	    snpunit[nin+  j*n_unit_lines] = 1.;
	    snpunit[nin+1+j*n_unit_lines] = -1.;
	  }
	  else if ( snprow[j] == 2 ) {
	    snpunit[nin+  j*n_unit_lines] = 1.;
	    snpunit[nin+1+j*n_unit_lines] = 1.;
	  }
	  else {
	    emmax_error("Invalid genotype %d\n",snprow[j]);
	  }
	}
      }
      nin += 2;
    }
    else {
      for(j=0; j < n; ++j) {
	if ( snprow[j] == (unsigned char)NA_GENO_CHAR ) {
	    snpunit[nin+j*n_unit_lines] = 0.;
	}
	else {
	  snpunit[nin+j*n_unit_lines] = ((double)snprow[j]-(aaf*2.))/sqrt(4*aaf*(1-aaf));
	}
      }
      ++nin;
    }

    // check if nin == n_unit_lines
    if ( nin >= n_unit_lines ) {
      if ( verbose ) {
	fprintf(stderr,"At SNP %d, intermediately computing kinship matrix with %d SNPs..\n",i,n_unit_lines);
      }
      // kin = kin + t(snpunit)%*%(snpunit)
#ifdef INTEL_COMPILER
      dgemm(&ct,&cn,&n,&n,&n_unit_lines,&one,snpunit,&n_unit_lines,snpunit,&n_unit_lines,&one,kin,&n);
#else
      cblas_dgemm(CblasColMajor, CblasTrans, CblasNoTrans, n, n, n_unit_lines, 1.0, snpunit, n_unit_lines, snpunit, n_unit_lines, 1.0, kin, n);
#endif
      memset(snpunit, 0, sizeof(double)*n_unit_lines*n);

      n_sum_nin += nin;
      nin = 0;
    }
  }
  close_file(&tpedh);

  if ( verbose ) {
    fprintf(stderr,"Succesfully finished reading TPED file\n");
  }
  if ( nin > 0 ) {
    if ( verbose ) {
      fprintf(stderr,"Computing kinship matrix with the remaining %d SNPs..\n",nin);
    }
#ifdef INTEL_COMPILER
    dgemm(&ct,&cn,&n,&n,&n_unit_lines,&one,snpunit,&n_unit_lines,snpunit,&n_unit_lines,&one,kin,&n);
#else
    cblas_dgemm(CblasColMajor, CblasTrans, CblasNoTrans, n, n, n_unit_lines, 1.0, snpunit, n_unit_lines, snpunit, n_unit_lines, 1.0, kin, n);
#endif
    n_sum_nin += nin;
  }

  if ( ibs_flag == 1 ) {
    if ( rand_fill_flag == 1 ) {
      suffix = "rIBS.kinf";
    }
    else {
      suffix = "aIBS.kinf";
    }
  }
  else {
    if ( rand_fill_flag == 1 ) {
      suffix = "rBN.kinf";
    }
    else {
      suffix = "aBN.kinf";
    }
  }
  kinsh = open_file_with_suffix( tpedf, suffix, 0, 1 );

  if ( verbose ) fprintf(stderr,"Printing the kinship matrix to file %s.%s\n",tpedf,suffix);

  for(i=0; i < n; ++i) {
    for(j=0; j < n; ++j) {
      if ( j > 0 ) fprintf(kinsh.fp,"\t");
      f = (double)kin[i+j*n]/(double)(n_sum_nin);
      if ( ibs_flag == 1 ) {
	fprintf(kinsh.fp,"%-.*lf",ndigits,0.5*f+0.5);
      }
      else {
	fprintf(kinsh.fp,"%-.*lf",ndigits,f);
      }
    }
    fprintf(kinsh.fp,"\n");
  }
  close_file(&kinsh);
  free(snprow);
  free(kin);
  free(snpunit);
  free(lbuf);
  free(tped_headers);
  free(tfam_headers);
  return 0;
}

max_stars_repo_name: slowkoni/EPI-EMMAX

max_issues_repo_name: slowkoni/EPI-EMMAX

max_forks_repo_name: slowkoni/EPI-EMMAX


text:
#pragma once

#include "CesiumGltfWriter/Library.h"

#include <CesiumJsonWriter/ExtensionWriterContext.h>

#include <gsl/span>

// forward declarations
namespace CesiumGltf {
struct Model;
}

namespace CesiumGltfWriter {

/**
 * @brief The result of writing a glTF with
 * {@link GltfWriter::writeGltf} or {@link GltfWriter::writeGlb}
 */
struct CESIUMGLTFWRITER_API GltfWriterResult {
  /**
   * @brief The final generated std::vector<std::byte> of the glTF or glb.
   */
  std::vector<std::byte> gltfBytes;

  /**
   * @brief Errors, if any, that occurred during the write process.
   */
  std::vector<std::string> errors;

  /**
   * @brief Warnings, if any, that occurred during the write process.
   */
  std::vector<std::string> warnings;
};

/**
 * @brief Options for how to write a glTF.
 */
struct CESIUMGLTFWRITER_API GltfWriterOptions {
  /**
   * @brief If the glTF JSON should be pretty printed. Usable with glTF or GLB
   * (not advised).
   */
  bool prettyPrint = false;

  /**
   * @brief Byte alignment of the GLB binary chunk. When using 64-bit types in
   * EXT_mesh_features or EXT_feature_metadata this value should be set to 8.
   */
  size_t binaryChunkByteAlignment = 4;
};

/**
 * @brief Writes glTF.
 */
class CESIUMGLTFWRITER_API GltfWriter {
public:
  /**
   * @brief Constructs a new instance.
   */
  GltfWriter();

  /**
   * @brief Gets the context used to control how glTF extensions are written.
   */
  CesiumJsonWriter::ExtensionWriterContext& getExtensions();

  /**
   * @brief Gets the context used to control how glTF extensions are written.
   */
  const CesiumJsonWriter::ExtensionWriterContext& getExtensions() const;

  /**
   * @brief Serializes the provided model into a glTF JSON byte vector.
   *
   * @details Ignores internal data such as {@link CesiumGltf::BufferCesium}
   * and {@link CesiumGltf::ImageCesium} when serializing the glTF. Internal
   * data must either be converted to data uris or saved as external files. The
   * buffer.uri and image.uri fields must be set accordingly prior to calling
   * this function.
   *
   * @param model The model.
   * @param options Options for how to write the glTF.
   * @return The result of writing the glTF.
   */
  GltfWriterResult writeGltf(
      const CesiumGltf::Model& model,
      const GltfWriterOptions& options = GltfWriterOptions()) const;

  /**
   * @brief Serializes the provided model into a glb byte vector.
   *
   * @details The first buffer object implicitly refers to the GLB binary chunk
   * and should not have a uri. Ignores internal data such as
   * {@link CesiumGltf::BufferCesium} and {@link CesiumGltf::ImageCesium}.
   *
   * @param model The model.
   * @param bufferData The buffer data to store in the GLB binary chunk.
   * @param options Options for how to write the glb.
   * @return The result of writing the glb.
   */
  GltfWriterResult writeGlb(
      const CesiumGltf::Model& model,
      const gsl::span<const std::byte>& bufferData,
      const GltfWriterOptions& options = GltfWriterOptions()) const;

private:
  CesiumJsonWriter::ExtensionWriterContext _context;
};

} // namespace CesiumGltfWriter

ext: h

max_stars_repo_name: yieryi/cesium-native

max_issues_repo_name: yieryi/cesium-native

max_forks_repo_name: yieryi/cesium-native


text:
/**
 *
 * @file
 *
 *  PLASMA is a software package provided by:
 *  University of Tennessee, US,
 *  University of Manchester, UK.
 *
 **/
#ifndef PLASMA_CORE_LAPACK_H
#define PLASMA_CORE_LAPACK_H

#if defined(HAVE_MKL) || defined(PLASMA_WITH_MKL)
    #define MKL_Complex16 double _Complex
    #define MKL_Complex8  float _Complex

    #include <mkl_cblas.h>
    #include <mkl_lapacke.h>

    // MKL LAPACKE doesn't provide LAPACK_GLOBAL macro, so define it here.
    // MKL provides all 3 name manglings (foo, foo_, FOO); pick foo_.
    #ifndef LAPACK_GLOBAL
    #define LAPACK_GLOBAL(lcname,UCNAME)  lcname##_
    #endif
#elif defined(HAVE_ESSL) || defined(PLASMA_WITH_ESSL)
    // GCC + ESSL(BLAS) + LAPACKE/CBLAS from LAPACK
    #include <cblas.h>
    #include <lapacke.h>

    #ifndef LAPACK_GLOBAL
    #define LAPACK_GLOBAL(lcname,UCNAME) lcname##_
    #endif
#else
    #include <cblas.h>
    #include <lapacke.h>

    // Original  cblas.h does: enum CBLAS_ORDER {...};
    // Intel mkl_cblas.h does: typedef enum {...} CBLAS_ORDER;
    // LAPACK    cblas.h does: typedef enum {...} CBLAS_ORDER;
    // OpenBLAS  cblas.h does: typedef enum CBLAS_ORDER {...} CBLAS_ORDER;
    // We use (CBLAS_ORDER), so add these typedefs for original cblas.h
    #ifdef CBLAS_ADD_TYPEDEF
    typedef enum CBLAS_ORDER CBLAS_ORDER;
    typedef enum CBLAS_TRANSPOSE CBLAS_TRANSPOSE;
    typedef enum CBLAS_UPLO CBLAS_UPLO;
    typedef enum CBLAS_DIAG CBLAS_DIAG;
    typedef enum CBLAS_SIDE CBLAS_SIDE;
    #endif
#endif

#ifndef lapack_int
#define lapack_int int
#endif

#include "core_lapack_s.h"
#include "core_lapack_d.h"
#include "core_lapack_c.h"
#include "core_lapack_z.h"

#endif // PLASMA_CORE_LAPACK_H


text:
#include <time.h>
#include <math.h>
#include <stdio.h>
#include <string.h>
#include <unistd.h>
#include <getopt.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_vector.h>
#include "../include/type.h"
#include "../include/util.h"
#include "../include/myfunc.h"
#include "../include/global.h"

double rho = 0.9;
double cutoff = 0.0001;
clock_t dur = 0;
// char *optarg;
// int optind = 0;

int main(int argc, char *argv[]) {
    char *inputf = "input.txt", *outputf = (char*)malloc(sizeof(char)*MAX_CHAR);
    char str_line[MAX_LINE], *title_element_list[100];
    int nthread = 1, batch_size = 1, opt, row_num = 0, i = 0;

    while ((opt = getopt(argc, argv, "t:bc:o:i:")) != -1) {
        switch (opt) {
        case 't':
            nthread = atoi(optarg);
            break;
        case 'b':
            batch_size = 0;
            break;
        case 'c':
            cutoff = atof(optarg);
            break;
        case 'o':
            // strcpy(outputf, optarg);
            // strcat(outputf, "/rMATS_Result_P.txt");
            outputf = optarg;
            break;
        case 'i':
            inputf = optarg;
            break;
        default: /* '?' */
            fprintf(stderr, "Usage: %s [-t nthread] [-c] cutoff [-o] output folder [-i] input file\n", argv[0]);
            exit(EXIT_FAILURE);
        }
    }

    printf("number of thread=%d; input file=%s; output folder=%s; cutoff=%g;\n",\
           nthread, inputf, outputf, cutoff);

    if (optind > argc) {
        fprintf(stderr, "Expected argument after options\n");
        exit(EXIT_FAILURE);
    }

    if (nthread > MAXTHREAD) {
        printf("nthread (%d) should be less than MAXTHREAD (%d)\n", nthread, MAXTHREAD);
        exit(0);
    }

    diff_list_node *list = (diff_list_node*)malloc(sizeof(diff_list_node));
    list->data = NULL, list->next = NULL, list->end = list;
    row_num = parse_file(inputf, list, title_element_list);

    clock_t start, diff;
    start = clock();

    odiff *datum[row_num];
    diff_list_node *node = list;
    double *prob[row_num];
    for (i = 0; i < row_num; ++i) {
        node = node->next;
        datum[i] = node->data;
    }
    mp_threadpool(nthread, row_num, thread_wrapper_for_LT, (void**)datum, (void**)prob);

    diff = clock() - start;
    int msec = diff * 1000 / CLOCKS_PER_SEC;
    printf("Total Wallclock time taken %d seconds %d milliseconds\n", msec/1000, msec%1000);
    printf("Wallclock time per thread taken %d seconds %d milliseconds\n", msec/1000/nthread, (msec%1000)/nthread);

    FILE *ifp = NULL, *ofp = NULL;
    if ((ifp = fopen(inputf, "r")) == NULL) {
        printf("Fail to open %s!", inputf);
        return -1;
    }
    if ((ofp = fopen(outputf, "w")) == NULL) {
        printf("Fail to open %s!", outputf);
        return -1;
    }

    // original comment: analyze the title of the inputed data file to find the
    //                   information of how much simulation are involved
    fgets(str_line, MAX_LINE, ifp);
    str_line[strlen(str_line)-1] = 0;
    fprintf(ofp, "%s\tPValue\n", str_line);
    for (i = 0; i < row_num; ++i) {
        fgets(str_line, MAX_LINE, ifp);
        str_line[strlen(str_line)-1] = 0;
        fprintf(ofp, "%s\t%.12g\n", str_line, *prob[i]);
    }
    fclose(ifp);
    fclose(ofp);

    msec = dur * 1000 / CLOCKS_PER_SEC;
    printf("Time for func(single thread): %d seconds %d milliseconds\n", msec/1000, msec%1000);

    return 0;
}

ext: c


text:
#ifndef COVAR_H_
#define COVAR_H_

#include <stdio.h>
#include <gsl/gsl_statistics_double.h>
#include <gsl/gsl_matrix_double.h>
#include <vector>

// \TODO: trim excess comments
// \TODO: split this into an implementation file

/*
 *  Okay, so here's what we need
 *  1) A data structure to store the data points (vector?)
 *  2) A way of packaging the data points (struct of floats?)
 *  3) A method to push new data points onto the data structure
 *
 *  So the above is an API towards collecting data. On the ROS side,
 *  1) As new messages from the ROS topic arrives, we need the
 *     data_structure.push() to happen
 *
 *  Next, we need a method for generating the covariance matrix. This
 *  should happen after the data collection period has expired. What needs
 *  to happen when we generate the covariance matrix?
 *  1) Convert the data_structure into a gsl_matrix
 *  2) Don't forget to allocate another gsl_matrix which will hold the
 *     resulting covariance matrix.
 *  3) Double for-loop to generate the variance-covariance value for each
 *     element of the matrix
 */

// the covariance matrix of our odometry topic tracks six variables
// as in, think of each data point as a six-dimensional point with the
// following axes: x, y, z, roll, pitch, yaw
#define ODOM_COVAR_MAT_VARSIZE 6

// struct odom_struct
// {
//     double x;
//     double y;
//     double z;
//     double roll;
//     double pitch;
//     double yaw;
// };

// Since we're traversing through the data points using a loop, maybe it'd
// be better to use an array of doubles rather than a struct
// double * data_point = new data_point[6];
// We establish the convention here:
// index: 0  1  2  3    4      5
// axis:  x, y, z, roll, pitch, yaw

// However, from a maintainability perspective, relying on a double
// pointer is bad - there is no way of guaranteeing the fixed six-element
// size that we want for data point
// \TODO: create a better container (i.e. class) for each data point that
//        has as a member method
//        1) the square bracket operator, for fetching an element

class covar
{
    private:
        gsl_matrix *covar_mat;
        std::vector<double*> *data_points;

    public:
        covar()
        {
            data_points = new std::vector<double*>;
            covar_mat = gsl_matrix_alloc(ODOM_COVAR_MAT_VARSIZE,
                                         ODOM_COVAR_MAT_VARSIZE);
        }

        ~covar()
        {
            delete covar_mat;
            covar_mat = 0;
            delete data_points;
            data_points = 0;
        }

        // We expect the input to be a double array with six elements
        bool insert(double* data_point)
        {
            // double * temp = new double[6];
            data_points->push_back(data_point);
            return true;
        }

        void print_dp()
        {
            for (int i = 0; i < data_points->size(); i++)
            {
                printf("%i: %f %f %f %f %f %f\n",
                            i,
                            (*data_points)[i][0],
                            (*data_points)[i][1],
                            (*data_points)[i][2],
                            (*data_points)[i][3],
                            (*data_points)[i][4],
                            (*data_points)[i][5]);
            }
        }

        void print_mat()
        {
            for (int i = 0; i < covar_mat->size1; i++)
            {
                printf("\n");
                for (int j = 0; j < covar_mat->size2; j++)
                {
                    printf("%f", covar_mat->data[covar_mat->size2*i+j]);
                    if(j < covar_mat->size2 - 1) printf(", ");
                }
            }
            printf("\n");
        }

        // \TODO: if covar_mat not available, deny get_mat
        bool get_mat(double* array)
        {
            for(int i = 0; i < 36; i++)
                array[i] = gsl_matrix_get(covar_mat,
                                          i/6,  //row
                                          i%6); //column

            return true;
        }

        bool generate_covar()
        {
            /* no point in generating a covariance matrix if we don't have
             * the data points to do so
             */
            int i, j;
            int rsize = data_points->size();

            if(rsize <= 1) return false;

            gsl_vector_view a, b;
            gsl_matrix *A;
            A = gsl_matrix_alloc(rsize, ODOM_COVAR_MAT_VARSIZE);
            /* the covariance matrix should already be instantiated by the
             * class constructor
             */
            // covar_mat = gsl_matrix_alloc(ODOM_COVAR_MAT_VARSIZE,
            //                              ODOM_COVAR_MAT_VARSIZE);

            // convert data_points into gsl_matrix A
            for ( i = 0; i < rsize; i++)
                for ( j = 0; j < ODOM_COVAR_MAT_VARSIZE; j++)
                    gsl_matrix_set (A, i, j, (*data_points)[i][j]);

            // generate the var-covar values based on A
            // insert these values into covar_mat
            for (i = 0; i < A->size2; i++) {
                for (j = 0; j < A->size2; j++) {
                a = gsl_matrix_column (A, i);
                b = gsl_matrix_column (A, j);
                double cov = gsl_stats_covariance(a.vector.data,
                                                    a.vector.stride,
                                                    b.vector.data,
                                                    b.vector.stride,
                                                    rsize);
                gsl_matrix_set (covar_mat, i, j, cov);
                }
            }

            return true;
        }
};

#endif // COVAR_H_


text:
#include <gsl/gsl_sf_erf.h>
#include <gsl/gsl_vector.h>
#include <math.h>
#include <time.h>
#include <stdio.h>
#include <stdlib.h>
//#include <gsl/gsl_ran_gaussian.h>

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

    double a = gsl_sf_erf(1);
    printf("The random number we get is %f.", a);
    exit(0);
}

max_stars_repo_name: ShreyaR/WorkerPoolSelection

max_issues_repo_name: ShreyaR/WorkerPoolSelection

max_forks_repo_name: ShreyaR/WorkerPoolSelection


text:
#include "transform_animation.h"
#include "entities/entities.h"
#include <gsl/gsl_poly.h>
#include <gsl/gsl_math.h>
#include "utility/logging.h"

// from x^0 ... x^(len)
static const double drop_in_coeffients[] = { 0.0, -2.0, -3.0, 1.0 };

static void no_op_release(void *data) { }
static void simple_release(void *data) { if (data != NULL) { free(data); } }

typedef struct {
  double duration;
  double scale;
  uint8_t flags;
} drop_in_data;

static void drop_in_get_transform(
    double t, ANIMATION_PLAYBACK_MODE mode, int reverse, void *pdata, double *m_out) {

  drop_in_data *data = pdata;

  double clamped_t;

  if (mode == PLAY_ONCE) {
    if (reverse) {
      clamped_t = GSL_MAX(t, -data->duration);
      clamped_t = data->duration + clamped_t;
    } else {
      clamped_t = GSL_MIN(t, data->duration);
    }
  } else if (mode == LOOP) {
    double full_iterations = floor(t / data->duration);
    clamped_t = t - (full_iterations * data->duration);
  } else if (mode == OSCILLATE) {
    int full_iterations = floor(t / data->duration);
    if (GSL_IS_EVEN(full_iterations)) { // same as loop
      clamped_t = t - (full_iterations * data->duration);
    } else { // invert it
      clamped_t = t - (full_iterations * data->duration);
      clamped_t = data->duration - clamped_t;
    }
  }

  uint8_t flags = data->flags;
  double translation[4];
  double value = data->scale * gsl_poly_eval(drop_in_coeffients, 4, clamped_t - 0.6);

  gsl_matrix_view out_view = gsl_matrix_view_array(m_out, 4, 4);
  gsl_matrix_set_identity(&out_view.matrix);
  gsl_vector_view tview = gsl_matrix_column(&out_view.matrix, 3);
  if (flags & XFORM_TRANSLATION_X) { gsl_vector_set(&tview.vector, 0, value); }
  if (flags & XFORM_TRANSLATION_Y) { gsl_vector_set(&tview.vector, 1, value); }
  if (flags & XFORM_TRANSLATION_Z) { gsl_vector_set(&tview.vector, 2, value); }
}

animation_channel build_drop_in_channel(uint8_t apply_to, double scale) {
  drop_in_data *data = calloc(1, sizeof(drop_in_data));
  data->duration = 5;
  data->scale = scale;
  data->flags = apply_to;
  animation_channel retval = {
    .data=data,
    .get_transform=drop_in_get_transform,
    .release_resources=simple_release
  };
  return retval;
}

static const uint8_t CONSTANT_TRANSLATION = (1<<0);
static const uint8_t CONSTANT_X_AXIS =      (1<<1);
static const uint8_t CONSTANT_Y_AXIS =      (1<<2);
static const uint8_t CONSTANT_Z_AXIS =      (1<<3);

typedef struct {
  gsl_matrix *basis;
  uint8_t constant_component_flags;
  double *frame_data;
  size_t entries_per_frame;
  double duration; // seconds
} optimized_keyframe_data;

static void optimized_keyframe_get_transform(
    double t, ANIMATION_PLAYBACK_MODE mode, int reverse, void *data, double *m_out) {
  // look through the collection of keyframes and determine indices for
  // the left and right frames bounding the provided `t`.
  // doublesPerFrame, flat array for frame data store
  // pack into (time, then optional translation, basis vectors, etc)
}

animation_channel build_optimized_keyframe_channel(
    keyframe *keyframes, size_t nframes) {
  animation_channel retval = {
    .data=NULL,
    .get_transform=optimized_keyframe_get_transform,
    .release_resources=simple_release
  };
  // TODO
  return retval;
}

void animate_entity(double deltaT, entity_info* e) {
  animation_state *state = &e->animation_state;

  // update time trackers
  double scaledT = deltaT * state->playback_speed;
  state->current_cycle_time += scaledT;
  state->time_in_transition += scaledT;
  state->transition_cycle_time += scaledT;
}

transform_animation_type const transform_animation = {
  .build_optimized_keyframe_channel=build_optimized_keyframe_channel,
  .build_drop_in_channel=build_drop_in_channel,
  .animate_entity=animate_entity
};


text:
#include <stdlib.h>            // malloc
#include <stdio.h>             // printf
#include <math.h>              // fabs, sqrt, etc.
#include <time.h>              // time
#include <unistd.h>            // getpid
#include <gsl/gsl_rng.h>       // GNU Scientific Library
#include <gsl/gsl_cdf.h>       // GNU Scientific Library
#include <gsl/gsl_randist.h>   // GNU Scientific Library

double dvonmises(double* x, double mu, double nu,
               double kappa1, double kappa2, double lambda) {
  double phi = x[0];
  double psi = x[1];
  double kernal = exp( kappa1*cos(phi - mu) + kappa2*cos(psi-nu) +
  lambda*sin(phi-mu)*sin(psi-nu));
  return kernal;
}

void rvonmises2(int* ptrn, double* ptrmu, double* ptrnu,
              double* ptrkappa1, double* ptrkappa2, double* ptrlambda,
              double* phiout, double* psiout)
{
  // Assign the variables
  int n = *ptrn;
  double mu = *ptrmu;
  double nu = *ptrnu;
  double kappa1 = *ptrkappa1;
  double kappa2 = *ptrkappa2;
  double lambda = *ptrlambda;
  double modeh;
  // If Unimodal...
   if(lambda*lambda<(kappa1*kappa2)) {
     double x[] = {mu, nu};
     modeh = dvonmises(x,mu,nu,kappa1,kappa2,lambda);
   }
   //If Bimodal...
   if(lambda*lambda>(kappa1*kappa2)) {
     double psi =
     acos(kappa2/fabs(lambda)*sqrt((lambda*lambda+kappa1*kappa1)/(lambda*lambda+kappa2*kappa2)));
     double phi =
     acos(kappa1/fabs(lambda)*sqrt((lambda*lambda+kappa2*kappa2)/(lambda*lambda+kappa1*kappa1)));
     //Caclulate the 4 possible maximum if lambda >0
     if(lambda>0.0){
       double x1[] = {mu+psi , nu+phi};
       double x2[] = {mu-psi , nu-phi};
       double mode1 = dvonmises(x1,mu,nu,kappa1,kappa2,lambda);
       double mode2 = dvonmises(x2,mu,nu,kappa1,kappa2,lambda);
       modeh = fmax(mode1,mode2);
     }
     if(lambda<0.0){
       double x1[] = {mu-phi , nu+psi};
       double x2[] = {mu+psi , nu-psi};
       double mode1 = dvonmises(x1,mu,nu,kappa1,kappa2,lambda);
       double mode2 = dvonmises(x2,mu,nu,kappa1,kappa2,lambda);
       modeh = fmax(mode1,mode2);
     }
   }
  //implement the rejection sampler now that we can use the mode to get alpha
  double pi = acos(-1);
  double height_of_sample_density = 1/(4*pi*pi);
  double alpha = height_of_sample_density/modeh;
  //printf("%s","mode height\n");
  //printf("%f\n",modeh);

  gsl_rng* r = gsl_rng_alloc (gsl_rng_taus);
  int seed = time(NULL)*getpid();
  gsl_rng_set(r,seed);

  int i = 0;
  while( i < n){
    //sample from the sample envelope: the uniform
    double two_uniforms[0];
    //evaluate the density at those uniforms
    two_uniforms[0] = gsl_ran_flat(r, -pi, pi);
    two_uniforms[1] = gsl_ran_flat(r, -pi, pi);
    double f_x = dvonmises(two_uniforms,mu,nu,kappa1,kappa2,lambda);
    double e_x = height_of_sample_density/alpha;

    // Test if U <= f_x/e_x and store the uniforms if they work
    double U = gsl_ran_flat(r,0,1);
    if( U <= f_x/e_x){
      phiout[i] = two_uniforms[0];
      psiout[i] = two_uniforms[1];
      i +=1;
    }
  }
}

int main(int argc, char* argv[]) {
  if ( argc != 7 ) {
    fprintf(stderr,"usage: %s nreps mu nu kappa1 kappa2 lambda\n",argv[0]);

    return 1;
  }
  if (argv[7]=0){
    fprintf(stderr,"lambda cannot be 0");
    return 1;
  }
  // Simulation parameters
  int n = atoi(argv[1]);
  double mu = atof(argv[2]);
  double nu = atof(argv[3]);
  double kappa1 = atof(argv[4]);
  double kappa2 = atof(argv[5]);
  double lambda = atof(argv[6]);
//because our rvonmises funciton is void, we must allocate memor to psi and phi
//so that we can get results and call them later
  double* phiout = (double*) malloc(n*sizeof(double));
  double* psiout = (double*) malloc(n*sizeof(double));
  rvonmises2(&n,&mu,&nu,&kappa1,&kappa2,&lambda,phiout,psiout);
  printf("%s %s\n","phi","psi");
  for(int i = 0; i < n; i++){
    printf("%f %f\n",phiout[i],psiout[i]);
  }
  free(phiout);
  free(psiout);
  return 0;//to get it to work, comment this out
}

max_stars_repo_name: bryanwhiting/randomCsampler

max_issues_repo_name: bryanwhiting/randomCsampler

max_forks_repo_name: bryanwhiting/randomCsampler


text:
#include <math.h>
#include <stdlib.h>
#if !defined(__APPLE__)
#include <malloc.h>
#endif
#include <stdio.h>
#include <assert.h>
#include <time.h>
#include <string.h>

#include <fftw3.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_erf.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_legendre.h>
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_sf_expint.h>
#include <gsl/gsl_deriv.h>
#include <gsl/gsl_interp2d.h>
#include <gsl/gsl_spline2d.h>

#include "../cosmolike_core/theory/basics.c"
#include "../cosmolike_core/theory/structs.c"
#include "../cosmolike_core/theory/parameters.c"
#include "../cosmolike_core/emu17/P_cb/emu.c"
#include "../cosmolike_core/theory/recompute.c"
#include "../cosmolike_core/theory/cosmo3D.c"
#include "../cosmolike_core/theory/redshift_spline.c"
#include "../cosmolike_core/theory/halo.c"
#include "../cosmolike_core/theory/HOD.c"
#include "../cosmolike_core/theory/pt.c"
#include "../cosmolike_core/theory/cosmo2D_fourier.c"
#include "../cosmolike_core/theory/IA.c"
#include "../cosmolike_core/theory/cluster.c"
#include "../cosmolike_core/theory/BAO.c"
#include "../cosmolike_core/theory/external_prior.c"
#include "../cosmolike_core/theory/GRS.c"
#include "init_WFIRST_forecasts.c"
#include "like_grs.c"

double C_shear_tomo_sys(double ell,int z1,int z2);
double C_cgl_tomo_sys(double ell_Cluster,int zl,int nN, int zs);
double C_gl_tomo_sys(double ell,int zl,int zs);
void set_data_shear(int Ncl, double *ell, double *data, int start);
void set_data_ggl(int Ncl, double *ell, double *data, int start);
void set_data_clustering(int Ncl, double *ell, double *data, int start);
void set_data_cluster_N(double *data, int start);
void set_data_cgl(double *ell_Cluster, double *data, int start);
double log_L_3x2pt_clusterN_clusterWL_GRS_SN();
double log_L_3x2pt_clusterN_clusterWL_GRS();
void compute_data_vector(char *details, double OMM, double S8, double NS, double W0,double WA, double OMB, double H0, double MGSigma, double MGmu, double B1, double B2, double B3, double B4,double B5, double B6, double B7, double B8, double B9, double B10, double SP1, double SP2, double SP3, double SP4, double SP5, double SP6, double SP7, double SP8, double SP9, double SP10, double SPS1, double CP1, double CP2, double CP3, double CP4, double CP5, double CP6, double CP7, double CP8, double CP9, double CP10, double CPS1, double M1, double M2, double M3, double M4, double M5, double M6, double M7, double M8, double M9, double M10, double A_ia, double beta_ia, double eta_ia, double eta_ia_highz, double LF_alpha, double LF_P, double LF_Q, double LF_red_alpha, double LF_red_P, double LF_red_Q,double mass_obs_norm, double mass_obs_slope, double mass_z_slope, double mass_obs_scatter_norm, double mass_obs_scatter_mass_slope, double mass_obs_scatter_z_slope);
double log_multi_like(double OMM, double S8, double NS, double W0,double WA, double OMB, double H0, double MGSigma, double MGmu, double B1, double B2, double B3, double B4,double B5, double B6, double B7, double B8, double B9, double B10, double SP1, double SP2, double SP3, double SP4, double SP5, double SP6, double SP7, double SP8, double SP9, double SP10, double SPS1, double CP1, double CP2, double CP3, double CP4, double CP5, double CP6, double CP7, double CP8, double CP9, double CP10, double CPS1, double M1, double M2, double M3, double M4, double M5, double M6, double M7, double M8, double M9, double M10, double A_ia, double beta_ia, double eta_ia, double eta_ia_highz, double LF_alpha, double LF_P, double LF_Q, double LF_red_alpha, double LF_red_P, double LF_red_Q,double mass_obs_norm, double mass_obs_slope, double mass_z_slope, double mass_obs_scatter_norm, double mass_obs_scatter_mass_slope, double mass_obs_scatter_z_slope, double GRSB1, double GRSB2, double GRSB3, double GRSB4, double GRSB5, double GRSB6, double GRSB7, double SIGMAP1, double SIGMAP2, double SIGMAP3, double SIGMAP4, double SIGMAP5, double SIGMAP6, double SIGMAP7,double SIGMAZ, double PSHOT, double KSTAR);
void write_vector_wrapper(char *details, input_cosmo_params ic, input_nuisance_params in);
double log_like_wrapper(input_cosmo_params ic, input_nuisance_params in, input_nuisance_params_grs ingr);
int get_N_tomo_shear(void);
int get_N_tomo_clustering(void);
int get_N_ggl(void);
int get_N_ell(void);

int get_N_tomo_shear(void){
  return tomo.shear_Nbin;
}
int get_N_tomo_clustering(void){
  return tomo.clustering_Nbin;
}
int get_N_ggl(void){
  return tomo.ggl_Npowerspectra;
}
int get_N_ell(void){
  return like.Ncl;
}

double C_shear_tomo_sys(double ell, int z1, int z2)
{
  double C;
  // C= C_shear_tomo_nointerp(ell,z1,z2);
  // if(like.IA==1) C+=C_II_nointerp(ell,z1,z2)+C_GI_nointerp(ell,z1,z2);

  if(like.IA!=1) C= C_shear_tomo_nointerp(ell,z1,z2);
  //if(like.IA==1) C= C_shear_shear_IA(ell,z1,z2);
  if(like.IA==1) C = C_shear_tomo_nointerp(ell,z1,z2)+C_II_nointerp(ell,z1,z2)+C_GI_nointerp(ell,z1,z2);
  if(like.IA==2) C += C_II_lin_nointerp(ell,z1,z2)+C_GI_lin_nointerp(ell,z1,z2);
  if(like.shearcalib==1) C *=(1.0+nuisance.shear_calibration_m[z1])*(1.0+nuisance.shear_calibration_m[z2]);
  //printf("%le %d %d %le\n",ell,z1,z2,C_shear_tomo_nointerp(ell,z1,z2)+C_II_JB_nointerp(ell,z1,z2)+C_GI_JB_nointerp(ell,z1,z2));
return C;
}

double C_gl_tomo_sys(double ell,int zl,int zs)
{
  double C;
  // C=C_gl_tomo_nointerp(ell,zl,zs);
  // if(like.IA==1) C += C_gI_nointerp(ell,zl,zs);

  if(like.IA!=1) C=C_gl_tomo_nointerp(ell,zl,zs);
  if(like.IA==1) C = C_ggl_IA(ell,zl,zs);
  if(like.IA==2) C += C_gI_lin_nointerp(ell,zl,zs);
  if(like.shearcalib==1) C *=(1.0+nuisance.shear_calibration_m[zs]);
return C;
}

double C_cgl_tomo_sys(double ell_Cluster, int zl,int nN, int zs)
{
  double C;
  C=C_cgl_tomo_nointerp(ell_Cluster,zl,nN,zs);
  //if(like.IA!=0) C +=
  if(like.shearcalib==1) C *=(1.0+nuisance.shear_calibration_m[zs]);
return C;
}

void set_data_shear(int Ncl, double *ell, double *data, int start)
{
  int i,z1,z2,nz;
  double a;
  for (nz = 0; nz < tomo.shear_Npowerspectra; nz++){
    z1 = Z1(nz); z2 = Z2(nz);
    for (i = 0; i < Ncl; i++){
      if (ell[i] < like.lmax_shear){ data[Ncl*nz+i] = C_shear_tomo_sys(ell[i],z1,z2);}
      else {data[Ncl*nz+i] = 0.;}
    }
  }
}

void set_data_ggl(int Ncl, double *ell, double *data, int start)
{
  int i, zl,zs,nz;
  for (nz = 0; nz < tomo.ggl_Npowerspectra; nz++){
    zl = ZL(nz); zs = ZS(nz);
    for (i = 0; i < Ncl; i++){
      if (test_kmax(ell[i],zl)){
        data[start+(Ncl*nz)+i] = C_gl_tomo_sys(ell[i],zl,zs);
      }
      else{
        data[start+(Ncl*nz)+i] = 0.;
      }
    }
  }
}

void set_data_clustering(int Ncl, double *ell, double *data, int start){
  int i, nz;
  for (nz = 0; nz < tomo.clustering_Npowerspectra; nz++){
    //printf("%d %e %e\n",nz, gbias.b[nz][1],pf_photoz(gbias.b[nz][1],nz));
    for (i = 0; i < Ncl; i++){
      if (test_kmax(ell[i],nz)){data[start+(Ncl*nz)+i] = C_cl_tomo_nointerp(ell[i],nz,nz);}
      else{data[start+(Ncl*nz)+i] = 0.;}
      //printf("%d %d %le %le\n",nz,nz,ell[i],data[Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra + nz)+i]);
    }
  }
}

void set_data_cluster_N(double *data, int start){
  int nN, nz;
  for (nz = 0; nz < tomo.cluster_Nbin; nz++){
    for (nN = 0; nN < Cluster.N200_Nbin; nN++){
      data[start+Cluster.N200_Nbin*nz+nN] = N_N200(nz, nN);
    }
  }
}

void set_data_cgl(double *ell_Cluster, double *data, int start)
{
  int zl,zs,nN,nz,i,j;
  for(nN = 0; nN < Cluster.N200_Nbin; nN++){
    for (nz = 0; nz < tomo.cgl_Npowerspectra; nz++){
      zl = ZC(nz); zs = ZSC(nz);
      for (i = 0; i < Cluster.lbin; i++){
        j = start;
        j += (nz*Cluster.N200_Nbin+nN)*Cluster.lbin +i;
        data[j] = C_cgl_tomo_sys(ell_Cluster[i],zl,nN,zs);
      }
    }
  }
}

int set_cosmology_params(double OMM, double S8, double NS, double W0,double WA, double OMB, double H0, double MGSigma, double MGmu)
{
  cosmology.Omega_m=OMM;
  cosmology.Omega_v= 1.0-cosmology.Omega_m;
  cosmology.sigma_8=S8;
  cosmology.n_spec= NS;
  cosmology.w0=W0;
  cosmology.wa=WA;
  cosmology.omb=OMB;
  cosmology.h0=H0;
  cosmology.MGSigma=MGSigma;
  cosmology.MGmu=MGmu;

  if (cosmology.Omega_m < 0.05 || cosmology.Omega_m > 0.6) return 0;
  if (cosmology.omb < 0.04 || cosmology.omb > 0.055) return 0;
  if (cosmology.sigma_8 < 0.5 || cosmology.sigma_8 > 1.1) return 0;
  if (cosmology.n_spec < 0.84 || cosmology.n_spec > 1.06) return 0;
  if (cosmology.w0 < -2.1 || cosmology.w0 > -0.0) return 0;
  if (cosmology.wa < -2.6 || cosmology.wa > 2.6) return 0;
  if (cosmology.h0 < 0.4 || cosmology.h0 > 0.9) return 0;
  //CH BEGINS
  //CH: to use for running planck15_BA0_w0_wa prior alone)
  //printf("like_fourier.c from WFIRST_forecasts: cosmology bounds set for running with planck15_BA0_w0_wa prior\n");
  //if (cosmology.Omega_m < 0.05 || cosmology.Omega_m > 0.6) return 0;
  //if (cosmology.omb < 0.01 || cosmology.omb > 0.1) return 0;
  //if (cosmology.sigma_8 < 0.5 || cosmology.sigma_8 > 1.1) return 0;
  //if (cosmology.n_spec < 0.84 || cosmology.n_spec > 1.06) return 0;
  //if (cosmology.w0 < -2.1 || cosmology.w0 > 1.5) return 0;
  //if (cosmology.wa < -5.0 || cosmology.wa > 2.6) return 0;
  //if (cosmology.h0 < 0.3 || cosmology.h0 > 0.9) return 0;
  //CH ENDS
  return 1;
}

void set_nuisance_shear_calib(double M1, double M2, double M3, double M4, double M5, double M6, double M7, double M8, double M9, double M10)
{
  nuisance.shear_calibration_m[0] = M1;
  nuisance.shear_calibration_m[1] = M2;
  nuisance.shear_calibration_m[2] = M3;
  nuisance.shear_calibration_m[3] = M4;
  nuisance.shear_calibration_m[4] = M5;
  nuisance.shear_calibration_m[5] = M6;
  nuisance.shear_calibration_m[6] = M7;
  nuisance.shear_calibration_m[7] = M8;
  nuisance.shear_calibration_m[8] = M9;
  nuisance.shear_calibration_m[9] = M10;
}

int set_nuisance_shear_photoz(double SP1,double SP2,double SP3,double SP4,double SP5,double SP6,double SP7,double SP8,double SP9,double SP10,double SPS1)
{
  int i;
  nuisance.bias_zphot_shear[0]=SP1;
  nuisance.bias_zphot_shear[1]=SP2;
  nuisance.bias_zphot_shear[2]=SP3;
  nuisance.bias_zphot_shear[3]=SP4;
  nuisance.bias_zphot_shear[4]=SP5;
  nuisance.bias_zphot_shear[5]=SP6;
  nuisance.bias_zphot_shear[6]=SP7;
  nuisance.bias_zphot_shear[7]=SP8;
  nuisance.bias_zphot_shear[8]=SP9;
  nuisance.bias_zphot_shear[9]=SP10;

  for (i=0;i<tomo.shear_Nbin; i++){
    nuisance.sigma_zphot_shear[i]=SPS1;
    if (nuisance.sigma_zphot_shear[i]<0.001) return 0;
  }
  return 1;
}

int set_nuisance_clustering_photoz(double CP1,double CP2,double CP3,double CP4,double CP5,double CP6,double CP7,double CP8,double CP9,double CP10,double CPS1)
{
  int i;
  nuisance.bias_zphot_clustering[0]=CP1;
  nuisance.bias_zphot_clustering[1]=CP2;
  nuisance.bias_zphot_clustering[2]=CP3;
  nuisance.bias_zphot_clustering[3]=CP4;
  nuisance.bias_zphot_clustering[4]=CP5;
  nuisance.bias_zphot_clustering[5]=CP6;
  nuisance.bias_zphot_clustering[6]=CP7;
  nuisance.bias_zphot_clustering[7]=CP8;
  nuisance.bias_zphot_clustering[8]=CP9;
  nuisance.bias_zphot_clustering[9]=CP10;

  for (i=0;i<tomo.clustering_Nbin; i++){
    nuisance.sigma_zphot_clustering[i]=CPS1;
    if (nuisance.sigma_zphot_clustering[i]<0.001) return 0;
  }
  return 1;
}

int set_nuisance_ia(double A_ia, double beta_ia, double eta_ia, double eta_ia_highz, double LF_alpha, double LF_P, double LF_Q, double LF_red_alpha, double LF_red_P, double LF_red_Q)
{
  nuisance.A_ia=A_ia;
  nuisance.beta_ia=beta_ia;
  nuisance.eta_ia=eta_ia;
  nuisance.eta_ia_highz=eta_ia_highz;
  nuisance.LF_alpha=LF_alpha;
  nuisance.LF_P=LF_P;
  nuisance.LF_Q=LF_Q;
  nuisance.LF_red_alpha=LF_red_alpha;
  nuisance.LF_red_P=LF_red_P;
  nuisance.LF_red_Q=LF_red_Q;
  if (nuisance.A_ia < 0.0 || nuisance.A_ia > 10.0) return 0;
  if (nuisance.beta_ia < -1.0 || nuisance.beta_ia > 3.0) return 0;
  if (nuisance.eta_ia < -3.0 || nuisance.eta_ia> 3.0) return 0;
  if (nuisance.eta_ia_highz < -1.0 || nuisance.eta_ia_highz> 1.0) return 0;
  // if(like.IA!=0){
  //  if (check_LF()) return 0;
  // }
return 1;
}

int set_nuisance_cluster_Mobs(double cluster_Mobs_lgN0,  double cluster_Mobs_alpha, double cluster_Mobs_beta, double cluster_Mobs_sigma0, double cluster_Mobs_sigma_qm, double cluster_Mobs_sigma_qz)
{
  //  nuisance.cluster_Mobs_lgM0 = mass_obs_norm;  //fiducial : 1.72+log(1.e+14*0.7); could use e.g. sigma = 0.2 Gaussian prior
  //  nuisance.cluster_Mobs_alpha = mass_obs_slope; //fiducial: 1.08; e.g. sigma = 0.1 Gaussian prior
  //  nuisance.cluster_Mobs_beta = mass_z_slope; //fiducial: 0.0; e.g. sigma = 0.1 Gaussian prior
  //  nuisance.cluster_Mobs_sigma = mass_obs_scatter; //fiducial 0.25; e.g. sigma = 0.05 Gaussian prior

  // fiducial values and priors from Murata et al. (2018) except for redshift-related parameters
  nuisance.cluster_Mobs_lgN0 = cluster_Mobs_lgN0; //fiducial: 3.207, flat prior [0.5, 5.0]
  nuisance.cluster_Mobs_alpha = cluster_Mobs_alpha; //fiducial: 0.993, flat prior [0.0, 2.0]
  nuisance.cluster_Mobs_beta = cluster_Mobs_beta; //fiducial: 0.0, flat prior [-1.5, 1.5]
  nuisance.cluster_Mobs_sigma0 = cluster_Mobs_sigma0; //fiducial: 0.456, flat prior [0.0, 1.5]
  nuisance.cluster_Mobs_sigma_qm = cluster_Mobs_sigma_qm; //fiducial: 0.0, flat prior [-1.5, 1.5]
  nuisance.cluster_Mobs_sigma_qz = cluster_Mobs_sigma_qz; //fiducial: 0.0, flat prior [-1.5, 1.5]

  if (nuisance.cluster_Mobs_lgN0 < 0.5 || nuisance.cluster_Mobs_lgN0 > 5.0) return 0;
  if (nuisance.cluster_Mobs_alpha < 0.00001 || nuisance.cluster_Mobs_alpha > 2.0) return 0;
  if (nuisance.cluster_Mobs_beta < -1.5 || nuisance.cluster_Mobs_beta > 1.5) return 0;
  if (nuisance.cluster_Mobs_sigma0 < 0.00001|| nuisance.cluster_Mobs_sigma0 > 1.5) return 0;
  if (nuisance.cluster_Mobs_sigma_qm < -1.5 && nuisance.cluster_Mobs_sigma_qm > 1.5) return 0;
  if (nuisance.cluster_Mobs_sigma_qz < -1.5 && nuisance.cluster_Mobs_sigma_qz > 1.5)return 0;

return 1;
}

int set_nuisance_gbias(double B1, double B2, double B3, double B4,double B5, double B6, double B7, double B8,double B9, double B10)
{
  int i;
  gbias.b[0] = B1;
  gbias.b[1] = B2;
  gbias.b[2] = B3;
  gbias.b[3] = B4;
  gbias.b[4] = B5;
  gbias.b[5] = B6;
  gbias.b[6] = B7;
  gbias.b[7] = B8;
  gbias.b[8] = B9;
  gbias.b[9] = B10;
  if(like.bias==1){
    for (i = 0; i < 10; i++){
      if (gbias.b[i] < 0.8 || gbias.b[i] > 3.0) return 0;
    }
  }
  return 1;
}

double log_L_clusterMobs_WFIRST()
{
  double log_L = 0.;
  int i;

  log_L -=  pow((nuisance.cluster_Mobs_lgN0 - prior.cluster_Mobs_lgM0[0])/ prior.cluster_Mobs_lgM0[1],2.0);
  log_L -=  pow((nuisance.cluster_Mobs_alpha - prior.cluster_Mobs_alpha[0])/ prior.cluster_Mobs_alpha[1],2.0);
  log_L -=  pow((nuisance.cluster_Mobs_beta - prior.cluster_Mobs_beta[0])/ prior.cluster_Mobs_beta[1],2.0);
  log_L -=  pow((nuisance.cluster_Mobs_sigma0 - prior.cluster_Mobs_sigma0[0])/ prior.cluster_Mobs_sigma0[1],2.0);
  log_L -=  pow((nuisance.cluster_Mobs_sigma_qm - prior.cluster_Mobs_sigma_qm[0])/ prior.cluster_Mobs_sigma_qm[1],2.0);
  log_L -=  pow((nuisance.cluster_Mobs_sigma_qz - prior.cluster_Mobs_sigma_qz[0])/ prior.cluster_Mobs_sigma_qz[1],2.0);
  log_L = 0.5*log_L;

  return log_L;
}

double log_L_3x2pt_clusterN_clusterWL_GRS_SN()
{
  double log_L = 0.;
  int n_param = 7;
  double param_fid[n_param], param_diff[n_param];
  int c, r;
  double table[7][7]={{5.49220273e+05,2.33044116e+05,1.47504189e+05,-5.92824792e+03,-2.08565457e+03,-4.86800735e+05,4.48613425e+04},{2.33044116e+05,2.30971497e+05,3.54103600e+04,-6.40553809e+03,-2.59356689e+03,-8.01835609e+04,3.04467069e+04},{1.47504189e+05,3.54103600e+04,6.49198594e+04,-4.10636480e+02,2.04259100e+01,-1.74504277e+05,1.95029955e+04},{-5.92824792e+03,-6.40553809e+03,-4.10636480e+02,3.86909977e+03,6.54697110e+02,7.12569941e+01,-9.99312405e+02},{-2.08565457e+03,-2.59356689e+03,2.04259100e+01,6.54697110e+02,1.57299983e+02,-1.33066811e+02,-2.86375759e+02},{-4.86800735e+05,-8.01835609e+04,-1.74504277e+05,7.12569941e+01,-1.33066811e+02,9.98008316e+05,-1.24566696e+05},{4.48613425e+04,3.04467069e+04,1.95029955e+04,-9.99312405e+02,-2.86375759e+02 -1.24566696e+05,5.51799238e+04}};

  param_diff[0] = cosmology.Omega_m-prior.Omega_m;
  param_diff[1] = cosmology.sigma_8-prior.sigma_8;
  param_diff[2] = cosmology.n_spec-prior.n_spec;
  param_diff[3] = cosmology.w0-prior.w0;
  param_diff[4] = cosmology.wa-prior.wa;
  param_diff[5] = cosmology.omb-prior.omb;
  param_diff[6] = cosmology.h0-prior.h0;

  log_L = -0.5*do_matrix_mult_invcov(n_param,table, param_diff);

  return log_L;
}

double log_L_3x2pt_clusterN_clusterWL_GRS()
{
  double log_L = 0.;
  int n_param = 7;
  double param_fid[n_param], param_diff[n_param];
  int c, r;
  double table[7][7]={{5.43798917e+05,2.33044116e+05,1.47504189e+05,-9.33646911e+03,-2.56808086e+03 -4.86800735e+05,4.48613258e+04},{2.33044116e+05,2.30971497e+05,3.54103600e+04,-6.40553809e+03,-2.59356689e+03,-8.01835609e+04,3.04467069e+04},{1.47504189e+05,3.54103600e+04,6.49198594e+04 -4.10636480e+02,2.04259100e+01,-1.74504277e+05,1.95029955e+04},{-9.33646911e+03,-6.40553809e+03,-4.10636480e+02,1.46051254e+03,3.62655979e+02,7.12569941e+01,-9.99337099e+02},{-2.56808086e+03,-2.59356689e+03,2.04259100e+01,3.62655979e+02,1.09984619e+02,-1.33066811e+02,-2.86374538e+02},{-4.86800735e+05,-8.01835609e+04,-1.74504277e+05,7.12569941e+01,-1.33066811e+02,9.98008316e+05 -1.24566696e+05},{4.48613258e+04,3.04467069e+04,1.95029955e+04,-9.99337099e+02,-2.86374538e+02,-1.24566696e+05,5.51799219e+04}};

  param_diff[0] = cosmology.Omega_m-prior.Omega_m;
  param_diff[1] = cosmology.sigma_8-prior.sigma_8;
  param_diff[2] = cosmology.n_spec-prior.n_spec;
  param_diff[3] = cosmology.w0-prior.w0;
  param_diff[4] = cosmology.wa-prior.wa;
  param_diff[5] = cosmology.omb-prior.omb;
  param_diff[6] = cosmology.h0-prior.h0;

  log_L = -0.5*do_matrix_mult_invcov(n_param,table, param_diff);

  return log_L;
}

double log_multi_like(double OMM, double S8, double NS, double W0,double WA, double OMB, double H0, double MGSigma, double MGmu, double B1, double B2, double B3, double B4,double B5, double B6, double B7, double B8, double B9, double B10, double SP1, double SP2, double SP3, double SP4, double SP5, double SP6, double SP7, double SP8, double SP9, double SP10, double SPS1, double CP1, double CP2, double CP3, double CP4, double CP5, double CP6, double CP7, double CP8, double CP9, double CP10, double CPS1, double M1, double M2, double M3, double M4, double M5, double M6, double M7, double M8, double M9, double M10, double A_ia, double beta_ia, double eta_ia, double eta_ia_highz, double LF_alpha, double LF_P, double LF_Q, double LF_red_alpha, double LF_red_P, double LF_red_Q,double mass_obs_norm, double mass_obs_slope, double mass_z_slope, double mass_obs_scatter_norm, double mass_obs_scatter_mass_slope, double mass_obs_scatter_z_slope, double GRSB1, double GRSB2, double GRSB3, double GRSB4, double GRSB5, double GRSB6, double GRSB7, double SIGMAP1, double SIGMAP2, double SIGMAP3, double SIGMAP4, double SIGMAP5, double SIGMAP6, double SIGMAP7,double SIGMAZ, double PSHOT, double KSTAR)
{
  int i,j,k,m=0,l;
  static double *pred;
  static double *ell;
  static double *ell_Cluster;
  static double darg;
  double chisqr,a,log_L_prior=0.0;

  if(ell==0){
    pred= create_double_vector(0, like.Ndata-1);
    ell= create_double_vector(0, like.Ncl-1);
    darg=(log(like.lmax)-log(like.lmin))/like.Ncl;
    for (l=0;l<like.Ncl;l++){
      ell[l]=exp(log(like.lmin)+(l+0.5)*darg);
    }
    ell_Cluster= create_double_vector(0, Cluster.lbin-1);
    darg=(log(Cluster.l_max)-log(Cluster.l_min))/Cluster.lbin;
    for (l=0;l<Cluster.lbin;l++){
      ell_Cluster[l]=exp(log(Cluster.l_min)+(l+0.5)*darg);
    }
  }
  if (set_cosmology_params(OMM,S8,NS,W0,WA,OMB,H0,MGSigma,MGmu)==0){
    printf("Cosmology out of bounds\n");
    return -1.0e12;
  }
  set_nuisance_shear_calib(M1,M2,M3,M4,M5,M6,M7,M8,M9,M10);
  if (set_nuisance_shear_photoz(SP1,SP2,SP3,SP4,SP5,SP6,SP7,SP8,SP9,SP10,SPS1)==0){
    printf("Shear photo-z sigma too small\n");
    return -1.0e12;
  }
  if (set_nuisance_clustering_photoz(CP1,CP2,CP3,CP4,CP5,CP6,CP7,CP8,CP9,CP10,CPS1)==0){
    printf("Clustering photo-z sigma too small\n");
    return -1.0e12;
  }
  if (set_nuisance_ia(A_ia,beta_ia,eta_ia,eta_ia_highz,LF_alpha,LF_P,LF_Q,LF_red_alpha,LF_red_P,LF_red_Q)==0){
    printf("IA parameters out of bounds\n");
    return -1.0e12;
  }
  if (set_nuisance_gbias(B1,B2,B3,B4,B5,B6,B7,B8,B9,B10)==0){
    printf("Bias out of bounds\n");
    return -1.0e12;
  }
  if (set_nuisance_cluster_Mobs(mass_obs_norm, mass_obs_slope, mass_z_slope, mass_obs_scatter_norm, mass_obs_scatter_mass_slope, mass_obs_scatter_z_slope)==0){
    printf("Mobs out of bounds\n");
    return -1.0e12;
  }

  //printf("like %le %le %le %le %le %le %le %le\n",cosmology.Omega_m, cosmology.Omega_v,cosmology.sigma_8,cosmology.n_spec,cosmology.w0,cosmology.wa,cosmology.omb,cosmology.h0);
  // printf("like %le %le %le %le\n",gbias.b[0][0], gbias.b[1][0], gbias.b[2][0], gbias.b[3][0]);
  // for (i=0; i<10; i++){
  //   printf("nuisance %le %le %le\n",nuisance.shear_calibration_m[i],nuisance.bias_zphot_shear[i],nuisance.sigma_zphot_shear[i]);
  // }

  log_L_prior=log_L_3x2pt_clusterN_clusterWL_GRS();
  //log_L_prior+=log_L_3x2pt_clusterN_clusterWL_GRS_SN();
  // if(like.wlphotoz!=0) log_L_prior+=log_L_wlphotoz();
  // if(like.clphotoz!=0) log_L_prior+=log_L_clphotoz();
  // if(like.shearcalib==1) log_L_prior+=log_L_shear_calib();
  // if(strcmp(like.ext_data,"WFIRST_SN")==0) log_L_prior+=log_L_SN_WFIRST_w0wa();
  // if(like.GRS==1) log_L_prior+=log_like_GRS(OMM, S8, NS, W0,WA, OMB, H0, MGSigma, MGmu, GRSB1, GRSB2, GRSB3, GRSB4, GRSB5, GRSB6, GRSB7, SIGMAP1, SIGMAP2, SIGMAP3, SIGMAP4, SIGMAP5, SIGMAP6, SIGMAP7,SIGMAZ, PSHOT, KSTAR);
  // if(like.clusterMobs==1) log_L_prior+=log_L_clusterMobs_WFIRST();

  //printf("%d %d %d %d\n",like.BAO,like.wlphotoz,like.clphotoz,like.shearcalib);
  // printf("logl %le %le %le %le\n",log_L_shear_calib(),log_L_wlphotoz(),log_L_clphotoz(),log_L_clusterMobs());
  int start=0;

  // if(like.shear_shear==1) {
  //   set_data_shear(like.Ncl, ell, pred, start);
  //   start=start+like.Ncl*tomo.shear_Npowerspectra;
  // }
  // if(like.shear_pos==1){
  //   set_data_ggl(like.Ncl, ell, pred, start);
  //   start=start+like.Ncl*tomo.ggl_Npowerspectra;
  // }
  // if(like.pos_pos==1){
  //   set_data_clustering(like.Ncl,ell,pred, start);
  //   start=start+like.Ncl*tomo.clustering_Npowerspectra;
  // }
  // if(like.clusterN==1){
  //   set_data_cluster_N(pred,start);
  //   start= start+tomo.cluster_Nbin*Cluster.N200_Nbin;
  // }
  // if(like.clusterWL==1){
  //   set_data_cgl(ell_Cluster,pred, start);
  // }
  chisqr=0.0;
  // for (i=0; i<like.Ndata; i++){
  //   for (j=0; j<like.Ndata; j++){
  //     a=(pred[i]-data_read(1,i))*invcov_read(1,i,j)*(pred[j]-data_read(1,j));
  //     chisqr=chisqr+a;
  //   }
    // if (fabs(data_read(1,i)) < 1.e-25){
    //    printf("%d %le %le %le\n",i,data_read(1,i),pred[i],invcov_read(1,i,i));
    // }
  // }
  // if (chisqr<0.0){
  //   printf("error: chisqr = %le\n",chisqr);
  //   //exit(EXIT_FAILURE);
  // }
  // if (like.GRS == 1){
  //   printf("like_grs activated!\n");
  //   log_L_GRS = log_like_GRS(cosmology.Omega_m, cosmology.sigma_8, cosmology.n_spec, cosmology.w0,cosmology.wa, cosmology.omb, cosmology.h0,-42.,-42.,-42.,-42.,-42.,-42.,-42.,-42.,-42.,-42.,-42.,-42.,-42.,-42.,-42.,-42.,-42.);
  // }
 //printf("%le\n",chisqr);
//  return -0.5*chisqr+log_L_prior+log_L_GRS;
  return -0.5*chisqr+log_L_prior;
}

void compute_data_vector(char *details, double OMM, double S8, double NS, double W0,double WA, double OMB, double H0, double MGSigma, double MGmu, double B1, double B2, double B3, double B4,double B5, double B6, double B7, double B8, double B9, double B10, double SP1, double SP2, double SP3, double SP4, double SP5,double SP6, double SP7, double SP8, double SP9, double SP10, double SPS1, double CP1, double CP2, double CP3, double CP4, double CP5, double CP6, double CP7, double CP8, double CP9, double CP10, double CPS1, double M1, double M2, double M3, double M4, double M5, double M6, double M7, double M8, double M9, double M10, double A_ia, double beta_ia, double eta_ia, double eta_ia_highz, double LF_alpha, double LF_P, double LF_Q, double LF_red_alpha, double LF_red_P, double LF_red_Q, double mass_obs_norm, double mass_obs_slope, double mass_z_slope, double mass_obs_scatter_norm, double mass_obs_scatter_mass_slope, double mass_obs_scatter_z_slope)
{

  int i,j,k,m=0,l;
  static double *pred;
  static double *ell;
  static double *ell_Cluster;
  static double darg;
  double chisqr,a,log_L_prior=0.0;

  if(ell==0){
    pred= create_double_vector(0, like.Ndata-1);
    ell= create_double_vector(0, like.Ncl-1);
    darg=(log(like.lmax)-log(like.lmin))/like.Ncl;
    for (l=0;l<like.Ncl;l++){
      ell[l]=exp(log(like.lmin)+(l+0.5)*darg);
    }
    ell_Cluster= create_double_vector(0, Cluster.lbin-1);
    darg=(log(Cluster.l_max)-log(Cluster.l_min))/Cluster.lbin;
    for (l=0;l<Cluster.lbin;l++){
      ell_Cluster[l]=exp(log(Cluster.l_min)+(l+0.5)*darg);
    }
  }
// for (l=0;l<like.Ncl;l++){
//   printf("%d %le\n",i,ell[l]);
// }
  set_cosmology_params(OMM,S8,NS,W0,WA,OMB,H0,MGSigma,MGmu);
  set_nuisance_shear_calib(M1,M2,M3,M4,M5,M6,M7,M8,M9,M10);
  set_nuisance_shear_photoz(SP1,SP2,SP3,SP4,SP5,SP6,SP7,SP8,SP9,SP10,SPS1);
  set_nuisance_clustering_photoz(CP1,CP2,CP3,CP4,CP5,CP6,CP7,CP8,CP9,CP10,CPS1);
  set_nuisance_ia(A_ia,beta_ia,eta_ia,eta_ia_highz,LF_alpha,LF_P,LF_Q,LF_red_alpha,LF_red_P,LF_red_Q);
  set_nuisance_gbias(B1,B2,B3,B4,B5,B6,B7,B8,B9,B10);
  set_nuisance_cluster_Mobs(mass_obs_norm,mass_obs_slope,mass_z_slope,mass_obs_scatter_norm,mass_obs_scatter_mass_slope,mass_obs_scatter_z_slope);

  int start=0;
  if(like.shear_shear==1) {
    set_data_shear(like.Ncl, ell, pred, start);
    start=start+like.Ncl*tomo.shear_Npowerspectra;
  }
  if(like.shear_pos==1){
    //printf("ggl\n");
    set_data_ggl(like.Ncl, ell, pred, start);
    start=start+like.Ncl*tomo.ggl_Npowerspectra;
  }
  if(like.pos_pos==1){
    //printf("clustering\n");
    set_data_clustering(like.Ncl,ell,pred, start);
    start=start+like.Ncl*tomo.clustering_Npowerspectra;
  }

  if(like.clusterN==1){
    set_data_cluster_N(pred,start);
    start= start+tomo.cluster_Nbin*Cluster.N200_Nbin;
  }
  if(like.clusterWL==1){
    set_data_cgl(ell_Cluster,pred, start);
  }
  FILE *F;
  char filename[300];
  if (strstr(details,"FM") != NULL){
    sprintf(filename,"%s",details);
  }
  else {sprintf(filename,"datav/%s_%s_%s",survey.name,like.probes,details);}
  F=fopen(filename,"w");
  for (i=0;i<like.Ndata; i++){
    fprintf(F,"%d %le\n",i,pred[i]);
    //printf("%d %le\n",i,pred[i]);
  }
  fclose(F);
}

void write_vector_wrapper(char *details, input_cosmo_params ic, input_nuisance_params in)
{
  compute_data_vector(details, ic.omega_m, ic.sigma_8, ic.n_s, ic.w0, ic.wa, ic.omega_b, ic.h0, ic.MGSigma, ic.MGmu,
    in.bias[0], in.bias[1], in.bias[2], in.bias[3],in.bias[4], in.bias[5], in.bias[6], in.bias[7],in.bias[8], in.bias[9],
    in.source_z_bias[0], in.source_z_bias[1], in.source_z_bias[2], in.source_z_bias[3], in.source_z_bias[4],
    in.source_z_bias[5], in.source_z_bias[6], in.source_z_bias[7], in.source_z_bias[8], in.source_z_bias[9],
    in.source_z_s,
    in.lens_z_bias[0], in.lens_z_bias[1], in.lens_z_bias[2], in.lens_z_bias[3], in.lens_z_bias[4],
    in.lens_z_bias[5], in.lens_z_bias[6], in.lens_z_bias[7], in.lens_z_bias[8], in.lens_z_bias[9],
    in.lens_z_s,
    in.shear_m[0], in.shear_m[1], in.shear_m[2], in.shear_m[3], in.shear_m[4],
    in.shear_m[5], in.shear_m[6], in.shear_m[7], in.shear_m[8], in.shear_m[9],
    in.A_ia, in.beta_ia, in.eta_ia, in.eta_ia_highz,
    in.lf[0], in.lf[1], in.lf[2], in.lf[3], in.lf[4], in.lf[5],
    in.m_lambda[0], in.m_lambda[1], in.m_lambda[2], in.m_lambda[3],
    in.m_lambda[4], in.m_lambda[5]);
}

double log_like_wrapper(input_cosmo_params ic, input_nuisance_params in,input_nuisance_params_grs ingr)
{
  double like = log_multi_like(ic.omega_m, ic.sigma_8, ic.n_s, ic.w0, ic.wa, ic.omega_b, ic.h0, ic.MGSigma, ic.MGmu,
    in.bias[0], in.bias[1], in.bias[2], in.bias[3],in.bias[4], in.bias[5], in.bias[6], in.bias[7],in.bias[8], in.bias[9],
    in.source_z_bias[0], in.source_z_bias[1], in.source_z_bias[2], in.source_z_bias[3], in.source_z_bias[4],
    in.source_z_bias[5], in.source_z_bias[6], in.source_z_bias[7], in.source_z_bias[8], in.source_z_bias[9],
    in.source_z_s,
    in.lens_z_bias[0], in.lens_z_bias[1], in.lens_z_bias[2], in.lens_z_bias[3], in.lens_z_bias[4],
    in.lens_z_bias[5], in.lens_z_bias[6], in.lens_z_bias[7], in.lens_z_bias[8], in.lens_z_bias[9],
    in.lens_z_s,
    in.shear_m[0], in.shear_m[1], in.shear_m[2], in.shear_m[3], in.shear_m[4],
    in.shear_m[5], in.shear_m[6], in.shear_m[7], in.shear_m[8], in.shear_m[9],
    in.A_ia, in.beta_ia, in.eta_ia, in.eta_ia_highz,
    in.lf[0], in.lf[1], in.lf[2], in.lf[3], in.lf[4], in.lf[5],
    in.m_lambda[0], in.m_lambda[1], in.m_lambda[2], in.m_lambda[3],
    in.m_lambda[4], in.m_lambda[5],ingr.grsbias[0],ingr.grsbias[1],ingr.grsbias[2],ingr.grsbias[3],ingr.grsbias[4],ingr.grsbias[5],ingr.grsbias[6],ingr.grssigmap[0],ingr.grssigmap[1],ingr.grssigmap[2],ingr.grssigmap[3],ingr.grssigmap[4],ingr.grssigmap[5],ingr.grssigmap[6],ingr.grssigmaz,ingr.grspshot,ingr.grskstar);
  return like;
}

void save_zdistr_sources(int zs){
  double z,dz =(redshift.shear_zdistrpar_zmax-redshift.shear_zdistrpar_zmin)/300.0;
  printf("Printing redshift distribution n(z) for source redshift bin %d\n",zs);

   FILE *F1;
   char filename[300];
   sprintf(filename,"zdistris/pessi_zdist_sources_bin%d.txt",zs);
   F1 = fopen(filename,"w");
   for (z =redshift.shear_zdistrpar_zmin; z< redshift.shear_zdistrpar_zmax; z+= dz){
      fprintf(F1,"%e %e\n", z, zdistr_photoz(z,zs));
   }
}

void save_zdistr_lenses(int zl){
   double z,dz =(redshift.clustering_zdistrpar_zmax-redshift.clustering_zdistrpar_zmin)/300.0;
  printf("Printing redshift distribution n(z) and bias b(z) for lens redshift bin %d\n",zl);

   FILE *F1;
   char filename[300];
   sprintf(filename,"zdistris/pessi_zdist_lenses_bin%d.txt", zl);
   F1 = fopen(filename,"w");
   for (z =redshift.clustering_zdistrpar_zmin; z< redshift.clustering_zdistrpar_zmax; z+= dz){
      fprintf(F1,"%e %e\n", z, pf_photoz(z,zl));
   }
}

int main(int argc, char** argv)
{
  clock_t begin, end;
  double time_spent,loglike=0.0;
  int i;
/* here, do your time-consuming job */

  init_cosmo_runmode("halofit");
  init_binning_fourier(25,30.0,15000.0,4000.0,21.0,10,10);
  if(strcmp(argv[1],"opti")==0) init_priors("photo_opti","shear_opti","none","none");
  if(strcmp(argv[1],"pessi")==0) init_priors("photo_pessi","shear_pessi","none","none");
  init_survey(argv[2]);
  if(strcmp(argv[2],"WFIRST")==0) init_galaxies("zdistris/zdistri_WFIRST_LSST_lensing_fine_bin","zdistris/zdistri_WFIRST_LSST_clustering_fine_bin", "gaussian", "gaussian", "SN10");
  init_clusters();
  init_IA("none", "none");
  init_probes(argv[3]);

  if(strcmp(argv[1],"opti")==0) compute_data_vector(argv[1],0.3156,0.831,0.9645,-1.,0.,0.0491685,0.6727,0.,0.,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.1,2.2,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.01,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.01,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,5.92,1.1,-0.47,0.0,0.0,0.0,0.0,0.0,0.0,0.0,3.207,0.993,0.0,0.456,0.0,0.0);
  if(strcmp(argv[1],"pessi")==0) compute_data_vector(argv[1],0.3156,0.831,0.9645,-1.,0.,0.0491685,0.6727,0.,0.,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.1,2.2,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.05,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.05,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,5.92,1.1,-0.47,0.0,0.0,0.0,0.0,0.0,0.0,0.0,3.207,0.993,0.0,0.456,0.0,0.0);

  // init_data_inv("cov/WFIRST_3x2pt_inv","datav/WFIRST_all_2pt_fid_opti");

  // begin = clock();
  // loglike=log_multi_like(0.3156,0.831,0.9645,-1.,0.,0.0491685,0.6727,0.,0.,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.1,2.2,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.05,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,-0.0005,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,5.92,1.1,-0.47,0.0,0.0,0.0,0.0,0.0,0.0,0.0,3.207,0.993,0.0,0.456,0.0,0.0);
  // printf("%le\n",loglike);
  // // printf("knonlin %le\n",nonlinear_scale_computation(1.0));
  // // printf("knonlin %le\n",nonlinear_scale_computation(0.5));
  // end = clock();
  // time_spent = (double)(end - begin) / CLOCKS_PER_SEC;
  // printf("timespent %le\n",time_spent);

  //CH BEGINS
  //for testing Planck15_BAO_w0wa prior alone
  //compute_data_vector("fid",3.50989e-01,8.04675e-01,9.64061e-01,-5.05518e-01,-1.46884e+00,5.46245e-02,6.39839e-01,0.,0.,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.1,2.2,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.05,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.01,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,5.92,1.1,-0.47,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.72+log(1.e+14*0.7),1.08,0.0,0.25,0.9,0.9,0.9,0.9);
  //expect 0.0 for the following
  //log_multi_like(3.50989e-01,8.04675e-01.1,9.64061e-01,-5.05518e-01,-1.46884e+00,5.46245e-02,6.39839e-01,0.,0.,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.1,2.2,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.05,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.01,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,5.92,1.1,-0.47,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.72+log(1.e+14*0.7),1.08,0.0,0.25,0.9,0.9,0.9,0.9);
  //expect -13.195605 for the following
  //log_multi_like(0.35449914, 0.81272201,0.97370111, -0.51057289,-1.48353327,0.05517077,0.64623714,0.,0.,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.1,2.2,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.05,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.01,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,5.92,1.1,-0.47,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.72+log(1.e+14*0.7),1.08,0.0,0.25,0.9,0.9,0.9,0.9);
  //CH ENDS

  return 0;
}


text:
/*
    This code is written by <albanese@fbk.it>.
    (C) 2010 mlpy Developers.

    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, see <http://www.gnu.org/licenses/>.
*/

#include <stdlib.h>
#include <math.h>
#include <float.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

#define MIN( A , B ) ((A) < (B) ? (A) : (B))

#define INIT_STD 0
#define INIT_PLUSPLUS 1

void init_std(double *data,      /* data points (nn points x pp dimensions) */
	      double *means,     /* means (kk clusters x pp dimensions) */
	      int nn,            /* number od data points */
	      int pp,            /* number of dimensions */
	      int kk,            /* number of clusters */
	      unsigned long seed /* random seed for init */
	      )
{
  int n, p, k;
  int *ridx;
  const gsl_rng_type * T;
  gsl_rng * r;

  T = gsl_rng_default;
  r = gsl_rng_alloc (T);
  gsl_rng_set (r, seed);

  ridx = (int *) malloc (nn * sizeof(int));

  for (n=0; n<nn; n++)
    ridx[n] = n;

  gsl_ran_shuffle (r, ridx, nn, sizeof (int));

  for (k=0; k<kk; k++)
    for (p=0; p<pp; p++)
      means[p + (k * pp)] = data[p + (ridx[k] * pp)];

  free(ridx);
}

/* for init_plus */
void
dist_min(double *a, double *b, int nn)
{
  int n;

  for (n=0; n<nn; n++)
    a[n] = MIN (a[n], b[n]);
}

/* for init_plus */
int
idx_max(double *a, int nn)
{
  int n, idx = 0;
  double max = -DBL_MAX;

  for (n=0; n<nn; n++)
    if (a[n] > max)
      {
	max = a[n];
	idx = n;
      }

  return idx;
}

void
init_plus(double *data,      /* data points (nn points x pp dimensions) */
	  double *means,     /* means (kk clusters x pp dimensions) */
	  int nn,            /* number od data points */
	  int pp,            /* number of dimensions */
	  int kk,            /* number of clusters */
	  unsigned long seed /* random seed for init */
	  )
{
  int n, p, k;
  double *dist, *distk;
  int sidx;

  const gsl_rng_type *T;
  gsl_rng *r;

  T = gsl_rng_default;
  r = gsl_rng_alloc (T);
  gsl_rng_set (r, seed);

  dist = (double *) malloc (nn * sizeof(double));
  distk = (double *) malloc (nn * sizeof(double));

  /* first mean (randomly selected) */
  sidx = (int) gsl_rng_uniform_int (r, nn);
  gsl_rng_free(r);
  for (p=0; p<pp; p++)
    means[p] = data[p + (sidx * pp)];

  /* initialize dist */
  for (n=0; n<nn; n++) /* for each data point */
    dist[n] = DBL_MAX;

  for (k=0; k<kk-1; k++)
    {
      /* for each data point x compute distance from mean k */
      for (n=0; n<nn; n++)
	{
	  distk[n] = 0.0;
	  for (p=0; p<pp; p++)
	    distk[n] += pow(data[p + (n * pp)] - means[p + (k * pp)], 2);
	}

      /* for each data point x, compute the distance between x and */
      /* the nearest center that has already been chosen */
      dist_min (dist, distk, nn);

      /* add one new data point as a new center, using */
      sidx = idx_max (dist, nn);

      for (p=0; p<pp; p++)
	means[p + ((k+1) * pp)] = data[p + (sidx * pp)];
    }

  free(dist);
  free(distk);
}

/* assignment step */
int
a_step(double *data,  /* data points (nn x pp) */
       double *means, /* means (kk x pp) */
       int *cls,      /* cluster assignement for each data point (nn) */
       int *nelems,   /* number of elements of each cluster (kk) */
       int nn,        /* number od data points */
       int pp,        /* number of dimensions */
       int kk         /* number of clusters */
       )
{
  int n, p, k, kn = 0;
  double dist, dmin;
  int changed = 0;

  for (k=0; k<kk; k++)
    nelems[k] = 0;

  for (n=0; n<nn; n++) /* for each data point */
    {
      dmin = DBL_MAX;
      for (k=0; k<kk; k++) /* for each cluster */
	{
	  /* compute distance */
	  dist = 0.0;
	  for (p=0; p<pp; p++)
	    dist += pow(data[p + (n * pp)] - means[p + (k * pp)], 2);

	  /* remember the cluster k if dist < dmin */
	  if (dist < dmin)
	    {
	      dmin = dist;
	      kn = k;
	    }
	}

      /* if the cluster assignement change */
      if (kn != cls[n])
	changed++;

      /* update clusters and number of elements of each cluster */
      cls[n] = kn;
      nelems[kn]++;
    }

  return changed;
}

/* update step */
int
u_step(double *data,  /* data points (nn x pp) */
       double *means, /* means (kk x pp) */
       int *cls,      /* cluster assignement for each data point (nn) */
       int *nelems,   /* number of elements of each cluster (kk) */
       int nn,        /* number od data points */
       int pp,        /* number of dimensions */
       int kk         /* number of clusters */
       )
{
  int n, p, k;

  /* reset means */
  for (k=0; k<kk; k++)
    for (p=0; p<pp; p++)
      means[p + (k * pp)] = 0.0;

  for (n=0; n<nn; n++) /* for each data point */
    for (p=0; p<pp; p++) /* for each dimension */
      means[p + (cls[n] * pp)] += data[p + (n * pp)];

  for (k=0; k<kk; k++)
    if (nelems[k] > 0)
      for (p=0; p<pp; p++)
	means[p + (k * pp)] /= nelems[k];

  return 1;
}

int
km(double *data,  /* data points (nn x pp) */
   double *means, /* initialized means (kk x pp) */
   int *cls,      /* cluster assignement for each data point (nn) */
   int nn,        /* number od data points */
   int pp,        /* number of dimensions */
   int kk         /* number of clusters */
   )
{
  int n;
  int ret, steps = 0;
  int changed = -1;
  int *nelems; /* number of elements of each cluster (kk) */

  nelems = (int *) malloc (kk * sizeof(int));

  /* init cls */
  for (n=0; n<nn; n++)
    cls[n] = -1.0;

  /* k-means algorithm  */
  while (changed != 0)
    {
      changed = a_step (data, means, cls, nelems, nn, pp, kk);
      ret = u_step (data, means, cls, nelems, nn, pp, kk);
      steps++;
    }

  free(nelems);

  return steps;
}


text:
/**
 *
 * @file qwrapper_slaset.c
 *
 *  PLASMA core_blas quark wrapper
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Hatem Ltaief
 * @date 2010-11-15
 * @generated s Tue Jan  7 11:44:58 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 **/
void QUARK_CORE_slaset(Quark *quark, Quark_Task_Flags *task_flags,
                       PLASMA_enum uplo, int M, int N,
                       float alpha, float beta,
                       float *A, int LDA)
{
    DAG_CORE_LASET;
    QUARK_Insert_Task(quark, CORE_slaset_quark, task_flags,
        sizeof(PLASMA_enum),                &uplo,  VALUE,
        sizeof(int),                        &M,     VALUE,
        sizeof(int),                        &N,     VALUE,
        sizeof(float),         &alpha, VALUE,
        sizeof(float),         &beta,  VALUE,
        sizeof(float)*LDA*N,    A,      OUTPUT,
        sizeof(int),                        &LDA,   VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_slaset_quark = PCORE_slaset_quark
#define CORE_slaset_quark PCORE_slaset_quark
#endif
void CORE_slaset_quark(Quark *quark)
{
    PLASMA_enum uplo;
    int M;
    int N;
    float alpha;
    float beta;
    float *A;
    int LDA;

    quark_unpack_args_7(quark, uplo, M, N, alpha, beta, A, LDA);
    LAPACKE_slaset_work(
        LAPACK_COL_MAJOR,
        lapack_const(uplo),
        M, N, alpha, beta, A, LDA);
}


text:
/* specfunc/mathieu_radfunc.c
 *
 * Copyright (C) 2002 Lowell Johnson
 *
 * 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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Author:  L. Johnson */

#include <config.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_sf_mathieu.h>

int gsl_sf_mathieu_Mc_e(int kind, int order, double qq, double zz,
                      gsl_sf_result *result)
{
  int even_odd, kk, mm, status;
  double maxerr = 1e-14, amax, pi = M_PI, fn, factor;
  double coeff[GSL_SF_MATHIEU_COEFF], fc;
  double j1c, z2c, j1pc, z2pc;
  double u1, u2;
  gsl_sf_result aa;

  /* Check for out of bounds parameters. */
  if (qq <= 0.0)
  {
      GSL_ERROR("q must be greater than zero", GSL_EINVAL);
  }
  if (kind < 1 || kind > 2)
  {
      GSL_ERROR("kind must be 1 or 2", GSL_EINVAL);
  }

  mm = 0;
  amax = 0.0;
  fn = 0.0;
  u1 = sqrt(qq)*exp(-1.0*zz);
  u2 = sqrt(qq)*exp(zz);

  even_odd = 0;
  if (order % 2 != 0)
      even_odd = 1;

  /* Compute the characteristic value. */
  status = gsl_sf_mathieu_a_e(order, qq, &aa);
  if (status != GSL_SUCCESS)
  {
      return status;
  }

  /* Compute the series coefficients. */
  status = gsl_sf_mathieu_a_coeff(order, qq, aa.val, coeff);
  if (status != GSL_SUCCESS)
  {
      return status;
  }

  if (even_odd == 0)
  {
      for (kk=0; kk<GSL_SF_MATHIEU_COEFF; kk++)
      {
          amax = GSL_MAX(amax, fabs(coeff[kk]));
          if (fabs(coeff[kk])/amax < maxerr)
              break;

          j1c = gsl_sf_bessel_Jn(kk, u1);
          if (kind == 1)
          {
              z2c = gsl_sf_bessel_Jn(kk, u2);
          }
          else /* kind = 2 */
          {
              z2c = gsl_sf_bessel_Yn(kk, u2);
          }

          fc = pow(-1.0, 0.5*order+kk)*coeff[kk];
          fn += fc*j1c*z2c;
      }

      fn *= sqrt(pi/2.0)/coeff[0];
  }
  else
  {
      for (kk=0; kk<GSL_SF_MATHIEU_COEFF; kk++)
      {
          amax = GSL_MAX(amax, fabs(coeff[kk]));
          if (fabs(coeff[kk])/amax < maxerr)
              break;

          j1c = gsl_sf_bessel_Jn(kk, u1);
          j1pc = gsl_sf_bessel_Jn(kk+1, u1);
          if (kind == 1)
          {
              z2c = gsl_sf_bessel_Jn(kk, u2);
              z2pc = gsl_sf_bessel_Jn(kk+1, u2);
          }
          else /* kind = 2 */
          {
              z2c = gsl_sf_bessel_Yn(kk, u2);
              z2pc = gsl_sf_bessel_Yn(kk+1, u2);
          }
          fc = pow(-1.0, 0.5*(order-1)+kk)*coeff[kk];
          fn += fc*(j1c*z2pc + j1pc*z2c);
      }

      fn *= sqrt(pi/2.0)/coeff[0];
  }

  result->val = fn;
  result->err = 2.0*GSL_DBL_EPSILON;
  factor = fabs(fn);
  if (factor > 1.0)
      result->err *= factor;

  return GSL_SUCCESS;
}

int gsl_sf_mathieu_Ms_e(int kind, int order, double qq, double zz,
                      gsl_sf_result *result)
{
  int even_odd, kk, mm, status;
  double maxerr = 1e-14, amax, pi = M_PI, fn, factor;
  double coeff[GSL_SF_MATHIEU_COEFF], fc;
  double j1c, z2c, j1mc, z2mc, j1pc, z2pc;
  double u1, u2;
  gsl_sf_result aa;

  /* Check for out of bounds parameters. */
  if (qq <= 0.0)
  {
      GSL_ERROR("q must be greater than zero", GSL_EINVAL);
  }
  if (kind < 1 || kind > 2)
  {
      GSL_ERROR("kind must be 1 or 2", GSL_EINVAL);
  }

  /* Handle the trivial cases where order = 0. */
  if (order == 0)
  {
      result->val = 0.0;
      result->err = 0.0;
      return GSL_SUCCESS;
  }

  mm = 0;
  amax = 0.0;
  fn = 0.0;
  u1 = sqrt(qq)*exp(-1.0*zz);
  u2 = sqrt(qq)*exp(zz);

  even_odd = 0;
  if (order % 2 != 0)
      even_odd = 1;

  /* Compute the characteristic value. */
  status = gsl_sf_mathieu_b_e(order, qq, &aa);
  if (status != GSL_SUCCESS)
  {
      return status;
  }

  /* Compute the series coefficients. */
  status = gsl_sf_mathieu_b_coeff(order, qq, aa.val, coeff);
  if (status != GSL_SUCCESS)
  {
      return status;
  }

  if (even_odd == 0)
  {
      for (kk=0; kk<GSL_SF_MATHIEU_COEFF; kk++)
      {
          amax = GSL_MAX(amax, fabs(coeff[kk]));
          if (fabs(coeff[kk])/amax < maxerr)
              break;

          j1mc = gsl_sf_bessel_Jn(kk, u1);
          j1pc = gsl_sf_bessel_Jn(kk+2, u1);
          if (kind == 1)
          {
              z2mc = gsl_sf_bessel_Jn(kk, u2);
              z2pc = gsl_sf_bessel_Jn(kk+2, u2);
          }
          else /* kind = 2 */
          {
              z2mc = gsl_sf_bessel_Yn(kk, u2);
              z2pc = gsl_sf_bessel_Yn(kk+2, u2);
          }

          fc = pow(-1.0, 0.5*order+kk+1)*coeff[kk];
          fn += fc*(j1mc*z2pc - j1pc*z2mc);
      }

      fn *= sqrt(pi/2.0)/coeff[0];
  }
  else
  {
      for (kk=0; kk<GSL_SF_MATHIEU_COEFF; kk++)
      {
          amax = GSL_MAX(amax, fabs(coeff[kk]));
          if (fabs(coeff[kk])/amax < maxerr)
              break;

          j1c = gsl_sf_bessel_Jn(kk, u1);
          j1pc = gsl_sf_bessel_Jn(kk+1, u1);
          if (kind == 1)
          {
              z2c = gsl_sf_bessel_Jn(kk, u2);
              z2pc = gsl_sf_bessel_Jn(kk+1, u2);
          }
          else /* kind = 2 */
          {
              z2c = gsl_sf_bessel_Yn(kk, u2);
              z2pc = gsl_sf_bessel_Yn(kk+1, u2);
          }

          fc = pow(-1.0, 0.5*(order-1)+kk)*coeff[kk];
          fn += fc*(j1c*z2pc - j1pc*z2c);
      }

      fn *= sqrt(pi/2.0)/coeff[0];
  }

  result->val = fn;
  result->err = 2.0*GSL_DBL_EPSILON;
  factor = fabs(fn);
  if (factor > 1.0)
      result->err *= factor;

  return GSL_SUCCESS;
}

int gsl_sf_mathieu_Mc_array(int kind, int nmin, int nmax, double qq,
                            double zz, gsl_sf_mathieu_workspace *work,
                            double result_array[])
{
  int even_odd, order, ii, kk, mm, status;
  double maxerr = 1e-14, amax, pi = M_PI, fn;
  double coeff[GSL_SF_MATHIEU_COEFF], fc;
  double j1c, z2c, j1pc, z2pc;
  double u1, u2;
  double *aa = work->aa;

  /* Initialize the result array to zeroes. */
  for (ii=0; ii<nmax-nmin+1; ii++)
      result_array[ii] = 0.0;

  /* Check for out of bounds parameters. */
  if (qq <= 0.0)
  {
      GSL_ERROR("q must be greater than zero", GSL_EINVAL);
  }
  if (kind < 1 || kind > 2)
  {
      GSL_ERROR("kind must be 1 or 2", GSL_EINVAL);
  }

  mm = 0;
  amax = 0.0;
  u1 = sqrt(qq)*exp(-1.0*zz);
  u2 = sqrt(qq)*exp(zz);

  /* Compute all eigenvalues up to nmax. */
  gsl_sf_mathieu_a_array(0, nmax, qq, work, aa);

  for (ii=0, order=nmin; order<=nmax; ii++, order++)
  {
      fn = 0.0;
      even_odd = 0;
      if (order % 2 != 0)
          even_odd = 1;

      /* Compute the series coefficients. */
      status = gsl_sf_mathieu_a_coeff(order, qq, aa[order], coeff);
      if (status != GSL_SUCCESS)
      {
          return status;
      }

      if (even_odd == 0)
      {
          for (kk=0; kk<GSL_SF_MATHIEU_COEFF; kk++)
          {
              amax = GSL_MAX(amax, fabs(coeff[kk]));
              if (fabs(coeff[kk])/amax < maxerr)
                  break;

              j1c = gsl_sf_bessel_Jn(kk, u1);
              if (kind == 1)
              {
                  z2c = gsl_sf_bessel_Jn(kk, u2);
              }
              else /* kind = 2 */
              {
                  z2c = gsl_sf_bessel_Yn(kk, u2);
              }

              fc = pow(-1.0, 0.5*order+kk)*coeff[kk];
              fn += fc*j1c*z2c;
          }

          fn *= sqrt(pi/2.0)/coeff[0];
      }
      else
      {
          for (kk=0; kk<GSL_SF_MATHIEU_COEFF; kk++)
          {
              amax = GSL_MAX(amax, fabs(coeff[kk]));
              if (fabs(coeff[kk])/amax < maxerr)
                  break;

              j1c = gsl_sf_bessel_Jn(kk, u1);
              j1pc = gsl_sf_bessel_Jn(kk+1, u1);
              if (kind == 1)
              {
                  z2c = gsl_sf_bessel_Jn(kk, u2);
                  z2pc = gsl_sf_bessel_Jn(kk+1, u2);
              }
              else /* kind = 2 */
              {
                  z2c = gsl_sf_bessel_Yn(kk, u2);
                  z2pc = gsl_sf_bessel_Yn(kk+1, u2);
              }
              fc = pow(-1.0, 0.5*(order-1)+kk)*coeff[kk];
              fn += fc*(j1c*z2pc + j1pc*z2c);
          }

          fn *= sqrt(pi/2.0)/coeff[0];
      }

      result_array[ii] = fn;
  } /* order loop */

  return GSL_SUCCESS;
}

int gsl_sf_mathieu_Ms_array(int kind, int nmin, int nmax, double qq,
                            double zz, gsl_sf_mathieu_workspace *work,
                            double result_array[])
{
  int even_odd, order, ii, kk, mm, status;
  double maxerr = 1e-14, amax, pi = M_PI, fn;
  double coeff[GSL_SF_MATHIEU_COEFF], fc;
  double j1c, z2c, j1mc, z2mc, j1pc, z2pc;
  double u1, u2;
  double *bb = work->bb;

  /* Initialize the result array to zeroes. */
  for (ii=0; ii<nmax-nmin+1; ii++)
      result_array[ii] = 0.0;

  /* Check for out of bounds parameters. */
  if (qq <= 0.0)
  {
      GSL_ERROR("q must be greater than zero", GSL_EINVAL);
  }
  if (kind < 1 || kind > 2)
  {
      GSL_ERROR("kind must be 1 or 2", GSL_EINVAL);
  }

  mm = 0;
  amax = 0.0;
  u1 = sqrt(qq)*exp(-1.0*zz);
  u2 = sqrt(qq)*exp(zz);

  /* Compute all eigenvalues up to nmax. */
  gsl_sf_mathieu_b_array(0, nmax, qq, work, bb);

  for (ii=0, order=nmin; order<=nmax; ii++, order++)
  {
      fn = 0.0;
      even_odd = 0;
      if (order % 2 != 0)
          even_odd = 1;

      /* Handle the trivial cases where order = 0. */
      if (order == 0)
      {
          result_array[ii] = 0.0;
          continue;
      }

      /* Compute the series coefficients. */
      status = gsl_sf_mathieu_b_coeff(order, qq, bb[order], coeff);
      if (status != GSL_SUCCESS)
      {
          return status;
      }

      if (even_odd == 0)
      {
          for (kk=0; kk<GSL_SF_MATHIEU_COEFF; kk++)
          {
              amax = GSL_MAX(amax, fabs(coeff[kk]));
              if (fabs(coeff[kk])/amax < maxerr)
                  break;

              j1mc = gsl_sf_bessel_Jn(kk, u1);
              j1pc = gsl_sf_bessel_Jn(kk+2, u1);
              if (kind == 1)
              {
                  z2mc = gsl_sf_bessel_Jn(kk, u2);
                  z2pc = gsl_sf_bessel_Jn(kk+2, u2);
              }
              else /* kind = 2 */
              {
                  z2mc = gsl_sf_bessel_Yn(kk, u2);
                  z2pc = gsl_sf_bessel_Yn(kk+2, u2);
              }

              fc = pow(-1.0, 0.5*order+kk+1)*coeff[kk];
              fn += fc*(j1mc*z2pc - j1pc*z2mc);
          }

          fn *= sqrt(pi/2.0)/coeff[0];
      }
      else
      {
          for (kk=0; kk<GSL_SF_MATHIEU_COEFF; kk++)
          {
              amax = GSL_MAX(amax, fabs(coeff[kk]));
              if (fabs(coeff[kk])/amax < maxerr)
                  break;

              j1c = gsl_sf_bessel_Jn(kk, u1);
              j1pc = gsl_sf_bessel_Jn(kk+1, u1);
              if (kind == 1)
              {
                  z2c = gsl_sf_bessel_Jn(kk, u2);
                  z2pc = gsl_sf_bessel_Jn(kk+1, u2);
              }
              else /* kind = 2 */
              {
                  z2c = gsl_sf_bessel_Yn(kk, u2);
                  z2pc = gsl_sf_bessel_Yn(kk+1, u2);
              }

              fc = pow(-1.0, 0.5*(order-1)+kk)*coeff[kk];
              fn += fc*(j1c*z2pc - j1pc*z2c);
          }

          fn *= sqrt(pi/2.0)/coeff[0];
      }

      result_array[ii] = fn;
  } /* order loop */

  return GSL_SUCCESS;
}

/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/

#include "eval.h"

double gsl_sf_mathieu_Mc(int kind, int order, double qq, double zz)
{
	EVAL_RESULT(gsl_sf_mathieu_Mc_e(kind, order, qq, zz, &result));
}

double gsl_sf_mathieu_Ms(int kind, int order, double qq, double zz)
{
	EVAL_RESULT(gsl_sf_mathieu_Ms_e(kind, order, qq, zz, &result));
}

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#ifndef _INC_MAXMULTIMIN_
#define _INC_MAXMULTIMIN_

#include <pthread.h> // for debug info
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_rng.h>

struct estimate_thetas_params;

/**
 * compute the gradient matrix for the length setting theta values
 * dC/dTheta = (C-nugget) * (1/2)*(x_i - x_j)^(alpha) * alpha / (thetaLength)
 *
 * this is a fn-ptr which will be set when the cov function is setup
 * the different target fns are in emulator.c called derivative_l_<covfnname>
 *
 * @param covsub the covariance matrix with the nugget term subtracted
 * @param thetaLength the current value of the length scale we're differentiating wrt
 * @param index, the direction we're looking in
 * @return dCdTheta the matrix of the derivative of C wrt index
 *
 */
void (*makeGradMatLength)(gsl_matrix *dCdTheta,  gsl_matrix* xmodel,
														double thetaLength, int index, int nmodel_points, int nparams);

void maxWithMultiMin(struct estimate_thetas_params *params);

double evalFnMulti(const gsl_vector* theta_vec, void* params_in);

void gradFnMulti(const gsl_vector* theta_vec, void* params_in, gsl_vector * grad_vec);

double getGradientCn(gsl_matrix * dCdtheta, gsl_matrix *cinverse,
										 gsl_vector* training_vector,int nmodel_points, int nthetas);

void evalFnGradMulti(const gsl_vector* theta_vec, void* params,
										 double* fnval, gsl_vector * grad_vec);

double estimateSigma(gsl_matrix* cinverse, void* params_in);

double estimateSigmaFull(gsl_vector *thetas, void* params_in);

int doOptimizeMultiMin( double(*fn)(const gsl_vector*, void*),
												void(*gradientFn)(const gsl_vector*,void*, gsl_vector*),
												void(*fnGradFn)(const gsl_vector*, void*, double*, gsl_vector*),
												gsl_vector *thetaInit, gsl_vector* thetaFinal, void* args);

void set_random_init_value(gsl_rng* rand, gsl_vector* x, gsl_matrix* ranges,int  nthetas);

#endif

max_stars_repo_name: jackdawjackdaw/emulator

max_issues_repo_name: MADAI/MADAIEmulator

max_forks_repo_name: MADAI/MADAIEmulator