text:
/*
** read nifti-1 or nifti-2 to vista format
**
** G.Lohmann, MPI-KYB, 2016
*/
#include <viaio/Vlib.h>
#include <viaio/VImage.h>
#include <viaio/mu.h>
#include <viaio/option.h>

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <strings.h>
#include <math.h>
#include <ctype.h>

#include <gsl/gsl_math.h>

#include <nifti/nifti2.h>
#include <nifti/nifti1_io.h>

#define MIN_HEADER_SIZE 348
#define NII_HEADER_SIZE 352

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

/*extern int isinff(float x);
  extern int isnanf(float x);
*/
extern char *VReadGzippedData(char *filename,size_t *len);
extern char *VReadUnzippedData(char *filename,VBoolean nofail,size_t *size);
extern char *VReadDataContainer(char *filename,VBoolean nofail,size_t *size);
extern FILE *VReadInputFile (char *filename,VBoolean nofail);
extern int  CheckGzip(char *filename);

extern void VByteSwapNiftiHeader(nifti_1_header *hdr);
extern void VByteSwapNifti2Header(nifti_2_header *hdr);
extern void VByteSwapData(char *data,size_t ndata,size_t nsize);
extern float strtof(const char *str, char **endptr);

int NiftiVersion(char *databuffer,int *swap)
{
  int *k = (int *)(&databuffer[0]);
  int type = *k;
  *swap = 0;
  if (type == 348) return 1;
  else if (type == 540) return 2;
  else {
    nifti_swap_4bytes(1,k);
    type = *k;
    *swap = 1;
    if (type == 348) return 1;
    else if (type == 540) return 2;
    else return 0;
  }
  return 0;
}

/* read repetition time from description field instead of pixdim  */
float TRfromString(const char *str)
{
  float f_tr   = -1;    /* return value on failure */
  float factor = 1000;  /* expect seconds as default */
  char* c_tr;

  if( ( c_tr = strstr( str, "TR=" ) ) || ( c_tr = strstr( str, "TR =" ) ) ) {
      while( *c_tr++ != '=' )
    ;

      while( *c_tr && isspace( *c_tr ) )
	c_tr++;

      f_tr = strtof( c_tr, &c_tr );

      if( f_tr > 0.0 ) {
	while( *c_tr && isspace( *c_tr ) )
	  c_tr++;

	if( strncasecmp( c_tr, "ms", 2 ) == 0 )
	  factor = 1;
      }
  }
  return f_tr * factor;
}

/* set dimension in geo info */
void VUpdateGeoinfo(VAttrList geolist,int dimtype,VLong tr)
{
  double *D = VGetGeoDim(geolist,NULL);
  D[0] = (double)dimtype;
  VSetGeoDim(geolist,D);

  double *E = VGetGeoPixdim(geolist,NULL);
  if (tr > 0) {
    E[4] = (double)tr;
    VSetGeoPixdim(geolist,E);
  }
}

/*
 *  VIniImage
 *
 *  Allocate Image structure and fill voxel values with databuffer.
 *  Returns a pointer to the image if successful, zero otherwise.
 *  Useful for converting from 4D nifti to list of 3D vista images
 */
VImage VIniImage (int nbands, int nrows, int ncolumns, VRepnKind pixel_repn, char *databuffer)
{
  size_t row_size = ncolumns * VRepnSize (pixel_repn);
  size_t row_index_size = nbands * nrows * sizeof (char *);
  size_t band_index_size = nbands * sizeof (char **);
  size_t pixel_size=0;
  VImage image;
  int band=0, row=0;
  char *p = NULL;

  /* Check parameters: */
  if (nbands < 1) VError("VIniImage: Invalid number of bands: %d", (int) nbands);
  if (nrows < 1)  VError ("VIniImage: Invalid number of rows: %d", (int) nrows);
  if (ncolumns < 1) VError ("VIniImage: Invalid number of columns: %d",(int) ncolumns);

#define AlignUp(v, b) ((((v) + (b) - 1) / (b)) * (b))

  /* Initialize VImage data struct */
  pixel_size = VRepnSize (pixel_repn);
  p = VMalloc(AlignUp (sizeof (VImageRec) + row_index_size + band_index_size, pixel_size));

  image = (VImage) p;
  image->nbands = nbands;
  image->nrows  = nrows;
  image->ncolumns = ncolumns;
  image->flags    = VImageSingleAlloc;
  image->pixel_repn = pixel_repn;
  image->attributes = VCreateAttrList ();
  image->band_index = (VPointer **) (p += sizeof (VImageRec));
  image->row_index  = (VPointer *) (p += band_index_size);
  image->data = (VPointer) AlignUp((long) (databuffer+4), pixel_size);

  image->nframes = nbands;
  image->nviewpoints = image->ncolors = image->ncomponents = 1;

  /* Initialize the indices: */
  for (band = 0; band < nbands; band++)
    image->band_index[band] = image->row_index + band * nrows;
  for (row = 0, p = image->data; row < nbands * nrows; row++, p += row_size)
    image->row_index[row] = p;

  return image;

#undef AlignUp
}

#define xgetval(data,index,Datatype) \
{ \
  Datatype *k = (Datatype *)(&data[index]); \
  u = (float)(*k); \
}

float VGetValue(char *data,size_t index,int datatype)
{
  float u=0;
  switch(datatype) {
  case DT_BINARY:
    xgetval(data,index,VBit);
    break;
  case DT_UNSIGNED_CHAR:
    xgetval(data,index,VUByte);
    break;
  case DT_SIGNED_SHORT:
    xgetval(data,index,short);
    break;
  case DT_SIGNED_INT:
    xgetval(data,index,int);
    break;
  case DT_FLOAT:
    xgetval(data,index,float);
    break;
  case DT_DOUBLE:
    xgetval(data,index,double);
    break;
  case DT_INT8:
    xgetval(data,index,VSByte);
    break;
  case DT_UINT16:
    xgetval(data,index,unsigned short);
    break;
  case DT_UINT32:
    xgetval(data,index,unsigned int);
    break;
  case DT_INT64:
    xgetval(data,index,long);
    break;
  case DT_UINT64:
    xgetval(data,index,unsigned long);
    break;
  default:
    VError(" unknown datatype %d",datatype);
  }
  if (isinf(u) || isnan(u))
    u = 0;
  return u;
}

void VCleanData(VImage src)
{
  int b,r,c;
  double u=0;
  for (b=0; b<VImageNBands(src); b++) {
    for (r=0; r<VImageNRows(src); r++) {
      for (c=0; c<VImageNColumns(src); c++) {
	u = VGetPixel(src,b,r,c);
	if (isinf(u) || isnan(u)) {
	  u = 0;
	  VSetPixel(src,b,r,c,u);
	}
      }
    }
  }
}

/* get image statistics for re-scaling parameters */
void VDataStats(char *data,size_t ndata,size_t nsize,int datatype,float *xmin,float *xmax)
{
  size_t i,n;
  float u=0;
  float zmin = VRepnMaxValue(VFloatRepn);
  float zmax = VRepnMinValue(VFloatRepn);

  n=0;
  for (i=0; i<ndata; i+= nsize) {
    u = VGetValue(data,i,datatype);
    if (fabs(u) < TINY) continue;
    if (u < zmin) zmin = u;
    if (u > zmax) zmax = u;
    n++;
  }
  if (n < 1) VError(" no non-zero data points found");

  *xmin = (zmin+TINY);
  *xmax = (zmax-TINY);
}

/* list of 3D images */
void Nii2Vista3DList(char *data,size_t nsize,size_t nslices,size_t nrows,size_t ncols,size_t nt,
		     VRepnKind dst_repn,int datatype,float scl_slope,float scl_inter,
		     VString voxelstr,VLong tr,VAttrList out_list)
{
  int slice,row,col;
  float u=0;
  size_t i;
  size_t add=0;
  size_t nrnc = nrows*ncols;
  size_t npix = nrnc*nslices;
  size_t ndata = nt * npix * nsize;

  if (nsize == 1) add=4;
  if (nsize == 2) add=2;
  if (nsize == 4) add=1;

  VImage *dst = (VImage *) VCalloc(nt,sizeof(VImage));
  if (nt > 1) VAppendAttr(out_list,"nimages",NULL,VLongRepn,(VLong) nt);

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

    /* rescale if needed */
    if (fabs(scl_slope) > 0 || fabs(scl_inter) > 0) {

      VImage tmp = VCreateImage(nslices,nrows,ncols,dst_repn);
      VSetAttr(VImageAttrList(tmp),"voxel",NULL,VStringRepn,voxelstr);
      if (tr > 0) VSetAttr(VImageAttrList(tmp),"repetition_time",NULL,VLongRepn,(VLong)tr);
      for (slice=0; slice<nslices; slice++) {
	for (row=0; row<nrows; row++) {
	  for (col=0; col<ncols; col++) {
	    const size_t src_index = (col + row*ncols + slice*nrnc + i*npix)*nsize;
	    if (src_index >= ndata || src_index < 0) continue;
	    u = VGetValue(data,src_index,datatype);
	    u = scl_slope*u + scl_inter;
	    VSetPixel(tmp,slice,row,col,(double)u);
	  }
	}
      }
      VAppendAttr(out_list,"image",NULL,VImageRepn,tmp);
    }

    /* otherwise just copy */
    else {
      const size_t index = (i*npix-add)*nsize;
      dst[i] = VIniImage(nslices,nrows,ncols,dst_repn,&data[index]);
      VSetAttr(VImageAttrList(dst[i]),"voxel",NULL,VStringRepn,voxelstr);
      if (tr > 0) VSetAttr(VImageAttrList(dst[i]),"repetition_time",NULL,VLongRepn,(VLong)tr);
      VAppendAttr(out_list,"image",NULL,VImageRepn,dst[i]);
      VCleanData(dst[i]);
    }
  }
}

/* 4D time series data */
void Nii2Vista4D(char *data,size_t nsize,size_t nslices,size_t nrows,size_t ncols,size_t nt,
		 VRepnKind dst_repn,int datatype,VBoolean do_scaling,float scl_slope,float scl_inter,
		 VString voxelstr,double *slicetime,VLong tr,VAttrList out_list)
{
  size_t slice,row,col,ti;
  size_t nrnc = nrows*ncols;
  size_t npix = nrnc*nslices;
  size_t ndata = nt * npix * nsize;
  size_t add=0;
  float u=0;

  /* rescale to 16bit integer if needed */
  float xmin=0,xmax=0,umin=0,umax=0;
  if (do_scaling && dst_repn != VShortRepn) {
    dst_repn = VShortRepn;
    VDataStats(data,ndata,nsize,datatype,&xmin,&xmax);
    fprintf(stderr," data range: [%f, %f]\n",xmin,xmax);
    umin = 0;
    umax = VRepnMaxValue(VShortRepn);
  }

  VImage *dst = (VImage *) VCalloc(nslices,sizeof(VImage));

  for (slice=0; slice<nslices; slice++) {
    dst[slice] = VCreateImage(nt,nrows,ncols,dst_repn);
    if (!dst[slice]) VError(" err allocating image");
    VFillImage(dst[slice],VAllBands,0);

    VSetAttr(VImageAttrList(dst[slice]),"voxel",NULL,VStringRepn,voxelstr);
    if (tr > 0) VSetAttr(VImageAttrList(dst[slice]),"repetition_time",NULL,VLongRepn,(VLong)tr);
    if (slicetime != NULL)
      VSetAttr(VImageAttrList(dst[slice]),"slice_time",NULL,VShortRepn,(VShort)slicetime[slice]);

    for (ti=0; ti<nt; ti++) {
      for (row=0; row<nrows; row++) {
	for (col=0; col<ncols; col++) {
	  const size_t src_index = (col + row*ncols + slice*nrnc + ti*npix + add)*nsize;
	  if (src_index >= ndata || src_index < 0) continue;

	  if (!do_scaling) {
	    u = VGetValue(data,src_index,datatype);
	    if (fabs(scl_slope) > 0 || fabs(scl_inter) > 0) {
	      u = scl_slope*u + scl_inter;
	    }
	    if (VPixelRepn(dst[slice]) == VShortRepn) {
	      VPixel(dst[slice],ti,row,col,VShort) = u;
	    }
	    else if (VPixelRepn(dst[slice]) == VFloatRepn) {
	      VPixel(dst[slice],ti,row,col,VFloat) = u;
	    }
	    else {
	      VSetPixel(dst[slice],ti,row,col,(double)u);
	    }
	  }
	  else {
	    u = VGetValue(data,src_index,datatype);
	    u = umax * (u-xmin)/(xmax-xmin);
	    if (u < umin) u = umin;
	    if (u > umax) u = umax;
	    VPixel(dst[slice],ti,row,col,VShort) = (VShort)u;
	  }
	}
      }
    }
    VCleanData(dst[slice]);
    VAppendAttr(out_list,"image",NULL,VImageRepn,dst[slice]);
  }
}

