text:
/*
** BCM - bipartite connectivity mapping
**
** G.Lohmann, MPI-KYB, Aug 2017
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

#include "viaio/Vlib.h"
#include "viaio/VImage.h"
#include "viaio/mu.h"
#include "viaio/option.h"

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_math.h>

extern void gsl_sort_vector_index(gsl_permutation *,gsl_vector *);

#define SQR(x) ((x) * (x))
#define ABS(x) ((x) > 0 ? (x) : -(x))

#ifdef _OPENMP
#include <omp.h>
#endif /*_OPENMP*/

extern gsl_matrix_float *DataMatrix(VImage *src,int nslices,VImage roi,int);
extern size_t NumVoxels(VImage src);
extern VImage VoxelMap(VImage roi);

typedef struct SpointStruct{
  VShort x;
  VShort y;
  VShort z;
} SPoint;

VDictEntry MetricDict[] = {
  { "pearson", 0 },
  { "spearman", 1 },
  { "MI", 2 },
  { NULL }
};

VImage VCallocImage(int nslices,int nrows,int ncols,VRepnKind repn,VImage ref)
{
  VImage image = VCreateImage(nslices,nrows,ncols,repn);
  if (image == NULL) VError(" error allocating image");
  VFillImage(image,VAllBands,0);
  if (ref != NULL) VCopyImageAttrs (ref,image);
  return image;
}

void ShowImage(gsl_matrix_float *A,char *filename)
{
  int i,j;
  int n=A->size1;
  int m=A->size2;
  VImage tstimage = VCallocImage(1,n,m,VFloatRepn,NULL);
  for (i=0; i <n; i++) {
    for (j=0; j<m; j++) {
      VPixel(tstimage,0,i,j,VFloat) = gsl_matrix_float_get(A,i,j);
    }
  }
  VAttrList out_list = VCreateAttrList();
  VAppendAttr(out_list,"image",NULL,VImageRepn,tstimage);
  FILE *out_file = fopen(filename,"w");
  VWriteFile (out_file, out_list);
}

VImage GetRoi(SPoint addr,int nslices,int nrows,int ncols,int radius)
{
  int i,b,r,c;
  int rad2 = radius*radius+1;

  int b0 = addr.z;
  int r0 = addr.y;
  int c0 = addr.x;

  i=0;
  for (b=b0-radius; b<=b0+radius; b++) {
    for (r=r0-radius; r<=r0+radius; r++) {
      for (c=c0-radius; c<=c0+radius; c++) {
	if (SQR(b-b0) + SQR(r-r0) + SQR(c-c0) > rad2) continue;
	i++;
      }
    }
  }
  VImage roi = VCreateImage(nslices,nrows,ncols,VBitRepn);
  VFillImage(roi,VAllBands,0);
  for (b=b0-radius; b<=b0+radius; b++) {
    if (b < 0 || b >= nslices) VError(" b %d",b);
    for (r=r0-radius; r<=r0+radius; r++) {
      if (r < 0 || r >= nrows) VError(" r %d",r);
      for (c=c0-radius; c<=c0+radius; c++) {
	if (c < 0 || c >= ncols) VError(" c %d",c);
	if (SQR(b-b0) + SQR(r-r0) + SQR(c-c0) > rad2) continue;
	VPixel(roi,b,r,c,VBit) = 1;
      }
    }
  }
  return roi;
}

double Pearson(const float *data1,const float *data2,int n)
{
  int i;
  double corr=0;
  for (i=0; i<n; i++) {
    const double u = (double)data1[i];
    const double v = (double)data2[i];
    corr += u*v;
  }
  corr /= (double)n;
  return corr;
}

double Spearman(const float *data1,const float *data2,int n)
{
  int i;
  double nx = (double)n;
  double kx = nx*(nx*nx-1.0);
  double sxy=0.0;
  for (i=0; i<n; i++) {
    const double u = (double)data1[i];
    const double v = (double)data2[i];
    const double d = (u-v);
    sxy += d*d;
  }
  double rho = 1.0 - 6.0*sxy/kx;
  return rho;
}

/* if X,Y are normalized Gaussians, then I=-0.5*log(1-r^2), with r covariance */
double MutualInformation(const float *data1,const float *data2,int nt)
{
  int i;
  double covar,mi,u;
  double tiny=1.0e-10;

  covar = 0;
  for (i=0; i<nt; i++) {
    const double u = data1[i];
    const double v = data2[i];
    covar += u*v;
  }
  const double nx= nt;
  covar /= nx;

  mi = 0;
  u  = 1.0 - covar*covar;
  if (u > tiny) mi = -0.5 * log(u);
  if (mi < 0) {
    VWarning("MI: %f",mi);
    mi = 0;
  }
  return mi;
}

double Correlation(const float *data1,const float *data2,int n,int metric)
{
  double corr=0.0;
  if (metric == 0) corr = Pearson(data1,data2,n);
  if (metric == 1) corr = Spearman(data1,data2,n);
  if (metric == 2) corr = MutualInformation(data1,data2,n);
  if (gsl_isnan(corr) || gsl_isinf(corr)) corr = 0;
  return corr;
}

void VectorNormalize(gsl_vector_float *x)
{
  int i;
  double nx = (float)x->size;
  double sum1=0,sum2=0;
  for (i=0; i<x->size; i++) {
    sum1 += x->data[i];
    sum2 += x->data[i]*x->data[i];
  }
  double mean = sum1/nx;
  double sigma = sqrt((double)((sum2 - nx * mean * mean) / (nx - 1.0)));

  for (i=0; i<x->size; i++) {
    x->data[i] = (x->data[i] - mean)/sigma;
  }
}

/* bipartite eigenvector centrality mapping */
void VBiadjacencyECM(gsl_matrix_float *C,gsl_vector_float *xa,gsl_vector_float *xb)
{
  int na = C->size1;
  int nb = C->size2;
  int i,iter,maxiter=100;
  double d=0,sum=0,sum_old=0;

  /* ini */
  gsl_vector_float *ya = gsl_vector_float_calloc(na);
  gsl_vector_float *yb = gsl_vector_float_calloc(nb);
  double nx = (double)(na+nb);
  for (i=0; i<na; i++) xa->data[i] = 1.0/nx;
  for (i=0; i<nb; i++) xb->data[i] = 1.0/nx;

  /* iterations */
  for (iter=0; iter<maxiter; iter++) {
    gsl_blas_sgemv (CblasNoTrans,1.0,C,xb,0.0,ya);
    gsl_blas_sgemv (CblasTrans,1.0,C,xa,0.0,yb);

    sum = 0;
    for (i=0; i<na; i++) sum += ya->data[i]*ya->data[i];
    for (i=0; i<nb; i++) sum += yb->data[i]*yb->data[i];
    sum = sqrt(sum);

    for (i=0; i<na; i++) xa->data[i] = ya->data[i]/sum;
    for (i=0; i<nb; i++) xb->data[i] = yb->data[i]/sum;

    d = fabs(sum-sum_old);
    if (iter > 0) fprintf(stderr," %5d   %f\n",(int)iter,d);
    if (d < 1.0e-8 && iter > 3) break;
    if (iter > 5 && sum > sum_old) break;
    sum_old = sum;
  }

  gsl_vector_float_free(ya);
  gsl_vector_float_free(yb);
}

/* bipartite degree centrality mapping */
void VBiadjacencyDegreeMap(gsl_matrix_float *C,gsl_vector_float *xa,gsl_vector_float *xb)
{
  size_t na = C->size1;
  size_t nb = C->size2;
  size_t i;

  gsl_vector_float *ya = gsl_vector_float_calloc(na);
  gsl_vector_float *yb = gsl_vector_float_calloc(nb);
  for (i=0; i<na; i++) xa->data[i] = 1.0/(double)na;
  for (i=0; i<nb; i++) xb->data[i] = 1.0/(double)nb;

  gsl_blas_sgemv (CblasNoTrans,1.0,C,xb,0.0,ya);
  gsl_blas_sgemv (CblasTrans,1.0,C,xa,0.0,yb);

  gsl_vector_float_memcpy(xa,ya);
  gsl_vector_float_memcpy(xb,yb);

  gsl_vector_float_free(ya);
  gsl_vector_float_free(yb);
}

/* project onto ROI */
void VNetworkProjection(gsl_matrix_float *C,gsl_vector_float *x)
{
  size_t na = C->size1;
  size_t nb = C->size2;
  size_t i,n1=0,n2=0;

  if (x->size == na) {
    n1 = na;
    n2 = nb;
  }
  else {
    n1 = nb;
    n2 = na;
  }
  float nx = (float)n2;

  fprintf(stderr," network projection\n");

  /* project bipartite graph onto ROI */
  gsl_matrix_float *A = gsl_matrix_float_calloc(n1,n1);
  if (!A) VError(" err allocating projection network");
  int progress=0;

#pragma omp parallel for schedule(guided) firstprivate(C,A)
  for (i=0; i<n1; i++) {
    size_t j=0,k=0;
    if (i%10 == 0) fprintf(stderr," %d0  of  %lu\r",(int)(++progress),n1);
    for (j=0; j<n1; j++) {

      float w=0;
      for (k=0; k<n2; k++) {
	float u = gsl_matrix_float_get(C,i,k);
	float v = gsl_matrix_float_get(C,j,k);
	w += u*v;
      }
#pragma omp critical
      {
	gsl_matrix_float_set(A,i,j,w/nx);
      }
    }
  }
  fprintf(stderr,"\n");

  /* Hubs of projected network */
  float w=0;
  size_t j=0;
  for (i=0; i<n1; i++) {
    w=0;
    for (j=0; j<n1; j++) {
      w += gsl_matrix_float_get(A,i,j);
    }
    x->data[i] = w/nx;
  }
  gsl_matrix_float_free(A);
}

/* make matrix positive */
void VThresholdMatrix(gsl_matrix_float *C,float threshold)
{
  int i,j;
  size_t nneg=0,npos=0;
  float u=0,tiny=1.0e-6;
  double sum1=0,sum2=0;
  double nx = (double)(C->size1 * C->size2);
  for (i=0; i<C->size1; i++) {
    for (j=0; j<C->size2; j++) {
      u = gsl_matrix_float_get(C,i,j);
      sum1 += u;
      sum2 += u*u;
      if (u < threshold+tiny) nneg++;
      else npos++;
      if (u < threshold+tiny) gsl_matrix_float_set(C,i,j,0.0);
    }
  }
  if (npos < 1) VWarning(" Subthreshold correlations only");
  double mean = sum1/nx;
  double sigma = sqrt((double)((sum2 - nx * mean * mean) / (nx - 1.0)));
  fprintf(stderr," matrix mean,std: %f %f,  npos: %lu, nneg: %lu\n",mean,sigma,npos,nneg);
}

