text:
/*
** read data matrix, map of voxel addresses
**
** G.Lohmann, Jan 2013
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_sort.h>

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

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

VImage *VImagePointer(VAttrList list,int *nt)
{
  VImage tmp;
  VAttrListPosn posn;
  int i,ntimesteps,nrows,ncols,nslices;

  /* get image dimensions, read functional data */
  nslices = 0;
  for (VFirstAttr (list, & posn); VAttrExists (& posn); VNextAttr (& posn)) {
    if (VGetAttrRepn (& posn) != VImageRepn) continue;
    VGetAttrValue (& posn, NULL,VImageRepn, & tmp);
    if (VPixelRepn(tmp) != VShortRepn) continue;
    nslices++;
  }
  /* VDestroyImage(tmp); */
  if (nslices < 1) VError(" no slices");

  VImage *src = (VImage *) VCalloc(nslices,sizeof(VImage));
  i = ntimesteps = nrows = ncols = 0;
  for (VFirstAttr (list, & posn); VAttrExists (& posn); VNextAttr (& posn)) {
    if (VGetAttrRepn (& posn) != VImageRepn) continue;
    VGetAttrValue (& posn, NULL,VImageRepn, & src[i]);
    if (VPixelRepn(src[i]) != VShortRepn) continue;
    if (VImageNBands(src[i]) > ntimesteps) ntimesteps = VImageNBands(src[i]);
    if (VImageNRows(src[i])  > nrows)  nrows = VImageNRows(src[i]);
    if (VImageNColumns(src[i]) > ncols) ncols = VImageNColumns(src[i]);
    i++;
  }
  nslices = i;
  if (nslices < 1) return NULL;
  if (ntimesteps < 2) return NULL;
  *nt = ntimesteps;
  return src;
}

/* read functional data */
void VReadImagePointer(VAttrList list,VImage *src)
{
  VAttrListPosn posn;

  int i = 0;
  for (VFirstAttr (list, & posn); VAttrExists (& posn); VNextAttr (& posn)) {
    if (VGetAttrRepn (& posn) != VImageRepn) continue;
    VGetAttrValue (& posn, NULL,VImageRepn, & src[i]);
    if (VPixelRepn(src[i]) != VShortRepn) continue;
    i++;
  }
}

/* check if mask covers data. If not, set surplus voxels to zero */
long VMaskCoverage(VAttrList list,VImage mask)
{
  VImage src=NULL;
  VAttrListPosn posn;
  int b,r,c;
  long count=0;

  b = 0;
  for (VFirstAttr (list, & posn); VAttrExists (& posn); VNextAttr (& posn)) {
    if (VGetAttrRepn (& posn) != VImageRepn) continue;
    VGetAttrValue (& posn, NULL,VImageRepn, & src);
    if (VPixelRepn(src) != VShortRepn) continue;
    if (VImageNRows(src) < 3) continue;

    if (VImageNRows(src) != VImageNRows(mask))
      VError(" inconsistent nrows: %d,  mask: %d",VImageNRows(src),VImageNRows(mask));
    if (VImageNColumns(src) != VImageNColumns(mask))
      VError(" inconsistent ncols: %d,  mask: %d",VImageNColumns(src),VImageNColumns(mask));

    for (r=0; r<VImageNRows(src); r++) {
      for (c=0; c<VImageNColumns(src); c++) {
	if (b >= VImageNBands(mask)) VError("VMaskCoverage, illegal addr, band= %d",b);
	float u = VGetPixel(mask,b,r,c);
	int j = (int)VPixel(src,0,r,c,VShort);
	if (j == 0 && u > 0.3) {
	  VSetPixel(mask,b,r,c,0);
	  count++;
	}
      }
    }
    b++;
  }
  return count;
}

/* check if data matrix contains zero voxels */
void VCheckMatrix(gsl_matrix_float *X)
{
  long i,j,nvox=X->size1,nt=X->size2;
  double sum1,sum2,nx,mean,var,tiny=1.0e-4;

  long count = 0;
  nx = (double)nt;
  for (i=0; i<nvox; i++) {
    const float *arr1 = gsl_matrix_float_const_ptr(X,i,0);
    sum1 = sum2 = 0;
    for (j=0; j<nt; j++) {
      const double u = (double)arr1[j];
      sum1 += u;
      sum2 += u*u;
    }
    mean = sum1/nx;
    var = (sum2 - nx * mean * mean) / (nx - 1.0);
    if (var < tiny) count++;
  }
  if (count > 0)
    VWarning(" number of empty voxels: %ld",count);
}

VImage VoxelMap(VImage mask,size_t *nvoxels)
{
  int b,r,c;
  int nslices = VImageNBands(mask);
  int nrows = VImageNRows(mask);
  int ncols = VImageNColumns(mask);

  /* count number of non-zero voxels */
  size_t nvox = 0;
  for (b=0; b<nslices; b++) {
    for (r=0; r<nrows; r++) {
      for (c=0; c<ncols; c++) {
	float u = VGetPixel(mask,b,r,c);
	if (ABS(u) < 0.5) continue;
	nvox++;
      }
    }
  }
  (*nvoxels) = nvox;

  /* voxel addresses */
  VImage map = VCreateImage(1,4,nvox,VShortRepn);
  if (map == NULL) VError(" error allocating addr map");
  VFillImage(map,VAllBands,0);
  VImage xmap=NULL;
  VExtractAttr (VImageAttrList(mask),"map",NULL,VImageRepn,&xmap,FALSE);
  VCopyImageAttrs (mask,map);
  VSetAttr(VImageAttrList(map),"nvoxels",NULL,VLongRepn,(VLong)nvox);
  VSetAttr(VImageAttrList(map),"nslices",NULL,VLongRepn,(VLong)nslices);
  VSetAttr(VImageAttrList(map),"nrows",NULL,VLongRepn,(VLong)nrows);
  VSetAttr(VImageAttrList(map),"ncols",NULL,VLongRepn,(VLong)ncols);

  VPixel(map,0,3,0,VShort) = nslices;
  VPixel(map,0,3,1,VShort) = nrows;
  VPixel(map,0,3,2,VShort) = ncols;

  size_t i = 0;
  for (b=0; b<nslices; b++) {
    for (r=0; r<nrows; r++) {
      for (c=0; c<ncols; c++) {
	float u = VGetPixel(mask,b,r,c);
	if (ABS(u) < 0.5) continue;
	VPixel(map,0,0,i,VShort) = b;
	VPixel(map,0,1,i,VShort) = r;
	VPixel(map,0,2,i,VShort) = c;
	i++;
      }
    }
  }
  return map;
}

/* two-pass formula, correct for round-off error */
void VNormalize(float *data,int nt,VBoolean stddev)
{
  int j;
  float ave,var,sd,nx,s,u,tiny=1.0e-6;

  nx = (float)nt;
  ave = 0;
  for (j=0; j<nt; j++) ave += data[j];
  ave /= nx;

  var = u = 0;
  for (j=0; j<nt; j++) {
    s = data[j]-ave;
    u   += s;
    var += s*s;
  }
  var=(var-u*u/nx)/(nx-1);
  sd = sqrt(var);
  if (sd < tiny) {
    for (j=0; j<nt; j++) data[j] = 0;
  }
  else {
    if (stddev==FALSE) sd = 1.0;  /* only subtract mean !!!! */
    for (j=0; j<nt; j++) {
      u = data[j];
      data[j] = (u-ave)/sd;
    }
  }
}

/*  copy data to input matrix X, contains fMRI time series */
void VDataMatrix(VImage *src,int first,int len,VImage map,gsl_matrix_float *X)
{
  long i,j,k,nvox;
  int b,r,c;
  float *data=NULL;

  data = (float *) VCalloc(len,sizeof(float));

  gsl_matrix_float_set_zero (X);
  nvox = VImageNColumns(map);
  if (nvox != (long)X->size1) VError(" err, %ld %ld",(long)X->size1,nvox);
  if (len  != (long)X->size2) VError(" err, %ld %ld",(long)X->size2,len);

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

    b = VPixel(map,0,0,i,VShort);
    r = VPixel(map,0,1,i,VShort);
    c = VPixel(map,0,2,i,VShort);

    if ((first+len) > VImageNBands(src[b])) VError(" illegal len addr, %d  %d %d %d",
						   b,first,len,VImageNBands(src[b]));
    if (r >= VImageNRows(src[b])) VError(" illegal row addr");
    if (c >= VImageNColumns(src[b])) VError(" illegal column addr");

    k = 0;
    for (j=first; j<first+len; j++) {
      data[k] = (float) VGetPixel(src[b],j,r,c);
      k++;
    }
    VNormalize(data,len,TRUE);

    float *ptr = gsl_matrix_float_ptr(X,i,0);
    for (j=0; j<len; j++) *ptr++ = data[j];
  }
  VFree(data);
}

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
#pragma once

/***************************************************
 *************** Auto-Generated File ***************
 ***************************************************/

#include <cstring>
#include <string>

#include <gsl/gsl_const_mksa.h>
#include <gsl/gsl_const_num.h>

#include "../../Utils/BinarySearch.h"

enum UnitType {
	NONE,
	{UnitTypes}
};

constexpr const int numUnits = {numUnits};
constexpr const char* unitNames[numUnits] = {
	{unitNames}
};
constexpr const int unitIndices[numUnits] = {
	{unitIndices}
};
constexpr const char * unitAbbreviations[numUnits] = {
	{unitAbbreviations}
};
constexpr const UnitType unitTypes[numUnits] = {
	{unitTypes}
};
constexpr const double unitConversions[numUnits] = {
	{unitConversions}
};

/*
Returns index of the unit in the unitNames array.

Uses binary search under the hood to search for the index.

Parameters
----------
name: The name of the unit

Returns
-------
The index or -1 if the provided name is not a unit.
*/
CONSTEXPR_BINARY_SEARCH(getUnitIndex, unitNames, numUnits)