/* copy nifti header infos to geolist in vista header */
double *VGetNiftiHeader(VAttrList geolist,nifti_2_header hdr,VLong tr)
{

  /* units in lipsia: mm and sec */
  char xyzt = hdr.xyzt_units;
  int spaceunits = XYZT_TO_SPACE(xyzt);
  int timeunits  = XYZT_TO_TIME(xyzt);
  double xscale  = 1.0;
  double tscale  = 1.0;

  if (spaceunits == NIFTI_UNITS_MICRON) xscale=1000.0;
  if (timeunits == NIFTI_UNITS_SEC) tscale=1000.0;

  /* dim info */
  VSetAttr(geolist,"dim_info",NULL,VShortRepn,(VShort)hdr.dim_info);

  /* dim */
  /* hdr.dim[0]==3 means dim=3,  hdr.dim[0]==4 means dim=4 (timeseries) */
  float *E = VCalloc(8,sizeof(float));
  VAttrList elist = VCreateAttrList();
  VBundle ebundle = VCreateBundle ("bundle",elist,8*sizeof(float),(VPointer)E);
  VSetAttr(geolist,"dim",NULL,VBundleRepn,ebundle);
  int i;
  for (i=0; i<8; i++) E[i] = hdr.dim[i];
  for (i=5; i<8; i++) E[i] = 0;

  /* pixdim */
  float *D = VCalloc(8,sizeof(float));
  for (i=0; i<8; i++) D[i] = hdr.pixdim[i];
  for (i=1; i<=3; i++) D[i] *= xscale;
  D[4] *= tscale;
  if (tr > 0) D[4] = (double)tr;

  /* override if TR info is in description field */
  float xtr = TRfromString(hdr.descrip);
  if (xtr > 0) {
    D[4] = (double)(xtr*tscale);
    fprintf(stderr," reading TR from description field, TR= %.4f ms\n",D[4]);
  }
  if (fabs(D[0]) < 0.0001) D[0] = 1.0;  /* if not specified, assume pixdim[0] = 1 */

  for (i=5; i<8; i++) D[i] = 0;
  VAttrList dlist = VCreateAttrList();
  VBundle dbundle = VCreateBundle ("bundle",dlist,8*sizeof(float),(VPointer)D);
  VSetAttr(geolist,"pixdim",NULL,VBundleRepn,dbundle);

  /* sform */
  VImage sform = VCreateImage(1,4,4,VFloatRepn);
  VFillImage (sform,VAllBands,0);
  int j;
  for (j=0; j<4; j++) {
    VPixel(sform,0,0,j,VFloat) = hdr.srow_x[j];
    VPixel(sform,0,1,j,VFloat) = hdr.srow_y[j];
    VPixel(sform,0,2,j,VFloat) = hdr.srow_z[j];
  }
  VSetAttr(geolist,"sform_code",NULL,VShortRepn,(VShort)hdr.sform_code);
  VSetAttr(geolist,"sform",NULL,VImageRepn,sform);

  /* qform */
  size_t dim=6;
  float *Q = VCalloc(dim,sizeof(float));
  VAttrList qlist = VCreateAttrList();
  VBundle qbundle = VCreateBundle ("bundle",qlist,dim*sizeof(float),(VPointer)Q);
  Q[0] = hdr.quatern_b;
  Q[1] = hdr.quatern_c;
  Q[2] = hdr.quatern_d;
  Q[3] = hdr.qoffset_x;
  Q[4] = hdr.qoffset_y;
  Q[5] = hdr.qoffset_z;

  VShort qform_code = hdr.qform_code;
  if (hdr.sform_code==0 && qform_code==0) qform_code = 1;    /* if both codes are unspecified, assume scanner coord */
  VSetAttr(geolist,"qform_code",NULL,VShortRepn,(VShort)qform_code);
  VSetAttr(geolist,"qform",NULL,VBundleRepn,qbundle);

  /* MRI encoding directions */
  int freq_dim=0,phase_dim=0,slice_dim=0;
  if (hdr.dim_info != 0) {
    freq_dim  = DIM_INFO_TO_FREQ_DIM (hdr.dim_info );
    phase_dim = DIM_INFO_TO_PHASE_DIM(hdr.dim_info );
    slice_dim = DIM_INFO_TO_SLICE_DIM(hdr.dim_info );
    VSetAttr(geolist,"freq_dim",NULL,VShortRepn,(VShort)freq_dim);
    VSetAttr(geolist,"phase_dim",NULL,VShortRepn,(VShort)phase_dim);
    VSetAttr(geolist,"slice_dim",NULL,VShortRepn,(VShort)slice_dim);
  }
  if (slice_dim == 0) return NULL;

  /* slicetiming information */
  int slice_start = hdr.slice_start;
  int slice_end = hdr.slice_end;
  int slice_code = hdr.slice_code;
  double slice_duration = hdr.slice_duration;
  if (slice_code == 0 || slice_duration < 1.0e-10) return NULL;
  fprintf(stderr," slice duration: %.2f ms\n", slice_duration);
  VSetAttr(geolist,"slice_start",NULL,VShortRepn,(VShort)slice_start);
  VSetAttr(geolist,"slice_end",NULL,VShortRepn,(VShort)slice_end);
  VSetAttr(geolist,"slice_code",NULL,VShortRepn,(VShort)slice_code);
  VSetAttr(geolist,"slice_duration",NULL,VFloatRepn,(VFloat)slice_duration);

  int nslices = (int)E[3];
  double *slicetimes = (double *) VCalloc(nslices,sizeof(double));
  VGetSlicetimes(slice_start,slice_end,slice_code,slice_duration,slicetimes);
  return slicetimes;
}

VAttrList Nifti1_to_Vista(char *databuffer,VLong tr,VBoolean attrtype,VBoolean do_scaling,VBoolean *ok)
{
  /* read header */
  nifti_1_header hdr1;
  nifti_2_header hdr;
  int i=0,swap=0;
  size_t header_size=0;

  int nifti_version = NiftiVersion(databuffer,&swap);
  if (nifti_version < 1) VError(" unknown nifti version");
  /* fprintf(stderr," nifti version: %d\n",nifti_version); */

  if (nifti_version == 1) {  /* nifti-1 */
    header_size = 348;
    memcpy(&hdr1,databuffer,header_size);
    swap = NIFTI_NEEDS_SWAP(hdr1);
    if (swap == 1) VByteSwapNiftiHeader(&hdr1);

    hdr.datatype = hdr1.datatype;
    for (i=0; i<8; i++) hdr.dim[i]=hdr1.dim[i];
    for (i=0; i<8; i++) hdr.pixdim[i]=hdr1.pixdim[i];
    hdr.bitpix=hdr1.bitpix;
    hdr.vox_offset=hdr1.vox_offset;
    hdr.intent_p1 = hdr1.intent_p1;
    hdr.intent_p2 = hdr1.intent_p2;
    hdr.intent_p3 = hdr1.intent_p3;
    hdr.scl_slope = hdr1.scl_slope;
    hdr.scl_inter = hdr1.scl_inter;
    hdr.slice_duration = hdr1.slice_duration;
    hdr.slice_start = hdr1.slice_start;
    hdr.slice_end = hdr1.slice_end;
    for (i=0; i<80; i++) hdr.descrip[i] = hdr1.descrip[i];
    hdr.qform_code = hdr1.qform_code;
    hdr.sform_code = hdr1.sform_code;
    hdr.quatern_b = hdr1.quatern_b;
    hdr.quatern_c = hdr1.quatern_c;
    hdr.quatern_d = hdr1.quatern_d;
    hdr.qoffset_x = hdr1.qoffset_x;
    hdr.qoffset_y = hdr1.qoffset_y;
    hdr.qoffset_z = hdr1.qoffset_z;
    for (i=0; i<4; i++) hdr.srow_x[i]=hdr1.srow_x[i];
    for (i=0; i<4; i++) hdr.srow_y[i]=hdr1.srow_y[i];
    for (i=0; i<4; i++) hdr.srow_z[i]=hdr1.srow_z[i];
    hdr.xyzt_units = hdr1.xyzt_units;
    hdr.dim_info = hdr1.dim_info;
  }
  else if (nifti_version == 2) {  /* nifti-2 */
    header_size = 540;
    memcpy(&hdr,databuffer,header_size);
    if (swap) {
      VByteSwapNifti2Header(&hdr);
    }
  }
  else {
    VError(" not a nifti file");
  }

  /* get data type */
  VRepnKind dst_repn = DT_UNKNOWN;
  int datatype = (int)hdr.datatype;

  switch(datatype) {
  case DT_UNKNOWN:
    VError(" unknown data type");
    break;
  case DT_BINARY:
    dst_repn = VBitRepn;
    break;
  case DT_UNSIGNED_CHAR:
    dst_repn = VUByteRepn;
    break;
  case DT_SIGNED_SHORT:
    dst_repn = VShortRepn;
    break;
  case DT_SIGNED_INT:
    dst_repn = VIntegerRepn;
    break;
  case DT_FLOAT:
    dst_repn = VFloatRepn;
    break;
  case DT_DOUBLE:
    dst_repn = VDoubleRepn;
    break;
  case DT_INT8:
    dst_repn = VSByteRepn;
    break;
  case DT_UINT16:
    dst_repn = VUShortRepn;
    break;
  case DT_UINT32:
    dst_repn = VUIntegerRepn;
    break;
  case DT_INT64:
    dst_repn = VLongRepn;
    break;
  case DT_UINT64:
    dst_repn = VULongRepn;
    break;
  default:
    VError(" unknown data type %d",datatype);
  }

  /* get scaling parameters */
  double scl_slope = (double)hdr.scl_slope;
  double scl_inter = (double)hdr.scl_inter;
  double tiny=1.0e-6;
  if (fabs(scl_slope) < tiny && fabs(scl_inter) < tiny) {
    scl_slope = 1.0;
    scl_inter = 0;
  }
  else if (fabs(scl_slope-1.0) < tiny && fabs(scl_inter) < tiny) {
    scl_slope = 1.0;
    scl_inter = 0;
  }
  else {
    dst_repn = VFloatRepn;
    do_scaling = FALSE;  /* no automatic scaling */
  }

  /* number of values stored at each time point */
  if (hdr.dim[5] > 1) VError("data type not supported, dim[5]= %d\n",hdr.dim[5]);

  /* image size */
  /* size_t dimtype = (size_t)hdr.dim[0]; */
  size_t ncols   = (size_t)hdr.dim[1];
  size_t nrows   = (size_t)hdr.dim[2];
  size_t nslices = (size_t)hdr.dim[3];
  size_t nt      = (size_t)hdr.dim[4];

  /* fill data container */
  size_t bytesize = 8;
  if (dst_repn == VBitRepn) bytesize = 1;
  size_t nsize   = hdr.bitpix/bytesize;
  size_t npixels = nslices * nrows * ncols;
  size_t ndata   = nt * npixels * nsize;
  size_t vox_offset = (size_t)hdr.vox_offset;
  size_t startdata = MIN_HEADER_SIZE;
  if (vox_offset > 0) startdata = vox_offset;
  char *data = &databuffer[startdata];

  /* byte swap image data, if needed */
  if (swap == 1) {
    VByteSwapData(data,ndata,nsize);
  }

  /* repetition time (may be wrong in some cases) */
  char xyzt = hdr.xyzt_units;
  int tcode = XYZT_TO_TIME(xyzt);
  float factor = 1.0;
  float xtr=0;
  if (tcode == NIFTI_UNITS_MSEC) factor = 1.0;
  if (tcode == NIFTI_UNITS_SEC) factor = 1000.0;

  if (nt > 1)  {
    if (tr == 0) tr = (short)(factor*hdr.pixdim[4]);

    /* override if TR info is in description field */
    xtr = TRfromString(hdr.descrip);
    if (xtr > 0) {
      tr = (short) (xtr*factor);
      fprintf(stderr," reading TR from description field, TR= %ld ms\n",tr);
    }
    /* fprintf(stderr," nt=%ld,  TR= %ld milliseconds\n",nt,tr); */
    if (tr < 0.1) VWarning(" implausible TR (%d ms), use 'vnifti' for data conversion to specify TR on the command line",tr);
  }

  /* voxel reso */
  int blen=512;
  char *voxelstr = (char *) VCalloc(blen,sizeof(char));
  memset(voxelstr,0,blen);
  sprintf(voxelstr,"%f %f %f",hdr.pixdim[1],hdr.pixdim[2],hdr.pixdim[3]);

  /* geometry information */
  VAttrList geolist = VCreateAttrList();
  double *slicetime = VGetNiftiHeader(geolist,hdr,tr);

  /* read nii image into vista attrlist */
  VAttrList out_list = VCreateAttrList();
  VAppendAttr(out_list,"geoinfo",NULL,VAttrListRepn,geolist);

  (*ok) = FALSE;
  if (nt <= 1) {      /* output one 3D image */
    Nii2Vista3DList(data,nsize,nslices,nrows,ncols,nt,dst_repn,datatype,scl_slope,scl_inter,voxelstr,tr,out_list);
    VUpdateGeoinfo(geolist,(int)3,tr);
  }
  else if (attrtype == FALSE) {        /* output list of 3D images */
    Nii2Vista3DList(data,nsize,nslices,nrows,ncols,nt,dst_repn,datatype,scl_slope,scl_inter,voxelstr,tr,out_list);
    VUpdateGeoinfo(geolist,(int)3,tr);
  }
  else if (attrtype == TRUE && nt > 1) {   /* output one 4D image */
    Nii2Vista4D(data,nsize,nslices,nrows,ncols,nt,dst_repn,datatype,do_scaling,scl_slope,scl_inter,voxelstr,slicetime,tr,out_list);
    VUpdateGeoinfo(geolist,(int)4,tr);
    (*ok) = TRUE;
  }

  return out_list;
}

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
/* rng/ran1.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_rng.h>

/* This is an implementation of the algorithm used in Numerical
   Recipe's ran1 generator.  It is MINSTD with a 32-element
   shuffle-box. */