int main (int argc,char *argv[])
{
  static SPoint addr1;
  static SPoint addr2;
  static VString roi1_filename = "";
  static VString roi2_filename = "";
  static VShort type=1;
  static VFloat threshold = 0.0;
  static VShort radius=5;
  static VShort   metric = 0;
  static VBoolean plotboth=FALSE;
  static VBoolean normalize=TRUE;
  static VOptionDescRec  options[] = {
    {"roi1",VStringRepn,1,(VPointer) &roi1_filename,VOptionalOpt,NULL,"ROI 1"},
    {"roi2",VStringRepn,1,(VPointer) &roi2_filename,VOptionalOpt,NULL,"ROI 2"},
    {"seed1",VShortRepn,3,(VPointer) &addr1,VOptionalOpt,NULL,"Voxel address of seed point (x,y,z)"},
    {"seed2",VShortRepn,3,(VPointer) &addr2,VOptionalOpt,NULL,"Voxel address of seed point (x,y,z)"},
    {"type", VShortRepn,1,(VPointer) &type,VOptionalOpt,NULL,"Type, 0:ECM, 1: DCM, 2: project onto ROI1, 3: project onto ROI2"},
    {"threshold", VFloatRepn,1,(VPointer) &threshold,VOptionalOpt,NULL,"Matrix threshold"},
    {"metric",VShortRepn,1,(VPointer) &metric,VOptionalOpt,MetricDict,"Correlation metric"},
    {"radius", VShortRepn,1,(VPointer) &radius,VOptionalOpt,NULL,"Radius around seed voxel"},
    {"both", VBooleanRepn,1,(VPointer) &plotboth,VOptionalOpt,NULL,"Whether to plot results of both ROIs"},
    {"normalize", VBooleanRepn,1,(VPointer) &normalize,VOptionalOpt,NULL,"Whether to normalize results"},
  };
  FILE *out_file=NULL;
  VString in_file=NULL;
  VImage roi1=NULL,roi2=NULL;
  int i;
  char *prg = GetLipsiaName("vbcm");

  /* Parse command line arguments and identify files: */
  VParseFilterCmdX (VNumber (options), options, argc, argv,&in_file,&out_file);
  fprintf(stderr," type= %d\n",type);

  /* read functional data */
  VAttrList list = VReadAttrList(in_file,0L,TRUE,FALSE);
  if (list == NULL) VError(" error reading input file %s",in_file);

  /* get pointers to image data */
  int nrows=0,ncols=0,nt=0;
  int nslices = VAttrListNumImages(list);
  VImage *src = VAttrListGetImages(list,nslices);
  VImageDimensions(src,nslices,&nt,&nrows,&ncols);
  fprintf(stderr," image dims:  %d  %d  %d,  nt: %d\n",nslices,nrows,ncols,nt);

  /* omp-stuff */
#ifdef _OPENMP
  int num_procs=omp_get_num_procs();
  printf("using %d cores\n",(int)num_procs);
  omp_set_num_threads(num_procs);
#endif /* _OPENMP */

  /* ROI 1 */
  if (strlen(roi1_filename) > 0) {
    VAttrList list1 = VReadAttrList(roi1_filename,0L,TRUE,FALSE);
    if (list1 == NULL) VError(" error reading roi file %s",roi1_filename);
    roi1 = VReadImage(list1);
    if (roi1 == NULL) VError(" err reading %s",roi1_filename);
  }

  /* ROI 2 */
  if (strlen(roi2_filename) > 0) {
    VAttrList list2 = VReadAttrList(roi2_filename,0L,TRUE,FALSE);
    if (list2 == NULL) VError(" error reading roi file %s",roi2_filename);
    roi2 = VReadImage(list2);
    if (roi2 == NULL) VError(" err reading %s",roi2_filename);
  }

  /* ROIs from seed voxels */
  if (roi1 == NULL) roi1 = GetRoi(addr1,nslices,nrows,ncols,radius);
  if (roi2 == NULL) roi2 = GetRoi(addr2,nslices,nrows,ncols,radius);

  /* voxel maps */
  VImage map1  = VoxelMap(roi1);
  VImage map2  = VoxelMap(roi2);
  size_t nvox1 = NumVoxels(roi1);
  size_t nvox2 = NumVoxels(roi2);
  fprintf(stderr," nvox1: %lu,  nvox2: %lu\n",nvox1,nvox2);

  /* read data matrices */
  gsl_matrix_float *X1 = DataMatrix(src,nslices,roi1,(int)metric);
  gsl_matrix_float *X2 = DataMatrix(src,nslices,roi2,(int)metric);

  /* correlation matrix */
  fprintf(stderr," CorrMatrix...\n");
  int rad2 = 3*3;     /* exclusion radius */
  gsl_matrix_float *A = gsl_matrix_float_calloc(nvox1,nvox2);
  if (A == NULL) VError(" err allocating corr matrix, %d x %d",nvox1,nvox2);
  int progress=0;

#pragma omp parallel for schedule(guided) firstprivate(A)
  for (i=0; i<nvox1; i++) {
    if (i%100 == 0) fprintf(stderr," %d00\r",(int)(++progress));
    int bi = VPixel(map1,0,0,i,VShort);
    int ri = VPixel(map1,0,1,i,VShort);
    int ci = VPixel(map1,0,2,i,VShort);

    int j=0;
    float u=0;
    float *tmp = (float *) VCalloc(nvox2,sizeof(float));
    for (j=0; j<nvox2; j++) tmp[j] = -1.0;
    const float *data1 = gsl_matrix_float_const_ptr(X1,i,0);
    for (j=0; j<nvox2; j++) {
      int bj = VPixel(map2,0,0,j,VShort);
      int rj = VPixel(map2,0,1,j,VShort);
      int cj = VPixel(map2,0,2,j,VShort);
      int di = SQR(bi-bj) + SQR(ri-rj) + SQR(ci-cj);
      if (di < rad2) continue;

      const float *data2 = gsl_matrix_float_const_ptr(X2,j,0);
      tmp[j] = Correlation(data1,data2,nt,(int)metric);
    }

#pragma omp critical
    {
      for (j=0; j<nvox2; j++) {
	u = tmp[j];
	gsl_matrix_float_set(A,i,j,u);
      }
    }
    VFree(tmp);
  }
  fprintf(stderr,"\n");

  /*
  ShowImage(A,"test.v");
  exit(0);
  */

  /* Biadjacency ECM */
  fprintf(stderr," Biadjacency Connectivity...\n");
  gsl_vector_float *xa = gsl_vector_float_calloc(nvox1);
  gsl_vector_float *xb = gsl_vector_float_calloc(nvox2);

  fprintf(stderr," type= %d\n",type);

  VThresholdMatrix(A,threshold);

  switch(type) {
  case 0:
    VBiadjacencyECM(A,xa,xb);
    break;
  case 1:
    VBiadjacencyDegreeMap(A,xa,xb);
    break;
  case 2:
    VNetworkProjection(A,xa);
    break;
  case 3:
    VNetworkProjection(A,xb);
    break;
  default:
    VError(" unknown type");
  }

  /* normalize */
  if (normalize) {
    VectorNormalize(xa);
    VectorNormalize(xb);
  }

  /* output */
  VImage dest = VCallocImage(nslices,nrows,ncols,VFloatRepn,src[0]);
  int b,r,c;
  for (i=0; i<nvox1; i++) {
    b = VPixel(map1,0,0,i,VShort);
    r = VPixel(map1,0,1,i,VShort);
    c = VPixel(map1,0,2,i,VShort);
    VPixel(dest,b,r,c,VFloat) = 100.0*xa->data[i];
  }
  if (plotboth) {
    for (i=0; i<nvox2; i++) {
      b = VPixel(map2,0,0,i,VShort);
      r = VPixel(map2,0,1,i,VShort);
      c = VPixel(map2,0,2,i,VShort);
      VPixel(dest,b,r,c,VFloat) = 100.0*xb->data[i];
    }
  }

  /* write output to disk */
  VAttrList out_list = VCreateAttrList();
  VAttrList geolist = VGetGeoInfo(list);
  if (geolist != NULL) {
    double *D = VGetGeoDim(geolist,NULL);
    D[0] = 3;  /* 3D */
    D[4] = 1;
    VSetGeoDim(geolist,D);
    VSetGeoInfo(geolist,out_list);
  }
  VHistory(VNumber(options),options,prg,&list,&out_list);
  VAppendAttr(out_list,"image",NULL,VImageRepn,dest);
  VWriteFile (out_file, out_list);
  fclose(out_file);
  fprintf (stderr, "%s: done.               \n", argv[0]);
  exit(0);
}


text:
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <gsl/gsl_multifit.h>

//DEFINES
#define SUCCESS 1
#define FAILURE -1
#define PROGRESS "===================="
#define PRG_WIDTH 20

//FUNCTIONS DECLARATION
int rows_number(char *);
double mean(double *, int);
void cumsum(double *, double *, int);
void slice_vec(double *, double *, int, int);
void polynomialFit(int, int, double *, double *, double *);

//MAIN
int main(int argc, char **argv)
{
    //inputs
    int input_args = 6;
    if(argc < input_args){
        printf("Not enough input arguments!\n");
        return FAILURE;
    }else if(argc > input_args){
        printf("Too many input arguments!\n");
        return FAILURE;
    }
    char file_name[255];
    memset(file_name, 0x00, sizeof(file_name));
    sprintf(file_name, "%s", argv[1]);
    int min_win = atoi(argv[2]);
    int ord = atoi(argv[3]);
    int rev_seg = atoi(argv[4]);
    char path_tot[255];
    memset(path_tot, 0x00, sizeof(path_tot));
    sprintf(path_tot, "%s/dfa_mw=%d_pol=%d_bck=%d.txt", argv[5], min_win, ord, rev_seg);
    //file length
    int N = rows_number(file_name);
    if(N == FAILURE){
        printf("Cannot open input file.\n");
        return FAILURE;
    }
    //time series
    double *pn;
    pn = calloc(N, sizeof(double));
    if(!pn){
        printf("MALLOC ERROR (pn)\n");
        return FAILURE;
    }
    FILE *f;
    f = fopen(file_name, "r");
    if(!f){
        printf("Cannot open input file.\n");
        return FAILURE;
    }
    for(int i = 0; i < N; i++)
        fscanf(f, "%lf", pn+i);
    fclose(f);
    printf("Data length: <%d>\n", N);
    //time vector
    double *t;
    t = calloc(N, sizeof(double));
    if(!t){
        printf("MALLOC ERROR (t)\n");
        return FAILURE;
    }
    for(int i = 0; i < N; i++)
        t[i] = (double)(i+1);
    //time series minus its mean
    double a_ave = mean(pn, N);
    double *pn_nomean;
    pn_nomean = calloc(N, sizeof(double));
    if(!pn_nomean){
        printf("MALLOC ERROR (pn_nomean)\n");
        return FAILURE;
    }
    for(int i = 0; i < N; i++)
        pn_nomean[i] = pn[i] - a_ave;
    //cumulative sum
    double *y;
    y = calloc(N, sizeof(double));
    if(!y){
        printf("MALLOC ERROR (y)\n");
        return FAILURE;
    }
    printf("Data integration...");
    cumsum(pn_nomean, y, N);
    printf("done!\n");
    //defining parameters
    int max_win = 5;
    int end_dfa = N / max_win;
    int *s;
    int range_dfa = end_dfa - min_win + 1;
    s = calloc(range_dfa, sizeof(int));
    if(!s){
        printf("MALLOC ERROR (s)\n");
        return FAILURE;
    }
    for(int i = 0; i < range_dfa; i++)
        s[i] = i + min_win;
    //fluctuations vector and other arrays
    double *F, *F_nu1, *F_nu2, *t_fit, *y_fit, *diff_vec;
    int F_len = N / min_win;
    F = calloc(range_dfa, sizeof(double));
    if(!F){
        printf("MALLOC ERROR (F)\n");
        return FAILURE;
    }
    F_nu1 = calloc(F_len, sizeof(double));
    if(!F_nu1){
        printf("MALLOC ERROR (F_nu1)\n");
        return FAILURE;
    }
    F_nu2 = calloc(F_len, sizeof(double));
    if(!F_nu2){
        printf("MALLOC ERROR (F_nu2)\n");
        return FAILURE;
    }
    t_fit = calloc(end_dfa, sizeof(double));
    if(!t_fit){
        printf("MALLOC ERROR (t_fit)\n");
        return FAILURE;
    }
    y_fit = calloc(end_dfa, sizeof(double));
    if(!y_fit){
        printf("MALLOC ERROR (y_fit)\n");
        return FAILURE;
    }
    diff_vec = calloc(end_dfa, sizeof(double));
    if(!diff_vec){
        printf("MALLOC ERROR (diff_vec)\n");
        return FAILURE;
    }
    //computation
    int start_lim, end_lim;
    double *fit_coeffs;
    fit_coeffs = calloc(ord+1, sizeof(double));
    if(!fit_coeffs){
        printf("MALLOC ERROR (fit_coeffs)\n");
        return FAILURE;
    }
    for(int i = 0; i < range_dfa; i++){
        int N_s = N / s[i];
        double perc = i * 100 / (double)range_dfa;
        int prg = (i * PRG_WIDTH) / range_dfa;
        printf("Computing fluctuations => [%.*s%*s] %.2lf%%\r", prg, PROGRESS, PRG_WIDTH-prg, "", perc);
        fflush(stdout);
        for(int v = 0; v < N_s; v++){
            start_lim = v * s[i];
            end_lim = (v + 1) * s[i] - 1;
            slice_vec(t, t_fit, start_lim, end_lim);
            slice_vec(y, y_fit, start_lim, end_lim);
            polynomialFit(s[i], ord+1, t_fit, y_fit, fit_coeffs);
            for(int j = 0; j < s[i]; j++){
                for(int k = 0; k < ord+1; k++)
                    y_fit[j] -= fit_coeffs[k] * pow(t_fit[j], k);
                diff_vec[j] = pow(y_fit[j], 2.0);
            }
            F_nu1[v] = mean(diff_vec, s[i]);
        }
        if(rev_seg == 1){
            for(int v = 0; v < N_s; v++){
                start_lim = v * s[i] + (N - N_s * s[i]);
                end_lim = (v + 1) * s[i] + (N - N_s * s[i]);
                slice_vec(t, t_fit, start_lim, end_lim);
                slice_vec(y, y_fit, start_lim, end_lim);
                polynomialFit(s[i], ord+1, t_fit, y_fit, fit_coeffs);
                for(int j = 0; j < s[i]; j++){
                    for(int k = 0; k < ord+1; k++)
                        y_fit[j] -= fit_coeffs[k] * pow(t_fit[j], k);
                    diff_vec[j] = pow(y_fit[j], 2.0);
                }
                F_nu2[v] = mean(diff_vec, s[i]);
            }
            F[i] = sqrt((mean(F_nu1, N_s) + mean(F_nu2, N_s)) / (double)2);
        }else{
            F[i] = sqrt(mean(F_nu1, N_s));
        }
    }
    printf("Computing fluctuations => [%s] 100.00%%\r", PROGRESS);
    fflush(stdout);
    printf("\n");
    free(F_nu1); free(F_nu2);
    free(t_fit); free(y_fit);
    free(diff_vec); free(fit_coeffs);
    //HURST EXPONENT
    printf("Computing hurst exponent...");
    double *log_s, *log_F, *H_fit;
    log_s = calloc(range_dfa, sizeof(double));
	if(!log_s){
        printf("MALLOC ERROR (log_s)\n");
        return FAILURE;
    }
    log_F = calloc(range_dfa, sizeof(double));
    if(!log_F){
        printf("MALLOC ERROR (log_F)\n");
        return FAILURE;
    }
    H_fit = calloc(2, sizeof(double));
    if(!H_fit){
        printf("MALLOC ERROR (H_fit)\n");
        return FAILURE;
    }
    for(int i = 0; i < range_dfa; i++){
        log_s[i] = log(s[i]);
        log_F[i] = log(F[i]);
    }
    polynomialFit(range_dfa, 2, log_s, log_F, H_fit);
    printf("done\n");
    printf("Output file...");
    f = fopen(path_tot, "w");
    if(!f){
        printf("Cannot open output file.\n");
        return FAILURE;
    }
    for(int i = 0; i < range_dfa; i++)
        fprintf(f, "%d %lf %lf %lf\n", s[i], F[i], H_fit[0]+log_s[i]*H_fit[1], H_fit[1]);
    fclose(f);
    printf("done\n");
    free(F); free(log_s); free(log_F); free(H_fit);
    return 0;
}