CONSTEXPR_BINARY_SEARCH(getAbbrIndex, unitAbbreviations, numUnits)

constexpr UnitType getUnitType(const char* name){
    int index = getUnitIndex(name);
    if (index != -1){
		return NONE;
	}
	return unitTypes[index];
}
inline UnitType getUnitType(const std::string& name){
    int index = getUnitIndex(name.c_str());
    if (index != -1){
		return NONE;
	}
	return unitTypes[index];
}
inline UnitType getUnitType(int index){
    if (index < 0 || index >= numUnits){
		return NONE;
	}
	return unitTypes[index];
}

constexpr double getUnitConversion(const char* name){
    int index = getUnitIndex(name);
    if (index != -1){
		return NONE;
	}
	return unitConversions[index];
}
inline double getUnitConversion(const std::string& name){
    int index = getUnitIndex(name.c_str());
    if (index != -1){
		return NONE;
	}
	return unitConversions[index];
}
inline double getUnitConversion(int index){
    if (index < 0 || index >= numUnits){
		return NONE;
	}
	return unitConversions[index];
}

ext: h

max_stars_repo_name: antoniojkim/CalcPlusPlus

max_issues_repo_name: antoniojkim/CalcPlusPlus

max_forks_repo_name: antoniojkim/CalcPlusPlus


text:
// Test driver for the pyramid class.

#include <math.h>
#include <stdlib.h>

#include <glib.h>
#include <gsl/gsl_math.h>

#include <asf_meta.h>
#include <float_image.h>

#include "meta_read_wrapper.h"
#include "pyramid.h"

// Base name of .img/.meta pair we will use for testing.
//#define TEST_IMG_BASENAME "E116306133G1S009"
#define TEST_IMG_BASENAME "E116306133G1S009_subimage"

// Directory to use as a scratch directory for pyramid layers.
#define SCRATCH_DIR "/tmp/"

// When we make JPEG images we give this mask value to the routine
// which does it.  This make should have any effect for the image we
// are using.
#define TEST_IMG_BOGUS_MASK_VALUE -999999.9

int
main (void)
{
  g_print ("This test may take a few moments...\n");

  GString *meta_name = g_string_new (TEST_IMG_BASENAME);
  g_string_append (meta_name, ".meta");

  GString *data_name = g_string_new (TEST_IMG_BASENAME);
  g_string_append (data_name, ".img");

  meta_parameters *md = meta_read_wrapper (meta_name->str);

  g_string_free (meta_name, TRUE);

  ssize_t tisx = md->general->sample_count;
  ssize_t tisy = md->general->line_count;

  FloatImage *tiafi
    = float_image_new_from_file (tisx, tisy, data_name->str, 0,
				 FLOAT_IMAGE_BYTE_ORDER_BIG_ENDIAN);

  g_string_free (data_name, TRUE);

  float_image_export_as_jpeg (tiafi, "test_pyramid_base_layer_direct.jpg",
			      GSL_MAX (tisx, tisy), TEST_IMG_BOGUS_MASK_VALUE);

  meta_free (md);

  Pyramid *p = pyramid_new (TEST_IMG_BASENAME, SCRATCH_DIR);

  // "Testing" currently consists of creating jpeg images of each
  // layer and checking to see if they look decent.

  size_t ii;
  for ( ii = 0 ; ii < p->layers->len ; ii++ ) {

    pyramid_layer *cl = g_ptr_array_index (p->layers, ii);

    float *region;
    size_t rstart_x, rstart_y, rw, rh;
    gboolean unowned_memory;

    pyramid_get_region (p, 0, 0, tisx, tisy, pow (2, ii), &region, &rstart_x,
			&rstart_y, &rw, &rh, &unowned_memory);
    g_assert (rstart_x == 0 & rstart_y == 0);
    g_assert (rstart_x + rw * pow (2.0, ii) >= tisx);
    g_assert (rstart_y + rh * pow (2.0, ii) >= tisy);

    FloatImage *layer_image = float_image_new_from_memory (rw, rh, region);

    if ( unowned_memory ) {
      g_free (region);
    }

    GString *layer_file_name = g_string_new ("");
    g_string_append_printf (layer_file_name, "layer_%d.jpg", ii);

    ssize_t max_dimension = rw > rh ? rw : rh;

    int return_code
      = float_image_export_as_jpeg (layer_image, layer_file_name->str,
				    max_dimension, TEST_IMG_BOGUS_MASK_VALUE);
    g_assert (return_code == 0);

    float_image_unref (layer_image);

    g_string_free (layer_file_name, TRUE);
  }

  g_print ("JPEG images of pyramid layers formed, now look at them and\n"
	   "'touch test_pyramid_stamp' if they are ok.\n");

  exit (EXIT_SUCCESS);
}

ext: c

max_stars_repo_name: glshort/MapReady

max_issues_repo_name: glshort/MapReady

max_forks_repo_name: glshort/MapReady


text:
/* -*- linux-c -*- */
/* fewbody_hier.c

   Copyright (C) 2002-2004 John M. Fregeau

   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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

#include <stdio.h>
#include <stddef.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include "fewbody.h"

/* allocate memory for a hier_t */
void fb_malloc_hier(fb_hier_t *hier)
{
	int i;

	hier->hi = (int *) malloc(hier->nstarinit * sizeof(int)) - 1;
	fb_create_indices(hier->hi, hier->nstarinit);
	hier->narr = (int *) malloc((hier->nstarinit-1)*sizeof(int)) - 2;
	/* change from malloc to calloc to get rid of harmless valgrind errors...
	   the errors occur in fb_normalize(), where uninitialized memory is normalized -
	   these are clearly harmless */
	hier->hier = (fb_obj_t *) calloc(hier->hi[hier->nstarinit] + 1, sizeof(fb_obj_t));
	for (i=0; i<hier->hi[hier->nstarinit]+1; i++) {
		hier->hier[i].ncoll = 0;
		hier->hier[i].id = (long *) malloc(hier->nstarinit * sizeof(long));
	}
	hier->obj = (fb_obj_t **) malloc(hier->nstarinit * sizeof(fb_obj_t *));
}

/* initialize to a flat hier */
void fb_init_hier(fb_hier_t *hier)
{
	int i;

	hier->nobj = hier->nstar;

	for (i=2; i<=hier->nstar; i++) {
		hier->narr[i] = 0;
	}

	for (i=0; i<hier->nstar; i++) {
		hier->obj[i] = &(hier->hier[hier->hi[1]+i]);
	}
}

/* free memory */
void fb_free_hier(fb_hier_t hier)
{
	int i;

	for (i=0; i<hier.hi[hier.nstarinit]+1; i++) {
		free(hier.hier[i].id);
	}
	free(hier.hi+1);
	free(hier.narr+2);
	free(hier.hier);
	free(hier.obj);
}

/* trickle down hier */
void fb_trickle(fb_hier_t *hier, double t)
{
	int i, j;

	for (i=hier->nstar; i>=2; i--) {
		for (j=0; j<hier->narr[i]; j++) {
			fb_downsync(&(hier->hier[hier->hi[i] + j]), t);
		}
	}
}

/* trickle up hier */
void fb_elkcirt(fb_hier_t *hier, double t)
{
	int i, j;

	for (i=2; i<=hier->nstar; i++) {
		for (j=0; j<hier->narr[i]; j++) {
			fb_upsync(&(hier->hier[hier->hi[i] + j]), t);
		}
	}
}

/* created the index array */
int fb_create_indices(int *hi, int nstar)
{
	int i, j=0;

	for (i=1; i<=nstar; i++) {
		hi[i] = j;
		j += nstar/i;
	}

	return(0);
}

/* how many members are in the immediate hierarchy */
int fb_n_hier(fb_obj_t *obj)
{
	if (obj == NULL) {
		return(0);
	} else if ((obj->obj[0] == NULL) && (obj->obj[1] == NULL)) {
		return(1);
	} else if ((obj->obj[0]->obj[0] == NULL) && (obj->obj[0]->obj[1] == NULL) && (obj->obj[1]->obj[0] == NULL) && (obj->obj[1]->obj[1] == NULL)) {
		return(2);
	} else if (((obj->obj[0]->obj[0] == NULL) && (obj->obj[0]->obj[1] == NULL)) || ((obj->obj[1]->obj[0] == NULL) && (obj->obj[1]->obj[1] == NULL))) {
		return(3);
	} else {
		return(4);
	}
}

/* print the hierarchy information */
char *fb_sprint_hier(fb_hier_t hier, char string[FB_MAX_STRING_LENGTH])
{
	int i;

	snprintf(string, FB_MAX_STRING_LENGTH, "nstar=%d nobj=%d: ", hier.nstar, hier.nobj);
	for (i=0; i<hier.nobj; i++) {
		snprintf(&(string[strlen(string)]), FB_MAX_STRING_LENGTH-strlen(string), " %s", hier.obj[i]->idstring);
	}

	return(string);
}