static inline unsigned long int ran1_get (void *vstate);
static double ran1_get_double (void *vstate);
static void ran1_set (void *state, unsigned long int s);

static const long int m = 2147483647, a = 16807, q = 127773, r = 2836;

#define N_SHUFFLE 32
#define N_DIV (1 + 2147483646/N_SHUFFLE)

typedef struct
  {
    unsigned long int x;
    unsigned long int n;
    unsigned long int shuffle[N_SHUFFLE];
  }
ran1_state_t;

static inline unsigned long int
ran1_get (void *vstate)
{
  ran1_state_t *state = (ran1_state_t *) vstate;

  const unsigned long int x = state->x;

  const long int h = x / q;
  const long int t = a * (x - h * q) - h * r;

  if (t < 0)
    {
      state->x = t + m;
    }
  else
    {
      state->x = t;
    }

  {
    unsigned long int j = state->n / N_DIV;
    state->n = state->shuffle[j];
    state->shuffle[j] = state->x;
  }

  return state->n;
}

static double
ran1_get_double (void *vstate)
{
  float x_max = 1 - 1.2e-7f ; /* Numerical Recipes version of 1-FLT_EPS */

  float x = ran1_get (vstate) / 2147483647.0f ;

  if (x > x_max)
    return x_max ;

  return x ;
}

static void
ran1_set (void *vstate, unsigned long int s)
{
  ran1_state_t *state = (ran1_state_t *) vstate;
  int i;

  if (s == 0)
    s = 1;      /* default seed is 1 */

  for (i = 0; i < 8; i++)
    {
      long int h = s / q;
      long int t = a * (s - h * q) - h * r;
      if (t < 0)
        t += m;
      s = t;
    }

  for (i = N_SHUFFLE - 1; i >= 0; i--)
    {
      long int h = s / q;
      long int t = a * (s - h * q) - h * r;
      if (t < 0)
        t += m;
      s = t;
      state->shuffle[i] = s;
    }

  state->x = s;
  state->n = s;

  return;
}

static const gsl_rng_type ran1_type =
{"ran1",                        /* name */
 2147483646,                    /* RAND_MAX */
 1,                             /* RAND_MIN */
 sizeof (ran1_state_t),
 &ran1_set,
 &ran1_get,
 &ran1_get_double};

const gsl_rng_type *gsl_rng_ran1 = &ran1_type;

ext: c

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
/* Header file for ridgeRegressionFunctions.c */

#include "depends.h"
#ifdef HAVE_GSL_HEADER

/* includes */
#if _CUDA_
#include <cuda.h>
#include <cuda_runtime_api.h>
#include <cublas.h>
#include <cula_blas.h>
#include <cula.h>
#include "cudaOnlyFunctions.h"
#endif

#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <math.h>

#include "computeLinearRidge.h"

#include "ReadInData.h"

/* Compute linear ridge coefficeints */
int computeLinearRidge(GSL_TYPE(vector) * ahat,
		       GSL_TYPE(vector) * B,
		       GSL_TYPE(vector) * D2,
		       GSL_TYPE(matrix) * V,
		       PREC lambda);

/* linear generalized ridge regression */
GSL_TYPE(vector) * computeLinearGeneralizedRidge(GSL_TYPE(vector) * beta,
					   GSL_TYPE(matrix) * pred,
					   GSL_TYPE(vector) * pheno,
					   GSL_TYPE(vector) * shrinkage,
					   int intercept_flag);

int computeLogisticRidge(GSL_TYPE(vector) * beta,
				GSL_TYPE(matrix) * pred,
				GSL_TYPE(vector) * pheno,
				GSL_TYPE(vector) * shrinkage,
				int intercept_flat,
				int Doff_flat,
				PREC * DofF);

PREC objectiveFunction(GSL_TYPE(vector) * beta,
				GSL_TYPE(matrix) * X,
				GSL_TYPE(vector) * pheno,
				GSL_TYPE(vector) * shrinkage,
				int intercept_flat);

int updateBeta(GSL_TYPE(vector) * beta,
			GSL_TYPE(matrix) * X,
			GSL_TYPE(vector) * pheno,
			GSL_TYPE(matrix) * kI,
			int intercept_flag,
			int DofF_flag,
			GSL_TYPE(matrix) * invtXWX_return,
			GSL_TYPE(matrix) * W_return);

int getProb(GSL_TYPE(vector) * p, GSL_TYPE(vector) * XB);

/* NB my_gsl_solve only works for matrices of type double
  due to the linalg functions only having been written for this type */
int my_gsl_solve(gsl_matrix * X,
	gsl_matrix * solvedX);

int compute_XB_and_p(GSL_TYPE(matrix) * X,
	GSL_TYPE(vector) * B,
	GSL_TYPE(vector) * XB,
	GSL_TYPE(vector) * p);

int chooseHowManyK(GSL_TYPE(vector) * D);

int returnToOriginalScaleLinear(GSL_TYPE(vector) * betaOut,
				GSL_TYPE(vector) * Bridge,
				GSL_TYPE(vector) * means,
				GSL_TYPE(vector) * scales,
				PREC y_mean,
				int intercept_flag);

/* return to original scale - generalized linear ridge regression */
int returnToOriginalScaleGenLinear(GSL_TYPE(vector) * Bridge,
					    GSL_TYPE(vector) * betaOut,
					    GSL_TYPE(matrix) * pred,
					    GSL_TYPE(vector) * pheno,
					    GSL_TYPE(vector) * scales,
					    int intercept_flag);

/* Convert trans */
char setTrans(CBLAS_TRANSPOSE_t Trans);

/* prepare for linear ridge
	used when we are going to call linear rr on a range of different shrinkage parameters */
int prepareForLinearRidge(GSL_TYPE(matrix) * X,
			  GSL_TYPE(vector) * y,
			  GSL_TYPE(matrix) * U,
			  GSL_TYPE(matrix) * V,
			  GSL_TYPE(vector) * D,
			  GSL_TYPE(vector) * D2,
			  GSL_TYPE(matrix) * Z,
			  GSL_TYPE(vector) * ahat);

/* compute DofF for ridge model */
int computeDofF(GSL_TYPE(vector) * D2,
		PREC Kr,
		PREC * DofF);

/* compute (A + BDC')^(-1) */
int invert_sum_of_matrices(const PREC Ainv,
			   const GSL_TYPE(matrix) * B,
			   const GSL_TYPE(vector) * Dinv,
			   const GSL_TYPE(matrix) * tC,
			   GSL_TYPE(matrix) * out);

/* ordinary regression */
gsl_vector * my_gsl_linear_fit(gsl_matrix * X,
			       gsl_vector * y,
			       int NROW,
			       int NCOL);

#endif

ext: h

max_stars_repo_name: dfrankow/ridge-cran

max_issues_repo_name: dfrankow/ridge-cran

max_forks_repo_name: dfrankow/ridge-cran


text:
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/block/gsl_block.h>
#include <gsl/vector/gsl_vector.h>

#define BASE_DOUBLE
#include <gsl/templates_on.h>
#include <gsl/vector/file_source.c>
#include <gsl/templates_off.h>
#undef BASE_DOUBLE

max_stars_repo_name: MontyThibault/centre-of-mass-awareness

max_issues_repo_name: MontyThibault/centre-of-mass-awareness

max_forks_repo_name: MontyThibault/centre-of-mass-awareness


text:
/**
* Copyright 2016 BitTorrent 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 <scraps/config.h>

#include <scraps/Byte.h>
#include <scraps/Temp.h>
#include <scraps/ByteArray.h>

#include <gsl.h>

#include <iterator>

namespace scraps {

constexpr signed char DecToHex(int c) {
    if (c >= 0 && c < 16) {
        return "0123456789abcdef"[c];
    }
    return -1;
}

constexpr signed char DecToHex(const Byte& c) {
    return DecToHex(c.value());
}

template <typename CharT>
constexpr int8_t HexToDec(CharT c) {
    if (c >= '0' && c <= '9') { return c - '0'; }
    if (c >= 'a' && c <= 'f') { return c - 'a' + 10; }
    if (c >= 'A' && c <= 'F') { return c - 'A' + 10; }
    return -1;
}

/**
 * Returns a hex string representation of the input span. The output does not
 * include an "0x" prefix.
 *
 * example: ToHex(std::array<uint8_t, 1>{0xAB}) == "AB";
 */
template <typename T, std::ptrdiff_t... BytesDimension>
std::string ToHex(const gsl::span<T, BytesDimension...> range) {
    std::string ret;

    static_assert(sizeof(T) == 1, "Input span type too large");

    ret.reserve(range.size() * 2 + 2);

    for (auto& b : range) {
        ret += DecToHex(b >> 4);
        ret += DecToHex(b & 0x0F);
    }

    return ret;
}

template <typename T, size_t N>
std::string ToHex(const std::array<T, N>& in) {
    return ToHex(gsl::span<const T, N>{in});
}

template <size_t N>
std::string ToHex(const ByteArray<N>& in) {
    return ToHex(gsl::span<const unsigned char, N>{in.bytes});
}

/**
 * Fills a byte range from a range of hex characters. Returns true if the input
 * range is successfully converted to the output range. A prefix of "0x" or "0X"
 * is optional.
 */
template <typename HexCharT, std::ptrdiff_t HexExtent, typename ByteType, std::ptrdiff_t... BytesDimension>
constexpr bool ToBytes(const gsl::basic_string_span<const HexCharT, HexExtent> hex, const gsl::span<ByteType, BytesDimension...> bytes) {
    auto prefixSize = 0;
    if (hex.size() > 2 && hex[0] == '0' && (hex[1] == 'x' || hex[1] == 'X')) {
        prefixSize += 2;
    }
    if (hex.size() - prefixSize != bytes.size() * 2) {
        return false;
    }
    for (decltype(bytes.size()) i = 0; i < bytes.size(); ++i) {
        auto hi = HexToDec(hex[i * 2 + prefixSize]);
        auto lo = HexToDec(hex[i * 2 + 1 + prefixSize]);
        if (hi < 0 || lo < 0) {
            return false;
        }
        bytes[i] = ByteType{static_cast<uint8_t>((hi << 4) | lo)};
    }
    return true;
}

template <typename CharT, typename T, size_t N>
constexpr bool ToBytes(const std::basic_string<CharT>& in, std::array<T, N>& out) {
    return ToBytes(gsl::basic_string_span<const CharT>{in}, gsl::span<T, N>{out});
}

} // namespace scraps

max_stars_repo_name: marakew/scraps

max_issues_repo_name: marakew/scraps

max_forks_repo_name: marakew/scraps


text:
/**
 * @file binauralprocessing.h
 * @brief Binaural processing
 * @author Kenichi Kumatani
 */
#ifndef BINAURALPROCESSING_H
#define BINAURALPROCESSING_H

#include <stdio.h>
#include <assert.h>
#include <float.h>

#include <gsl/gsl_block.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_fft_complex.h>
#include <common/refcount.h>
#include "common/jexception.h"

#include "stream/stream.h"
#include "postfilter/spectralsubtraction.h"
#include "beamformer/spectralinfoarray.h"
#include "beamformer/beamformer.h"

class BinaryMaskFilter : public VectorComplexFeatureStream {
public:
  BinaryMaskFilter( unsigned chanX, VectorComplexFeatureStreamPtr &srcL, VectorComplexFeatureStreamPtr &srcR,
		       unsigned M, float threshold, float alpha,
		       float dEta = 0.01, const String& nm = "BinaryMaskFilter" );
  ~BinaryMaskFilter();
  virtual const gsl_vector_complex* next(int frame_no = -5);
  virtual void reset();
  void set_threshold( float threshold ){ threshold_ = threshold; }
  void set_thresholds( const gsl_vector *thresholds );
  double threshold(){ return threshold_; }
  gsl_vector *thresholds(){ return threshold_per_freq_; }

#ifdef ENABLE_LEGACY_BTK_API
  void setThreshold( float threshold ){ set_threshold(threshold); }
  void setThresholds( const gsl_vector *thresholds ){ set_thresholds(thresholds); }
  double getThreshold(){ return threshold(); }
  gsl_vector *getThresholds(){ return thresholds(); }
#endif

protected:
  VectorComplexFeatureStreamPtr srcL_; /* left channel */
  VectorComplexFeatureStreamPtr srcR_; /* right channel */
  unsigned          chanX_;  /* want to extract the index of a channel */
  gsl_vector_float *prevMu_; /* binary mask at a previous frame */
  float             alpha_;  /* forgetting factor */
  float             dEta_;        /* flooring value */
  float             threshold_;   /* threshold for the ITD */
  gsl_vector        *threshold_per_freq_;
};