//FUNCTIONS
int rows_number(char *file_name)
{
    FILE *f;
    int stop;
    int lines = 0;
    f = fopen(file_name,"r");
    if(!f){
        printf("Cannot open file %s\n", file_name);
        return FAILURE;
    }
    while(!feof(f)){
        stop = fgetc(f);
        if(stop == '\n')
            lines++;
    }
    fclose(f);
    return lines;
}

double mean(double *vec, int L)
{
    double avg = 0.0;
    for(int i = 0; i < L; i++)
        avg += vec[i];
    avg /= (double)L;
    return avg;
}

void cumsum(double *vec, double *sum_vec, int L)
{
    sum_vec[0] = vec[0];
    for(int i = 1; i < L; i++)
        sum_vec[i] = sum_vec[i-1] + vec[i];
}

void slice_vec(double *all_vec, double *sliced_vec, int start, int end)
{
    for(int i = 0; i <= (end-start); i++)
        sliced_vec[i] = all_vec[start+i];
}

void polynomialFit(int obs, int degree, double *dx, double *dy, double *store)
{
    gsl_multifit_linear_workspace *ws;
    gsl_matrix *cov, *X;
    gsl_vector *y, *c;
    double chisq;
    int i, j;
    //alloc
    X = gsl_matrix_alloc(obs, degree);
    y = gsl_vector_alloc(obs);
    c = gsl_vector_alloc(degree);
    cov = gsl_matrix_alloc(degree, degree);
    //computation
    for(i = 0; i < obs; i++){
        for(j = 0; j < degree; j++)
            gsl_matrix_set(X, i, j, pow(dx[i], j));
        gsl_vector_set(y, i, dy[i]);
    }
    ws = gsl_multifit_linear_alloc(obs, degree);
    gsl_multifit_linear(X, y, c, cov, &chisq, ws);
    for(i = 0; i < degree; i++)
        store[i] = gsl_vector_get(c, i);
    //free memory
    gsl_multifit_linear_free(ws);
    gsl_matrix_free(X);
    gsl_matrix_free(cov);
    gsl_vector_free(y);
    gsl_vector_free(c);
}

ext: c


text:
#pragma once

#include <winrt\Windows.Foundation.h>
#include <d3d11.h>
#include <cstdint>
#include <array>
#include <limits>
#include <gsl\gsl>

namespace Library
{
	class Game;

	class RenderStateHelper final
	{
	public:
		RenderStateHelper(Game& game);
		RenderStateHelper(const RenderStateHelper&) = delete;
		RenderStateHelper(RenderStateHelper&&) = default;
		RenderStateHelper& operator=(const RenderStateHelper&) = delete;
		RenderStateHelper& operator=(RenderStateHelper&&) = default;
		~RenderStateHelper() = default;

		static void ResetAll(gsl::not_null<ID3D11DeviceContext*> deviceContext);
		static void ResetRasterizerState(gsl::not_null<ID3D11DeviceContext*> deviceContext);
		static void ResetBlendState(gsl::not_null<ID3D11DeviceContext*> deviceContext);
		static void ResetDepthStencilState(gsl::not_null<ID3D11DeviceContext*> deviceContext);

		ID3D11RasterizerState* RasterizerState();
		ID3D11BlendState* BlendState();
		ID3D11DepthStencilState* DepthStencilState();

		void SaveRasterizerState();
		void RestoreRasterizerState() const;

		void SaveBlendState();
		void RestoreBlendState() const;

		void SaveDepthStencilState();
		void RestoreDepthStencilState() const;

		void SaveAll();
		void RestoreAll() const;
		void ClearAll();

	private:
		Game& mGame;

		winrt::com_ptr<ID3D11RasterizerState> mRasterizerState;
		winrt::com_ptr<ID3D11BlendState> mBlendState;
		std::array<float, 4> mBlendFactor;
		std::uint32_t mSampleMask{ std::numeric_limits<std::uint32_t>::max() };
		winrt::com_ptr<ID3D11DepthStencilState> mDepthStencilState;
		std::uint32_t mStencilRef{ std::numeric_limits<std::uint32_t>::max() };
	};
}

ext: h


text:
#ifndef UTILSH
#define UTILSH

#include <gsl/gsl_vector.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_permutation.h>
#include <math.h>
#include <time.h>
#include <assert.h>
#include <sys/stat.h>
#include <sys/types.h>
#include <stdlib.h>
#include <stdarg.h>

time_t TIME;
char line1[80];

void outlog(char* fmt, ...);

typedef struct int_vector
{
    int size;
    int* val;
} int_vector;

int_vector* new_int_vector(int size);
void delete_int_vector(int_vector* iv);

void ivappend(int_vector* v, int val);

static inline int ivget(int_vector* v, int n)
{ return(v->val[n]); };

static inline void ivset(int_vector* v, int n, int val)
{ v->val[n] = val; };

static inline void vset(gsl_vector* v, int i, double x)
{ gsl_vector_set(v, i, x); };

static inline double vget(const gsl_vector* v, int i)
{ return(gsl_vector_get(v, i)); };

static inline void vinc(gsl_vector* v, int i, double x)
{ vset(v, i, vget(v, i) + x); };

static inline double log_sum(double log_a, double log_b)
{
  if (log_a < log_b)
      return(log_b+log(1 + exp(log_a-log_b)));
  else
      return(log_a+log(1 + exp(log_b-log_a)));
};

static inline double lgam(double x)
{ return(gsl_sf_lngamma(x)); }

double runif();

double rgauss(double mean, double stdev);

int sample_from_log(gsl_vector* log_prob);

void init_random_number_generator();

void vct_fscanf(const char* filename, gsl_vector* v);

void iv_permute(int_vector* iv);
void iv_permute_from_perm(int_vector* iv, gsl_permutation* p);
int_vector* iv_copy(int_vector* iv);

double sum(gsl_vector* v);

void print_vector(gsl_vector* v);

gsl_permutation* rpermutation(int size);

void write_vect(gsl_vector* vect, char* name, FILE* file);

gsl_vector* read_vect(char* name, int size, FILE* file);

int read_int(char* name, FILE* file);

void write_int(int x, char* name, FILE* file);

void write_double(double x, char* name, FILE* file);

double read_double(char* name, FILE* file);

int directory_exist(const char *dname);

void make_directory(char* name);

int rbernoulli(double p);

double rbeta(double a,
             double b);

double rgamma(double shape,
              double scale);

double log_dgamma(double x,
              double shape,
              double scale);

void resize(gsl_vector * vec, size_t newsize);

#endif

max_stars_repo_name: thatindiandude/PGM-Project

max_issues_repo_name: thatindiandude/PGM-Project

max_forks_repo_name: thatindiandude/PGM-Project


text:
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <string.h>
#include <math.h>
#include <assert.h>
#include <alloca.h>

#include <gsl/gsl_cdf.h>

#include "limits.h"
#include "stattest.h"
#include "cat.h"
#include "fisher.h"
#include "min2.h"
#include "bvr.h"

typedef unsigned int count_t;
typedef const count_t COUNT_T;
typedef double       prob_t;

/**
  * This class is stateful in the interest of efficiency
  * because some methods share preliminary computations:
  * 1) statistics require expectation, and
  * 2) expectation requires marginals
  * (which may also be used incidentally for other purposes), so this
  * variable captures the progress towards complete calculations
  * to enable implicit just-in-time calculation.
  */
enum State {
	Nothing = 0,
	Marginals,
	Expectation
};

struct CatCovars {

	/**
	 * Buffers are sized to accomodate these dimensions at construction
	 * and are not changeable thereafter.
	 */
	unsigned int ROW_CAPACITY, COL_CAPACITY;

	/**
	  * "decl(ared)_" because the bounds says nothing about the row/col
	  * content--empty rows/columns may exist! I'm just driving this point
	  * home in these members' names because it has already burned me once!
	  */
	unsigned int decl_rows, decl_cols, decl_cells;

	/**
	  * Count of samples pushed since last cat_clear.
	  */
	count_t sample_count;

	// Calculation state ///////////////////////////////////////////////////

	prob_t minimumExpected;

	/**
	 * This is used strictly by badCellExists
	 */
	count_t REQUIRED_CELL_MINIMUM;

	enum State calculated;

	// Buffers /////////////////////////////////////////////////////////////

	size_t SIZEOF_PROB_BUF;

	prob_t  *expect;

	size_t SIZEOF_COUNT_BUF;

	/**
	  * The counts matrix is stored packed against the head of this
	  * (heap-allocated) buffer. In other words, though this buffer
	  * is -expected- to be larger than required (since the class
	  * is designed to be iteratively re-used), only the first
	  * (decl_rows*decl_cols) of the buffer are used for any given
	  * matrix.
	  */
	count_t *counts;

	/**
	  * Following are derivative data...computed from counts.
	  */
	count_t *rmarg;
	count_t *cmarg;
};

static void _coordsOfOffset(
	struct CatCovars *co, unsigned int off, unsigned int *r, unsigned int *c ) {
	*r = off / co->decl_cols;
	*c = off % co->decl_cols;
}

static count_t _count( struct CatCovars *co, coord_t r, coord_t c ) {
	const int C = co->decl_cols;
	assert( r < co->decl_rows && c < C );
	return co->counts[ r*C + c ];
}