/* print the hierarchy information in a nice human-readable format */
char *fb_sprint_hier_hr(fb_hier_t hier, char string[FB_MAX_STRING_LENGTH])
{
	int i;

	/* zero out string so strlen(string) returns 0 */
	string[0] = '\0';

	/* loop through objects */
	for (i=0; i<hier.nobj; i++) {
		/* put a dash inbetween labels */
		if (i > 0) {
			snprintf(&(string[strlen(string)]), FB_MAX_STRING_LENGTH-strlen(string), "-");
		}
		/* a name for each type of hierarchy */
		if (hier.obj[i]->n == 1) {
			snprintf(&(string[strlen(string)]), FB_MAX_STRING_LENGTH-strlen(string), "single");
		} else if (hier.obj[i]->n == 2) {
			snprintf(&(string[strlen(string)]), FB_MAX_STRING_LENGTH-strlen(string), "binary");
		} else if (hier.obj[i]->n == 3) {
			snprintf(&(string[strlen(string)]), FB_MAX_STRING_LENGTH-strlen(string), "triple");
		} else if (hier.obj[i]->n == 4) {
			snprintf(&(string[strlen(string)]), FB_MAX_STRING_LENGTH-strlen(string), "quadruple");
		} else if (hier.obj[i]->n == 5) {
			snprintf(&(string[strlen(string)]), FB_MAX_STRING_LENGTH-strlen(string), "quintuple");
		} else if (hier.obj[i]->n == 6) {
			snprintf(&(string[strlen(string)]), FB_MAX_STRING_LENGTH-strlen(string), "sextuple");
		} else if (hier.obj[i]->n == 7) {
			snprintf(&(string[strlen(string)]), FB_MAX_STRING_LENGTH-strlen(string), "septuple");
		} else if (hier.obj[i]->n == 8) {
			snprintf(&(string[strlen(string)]), FB_MAX_STRING_LENGTH-strlen(string), "octuple");
		} else if (hier.obj[i]->n == 9) {
			snprintf(&(string[strlen(string)]), FB_MAX_STRING_LENGTH-strlen(string), "nontuple");
		} else if (hier.obj[i]->n == 10) {
			snprintf(&(string[strlen(string)]), FB_MAX_STRING_LENGTH-strlen(string), "decituple");
		} else {
			snprintf(&(string[strlen(string)]), FB_MAX_STRING_LENGTH-strlen(string), "%d", hier.obj[i]->n);
			snprintf(&(string[strlen(string)]), FB_MAX_STRING_LENGTH-strlen(string), "tuple");
		}
	}

	return(string);
}

/* merge the object's properties up---calculate the binary's properties from the
   stars' properties */
void fb_upsync(fb_obj_t *obj, double t)
{
	int i;
	double m0, m1, x0[3], x1[3], xrel[3], v0[3], v1[3], vrel[3], E, l0[3], l1[3], l[3];
	double ypp[3], psi, ecc_anom, L[3], A[3];

	/* a little bit of paranoia here */
	obj->ncoll = 0;
	obj->id[0] = -1;

	/* set time of upsync */
	obj->t = t;

	/* create idstring */
	snprintf(obj->idstring, FB_MAX_STRING_LENGTH, "[%s %s]", obj->obj[0]->idstring, obj->obj[1]->idstring);

	/* update number of stars in hierarchy */
	obj->n = obj->obj[0]->n + obj->obj[1]->n;

	/* update mass */
	m0 = obj->obj[0]->m;
	m1 = obj->obj[1]->m;
	obj->m = m0 + m1;

	/* update position and velocity */
	for (i=0; i<3; i++) {
		obj->x[i] = (m0*(obj->obj[0]->x[i]) + m1*(obj->obj[1]->x[i])) / (obj->m);
		obj->v[i] = (m0*(obj->obj[0]->v[i]) + m1*(obj->obj[1]->v[i])) / (obj->m);
	}

	/* update semimajor axis */
	for (i=0; i<3; i++) {
		x0[i] = obj->obj[0]->x[i] - obj->x[i];
		x1[i] = obj->obj[1]->x[i] - obj->x[i];
		xrel[i] = obj->obj[0]->x[i] - obj->obj[1]->x[i];
		v0[i] = obj->obj[0]->v[i] - obj->v[i];
		v1[i] = obj->obj[1]->v[i] - obj->v[i];
		vrel[i] = obj->obj[0]->v[i] - obj->obj[1]->v[i];
	}

	E = 0.5*(m0*fb_dot(v0, v0) + m1*fb_dot(v1, v1)) - m0*m1/fb_mod(xrel);

	if (E >= 0.0) {
		fprintf(stderr, "fb_upsync: E = %g >= 0!\n", E);
		exit(1);
	}

	obj->a = -m0 * m1 / (2.0 * E);

	/* update angular momentum, Runge-Lenz vector, and eccentricity (using the Runge-Lenz vector) */
	fb_cross(x0, v0, l0);
	fb_cross(x1, v1, l1);

	for (i=0; i<3; i++) {
		L[i] = m0 * l0[i] + m1 * l1[i];
		l[i] = L[i] * (m0+m1)/(m0*m1);
	}

	/* -A = l x v + G M \hat r */
	fb_cross(vrel, l, A);
	for (i=0; i<3; i++) {
		A[i] -= (m0+m1) * xrel[i]/fb_mod(xrel);
	}

	obj->e = fb_mod(A)/(m0+m1);

	/* must be careful with circular orbits, which have a null A vector */
	if (obj->e == 0.0) {
		for (i=0; i<3; i++) {
			A[i] = xrel[i];
		}
	}

	/* define coord system */
	for (i=0; i<3; i++) {
		obj->Lhat[i] = L[i]/fb_mod(L); /* z-direction */
		obj->Ahat[i] = A[i]/fb_mod(A); /* x-direction */
	}
	fb_cross(obj->Lhat, obj->Ahat, ypp);

	/* and set mean anomaly */
	psi = atan2(fb_dot(x0, ypp), fb_dot(x0, obj->Ahat));
	ecc_anom = acos((obj->e + cos(psi)) / (1.0 + obj->e * cos(psi)));
	if (psi < 0.0) {
		ecc_anom = 2.0 * FB_CONST_PI - ecc_anom;
	}
	obj->mean_anom = ecc_anom - obj->e * sin(ecc_anom);
}

/* randomly orient binary */
void fb_randorient(fb_obj_t *obj, gsl_rng *rng)
{
	int i;
	double theta, phi, omega, xp[3], yp[3], zp[3];

	/* random angles */
	theta = acos(2.0 * gsl_rng_uniform(rng) - 1.0);
	phi = 2.0 * FB_CONST_PI * gsl_rng_uniform(rng);
	omega = 2.0 * FB_CONST_PI * gsl_rng_uniform(rng);
	obj->mean_anom = 2.0 * FB_CONST_PI * gsl_rng_uniform(rng);

	/* set up coordinate transformations */
	zp[0] = -sin(theta) * sin(phi);
	zp[1] = sin(theta) * cos(phi);
	zp[2] = cos(theta);

	xp[0] = cos(phi);
	xp[1] = sin(phi);
	xp[2] = 0.0;

	fb_cross(zp, xp, yp);

	for (i=0; i<3; i++) {
		obj->Lhat[i] = zp[i];
		obj->Ahat[i] = xp[i] * cos(omega) + yp[i] * sin(omega);
	}
}

/* merge the object's properties down---calculate the objs' properties from the
   binary's properties */
void fb_downsync(fb_obj_t *obj, double t)
{
	int i;
	double xpp[3], ypp[3], zpp[3], er[3], epsi[3];
	double a, e, m0, m1, omega, mean_anom, ecc_anom, psi, r0, r1, L, psidot, E, r0dot, r0dot2, r1dot;

	/* we're assuming here that the objs' masses are already set */
	a = obj->a;
	e = obj->e;
	m0 = obj->obj[0]->m;
	m1 = obj->obj[1]->m;

	/* angular frequency */
	omega = sqrt(obj->m/fb_cub(a));
	/* mean anomaly, between 0 and 2PI */
	mean_anom = (obj->mean_anom + omega * (t - obj->t)) / (2.0 * FB_CONST_PI);
	mean_anom = (mean_anom - floor(mean_anom)) * 2.0 * FB_CONST_PI;
	/* eccentric anomaly, from solving the Kepler equation */
	ecc_anom = fb_kepler(e, mean_anom);
	/* true anomaly, between 0 and 2PI */
	psi = acos((cos(ecc_anom) - e) / (1.0 - e * cos(ecc_anom)));
	/* this step is necessary because acos() returns a value between 0 and PI */
	if (ecc_anom > FB_CONST_PI) {
		psi = 2.0 * FB_CONST_PI - psi;
	}

	/* define vectors, based on angular momentum and Runge-Lenz vectors */
	for (i=0; i<3; i++) {
		zpp[i] = obj->Lhat[i];
		xpp[i] = obj->Ahat[i];
	}
	fb_cross(zpp, xpp, ypp);

	for (i=0; i<3; i++) {
		er[i] = xpp[i] * cos(psi) + ypp[i] * sin(psi);
		epsi[i] = ypp[i] * cos(psi) - xpp[i] * sin(psi);
	}

	/* determine binary params */
	r0 = a * (1.0 - fb_sqr(e)) / ((1.0 + e*cos(psi)) * (1.0 + m0/m1));
	r1 = a * (1.0 - fb_sqr(e)) / ((1.0 + e*cos(psi)) * (1.0 + m1/m0));

	L = m0 * m1 * sqrt((m0+m1)*a*(1.0-fb_sqr(e))) / (m0 + m1);
	psidot = L / (m0 * fb_sqr(r0) + m1 * fb_sqr(r1));
	E = -m0 * m1 / (2.0 * a);

	/* must be careful for circular orbits */
	r0dot2 = (E + m0*m1/(r0+r1)) * 2.0 / (m0 * (1.0+m0/m1)) - fb_sqr(r0*psidot);
	if (ecc_anom > FB_CONST_PI) {
		r0dot = -(r0dot2<=0.0?0.0:sqrt(r0dot2));
	} else {
		r0dot = (r0dot2<=0.0?0.0:sqrt(r0dot2));
	}
	r1dot = m0 * r0dot / m1;

	for (i=0; i<3; i++) {
		obj->obj[0]->x[i] = r0 * er[i] + obj->x[i];
		obj->obj[0]->v[i] = r0dot * er[i] + r0 * psidot * epsi[i] + obj->v[i];
		obj->obj[1]->x[i] = -r1 * er[i] + obj->x[i];
		obj->obj[1]->v[i] = -r1dot * er[i] - r1 * psidot * epsi[i] + obj->v[i];
	}
}