typedef Inherit<BinaryMaskFilter, VectorComplexFeatureStreamPtr> BinaryMaskFilterPtr;

/**
   @class Implementation of binary masking based on C. Kim's Interspeech2010 paper
   @brief binary mask two inputs based of the threshold of the interaural time delay (ITD)
   @usage
 */
class KimBinaryMaskFilter : public BinaryMaskFilter {
public:
  KimBinaryMaskFilter( unsigned chanX, VectorComplexFeatureStreamPtr &srcL, VectorComplexFeatureStreamPtr &srcR,
		       unsigned M, float threshold, float alpha,
		       float dEta = 0.01, float dPowerCoeff = 1/15.0, const String& nm = "KimBinaryMaskFilter" );
  ~KimBinaryMaskFilter();
  virtual const gsl_vector_complex* next(int frame_no = -5);
  virtual void reset();

  virtual const gsl_vector_complex* masking1( const gsl_vector_complex* ad_X_L, const gsl_vector_complex* ad_X_R, float threshold );

protected:
  float             dpower_coeff_; /* power law non-linearity */
};

typedef Inherit<KimBinaryMaskFilter, BinaryMaskFilterPtr> KimBinaryMaskFilterPtr;

/**
   @class Implementation of estimating the threshold for C. Kim's ITD-based binary masking
   @brief binary mask two inputs based of the threshold of the interaural time delay (ITD)
   @usage
 */
class KimITDThresholdEstimator : public KimBinaryMaskFilter {
public:
  KimITDThresholdEstimator( VectorComplexFeatureStreamPtr &srcL, VectorComplexFeatureStreamPtr &srcR, unsigned M, float minThreshold = 0, float maxThreshold = 0, float width = 0.02, float minFreq= -1, float maxFreq=-1, int sampleRate=-1, float dEta = 0.01, float dPowerCoeff = 1/15.0, const String& nm = "KimITDThresholdEstimator" );
  ~KimITDThresholdEstimator();
  virtual const gsl_vector_complex* next(int frame_no = -5);
  virtual void   reset();
  virtual double calc_threshold();
  const gsl_vector* cost_function();

#ifdef ENABLE_LEGACY_BTK_API
  virtual double calcThreshold(){ return calc_threshold(); }
  const gsl_vector* getCostFunction(){ return cost_function(); }
#endif

protected:
  virtual void accumStats1( const gsl_vector_complex* ad_X_L, const gsl_vector_complex* ad_X_R );

  /* for restricting the search space */
  float min_threshold_;
  float max_threshold_;
  float width_;
  unsigned min_fbinX_;
  unsigned max_fbinX_;
  /* work space */
  double *cost_func_values_;
  double *sigma_T_;
  double *sigma_I_;
  double *mean_T_;
  double *mean_I_;
  unsigned int nCand_;
  unsigned int nSamples_;
  gsl_vector *buffer_; /* for returning values from the python */
  bool cost_func_computed_;
};

typedef Inherit<KimITDThresholdEstimator, KimBinaryMaskFilterPtr> KimITDThresholdEstimatorPtr;

/**
	@class
 */
class IIDBinaryMaskFilter : public BinaryMaskFilter {
public:
	IIDBinaryMaskFilter( unsigned chanX, VectorComplexFeatureStreamPtr &srcL, VectorComplexFeatureStreamPtr &srcR,
						unsigned M, float threshold, float alpha,
						float dEta = 0.01, const String& nm = "IIDBinaryMaskFilter" );
	~IIDBinaryMaskFilter();
	virtual const gsl_vector_complex* next(int frame_no = -5);
	virtual void reset();

	virtual const gsl_vector_complex* masking1( const gsl_vector_complex* ad_X_L, const gsl_vector_complex* ad_X_R, float threshold );
};

typedef Inherit<IIDBinaryMaskFilter, BinaryMaskFilterPtr> IIDBinaryMaskFilterPtr;

/**
 @class binary masking based on a difference between magnitudes of two beamformers' outputs.
 @brief set zero to beamformer's output with the smaller magnitude over frequency bins
 */
class IIDThresholdEstimator : public KimITDThresholdEstimator {
public:
	IIDThresholdEstimator( VectorComplexFeatureStreamPtr &srcL, VectorComplexFeatureStreamPtr &srcR, unsigned M,
							 float minThreshold = 0, float maxThreshold = 0, float width = 0.02,
							 float minFreq= -1, float maxFreq=-1, int sampleRate=-1, float dEta = 0.01, float dPowerCoeff = 0.5, const String& nm = "IIDThresholdEstimator" );
	~IIDThresholdEstimator();
	virtual const gsl_vector_complex* next(int frame_no = -5);
	virtual void reset();
	virtual double calc_threshold();

#ifdef ENABLE_LEGACY_BTK_API
	virtual double calcThreshold(){ return calc_threshold(); }
#endif

protected:
        virtual void accumStats1( const gsl_vector_complex* ad_X_L, const gsl_vector_complex* ad_X_R );

private:
        double *_Y4_T;
        double *_Y4_I;
        double _beta;
};

typedef Inherit<IIDThresholdEstimator, KimITDThresholdEstimatorPtr> IIDThresholdEstimatorPtr;

/**
 @class binary masking based on a difference between magnitudes of two beamformers' outputs at each frequency bin.
 @brief set zero to beamformer's output with the smaller magnitude at each frequency bin
 */
class FDIIDThresholdEstimator : public BinaryMaskFilter {
public:
  FDIIDThresholdEstimator( VectorComplexFeatureStreamPtr &srcL, VectorComplexFeatureStreamPtr &srcR, unsigned M,
			   float minThreshold = 0, float maxThreshold = 0, float width = 1000,
			   float dEta = 0.01, float dPowerCoeff = 1/15.0, const String& nm = "FDIIDThresholdEstimator" );
  ~FDIIDThresholdEstimator();
  virtual const gsl_vector_complex* next(int frame_no = -5);
  virtual void   reset();
  virtual double calc_threshold();
  const gsl_vector* cost_function(unsigned freqX);

#ifdef ENABLE_LEGACY_BTK_API
  virtual double calcThreshold(){ return calc_threshold(); }
  const gsl_vector* getCostFunction( unsigned freqX ){ return cost_function(freqX); }
#endif

protected:
  virtual void accumStats1( const gsl_vector_complex* ad_X_L, const gsl_vector_complex* ad_X_R );

  /* for restricting the search space */
  float min_threshold_;
  float max_threshold_;
  float width_;
  float dpower_coeff_;

  /* work space */
  double **cost_func_values_;
  double **_Y4;
  double **_sigma;
  double **_mean;
  double _beta;
  unsigned int nCand_;
  unsigned int nSamples_;
  gsl_vector *buffer_; /* for returning values from the python */
  bool cost_func_computed_;
};

typedef Inherit<FDIIDThresholdEstimator, BinaryMaskFilterPtr> FDIIDThresholdEstimatorPtr;

#endif


text:
#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include <gsl/gsl_blas.h>
#include "my_timer.h"

#define N 5

int cal_index_offset(int ndim, int *darr, int *position) {
    int size = 1;
    int offset = 0;
    for(int i = ndim-1; i >= 0; i--) {
        offset += position[i] * size;
        size *= darr[i];
    }
    return offset;
}

void print_mat_recursive(double *mat, int ndim, int *darr, int pdim, int *position) {
    if(pdim >= ndim) {
        int index_offset = cal_index_offset(ndim, darr, position);
        printf("%lf ", mat[index_offset]);
        return;
    }

    printf("[ ");
    for(int i = 0; i < darr[pdim]; i++) {
        position[pdim] = i;
        print_mat_recursive(mat, ndim, darr, pdim+1, position);
    }
    printf("]");
}

void print_mat(double *mat, int ndim, ...) {
    int position[ndim];
    int darr[ndim];

    va_list args;
    va_start(args, ndim);
    for(int i = 0; i < ndim; i++) {
        darr[i] = va_arg(args, int);
    }
    va_end(args);

    print_mat_recursive(mat, ndim, darr, 0, position);
    printf("\n");
}

void mm_brute_force(double *a, double *b, double *c, int p, int q, int r) {
    for(int i = 0; i < p; i++) {
        for(int k = 0; k < q; k++) {
            for(int j = 0; j < r; j++) {
                c[i * r + j] += a[i * q + k] * b[k * r + j];
            }
        }
    }
}

void mm_dgemm(double *a, double *b, double *c, int p, int q, int r) {
    int l = p;
    int m = q;
    int n = r;

    int lda = m;
    int ldb = n;
    int ldc = n;

    double alpha = 1.0;
    double beta = 0.0;

    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
            l, n, m, alpha, a, lda, b, ldb, beta, c, ldc);
}

void example_dgemm() {
    int l, m, n;
    l = m = n = N;

    double *mat_a = (double*)malloc(sizeof(double) * l * m);
    double *mat_b = (double*)malloc(sizeof(double) * m * n);
    double *mat_c = (double*)malloc(sizeof(double) * l * n);

    for(int i = 0; i < (l*m); i++) {
        mat_a[i] = i + 1;
    }

    for(int i = 0; i < (m*n); i++) {
        mat_b[i] = i + 1;
    }

    for(int i = 0; i < (l*n); i++) {
        mat_c[i] = 0;
    }

    my_timer_t mt;
    my_timer_start(&mt);
    mm_dgemm(mat_a, mat_b, mat_c, l, m, n);
    my_timer_end(&mt);
    my_timer_print(&mt);
    printf("\n");

    print_mat(mat_a, 2, l, m);
    printf("\n");
    print_mat(mat_b, 2, m, n);
    printf("\n");
    print_mat(mat_c, 2, l, n);
    printf("\n");

    free(mat_a);
    free(mat_b);
    free(mat_c);
}

int main(void) {
    example_dgemm();
    return 0;
}


text:
#pragma once

#include "type_traits.h"

#include <gsl/gsl>

#include <system_error>
#include <variant>

namespace arcana
{
    template<typename Type>
    struct is_expected;

    template<typename T, typename E>
    class basic_expected;

    namespace internal
    {
        template<typename T>
        struct expected_error_traits;

        template<>
        struct expected_error_traits<std::error_code>
        {
            template<typename ErrorT>
            using is_convertible = std::bool_constant<
                std::is_error_code_enum<std::decay_t<ErrorT>>::value ||
                std::is_error_condition_enum<std::decay_t<ErrorT>>::value
            >;

            template<typename ErrorT>
            using enable_conversion = std::enable_if<
                is_convertible<ErrorT>::value
            >;

            template<typename ErrorT>
            static std::error_code convert(ErrorT errorEnum)
            {
                using std::make_error_code;

                return make_error_code(errorEnum);
            }
        };

        template<>
        struct expected_error_traits<std::exception_ptr>
        {
            template<typename ErrorT>
            using enable_conversion = std::enable_if<
                std::is_same<std::error_code, ErrorT>::value || // support conversion from error_codes
                expected_error_traits<std::error_code>::template is_convertible<ErrorT>::value // and also enumerations that are error codes
            >;

            static std::exception_ptr convert(std::error_code code)
            {
                return std::make_exception_ptr(std::system_error(code));
            }

            template<typename ErrorT, typename = typename expected_error_traits<std::error_code>::template enable_conversion<ErrorT>::type>
            static std::exception_ptr convert(ErrorT errorEnum)
            {
                return convert(
                    expected_error_traits<std::error_code>::convert(errorEnum));
            }
        };
    }

    class bad_expected_access : public std::exception
    {
    public:
        const char* what() const override
        {
            return "tried accessing value()/error() of an expected when it wasn't set";
        }
    };

    template<typename E>
    class unexpected
    {
    public:
        unexpected() = delete;

        constexpr explicit unexpected(const E& e)
            : m_e{ e }
        {}

        constexpr explicit unexpected(E&& e)
            : m_e{ std::move(e) }
        {}

        constexpr const E& value() const &
        {
            return m_e;
        }

        constexpr E& value() &
        {
            return m_e;
        }

        constexpr E&& value() &&
        {
            return std::move(m_e);
        }
    private:
        E m_e;
    };

    template<typename T>
    inline auto make_unexpected(T&& error)
    {
        return unexpected<std::decay_t<T>>{std::forward<T>(error)};
    }

    template<typename T, typename E>
    class basic_expected
    {
        using traits = internal::expected_error_traits<E>;

        template<typename E2,
            typename = typename traits::template enable_conversion<E2>::type>
        std::variant<E, T> convert_variant(const basic_expected<T, E2>& exp)
        {
            if (exp.has_error())
                return { traits::convert(exp.error()) };
            else
                return { exp.value() };
        }
    public:
        using value_type = T;
        using error_type = E;

        basic_expected(const value_type& value)
            : m_data{ value }
        {}

        basic_expected(value_type&& value)
            : m_data{ std::move(value) }
        {}

        template<typename ErrorT,
            typename = typename traits::template enable_conversion<ErrorT>::type>
        basic_expected(const unexpected<ErrorT>& errorEnum)
            : m_data{ traits::convert(errorEnum.value()) }
        {}

        basic_expected(const unexpected<error_type>& ec)
            : m_data{ ec.value() }
        {
            GSL_CONTRACT_CHECK("you should never build an basic_expected<T,E> with a non-error", static_cast<bool>(error()));
        }