static prob_t _expected( struct CatCovars *co, coord_t r, coord_t c ) {
	const int C = co->decl_cols;
	assert( r < co->decl_rows && c < C );
	return co->expect[ r*C + c ];
}

static void _recalcSampleCount( struct CatCovars *co ) {
	int n = co->decl_cells;
	co->sample_count = 0;
	while( n-- > 0 ) co->sample_count += co->counts[n];
}

/**
  * This starts searching at the given coordinates and continues linearly
  * through the matrix until it finds a violating cell at which point it
  * turns the linear index into 2D coordinates.
  */
static bool _badCellExists( struct CatCovars *co,
		unsigned int *offset ) {

	unsigned int i = *offset;
	while( i < co->decl_cells ) {
		if( co->counts[i] < co->REQUIRED_CELL_MINIMUM ) {
			*offset = i;
			return true;
		}
		i++;
	}
	return false;
}

/**
  * If after accumulation there are fewer than 2 rows or fewer than 2
  * columns, the table is degenerate. Rows or columns with zero marginal
  * sum are OK! For a contingency table that's bad but not -degenerate-.
  */
static bool _immediatelyDegenerate( struct CatCovars *co ) {
	return ! ( co->decl_rows >= 2 && co->decl_cols >= 2 );
}

#if 0
static bool _eventuallyDegenerate( struct CatCovars *co ) {
	unsigned nonemprows = 0;
	unsigned nonempcols = 0;
	unsigned i;
	for(i = 0; i < co->decl_rows; i++ ) if( co->rmarg[i] > 0 ) nonemprows++;
	for(i = 0; i < co->decl_cols; i++ ) if( co->cmarg[i] > 0 ) nonempcols++;
	return  !  ( nonemprows >= 2 && nonempcols >= 2 );
}
#endif

static void _cullRow( struct CatCovars *co, unsigned int r ) {

	// An over-writing copy suffices since matrix is
	// stored row-major. mergeCols is more intricate.

	if( r+1 < co->decl_rows ) { // ...so 0 < rows-U-1
		memmove(
			co->counts + (r  )*co->decl_cols,
			co->counts + (r+1)*co->decl_cols,
			(co->decl_rows - r - 1)*co->decl_cols*sizeof(count_t) );
	} else {
		// ...it's the last row, and there's nothing to move.
	}
	co->decl_rows  -= 1;
	co->decl_cells -= co->decl_cols;

	co->sample_count = 0;
	co->calculated = Nothing; // ...to trigger recomputation.
}

static void _cullCol( struct CatCovars *co, unsigned int c ) {

	const unsigned int N = co->decl_cells;
	// Shrink. Because matrices are stored row-major and calculations
	// throughout this class assume the matrix is packed at the head of
	// the buffer, there is actually a fairly complex shifting of all
	// but the first c entries that has to occurr...
	// Basically, traverse the matrix linearly (storage order) shifting
	// "left" all entries that are NOT in the c'th column...
	for(unsigned int src = 0, dst = 0 ; src < N; src++ ) {
		if( src % co->decl_cols != c ) {
			co->counts[dst++] = co->counts[src];
		}
	}
	co->decl_cols -= 1;
	co->decl_cells -= co->decl_rows;

	co->sample_count = 0;
	co->calculated = Nothing; // ...to trigger recomputation.
}

/**
  * Following two MUST be called in order before any of
  * the statistics are computed. In other words,
  *    !!!! THESE ARE NOT CALLED IMPLICITLY! !!!!
  * (Statistical methods which -ought- to be const can't be unless
  *  'expect', '[rc]marg' and friends are made mutable->)
  */
static void _calc_marginals( struct CatCovars *co ) {

	const int R = co->decl_rows;
	const int C = co->decl_cols;
	co->sample_count = 0;

	memset( co->rmarg, 0, R*sizeof(count_t) );
	memset( co->cmarg, 0, C*sizeof(count_t) );

	// Calculate the marginals...

	for(unsigned int i = 0; i < R; i++ ) {
		for(unsigned int j = 0; j < C; j++ ) {
			const count_t C = _count( co, i,j);
			co->rmarg[i] += C;
			co->cmarg[j] += C;
			co->sample_count    += C;
		}
	}

	co->calculated = Marginals;
}

/**
 * Calculates the expected probabilities of a contingency table assuming
 * independence of the two variables.
 * "Since the Chi-square and G^2 statistics rest on the normal approximation
 *  to the hypergeometric distribution, these approximations get strained
 *  when there are zero frequency cells or when the expected counts are
 *  'very small'." ...from Yaramakala's Masters Thesis, p.43
 *
 * "...for an independence test to be reliable, all cells have non-zero
 *  expected values and at least 80% of the cells should have expected
 *  values greater than 5." [Everitt 1977, Upton 1978]
 */
static void _calc_expectation( struct CatCovars *co ) {

	const int R = co->decl_rows;
	const int C = co->decl_cols;

	if( co->calculated < Marginals )
		_calc_marginals( co );

	if( co->sample_count > 0 ) {
		co->minimumExpected
			= (prob_t)(co->rmarg[0] * co->cmarg[0])
			/ (prob_t)co->sample_count;
		for(unsigned int i = 0; i < R; i++ ) {
			for(unsigned int j = 0; j < C; j++ ) {
				const prob_t E
					= (prob_t)(co->rmarg[i] * co->cmarg[j])
					/ (prob_t)co->sample_count;
				co->expect[i*co->decl_cols + j] = E;
				if( co->minimumExpected > E )
					co->minimumExpected = E;
			}
		}
	}

	co->calculated = Expectation;
}

/**
 * This function cannot fail. Valid results always exist.
 * Returns the indices of the minimal row and columns marginals.
 */
#if 0
static void _minimalMarginals( struct CatCovars *co, unsigned int *rm, unsigned int *cm ) {

	unsigned int MIN_DIM
		= co->decl_rows < co->decl_cols ? co->decl_rows : co->decl_cols;
	unsigned int i, ri = 0, ci = 0;
	count_t mr, mc;

	if( co->calculated < Marginals )
		_calc_marginals( co );

	mr = co->rmarg[ri];
	mc = co->cmarg[ci];

	// Iterate over marginals at indices [0,min(rows,cols) )

	for(i = 0; i < MIN_DIM; i++ ) {
		if( mr > co->rmarg[i] ) {
			mr = co->rmarg[i];
			ri = i;
		}
		if( mc > co->cmarg[i] ) {
			mc = co->cmarg[i];
			ci = i;
		}
	}

	// ...continue iteration over [ min(rows,cols), max(rows,cols) ).

	if( co->decl_cols > co->decl_rows )
		for(; i < co->decl_cols; i++ ) {
			if( mc > co->cmarg[i] ) {
				mc = co->cmarg[i];
				ci = i;
			}
		}
	else
	if( co->decl_rows > co->decl_cols )
		for(; i < co->decl_rows; i++ ) {
			if( mr > co->rmarg[i] ) {
				mr = co->rmarg[i];
				ri = i;
			}
		}

	*rm = ri;
	*cm = ci;
}
#endif

/**
 * A pretty printer...primarily for debugging.
 */
#ifdef _DEBUG
static void dbg_dump( struct CatCovars *co, FILE *fp, const char *prefix ) {

	if( prefix==NULL ) prefix="";

	fprintf( fp, "%s %d(r) x %d(c) (capacity = %dx%d)\n",
		prefix,
		co->decl_rows,
		co->decl_cols,
		co->ROW_CAPACITY,
		co->COL_CAPACITY );

	fprintf( fp, "%s cell counts...\n", prefix );
	for(unsigned int r = 0; r < co->decl_rows; r++ ) {
		fputs( prefix, fp );
		for(unsigned int c = 0; c < co->decl_cols; c++ ) {
			fprintf( fp, "\t%d", _count(co,r,c) );
		}
		fputc( '\n', fp );
	}

	if( co->calculated >= Expectation ) {
		fprintf( fp, "%s %d total samples, expectation...\n", prefix, co->sample_count );
		for(unsigned int r = 0; r < co->decl_rows; r++ ) {
			fputs( prefix, fp );
			for(unsigned int c = 0; c < co->decl_cols; c++ ) {
				fprintf( fp, "\t%.3e", _expected( co, r, c) );
			}
			fputc( '\n', fp );
		}
	}
}
#endif

/***************************************************************************
  * Publics
  */

void cat_destroy( void *pv ) {

	if( pv ) {
		struct CatCovars *co = (struct CatCovars *)pv;
		if( co->expect )
			free( co->expect );
		free( pv );
	}
}

/**
  * Pre-allocate a set of working buffers large enough for all anticipated
  * calculations (max feature length) and a struct to wrap them.
  */
void *cat_create( unsigned int rcap, unsigned int ccap ) {

	struct CatCovars *co
		= calloc( 1, sizeof(struct CatCovars) );

	if( co ) {
		co->ROW_CAPACITY = rcap;
		co->COL_CAPACITY = ccap;
		co->SIZEOF_COUNT_BUF
			= (rcap*ccap+rcap+ccap)*sizeof(count_t);
		co->SIZEOF_PROB_BUF
			= rcap*ccap*sizeof(prob_t);

		// Allocate one large buffer and partition it up.
		co->expect
			= calloc( co->SIZEOF_PROB_BUF + co->SIZEOF_COUNT_BUF, sizeof(char) );
		co->counts
			= (count_t*)( co->expect + rcap*ccap );
		co->rmarg
			= co->counts + rcap*ccap;
		co->cmarg
			= co->rmarg + rcap;

		co->REQUIRED_CELL_MINIMUM = 5;

		// If -anything- failed clean up any successes.
		if( (NULL == co->counts) ||
			(NULL == co->expect) ) {
			cat_destroy( co );
			return NULL;
		}
		return co;
	}
	return NULL;
}

void cat_clear( void *pv, coord_t nr, coord_t nc ) {

	struct CatCovars *co = (struct CatCovars *)pv;

	assert( nr <= co->ROW_CAPACITY && nc <= co->COL_CAPACITY );

	// Set dimensions and cell count.

	co->decl_rows  = nr;
	co->decl_cols  = nc;
	co->decl_cells = nr*nc;

	co->sample_count    = 0;
	co->minimumExpected = 0.0;
	co->calculated = Marginals;

	// OPTIMIZATION: Only clearing the part of capacity intended for use!

	memset( co->counts, 0, co->decl_cells*sizeof(count_t) );
	memset( co->rmarg,  0, co->decl_rows*sizeof(count_t) );
	memset( co->cmarg,  0, co->decl_cols*sizeof(count_t) );
	memset( co->expect, 0, co->decl_cells*sizeof(prob_t) );
}

void cat_setMinCellCount( void *pv, unsigned n ) {
	struct CatCovars *co = (struct CatCovars *)pv;
	co->REQUIRED_CELL_MINIMUM = n;
}

/**
 * Only this one is public in order that we may count the samples
 * as they're added...
 */
void cat_push( void *pv, coord_t r, coord_t c ) {

	struct CatCovars *co = (struct CatCovars *)pv;

	assert( r < co->decl_rows );
	assert( c < co->decl_cols );

	co->counts[ r*co->decl_cols + c ] += 1;

	// size and marginals are kept up to date...
	co->sample_count++;
	co->rmarg[r] += 1;
	co->cmarg[c] += 1;
}

size_t cat_size( void *pv ) {
	struct CatCovars *co = (struct CatCovars *)pv;
	return co->sample_count;
}

bool cat_complete( void *pv ) {
	struct CatCovars *co = (struct CatCovars *)pv;
	return ! _immediatelyDegenerate( co );
}

bool cat_is2x2( void *pv ) {
	struct CatCovars *co = (struct CatCovars *)pv;
	return (2 == co->decl_rows) && (2 == co->decl_cols);
}

/**
  * <bits> will ALWAYS be non-zero on entry, but it may only contain
  * a single set bit.
  */
static int _offsetOfLeastLoss( unsigned bits, COUNT_T *loss, count_t *leastLoss ) {
	int v = __builtin_ctz( bits );
	int l = loss[v];
	bits ^= (1<<v);
	while( bits ) {
		const int i = __builtin_ctz( bits );
		bits ^= (1<<i);
		if( l > loss[i] ) {
			l = loss[i];
			v = i;
		}
	}
	*leastLoss = l;
	return v;
}