/* copy one object to another, being careful about any pointers */
void fb_objcpy(fb_obj_t *obj1, fb_obj_t *obj2)
{
	int i;
	long *longptr;

	for (i=0; i<obj2->ncoll; i++) {
		obj1->id[i] = obj2->id[i];
	}

	/* prevent any dangling pointers */
	longptr = obj1->id;
	*obj1 = *obj2;
	obj1->id = longptr;
}

max_stars_repo_name: gnodvi/cosmos

max_issues_repo_name: gnodvi/cosmos

max_forks_repo_name: gnodvi/cosmos


text:
/* file: del2mat.h */

#ifndef DEL2MAT_H
#define DEL2MAT_H

#include <petsc.h>
#include <petscvec.h>
#include <petscmat.h>

/* external Fortran 90 subroutine */
#define Del2Apply del2apply_
extern void Del2Apply(int*,double*,const double*,double*);

/* user data structure and routines
 * defining the matrix-free operator */

typedef struct {
  PetscInt    N;
  PetscScalar *F;
} Del2Mat;

#undef  __FUNCT__
#define __FUNCT__ "Del2Mat_mult"
/* y <- A * x */
PetscErrorCode Del2Mat_mult(Mat A, Vec x, Vec y)
{
  Del2Mat *ctx;
  const PetscScalar *xx;
  PetscScalar *yy;
  PetscErrorCode ierr;
  PetscFunctionBegin;
  ierr = MatShellGetContext(A,(void**)&ctx);CHKERRQ(ierr);
  /* get raw vector arrays */
  ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr);
  ierr = VecGetArray(y,&yy);CHKERRQ(ierr);
  /* call external Fortran subroutine */
  Del2Apply(&ctx->N,ctx->F,xx,yy);
  /* restore raw vector arrays */
  ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr);
  ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}