        basic_expected(unexpected<error_type>&& ec)
            : m_data{ std::move(ec).value() }
        {
            GSL_CONTRACT_CHECK("you should never build an basic_expected<T,E> with a non-error", static_cast<bool>(error()));
        }

        template<typename ErrorT,
            typename = typename traits::template enable_conversion<ErrorT>::type>
        basic_expected(const basic_expected<value_type, ErrorT>& other)
            : m_data{ convert_variant(other) }
        {}

        template<typename ErrorT,
            typename = typename traits::template enable_conversion<ErrorT>::type>
        basic_expected& operator=(const basic_expected<value_type, ErrorT>& other)
        {
            m_data = convert_variant(other);
        }

        //
        // Returns the expected value or throws bad_expected_access if it is in an error state.
        //
        const value_type& value() const
        {
            return const_cast<basic_expected*>(this)->value();
        }

        value_type& value()
        {
            auto val = std::get_if<value_type>(&m_data);
            if (val == nullptr)
            {
                throw bad_expected_access();
            }
            return *val;
        }

        //
        // Returns the value contained or the default passed in if the expected is in an error state.
        //
        template<typename U>
        value_type value_or(U&& def) const
        {
            if (!has_value())
            {
                return std::forward<U>(def);
            }

            return value();
        }

        //
        // Returns the error or throws bad_expected_access if it is not in an error state.
        //
        const error_type& error() const
        {
            auto err = std::get_if<error_type>(&m_data);
            if (err == nullptr)
            {
                throw bad_expected_access();
            }
            return *err;
        }

        //
        // Returns whether or not the expected contains a value.
        //
        bool has_value() const noexcept
        {
            return std::holds_alternative<value_type>(m_data);
        }

        //
        // Returns whether or not the expected is in an error state.
        //
        bool has_error() const noexcept
        {
            return std::holds_alternative<error_type>(m_data);
        }

        //
        // Converts to true if the value is valid.
        //
        explicit operator bool() const noexcept
        {
            return has_value();
        }

        //
        // Shorthand access operators
        //
        const value_type& operator*() const
        {
            return value();
        }

        value_type& operator*()
        {
            return value();
        }

        //
        // Shorthand access operators
        //
        const value_type* operator->() const
        {
            return &value();
        }

        value_type* operator->()
        {
            return &value();
        }

    private:
        std::variant<error_type, value_type> m_data;
    };

    template<typename E>
    class basic_expected<void, E>
    {
    public:
        using value_type = void;
        using error_type = E;

        using traits = internal::expected_error_traits<error_type>;

        template<typename ErrorT,
            typename = typename traits::template enable_conversion<ErrorT>::type>
        basic_expected(const unexpected<ErrorT>& errorEnum)
            : m_error{ traits::convert(errorEnum.value()) }
        {}

        basic_expected(const unexpected<error_type>& ec)
            : m_error{ ec.value() }
        {
            GSL_CONTRACT_CHECK("you should never build an basic_expected<T,E> with a non-error", static_cast<bool>(error()));
        }

        basic_expected(unexpected<error_type>&& ec)
            : m_error{ std::move(ec).value() }
        {
            GSL_CONTRACT_CHECK("you should never build an basic_expected<T,E> with a non-error", static_cast<bool>(error()));
        }

        template<typename ErrorT,
            typename = typename traits::template enable_conversion<ErrorT>::type>
        basic_expected(const basic_expected<value_type, ErrorT>& other)
            : m_error{ other.has_error() ? traits::convert(other.error()) : nullptr }
        {}

        template<typename ErrorT,
            typename = typename traits::template enable_conversion<ErrorT>::type>
        basic_expected& operator=(const basic_expected<value_type, ErrorT>& other)
        {
            m_error = other.has_error() ? traits::convert(other.error()) : nullptr;
        }

        //
        // Returns the error, throws bad_expected_access if it is not in an error state
        //
        const error_type& error() const
        {
            if (!has_error())
            {
                throw bad_expected_access();
            }
            return m_error;
        }

        //
        // Returns whether or not the expected is in an error state.
        //
        bool has_error() const noexcept
        {
            return static_cast<bool>(m_error);
        }

        //
        // Converts to true if the value is valid.
        //
        explicit operator bool() const noexcept
        {
            return !has_error();
        }

        //
        // Creates an expected<void> that doesn't have an error
        //
        static basic_expected make_valid()
        {
            return basic_expected{};
        }

    private:
        basic_expected() = default;

        error_type m_error;
    };

    template<typename Type>
    struct is_expected : public std::false_type
    {};

    template<typename Type, typename Error>
    struct is_expected<basic_expected<Type, Error>> : public std::true_type
    {};

    template<typename ResultT, typename ErrorT>
    struct as_expected
    {
        using type = basic_expected<ResultT, ErrorT>;
        using value_type = typename type::value_type;
        using error_type = typename type::error_type;
    };

    template<typename ResultT, typename ErrorT, typename DesiredErrorT>
    struct as_expected<basic_expected<ResultT, ErrorT>, DesiredErrorT>
    {
        static_assert(std::is_same<ErrorT, DesiredErrorT>::value, "the error type of the expected and the desired error type don't match");

        using type = basic_expected<ResultT, ErrorT>;
        using value_type = typename type::value_type;
        using error_type = typename type::error_type;
    };

    template<typename ResultT, typename ErrorT>
    using expected = basic_expected<ResultT, ErrorT>;

    template<typename ErrorT>
    struct error_priority;

    template<>
    struct error_priority<std::error_code> : std::integral_constant<int, 0>
    {
        using type = std::error_code;
    };

    template<>
    struct error_priority<std::exception_ptr> : std::integral_constant<int, 1>
    {
        using type = std::exception_ptr;
    };

    template<typename Left, typename Right>
    struct largest_error
    {
        using type = typename typename largest_integral_constant<error_priority<Left>, error_priority<Right>>::type::type;
    };

    namespace internal
    {
        template<typename Expected, typename OrType, bool IsExpected>
        struct expected_error_or_impl
        {
            using type = OrType;
        };

        template<typename Expected, typename OrType>
        struct expected_error_or_impl<Expected, OrType, true>
        {
            using type = typename Expected::error_type;
        };
    }

    template<typename Expected, typename OrType>
    struct expected_error_or : public internal::expected_error_or_impl<Expected, OrType, is_expected<Expected>::value>
    {};
}


text:
#ifndef GEMMFORGE_REFERENCE_GEMM_H
#define GEMMFORGE_REFERENCE_GEMM_H

#include "typedef.h"
#include <iostream>

#define GEMMFORGE 1
#define OPENBLAS 2

#define CPU_BACKEND CONCRETE_CPU_BACKEND

#if CPU_BACKEND == OPENBLAS
#include <cblas.h>
#endif

namespace gemmforge {
  namespace reference {
    enum class LayoutType {
      Trans, NoTrans
    };

    void singleGemm(LayoutType TypeA,
                    LayoutType TypeB,
                    int M, int N, int K,
                    real Alpha, real *A, int Lda,
                    real *B, int Ldb,
                    real Beta, real *C, int Ldc);

    real *findData(real *Data, unsigned Stride, unsigned BlockId);
    real *findData(real **Data, unsigned Stride, unsigned BlockId);

    template<typename AT, typename BT, typename CT>
    void gemm(LayoutType TypeA,
              LayoutType TypeB,
              int M, int N, int K,
              real Alpha, AT A, int Lda,
              BT B, int Ldb,
              real Beta, CT C, int Ldc,
              unsigned OffsetA,
              unsigned OffsetB,
              unsigned OffsetC,
              unsigned NumElements) {

      for (unsigned Index = 0; Index < NumElements; ++Index) {
        real *MatrixA = findData(A, OffsetA, Index);
        real *MatrixB = findData(B, OffsetB, Index);
        real *MatrixC = findData(C, OffsetC, Index);

#if CPU_BACKEND == GEMMFORGE
        singleGemm(TypeA, TypeB,
                   M, N, K,
                   Alpha, MatrixA, Lda,
                   MatrixB, Ldb,
                   Beta, MatrixC, Ldc);

#elif CPU_BACKEND == OPENBLAS
  CBLAS_LAYOUT Layout = CblasColMajor;
  CBLAS_TRANSPOSE TransA = TypeA == LayoutType::Trans ? CblasTrans : CblasNoTrans;
  CBLAS_TRANSPOSE TransB = TypeB == LayoutType::Trans ? CblasTrans : CblasNoTrans;

#if REAL_SIZE == 4
  cblas_sgemm(Layout, TransA, TransB,
              M, N, K,
              Alpha, MatrixA, Lda,
              MatrixB, Ldb,
              Beta, MatrixC, Ldc);
#elif REAL_SIZE == 8
  cblas_dgemm(Layout, TransA, TransB,
              M, N, K,
              Alpha, MatrixA, Lda,
              MatrixB, Ldb,
              Beta, MatrixC, Ldc);
#endif

#else
#error "Chosen reference CPU-GEMM impl. is not supported"
#endif
      }
    }
  }
}

#endif //GEMMFORGE_REFERENCE_GEMM_H

max_stars_repo_name: ravil-mobile/gemmforge

max_issues_repo_name: ravil-mobile/gemmforge

max_forks_repo_name: ravil-mobile/gemmforge


text:
//
// SPDX-License-Identifier: MIT
// Copyright (C) 2020 - 2021 by the ryujin authors
//

#pragma once

#ifdef DEAL_II_WITH_GSL
#include <gsl/gsl_spline.h>

namespace ryujin
{
  /**
   * A cubic spline class implemented as a thin wrapper around the GSL
   * spline library functions.
   *
   * Usage:
   * @code
   * std::vector<double> xs{0.0, 0.2, 0.4, 0.6, 0.8, 1.0};
   * std::vector<double> ys{1.0, 0.2, 5.0, 2.0, 1.0, 10.0};
   * CubicSpline spline(xs, ys);
   *
   * spline.eval(0.5);
   * @endcode
   *
   * @ingroup Mesh
   */
  class CubicSpline
  {
  public:
    /**
     * Constructor.
     *
     * @pre The supplied vectors @p x and @p y must have the same size and
     * must contain at least two elements. The vector @p x must be sorted
     */
    CubicSpline(const std::vector<double> &x,
                const std::vector<double> &y) noexcept
        : x_(x)
        , y_(y)
    {
      AssertNothrow(x_.size() == y_.size(), dealii::ExcInternalError());
      AssertNothrow(x_.size() >= 2, dealii::ExcInternalError());
      AssertNothrow(std::is_sorted(x_.begin(), x_.end()),
                    dealii::ExcInternalError());

      spline = gsl_spline_alloc(gsl_interp_cspline, x_.size());
      gsl_spline_init(spline, x_.data(), y_.data(), x_.size());

      accel = gsl_interp_accel_alloc();
    }

    /**
     * Copy constructor.
     */
    CubicSpline(const CubicSpline &copy)
        : CubicSpline(copy.x_, copy.y_)
    {
    }

    /**
     * The copy assignment operator is deleted.
     */
    CubicSpline &operator=(const CubicSpline &) = delete;

    /**
     * Destructor.
     */
    ~CubicSpline()
    {
      gsl_interp_accel_free(accel);
      gsl_spline_free(spline);
    }

    /**
     * Evaluate the cubic spline at a given point @p x.
     *
     * @pre The point @p x must lie within the interval described by the
     * largest and smallest support point supplied to the constructor.
     */
    inline double eval(double x) const
    {
      return gsl_spline_eval(spline, x, accel);
    }

  private:
    const std::vector<double> x_;
    const std::vector<double> y_;
    gsl_spline *spline;
    mutable gsl_interp_accel *accel;
  };
} // namespace ryujin

#endif

max_stars_repo_name: Fanxiaotsing/ryujin

max_issues_repo_name: Fanxiaotsing/ryujin

max_forks_repo_name: Fanxiaotsing/ryujin


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_eigen.h>
#include "shared/libklp_matrix_header.h"
#include "shared/libtpl_header.h"
#include "population_constants.h"
#include "population_params.h"
#include "population_initializers.h"
#include "population_functions.h"

void population_from_row_ordered_transition_matrix(const KLP_PARAMS klp_params, const POPULATION_PARAMS parameters, TRANSITION_MATRIX row_transition_matrix) {
  if (parameters.equilibrium) {
    equilibrium_from_row_ordered_transition_matrix(klp_params, parameters, row_transition_matrix);
  } else {
    population_proportion_from_row_ordered_transition_matrix(klp_params, parameters, row_transition_matrix);
  }
}

void equilibrium_from_row_ordered_transition_matrix(const KLP_PARAMS klp_params, const POPULATION_PARAMS parameters, TRANSITION_MATRIX row_transition_matrix) {
  EIGENSYSTEM eigensystem;

  eigensystem = eigensystem_from_row_ordered_transition_matrix(row_transition_matrix);
  print_equilibrium(eigensystem, klp_params, parameters);

}

void population_proportion_from_row_ordered_transition_matrix(const KLP_PARAMS klp_params, const POPULATION_PARAMS parameters, TRANSITION_MATRIX row_transition_matrix) {
  EIGENSYSTEM eigensystem;

  eigensystem = eigensystem_from_row_ordered_transition_matrix(row_transition_matrix);
  print_population_proportion(eigensystem, klp_params, parameters);
}