/**
 * This removes all rows and any columns containing cells with
 * counts LESS THAN OR EQUAL TO REQUIRED_CELL_MINIMUM.
 *
 * This method assumes we start with a table larger in at least
 * one dimention than 2x2 and halts when/if we reach a 2x2.
 */
unsigned int cat_cullBadCells( void *pv, char *log, int buflen ) {

	unsigned int culled = 0;
	struct CatCovars *co = (struct CatCovars *)pv;
	const char * const EOL = log + buflen - 2; // leave room for "+\0".
	// BEGIN log maintenance...
	const int MIN_LOG_RECORD_LEN = 3; // "R99"
	// ...since tables will never exceed 99 rows or columns.
	char *pc = log;
	// ...END log maintenance.
	unsigned int victim;
	unsigned int raw_mem_offset = 0; // linear offsets of "bad" cell.

#ifdef _DEBUG
	const bool SHOW_CULLING = getenv("SHOW_CULLING") != NULL;
	if( SHOW_CULLING ) {
		fprintf( stderr, "CULL log: original matrix\n" );
		dbg_dump( co, stderr, "CULL" );
	}
#endif

	assert( ! _immediatelyDegenerate( co ) ); // so BOTH dims >= 2

	while( (co->decl_rows > 2 || co->decl_cols > 2) && _badCellExists( co, &raw_mem_offset ) ) {

		unsigned int r, c;
		char cullty = '?';

		/**
		  * It is always necessary to cull a row or column.
		  * If |columns| > 3, calculate "best" column to cull.
		  * If |rows| > 3, calculate "best" row to cull.
		  * "Best" victim is row OR column among the candidates for
		  * culling with minimal marginal.
		  */

		unsigned row_bitfield = 0;
		unsigned col_bitfield = 0;
		count_t minRowCost, minColCost;

		assert( sizeof(row_bitfield)*8 >= MAX_CATEGORY_COUNT );

		// Find ALL rows and columns that are candidates for culling by
		// virtue of containing a bad cell.

		do {
			_coordsOfOffset( co, raw_mem_offset, &r, &c );
			row_bitfield |= (1<<r);
			col_bitfield |= (1<<c);
			++ raw_mem_offset;
			// Add the row and column to the list of candidate victims.
		} while( _badCellExists( co, &raw_mem_offset ) );

		if( co->calculated < Marginals )
			_calc_marginals( co );

		if( co->decl_rows > 2 ) {

			r = _offsetOfLeastLoss( row_bitfield, co->rmarg, &minRowCost );

			if( co->decl_cols > 2 ) { // EITHER a col OR row may be culled

				c = _offsetOfLeastLoss( col_bitfield, co->cmarg, &minColCost );
				cullty = minRowCost < minColCost ? 'R' : 'C';

			} else                    // ONLY a row may be culled
				cullty = 'R';

		} else {                      // ONLY a column may be culled

			c = _offsetOfLeastLoss( col_bitfield, co->cmarg, &minColCost );
			cullty = 'C';
		}

		if( cullty == 'R' )
			_cullRow( co, (victim = r) );
		else {
			assert( cullty == 'C' );
			_cullCol( co, (victim = c) );
		}

		// Log it, if possible...

		if( pc ) {

			pc += sprintf( pc, "%c%d", cullty, victim );

			// Following is a bit of overkill to make sure we can
			// give some indication that we ran out of log space.

			if( EOL - pc < MIN_LOG_RECORD_LEN ) {
				if( EOL - pc > 0 ) strcpy( pc, "+");
				pc = NULL; // no more logging.
				// Really this should never be executed if log buffer
				// is always sufficient.
			}
		}

		// Unlike merging, decimating CAN leave a 0 cell at a lower
		// linear index, so reset is required...
		culled++;
		raw_mem_offset = 0;

#ifdef _DEBUG
		if( SHOW_CULLING ) {
			fprintf( stderr,
				"CULL After %d selected %c%d (from R:%08X,C:%08X) costing %d samples\n",
				culled, cullty, victim,
				row_bitfield, col_bitfield,
			   	cullty=='R' ? minRowCost : minColCost );
			dbg_dump( co, stderr, "CULL" );
		}
#endif
	}

	if( culled > 0 )
		_recalcSampleCount( co );

	return culled;
}

#ifdef HAVE_G_TEST
int cat_g( void *pv, struct Statistic *result ) {

	struct CatCovars *co = (struct CatCovars *)pv;
	const int R = co->decl_rows;
	const int C = co->decl_cols;
	double g = 0.0;

	abort(); // unfinished.

	if( co->calculated < Expectation )
		_calc_expectation( co );

	for(unsigned int i = 0; i < R; i++ ) {
		for(unsigned int j = 0; j < C; j++ ) {
			const double O = _count( co, i, j );
			const double E = _expected( co, i, j );
			if( E > 0 && O > 0 )
				g += ( O * log( O / E ) );
		}
	}

	g *= 2.0; // TODO: verify!

	result->name
		= "G";
	result->sample_count
		= co->sample_count;
	result->probability
		= gsl_cdf_chisq_Q( g, (R-1)*(C-1) );

	return 0;
}
#endif

int cat_chi_square( void *pv, struct Statistic *result ) {

	struct CatCovars *co = (struct CatCovars *)pv;
	const int R = co->decl_rows;
	const int C = co->decl_cols;

	unsigned int n_empty = 0;
	double chi = 0.0;

	if( co->calculated < Expectation )
		_calc_expectation( co );

	for(unsigned int i = 0; i < R; i++ ) {
		for(unsigned int j = 0; j < C; j++ ) {
			const double O = _count( co, i, j );
			const double E = _expected( co, i, j );
			if( E > 0 && O > 0 )
				chi += ((O-E)*(O-E))/E;
			else
			if( O == 0 )
				n_empty++;
		}
	}

	result->name
		= "Chi-square";
	result->sample_count
		= co->sample_count;
	result->probability
		= gsl_cdf_chisq_Q( chi, (R-1)*(C-1) );

	result->extra_value[0] = R;
	result->extra_value[1] = C;
	result->extra_value[2] = co->minimumExpected;
	result->extra_value[3] = n_empty;

	return 0;
}

int cat_fisher_exact( void *pv, struct Statistic *result ) {

	struct CatCovars *co = (struct CatCovars *)pv;

	// OPTIMIZE: does not make sense to call fisher exact without
	// NULL stest_t argument. Keeping this signature for consistency
	// with all other stats, but it might make sense to change the
	// signature of just this function to return P-value (unless I
	// ever added log-likelihood ratio calculation which is what
	// ought to be returned as the function value).

	/**
	  * Given the table:
	  *    a | b
	  *    --+--
	  *    c | d
	  *
	  * ...we need to pass fexact_prob( x, m, n, k )
	  * m == a+c
	  * n == b+d
	  * k == a+b
	  * x == a
	  */
	result->name
		= "Fisher_Exact";
	result->sample_count
		= co->sample_count;
	result->probability
		= fexact_prob(
			_count(co,0,0),
			_count(co,0,0) + _count(co,1,0),   // a+c
			_count(co,0,1) + _count(co,1,1),   // b+d
			_count(co,0,0) + _count(co,0,1) ); // a+b

	// TODO: Don't need these now.
	result->extra_value[0] = co->decl_rows;
	result->extra_value[1] = co->decl_cols;

	return 0;
}

/**
  */
double cat_spearman_rho( void *pv ) {

	struct CatCovars *co = (struct CatCovars *)pv;
	double rho = nan("nan");
	const unsigned r0 = _count(co,0,0) + _count(co,0,1);
	const unsigned r1 = _count(co,1,0) + _count(co,1,1);
	const unsigned c0 = _count(co,0,0) + _count(co,1,0);
	const unsigned c1 = _count(co,0,1) + _count(co,1,1);

	// Continuing is pointless if any of these are 0.

	if( r0 > 0 && r1 > 0 && c0 > 0 && c1 > 0 ) {

		const float MU    = (co->sample_count+1.0)/2.0;
		const float r0_mu = (r0+1.0)/2.0;
		const float r1_mu = MEAN_RANK_OF_TIES( r0, co->sample_count );
		const float c0_mu = (c0+1.0)/2.0;
		const float c1_mu = MEAN_RANK_OF_TIES( c0, co->sample_count );

		/**
		  * Following is a heavily compacted form of the Pearson rho,
		  * compacted because thanks to all the ties there are only two
		  * "versions" of (x_i-x_mu) and (y_i-x_mu), and the general
		  * sums in the Pearson quotient just become multiples of these.
		  */
		const double DENOM
			= sqrt(
			  (r0*(r0_mu-MU)*(r0_mu-MU) + r1*(r1_mu-MU)*(r1_mu-MU))
			* (c0*(c0_mu-MU)*(c0_mu-MU) + c1*(c1_mu-MU)*(c1_mu-MU)) );

		if( isnormal(DENOM) ) {
			rho =
				( _count(co,0,0)*(r0_mu-MU)*(c0_mu-MU)
				+ _count(co,0,1)*(r0_mu-MU)*(c1_mu-MU)
				+ _count(co,1,0)*(r1_mu-MU)*(c0_mu-MU)
				+ _count(co,1,1)*(r1_mu-MU)*(c1_mu-MU) )
				/
				DENOM;
		}

	}
	return rho;
}

#ifdef _UNITTEST_CAT_

#include <err.h>

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

	if( argc >= 3 ) {

		const int MERGE_LOG_LEN = 128;
		char merge_log[ MERGE_LOG_LEN ];

		struct Statistic result;

		const unsigned int MAXROW
			= atoi(argv[1]);
		const unsigned int MAXCOL
			= atoi(argv[2]);

		FILE *fp
			= argc > 3
			? fopen( argv[3], "r" )
			: stdin;

		int err = 0;
		char *line = NULL;
		size_t n = 0;

		void *accum
			= cat_create( MAXROW, MAXCOL);

		cat_setMinCellCount( accum, 5 );
		cat_clear( accum, MAXROW, MAXCOL );   // ...as per command line.

		while( getline( &line, &n, fp ) > 0 ) {
			int cat[2];
			if( line[0] == '#' ) {
				fputs( line, stdout );
				continue;
			}
			if( 2 == sscanf( line, "%d\t%d\n", cat+0, cat+1 ) ) {
				if( 0 <= cat[0] && 0 <= cat[1] ) {
					cat_push( accum, cat[0], cat[1] );
				}
			} else {
				fprintf( stderr, "Failure parsing line: %s", line );
			}
		}
		free( line );
		fclose( fp );

		puts( "Original..." );
		dbg_dump( accum, stdout, ""  );

		if( cat_cullBadCells( accum, merge_log, MERGE_LOG_LEN ) > 0 ) {
			fprintf( stdout, "Culled %s, leaving...\n", merge_log );
			dbg_dump( accum, stdout, "" );
		}

		memset( &result, 0, sizeof(struct Statistic) );
		err = cat_is2x2( accum )
			? cat_fisher_exact( accum, &result )
			: cat_chi_square( accum, &result );

		if( ! err ) {
			printf( "p-value=%f", result.probability );
			if( cat_is2x2( accum ) )
				printf( ", spearman=%f\n", cat_spearman_rho( accum ) );
			else
				fputc( '\n', stdout );
		} else
			printf( "error\n" );

		cat_destroy( accum );

	} else
		err( -1, "%s <max row> <max col> [ <2-column file> ]", argv[0] );
	return 0;
}
#endif

max_stars_repo_name: IlyaLab/kramtools

max_issues_repo_name: IlyaLab/kramtools

max_forks_repo_name: IlyaLab/kramtools


text:
/*
GENETIC - A simple genetic algorithm.

Copyright 2014, Javier Burguete Tolosa.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

1. Redistributions of source code must retain the above copyright notice, this
	list of conditions and the following disclaimer.

2. Redistributions in binary form must reproduce the above copyright notice,
	this list of conditions and the following disclaimer in the documentation
	and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY Javier Burguete Tolosa ``AS IS'' AND ANY EXPRESS
OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
EVENT SHALL Javier Burguete Tolosa OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

/**
 * \file selection.c
 * \brief Source file to define parents selection functions.
 * \author Javier Burguete Tolosa.
 * \copyright Copyright 2014 Javier Burguete Tolosa. All rights reserved.
 */