#undef  __FUNCT__
#define __FUNCT__ "Del2Mat_diag"
/*D_i <- A_ii */
PetscErrorCode Del2Mat_diag(Mat A, Vec D)
{
  PetscErrorCode ierr;
  PetscFunctionBegin;
  ierr = VecSet(D,6.0);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#endif


text:
/**
 * \author Sylvain Marsat, University of Maryland - NASA GSFC
 *
 * \brief C header for the computation of the Fourier-domain overlaps, likelihoods.
 *
 *
 */

#ifndef _LIKELIHOOD_H
#define _LIKELIHOOD_H

#define _XOPEN_SOURCE 500

#ifdef __GNUC__
#define UNUSED __attribute__ ((unused))
#else
#define UNUSED
#endif

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>
#include <time.h>
#include <unistd.h>
#include <getopt.h>
#include <stdbool.h>
#include <string.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_complex.h>

#include "constants.h"
#include "struct.h"
#include "waveform.h"
#include "wip.h"

#if defined(__cplusplus)
extern "C" {
#elif 0
} /* so that editors will match preceding brace */
#endif

/****** Prototypes: utilities *******/

/* Function to evaluate a Noise function  */
void EvaluateNoise(
  gsl_vector* noisevalues,                         /* Output: vector of the noise values */
  gsl_vector* freq,                                /* Input: vector of frequencies on which to evaluate */
  ObjectFunction * Snoise,                         /* Noise function */
  double fLow,                                     /* Lower bound of the frequency window for the detector */
  double fHigh);                                   /* Upper bound of the frequency window for the detector */

/* Function building a frequency vector with logarithmic sampling */
void SetLogFrequencies(
  gsl_vector* freqvector,    /* Output pointer to gsl_vector, already allocated */
  const double fmin,         /* Lower bound of the frequency interval */
  const double fmax,         /* Upper bound of the frequency interval */
  const int nbpts);          /* Number of points */
/* Function building a frequency vector with linear sampling */
void SetLinearFrequencies(
  gsl_vector* freqvector,    /* Output pointer to gsl_vector, already allocated */
  const double fmin,         /* Lower bound of the frequency interval */
  const double fmax,         /* Upper bound of the frequency interval */
  const int nbpts);          /* Number of points */

/* Function determining logarithmically spaced frequencies from a waveform given as a list of modes, a minimal frequency and a number of points */
void ListmodesSetFrequencies(
  struct tagListmodesCAmpPhaseFrequencySeries *list,     /* Waveform, list of modes in amplitude/phase form */
  double fLow,                                           /* Additional lower frequency limit */
  double fHigh,                                          /* Additional upper frequency limit */
  int nbpts,                                             /* Number of frequency samples */
  int tagsampling,                                       /* Tag for linear (0) or logarithmic (1) sampling */
  gsl_vector* freqvector);                               /* Output: vector of frequencies, already allocated with nbpts */

/****** Prototypes: overlaps and likelihood computation *******/

/***************************** Functions for overlaps using linear integration ******************************/

/* Function computing the overlap (h1|h2) between two given modes, for a given noise function - uses simple trapeze integration on logarithmically sampled frequencies  */
double FDSinglemodeLogLinearOverlap(
  struct tagCAmpPhaseFrequencySeries *freqseries1, /* First mode h1, in amplitude/phase form */
  struct tagCAmpPhaseFrequencySeries *freqseries2, /* Second mode h2, in amplitude/phase form */
  ObjectFunction * Snoise,                         /* Noise function */
  double fLow,                                     /* Lower bound of the frequency window for the detector */
  double fHigh);                                   /* Upper bound of the frequency window for the detector */

/* Function computing the overlap (h1|h2) between two waveforms given as lists of mode contributions (factos sYlm already included), for a given noise function - uses simple trapeze integration on logarithmically sampled frequencies - generates the frequency series in Re/Im form by summing the mode contributions first, then computes the overlap */
double FDListmodesLogLinearOverlap(
  struct tagListmodesCAmpPhaseFrequencySeries *list1,    /* First waveform, list of modes in amplitude/phase form */
  struct tagListmodesCAmpPhaseFrequencySeries *list2,    /* First waveform, list of modes in amplitude/phase form */
  ObjectFunction * Snoise,                         /* Noise function */
  double fLow,                                           /* Lower bound of the frequency window for the detector */
  double fHigh,                                          /* Upper bound of the frequency window for the detector */
  double fstartobs1,                                     /* Starting frequency for the 22 mode of wf 1 - as determined from a limited duration of the observation - set to 0 to ignore */
  double fstartobs2);                                    /* Starting frequency for the 22 mode of wf 2 - as determined from a limited duration of the observation - set to 0 to ignore */

/* Function computing the overlap (h1|h2) between h1 given in Re/Im form and h2 given as a list of mode contributions (factos sYlm already included), for a given vector of noise values - uses simple trapeze integration - generates the frequency series of h2 in Re/Im form by summing the mode contributions first, then computes the overlap */
double FDOverlapReImvsListmodesCAmpPhase(
  struct tagReImFrequencySeries *freqseries1,            /* First waveform, in Re/Im form */
  struct tagListmodesCAmpPhaseFrequencySeries *list2,    /* Second waveform, list of modes in amplitude/phase form */
  gsl_vector* noisevalues,                               /* Vector for the noise values on the freq of h1 */
  double fLow,                                           /* Minimal frequency - set to 0 to ignore */
  double fHigh,                                          /* Maximal frequency - set to 0 to ignore */
  double fstartobs2);                                    /* Starting frequency for the 22 mode of wf 2 - as determined from a limited duration of the observation - set to 0 to ignore */

/* Function computing the overlap (h1|h2) between two waveforms given as Re/Im frequency series (common freq values), for a given vector of noise values - uses simple trapeze integration */
double FDOverlapReImvsReIm(
  struct tagReImFrequencySeries *h1,  /* First waveform, frequency series in Re/Im form */
  struct tagReImFrequencySeries *h2,  /* Second waveform, frequency series in Re/Im form */
  gsl_vector* noisevalues);           /* Vector for the noise values on common freq of the freqseries */

/***************************** Functions for overlaps using amplitude/phase (wip) ******************************/

/* Function computing the overlap (h1|h2) between two given modes in amplitude/phase form, for a given noise function - uses the amplitude/phase representation (wip) */
double FDSinglemodeWIPOverlap(
  struct tagCAmpPhaseFrequencySeries *freqseries1, /* First mode h1, in amplitude/phase form */
  struct tagCAmpPhaseFrequencySeries *freqseries2, /* Second mode h2, in amplitude/phase form */
  ObjectFunction * Snoise,                         /* Noise function */
  double fLow,                                     /* Lower bound of the frequency window for the detector */
  double fHigh);                                   /* Upper bound of the frequency window for the detector */

/* Function computing the overlap (h1|h2) between two waveforms given as list of modes, for a given noise function - no cross-products between modes are taken into account - two additional parameters for the starting 22-mode frequencies (then properly scaled for the other modes) for a limited duration of the observations */
double FDListmodesWIPOverlap(
  struct tagListmodesCAmpPhaseFrequencySeries *listh1, /* First waveform, list of modes in amplitude/phase form */
  struct tagListmodesCAmpPhaseFrequencySeries *listh2, /* Second waveform, list of modes in amplitude/phase form */
  ObjectFunction * Snoise,                         /* Noise function */
  double fLow,                                         /* Lower bound of the frequency window for the detector */
  double fHigh,                                        /* Upper bound of the frequency window for the detector */
  double fstartobs1,                                   /* Starting frequency for the 22 mode of wf 1 - as determined from a limited duration of the observation - set to 0 to ignore */
  double fstartobs2);                                  /* Starting frequency for the 22 mode of wf 2 - as determined from a limited duration of the observation - set to 0 to ignore */

/************** Functions for overlap/likelihood allowing to switch between wip and loglinear integration *****************/

/* Wrapping of FDListmodesWIPOverlap or FDListmodesLogLinearOverlap according to tagint */
double FDListmodesOverlap(
  struct tagListmodesCAmpPhaseFrequencySeries *listh1, /* First mode h1, list of modes in amplitude/phase form */
  struct tagListmodesCAmpPhaseFrequencySeries *listh2, /* Second mode h2, list of modes in amplitude/phase form */
  ObjectFunction * Snoise,                         /* Noise function */
  double fLow,                                         /* Lower bound of the frequency window for the detector */
  double fHigh,                                        /* Upper bound of the frequency window for the detector */
  double fstartobs1,                                   /* Starting frequency for the 22 mode of wf 1 - as determined from a limited duration of the observation - set to 0 to ignore */
  double fstartobs2,                                   /* Starting frequency for the 22 mode of wf 2 - as determined from a limited duration of the observation - set to 0 to ignore */
  int tagint);                                         /* Tag choosing the integrator: 0 for wip, 1 for log linear integration */

/* Function computing the log likelihood (h|s) - 1/2 (h|h) - 1/2 (s|s), with s the signal, h the template, and where we keep the constant term (s|s) - passed to the function as a parameter - all the cross-products between modes are taken into account - two additionals parameters for the starting 22-mode frequencies (then properly scaled for the other modes) for a limited duration of the observations */
double FDLogLikelihood(
  struct tagListmodesCAmpPhaseFrequencySeries *lists,  /* Input: list of modes for the signal s, in Frequency-domain amplitude and phase form */
  struct tagListmodesCAmpPhaseFrequencySeries *listh,  /* Input: list of modes for the template, in Frequency-domain amplitude and phase form */
  ObjectFunction * Snoise,                         /* Noise function */
  double fLow,                                         /* Lower bound of the frequency window for the detector */
  double fHigh,                                        /* Upper bound of the frequency window for the detector */
  double ss,                                           /* Inner product (s|s), constant to be computed elsewhere and passed as an argument */
  double hh,                                           /* Inner product (h|h), constant to be computed elsewhere and passed as an argument */
  double fstartobss,                                   /* Starting frequency for the 22 mode of s - as determined from a limited duration of the observation - set to 0 to ignore */
  double fstartobsh,                                   /* Starting frequency for the 22 mode of h - as determined from a limited duration of the observation - set to 0 to ignore */
  int tagint);                                         /* Tag choosing the integrator: 0 for wip, 1 for log linear integration */

/************** Functions for overlap/likelihood specific to loglinear integration *****************/

/* Function computing the log likelihood -1/2(h-s|h-s), with s the signal, h the template (both given as frequency series in Re/Im form) - h, s assumed to be given on the same set of frequencies - same for the vector of noise values - fstartobs for s, h has already been taken into account */
double FDLogLikelihoodReIm(
  struct tagReImFrequencySeries *s,    /* First waveform (injection), frequency series in Re/Im form */
  struct tagReImFrequencySeries *h,    /* Second waveform (template), frequency series in Re/Im form */
  gsl_vector* noisevalues);            /* Vector for the noise values on common freq of the freqseries */

/***************************** Functions for overlaps using amplitude/phase (Fresnel) ******************************/

void ComputeIntegrandValues(
  CAmpPhaseFrequencySeries** integrand,     /* Output: values of the integrand on common frequencies (initialized in the function) */
  CAmpPhaseFrequencySeries* freqseries1,    /* Input: frequency series for wf 1 */
  CAmpPhaseSpline* splines2,                /* Input: splines in matrix form for wf 2 */
  ObjectFunction * Snoise,                         /* Noise function */
  double fLow,                              /* Lower bound of the frequency - 0 to ignore */
  double fHigh);                            /* Upper bound of the frequency - 0 to ignore */
/* Function computing the integrand values, combining the three channels A, E and T */
int ComputeIntegrandValues3Chan(
  CAmpPhaseFrequencySeries** integrand,     /* Output: values of the integrand on common frequencies (initialized in the function) */
  CAmpPhaseFrequencySeries* freqseries1chan1,    /* Input: frequency series for wf 1, channel 1 */
  CAmpPhaseFrequencySeries* freqseries1chan2,    /* Input: frequency series for wf 1, channel 2 */
  CAmpPhaseFrequencySeries* freqseries1chan3,    /* Input: frequency series for wf 1, channel 3 */
  CAmpPhaseSpline* splines2chan1,                /* Input: splines in matrix form for wf 2, channel 1 */
  CAmpPhaseSpline* splines2chan2,                /* Input: splines in matrix form for wf 2, channel 2 */
  CAmpPhaseSpline* splines2chan3,                /* Input: splines in matrix form for wf 2, channel 3 */
  ObjectFunction * Snoise1,                         /* Noise function */
  ObjectFunction * Snoise2,                         /* Noise function */
  ObjectFunction * Snoise3,                         /* Noise function */
  double fLow,                              /* Lower bound of the frequency - 0 to ignore */
  double fHigh);                            /* Upper bound of the frequency - 0 to ignore */

/* Function computing the overlap (h1|h2) between two given modes in amplitude/phase form, one being already interpolated, for a given noise function - uses the amplitude/phase representation (Fresnel) */
double FDSinglemodeFresnelOverlap(
  struct tagCAmpPhaseFrequencySeries *freqseries1, /* First mode h1, in amplitude/phase form */
  struct tagCAmpPhaseSpline *splines2,             /* Second mode h2, already interpolated in matrix form */
  ObjectFunction * Snoise,                         /* Noise function */
  double fLow,                                     /* Lower bound of the frequency window for the detector */
  double fHigh);                                   /* Upper bound of the frequency window for the detector */
/* Function computing the overlap (h1|h2) between two given modes in amplitude/phase form for each channel A, E, T, one being already interpolated, for a given noise function - uses the amplitude/phase representation (Fresnel) */
double FDSinglemodeFresnelOverlap3Chan(
  struct tagCAmpPhaseFrequencySeries *freqseries1chan1, /* First mode h1 for channel 1, in amplitude/phase form */
  struct tagCAmpPhaseFrequencySeries *freqseries1chan2, /* First mode h1 for channel 2, in amplitude/phase form */
  struct tagCAmpPhaseFrequencySeries *freqseries1chan3, /* First mode h1 for channel 3, in amplitude/phase form */
  struct tagCAmpPhaseSpline *splines2chan1,             /* Second mode h2 for channel 1, already interpolated in matrix form */
  struct tagCAmpPhaseSpline *splines2chan2,             /* Second mode h2 for channel 2, already interpolated in matrix form */
  struct tagCAmpPhaseSpline *splines2chan3,             /* Second mode h2 for channel 3, already interpolated in matrix form */
  ObjectFunction * Snoise1,                         /* Noise function */
  ObjectFunction * Snoise2,                         /* Noise function */
  ObjectFunction * Snoise3,                         /* Noise function */
  double fLow,                                      /* Lower bound of the frequency window for the detector */
  double fHigh);                                    /* Upper bound of the frequency window for the detector */

/* Function computing the overlap (h1|h2) between two waveforms given as list of modes, one being already interpolated, for a given noise function - two additional parameters for the starting 22-mode frequencies (then properly scaled for the other modes) for a limited duration of the observations */
double FDListmodesFresnelOverlap(
  struct tagListmodesCAmpPhaseFrequencySeries *listh1, /* First waveform, list of modes in amplitude/phase form */
  struct tagListmodesCAmpPhaseSpline *listsplines2,    /* Second waveform, list of modes already interpolated in matrix form */
  ObjectFunction * Snoise,                         /* Noise function */
  double fLow,                                         /* Lower bound of the frequency window for the detector */
  double fHigh,                                        /* Upper bound of the frequency window for the detector */
  double fstartobs1,                                   /* Starting frequency for the 22 mode of wf 1 - as determined from a limited duration of the observation - set to 0 to ignore */
  double fstartobs2);                                  /* Starting frequency for the 22 mode of wf 2 - as determined from a limited duration of the observation - set to 0 to ignore */
/* Function computing the overlap (h1|h2) between two waveforms given as list of modes for each non-correlated channel 1,2,3, one being already interpolated, for a given noise function - two additional parameters for the starting 22-mode frequencies (then properly scaled for the other modes) for a limited duration of the observations */
double FDListmodesFresnelOverlap3Chan(
  struct tagListmodesCAmpPhaseFrequencySeries *listh1chan1, /* First waveform channel channel 1, list of modes in amplitude/phase form */
  struct tagListmodesCAmpPhaseFrequencySeries *listh1chan2, /* First waveform channel channel 2, list of modes in amplitude/phase form */
  struct tagListmodesCAmpPhaseFrequencySeries *listh1chan3, /* First waveform channel channel 3, list of modes in amplitude/phase form */
  struct tagListmodesCAmpPhaseSpline *listsplines2chan1,    /* Second waveform channel channel 1, list of modes already interpolated in matrix form */
  struct tagListmodesCAmpPhaseSpline *listsplines2chan2,    /* Second waveform channel channel 2, list of modes already interpolated in matrix form */
  struct tagListmodesCAmpPhaseSpline *listsplines2chan3,    /* Second waveform channel channel 3, list of modes already interpolated in matrix form */
  ObjectFunction * Snoise1,                         /* Noise function */
  ObjectFunction * Snoise2,                         /* Noise function */
  ObjectFunction * Snoise3,                         /* Noise function */
  double fLow,                                          /* Lower bound of the frequency window for the detector */
  double fHigh,                                         /* Upper bound of the frequency window for the detector */
  double fstartobs1,                                    /* Starting frequency for the 22 mode of wf 1 - as determined from a limited duration of the observation - set to 0 to ignore */
  double fstartobs2);                                    /* Starting frequency for the 22 mode of wf 2 - as determined from a limited duration of the observation - set to 0 to ignore */

#if 0
{ /* so that editors will match succeeding brace */
#elif defined(__cplusplus)
}
#endif

#endif /* _LIKELIHOOD_H */

max_stars_repo_name: titodalcanton/flare

max_issues_repo_name: titodalcanton/flare

max_forks_repo_name: titodalcanton/flare


text:
/**
  * Copyright 2016 José Manuel Abuín Mosquera <josemanuel.abuin@usc.es>
  *
  * This file is part of Matrix Market Suite.
  *
  * Matrix Market Suite is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
  * the Free Software Foundation, either version 3 of the License, or
  * (at your option) any later version.
  *
  * Matrix Market Suite 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 Matrix Market Suite. If not, see <http://www.gnu.org/licenses/>.
  */

#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>

#include <cblas.h>

#include "DMxV.h"

void usageDMxDM(){

	fprintf(stderr, "\n");
	fprintf(stderr, "Usage: MM-Suite DMxDM [options] <input-matrix-A> <input-matrix-B>\n");
	fprintf(stderr, "\nInput/output options:\n\n");
	fprintf(stderr, "       -o STR        Output file name. Default: stdout\n");
	fprintf(stderr, "       -r            Input format is row per line. Default: False\n");
	fprintf(stderr, "\nPerformance options:\n\n");
	fprintf(stderr, "       -t INT        Number of threads to use in OpenBLAS. Default: 1\n");
	fprintf(stderr, "\n");

}

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

	int 			ret_code = 1;
	int 			option;

	unsigned long 		*IA;
	unsigned long 		*JA;
	double 			*valuesA;

	unsigned long 		MA;
	unsigned long 		NA;
	unsigned long long 	nzA;

	unsigned long 		*IB;
	unsigned long 		*JB;
	double 			*valuesB;

	unsigned long 		MB;
	unsigned long 		NB;
	unsigned long long 	nzB;

	char			*outputFileName = NULL;

	char			*inputMatrixFileA = NULL;
	char			*inputMatrixFileB= NULL;

	int			inputFormatRow = 0;

	int			numThreads = 1;

	while ((option = getopt(argc, argv,"ro:t:")) >= 0) {
		switch (option) {
			case 'o' :
				//free(outputFileName);

				outputFileName = (char *) malloc(sizeof(char)*strlen(optarg)+1);
				strcpy(outputFileName,optarg);

				break;

			case 'r':
				inputFormatRow = 1;
				break;

			case 't':
				numThreads = atoi(optarg);
				break;

			default: break;
		}

	}

	if ((optind + 2 > argc) || (optind + 3 <= argc)) {
		usageDMxDM();
		return 0;
	}

	openblas_set_num_threads(numThreads);

	if(outputFileName == NULL) {
		outputFileName = (char *) malloc(sizeof(char)*7);
		sprintf(outputFileName,"stdout");
	}

	inputMatrixFileA = (char *)malloc(sizeof(char)*strlen(argv[optind])+1);
	inputMatrixFileB = (char *)malloc(sizeof(char)*strlen(argv[optind+1])+1);

	strcpy(inputMatrixFileA,argv[optind]);
	strcpy(inputMatrixFileB,argv[optind+1]);

	//Read matrices
	if(inputFormatRow){

		if(!readDenseCoordinateMatrixRowLine(inputMatrixFileA,&IA,&JA,&valuesA,&MA,&NA,&nzA)){
			fprintf(stderr, "[%s] Can not read Matrix\n",__func__);
			return 0;
		}

		if(!readDenseCoordinateMatrixRowLine(inputMatrixFileB,&IB,&JB,&valuesB,&MB,&NB,&nzB)){
			fprintf(stderr, "[%s] Can not read Matrix\n",__func__);
			return 0;
		}

	}
	else {
		if(!readDenseCoordinateMatrix(inputMatrixFileA,&IA,&JA,&valuesA,&MA,&NA,&nzA)){
			fprintf(stderr, "[%s] Can not read Matrix\n",__func__);
			return 0;
		}

		if(!readDenseCoordinateMatrix(inputMatrixFileB,&IB,&JB,&valuesB,&MB,&NB,&nzB)){
			fprintf(stderr, "[%s] Can not read Matrix\n",__func__);
			return 0;
		}
	}

	/*
	void cblas_dgemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA,
                 const enum CBLAS_TRANSPOSE TransB, const int M, const int N,
                 const int K, const double alpha, const double *A,
                 const int lda, const double *B, const int ldb,
                 const double beta, double *C, const int ldc);
                 */

        double *result=(double *) malloc(MA * NB * sizeof(double));

	//cblas_dgemv(CblasColMajor,CblasNoTrans,M,N,1.0,values,N,vectorValues,1,0.0,result,1);
	cblas_dgemm(CblasRowMajor,CblasNoTrans,CblasNoTrans, MA, NB, NA, 1.0,valuesA,NA,valuesB,NB,0.0,result,NB);

	if(inputFormatRow){
		writeDenseCoordinateMatrixRowLine(outputFileName, result,MA,NB,MA * NB);
	}
	else{
		writeDenseCoordinateMatrix(outputFileName, result,MA,NB,MA * NB);
	}

	return ret_code;
}