EIGENSYSTEM eigensystem_from_row_ordered_transition_matrix(TRANSITION_MATRIX row_transition_matrix) {
  TRANSITION_MATRIX column_transition_matrix;
  EIGENSYSTEM eigensystem;

  column_transition_matrix = transpose_matrix(row_transition_matrix);
  eigensystem              = convert_transition_matrix_to_eigenvectors(column_transition_matrix);
  invert_matrix(&eigensystem);

  return eigensystem;
}

TRANSITION_MATRIX convert_structures_to_transition_matrix(const SOLUTION* all_structures, int num_structures) {
  // This code to convert the list of structures to a transition matrix does not ensure that there are only single b.p. moves.
  int i, j;
  double col_sum;
  TRANSITION_MATRIX transition_matrix = init_transition_matrix(num_structures, 'R');

  for (i = 0; i < num_structures; ++i) {
    col_sum = 0;

    for (j = 0; j < num_structures; ++j) {
      if (i != j) {
        T_COL_ORDER(transition_matrix, i, j) = \
                                               MIN(1, exp(-((double)all_structures[j].energy - (double)all_structures[i].energy) / VIENNA_RT));

        col_sum += T_COL_ORDER(transition_matrix, i, j);
#ifdef INSANE_DEBUG
        printf("%d\t%d\t%.4e\n", i, j, T_COL_ORDER(transition_matrix, i, j));
#endif
      }
    }

#ifdef INSANE_DEBUG
    printf("%d col_sum:\t%.4e\n\n", i, col_sum);
#endif
    T_COL_ORDER(transition_matrix, i, i) = -col_sum;
  }

  return transition_matrix;
}

EIGENSYSTEM convert_transition_matrix_to_eigenvectors(TRANSITION_MATRIX transition_matrix) {
  // Generates right eigenvectors from a column-ordered transition rate matrix, and saves them in a one-dimensional column-ordered array (eigensystem.vectors).
  int i, j;
  EIGENSYSTEM eigensystem;
  eigensystem                             = init_eigensystem(transition_matrix.row_length);
  gsl_matrix_view matrix_view             = gsl_matrix_view_array(transition_matrix.matrix, eigensystem.length, eigensystem.length);
  gsl_vector_complex* eigenvalues         = gsl_vector_complex_alloc(eigensystem.length);
  gsl_matrix_complex* eigenvectors        = gsl_matrix_complex_alloc(eigensystem.length, eigensystem.length);
  gsl_eigen_nonsymmv_workspace* workspace = gsl_eigen_nonsymmv_alloc(eigensystem.length);
  gsl_eigen_nonsymmv_params(1, workspace);
  gsl_eigen_nonsymmv(&matrix_view.matrix, eigenvalues, eigenvectors, workspace);
  gsl_eigen_nonsymmv_free(workspace);

  for (i = 0; i < eigensystem.length; ++i) {
    eigensystem.values[i]               = GSL_REAL(gsl_vector_complex_get(eigenvalues, i));
    gsl_vector_complex_view eigenvector = gsl_matrix_complex_column(eigenvectors, i);

    for (j = 0; j < eigensystem.length; ++j) {
      E_COL_ORDER(eigensystem.vectors, i, j, eigensystem.length) = GSL_REAL(gsl_vector_complex_get(&eigenvector.vector, j));
    }
  }

  free_transition_matrix(transition_matrix);
  gsl_vector_complex_free(eigenvalues);
  gsl_matrix_complex_free(eigenvectors);
  return eigensystem;
}

void invert_matrix(EIGENSYSTEM* eigensystem) {
  int i, j, signum;
  gsl_matrix* matrix_to_invert = gsl_matrix_alloc(eigensystem->length, eigensystem->length);
  gsl_matrix* inversion_matrix = gsl_matrix_alloc(eigensystem->length, eigensystem->length);
  gsl_permutation* permutation = gsl_permutation_alloc(eigensystem->length);

  for (i = 0; i < eigensystem->length; ++i) {
    for (j = 0; j < eigensystem->length; ++j) {
      gsl_matrix_set(matrix_to_invert, i, j, E_ROW_ORDER(eigensystem->vectors, i, j, eigensystem->length));
    }
  }

  gsl_linalg_LU_decomp(matrix_to_invert, permutation, &signum);
  gsl_linalg_LU_invert(matrix_to_invert, permutation, inversion_matrix);

  for (i = 0; i < eigensystem->length; ++i) {
    for (j = 0; j < eigensystem->length; ++j) {
      E_ROW_ORDER(eigensystem->inverse_vectors, i, j, eigensystem->length) = gsl_matrix_get(inversion_matrix, i, j);
    }
  }

  gsl_matrix_free(matrix_to_invert);
  gsl_matrix_free(inversion_matrix);
  gsl_permutation_free(permutation);
}

double probability_at_logtime(const EIGENSYSTEM eigensystem, const KLP_PARAMS klp_params, double timepoint, int start_index, int target_index) {
  // This function is hard-wired to only consider the kinetics for folding from a distribution where p_{0}(start_index) == 1.
  int i;
  double cumulative_probability = 0;

  if (klp_params.start_state < 0 || klp_params.end_state < 0) {
    if (klp_params.start_state < 0) {
      fprintf(stderr, "Error: the starting structure could not be found. Usually this means that you did not specify an explicit starting structure so the empty structure was used, but the energy band was not wide enough for RNAsubopt to sample it.\n");
    }

    if (klp_params.end_state < 0) {
      fprintf(stderr, "Error: the ending structure could not be found. Usually this means that you did explicitly provided a suboptimal ending structure, but the energy band was not wide enough for RNAsubopt to sample it.\n");
    }

    abort();
  }

  for (i = 0; i < eigensystem.length; ++i) {
    cumulative_probability +=
      E_ROW_ORDER(eigensystem.vectors, target_index, i, eigensystem.length) *
      E_ROW_ORDER(eigensystem.inverse_vectors, i, start_index, eigensystem.length) *
      exp(eigensystem.values[i] * pow(10, timepoint));
  }

  return cumulative_probability;
}

void find_key_structure_indices_in_structure_list(KLP_PARAMS* klp_params, const POPULATION_PARAMS parameters, const SOLUTION* all_structures, int num_structures) {
  int i;

  for (i = 0; i < num_structures; ++i) {
    if (klp_params->start_state == -1) {
      if (!strcmp(parameters.start_structure, all_structures[i].structure)) {
        klp_params->start_state = i;
      }
    }

    if (klp_params->end_state == -1) {
      if (!strcmp(parameters.end_structure, all_structures[i].structure)) {
        klp_params->end_state = i;
      }
    }
  }
}

void serialize_eigensystem(const EIGENSYSTEM eigensystem, const POPULATION_PARAMS parameters) {
  tpl_node* tpl;
  int i;
  double value, vector, inverse_vector;
  tpl = tpl_map("iA(f)A(f)A(f)", &eigensystem.length, &value, &vector, &inverse_vector);
  tpl_pack(tpl, 0);

  for (i = 0; i < eigensystem.length * eigensystem.length; ++i) {
    if (i < eigensystem.length) {
      value = eigensystem.values[i];
      tpl_pack(tpl, 1);
    }

    vector         = eigensystem.vectors[i];
    inverse_vector = eigensystem.inverse_vectors[i];

    tpl_pack(tpl, 2);
    tpl_pack(tpl, 3);
  }

  tpl_dump(tpl, TPL_FILE, parameters.filename);
  tpl_free(tpl);
}

EIGENSYSTEM deserialize_eigensystem(const POPULATION_PARAMS parameters) {
  tpl_node* tpl;
  int i, length;
  double value, vector, inverse_vector;
  EIGENSYSTEM eigensystem;

  tpl = tpl_map("iA(f)A(f)A(f)", &length, &value, &vector, &inverse_vector);
  tpl_load(tpl, TPL_FILE, parameters.filename);
  tpl_unpack(tpl, 0);

  eigensystem = init_eigensystem(length);

  for (i = 0; i < eigensystem.length * eigensystem.length; ++i) {
    if (i < eigensystem.length) {
      tpl_unpack(tpl, 1);
      eigensystem.values[i] = value;
    }

    tpl_unpack(tpl, 2);
    tpl_unpack(tpl, 3);

    eigensystem.vectors[i]         = vector;
    eigensystem.inverse_vectors[i] = inverse_vector;
  }

  tpl_free(tpl);

  return eigensystem;
}

void ensure_key_structures_and_energies_assigned(POPULATION_PARAMS* parameters) {
  int i, seq_length = strlen(parameters->sequence);

  if (parameters->start_structure == NULL) {
    parameters->start_structure = malloc((seq_length + 1) * sizeof(char));

    for (i = 0; i < seq_length; ++i) {
      parameters->start_structure[i] = '.';
    }

    parameters->start_structure[i] = '\0';
  }

  if (parameters->end_structure == NULL) {
    set_end_structure(parameters);
  }
}

void set_end_structure(POPULATION_PARAMS* parameters) {
  paramT* vienna_params;

  vienna_params             = init_vienna_params(*parameters);
  parameters->end_structure = malloc((strlen(parameters->sequence) + 1) * sizeof(char));
  fold_par(parameters->sequence, parameters->end_structure, vienna_params, 0, 0);
}

SOLUTION* sample_structures(const POPULATION_PARAMS parameters) {
  double energy_cap;
  paramT* vienna_params;

  vienna_params = init_vienna_params(parameters);
  energy_cap    = parameters.energy_cap ? (int)round(parameters.energy_cap * 100) : INF;

  return subopt_par(parameters.sequence, parameters.start_structure, vienna_params, energy_cap, 0, 0, NULL);
}

double estimate_equilibrium(const EIGENSYSTEM eigensystem, const KLP_PARAMS klp_params, const POPULATION_PARAMS parameters, int eigensystem_index) {
  int i, everything_in_equilibrium = 0;
  double equilibrium_time;

  equilibrium_time = soft_bound_for_population_proportion(eigensystem, klp_params, parameters, eigensystem_index, parameters.end_time, 1);

#ifdef INSANE_DEBUG
  printf("equilibrium_time: %f\n", equilibrium_time);
#endif

  while (equilibrium_time >= parameters.start_time && index_in_equilibrium_within_window(eigensystem, klp_params, parameters, eigensystem_index, equilibrium_time)) {
    equilibrium_time -= parameters.step_size;
  }

#ifdef INSANE_DEBUG
  printf("equilibrium_time after backtrack: %f\n", equilibrium_time);
#endif

  if (equilibrium_time < parameters.start_time) {
    return NEG_INF(parameters);
  } else {
    while (equilibrium_time <= parameters.end_time) {
      if (parameters.all_subpop_for_eq) {
        everything_in_equilibrium = 1;

        for (i = 0; i < eigensystem.length && everything_in_equilibrium; ++i) {
          everything_in_equilibrium = index_in_equilibrium_within_window(eigensystem, klp_params, parameters, i, equilibrium_time);
        }
      } else {
        everything_in_equilibrium = index_in_equilibrium_within_window(eigensystem, klp_params, parameters, eigensystem_index, equilibrium_time);
      }

      if (everything_in_equilibrium) {
        return equilibrium_time;
      } else {
        equilibrium_time += parameters.step_size;
      }
    }

    return POS_INF(parameters);
  }
}

double soft_bound_for_population_proportion(const EIGENSYSTEM eigensystem, const KLP_PARAMS klp_params, const POPULATION_PARAMS parameters, int eigensystem_index, double final_time, int sign) {
  double current_probability, final_probability, current_time, temp_time, last_time = final_time;

  if (parameters.soft_bounds) {
    last_time           = final_time;
    current_time        = parameters.start_time + (parameters.end_time - parameters.start_time) / 2.0;
    current_probability = probability_at_logtime(eigensystem, klp_params, current_time, klp_params.start_state, eigensystem_index);
    final_probability   = probability_at_logtime(eigensystem, klp_params, final_time, klp_params.start_state, eigensystem_index);

    while (fabs(current_time - last_time) > parameters.step_size) {
      current_probability = probability_at_logtime(eigensystem, klp_params, current_time, klp_params.start_state, eigensystem_index);

#ifdef INSANE_DEBUG
      printf("current_time: %+.4f\tlast_time: %+.4f\tp(%d, current): %.4f\tp(%d, target): %.4f\n", current_time, last_time, eigensystem_index, current_probability, eigensystem_index, final_probability);
#endif

      temp_time = current_time;

      if (fabs(current_probability - final_probability) > parameters.delta) {
        current_time += sign * fabs(last_time - current_time) / 2.0;
      } else {
        current_time -= sign * fabs(last_time - current_time) / 2.0;
      }

      last_time = temp_time;
    }

#ifdef INSANE_DEBUG
    printf("Bounded at %+f within %f of requested delta, %f\n\n", current_time, fabs(current_probability - final_probability), parameters.delta);
#endif

    return round(current_time / parameters.step_size) * parameters.step_size;
  } else {
    return final_time;
  }
}