#define _GNU_SOURCE
#include <math.h>
#include <gsl/gsl_rng.h>
//#include "config.h"
#include "entity.h"
#include "population.h"
#include "selection.h"

void (*selection_mutation) (Population *, Entity **, gsl_rng *);
  ///< Pointer to the function to apply mutation selection operations.
void (*selection_reproduction) (Population *, Entity **, Entity **, gsl_rng *);
  ///< Pointer to the function to apply reproduction selection operations.
void (*selection_adaptation) (Population *, Entity **, gsl_rng *);
  ///< Pointer to the function to apply adaptation selection operations.

/**
 * Function to select a single random entity.
 */
static void
selection_mutation_random (Population * population,     ///< Population.
                           Entity ** mother,
                           ///< Pointer to the mother entity.
                           gsl_rng * rng)
                           ///< GSL random numbers generator.
{
  *mother = population->entity
    + gsl_rng_uniform_int (rng, population->nsurvival);
}

/**
 * Function to select the best of two random entities.
 */
static void
selection_mutation_bestof2 (Population * population,    ///< Population.
                            Entity ** mother,
                            ///< Pointer to the mother entity.
                            gsl_rng * rng)
                            ///< GSL random numbers generator.
{
  unsigned int i, j;
  i = gsl_rng_uniform_int (rng, population->nsurvival);
  j = gsl_rng_uniform_int (rng, population->nsurvival);
  if (j < i)
    i = j;
  *mother = population->entity + i;
}

/**
 * Function to select the best of three random entities.
 */
static void
selection_mutation_bestof3 (Population * population,    ///< Population.
                            Entity ** mother,
                            ///< Pointer to the mother entity.
                            gsl_rng * rng)
                            ///< GSL random numbers generator.
{
  unsigned int i, j;
  i = gsl_rng_uniform_int (rng, population->nsurvival);
  j = gsl_rng_uniform_int (rng, population->nsurvival);
  if (j < i)
    i = j;
  j = gsl_rng_uniform_int (rng, population->nsurvival);
  if (j < i)
    i = j;
  *mother = population->entity + i;
}

/**
 * Function to select the best entity only.
 */
static void
selection_mutation_best (Population * population,       ///< Population.
                         Entity ** mother,
                         ///< Pointer to the mother entity.
                         gsl_rng * rng __attribute__ ((unused)))
  ///< GSL random numbers generator.
{
  *mother = population->entity;
}

/**
 * Function to select an entity based on linear probability distribution with
 *   respect to rank.
 */
static void
selection_mutation_linearrank (Population * population, ///< Population.
                               Entity ** mother,
                               ///< Pointer to the mother entity.
                               gsl_rng * rng)
                               ///< GSL random numbers generator.
{
  unsigned int i;
  i = (1.0 - sqrt (gsl_rng_uniform (rng))) * population->nsurvival;
  *mother = population->entity + i;
}

/**
 * Function to select a pair of random entities.
 */
static void
selection_reproduction_random (Population * population, ///< Population.
                               Entity ** mother,
                               ///< Pointer to the mother entity.
                               Entity ** father,
                               ///< Pointer to the father entity.
                               gsl_rng * rng)
                               ///< GSL random numbers generator.
{
  *mother = population->entity
    + gsl_rng_uniform_int (rng, population->nsurvival);
  do
    *father = population->entity
      + gsl_rng_uniform_int (rng, population->nsurvival);
  while (father == mother);
}

/**
 * Function to select a pair of entities. For each parent select the best of two
 *   random entities.
 */
static void
selection_reproduction_bestof2 (Population * population,        ///< Population.
                                Entity ** mother,
                                ///< Pointer to the mother entity.
                                Entity ** father,
                                ///< Pointer to the father entity.
                                gsl_rng * rng)
                                ///< GSL random numbers generator.
{
  unsigned int i, j, k;
  i = gsl_rng_uniform_int (rng, population->nsurvival);
  j = gsl_rng_uniform_int (rng, population->nsurvival);
  if (i < j)
    i = j;
  do
    j = gsl_rng_uniform_int (rng, population->nsurvival);
  while (i == j);
  k = gsl_rng_uniform_int (rng, population->nsurvival);
  if (i != k && k < j)
    j = k;
  *mother = population->entity + i;
  *father = population->entity + j;
}

/**
 * Function to select a pair of entities. For each parent select the best of
 *   three random entities.
 */
static void
selection_reproduction_bestof3 (Population * population,        ///< Population.
                                Entity ** mother,
                                ///< Pointer to the mother entity.
                                Entity ** father,
                                ///< Pointer to the father entity.
                                gsl_rng * rng)
                                ///< GSL random numbers generator.
{
  unsigned int i, j, k;
  i = gsl_rng_uniform_int (rng, population->nsurvival);
  j = gsl_rng_uniform_int (rng, population->nsurvival);
  if (i < j)
    i = j;
  j = gsl_rng_uniform_int (rng, population->nsurvival);
  if (i < j)
    i = j;
  do
    j = gsl_rng_uniform_int (rng, population->nsurvival);
  while (i == j);
  k = gsl_rng_uniform_int (rng, population->nsurvival);
  if (i != k && k < j)
    j = k;
  k = gsl_rng_uniform_int (rng, population->nsurvival);
  if (i != k && k < j)
    j = k;
  *mother = population->entity + i;
  *father = population->entity + j;
}

/**
 * Function to select a pair of entities using the best and a random entity.
 */
static void
selection_reproduction_best (Population * population,   ///< Population.
                             Entity ** mother,
                             ///< Pointer to the mother entity.
                             Entity ** father,
                             ///< Pointer to the father entity.
                             gsl_rng * rng)
                             ///< GSL random numbers generator.
{
  *mother = population->entity
    + gsl_rng_uniform_int (rng, population->nsurvival);
  *father = population->entity;
}

/**
 * Function to select a pair of entities based on linear probability
 *   distribution with respect to rank.
 */
static void
selection_reproduction_linearrank (Population * population,     ///< Population.
                                   Entity ** mother,
                                   ///< Pointer to the mother entity.
                                   Entity ** father,
                                   ///< Pointer to the father entity.
                                   gsl_rng * rng)
                                   ///< GSL random numbers generator.
{
  unsigned int i, j;
  i = (unsigned int) ((1.0 - sqrt (gsl_rng_uniform (rng)))
                      * population->nsurvival);
  *mother = population->entity + i;
  do
    j = (unsigned int) ((1.0 - sqrt (gsl_rng_uniform (rng)))
                        * population->nsurvival);
  while (i == j);
  *father = population->entity + j;
}

/**
 * Function to select the selection operations.
 */
void
selection_init (unsigned int mutation_type,
                ///< Type of mutation selection operations.
                unsigned int reproduction_type,
                ///< Type of reproduction selection operations.
                unsigned int adaptation_type)
                ///< Type of adaptation selection operations.
{
  switch (mutation_type)
    {
    case SELECTION_MUTATION_TYPE_RANDOM:
      selection_mutation = &selection_mutation_random;
      break;
    case SELECTION_MUTATION_TYPE_BESTOF2:
      selection_mutation = &selection_mutation_bestof2;
      break;
    case SELECTION_MUTATION_TYPE_BESTOF3:
      selection_mutation = &selection_mutation_bestof3;
      break;
    case SELECTION_MUTATION_TYPE_BEST:
      selection_mutation = &selection_mutation_best;
      break;
    default:
      selection_mutation = &selection_mutation_linearrank;
    }
  switch (reproduction_type)
    {
    case SELECTION_REPRODUCTION_TYPE_RANDOM:
      selection_reproduction = &selection_reproduction_random;
      break;
    case SELECTION_REPRODUCTION_TYPE_BESTOF2:
      selection_reproduction = &selection_reproduction_bestof2;
      break;
    case SELECTION_REPRODUCTION_TYPE_BESTOF3:
      selection_reproduction = &selection_reproduction_bestof3;
      break;
    case SELECTION_REPRODUCTION_TYPE_BEST:
      selection_reproduction = &selection_reproduction_best;
      break;
    default:
      selection_reproduction = &selection_reproduction_linearrank;
    }
  switch (adaptation_type)
    {
    case SELECTION_ADAPTATION_TYPE_RANDOM:
      selection_adaptation = &selection_mutation_random;
      break;
    case SELECTION_ADAPTATION_TYPE_BESTOF2:
      selection_adaptation = &selection_mutation_bestof2;
      break;
    case SELECTION_ADAPTATION_TYPE_BESTOF3:
      selection_adaptation = &selection_mutation_bestof3;
      break;
    case SELECTION_ADAPTATION_TYPE_BEST:
      selection_adaptation = &selection_mutation_best;
      break;
    default:
      selection_adaptation = &selection_mutation_linearrank;
    }
}

max_stars_repo_name: jburguete/genetic

max_issues_repo_name: jburguete/genetic

max_forks_repo_name: jburguete/genetic


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

#ifndef __GSL_MATRIX_COMPLEX_LONG_DOUBLE_H__
#define __GSL_MATRIX_COMPLEX_LONG_DOUBLE_H__

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

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector_complex_long_double.h>
#include <gsl/gsl_blas_types.h>

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

__BEGIN_DECLS

typedef struct
{
  size_t size1;
  size_t size2;
  size_t tda;
  long double * data;
  gsl_block_complex_long_double * block;
  int owner;
} gsl_matrix_complex_long_double ;

typedef struct
{
  gsl_matrix_complex_long_double matrix;
} _gsl_matrix_complex_long_double_view;

typedef _gsl_matrix_complex_long_double_view gsl_matrix_complex_long_double_view;

typedef struct
{
  gsl_matrix_complex_long_double matrix;
} _gsl_matrix_complex_long_double_const_view;

typedef const _gsl_matrix_complex_long_double_const_view gsl_matrix_complex_long_double_const_view;

/* Allocation */

GSL_FUN gsl_matrix_complex_long_double *
gsl_matrix_complex_long_double_alloc (const size_t n1, const size_t n2);

GSL_FUN gsl_matrix_complex_long_double *
gsl_matrix_complex_long_double_calloc (const size_t n1, const size_t n2);

GSL_FUN gsl_matrix_complex_long_double *
gsl_matrix_complex_long_double_alloc_from_block (gsl_block_complex_long_double * b,
                                           const size_t offset,
                                           const size_t n1, const size_t n2, const size_t d2);

GSL_FUN gsl_matrix_complex_long_double *
gsl_matrix_complex_long_double_alloc_from_matrix (gsl_matrix_complex_long_double * b,
                                            const size_t k1, const size_t k2,
                                            const size_t n1, const size_t n2);

GSL_FUN gsl_vector_complex_long_double *
gsl_vector_complex_long_double_alloc_row_from_matrix (gsl_matrix_complex_long_double * m,
                                                const size_t i);

GSL_FUN gsl_vector_complex_long_double *
gsl_vector_complex_long_double_alloc_col_from_matrix (gsl_matrix_complex_long_double * m,
                                                const size_t j);

GSL_FUN void gsl_matrix_complex_long_double_free (gsl_matrix_complex_long_double * m);

/* Views */

GSL_FUN _gsl_matrix_complex_long_double_view
gsl_matrix_complex_long_double_submatrix (gsl_matrix_complex_long_double * m,
                            const size_t i, const size_t j,
                            const size_t n1, const size_t n2);

GSL_FUN _gsl_vector_complex_long_double_view
gsl_matrix_complex_long_double_row (gsl_matrix_complex_long_double * m, const size_t i);

GSL_FUN _gsl_vector_complex_long_double_view
gsl_matrix_complex_long_double_column (gsl_matrix_complex_long_double * m, const size_t j);

GSL_FUN _gsl_vector_complex_long_double_view
gsl_matrix_complex_long_double_diagonal (gsl_matrix_complex_long_double * m);

GSL_FUN _gsl_vector_complex_long_double_view
gsl_matrix_complex_long_double_subdiagonal (gsl_matrix_complex_long_double * m, const size_t k);

GSL_FUN _gsl_vector_complex_long_double_view
gsl_matrix_complex_long_double_superdiagonal (gsl_matrix_complex_long_double * m, const size_t k);