text:
/* This library is used by compare_hashes_* tests.
 */

#ifndef COMPARISON_STATS_H
#define COMPARISON_STATS_H

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdint.h>
#include <gsl/gsl_sf_log.h>
#include <gsl/gsl_math.h>

struct comparison_stats {
	uint64_t *ones_count;
	uint64_t tot_count;
	uint hash_width;
	uint64_t (*calc_hash)(void *hash_info, unsigned char *str, unsigned length);
	void *hash_info;
};

typedef enum {
	report_type_short,
	report_type_long } report_type;

void basic_init_comparison_stats(struct comparison_stats *cs, uint width);
void process_comparison_stats(struct comparison_stats *cs, unsigned char *str, unsigned sz);
void report_comparison_stats(struct comparison_stats *cs, const char *name, report_type rtyp);
void free_comparison_stats(struct comparison_stats *cs);

#endif

max_stars_repo_name: jonkanderson/BitBalancedTableHash

max_issues_repo_name: jonkanderson/BitBalancedTableHash

max_forks_repo_name: jonkanderson/BitBalancedTableHash


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

void
cblas_srotg (float *a, float *b, float *c, float *s)
{
#define BASE float
#include "source_rotg.h"
#undef BASE
}

max_stars_repo_licenses: Net-SNMP Xnet

max_issues_repo_licenses: Net-SNMP Xnet

max_forks_repo_licenses: Net-SNMP Xnet


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

#include <config.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_movstat.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_statistics.h>

/*
gsl_movstat_qqr()
  Apply a moving q-quantile range to an input vector

Inputs: endtype - how to handle end points
        x       - input vector, size n
        q       - quantile \in [0,0.5]
        xqqr    - (output) vector of q-quantile ranges of x, size n
                  xqqr_i = q-quantile range of i-th window:
        w       - workspace
*/