int index_in_equilibrium_within_window(const EIGENSYSTEM eigensystem, const KLP_PARAMS klp_params, const POPULATION_PARAMS parameters, int eigensystem_index, double logtime) {
  int i;
  double current_proportion, future_proportion;

  current_proportion = probability_at_logtime(
                         eigensystem,
                         klp_params,
                         logtime,
                         klp_params.start_state,
                         eigensystem_index
                       );

  for (i = 1; i < parameters.window_size; ++i) {
    future_proportion = probability_at_logtime(
                          eigensystem,
                          klp_params,
                          logtime + parameters.step_size * i,
                          klp_params.start_state,
                          eigensystem_index
                        );

#ifdef INSANE_DEBUG
    printf(
      "current_time: %f, window_time: %f, index: %d, abs(p(window) - p(current)): %e, within epsilon?: %s\n",
      logtime,
      logtime + parameters.step_size * i,
      eigensystem_index,
      fabs(current_proportion - future_proportion),
      fabs(current_proportion - future_proportion) < parameters.epsilon ? "yes" : "no"
    );
#endif

    if (fabs(current_proportion - future_proportion) > parameters.epsilon) {
      return 0;
    }
  }

  return 1;
}

static int compare_tuple(const void* pa, const void* pb) {
  const double* a = *(const double**)pa;
  const double* b = *(const double**)pb;

  if (a[0] > b[0]) {
    return 1;
  } else if (a[0] < b[0]) {
    return -1;
  } else {
    return 0;
  }
}

int* array_of_indices_to_output(const EIGENSYSTEM eigensystem, const KLP_PARAMS klp_params, const POPULATION_PARAMS parameters, double logtime) {
  int i, output_size;
  int* id_list;
  double** final_proportions;

  switch (parameters.num_subpop_to_show) {
    case -1:
      output_size = 2;
      break;

    case 0:
      output_size = eigensystem.length;
      break;

    default:
      output_size = parameters.num_subpop_to_show;
  }

#ifdef INSANE_DEBUG
  printf("number of subpopulations to output: %d\n", output_size);
#endif

  id_list    = malloc((output_size + 1) * sizeof(int));
  id_list[0] = output_size;

  if (parameters.num_subpop_to_show == -1) {
    id_list[1] = klp_params.end_state;
    id_list[2] = klp_params.start_state;
  } else {
    final_proportions = malloc(eigensystem.length * sizeof(double*));

    for (i = 0; i < eigensystem.length; ++i) {
      final_proportions[i] = malloc(2 * sizeof(double));
    }

    for (i = 0; i < eigensystem.length; ++i) {
      final_proportions[i][0] = probability_at_logtime(eigensystem, klp_params, logtime, klp_params.start_state, i);
      final_proportions[i][1] = i;
    }

#ifdef INSANE_DEBUG
    printf("\nBefore calling qsort:\n\n");

    for (i = 0; i < eigensystem.length; ++i) {
      printf("final_proportions[%d]\t%.8f\n", (int)final_proportions[i][1], final_proportions[i][0]);
    }

#endif

    qsort(final_proportions, eigensystem.length, sizeof(final_proportions[0]), compare_tuple);

#ifdef INSANE_DEBUG
    printf("\nAfter calling qsort:\n\n");

    for (i = 0; i < eigensystem.length; ++i) {
      printf("final_proportions[%d]\t%.8f\n", (int)final_proportions[i][1], final_proportions[i][0]);
    }

    printf("\nIndices to keep:\n\n");
#endif

    for (i = 1; i <= id_list[0]; ++i) {
      id_list[i] = (int)final_proportions[eigensystem.length - i][1];

#ifdef INSANE_DEBUG
      printf("final_proportions[%d]\t%.8f\n", id_list[i], final_proportions[eigensystem.length - i][0]);
#endif
    }
  }

  return id_list;
}

void print_equilibrium(const EIGENSYSTEM eigensystem, const KLP_PARAMS klp_params, const POPULATION_PARAMS parameters) {
  double soft_right_bound, equilibrium_time;
  int i, list_length;
  int* id_list;

  soft_right_bound = soft_bound_for_population_proportion(eigensystem, klp_params, parameters, klp_params.end_state, parameters.end_time, 1);
  id_list          = array_of_indices_to_output(eigensystem, klp_params, parameters, soft_right_bound);
  list_length      = id_list++[0];

  printf("index\tlogtime\n");

  for (i = 0; i < list_length; ++i) {
    equilibrium_time = estimate_equilibrium(eigensystem, klp_params, parameters, id_list[i]);

    printf("%d\t", id_list[i]);

    if (equilibrium_time == NEG_INF(parameters)) {
      printf("-Infinity\n");
    } else if (equilibrium_time == POS_INF(parameters)) {
      printf("Infinity\n");
    } else {
      printf("%+.8f\n", equilibrium_time);
    }
  }
}

void print_population_proportion(const EIGENSYSTEM eigensystem, const KLP_PARAMS klp_params, const POPULATION_PARAMS parameters) {
  double current_time, soft_left_bound, soft_right_bound;
  int i, list_length;
  int* id_list;

  soft_left_bound  = soft_bound_for_population_proportion(eigensystem, klp_params, parameters, klp_params.start_state, parameters.start_time, -1);
  soft_right_bound = soft_bound_for_population_proportion(eigensystem, klp_params, parameters, klp_params.end_state, parameters.end_time, 1);
  id_list          = array_of_indices_to_output(eigensystem, klp_params, parameters, soft_right_bound);
  list_length      = id_list++[0];

  printf("logtime\t");

  for (i = 0; i < list_length; ++i) {
    printf("%d\t", id_list[i]);
  }

  printf("\n");

  for (current_time = soft_left_bound; current_time <= soft_right_bound; current_time += parameters.step_size) {
    printf("%+.8f\t", current_time);

    for (i = 0; i < list_length - 1; ++i) {
      printf("%+.8f\t", probability_at_logtime(eigensystem, klp_params, current_time, klp_params.start_state, id_list[i]));
    }

    printf("%+.8f", probability_at_logtime(eigensystem, klp_params, current_time, klp_params.start_state, id_list[i]));

    printf("\n");
  }
}

void print_array(char* title, double* matrix, int length) {
  int i;
  printf("%s\n", title);

  for (i = 0; i < length; ++i) {
    printf("%+.4f\n", matrix[i]);
  }

  printf("\n");
}

void print_square_matrix(char* title, double* matrix, int length) {
  int i, j;
  printf("%s\n", title);

  for (i = 0; i < length; ++i) {
    for (j = 0; j < length; ++j) {
      printf("%+.2e ", E_ROW_ORDER(matrix, i, j, length));
    }

    printf("\n");
  }

  printf("\n");
}

void print_eigenvalues(const EIGENSYSTEM eigensystem) {
  int i;

  for (i = 0; i < eigensystem.length; ++i) {
    printf("%+.8f\n", eigensystem.values[i]);
  }
}

max_stars_repo_name: evansenter/hermes

max_issues_repo_name: evansenter/hermes

max_forks_repo_name: evansenter/hermes


text:
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <assert.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_fit.h>
#include <gsl/gsl_roots.h>

#include "cosmocalc.h"

double fourierTransformTopHat(double y);
double tophatradnorm_linear_powspec_exact_nonorm_lnk_integ_funct_I0(double lnk, void *p);
double tophatradnorm_linear_powspec_exact_nonorm_k_integ_funct_I0(double k, void *p);

double convert_cmbnorm2sigma8_k_integ_funct_I0(double k, void *p)
{
  //k is in units of 1/Mpc
  //to get to units of h/Mpc
  // k -> k/h
  double ft = fourierTransformTopHat(k*8.0/cosmoData.h);
  double tf = transfer_function(k/cosmoData.h);

  /*
  fprintf(stderr,"Tf^2 = %g, W(k)^2 = %g (2/5/OmegaM)^2 = %g (k/H0)^4/k = %g, nspow = %g\n",
	  tf*tf,ft*ft,4.0/25.0/cosmoData.OmegaM/cosmoData.OmegaM,
	  pow(k/100.0,4.0)/k,pow(k*cosmoData.h/cosmoData.As_pivot,cosmoData.SpectralIndex-1.0)
	  );
  */
  return cosmoData.As*4.0/25.0/cosmoData.OmegaM/cosmoData.OmegaM
    *pow(k/cosmoData.h*DH,4.0)/k*pow(k/cosmoData.As_pivot,cosmoData.SpectralIndex-1.0)
    *ft*ft*tf*tf;
}

double convert_cmbnorm2sigma8(void)
{
  double I0,I1;
  double abserr;
  double epsrel,epsabs;
  gsl_integration_workspace *workspace;
  gsl_function F;

#define WORKSPACE_NUM 10000000
#define ABSERR 1e-8
#define RELERR 0.0
  workspace = gsl_integration_workspace_alloc((size_t) WORKSPACE_NUM);

  epsabs = 1e-20;
  epsrel = 1e-6;
  F.function = &convert_cmbnorm2sigma8_k_integ_funct_I0;
  gsl_integration_qags(&F,0.0,2.0*M_PI/8.0,epsabs,epsrel,(size_t) WORKSPACE_NUM,workspace,&I0,&abserr);
  gsl_integration_qagiu(&F,2.0*M_PI/8.0,epsabs,epsrel,(size_t) WORKSPACE_NUM,workspace,&I1,&abserr);

  gsl_integration_workspace_free(workspace);
#undef ABSERR
#undef RELERR
#undef WORKSPACE_NUM

  //fprintf(stderr,"I0 = %g, I1 = %g, g = %g\n",I0,I1,growth_function_exact_nonorm(1.0));

  return sqrt(I0 + I1)*growth_function_exact_nonorm(1.0);
}

double fourierTransformTopHat(double y)
{
  if(y < 1e-3)
    return 1.0;
  else
    return 3.0/y/y/y*(sin(y) - y*cos(y));
}

double tophatradnorm_linear_powspec_exact_nonorm_lnk_integ_funct_I0(double lnk, void *p)
{
  // topHatRad = (*((double*)p));
  double k = exp(lnk);
  double ft = fourierTransformTopHat(k*(*((double*)p)));
  double tf = transfer_function(k);
  return ft*ft*tf*tf*pow(k,cosmoData.SpectralIndex)*k*k/2.0/M_PI/M_PI*k;
}

double tophatradnorm_linear_powspec_exact_nonorm_k_integ_funct_I0(double k, void *p)
{
  // topHatRad = (*((double*)p));
  double ft = fourierTransformTopHat(k*(*((double*)p)));
  double tf = transfer_function(k);
  return ft*ft*tf*tf*pow(k,cosmoData.SpectralIndex)*k*k/2.0/M_PI/M_PI;
}

double tophatradnorm_linear_powspec_exact_nonorm(double topHatRad)
{
  double I0,I1;
  double abserr;
  double epsrel,epsabs;
  gsl_integration_workspace *workspace;
  gsl_function F;

#define WORKSPACE_NUM 10000000
#define ABSERR 1e-8
#define RELERR 0.0
  workspace = gsl_integration_workspace_alloc((size_t) WORKSPACE_NUM);

  F.params = &(topHatRad);
  if(topHatRad > 1e-4)
    {
      epsabs = 1e-20;
      epsrel = 1e-6;
      F.function = &tophatradnorm_linear_powspec_exact_nonorm_k_integ_funct_I0;
      gsl_integration_qags(&F,0.0,2.0*M_PI/topHatRad,epsabs,epsrel,(size_t) WORKSPACE_NUM,workspace,&I0,&abserr);
      gsl_integration_qagiu(&F,2.0*M_PI/topHatRad,epsabs,epsrel,(size_t) WORKSPACE_NUM,workspace,&I1,&abserr);
    }
  else
    {
      F.function = &tophatradnorm_linear_powspec_exact_nonorm_lnk_integ_funct_I0;
      gsl_integration_qagil(&F,log(2.0*M_PI/topHatRad),ABSERR,RELERR,(size_t) WORKSPACE_NUM,workspace,&I0,&abserr);
      gsl_integration_qagiu(&F,log(2.0*M_PI/topHatRad),ABSERR,RELERR,(size_t) WORKSPACE_NUM,workspace,&I1,&abserr);
    }

  gsl_integration_workspace_free(workspace);
#undef ABSERR
#undef RELERR
#undef WORKSPACE_NUM

  return I0 + I1;
}

double linear_powspec_exact(double k, double a)
{
  static int initFlag = 1;
  static int currCosmoNum;
  static double linear_powspec_norm = 1.0;
  double gf = growth_function(a);
  double tf = transfer_function(k);

  if(initFlag == 1 || currCosmoNum != cosmoData.cosmoNum)
    {
      initFlag = 0;
      currCosmoNum = cosmoData.cosmoNum;

      linear_powspec_norm = cosmoData.Sigma8*cosmoData.Sigma8/tophatradnorm_linear_powspec_exact_nonorm(8.0);
    }

  return tf*tf*pow(k,cosmoData.SpectralIndex)*gf*gf*linear_powspec_norm;
}