GSL_FUN _gsl_vector_complex_long_double_view
gsl_matrix_complex_long_double_subrow (gsl_matrix_complex_long_double * m,
                                 const size_t i, const size_t offset,
                                 const size_t n);

GSL_FUN _gsl_vector_complex_long_double_view
gsl_matrix_complex_long_double_subcolumn (gsl_matrix_complex_long_double * m,
                                    const size_t j, const size_t offset,
                                    const size_t n);

GSL_FUN _gsl_matrix_complex_long_double_view
gsl_matrix_complex_long_double_view_array (long double * base,
                             const size_t n1,
                             const size_t n2);

GSL_FUN _gsl_matrix_complex_long_double_view
gsl_matrix_complex_long_double_view_array_with_tda (long double * base,
                                      const size_t n1,
                                      const size_t n2,
                                      const size_t tda);

GSL_FUN _gsl_matrix_complex_long_double_view
gsl_matrix_complex_long_double_view_vector (gsl_vector_complex_long_double * v,
                              const size_t n1,
                              const size_t n2);

GSL_FUN _gsl_matrix_complex_long_double_view
gsl_matrix_complex_long_double_view_vector_with_tda (gsl_vector_complex_long_double * v,
                                       const size_t n1,
                                       const size_t n2,
                                       const size_t tda);

GSL_FUN _gsl_matrix_complex_long_double_const_view
gsl_matrix_complex_long_double_const_submatrix (const gsl_matrix_complex_long_double * m,
                                  const size_t i, const size_t j,
                                  const size_t n1, const size_t n2);

GSL_FUN _gsl_vector_complex_long_double_const_view
gsl_matrix_complex_long_double_const_row (const gsl_matrix_complex_long_double * m,
                            const size_t i);

GSL_FUN _gsl_vector_complex_long_double_const_view
gsl_matrix_complex_long_double_const_column (const gsl_matrix_complex_long_double * m,
                               const size_t j);

GSL_FUN _gsl_vector_complex_long_double_const_view
gsl_matrix_complex_long_double_const_diagonal (const gsl_matrix_complex_long_double * m);

GSL_FUN _gsl_vector_complex_long_double_const_view
gsl_matrix_complex_long_double_const_subdiagonal (const gsl_matrix_complex_long_double * m,
                                    const size_t k);

GSL_FUN _gsl_vector_complex_long_double_const_view
gsl_matrix_complex_long_double_const_superdiagonal (const gsl_matrix_complex_long_double * m,
                                      const size_t k);

GSL_FUN _gsl_vector_complex_long_double_const_view
gsl_matrix_complex_long_double_const_subrow (const gsl_matrix_complex_long_double * m,
                                       const size_t i, const size_t offset,
                                       const size_t n);

GSL_FUN _gsl_vector_complex_long_double_const_view
gsl_matrix_complex_long_double_const_subcolumn (const gsl_matrix_complex_long_double * m,
                                          const size_t j, const size_t offset,
                                          const size_t n);

GSL_FUN _gsl_matrix_complex_long_double_const_view
gsl_matrix_complex_long_double_const_view_array (const long double * base,
                                   const size_t n1,
                                   const size_t n2);

GSL_FUN _gsl_matrix_complex_long_double_const_view
gsl_matrix_complex_long_double_const_view_array_with_tda (const long double * base,
                                            const size_t n1,
                                            const size_t n2,
                                            const size_t tda);

GSL_FUN _gsl_matrix_complex_long_double_const_view
gsl_matrix_complex_long_double_const_view_vector (const gsl_vector_complex_long_double * v,
                                    const size_t n1,
                                    const size_t n2);

GSL_FUN _gsl_matrix_complex_long_double_const_view
gsl_matrix_complex_long_double_const_view_vector_with_tda (const gsl_vector_complex_long_double * v,
                                             const size_t n1,
                                             const size_t n2,
                                             const size_t tda);

/* Operations */

GSL_FUN void gsl_matrix_complex_long_double_set_zero (gsl_matrix_complex_long_double * m);
GSL_FUN void gsl_matrix_complex_long_double_set_identity (gsl_matrix_complex_long_double * m);
GSL_FUN void gsl_matrix_complex_long_double_set_all (gsl_matrix_complex_long_double * m, gsl_complex_long_double x);

GSL_FUN int gsl_matrix_complex_long_double_fread (FILE * stream, gsl_matrix_complex_long_double * m) ;
GSL_FUN int gsl_matrix_complex_long_double_fwrite (FILE * stream, const gsl_matrix_complex_long_double * m) ;
GSL_FUN int gsl_matrix_complex_long_double_fscanf (FILE * stream, gsl_matrix_complex_long_double * m);
GSL_FUN int gsl_matrix_complex_long_double_fprintf (FILE * stream, const gsl_matrix_complex_long_double * m, const char * format);

GSL_FUN int gsl_matrix_complex_long_double_memcpy(gsl_matrix_complex_long_double * dest, const gsl_matrix_complex_long_double * src);
GSL_FUN int gsl_matrix_complex_long_double_swap(gsl_matrix_complex_long_double * m1, gsl_matrix_complex_long_double * m2);
GSL_FUN int gsl_matrix_complex_long_double_tricpy(CBLAS_UPLO_t Uplo, CBLAS_DIAG_t Diag, gsl_matrix_complex_long_double * dest, const gsl_matrix_complex_long_double * src);