int
gsl_movstat_qqr(const gsl_movstat_end_t endtype, const gsl_vector * x, const double q,
                gsl_vector * xqqr, gsl_movstat_workspace * w)
{
  if (x->size != xqqr->size)
    {
      GSL_ERROR("x and xqqr vectors must have same length", GSL_EBADLEN);
    }
  else if (q < 0.0 || q > 0.5)
    {
      GSL_ERROR("q must be between 0 and 0.5", GSL_EDOM);
    }
  else
    {
      double qq = q;
      int status = gsl_movstat_apply_accum(endtype, x, gsl_movstat_accum_qqr, (void *) &qq, xqqr, NULL, w);
      return status;
    }
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
/*
 * Calculation of Integral
 */
#include <stdio.h>

#include <gsl/gsl_rng.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_vegas.h>

#include "timer.h"

double g (double *t, size_t dim, void *params);

int main (void)
{
    double res, err;
    double dist, dmin = 1.001, dmax = 4.;
    size_t dim = 6;
    int np = 20;
    double vegas[20], vegaserr[20];
    double time1, time2;

    double dstep = (dmax - dmin) / (np - 1);

    double xl[] = { 0., 0., 0., 0., 0., 0. };
    double xu[] = { 1., 1., 1., 1., 1., 1. };

    gsl_rng *r = gsl_rng_alloc (gsl_rng_taus2);
    unsigned long seed = 1UL;

    gsl_rng_set (r, seed);

    size_t calls = 1000000;

    gsl_monte_function G = { &g, dim, &dist };

    gsl_monte_vegas_state *sv = gsl_monte_vegas_alloc (dim);

    gsl_monte_vegas_init (sv);

    dist = dmin;

// Vegas algorithm calculation

    timer_start ();

    // printf ("# Dist      Energy     ErrEst    Dipole\n");
    for (int i = 0; i < np; i++)
    {
        gsl_monte_vegas_integrate (&G, xl, xu, dim, calls / 5, r, sv, &res,
            &err);

        do
        {
            gsl_monte_vegas_integrate (&G, xl, xu, dim, calls, r, sv, &res,
                &err);
        }

        while (fabs (gsl_monte_vegas_chisq (sv) - 1.0) > 0.2);

        //printf ("% .6f  % .6f % .6f  % .6f\n", dist, res, err,
        //    -2. / pow (dist, 3.));

        //fflush (stdout);
        vegas[i] = res;
        vegaserr[i] = err;

        dist += dstep;
    }

    time1 = timer_stop ();

    gsl_monte_vegas_free (sv);

    dist = dmin;
    printf ("# Dist      HMC        Vegas      Verr       Dipole    Difference\n");

    double t[6];

// Calculation using homemade Monte Carlo algorithm

    timer_start ();

    for (int i = 0; i < np; i++)
    {
        double sum = 0.;

        for (int j = 0; j < (int) calls; j++)
        {
            for (int m = 0; m < (int) dim; m++)
            {
                t[m] = gsl_rng_uniform (r);
            }
            sum += g (t, dim, &dist);
        }
        double energy = sum / calls;
        double difference = fabs (energy - vegas[i]);

        printf ("% .6f  % .6f % .6f  % .6f  % .6f  % .6f\n",
            dist, energy, vegas[i], vegaserr[i], -2./pow(dist, 3.), difference );
        dist += dstep;

     }

   time2 = timer_stop ();
    printf ("Time Vegas: %.2f Time HMC: %.2f\n", time1, time2);
    printf ("Speedup: %.2f\n", time1 / time2);

    gsl_rng_free (r);

    return 0;
}


text:
#include <cblas.h>
#include <stdio.h>
#include <stdlib.h>
#include <papi.h>

#ifndef MATLEN
#define MATLEN 128
#endif

#ifndef QUIET
#define PRINT( exp ) exp
#else
#define PRINT( exp )
#endif

int main( ) {
  int EventSet = PAPI_NULL;
  long long values[ 4 ], s, e;
  int retval;
  double *a, *b, *c;

  /* INICIALIZAÇÃO */

  PRINT( printf( "Inicializando Matriz: %dx%d\n", MATLEN, MATLEN ) );

  a = ( double* ) malloc( MATLEN * MATLEN * sizeof( double ) );
  b = ( double* ) malloc( MATLEN * MATLEN * sizeof( double ) );
  c = ( double* ) malloc( MATLEN * MATLEN * sizeof( double ) );

  /*
   * CONFIGURAÇÃO DO PAPI
   * Init PAPI library
   */
  retval = PAPI_library_init( PAPI_VER_CURRENT );
  if( retval != PAPI_VER_CURRENT ) {
    printf( "Erro em PAPI_library_init : retval = %d\n", retval );
    exit( 1 );
  }
  if( ( retval = PAPI_create_eventset( &EventSet ) ) != PAPI_OK ) {
    printf( "Erro em PAPI_create_eventset : retval = %d\n", retval );
    exit( 1 );
  }
  if( PAPI_add_event( EventSet, PAPI_L2_DCM ) != PAPI_OK ) {
    printf( "Erro em PAPI_L2_DCM\n" );
    exit( 1 );
  }
  if( PAPI_add_event( EventSet, PAPI_DP_OPS ) != PAPI_OK ) {
    printf( "Erro em PAPI_DP_OPS\n" );
    exit( 1 );
  }
  if( PAPI_add_event( EventSet, PAPI_TOT_CYC ) != PAPI_OK ) {
    printf( "Erro em PAPI_TOT_CYC\n" );
    exit( 1 );
  }
  if( PAPI_add_event( EventSet, PAPI_TOT_INS ) != PAPI_OK ) {
    printf( "Erro em PAPI_TOT_INS\n" );
    exit( 1 );
  }
  if( ( retval = PAPI_start( EventSet ) ) != PAPI_OK ) {
    printf( "Erro em PAPI_start" );
    exit( 1 );
  }
  s = PAPI_get_real_usec( );
  /* FUNÇÃO A SER AVALIADA */

  cblas_dgemm( CblasRowMajor, CblasNoTrans, CblasNoTrans,
               MATLEN, MATLEN, MATLEN, 1.0, a, MATLEN,
               b, MATLEN, 0.0, c, MATLEN );

  /* FIM DA FUNÇÃO A SER AVALIADA */
  e = PAPI_get_real_usec( );
  if( ( retval = PAPI_read( EventSet, &values[ 0 ] ) ) != PAPI_OK ) {
    printf( "Erro em PAPI_read" );
    exit( 1 );
  }
  if( ( retval = PAPI_stop( EventSet, NULL ) ) != PAPI_OK ) {
    printf( "Erro em PAPI_stop" );
    exit( 1 );
  }
  double cpi = ( double ) values[ 2 ] / ( double ) values[ 3 ];
  double icp = ( double ) values[ 3 ] / ( double ) values[ 2 ];
  double mflops = ( double ) values[ 1 ];
  // double mflops = ( double ) 2 * MATLEN * MATLEN * MATLEN;
  mflops = ( mflops / ( ( double ) ( e - s ) ) );
  /* EXIBINDO INFORMAÇÕES */
  PRINT(
    printf( "PAPI_L2_DCM = %lld\n", values[ 0 ] );
    printf( "PAPI_DP_OPS = %lld\n", values[ 1 ] );

    /* CPI */
    printf( "PAPI_TOT_CYC = %lld\n", values[ 2 ] );
    printf( "PAPI_TOT_INS = %lld\n", values[ 3 ] );
    printf( "CPI: %.2f\n", cpi );
    printf( "ICP: %.2f\n", icp );

    printf( "Wallclock time: %lld ms\n", e - s );
    printf( "MFLOPS: %g\n", mflops );
    printf( "Fim\n" );
    );
  /*       MAT BLk Time  DCM   MFLOPS CPI */
  printf( "%d, %lld, %lld, %.2f, %.2f\n", MATLEN, e - s, values[ 0 ], mflops, cpi );
  free( a );
  free( b );
  free( c );
  return( 0 );
}

max_stars_repo_name: lellisls/PAD

max_issues_repo_name: lellisls/PAD

max_forks_repo_name: lellisls/PAD


text:
#include "stdlib.h"
#include "stdio.h"
#include "/home/lillian/work/install_fpdebug/valgrind-3.7.0/fpdebug/fpdebug.h"
#include <gsl/gsl_sf.h>
int main(int argc, const char * argv[]) {
unsigned long int hexdouble;
double a;
sscanf(argv[1], "%lX", &hexdouble);
a = *(double*)(&hexdouble);
double result = gsl_sf_Si(a);
//printf("%.15f\n", result);
VALGRIND_PRINT_ERROR("result", &result);
return 0;
}

max_stars_repo_name: floatfeather/FpGenetic

max_issues_repo_name: floatfeather/FpGenetic

max_forks_repo_name: floatfeather/FpGenetic


text:
#pragma once

#include <gsl/gsl>
#include <memory>

struct SpeexResamplerState_;
typedef struct SpeexResamplerState_ SpeexResamplerState;

namespace Halley
{
	struct AudioResamplerResult
	{
		size_t nRead;
		size_t nWritten;
	};

	class AudioResampler
	{
	public:
		AudioResampler(int from, int to, int nChannels, float quality = 1.0f);
		~AudioResampler();

		AudioResamplerResult resample(gsl::span<const float> src, gsl::span<float> dst, size_t channel);
		AudioResamplerResult resampleInterleaved(gsl::span<const float> src, gsl::span<float> dst);
		AudioResamplerResult resampleInterleaved(gsl::span<const short> src, gsl::span<short> dst);
		AudioResamplerResult resampleNoninterleaved(gsl::span<const float> src, gsl::span<float> dst, const size_t numChannels);
		size_t numOutputSamples(size_t numInputSamples) const;

	private:
		std::unique_ptr<SpeexResamplerState, void(*)(SpeexResamplerState*)> resampler;
		size_t nChannels;
		int from;
		int to;
	};
}

max_stars_repo_name: sunhay/halley

max_issues_repo_name: sunhay/halley

max_forks_repo_name: sunhay/halley


text:
#include <assert.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <sys/time.h>

#include "../src/fastlambertw.h"

#include "testmacros.h"

#include "../src/config.h"
#ifdef HAVE_GSL_GSL_SF_LAMBERT_H
#include <gsl/gsl_sf_lambert.h>
#endif

static inline float
lambertwrange (void)
{
  if (drand48 () < 0.5)
    {
      return -0.36787944117144232f + 1.36787944117144232f * drand48 ();
    }
  else
    {
      return 100.0f * drand48 ();
    }
}

static inline float
lambertwf (float x)
{
  float w = (x < 5) ? 0 : log (x) - log (log (x)) + log (log (x)) / log (x);
  unsigned int n;

  for (n = 0; n < 20; ++n)
    {
      w = (w * w + exp (-w) * x) / (1.0 + w);
    }

  return w;
}

static inline float
lambertwexpxf (float x)
{
  return lambertwf (expf (x));
}

test_scalar (fastlambertw, lambertwf, lambertwrange (), 1e-4f, 100000000)
test_scalar (fasterlambertw, lambertwf, lambertwrange (), 1e-2f, 100000000)
test_scalar (fastlambertwexpx, lambertwexpxf, -3.0f + 6.0f * drand48 (), 1e-3f, 100000000)
test_scalar (fasterlambertwexpx, lambertwexpxf, -3.0f + 6.0f * drand48 (), 1e-2f, 100000000)

#ifdef HAVE_GSL_GSL_SF_LAMBERT_H
test_scalar (gsl_sf_lambert_W0, lambertwf, lambertwrange (), 1e-2f, 1000000)
#endif

test_vector (vfastlambertw, lambertwf, lambertwrange (), 1e-4f, 100000000)
test_vector (vfasterlambertw, lambertwf, lambertwrange (), 1e-2f, 100000000)
test_vector (vfastlambertwexpx, lambertwexpxf, -3.0f + 6.0f * drand48 (), 1e-3f, 100000000)
test_vector (vfasterlambertwexpx, lambertwexpxf, -3.0f + 6.0f * drand48 (), 1e-2f, 100000000)

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

  (void) argc;

  srand48 (69);

//  fprintf (stderr, "fastlambertw (%g) = %g, "
//                   "fastlambertw (%g) = %g, "
//                   "fasterlambertwexpx (%g) = %g (%g)\n",
//           -0.36787944117144232f,
//           fastlambertw (-0.36787944117144232f),
//           -0.36787944117144232f + 0.01f,
//           fastlambertw (-0.36787944117144232f + 0.01f),
//           -5.0f,
//           fasterlambertwexpx (-5.0f),
//           v4sf_index (vfasterlambertwexpx (v4sfl (-5.0f)), 0));

  strncpy (buf, argv[0], sizeof (buf) - 5);
  strncat (buf, ".out", 5);

  fclose (stderr);
  stderr = fopen (buf, "w");

  test_fastlambertw ();
  test_fasterlambertw ();
  test_vfastlambertw ();
  test_vfasterlambertw ();
#ifdef HAVE_GSL_GSL_SF_LAMBERT_H
  test_gsl_sf_lambert_W0 ();
#endif
  test_fastlambertwexpx ();
  test_fasterlambertwexpx ();
  test_vfastlambertwexpx ();
  test_vfasterlambertwexpx ();

  time_fastlambertw ();
  time_fasterlambertw ();
  time_vfastlambertw ();
  time_vfasterlambertw ();
#ifdef HAVE_GSL_GSL_SF_LAMBERT_H
  time_gsl_sf_lambert_W0 ();
#endif
  time_fastlambertwexpx ();
  time_fasterlambertwexpx ();
  time_vfastlambertwexpx ();
  time_vfasterlambertwexpx ();

  return 0;
}