double linear_powspec(double k, double a)
{
  static int initFlag = 1;
  static int currCosmoNum;
  static double linear_powspec_norm = 1.0;
  static gsl_spline *cosmocalc_linear_powspec_spline = NULL;
  static gsl_interp_accel *cosmocalc_linear_powspec_acc = NULL;
  static double c0,c1;

  double linear_powspec_table[COSMOCALC_LINEAR_POWSPEC_TABLE_LENGTH];
  double k_table[COSMOCALC_LINEAR_POWSPEC_TABLE_LENGTH];
  long i;
  double gf,cov00,cov01,cov11,sumsq,tf;

  if(initFlag == 1 || currCosmoNum != cosmoData.cosmoNum)
    {
      initFlag = 0;
      currCosmoNum = cosmoData.cosmoNum;

      for(i=0;i<COSMOCALC_LINEAR_POWSPEC_TABLE_LENGTH;++i)
	{
	  k_table[i] = log(K_MAX/K_MIN)/(COSMOCALC_LINEAR_POWSPEC_TABLE_LENGTH-1.0)*((double) i) + log(K_MIN);
	  tf = transfer_function(exp(k_table[i]));
	  linear_powspec_table[i] = log(tf*tf*pow(exp(k_table[i]),cosmoData.SpectralIndex));
	}

      //init the spline and accelerators
      if(cosmocalc_linear_powspec_spline != NULL)
	gsl_spline_free(cosmocalc_linear_powspec_spline);
      cosmocalc_linear_powspec_spline = gsl_spline_alloc(GSL_SPLINE_TYPE,(size_t) (COSMOCALC_LINEAR_POWSPEC_TABLE_LENGTH));
      gsl_spline_init(cosmocalc_linear_powspec_spline,k_table,linear_powspec_table,(size_t) (COSMOCALC_LINEAR_POWSPEC_TABLE_LENGTH));
      if(cosmocalc_linear_powspec_acc != NULL)
	gsl_interp_accel_reset(cosmocalc_linear_powspec_acc);
      else
	cosmocalc_linear_powspec_acc = gsl_interp_accel_alloc();

      gsl_fit_linear(k_table+COSMOCALC_LINEAR_POWSPEC_TABLE_LENGTH-COSMOCALC_LINEAR_POWSPEC_FIT_LENGTH,(size_t) 1,
		     linear_powspec_table+COSMOCALC_LINEAR_POWSPEC_TABLE_LENGTH-COSMOCALC_LINEAR_POWSPEC_FIT_LENGTH,(size_t) 1,
		     (size_t) COSMOCALC_LINEAR_POWSPEC_FIT_LENGTH,&c0,&c1,&cov00,&cov01,&cov11,&sumsq);

      linear_powspec_norm = cosmoData.Sigma8*cosmoData.Sigma8/tophatradnorm_linear_powspec_exact_nonorm(8.0);
    }

  gf = growth_function(a);

  if(k < K_MIN)
    {
      tf = transfer_function(k);
      return tf*tf*pow(k,cosmoData.SpectralIndex)*linear_powspec_norm*gf*gf;
    }
  else if(k < K_MAX)
    return exp(gsl_spline_eval(cosmocalc_linear_powspec_spline,log(k),cosmocalc_linear_powspec_acc))*linear_powspec_norm*gf*gf;
  else
    return exp(c0+c1*log(k))*linear_powspec_norm*gf*gf;
}

double tophatnorm_linear_powspec(double topHatRad)
{
#define NL_RTOPHAT_MIN 0.0001
#define NL_RTOPHAT_MAX 100.0
#define WORKSPACE_NUM 10000000
#define ABSERR 1e-8
#define RELERR 0.0

  static int initFlag = 1;
  static int currCosmoNum;
  static gsl_spline *spline = NULL;
  static gsl_interp_accel *accel = NULL;
  static double linear_powspec_norm;

  double I0,I1;
  double abserr,epsabs,epsrel;
  gsl_integration_workspace *workspace;
  gsl_function F;
  int i;
  double xtab[COSMOCALC_LINEAR_POWSPEC_NORM_TABLE_LENGTH];
  double ytab[COSMOCALC_LINEAR_POWSPEC_NORM_TABLE_LENGTH];
  double dlnr;
  double lnrmin;
  double lnr;

  if(initFlag == 1 || currCosmoNum != cosmoData.cosmoNum)
    {
      initFlag = 0;
      currCosmoNum = cosmoData.cosmoNum;

      workspace = gsl_integration_workspace_alloc((size_t) WORKSPACE_NUM);

      dlnr = log(NL_RTOPHAT_MAX/NL_RTOPHAT_MIN)/(COSMOCALC_LINEAR_POWSPEC_NORM_TABLE_LENGTH-1.0);
      lnrmin = log(NL_RTOPHAT_MIN);

      for(i=0;i<COSMOCALC_LINEAR_POWSPEC_NORM_TABLE_LENGTH;++i)
        {
          lnr = dlnr*i + lnrmin;
          xtab[i] = exp(lnr);

	  F.params = &(xtab[i]);
	  if(topHatRad > 1e-4)
	    {
	      epsabs = 1e-20;
	      epsrel = 1e-6;
	      F.function = &tophatradnorm_linear_powspec_exact_nonorm_k_integ_funct_I0;
	      gsl_integration_qags(&F,0.0,2.0*M_PI/topHatRad,epsabs,epsrel,(size_t) WORKSPACE_NUM,workspace,&I0,&abserr);
	      gsl_integration_qagiu(&F,2.0*M_PI/topHatRad,epsabs,epsrel,(size_t) WORKSPACE_NUM,workspace,&I1,&abserr);
	    }
	  else
	    {
	      F.function = &tophatradnorm_linear_powspec_exact_nonorm_lnk_integ_funct_I0;
	      gsl_integration_qagil(&F,log(2.0*M_PI/topHatRad),ABSERR,RELERR,(size_t) WORKSPACE_NUM,workspace,&I0,&abserr);
	      gsl_integration_qagiu(&F,log(2.0*M_PI/topHatRad),ABSERR,RELERR,(size_t) WORKSPACE_NUM,workspace,&I1,&abserr);
	    }

          xtab[i] = lnr;
          ytab[i] = log(I0+I1);
        }

      gsl_integration_workspace_free(workspace);

      if(spline != NULL)
        gsl_spline_free(spline);
      spline = gsl_spline_alloc(GSL_SPLINE_TYPE,(size_t) (COSMOCALC_LINEAR_POWSPEC_NORM_TABLE_LENGTH));
      gsl_spline_init(spline,xtab,ytab,(size_t) (COSMOCALC_LINEAR_POWSPEC_NORM_TABLE_LENGTH));
      if(accel != NULL)
        gsl_interp_accel_reset(accel);
      else
        accel = gsl_interp_accel_alloc();

      linear_powspec_norm = cosmoData.Sigma8*cosmoData.Sigma8/exp(gsl_spline_eval(spline,log(8.0),accel));

#undef ABSERR
#undef RELERR
#undef WORKSPACE_NUM
#undef NL_RTOPHAT_MIN
#undef NL_RTOPHAT_MAX
    }

  return exp(gsl_spline_eval(spline,log(topHatRad),accel))*linear_powspec_norm;
}

double linear_tophatnorm_scale_funct(double rad, void *p)
{
  double gf = ((double*)p)[0];

  return tophatnorm_linear_powspec(rad)*gf*gf-1.0;
}

double get_linear_tophatnorm_scale(double a)
{
  double gf = growth_function(a);
  double Rsigma,Rlow=0.001,Rhigh=10.0;
  int itr,maxItr=1000,status;

  //fprintf(stderr,"sigma2(Rlow) = %f, sigma2(Rhigh) = %f\n",tophatnorm_linear_powspec(Rlow)*gf*gf,tophatnorm_linear_powspec(Rhigh)*gf*gf);

#define ABSERR 1e-6
#define RELERR 1e-6
  const gsl_root_fsolver_type *T;
  gsl_root_fsolver *s;
  gsl_function F;

  F.function = &linear_tophatnorm_scale_funct;
  F.params = &gf;

  T = gsl_root_fsolver_brent;
  s = gsl_root_fsolver_alloc(T);
  gsl_root_fsolver_set(s,&F,Rlow,Rhigh);
  itr = 0;

  do
    {
      itr++;
      status = gsl_root_fsolver_iterate(s);
      Rsigma = gsl_root_fsolver_root(s);
      Rlow = gsl_root_fsolver_x_lower(s);
      Rhigh = gsl_root_fsolver_x_upper(s);
      status = gsl_root_test_interval(Rlow,Rhigh,ABSERR,RELERR);
    }
  while(status == GSL_CONTINUE && itr < maxItr);

#undef ABSERR
#undef RELERR

  gsl_root_fsolver_free(s);

  return Rsigma;
}

max_stars_repo_name: beckermr/cosmocalc

max_issues_repo_name: beckermr/cosmocalc

max_forks_repo_name: beckermr/cosmocalc


text:
#include <stdio.h>
#include <gsl/gsl_combination.h>

int
main (void)
{
  gsl_combination * c;
  size_t i;

  printf ("All subsets of {0,1,2,3} by size:\n") ;
  for (i = 0; i <= 4; i++)
    {
      c = gsl_combination_calloc (4, i);
      do
        {
          printf ("{");
          gsl_combination_fprintf (stdout, c, " %u");
          printf (" }\n");
        }
      while (gsl_combination_next (c) == GSL_SUCCESS);
      gsl_combination_free (c);
    }

  return 0;
}

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/*************************************************************
  Generalized-ICP Copyright (c) 2009 Aleksandr Segal.
  All rights reserved.

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

* Redistributions of source code must retain the above
  copyright notice, this list of conditions and the
  following disclaimer.
* 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.
* The names of the contributors may not be used to endorse
  or promote products derived from this software
  without specific prior written permission.

  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
  CONTRIBUTORS "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 THE
  COPYRIGHT OWNER 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.
*************************************************************/

#ifndef GICP_H_
#define GICP_H_

#include <ANN.h>
#include <vector>
#include <iostream>
//#include <gsl/gsl.h>
#include <semaphore.h>
#include "transform.h"

namespace dgc {
namespace gicp {
typedef double gicp_mat_t[3][3];

struct GICPPoint {
    double x, y, z;
    float range;
    gicp_mat_t C; // covariance matrix
};

class GICPPointSet {
public:
    GICPPointSet();
    ~GICPPointSet();
    void BuildKDTree();
    void ComputeMatrices();

    void SavePoints(const char *filename);
    void SaveMatrices(const char *filename);

    int NumPoints() { return (int)point_.size(); }
    void Clear(void);
    int Size() { return point_.size(); }
    inline void AppendPoint(GICPPoint const & pt) { point_.push_back(pt); }
    void SetMaxIteration(int iter) { max_iteration_ = iter; }
    void SetMaxIterationInner(int iter) { max_iteration_inner_ = iter; }
    void SetEpsilon(double eps) { epsilon_ = eps; }
    void SetSolveRotation(bool s) { solve_rotation_ = s; }
    void SetGICPEpsilon(double eps) { gicp_epsilon_ = eps; }
    void SetDebug(bool d) { debug_ = d; }

    GICPPoint & operator[](int i) { return point_[i]; }
    GICPPoint const& operator[](int i) const { return point_[i]; }

    // returns number of iterations it took to converge
    int AlignScan(GICPPointSet *scan, dgc_transform_t base_t, dgc_transform_t t, double max_match_dist, bool save_error_plot = 0);

private:
    std::vector <GICPPoint> point_;
    ANNpointArray kdtree_points_;
    ANNkd_tree *kdtree_;
    int max_iteration_;
    int max_iteration_inner_;
    double epsilon_;
    double epsilon_rot_;
    double gicp_epsilon_;
    bool debug_;
    bool solve_rotation_;
    bool matrices_done_;
    bool kdtree_done_;
    pthread_mutex_t mutex_;
};

}

}

#endif

max_stars_repo_name: raandoom/gicp

max_issues_repo_name: raandoom/gicp

max_forks_repo_name: raandoom/gicp


text:
#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

void mcmc
(
 int &mcID,
 int &seed,
 int &imh,
 int &vmh,
 int &neighborPast,
 int &n_neighbor,
 int &ptAccpt,
 double &pllh,
 double *jump,
 struct FXP &fxp,
 struct PAR &par_a,
 struct PAR &par_c,
 struct RND &rnd,
 struct VAR &var,
 struct DATA &data,
 struct DIAG &diag,
 struct HIST &hist,
 struct INCV &incv,
 struct POST &post,
 struct PRED &pred,
 gsl_rng *rng,
 FILE *outp
);

void logQuasiposterior
(
 //=======================================================
 double &llh,
 struct FXP &fxp,
 struct PAR &par,
 struct RND &rnd,
 struct VAR &var,
 struct DATA &data,
 struct DIAG &diag,
 struct INCV &incv,
 struct PRED &pred
//=======================================================
);

void logQuasiposterior
(
 //=======================================================
 double &llh,
 double *expDeltaIn,
 struct PAR &par,
 struct VAR &var
//=======================================================
);

void rollbackTest
(
 int &imh,
 int &vmh,
 int &n_neighbor,
 int &neighborPast,
 double &pllh,
 double &pllhbackup,
 struct FXP &fxp,
 struct PAR &par,
 struct VAR &var,
 struct DATA &data,
 struct DIAG &diag,
 struct HIST &hist,
 struct INCV &incv,
 struct PRED &pred
);