GSL_FUN int gsl_matrix_complex_long_double_swap_rows(gsl_matrix_complex_long_double * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_complex_long_double_swap_columns(gsl_matrix_complex_long_double * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_complex_long_double_swap_rowcol(gsl_matrix_complex_long_double * m, const size_t i, const size_t j);

GSL_FUN int gsl_matrix_complex_long_double_transpose (gsl_matrix_complex_long_double * m);
GSL_FUN int gsl_matrix_complex_long_double_transpose_memcpy (gsl_matrix_complex_long_double * dest, const gsl_matrix_complex_long_double * src);
GSL_FUN int gsl_matrix_complex_long_double_transpose_tricpy(CBLAS_UPLO_t Uplo_src, CBLAS_DIAG_t Diag, gsl_matrix_complex_long_double * dest, const gsl_matrix_complex_long_double * src);

GSL_FUN int gsl_matrix_complex_long_double_conjtrans_memcpy (gsl_matrix_complex_long_double * dest, const gsl_matrix_complex_long_double * src);

GSL_FUN int gsl_matrix_complex_long_double_equal (const gsl_matrix_complex_long_double * a, const gsl_matrix_complex_long_double * b);

GSL_FUN int gsl_matrix_complex_long_double_isnull (const gsl_matrix_complex_long_double * m);
GSL_FUN int gsl_matrix_complex_long_double_ispos (const gsl_matrix_complex_long_double * m);
GSL_FUN int gsl_matrix_complex_long_double_isneg (const gsl_matrix_complex_long_double * m);
GSL_FUN int gsl_matrix_complex_long_double_isnonneg (const gsl_matrix_complex_long_double * m);

GSL_FUN int gsl_matrix_complex_long_double_add (gsl_matrix_complex_long_double * a, const gsl_matrix_complex_long_double * b);
GSL_FUN int gsl_matrix_complex_long_double_sub (gsl_matrix_complex_long_double * a, const gsl_matrix_complex_long_double * b);
GSL_FUN int gsl_matrix_complex_long_double_mul_elements (gsl_matrix_complex_long_double * a, const gsl_matrix_complex_long_double * b);
GSL_FUN int gsl_matrix_complex_long_double_div_elements (gsl_matrix_complex_long_double * a, const gsl_matrix_complex_long_double * b);
GSL_FUN int gsl_matrix_complex_long_double_scale (gsl_matrix_complex_long_double * a, const gsl_complex_long_double x);
GSL_FUN int gsl_matrix_complex_long_double_scale_rows (gsl_matrix_complex_long_double * a, const gsl_vector_complex_long_double * x);
GSL_FUN int gsl_matrix_complex_long_double_scale_columns (gsl_matrix_complex_long_double * a, const gsl_vector_complex_long_double * x);
GSL_FUN int gsl_matrix_complex_long_double_add_constant (gsl_matrix_complex_long_double * a, const gsl_complex_long_double x);
GSL_FUN int gsl_matrix_complex_long_double_add_diagonal (gsl_matrix_complex_long_double * a, const gsl_complex_long_double x);

/***********************************************************************/
/* The functions below are obsolete                                    */
/***********************************************************************/
GSL_FUN int gsl_matrix_complex_long_double_get_row(gsl_vector_complex_long_double * v, const gsl_matrix_complex_long_double * m, const size_t i);
GSL_FUN int gsl_matrix_complex_long_double_get_col(gsl_vector_complex_long_double * v, const gsl_matrix_complex_long_double * m, const size_t j);
GSL_FUN int gsl_matrix_complex_long_double_set_row(gsl_matrix_complex_long_double * m, const size_t i, const gsl_vector_complex_long_double * v);
GSL_FUN int gsl_matrix_complex_long_double_set_col(gsl_matrix_complex_long_double * m, const size_t j, const gsl_vector_complex_long_double * v);
/***********************************************************************/

/* inline functions if you are using GCC */

GSL_FUN INLINE_DECL gsl_complex_long_double gsl_matrix_complex_long_double_get(const gsl_matrix_complex_long_double * m, const size_t i, const size_t j);
GSL_FUN INLINE_DECL void gsl_matrix_complex_long_double_set(gsl_matrix_complex_long_double * m, const size_t i, const size_t j, const gsl_complex_long_double x);

GSL_FUN INLINE_DECL gsl_complex_long_double * gsl_matrix_complex_long_double_ptr(gsl_matrix_complex_long_double * m, const size_t i, const size_t j);
GSL_FUN INLINE_DECL const gsl_complex_long_double * gsl_matrix_complex_long_double_const_ptr(const gsl_matrix_complex_long_double * m, const size_t i, const size_t j);

#ifdef HAVE_INLINE

INLINE_FUN
gsl_complex_long_double
gsl_matrix_complex_long_double_get(const gsl_matrix_complex_long_double * m,
                     const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      gsl_complex_long_double zero = {{0,0}};

      if (i >= m->size1)
        {
          GSL_ERROR_VAL("first index out of range", GSL_EINVAL, zero) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_VAL("second index out of range", GSL_EINVAL, zero) ;
        }
    }
#endif
  return *(gsl_complex_long_double *)(m->data + 2*(i * m->tda + j)) ;
}

INLINE_FUN
void
gsl_matrix_complex_long_double_set(gsl_matrix_complex_long_double * m,
                     const size_t i, const size_t j, const gsl_complex_long_double x)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_VOID("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_VOID("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  *(gsl_complex_long_double *)(m->data + 2*(i * m->tda + j)) = x ;
}

INLINE_FUN
gsl_complex_long_double *
gsl_matrix_complex_long_double_ptr(gsl_matrix_complex_long_double * m,
                             const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  return (gsl_complex_long_double *)(m->data + 2*(i * m->tda + j)) ;
}

INLINE_FUN
const gsl_complex_long_double *
gsl_matrix_complex_long_double_const_ptr(const gsl_matrix_complex_long_double * m,
                                   const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  return (const gsl_complex_long_double *)(m->data + 2*(i * m->tda + j)) ;
}

#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_MATRIX_COMPLEX_LONG_DOUBLE_H__ */

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
//  This random generator is a C++ wrapper for the GNU Scientific Library
//  Copyright (C) 2001 Torbjorn Vik

//  This program is free software; you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation; either version 2 of the License, or
//  (at your option) any later version.

//  This program is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.

//  You should have received a copy of the GNU General Public License
//  along with this program; if not, write to the Free Software
//  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
#ifndef __histogram_h
#define __histogram_h

#include <gsl/gsl_histogram.h>
#include <stdexcept>
#include <exception>

namespace gsl{
#ifndef __HP_aCC
using std::string;
using std::runtime_error;
#endif

//! Encapsulates the histogram object of gsl. Only uniformly spaced bins yet.
class histogram
{
public:
	histogram(int nBins, double xmin, double xmax)
	{
		h=gsl_histogram_calloc(nBins);
		if (!h)
		{
			throw runtime_error("Couldn't allocate memory for histogram");
		}
		gsl_histogram_set_ranges_uniform(h, xmin, xmax);
	}
	~histogram(){gsl_histogram_free(h);}

//@{ Updating and Accessing Methods
	int increment(double x){return gsl_histogram_increment(h, x);}
	int accumulate(double x, double weight){return gsl_histogram_accumulate(h, x, weight);}
	double get(int i) const {return gsl_histogram_get(h, i);}
	double& operator[](const uint & i)
	{
		const uint n = h->n;

		if (i >= n)
		{
			throw runtime_error("index lies outside valid range of 0 .. n - 1");
//			GSL_ERROR_VAL ("index lies outside valid range of 0 .. n - 1", GSL_EDOM, 0);
		}

		return h->bin[i];
	}
	const double& operator[](const uint & i) const //{return (*this)[i];/*gsl_histogram_get(h, i);*/}
	{
		const uint n = h->n;

		if (i >= n)
		{
			throw runtime_error("index lies outside valid range of 0 .. n - 1");
//			GSL_ERROR_VAL ("index lies outside valid range of 0 .. n - 1", GSL_EDOM, 0);
		}

		return h->bin[i];
	}

	void get_range(int i, double& xmin, double& xmax) const {gsl_histogram_get_range(h, i, &xmin, &xmax);}
//@}

//@{These functions return the maximum upper and minimum lower range limits
// and the number of bins of the histogram h. They provide a way of determining these values without
//    accessing the gsl_histogram struct directly.
	double max() const {return gsl_histogram_max(h);}
	double min() const {return gsl_histogram_min(h);}
	int bins()const {return gsl_histogram_bins(h);}
	int size()const {return gsl_histogram_bins(h);}

//@}

//@{ Histogram statistics
	double mean()const {return gsl_histogram_mean(h);} // not in gsl library ?
	double max_val() const {return gsl_histogram_max_val(h);}
	int max_bin() const {return gsl_histogram_max_bin(h);}
	double min_val() const {return gsl_histogram_min_val(h);}
	int min_bin() const {return gsl_histogram_min_bin(h);}
	double sum() const {return gsl_histogram_sum(h);}
//@}

//@{ Accessor for gsl compatibility
	gsl_histogram*       gslobj()       { return h;}
	const gsl_histogram* gslobj() const { return h;}
//@}

protected:
	gsl_histogram * h;
};
}

#endif // __histogram_h

max_stars_repo_name: entn-at/GlottDNN

max_issues_repo_name: entn-at/GlottDNN

max_forks_repo_name: entn-at/GlottDNN


text:
/* -*- Mode: C++; tab-width: 4; c-basic-offset: 4; indent-tabs-mode: nil -*- */
/*
 *     Copyright 2018 Couchbase, Inc.
 *
 *   Licensed under the Apache License, Version 2.0 (the "License");
 *   you may not use this file except in compliance with the License.
 *   You may obtain a copy of the License at
 *
 *       http://www.apache.org/licenses/LICENSE-2.0
 *
 *   Unless required by applicable law or agreed to in writing, software
 *   distributed under the License is distributed on an "AS IS" BASIS,
 *   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *   See the License for the specific language governing permissions and
 *   limitations under the License.
 */

#pragma once

#include "engine_error.h"
#include "types.h"

#include <gsl/gsl>

struct ServerDocumentIface {
    virtual ~ServerDocumentIface() = default;

    /**
     * This callback is called from the underlying engine right before
     * it is linked into the list of available documents (it is currently
     * not visible to anyone). The engine should have validated all
     * properties set in the document by the client and the core, and
     * assigned a new CAS number for the document (and sequence number if
     * the underlying engine use those).
     *
     * The callback may at this time do post processing of the document
     * content (it is allowed to modify the content data, but not
     * reallocate or change the size of the data in any way).
     *
     * Given that the engine MAY HOLD LOCKS when calling this function
     * the core is *NOT* allowed to acquire *ANY* locks (except for doing
     * some sort of memory allocation for a temporary buffer).
     *
     * @param cookie The cookie provided to the engine for the storage
     *               command which may (which may hold more context)
     * @param info the items underlying data
     * @return ENGINE_SUCCESS means that the underlying engine should
     *                        proceed to link the item. All other
     *                        error codes means that the engine should
     *                        *NOT* link the item
     */
    virtual ENGINE_ERROR_CODE pre_link(gsl::not_null<const void*> cookie,
                                       item_info& info) = 0;

    /**
     * This callback is called from the underlying engine right before
     * a particular document expires. The callback is responsible for
     * modifying the contents of the itm_info passed in. The updated
     * total size is available in itm_info.nbytes.
     *
     * @param itm_info info pertaining to the item that is to be expired.
     * @return true indicating that the info has been modified in itm_info.
     *         false indicating that there is no data available in itm_info.
     * @throws std::bad_alloc in case of memory allocation failure
     * @throws std::logic_error if the data has grown
     */
    virtual bool pre_expiry(item_info& itm_info) = 0;
};


text:
#include <gsl/gsl_blas.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>

// author:    Andrew Liew <liew@arch.ethz.ch>
// copyright: Copyright 2018, BLOCK Research Group - ETH Zurich
// license:   MIT License
// email:     liew@arch.ethz.ch

int i;
int j;

// ==========================================================================================================
// One vector
// ==========================================================================================================

double length_vector(gsl_vector *u)
{
    double L = gsl_blas_dnrm2(u);

    return L;
}

double length_vector_squared(gsl_vector *u)
{
    double a = 0.;

    for (i = 0; i < 3; i++)
    {
        a += gsl_pow_2(gsl_vector_get(u, i));
    }

    return a;
}

void normalize_vector(gsl_vector *u)
{
    gsl_vector_scale(u, 1./length_vector(u));
}

void scale_vector(gsl_vector *u, double a)
{
    gsl_vector_scale(u, a);
}

void vector_from_pointer(double *ptr, gsl_vector *w)
{
    gsl_vector_set(w, 0, *(ptr + 0));
    gsl_vector_set(w, 1, *(ptr + 1));
    gsl_vector_set(w, 2, *(ptr + 2));
}

// ==========================================================================================================
// Two vectors
// ==========================================================================================================

void add_vectors(const gsl_vector *u, const gsl_vector *v, gsl_vector *w)
{
    gsl_vector_memcpy(w, u);
    gsl_vector_add(w, v);
}

void cross_vectors(const gsl_vector *u, const gsl_vector *v, gsl_vector *w)
{
    double w1 = gsl_vector_get(u, 1) * gsl_vector_get(v, 2) - gsl_vector_get(u, 2) * gsl_vector_get(v, 1);
    double w2 = gsl_vector_get(u, 2) * gsl_vector_get(v, 0) - gsl_vector_get(u, 0) * gsl_vector_get(v, 2);
    double w3 = gsl_vector_get(u, 0) * gsl_vector_get(v, 1) - gsl_vector_get(u, 1) * gsl_vector_get(v, 0);

    gsl_vector_set(w, 0, w1);
    gsl_vector_set(w, 1, w2);
    gsl_vector_set(w, 2, w3);
}

double dot_vectors(const gsl_vector *u, const gsl_vector *v)
{
    double a = 0.;

    for (i = 0; i < 3; i++)
    {
        a += gsl_vector_get(u, i) * gsl_vector_get(v, i);
    }

    return a;
}

void subtract_vectors(const gsl_vector *u, const gsl_vector *v, gsl_vector *w)
{
    gsl_vector_memcpy(w, u);
    gsl_vector_sub(w, v);
}

max_stars_repo_name: philianeles/compas

max_issues_repo_name: philianeles/compas

max_forks_repo_name: philianeles/compas


text:
#pragma once

#include "halleystring.h"
#include "halley/data_structures/vector.h"
#include <optional>
#include <gsl/span>

namespace Halley {
    class FuzzyTextMatcher {
    public:
		class Score {
		public:
			int curJumps = 0;
			int jumpLength = 0;
			int sections = 0;
			int sectionLens = 0;
			int sectionPos = 0xFFFF;
			Vector<std::pair<uint16_t, uint16_t>> matchPositions;

			bool operator<(const Score& other) const;

			Score advance(int jumpLen, int sectionPos, int newSectionLen) const;
			void makeMatchPositions(const Vector<int>& breadcrumbs);
		};

        class Result {
        public:
        	Result() = default;
            Result(String str, String id, Score score);

            bool operator<(const Result& other) const;

        	const String& getString() const;
        	const String& getId() const;
        	gsl::span<const std::pair<uint16_t, uint16_t>> getMatchPositions() const;

        private:
        	String str;
        	String id;
        	Vector<std::pair<uint16_t, uint16_t>> matchPositions;
        	Score score;
        };

        FuzzyTextMatcher(bool caseSensitive, std::optional<size_t> resultsLimit);

    	void addStrings(Vector<String> strings);
    	void addString(String string, String id = "");
    	void clear();

    	Vector<Result> match(const String& query) const;

    private:
        struct Entry {
	        String string;
            String id;
        };

    	Vector<Entry> strings;
    	bool caseSensitive;
    	std::optional<size_t> resultsLimit;

    	std::optional<Result> match(const String& str, const String& id, const StringUTF32& query) const;
    };
}

max_stars_repo_name: amrezzd/halley

max_issues_repo_name: amrezzd/halley

max_forks_repo_name: amrezzd/halley


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

/* Jan/2001 Modified by Brian Gough. Minor changes for GSL */

#ifndef __GSL_NTUPLE_H__
#define __GSL_NTUPLE_H__

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

#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_histogram.h>

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

__BEGIN_DECLS

typedef struct {
    FILE * file;
    void * ntuple_data;
    size_t size;
} gsl_ntuple;

typedef struct {
  int (* function) (void * ntuple_data, void * params);
  void * params;
} gsl_ntuple_select_fn;

typedef struct {
  double (* function) (void * ntuple_data, void * params);
  void * params;
} gsl_ntuple_value_fn;

GSL_FUN gsl_ntuple *
gsl_ntuple_open (char * filename, void * ntuple_data, size_t size);

GSL_FUN gsl_ntuple *
gsl_ntuple_create (char * filename, void * ntuple_data, size_t size);

GSL_FUN int gsl_ntuple_write (gsl_ntuple * ntuple);
GSL_FUN int gsl_ntuple_read (gsl_ntuple * ntuple);

GSL_FUN int gsl_ntuple_bookdata (gsl_ntuple * ntuple);  /* synonym for write */

GSL_FUN int gsl_ntuple_project (gsl_histogram * h, gsl_ntuple * ntuple,
                        gsl_ntuple_value_fn *value_func,
                        gsl_ntuple_select_fn *select_func);

GSL_FUN int gsl_ntuple_close (gsl_ntuple * ntuple);

__END_DECLS

#endif /* __GSL_NTUPLE_H__ */

max_stars_repo_name: berkus/music-cs

max_issues_repo_name: berkus/music-cs

max_forks_repo_name: berkus/music-cs


text:
#ifndef __NODE_H__
#define __NODE_H__

//#include <pybind11/functional.h>
#include <pybind11/numpy.h>
//#include <pybind11/stl.h>
//#include <pybind11/pybind11.h>
#include <torch/script.h>
#include <pybind11/embed.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <iomanip>
#include <sstream>
#include <cmath>
#include <math.h>
#include <string>
#include <memory>
#include <mutex>
#include "state.h"

namespace py = pybind11;

struct Node{
public:
	int select(State *state, double cpuct);
	struct Node* next_node(int action);
	struct Node* make_move(int action);
	void backpropagate(int action, float value);
	void set_prior(torch::Tensor p, gsl_rng *r, const double *alpha, double dir_eps);
	void set_prior(State *state,    gsl_rng *r, const double *alpha, const double dir_eps);
	std::array<long int, SIZE>* counts();
	std::string repr();
	std::string print_u(State *state, double cpuct);
	bool is_null(int a);
	long long int nodeN;
	std::mutex node_mutex;
	std::array<double, SIZE> childP;
//	std::string name;

	//Node(std::string &name, int a);
	Node();
	~Node();
private:
	std::array<std::unique_ptr<Node>, SIZE> child;
	int child_cnt;
	std::array<long int, SIZE> childN;
	std::array<double, SIZE> childW;
};

#endif

max_stars_repo_name: SlavkaMichal/alphazero

max_issues_repo_name: SlavkaMichal/alphazero

max_forks_repo_name: SlavkaMichal/alphazero