max_stars_repo_name: Mister-Meeseeks/fastapprox

max_issues_repo_name: Mister-Meeseeks/fastapprox

max_forks_repo_name: Mister-Meeseeks/fastapprox


text:
///
/// @file
///
/// @brief This file contains math functions that are shared among all
/// components of the library.
///
/// @author Mirko Myllykoski (mirkom@cs.umu.se), Umeå University
///
/// @internal LICENSE
///
/// Copyright (c) 2019-2020, Umeå Universitet
///
/// Redistribution and use in source and binary forms, with or without
/// modification, are permitted provided that the following conditions are met:
///
/// 1. Redistributions of source code must retain the above copyright notice,
///    this list of conditions and the following disclaimer.
///
/// 2. Redistributions in binary form must reproduce the above copyright notice,
///    this list of conditions and the following disclaimer in the documentation
///    and/or other materials provided with the distribution.
///
/// 3. Neither the name of the copyright holder nor the names of its
///    contributors may 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 HOLDER 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.
///

#include <starneig_config.h>
#include <starneig/configuration.h>
#include "math.h"
#include "common.h"
#include "sanity.h"
#include <stddef.h>
#include <string.h>
#include <cblas.h>

int starneig_largers_factor(int a, int b)
{
    while (b != 0) {
        int t = b;
        b = a % b;
        a = t;
    }
    return a;
}

void starneig_init_local_q(int n, size_t ldA, double *A)
{
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            A[i*ldA+j] = 0.0;

    for (int i = 0; i < n; i++)
        A[i*ldA+i] = 1.0;
}

void starneig_small_left_gemm_update(int rbegin, int rend, int cbegin, int cend,
    size_t ldQ, size_t ldA, size_t ldT, double const *Q, double *A, double *T) {

    STARNEIG_SANITY_CHECK_INF(rbegin, rend, cbegin, cend, ldA, A, "A (in)");
    STARNEIG_SANITY_CHECK_INF(0, rend-rbegin, 0, rend-rbegin, ldQ, Q, "Q");

    int m = rend-rbegin;
    int n = cend-cbegin;
    int k = rend-rbegin;

    if (m == 0 || n == 0)
        return;

    starneig_copy_matrix(
        k, n, ldA, ldT, sizeof(double), A+cbegin*ldA+rbegin, T);

    cblas_dgemm(
        CblasColMajor, CblasTrans, CblasNoTrans, m, n, k, 1.0, Q, ldQ, T, ldT,
        0.0, A+cbegin*ldA+rbegin, ldA);

    STARNEIG_SANITY_CHECK_INF(rbegin, rend, cbegin, cend, ldA, A, "A (out)");
}

void starneig_small_right_gemm_update(
    int rbegin, int rend, int cbegin, int cend,
    size_t ldQ, size_t ldA, size_t ldT, double const *Q, double *A, double *T) {

    STARNEIG_SANITY_CHECK_INF(rbegin, rend, cbegin, cend, ldA, A, "A (in)");
    STARNEIG_SANITY_CHECK_INF(0, cend-cbegin, 0, cend-cbegin, ldQ, Q, "Q");

    int m = rend-rbegin;
    int n = cend-cbegin;
    int k = cend-cbegin;

    if (m == 0 || n == 0)
        return;

    starneig_copy_matrix(
        m, k, ldA, ldT, sizeof(double), A+cbegin*ldA+rbegin, T);

    cblas_dgemm(
        CblasColMajor, CblasNoTrans, CblasNoTrans, m, n, k, 1.0, T, ldT, Q, ldQ,
        0.0, A+cbegin*ldA+rbegin, ldA);

    STARNEIG_SANITY_CHECK_INF(rbegin, rend, cbegin, cend, ldA, A, "A (out)");
}

void starneig_small_gemm_updates(
    int begin, int end, int n, size_t ldlQ, size_t ldlZ, size_t ldQ, size_t ldZ,
    size_t ldA, size_t ldB, size_t ldhT, size_t ldvT, double const *lQ,
    double const *lZ, double *Q, double *Z, double *A, double *B,
    double *hT, double *vT)
{
    // apply the local transformation matrices lQ and lZ to Q and Z
    if (Q != NULL)
        starneig_small_right_gemm_update(
            0, n, begin, end, ldlQ, ldQ, ldvT, lQ, Q, vT);
    if (Z != NULL && Z != Q)
        starneig_small_right_gemm_update(
            0, n, begin, end, ldlZ, ldZ, ldvT, lZ, Z, vT);

    // apply the local transformation matrices lQ and lZ to A
    if (A != NULL) {
        starneig_small_right_gemm_update(
            0, begin, begin, end, ldlZ, ldA, ldvT, lZ, A, vT);
        starneig_small_left_gemm_update(
            begin, end, end, n, ldlQ, ldA, ldhT, lQ, A, hT);
    }

    // apply the local transformation matrices lQ and lZ to B
    if (B != NULL) {
        starneig_small_right_gemm_update(
            0, begin, begin, end, ldlZ, ldB, ldvT, lZ, B, vT);
        starneig_small_left_gemm_update(
            begin, end, end, n, ldlQ, ldB, ldhT, lQ, B, hT);
    }
}

void starneig_compute_complex_eigenvalue(
    int ldA, int ldB, double const *A, double const *B,
    double *real1, double *imag1, double *real2, double *imag2,
    double *beta1, double *beta2)
{
    if (B != NULL) {
        extern void dlag2_(double const *, int const *, double const *,
            int const *, double const *, double*, double *, double *, double *,
            double *);

        extern double dlamch_(char const *);
        const double safmin = dlamch_("S");

        double _real1, _real2, _beta1, _beta2, wi;
        dlag2_(
            A, &ldA, B, &ldB, &safmin, &_beta1, &_beta2, &_real1, &_real2, &wi);

        if (beta1) {
            *real1 = _real1;
            *real2 = _real2;
            *imag1 = wi;
            *imag2 = -wi;
            *beta1 = _beta1;
            *beta2 = _beta2;
        }
        else {
            *real1 = _real1/_beta1;
            *real2 = _real2/_beta2;
            *imag1 = wi/_beta1;
            *imag2 = -wi/_beta2;
        }
    }
    else {
        extern void dlanv2_(
            double *, double *, double *, double *, double *,
            double *, double *, double *, double *, double *);

        double a[] = { A[0], A[1], A[ldA], A[ldA+1] };
        double cs, ss;

        dlanv2_(
            &a[0], &a[2], &a[1], &a[3], real1, imag1, real2, imag2, &cs, &ss);

        if (beta1) {
            *beta1 = 1.0;
            *beta2 = 1.0;
        }
    }
}

max_stars_repo_name: NLAFET/StarNEig

max_issues_repo_name: NLAFET/StarNEig

max_forks_repo_name: NLAFET/StarNEig


text:
#ifndef DEGNOME
#define DEGNOME

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

//Degnomedia Rogerus
typedef struct Degnome Degnome;
struct Degnome {
	double* dna_array;
	double hat_size;

};

Degnome* Degnome_new(void);
void Degnome_mate(Degnome* location, Degnome* p1, Degnome* p2, gsl_rng* rng,
	int mutation_rate, int mutation_effect, int crossover_rate);
void Degnome_free(Degnome* q);

int chrom_size;

#endif

max_stars_repo_name: Darokrithia/PolygenSim

max_issues_repo_name: Darokrithia/PopGenSim

max_forks_repo_name: Darokrithia/PolygenSim


text:
/**
 * \author Sylvain Marsat, University of Maryland - NASA GSFC
 *
 * \brief C header for the initialization of the instrumental noise for LIGO/VIRGO detectors.
 *
 *
 */

#ifndef _LLVNOISE_H
#define _LLVNOISE_H

#define _XOPEN_SOURCE 500

#ifdef __GNUC__
#define UNUSED __attribute__ ((unused))
#else
#define UNUSED
#endif

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>
#include <time.h>
#include <unistd.h>
#include <getopt.h>
#include <stdbool.h>
#include <string.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_complex.h>

#include "constants.h"

/************************************************************************/
/****** Global variables storing min and max f for the noise PSD  *******/

extern double __LLVSimFD_LHONoise_fLow;
extern double __LLVSimFD_LHONoise_fHigh;
extern double __LLVSimFD_LLONoise_fLow;
extern double __LLVSimFD_LLONoise_fHigh;
extern double __LLVSimFD_VIRGONoise_fLow;
extern double __LLVSimFD_VIRGONoise_fHigh;

/**************************************************************************/
/****** Prototypes: functions loading and evaluating the noise PSD  *******/

/* Function parsing the environment variable $LLV_NOISE_DATA_PATH and trying to run LLVSimFD_Noise_Init in each */
int LLVSimFD_Noise_Init_ParsePath(void);
/* Function loading the noise data from a directory */
int LLVSimFD_Noise_Init(const char dir[]);

/* The noise functions themselves */
double NoiseSnLHO(const double f);
double NoiseSnLLO(const double f);
double NoiseSnVIRGO(const double f);

#if 0
{ /* so that editors will match succeeding brace */
#elif defined(__cplusplus)
}
#endif

#endif /* _LLVNOISE_H */