text:
#include <math.h>
#include <stdio.h>
#include <stddef.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_linalg.h>
#include <getopt.h>
#include <string.h>

#include "libhdr"
#include "Library_DFE_v1.8.h"

/******************************************************************************/
int get_average_egf_vecs_1_and_2(double s,int n1, int n2, double t2_real,
				 int t2_lower, int t2_upper, double *egf_vec,
				 char *buffer_p1_n2_t2_lower,
				 char *buffer_p2_n2_t2_lower,
				 char *buffer_p1_n2_t2_upper,
				 char *buffer_p2_n2_t2_upper)
{
  static double egf_vec1_lower[maxnd+1], egf_vec2_lower[maxnd+1],
    egf_vec1_upper[maxnd+1], egf_vec2_upper[maxnd+1],
    egf_vec1[maxnd+1], egf_vec2[maxnd+1];
  double total_density;
  int i, res;
  int n1d=2*n1;
  int n2d=2*n2;

  if (s<-1)
    /*
       The contribution of strongly selected mutations
       is assumed to be inversely proportional to s, as
       expected from mutation selection balance, with a
       contribution from phase 2 mutations only.
    */
    {
      for (i=0; i<=2*n2; i++)
	{
	  egf_vec1[i] = 0;
	  egf_vec2[i] = 0;
	}
      egf_vec2[1] = -2/s;
    }
  else
    {
      res =  get_binary_egf_vec(buffer_p1_n2_t2_lower, n2, t2_lower,
				s, egf_vec1_lower);
      if (res==0) return 0;
      res =  get_binary_egf_vec(buffer_p2_n2_t2_lower, n2, t2_lower,
				s, egf_vec2_lower);
      if (res==0) return 0;
      res =  get_binary_egf_vec(buffer_p1_n2_t2_upper, n2, t2_upper,
				s, egf_vec1_upper);
      if (res==0) return 0;
      res =  get_binary_egf_vec(buffer_p2_n2_t2_upper, n2, t2_upper,
				s, egf_vec2_upper);
      if (res==0) return 0;

      compute_weighted_average_egf_vec(t2_real, t2_lower, t2_upper,
				       egf_vec1_lower, egf_vec1_upper, egf_vec1, 2*n2);
      //   monitorinput();
      compute_weighted_average_egf_vec(t2_real, t2_lower, t2_upper,
				       egf_vec2_lower, egf_vec2_upper, egf_vec2, 2*n2);
    }

  total_density = compute_weighted_average(egf_vec, egf_vec1, egf_vec2, n1d, 2*n2);
  if (s==0)
    {
      total_density_s0 = total_density;
      //      printf("s=0: total_density %lf\n", total_density);
    }

  return (1);
}

/******************************************************************************/
get_upper_lower_int(double parm, int *lower, int *upper, int n_int_evaluated, int *int_evaluated_vec)
{
  int i;
  *lower = undefined_int;
  *upper = undefined_int;
  //   printf("get_upper_lower_int parm %lf n_int_evaluated %d\n",
  //      parm, n_int_evaluated);
  for (i=0; i<n_int_evaluated; i++)
    {
      //      printf("int_evaluated_vec[%d] %d\n", i, int_evaluated_vec[i]);
      if ((double)int_evaluated_vec[i] >= parm)
	{
	  *upper = int_evaluated_vec[i];
	  if (i != 0) *lower = int_evaluated_vec[i-1];
	  return;
	}
    }
}
/******************************************************************************/
get_upper_lower_double(double parm, double*lower, double *upper,
		       int n_int_evaluated, double *double_evaluated_vec)
{
  int i;
  *lower = undefined;
  *upper = undefined;
  //   printf("get_upper_lower_int parm %lf n_int_evaluated %d\n",
  //      parm, n_int_evaluated);
  for (i=0; i<n_int_evaluated; i++)
    {
      //      printf("int_evaluated_vec[%d] %d\n", i, int_evaluated_vec[i]);
      if (double_evaluated_vec[i] >= parm)
	{
	  *upper =double_evaluated_vec[i];
	  if (i != 0) *lower = double_evaluated_vec[i-1];
	  return;
	}
    }
}
/******************************************************************************/
get_s_ranges()
{
  int i, stat, ind, i1;
  double d1, d2, d3;
  FILE *inf, *fopen();
  inf = openforreadsilent(s_range_file);
  stat = fscanf(inf, "%d", &nranges);
  if (stat!=1) gabort("get_s_ranges: input error 1", 0);
  if (nranges >max_s_ranges)
    gabort("get_s_ranges: nranges >max_s_ranges", nranges);
  for (i=1; i<=nranges; i++)
    {
      stat = fscanf(inf, "%d %lf %lf %lf %d", &ind, &d1, &d2, &d3, &i1);
      if (stat!=5) gabort("get_s_ranges: input error 2", 0);
      if (ind!=i) gabort("get_s_ranges: ind!=i", i);
      s_ranges[i].s_lower = d1;
      s_ranges[i].s_upper = d2;
      s_ranges[i].s_step = d3;
      s_ranges[i].s_lower_index = i1;
    }

}
/******************************************************************************/
get_int_evaluated_vec(int *int_evaluated_vec, int *n_int_evaluated,
		      int *int_lower, int *int_step,
		      char *int_evaluated_vec_file)
{
  int i, stat;
  FILE *inf, *fopen();
  int int_input, first_int, second_int, int_upper;
  inf = openforreadsilent(int_evaluated_vec_file);
  i = 0;
  for (;;)
    {
      stat = fscanf(inf, "%d", &int_input);
      if (stat==EOF) break;
      if (stat!=1) gabort("get_int_evaluated_vec: read error 1", stat);
      if (*n_int_evaluated >= max_n_int_evaluated)
	gabort("Value of n_int_evaluated too large: get_int_evaluated_vec",
	       get_int_evaluated_vec);
      int_evaluated_vec[*n_int_evaluated] = int_input;
      //      printf("int_evaluated_vec[%d] %d\n",
      //         *n_int_evaluated, int_evaluated_vec[*n_int_evaluated]);
      (*n_int_evaluated)++;
      if (i==0)
	{
	  first_int = int_input;
	  *int_lower = first_int;
	}
      if (i==1)
	{
	  second_int = int_input;
	}
      int_upper = int_input;
      i++;
    }
  fclose (inf);
  *int_step = second_int - first_int;
  //printf("int_lower %d int_upper %d int_step %d\n",
  //  *int_lower, int_upper, *int_step);

}
/******************************************************************************/
scale_vector(double *vec, int n2, double total_density_s0)
{
  int i;
  double tot = 0;
  //   printf("scale_vector: n2 %d total_density_s0 %lf\n", n2, total_density_s0);
  //   monitorinput();
  for (i=1; i<n2; i++)
    {
      vec[i] /= total_density_s0;
      tot += vec[i];
    }
  vec[0] = 1.0 - tot;        // Put all "missing" density in the lost mutant class
}
/******************************************************************************/
int get_binary_egf_vec(char *buffer, int n, int t2, double s, double *egf_vec)
{
  int i, j, file_size_bytes, n_read, t2_read, offset, tab_num, tab_size,
    size_of_float;
  float s_read, freq_read;
  char *ptr, ch;
  // Each expected gene freq table has two integers (n, t2), a float (s), then
  //   2*n -1 floats (frequencies)
  tab_size = (2*sizeof(int)+sizeof(float) + sizeof(float)*(2*n -1));
  file_size_bytes = n_s_evaluated_file*tab_size;
  //   printf("tab_size %d file_size_bytes %d\n", tab_size, file_size_bytes);
  tab_num =  compute_s_tab_num(s);
  //   printf("tab_num assuming ranges %d\n", tab_num);
  //   monitorinput();
  if (tab_num>=n_s_evaluated_file) return 0;
  //   printf("s %lf s_lower %lf s_step %lf (s - s_lower)/s_step %lf tab_num %d\n",
  //        s, s_lower, s_step, (s - s_lower)/s_step, tab_num);
  //   monitorinput();
  offset = tab_num*tab_size;
  assign_int_from_read_buff(&n_read, buffer, offset);
  //   printf("n_read %d\n", n_read);
  if (n_read!=n)
    {
      // printf("n %d n_read %d\n", n, n_read);
      gabort("get_binary_egf_vec: read error 1", n_read);
    }
  offset += sizeof(int);
  assign_int_from_read_buff(&t2_read, buffer, offset);
  //printf("t2 %d t2_read %d\n", t2, t2_read);
  if (t2_read!=t2) gabort("get_binary_egf_vec: read error 2", t2_read);
  offset += sizeof(int);
  s_read = assign_float_from_read_buff(buffer, offset);
  //   printf("s_read %f\n", s_read);
  size_of_float = sizeof(float);
  for (i=1; i<2*n; i++)
    {
      offset += size_of_float;
      //    assign_float_from_read_buff - in line code for speed
      ptr = &freq_read;
      for (j=0; j<size_of_float; j++)
	{
	  //      printf("assign_float_from_read_buff: ch %d\n", ch);
	  *(ptr + j) = buffer[offset + j];
	}
      //      printf("freq_read %f\n", freq_read);
      egf_vec[i] = (double)freq_read;
    }
  //   dumpvector(egf_vec, 1, 2*n -1, "egf_vec:");
  //   monitorinput();
  return 1;
}
/******************************************************************************/
fold_vector_double(double *vec, int counter)
{
  int i, j;
  j = 1;
  for (i=counter - 1; i > counter/2; i--)
    {
      vec[j] += vec[i];
      vec[i] = 0;
      j++;

    }
}
/******************************************************************************/
fold_vector_int(int *vec, int counter)
{
  int i, j;
  j = 1;
  for (i=counter - 1; i > counter/2; i--)
    {
      vec[j] += vec[i];
      vec[i] = 0;
      j++;

    }
}
/******************************************************************************/
int nearest_n2_ind(int n2_real)
{
  int i, ind = undefined_int;
  double min = undefined, diff, res;
  for (i=0; i<n_n2_evaluated; i++)
    {
      diff = n2_evaluated_vec[i] - n2_real;
      if (diff < 0) diff = -diff;

      if (min==undefined)
	{
	  min = diff;
	  ind = i;
	}
      else
	{
	  if (min > diff)
	    {
	      min = diff;
	      ind = i;
	    }
	}
    }
  if (ind == undefined_int) gabort("nearest_n2: ind == undefined_int", ind);
  res = n2_evaluated_vec[ind];;
  //   printf("nearest_n2: n2_real %lf result %lf\n", n2_real, res);
  //   monitorinput();
  return ind;
}
/******************************************************************************/
read_phase1_phase2_file_into_buffer(int n1, int phase, int n2, int t2,
				    char *buffer,
				    int file_size_bytes)
{
  static char n1_str[max_str_len], n2_str[max_str_len], t2_str[max_str_len];
  char file_str[max_str_len];
  int n2d;
  FILE *inf, *fopen();

  //   printf("read_phase1_phase2_file_into_buffer: phase %d n2 %d t2 %d\n",
  //      phase, n2, t2);
  n2d = 2*n2;
  int_to_string(n1, n1_str, max_str_len);

  int_to_string(n2, n2_str, max_str_len);

  int_to_string(t2, t2_str, max_str_len);

  if (strlen(n1_str) + strlen(n2_str) + strlen(t2_str) >= max_str_len)
    gabort("read_phase1_phase2_file_into_buffer: string too long", max_str_len);
  if (phase==1)
    {
      strcpy(file_str, phase_1_dir);
      strcat(file_str, "/n1:");
      strcat(file_str, n1_str);
      strcat(file_str, ":");
    }
  else
    {
      strcpy(file_str, phase_2_dir);
      strcat(file_str, "/");
    }
  strcat(file_str, "n2:");
  strcat(file_str, n2_str);
  strcat(file_str, ":t2:");
  strcat(file_str, t2_str);

  inf = openforreadsilent(file_str);
  fread(buffer, 1, file_size_bytes, inf);
  fclose(inf);
}
/******************************************************************************/
read_const_pop_file_into_buffer(int n1, char *buffer, int file_size_bytes)
{
  char file_str[max_str_len];
  FILE *inf, *fopen();

  //   printf("read_const_pop_file_into_buffer: n1 %d file_size_bytes %d\n",
  //      n1, file_size_bytes);
  strcpy(file_str, const_pop_dir);
  strcat(file_str, const_pop_file);
  //   printf("read_const_pop_file_into_buffer: file_str %s\n", file_str);
  //   monitorinput();
  inf = openforreadsilent(file_str);
  fread(buffer, 1, file_size_bytes, inf);
  fclose(inf);
}
/******************************************************************************/
compute_weighted_average_egf_vec(double t2_real, int t2_lower, int t2_upper,
				 double *egf_vec_lower, double *egf_vec_upper,
				 double *egf_vec, int n2d)
{
  int i;
  for (i=0; i<=n2d; i++)
    {
      egf_vec[i] = egf_vec_lower[i] + (egf_vec_upper[i] - egf_vec_lower[i])*(t2_real - (double)t2_lower)/(double)(t2_upper - t2_lower);
    }
}
/******************************************************************************/
double compute_weighted_average(double *weighted_vec, double *phase1_mut_vec,
				double *phase2_mut_vec, int n1d, int n2d)
{
  int i;
  double total = 0;
  for (i=1; i<n2d; i++)
    {
      weighted_vec[i]= (n1d*phase1_mut_vec[i] + n2d*phase2_mut_vec[i])/
	(n1d + n2d);
      total += weighted_vec[i];
    }
  //   dumpvector(weighted_vec, 0, n2d, "weighted average vector:"); monitorinput();
  return(total);
}
/******************************************************************************/
int get_const_pop_egf_vec(double s, int n1, double *egf_vec, char *buffer_const_pop,
			  int assign_tds0_mode)
{
  double total_density;
  int i, res;
  if (s<-1)          // The contribution of strongly selected mutations
    // is assumed to be inversely proportional to s, as
    // expected from mutation selection balance, with a
    // contribution from phase 2 mutations only.
    {
      for (i=0; i<=2*n1; i++)
	{
	  egf_vec[i] = 0;
	}
      egf_vec[1] = -2/s;
      //      printf("get_const_pop_egf_vec s %lf\n", s);
      //      dumpvector(egf_vec, 0, 2*n1, "get_const_pop_egf_vec:");
      //      monitorinput();
    }
  else
    {
      res =  get_binary_egf_vec(buffer_const_pop, n1, 0, s, egf_vec);
      if (res==0) return 0;
      //      dumpvector(egf_vec, 0, 2*n1, "get_const_pop_egf_vec:");
      //      monitorinput();
    }
  total_density = compute_weighted_average(egf_vec, egf_vec, egf_vec, 2*n1, 2*n1);
  if (assign_tds0_mode)
    {
      total_density_s0 = total_density;
    }
  //   printf("assign_tds0_mode %d: total_density_s0 %lf\n", assign_tds0_mode, total_density_s0);
  return (1);
}
/******************************************************************************/
assign_int_from_read_buff(int *res, char *buffer, int start)
{
  int temp, i;
  char *ptr, ch;
  ptr = &temp;
  for (i=0; i<sizeof(int); i++)
    {
      ch = buffer[start+i];

      *(ptr + i) = ch;
    }

  *res = temp;
}
/******************************************************************************/
float assign_float_from_read_buff(char *buffer, int start)
{
  float temp;
  int i;
  char *ptr, ch;
  ptr = &temp;
  for (i=0; i<sizeof(float); i++)
    {
      ch = buffer[start+i];
      //      printf("assign_float_from_read_buff: ch %d\n", ch);
      *(ptr + i) = ch;
    }
  //   printf("assign_float_from_read_buff: temp %f\n", temp);
  return temp;
}
/******************************************************************************/
int compute_s_tab_num(double s)
{
  double s_lower, s_upper, s_step;
  int i, tab_num;
  for (i=1; i<=nranges; i++)
    {
      s_lower = s_ranges[i].s_lower;
      s_upper = s_ranges[i].s_upper;
      if ((s >= s_lower)&&(s <= s_upper))
	{
	  s_step = s_ranges[i].s_step;
	  tab_num = (int)(floor((s - s_lower)/s_step + 0.0000000001));
	  return(tab_num + s_ranges[i].s_lower_index - 1);
	}
    }
  gabort("compute_s_tab_num: failure to find tab_num", 0);
}
/******************************************************************************/
int compute_file_size_bytes(int n)
{
  int tab_size, file_size_bytes;
  tab_size = (2*sizeof(int)+sizeof(float) + sizeof(float)*(2*n - 1));
  file_size_bytes = n_s_evaluated_file*tab_size;
  return (file_size_bytes);
}
/******************************************************************************/
set_up_file_name(int n1, char *froot, char *file_name)
{
  int len;
  static char n1_str[max_str_len];
  int_to_string(n1, n1_str, max_str_len);
  //   printf("set_up_file_name n1 %d froot %s n1_str %s\n", n1, froot, n1_str);
  //   monitorinput();
  len = strlen(data_path);
  len += strlen(n1_str) + 1;
  len += strlen(froot);
  if (len >= max_str_len) gabort("set_up_file_name: string too long (1)", len);
  strcpy(file_name, data_path);
  strcat(file_name, "n1_");
  strcat(file_name, n1_str);
  strcat(file_name, "/");
  strcat(file_name, froot);
  //   printf("set_up_file_name n1 file_name %s\n", file_name);
  //   monitorinput();
}
/******************************************************************************/
get_data_path(char *data_path)
{
  FILE *inf, *fopen();
  int stat;
  inf = openforread("directory_config.dat");
  stat = fscanf(inf, "%s", data_path);
  if (stat!=1) gabort("get_data_path: read error", 1);
  //   printf("get_data_path: data_path %s\n", data_path); monitorinput();
}
/******************************************************************************/
get_s_evaluated_vec(double *s_evaluated_vec, int *n_s_evaluated,
		    int *n_s_evaluated_file, char *s_evaluated_vec_file)
{
  int stat;
  FILE *inf, *fopen();
  double s;
  *n_s_evaluated = 0;
  inf = openforreadsilent(s_evaluated_vec_file);
  s = -100.0;                                 // Set up large |s| values specially
  while (s < -1)
    {
      s_evaluated_vec[*n_s_evaluated] = s;
      (*n_s_evaluated)++;
      if (s > -20) s = s + 1.0; else s = s + 5.0;
    }
  *n_s_evaluated_file = 0;
  for (;;)
    {
      stat = fscanf(inf, "%lf", &s);
      if (stat==EOF) break;
      if (stat!=1) gabort("get_s_limits: read error 1", stat);
      if (*n_s_evaluated >= max_n_s_evaluated)
	gabort("Value of n_s_evaluated too large: get_s_evaluated_vec",
	       get_s_evaluated_vec);
      s_evaluated_vec[*n_s_evaluated] = s;
      (*n_s_evaluated)++;
      (*n_s_evaluated_file)++;
    }
  fclose (inf);
  //   printf("n_s_evaluated %d n_s_evaluated_file %d\n", *n_s_evaluated,
  //      *n_s_evaluated_file);
  //   monitorinput();
}
/******************************************************************************/
get_lower_upper(int ind, double *s_evaluated_vec, int n_s_evaluated,
		double *lower, double *upper)
{
  double x, temp_r;
  if ((ind<0)||(ind>=n_s_evaluated))
    gabort("get_lower_upper: ind out of range", ind);
  x = s_evaluated_vec[ind];
  if (ind<n_s_evaluated-1)
    *upper = (x + s_evaluated_vec[ind+1])/2.0;
  else
    *upper = x;
  if (ind!=0)
    *lower = (x + s_evaluated_vec[ind-1])/2.0;
  else
    *lower = x;
  if (x <=0.0)
    {
      temp_r = *lower;
      *lower = *upper;
      *upper = temp_r;
      if (*lower!=0.0) *lower = -*lower;
      if (*upper!=0.0) *upper = -*upper;
    }
}
/******************************************************************************/
dumpvector(double *v, int min, int max, char *s)
{
  int i, control = 0;
  //   printf("\n");
  printf("%s", s); printf("\n");
  for (i = min; i<max; i++)
    {
      printf("%3d %10.8lf ", i, v[i]);
      control++;
      if (control==5)
	{
	  control = 0;
	  printf("\n");
	}
    }
  printf("\n");
}
/******************************************************************************/
dumpmatrix(double **m, int s1,
	   int rows, int s2, int cols, char *s)
{
  int i, j;
  printf("\n%s\n", s);
  for (i = s1; i<=rows; i++)
    {
      for (j = s2; j<=cols; j++)
	{
	  printf(" %2.3f", m[i][j]);
	}
      printf("\n");
    }
  printf("\n");
}
/******************************************************************************/
double selfn(double q, double s, double h)
{
  double num, denom;
  num = s*q*(1-q)*(q + h*(1-2*q));
  denom = 1 + 2.0*h*s*q*(1-q) + s*q*q;
  if (denom== 0) denom = .00000000000001;
  return( q + (num/denom));
}
/******************************************************************************/
setuprow(double **a, double p, int row, int n)
{

  double z;         /*probability of zero failures*/
  double x, y, temp1, temp2;
  int j;

  temp1 = 1.0 - p;
  if (temp1 == 0.0)  temp1 = .000000000001;            /*prevent zero divide*/
  temp1 = p/temp1;       /*First prevent negative log*/
  if (temp1 <=0)   temp1 = .000000000001;
  x = log(temp1);    /*used for later multiplication of series*/
  temp2 = 1.0-p;
  if (temp2<=0)   temp2 = .000000000001;
  z = (double)n*log(temp2);
  if (z > 150)   z = 150;
  a[row][0] = exp(z);
  y = a[row][0];        /*add up for later*/
  for (j = 1;j<=n; j++)
    {

      z = z + log(((double)n + 1.0 - (double)j)/(double)j) + x;
      /*             cumulative binomial coeff.*/
      a[row][j] = exp(z);
      y = y + a[row][j];          /*add up*/
    }
  if (y == 0.0)   y = .0000000000001;
  for (j=0; j<=n; j++) a[row][j] = a[row][j]/y;
}
/******************************************************************************/
setupmatrix(double **a, int n1, int n2, double s, double h)
{
  int i;
  double p;

  for (i = 1; i<=n1 - 1; i++)             /*do variable rows*/
    {
      p = (double)i/(double)n1;
      p = selfn(p, s,h);
      setuprow(a, p, i, n2);
    }
  for (i = 0; i<=n2; i++)
    {              /*do constant edges*/
      a[0][i] = 0;
      a[n1][i] = 0;
    }
  a[0][0] = 1.0;                          /*do absorbing corners*/
  a[n1][n2] = 1.0;

}
/******************************************************************************/
matrixinvert(int N, double **a,double *inverter)
{

  int i, j, lotkin_signum;
  gsl_matrix  *lotkin_a,*lotkin_inv;
  gsl_vector *x, *lotkin_b, *lotkin_x;

  gsl_permutation *lotkin_perm;

  /* allocate a, x, b */
  lotkin_a = gsl_matrix_alloc(N-1, N-1);
  lotkin_inv = gsl_matrix_alloc(N-1, N-1);
  x = gsl_vector_alloc(N-1);
  lotkin_b = gsl_vector_alloc(N-1);

  lotkin_x = gsl_vector_alloc(N-1);

  gsl_matrix *identity=gsl_matrix_alloc(N-1, N-1);

  for (i = 0; i<=N-2; i++)
    {
      for (j = 0; j<=N-2; j++)
	{
          gsl_matrix_set( identity, i, i,1);
	}
      ;
    }

  for (i = 0; i<=N-2; i++)
    {
      for (j = 0; j<=N-2; j++)
	{
	  gsl_matrix_set(lotkin_a,i,j,gsl_matrix_get(identity,i,j)-a[i+1][j+1]);
	}

    }

  /* LU decomposition and forward&backward substition */

  //gsl_blas_dgemv(CblasNoTrans, 1.0, lotkin_a, x, 0.0, lotkin_b);

  lotkin_perm = gsl_permutation_alloc(N-1);
  gsl_linalg_LU_decomp(lotkin_a, lotkin_perm, &lotkin_signum);
  gsl_linalg_LU_solve(lotkin_a, lotkin_perm, lotkin_b, lotkin_x);
  gsl_linalg_LU_invert (lotkin_a,lotkin_perm, lotkin_inv);

  /*apparently the inversion adds +1 to the first element of the vector
    so I substract it*/

  gsl_matrix_set(lotkin_inv,0,0,gsl_matrix_get(lotkin_inv,0,0)-1);

  //printf("\n%s\n", "Equilibrium");

  for (i = 0; i<1; i++)
    {
      for (j = 0; j<=N-2; j++)
	{

	  // printf("%d: %2.3f  ",j,gsl_matrix_get(lotkin_inv, i, j));
	}
      //printf("\n");
    }

  for (i=0;i<=N-2;i++){

    inverter[i+1]=gsl_matrix_get(lotkin_inv,0,i);

  }

  gsl_matrix_free(lotkin_a);
  gsl_matrix_free(lotkin_inv);
  gsl_vector_free(x);
  gsl_vector_free(lotkin_b);
  gsl_vector_free(lotkin_x);
  gsl_matrix_free(identity);
  gsl_permutation_free(lotkin_perm);
}
/******************************************************************************/
tmiterate(double **a,double *mut1,
	  int t, int n1, int n2,int decay,double *sumf)
{

  int k, i, j = 0;
  double z;

  double * steadystatefreq = (double*) calloc (maxnd+1, sizeof(double));
  double * mut2 = (double*) calloc (maxnd+1, sizeof(double));

  for (k =1 ;k<=t; k++)
    {                 /*t iterations*/

      /* Increment steadystatefreq. distribution*/
      for (i = 0; i<=n1; i++)
	{
	  steadystatefreq[i] = steadystatefreq[i] + mut1[i];
	}

      /* Perform matrix multiplication*/
      for (i = 0; i<=n2; i++)
	{
	  z = 0;
	  for (j = 0; j<=n1; j++)
	    {
	      z = z + a[j][i]*mut1[j];
	    }
	  mut2[i] = z;
	}

      /* Copy result of multiplication*/
      for (i=0; i<=n2; i++) {mut1[i] = mut2[i];

	if (decay==0){
	  sumf[i]+=mut1[i];}
	if(decay==1){sumf[i]=mut1[i];}
      }

    }

  free(steadystatefreq);
  free(mut2);
}
/******************************************************************************/
tmiterate2(double **a,double *mut1,
	   int t, int n1, int n2,int decay,double **a1)
{

  double *  steadystatefreq = (double*) calloc (maxnd+1, sizeof(double));
  double * mut2 = (double*) calloc (maxnd+1, sizeof(double));

  int k, i, j = 0;
  double z;

  for (k = 1; k<=t; k++){                 /*t iterations*/

    /* Increment steadystatefreq. distribution*/
    for (i = 0; i<=n1; i++)
      {
	steadystatefreq[i] = steadystatefreq[i] + mut1[i];
      }

    /* Perform matrix multiplication*/
    for (i = 0; i<=n2; i++)
      {
	z = 0;
	for (j = 0; j<=n1; j++)
	  {
            z = z + a[j][i]*mut1[j];
	  }
	mut2[i] = z;
      }

    /* Copy result of multiplication*/
    for (i=0; i<=n2; i++) {mut1[i] = mut2[i];

      if (decay==0){

	a1[k][i]+=mut1[i];
	if (k>1){a1[k][i]+=a1[k-1][i];}
      }

      if(decay==1){

	a1[k][i]=mut1[i];
	//if (k>0){ a1[k][i]+=a1[k-1][i];}
      }
    }

    //dumpvector(a1[k],0,n2,"matrix");
  }

  free(steadystatefreq);
  free(mut2);

}
/******************************************************************************/
eqf_using_matrix_inversion(int n1,double s,double **a,double *egf_out)
{
  setupmatrix(a, n1, n1, s, H);
  matrixinvert(n1,a,egf_out);
}
/******************************************************************************/
vector_average(int n1,int n2,double *fv1,double *fv2, double *fv_out)
{
  int i;
  for (i=0;i<n2;i++)
    {

      fv_out[i]=((n1*fv1[i])+(n2*fv2[i]))/(n1+n2);

    }
}
/******************************************************************************/
vector_s_average(int n,double P1,double *fv1,double *fv2, double *fv_out)
{
  int i;
  for (i=0;i<n;i++)
    {

      fv_out[i]=(P1*fv1[i])+((1-P1)*fv2[i]);

    }
}
/******************************************************************************/
output_sfs_to_file_thanasis_format(int n, int *sfs1,int *sfs2,char *sfs_filename)
{

  FILE *file= fopen(sfs_filename, "w" );

  int i;

  //neutral sites first

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

    fprintf(file,"%d ",sfs1[i]);

  }

  fprintf(file,"\n");

  //selected sites second

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

    fprintf(file,"%d ",sfs2[i]);

  }

  fclose(file);

}
/******************************************************************************/
output_sfs_to_file_peter_format(int n,int sample1, int sample2, int *sfs1,int *sfs2,char *sfs_filename)
{

  FILE *file= fopen(sfs_filename, "w" );
  int i;
  fprintf(file,"1\n");
  fprintf(file,"%d\n",n);

  //selected sites first
  fprintf(file,"%d\n",sample2);

  for (i=1;i<n;i++){

    fprintf(file,"%d ",sfs2[i]);

  }

  fprintf(file,"\n");
  //neutral sites second

  fprintf(file,"%d\n",sample1);
  for (i=1;i<n;i++){

    fprintf(file,"%d ",sfs1[i]);

  }

  fclose(file);

}
/******************************************************************************/
output_sfs_to_file_peter_format2(int n,int sample1, int sample2, int *sfs1,int *sfs2,char *sfs_filename)
{

  FILE *file= fopen(sfs_filename, "w" );
  int i;
  fprintf(file,"name\n");

  fprintf(file,"0 0\n0 0\n");
  fprintf(file,"%d\n",n);
  for (i=0;i<n;i++){

    fprintf(file,"%d ",sfs1[i]);

  }

  fprintf(file,"0\n");
  //fprintf(file,"%d\n",sample1);
  for (i=0;i<n;i++){

    fprintf(file,"%d ",sfs2[i]);

  }
fprintf(file,"0");
  fclose(file);

}
/******************************************************************************/

get_sfs(int alleles,int *par,char *sfs_filename)
{

  int i;
  int read;
  // char line[80];
  FILE *file= fopen(sfs_filename, "r" );

  for (i =0 ;i<alleles; i++)
    {
      fscanf (file, "%d", &read);
      par[i]=read;
    }

  fclose(file);

}
/******************************************************************************/
get_sfs_peter1(int *total_neu,int *total_sel,int *nalleles,int *sfs_sel,int *sfs_neu,char *sfs_filename)
{
  int i;
  int gene_number,sum;

  FILE *file= fopen(sfs_filename, "r" );
  fscanf (file, "%d", &gene_number);
  fscanf (file, "%d", nalleles);

  //printf("%d %d\n",read_sfs,*nalleles);

  fscanf (file, "%d", total_sel);
  sum=0;
  for (i =1 ;i<*nalleles; i++)
    {
      fscanf (file, "%d", &sfs_sel[i]);
      //printf("%d\n",sfs_sel[i]);
      sum+=sfs_sel[i];
    }
  sfs_sel[0]=*total_sel-sum;
  sfs_sel[*nalleles]=0;
  fscanf (file, "%d", total_neu);
  sum=0;
  for (i =1 ;i<*nalleles; i++)
    {
      fscanf (file, "%d", &sfs_neu[i]);
      //printf("%d\n",sfs_neu[i]);
      sum+=sfs_neu[i];
    }
  sfs_neu[0]=*total_neu-sum;
  sfs_neu[*nalleles]=0;
  //printf("%d\n",sfs_neu[*nalleles]);
  //monitorinput();
  fclose(file);
}
/******************************************************************************/
get_sfs_peter2(int *total_neu,int *total_sel,int *nalleles,int *sfs_sel,int *sfs_neu,char *sfs_filename,int *sel_sites,int *sel_diff,int *neu_sites,int *neu_diff)
{

  int i;
  int gene_number,sum;
  char name[100];
  FILE *file= fopen(sfs_filename, "r" );

  fscanf (file, "%s", name);
  fscanf (file, "%d", sel_sites);
  fscanf (file, "%d", sel_diff);
  fscanf (file, "%d", neu_sites);
  fscanf (file, "%d", neu_diff);
  fscanf (file, "%d", nalleles);

  //printf("%s %d %d %d %d %d\n",name,sel_sites,sel_diff,neu_sites,neu_diff,*nalleles);

  fscanf (file, "%d", &sfs_sel[0]);
  //printf("%d\n",sfs_sel[0]);
  sum=0;
  for (i =1 ;i<=*nalleles; i++)
    {
      fscanf (file, "%d", &sfs_sel[i]);
      //printf("%d\n",sfs_sel[i]);
      sum+=sfs_sel[i];
    }

//printf("next\n");
  fscanf (file, "%d", &sfs_neu[0]);
  //printf("%d\n",sfs_neu[0]);
  sum=0;
  for (i =1 ;i<=*nalleles; i++)
    {
      fscanf (file, "%d", &sfs_neu[i]);
      //printf("%d\n",sfs_neu[i]);
      sum+=sfs_neu[i];
    }

  //monitorinput();
  fclose(file);
}

/******************************************************************************/
double sfsfold_f(int nalleles,double *sfs,double *sfs_folded,int index)
{

int i;

      for (i=0;i<=nalleles/2;i++)
	{
	  sfs_folded[index+i]=sfs[i]+sfs[index-i];

	  if(i==nalleles/2&&(!odd(nalleles)))
	  {
	  sfs_folded[index+i]=sfs[i];
	  }

	}

}

/******************************************************************************/
int binarysearchcumulative_from_zero(double *array, int size)
{
  int cur, pre, post;
  double r;
  r = uniform();
  /*   printf("Uniform %f\n", r);*/
  pre = 0;
  post = size;
  cur = (size + 1)/2;
  do
    {
      if (array[cur] < r) pre = cur; else post = cur;
      cur = (pre + post)/2;
      if ((cur==post) || (cur==pre))
	{
	  if (r < array[pre])
	    {
	      /*            printf("Returning pre %d\n", pre);*/
	      return(pre);
	    }
	  else
	    {
	      /*            printf("Returning post %d\n", post);*/
	      return(post);
	    }
	}
    } while (size > 0);
}
/******************************************************************************/
double randomly_select_allele_frequency(double *egf_vec, double *temp1, int nd)
{
  double u, freq;
  int ind;
  ind = binarysearchcumulative_from_zero(egf_vec, nd-1);
  //   printf("ind %d\n", ind);
  (temp1[ind])++;
  freq = (double)ind/(double)nd;
  //   printf("ind %d freq %lf\n", ind, freq); monitorinput();
  return(freq);
}
/******************************************************************************/
set_up_cumulative_allele_freq_vec(int nd, double *egf_vec, double *cum_vec,
				  char *str)
{
  int i;
  for (i=0; i<nd; i++)
    {
      cum_vec[i] = egf_vec[i];
    }
  for (i=1; i<=nd; i++)
    {
      cum_vec[i] += cum_vec[i-1];
    }
}
/******************************************************************************/
generate_simulated_data(double *egf_vec, int nsites, int sample_size, int nd,
			int *sfs)
{
  int i, n_seg;
  double allele_freq;
  struct acc a;
  initacc(&a);
  for (i=1; i<=nsites; i++)
    {
      allele_freq =
	randomly_select_allele_frequency(egf_vec, allele_freq_sampled, nd);
      n_seg = bnldev(allele_freq, sample_size);
      accum(&a, (double)n_seg);
      (sfs[n_seg])++;
    }
}

/******************************************************************************/
egf_scaling_s(int n1,double f0,double *v_s_in, double *v_s_out)
{

  double sumx,sumy;
  int i;
  sumx=0;
  sumy=0;

  for (i=1;i<n1;i++){

    sumx+=v_s_in[i];

  }

  for (i=1;i<n1;i++){
    v_s_out[i]/=sumx;
    sumy+=v_s_out[i];
  }

  v_s_out[0]=(1-sumy);
}
/******************************************************************************/
egf_scaling_f0(int n1,double f0,double *fv)
{

  int i;
  for (i=0;i<n1;i++){
    fv[i]*=(1-f0);
  }
  fv[0]+=f0;
}
/******************************************************************************/
egf_scaling(int N,double f0,double *FV0,double *FVS)
{
  int n2d=N*2;
  egf_scaling_s(n2d,f0,FV0, FVS);
  egf_scaling_f0(n2d,f0,FVS);
}
/******************************************************************************/
binomial_sampling(int n1,int alleles,int sample, double *invy,int *discrete)
{

  int s1,success;
  int i;
  int n1d=2*n1;
  double prob;

  const gsl_rng_type * T;
  T = gsl_rng_taus;
  gsl_rng  *rgen =gsl_rng_alloc(T);

  /*
    based on equilibrium frequency vector (use it as a probability vector),
    I randomly generate numbers from 0 to n1d-1 for sample sites.
  */

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

      gsl_ran_discrete_t *r= gsl_ran_discrete_preproc (n1d,invy);
      s1=gsl_ran_discrete (rgen, r);
      prob=(double)(s1)/n1d;

      success=gsl_ran_binomial(rgen,prob,alleles-1);
      discrete[success]++;

      gsl_ran_discrete_free(r);
    }
  gsl_rng_free (rgen);

}
/******************************************************************************/
int set_up_density_vec_equal_effects(double s, int n_s_evaluated, double *gamma_density_vec)
{
   double s_lower, s_upper, p_upper;
   int i, s_ind;
   s_lower = undefined;
   s_upper = undefined;
//   printf("set_up_gamma_density_vec_equal_effects s %lf\n", s);
//   monitorinput();
   for (i=0; i<n_s_evaluated; i++)
   {
      gamma_density_vec[i] = 0.0;
   }
   for (i=1; i<n_s_evaluated; i++)
   {
//      printf("%d %lf  ", i, s_evaluated_vec[i]);
//      printf("\n");
      if (s_evaluated_vec[i] >= s)
      {
         s_upper = s_evaluated_vec[i];
         s_ind = i;
         if (i != 0 ) s_lower = s_evaluated_vec[i-1];
         break;
      }
   }
   if ((s_lower==undefined)||(s_upper==undefined)) return 0;
   p_upper = (s - s_lower)/(s_upper - s_lower);
//   printf("set_up_gamma_density_vec_equal_effects: s_lower %lf s_upper %lf\n",
//      s_lower, s_upper);
//   printf("p_upper %lf\n", p_upper);
//   monitorinput();
   gamma_density_vec[s_ind] = p_upper;
   if (s_ind==1)
      gamma_density_vec[s_ind] += (1.0 - p_upper);
   else
      gamma_density_vec[s_ind-1] = (1.0 - p_upper);
//   for (i=1; i<n_s_evaluated; i++)
//   {
//      printf("s %lf gamma_density_vec[%d] %lf\n", s_evaluated_vec[i], i,
//         gamma_density_vec[i]);
//   }
//   monitorinput();
return (1);
}
/******************************************************************************/
int set_up_density_vec_step_effects(double alpha, double beta,
				    int n_s_evaluated, double *gamma_density_vec)
{
 if (alpha<0){alpha=-alpha;}
  if (beta<0){beta=-beta;}
  int i,s_ind;
  double lower, upper, area, total_area = 0.0, x, diff, step, inc_x, inc_y, gf;
  double mean;

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

      x = s_evaluated_vec[i];

      get_lower_upper(i, s_evaluated_vec, n_s_evaluated, &lower, &upper);

      /*calculate cumulative uniform distribution for lower and upper*/
      inc_x=gsl_cdf_flat_P(lower,alpha,beta);
      inc_y=gsl_cdf_flat_P(upper,alpha,beta);

      area=inc_y-inc_x;
      gamma_density_vec[i] = area;
      total_area+=area;
      //printf("x %f inc_x %f inc_y %f area %f\n",x,inc_x,inc_y,area);
    }
//printf("alpha:%f beta:%f total_area:%f\n",alpha,beta, total_area);
//monitorinput();
// area =gsl_cdf_flat_P(-s_evaluated_vec[0], alpha, beta);
// gamma_density_vec[0]+=1.0-area;

 //if (gamma_density_vec[0]>0){printf("%f\n",gamma_density_vec[0]);}
  return (1);
}
/******************************************************************************/
double compute_beta_densities(double alpha, double *s_evaluated_vec,
			      int n_s_evaluated,
			      double *gamma_density_vec, double beta)
{
  int i;
  double lower, upper, area, total_area = 0.0, x, diff, step, inc_x, inc_y, gf;

  if (alpha <= 0) return undefined;
  if (beta <= 0) return undefined;

  for (i=0; i<n_s_evaluated; i++)
    {
      x = s_evaluated_vec[i];
      get_lower_upper(i, s_evaluated_vec, n_s_evaluated, &lower, &upper);

      inc_x = gsl_cdf_beta_P (lower, alpha, beta);
      inc_y = gsl_cdf_beta_P (upper, alpha, beta);

      area = (inc_y - inc_x);
      gamma_density_vec[i] = area;

      //printf("x %f inc_x %f inc_y %f area %f\n",x,inc_x,inc_y,area);

    }
   area = gsl_cdf_beta_P(-s_evaluated_vec[0], alpha, beta);
   gamma_density_vec[0]+=1.0-area;

  return (1);
}
/******************************************************************************/
double compute_gamma_densities(double alpha, double *s_evaluated_vec,
			       int n_s_evaluated, double *gamma_density_vec, double beta)
{
  int i;
  double lower, upper, area, total_area = 0.0, x, diff, step, inc_x, inc_y, gf;

  if (alpha <= 0) return undefined;

  for (i=0; i<n_s_evaluated; i++)
    {
      x = s_evaluated_vec[i];
      get_lower_upper(i, s_evaluated_vec, n_s_evaluated, &lower, &upper);

      inc_x = gsl_cdf_gamma_P(lower, beta,1/alpha);
      inc_y = gsl_cdf_gamma_P(upper, beta,1/alpha);

      area = (inc_y - inc_x);
      gamma_density_vec[i] = area;
      total_area+=area;
      //printf("x%f 1/alpha %f beta %f area %f\n",x,1/alpha,beta,area);

    }

  area = gsl_cdf_gamma_P(-s_evaluated_vec[0], beta,1/alpha);
  gamma_density_vec[0]+=1.0-area;
  //printf("x%f inc_x %f inc_y %f area %f\n",x,inc_x,inc_y,area);

  return (1);

}

/******************************************************************************/
double calculate_ne(int n1,int n2,int t)
{
double w1,w2,n_e;
double N1=(double)n1;
double N2=(double)n2;
double t2=(double)t;
/*compute weighted Ne*/
w1=N1*pow((1-1/(2*N2)),t2);
w2=N2*(1-exp(-t2/(2*N2)));

n_e=(N1*w1+N2*w2)/(w1+w2);

return n_e;
}


text:
// --------------------------------------------------------
// Dragon
// Copyright(c) 2017 SeetaTech
// Written by Ting Pan
// --------------------------------------------------------

#ifndef DRAGON_UTILS_MATH_FUNCTIONS_H_
#define DRAGON_UTILS_MATH_FUNCTIONS_H_

#include <float.h>
#include <cstdint>
#include <climits>

#ifdef WITH_BLAS
extern "C" {
#include <cblas.h>
}
#else    // WITH_BLAS
typedef enum CBLAS_TRANSPOSE {CblasNoTrans=111, CblasTrans=112} CBLAS_TRANSPOSE;
#endif

#include "protos/dragon.pb.h"

namespace dragon {

namespace math {

/******************** Level-0 ********************/

template <typename T, class Context>
void Set(const int n, const T alpha, T* x);

template <typename T, class Context>
void RandomUniform(const int n, const float low, const float high, T *x);

template <typename T, class Context>
void RandomNormal(const int n, const float mu, const float sigma, T* x);

template <typename T, class Context>
void RandomTruncatedNormal(const int n,
                           const float mu,
                           const float sigma,
                           const float low,
                           const float high,
                           T* x);

template <typename T, class Context>
void RandomBernoulli(const int n, const float p, uint32_t* x);

/******************** Level-1 ********************/

template <typename T, class Context>
void Add(const int n, const T* a, const T* b, T* y);

template <typename T, class Context>
void Sub(const int n, const T* a, const T* b, T* y);

template <typename T, class Context>
void Mul(const int n, const T* a, const T* b, T* y);

template <typename T, class Context>
void Div(const int n, const T* a, const T* b, T* y);

template <typename T, class Context>
void Clip(const int n, const float low, const float high, T* x);

template <typename T, class Context>
void Exp(const int n, const T* x, T* y);

template <typename T, class Context>
void Log(const int n, const T* x, T* y);

template <typename T, class Context>
void Square(const int n, const T* x, T* y);

template <typename T, class Context>
void Sqrt(const int n, const T* x, T* y);

template <typename T, class Context>
void Pow(const int n, const float alpha, const T* x, T* y);

template <typename T, class Context>
void Inv(const int n, const float numerator, const T* x, T* y);

/******************** Level-2 ********************/

template <typename T, class Context>
void Scal(const int n, const float alpha, T* y);

template <typename T, class Context>
void Scale(const int n, const float alpha, const T* x, T* y);

template <typename T, class Context>
T StridedDot(const int n,
             const T* a,
             const int incx,
             const T* b,
             const int incy);

template <typename T, class Context>
float Dot(const int n, const T* a, const T* b);

template<typename T, class Context>
float ASum(const int n, const T *x);

template<typename T, class Context>
void AddScalar(const int n, const float alpha, T* y);

template<typename T, class Context>
void MulScalar(const int n, const float alpha, T* y);

template<typename T, class Context>
void Axpy(const int n, float alpha, const T* x, T *y);

template<typename T, class Context>
void Axpby(const int n, float alpha, const T* x, float beta, T *y);

/******************** Level-3 ********************/

template <typename T, class Context>
void Gemm(const CBLAS_TRANSPOSE transA,
          const CBLAS_TRANSPOSE transB,
          const int M,
          const int N,
          const int K,
          const float alpha,
          const T* A,
          const T* B,
          const float beta,
          T* C,
          TensorProto_DataType math_type = TensorProto_DataType_FLOAT);

template<typename T, class Context>
void Gemv(const CBLAS_TRANSPOSE transA,
          const int M,
          const int N,
          const float alpha,
          const T* A,
          const T* x,
          const float beta,
          T* y,
          TensorProto_DataType math_type = TensorProto_DataType_FLOAT);

}    // namespace math

}    // namespace dragon

#endif    // DRAGON_UTILS_MATH_FUNCTIONS_H_

ext: h

max_stars_repo_name: neopenx/Dragon

max_issues_repo_name: neopenx/Dragon

max_forks_repo_name: neopenx/Dragon


text:
/*
 # This file is part of the Astrometry.net suite.
 # Licensed under a 3-clause BSD style license - see LICENSE
 */

#include <assert.h>
#include <gsl/gsl_matrix_double.h>
#include <gsl/gsl_vector_double.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_errno.h>
#include <stdarg.h>

#include "os-features.h"
#include "gslutils.h"
#include "errors.h"

static void errhandler(const char * reason,
                       const char * file,
                       int line,
                       int gsl_errno) {
    ERROR("GSL error: \"%s\" in %s:%i (gsl errno %i = %s)",
          reason, file, line,
          gsl_errno,
          gsl_strerror(gsl_errno));
}

void gslutils_use_error_system() {
    gsl_set_error_handler(errhandler);
}

int gslutils_invert_3x3(const double* A, double* B) {
    gsl_matrix* LU;
    gsl_permutation *p;
    gsl_matrix_view mB;
    int rtn = 0;
    int signum;

    p = gsl_permutation_alloc(3);
    gsl_matrix_const_view mA = gsl_matrix_const_view_array(A, 3, 3);
    mB = gsl_matrix_view_array(B, 3, 3);
    LU = gsl_matrix_alloc(3, 3);

    gsl_matrix_memcpy(LU, &mA.matrix);
    if (gsl_linalg_LU_decomp(LU, p, &signum) ||
        gsl_linalg_LU_invert(LU, p, &mB.matrix)) {
        ERROR("gsl_linalg_LU_decomp() or _invert() failed.");
        rtn = -1;
    }

    gsl_permutation_free(p);
    gsl_matrix_free(LU);
    return rtn;
}

void gslutils_matrix_multiply(gsl_matrix* C,
                              const gsl_matrix* A, const gsl_matrix* B) {
    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, A, B, 0.0, C);
}

int gslutils_solve_leastsquares_v(gsl_matrix* A, int NB, ...) {
    int i, res;
    gsl_vector**  B = malloc(NB * sizeof(gsl_vector*));
    // Whoa, three-star programming!
    gsl_vector*** X = malloc(NB * sizeof(gsl_vector**));
    gsl_vector*** R = malloc(NB * sizeof(gsl_vector**));

    gsl_vector** Xtmp = malloc(NB * sizeof(gsl_vector*));
    gsl_vector** Rtmp = malloc(NB * sizeof(gsl_vector*));

    va_list va;
    va_start(va, NB);
    for (i=0; i<NB; i++) {
        B[i] = va_arg(va, gsl_vector*);
        X[i] = va_arg(va, gsl_vector**);
        R[i] = va_arg(va, gsl_vector**);
    }
    va_end(va);

    res = gslutils_solve_leastsquares(A, B, Xtmp, Rtmp, NB);
    for (i=0; i<NB; i++) {
        if (X[i])
            *(X[i]) = Xtmp[i];
        else
            gsl_vector_free(Xtmp[i]);
        if (R[i])
            *(R[i]) = Rtmp[i];
        else
            gsl_vector_free(Rtmp[i]);
    }
    free(Xtmp);
    free(Rtmp);
    free(X);
    free(R);
    free(B);
    return res;
}

int gslutils_solve_leastsquares(gsl_matrix* A, gsl_vector** B,
                                gsl_vector** X, gsl_vector** resids,
                                int NB) {
    int i;
    gsl_vector *tau, *resid = NULL;
    Unused int ret;
    int M, N;

    M = A->size1;
    N = A->size2;

    for (i=0; i<NB; i++) {
        assert(B[i]);
        assert(B[i]->size == M);
    }

    tau = gsl_vector_alloc(MIN(M, N));
    assert(tau);

    ret = gsl_linalg_QR_decomp(A, tau);
    assert(ret == 0);
    // A,tau now contains a packed version of Q,R.

    for (i=0; i<NB; i++) {
        if (!resid) {
            resid = gsl_vector_alloc(M);
            assert(resid);
        }
        X[i] = gsl_vector_alloc(N);
        assert(X[i]);
        ret = gsl_linalg_QR_lssolve(A, tau, B[i], X[i], resid);
        assert(ret == 0);
        if (resids) {
            resids[i] = resid;
            resid = NULL;
        }
    }

    gsl_vector_free(tau);
    if (resid)
        gsl_vector_free(resid);

    return 0;
}

ext: c

max_stars_repo_licenses: Net-SNMP Xnet

max_issues_repo_licenses: Net-SNMP Xnet

max_forks_repo_licenses: Net-SNMP Xnet


text:
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>
#include <math.h>

#include "gauss.h"
#include "ts_fac.h"
#include "ts_tr_fac.h"

static double v_func(double diff){
	double denom = gsl_cdf_ugaussian_P(diff);
	double nom = gsl_ran_ugaussian_pdf(diff);
	return nom / denom;
}

static double w_func(double diff, double v){
	return (v*(v+diff));
}
static void ts_tr_fac_up(ts_fac_t *self, int useless){
	ts_var_t *var  = self->vars[TS_TR_SUM_VAR];
	gauss_t msg    = ts_fg_msg_to_var(var, self);
	gauss_t val    = gauss_div(var->value, msg);
	double sqrt_pi = sqrt(gauss_pi(val));
	double normal  = gauss_tau(val) / sqrt_pi;
	double v       = v_func(normal);
	double w       = w_func(normal, v);
	double denom   = 1.0 - w;
	double new_pi  = gauss_pi(val) / denom;
	double new_tau = (gauss_tau(val) + (sqrt_pi * v)) / denom;
	gauss_t nv     = gauss_init_prec(new_tau, new_pi);
	ts_var_update_val(var, self, nv);
}

static void ts_tr_fac_down(ts_fac_t *self, int useless){

}

ts_fac_t* ts_tr_fac_new(ts_fg_t *fg, ts_var_t *sum_var){
	ts_fac_t *self = ts_fac_new(fg);
	self->vars[TS_TR_SUM_VAR] = sum_var;
	self->up_func   = ts_tr_fac_up;
	self->down_func = ts_tr_fac_down;
	ts_fg_msg_to_var(sum_var, self) = gauss_init_prec(0.0, 0.0);

	return self;
}

max_stars_repo_name: cosmicturtle/gamesetmatch

max_issues_repo_name: cosmicturtle/gamesetmatch

max_forks_repo_name: cosmicturtle/gamesetmatch


text:
/**
 * @file pflow.h
 * @brief Public header file for power flow application
 */

#ifndef PFLOW_H
#define PFLOW_H

#include <petsc.h>
#include "cJSON.h"

typedef struct _p_PFLOW *PFLOW;

PETSC_EXTERN PetscErrorCode PFLOWCreate(MPI_Comm,PFLOW*);
PETSC_EXTERN PetscErrorCode PFLOWDestroy(PFLOW*);
PETSC_EXTERN PetscErrorCode PFLOWReadMatPowerData(PFLOW,const char[]);
PETSC_EXTERN PetscErrorCode PFLOWSolve(PFLOW);
PETSC_EXTERN PetscErrorCode PFLOWPostSolve(PFLOW);

/* FUNCTIONS FOR PARSING DATA */
int PFLOWParserGetVal(cJSON *, const char *, double *, int);
double PFLOWParserGet(PFLOW,int,char *, const char *);
void PFLOWParserSet(PFLOW,int,char *, PetscScalar, const char *);

#endif

max_stars_repo_name: GMLC-TDC/Use-Cases

max_issues_repo_name: GMLC-TDC/Use-Cases

max_forks_repo_name: GMLC-TDC/Use-Cases


text:
/*
 * Copyright 2020 Makani Technologies LLC
 *
 * 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.
 */

#ifndef SIM_PHYSICS_ROTOR_DATABASE_H_
#define SIM_PHYSICS_ROTOR_DATABASE_H_

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

#include <string>

#include "common/macros.h"

class RotorDatabase {
 public:
  explicit RotorDatabase(const std::string &filename);
  virtual ~RotorDatabase();

  // Calculates the rotor's thrust [N], applied torque [N-m], or
  // generated power [W] based on the angular rate, freestream
  // velocity, and air density.  Torque and power are defined to be
  // negative during thrusting.
  double CalcThrust(double angular_rate, double freestream_vel,
                    double air_density) const;
  double CalcTorque(double angular_rate, double freestream_vel,
                    double air_density) const;
  double CalcPower(double angular_rate, double freestream_vel,
                   double air_density) const;

  // Looks-up the value of the rotor's thrust, torque, or power
  // coefficients [#] at a given angular rate and freestream velocity
  // in a 2-D look-up table.  Values are interpolated linearly between
  // points and the inputs are saturated to the limits of the table.
  double LookupThrustCoeff(double angular_rate, double freestream_vel) const;
  double LookupTorqueCoeff(double angular_rate, double freestream_vel) const;
  double LookupPowerCoeff(double angular_rate, double freestream_vel) const;

 private:
  // Returns true if the rotor database follows our sign conventions:
  //
  // - Power coefficient decreases with increasing angular rate for
  //   the low freestream velocity case, i.e. power is positive in
  //   generation.
  //
  // - Thrust coefficient increases with increasing angular rate for
  //   the low freestream velocity case.
  bool IsValid() const;

  // Diameter of the propeller [m].
  double diameter_;

  // Vector of angular rates [rad/s], which define the rows of the
  // look-up table.
  gsl_vector *angular_rates_;

  // Vector of freestream velocities [m/s], which define the columns
  // of the look-up table.
  gsl_vector *freestream_vels_;

  // Matrix of thrust coefficients [#] defined over a tensor grid of
  // freestream velocities and angular rates.  The thrust coefficient
  // is defined as:
  //
  //   thrust_coeff = thrust / (air_density * angular_rate_hz^2 * diameter^4)
  //
  // See Eq. 9.18, Houghton & Carpenter, 4th ed.
  gsl_matrix *thrust_coeffs_;

  // Matrix of power coefficients [#] defined over a tensor grid of
  // freestream velocities and angular rates.  The power coefficient
  // is defined as:
  //
  //   power_coeff = power / (air_density * angular_rate_hz^3 * diameter^5)
  //
  // See Eq. 9.22, Houghton & Carpenter, 4th ed.
  gsl_matrix *power_coeffs_;

  DISALLOW_COPY_AND_ASSIGN(RotorDatabase);
};

#endif  // SIM_PHYSICS_ROTOR_DATABASE_H_


text:
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_statistics_double.h>
#include "tz_error.h"
#include "tz_constant.h"
#include "tz_image_lib.h"
#include "tz_objdetect.h"
#include "tz_imatrix.h"
#include "tz_stack_math.h"
#include "tz_stack_bwdist.h"
#include "tz_voxel_linked_list.h"
#include "tz_stack_sampling.h"
#include "tz_matlabio.h"
#include "tz_voxel_graphics.h"
#include "tz_darray.h"
#include "tz_geo3d_scalar_field.h"

INIT_EXCEPTION_MAIN(e)

int main(int argc, char* argv[])
{
  char neuron_name[100];
  char matlab_path[100];

  if (argc == 1) {
    strcpy(neuron_name, "fly_neuron");
  } else {
    strcpy(neuron_name, argv[1]);
  }

  if (argc < 3) {
    strcpy(matlab_path, "/Applications/MATLAB74/bin/matlab");
  } else {
    strcpy(matlab_path, argv[2]);
  }

  char file_path[100];

  /*
  sprintf(file_path, "../data/%s/mask.tif", neuron_name);
  Stack *stack = Read_Stack(file_path);
  Stack_Not(stack, stack);

  Set_Matlab_Path(matlab_path);
  Stack *dist = Stack_Bwdist(stack);
  */

  sprintf(file_path, "../data/%s/dist.tif", neuron_name);
  Stack *dist = Read_Stack(file_path);

  //Stack *seeds = Stack_Local_Max(dist, NULL, STACK_LOCMAX_ALTER1);
  sprintf(file_path, "../data/%s/seeds.tif", neuron_name);
  Stack *seeds = Read_Stack(file_path);

  //sprintf(file_path, "../data/%s/grow_bundle.tif", neuron_name);
  sprintf(file_path, "../data/%s/bundle.tif", neuron_name);
  Stack *bundle = Read_Stack(file_path);
  Stack_And(seeds, bundle, seeds);

  sprintf(file_path, "../data/%s/bundle_seed.tif", neuron_name);
  Write_Stack(file_path, seeds);
  printf("%s created.\n", file_path);

  Voxel_List *list = Stack_To_Voxel_List(seeds);
  Pixel_Array *pa = Voxel_List_Sampling(dist, list);
  sprintf(file_path, "../data/%s/bundle_seeds.pa", neuron_name);
  Pixel_Array_Write(file_path, pa);
  printf("%s created.\n", file_path);

  int i;
  //double *pa_array = (double *) pa->array;
  uint16 *pa_array16 = (uint16 *) pa->array;
  double *pa_array = darray_malloc(pa->size);
  for (i = 0; i < pa->size; i++) {
    pa_array[i] = sqrt(pa_array16[i]);
  }

  Voxel_P *voxel_array = Voxel_List_To_Array(list, 1, NULL, NULL);
  Geo3d_Scalar_Field *field = Make_Geo3d_Scalar_Field(pa->size);

  for (i = 0; i < pa->size; i++) {
    field->points[i][0] = voxel_array[i]->x;
    field->points[i][1] = voxel_array[i]->y;
    field->points[i][2] = voxel_array[i]->z;
    field->values[i] = pa_array[i];
  }

  sprintf(file_path, "../data/%s/seeds", neuron_name);
  Write_Geo3d_Scalar_Field(file_path, field);
  printf("%s created.\n", file_path);

  printf("mean: %g, std: %g\n", gsl_stats_mean(pa_array, 1, pa->size),
	 sqrt(gsl_stats_variance(pa_array, 1, pa->size)));
  printf("%g\n", gsl_stats_max(pa_array, 1, pa->size));

  return 0;
}


text:
#pragma once

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

#include "../math/Constants.h"
#include "../math/Exponentials.h"
#include "../math/General.h"
#include "../math/Point2.h"
#include "../math/Point3.h"
#include "../math/Random.h"
#include "../math/TrigFuncs.h"
#include "../math/Vec2.h"
#include "../math/Vec3.h"
#include "../utils/Concepts.h"
/// @brief custom mathematics implementaitons
namespace math { } // namespace math

/// @brief utility types
namespace utils { } // namespace utils

using utils::concepts::ConstructibleFrom;		// NOLINT
using utils::concepts::Copyable;				// NOLINT
using utils::concepts::CopyOrMovable;			// NOLINT
using utils::concepts::DefaultConstructible;	// NOLINT
using utils::concepts::Derived;					// NOLINT
using utils::concepts::FloatingPoint;			// NOLINT
using utils::concepts::Integral;				// NOLINT
using utils::concepts::Movable;					// NOLINT
using utils::concepts::NotCopyable;				// NOLINT
using utils::concepts::NotCopyOrMovable;		// NOLINT
using utils::concepts::NotDefaultConstructible; // NOLINT
using utils::concepts::NotMovable;				// NOLINT
using utils::concepts::NotPassable;				// NOLINT
using utils::concepts::NotPointer;				// NOLINT
using utils::concepts::NotReference;			// NOLINT
using utils::concepts::NotSemiRegular;			// NOLINT
using utils::concepts::Numeric;					// NOLINT
using utils::concepts::Passable;				// NOLINT
using utils::concepts::Pointer;					// NOLINT
using utils::concepts::Reference;				// NOLINT
using utils::concepts::SemiRegular;				// NOLINT
using utils::concepts::SignedIntegral;			// NOLINT
using utils::concepts::SignedNumeric;			// NOLINT

using math::Constants;	  // NOLINT
using math::Exponentials; // NOLINT
using math::General;	  // NOLINT
using math::Point2;		  // NOLINT
using math::Point2Idx;	  // NOLINT
using math::Point3;		  // NOLINT
using math::Point3Idx;	  // NOLINT
using math::random_value; // NOLINT
using math::Trig;		  // NOLINT
using math::Vec2;		  // NOLINT
using math::Vec2Idx;	  // NOLINT
using math::Vec3;		  // NOLINT
using math::Vec3Idx;	  // NOLINT

using gsl::narrow_cast; // NOLINT
template<typename T>
using NotNull = gsl::not_null<T*>; // NOLINT
template<typename T>
using Owner = gsl::owner<T*>; // NOLINT

template<FloatingPoint T = float>
inline constexpr auto degrees_to_radians(T degrees) noexcept -> T {
	return degrees * Constants<T>::pi / narrow_cast<T>(180.0);
}

template<FloatingPoint T = float>
inline constexpr auto radians_to_degrees(T radians) noexcept -> T {
	return radians * narrow_cast<T>(180.0) / Constants<T>::pi;
}

template<SignedNumeric T = float>
inline constexpr auto clamp(T val, T min, T max) noexcept -> T {
	if(val < min) {
		return min;
	}
	else if(val > max) {
		return max;
	}
	else {
		return val;
	}
}

template<typename... Args>
inline constexpr auto ignore(Args&&... args) noexcept -> void {
	std::ignore = std::tuple<Args...>(args...);
}

// clang-format off
#ifndef _MSC_VER
	// NOLINTNEXTLINE
	#define IGNORE_UNUSED_MACROS_START \
		_Pragma("GCC diagnostic push")\
		_Pragma("GCC diagnostic ignored \"-Wunused-macros\"")
#else
	// NOLINTNEXTLINE
	#define IGNORE_UNUSED_MACROS_START
#endif
// clang-format on

// clang-format off
#ifndef _MSC_VER
	// NOLINTNEXTLINE
	#define IGNORE_UNUSED_MACROS_STOP \
		_Pragma("GCC diagnostic pop")
#else
	// NOLINTNEXTLINE
	#define IGNORE_UNUSED_MACROS_STOP
#endif
// clang-format on

IGNORE_UNUSED_MACROS_START

// clang-format off
#ifndef _MSC_VER
	// NOLINTNEXTLINE
	#define IGNORE_PADDING_START \
		_Pragma("GCC diagnostic push") \
		_Pragma("GCC diagnostic ignored \"-Wpadded\"")
#else
	// NOLINTNEXTLINE
	#define IGNORE_PADDING_START \
		_Pragma("warning( push )") \
		_Pragma("warning( disable : 4820 )")
#endif
// clang-format on

// clang-format off
#ifndef _MSC_VER
	// NOLINTNEXTLINE
	#define IGNORE_PADDING_STOP \
		_Pragma("GCC diagnostic pop")
#else
	// NOLINTNEXTLINE
	#define IGNORE_PADDING_STOP \
		_Pragma("warning( pop )")
#endif
// clang-format on

// clang-format off
#ifndef _MSC_VER
	// NOLINTNEXTLINE
	#define IGNORE_UNUSED_TEMPLATES_START \
		_Pragma("GCC diagnostic push")\
		_Pragma("GCC diagnostic ignored \"-Wunused-template\"")
#else
	// NOLINTNEXTLINE
	#define IGNORE_UNUSED_TEMPLATES_START
#endif
// clang-format on

// clang-format off
#ifndef _MSC_VER
	// NOLINTNEXTLINE
	#define IGNORE_UNUSED_TEMPLATES_STOP \
		_Pragma("GCC diagnostic pop")
#else
	// NOLINTNEXTLINE
	#define IGNORE_UNUSED_TEMPLATES_STOP
#endif
// clang-format on

IGNORE_UNUSED_MACROS_STOP


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/time.h>
#include <math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

void read_matrix(int** index, int** matrix, double scaling, int N_kw, char* input_fileName)
{
    FILE *fp = fopen(input_fileName, "r");
    fscanf(fp, "%*[^\n]\n");

    for (int ii = 0; ii < N_kw; ii++)
        fscanf(fp, "%d,", &((*index)[ii]));
    fscanf(fp, "%*[^\n]\n");

    int tmp;
    for (int ii = 0; ii < N_kw; ii++)
    {
        for (int jj = 0; jj < N_kw; jj++)
        {
            fscanf(fp, "%d,", &tmp);
            (*matrix)[ii*N_kw + jj] = (int) (tmp * scaling);
        }
        fscanf(fp, "%*[^\n]\n");
    }
    fclose(fp);
}

void pad_matrix(int** matrix_padded, int** matrix, int N_kw, int N_doc, double sec_budget, int fixed_padding)
{
    // Initialising RNG
    const gsl_rng_type * T;
    gsl_rng * r;
    gsl_rng_env_setup();
    T = gsl_rng_default;
    r = gsl_rng_alloc(T);

    // perform padding on the keywords
    int ii, jj;
    #pragma omp parallel for private(ii)
    for (int ii = 0; ii < N_kw; ii++)
        (*matrix_padded)[ii*N_kw + ii] = 2 * (*matrix)[ii*N_kw + ii] + fixed_padding + (int) gsl_ran_laplace(r, 2.0/sec_budget);

    // perform padding
    #pragma omp parallel for private(ii, jj)
    for (ii = 0; ii < N_kw; ii++)
    {
        for (jj = 0; jj < N_kw; jj++)
        {
            if (ii > jj)
            {
                int n1 = (*matrix_padded)[ii*N_kw + ii] - (*matrix)[ii*N_kw + ii];
                int n2 = (*matrix_padded)[jj*N_kw + jj] - (*matrix)[jj*N_kw + jj];
                int x1 = gsl_ran_hypergeometric(r, n1, 2*N_doc-n1, (*matrix_padded)[jj*N_kw + jj]);
                int x2 = gsl_ran_hypergeometric(r, n2, 2*N_doc-n2, (*matrix_padded)[ii*N_kw + ii]);

                (*matrix_padded)[ii*N_kw + jj] = (*matrix)[ii*N_kw + jj] + x1 + x2;
                (*matrix_padded)[jj*N_kw + ii] = (*matrix_padded)[ii*N_kw + jj];
            }
        }
    }
    gsl_rng_free(r);
}

void observe_matrix(gsl_matrix* matrix_obs, int** matrix_padded, int N_kw)
{
    // perform observed count generation
    for (int ii = 0; ii < N_kw; ii++)
        for (int jj = 0; jj < N_kw; jj++)
            gsl_matrix_set(matrix_obs, ii, jj, (double) ((*matrix_padded)[ii*N_kw + jj]));
}

void permutation_generation(int* idx1, int* idx2, int** permutation_tmp, int** permutation, int** permutation_inv, gsl_matrix* matrix_obs, int** matrix, int fixed_padding, int N_kw, int N_obs, int N_doc, double sec_budget)
{
    double N1_lower, N1_upper, N2_lower, N2_upper;
    int check = -1;
    int count = 0;
    *idx1 = rand() % N_obs;
    *idx2 = -1;
    int idx_old = (*permutation)[*idx1];
    int idx_new = 0;

    do {
        idx_new = rand() % N_kw;
        count++;

        if ((*permutation_inv)[idx_new] >= 0)
        {
            *idx2 = (*permutation_inv)[idx_new];
            N1_lower = 2 * (*matrix)[idx_new*N_kw + idx_new] - 4 * sqrt((*matrix)[idx_new*N_kw + idx_new] * (N_doc - (*matrix)[idx_new*N_kw + idx_new]) / N_doc) - 4 / sec_budget + fixed_padding;
            N1_upper = 2 * (*matrix)[idx_new*N_kw + idx_new] + 4 * sqrt((*matrix)[idx_new*N_kw + idx_new] * (N_doc - (*matrix)[idx_new*N_kw + idx_new]) / N_doc) + 4 / sec_budget + fixed_padding;

            N2_lower = 2 * (*matrix)[idx_old*N_kw + idx_old] - 4 * sqrt((*matrix)[idx_old*N_kw + idx_old] * (N_doc - (*matrix)[idx_old*N_kw + idx_old]) / N_doc) - 4 / sec_budget + fixed_padding;
            N2_upper = 2 * (*matrix)[idx_old*N_kw + idx_old] + 4 * sqrt((*matrix)[idx_old*N_kw + idx_old] * (N_doc - (*matrix)[idx_old*N_kw + idx_old]) / N_doc) + 4 / sec_budget + fixed_padding;

            if ((gsl_matrix_get(matrix_obs, *idx1, *idx1) > N1_lower) && (gsl_matrix_get(matrix_obs, *idx1, *idx1) < N1_upper))
                if ((gsl_matrix_get(matrix_obs, *idx2, *idx2) > N2_lower) && (gsl_matrix_get(matrix_obs, *idx2, *idx2) > N2_lower))
                    check = 1;
        }
        else
        {
            *idx2 = -1;
            N1_lower = 2 * (*matrix)[idx_new*N_kw + idx_new] - 4 * sqrt((*matrix)[idx_new*N_kw + idx_new] * (N_doc - (*matrix)[idx_new*N_kw + idx_new]) / N_doc) - 4 / sec_budget + fixed_padding;
            N1_upper = 2 * (*matrix)[idx_new*N_kw + idx_new] + 4 * sqrt((*matrix)[idx_new*N_kw + idx_new] * (N_doc - (*matrix)[idx_new*N_kw + idx_new]) / N_doc) + 4 / sec_budget + fixed_padding;

            if ((gsl_matrix_get(matrix_obs, *idx1, *idx1) > N1_lower) && (gsl_matrix_get(matrix_obs, *idx1, *idx1) < N1_upper))
                check = 1;

        }

    } while ((check < 0) && (count < 200));

    if (count == 200)
        *idx2 = *idx1;

    if (*idx1 != *idx2)
    {
        (*permutation_tmp)[*idx1] = idx_new;
        if ((*permutation_inv)[idx_new] >= 0)
        {
            *idx2 = (*permutation_inv)[idx_new];
            (*permutation_tmp)[*idx2] = idx_old;
        }
    }
}

max_stars_repo_name: RethinkingSSE/Attacks-on-SSE

max_issues_repo_name: RethinkingSSE/Attacks-on-SSE

max_forks_repo_name: RethinkingSSE/Attacks-on-SSE


text:
/* multifit_nlinear/fdfvv.c
 *
 * Copyright (C) 2015 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 <gsl/gsl_math.h>
#include <gsl/gsl_multifit_nlinear.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>

/*
fdfvv()
  Compute approximate second directional derivative using
finite differences.

See Eq. 19 of:

M. K. Transtrum, J. P. Sethna, Improvements to the Levenberg
Marquardt algorithm for nonlinear least-squares minimization,
arXiv:1201.5885, 2012.

Inputs: h     - step size for finite difference
        x     - parameter vector, size p
        v     - geodesic velocity, size p
        f     - vector of function values f_i(x), size n
        J     - Jacobian matrix J(x), n-by-p
        swts  - data weights
        fdf   - fdf struct
        fvv   - (output) approximate second directional derivative
                vector D_v^2 f(x)
        work  - workspace, size p

Return: success or error
*/

static int
fdfvv(const double h, const gsl_vector *x, const gsl_vector *v,
      const gsl_vector *f, const gsl_matrix *J, const gsl_vector *swts,
      gsl_multifit_nlinear_fdf *fdf, gsl_vector *fvv, gsl_vector *work)
{
  int status;
  const size_t n = fdf->n;
  const size_t p = fdf->p;
  const double hinv = 1.0 / h;
  size_t i;

  /* compute work = x + h*v */
  for (i = 0; i < p; ++i)
    {
      double xi = gsl_vector_get(x, i);
      double vi = gsl_vector_get(v, i);

      gsl_vector_set(work, i, xi + h * vi);
    }

  /* compute f(x + h*v) */
  status = gsl_multifit_nlinear_eval_f (fdf, work, swts, fvv);
  if (status)
    return status;

  for (i = 0; i < n; ++i)
    {
      double fi = gsl_vector_get(f, i);    /* f_i(x) */
      double fip = gsl_vector_get(fvv, i); /* f_i(x + h*v) */
      gsl_vector_const_view row = gsl_matrix_const_row(J, i);
      double u, fvvi;

      /* compute u = sum_{ij} J_{ij} D v_j */
      gsl_blas_ddot(&row.vector, v, &u);

      fvvi = (2.0 * hinv) * ((fip - fi) * hinv - u);

      gsl_vector_set(fvv, i, fvvi);
    }

  return status;
}

/*
gsl_multifit_nlinear_fdfvv()
  Compute approximate second directional derivative
using finite differences

Inputs: h    - step size for finite difference
        x    - parameter vector, size p
        v    - geodesic velocity, size p
        f    - function values f_i(x), size n
        J    - Jacobian matrix J(x), n-by-p
        swts - sqrt data weights (set to NULL if not needed)
        fdf  - fdf
        fvv  - (output) approximate (weighted) second directional derivative
               vector, size n, sqrt(W) fvv
        work - workspace, size p

Return: success or error
*/

int
gsl_multifit_nlinear_fdfvv(const double h, const gsl_vector *x, const gsl_vector *v,
                           const gsl_vector *f, const gsl_matrix *J,
                           const gsl_vector *swts, gsl_multifit_nlinear_fdf *fdf,
                           gsl_vector *fvv, gsl_vector *work)
{
  return fdfvv(h, x, v, f, J, swts, fdf, fvv, work);
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


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

/* Author: MJB */

#ifndef __GSL_MONTE_MISER_H__
#define __GSL_MONTE_MISER_H__

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

#include <gsl/gsl_rng.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_plain.h>

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

__BEGIN_DECLS

typedef struct {
  size_t min_calls;
  size_t min_calls_per_bisection;
  double dither;
  double estimate_frac;
  double alpha;
  size_t dim;
  int estimate_style;
  int depth;
  int verbose;
  double * x;
  double * xmid;
  double * sigma_l;
  double * sigma_r;
  double * fmax_l;
  double * fmax_r;
  double * fmin_l;
  double * fmin_r;
  double * fsum_l;
  double * fsum_r;
  double * fsum2_l;
  double * fsum2_r;
  size_t * hits_l;
  size_t * hits_r;
} gsl_monte_miser_state;

GSL_FUN int gsl_monte_miser_integrate(gsl_monte_function * f,
                              const double xl[], const double xh[],
                              size_t dim, size_t calls,
                              gsl_rng *r,
                              gsl_monte_miser_state* state,
                              double *result, double *abserr);

GSL_FUN gsl_monte_miser_state* gsl_monte_miser_alloc(size_t dim);

GSL_FUN int gsl_monte_miser_init(gsl_monte_miser_state* state);

GSL_FUN void gsl_monte_miser_free(gsl_monte_miser_state* state);

typedef struct {
  double estimate_frac;
  size_t min_calls;
  size_t min_calls_per_bisection;
  double alpha;
  double dither;
} gsl_monte_miser_params;

GSL_FUN void gsl_monte_miser_params_get (const gsl_monte_miser_state * state,
				 gsl_monte_miser_params * params);

GSL_FUN void gsl_monte_miser_params_set (gsl_monte_miser_state * state,
				 const gsl_monte_miser_params * params);

__END_DECLS

#endif /* __GSL_MONTE_MISER_H__ */

max_stars_repo_name: mgreter/astrometrylib

max_issues_repo_name: zzpwahaha/VimbaCamJILA

max_forks_repo_name: zzpwahaha/VimbaCamJILA


text:
// Copyright (c) Microsoft Corporation.
// Licensed under the MIT license.

// clang-format off

#pragma once

#pragma warning(push)

// C
#include <climits>
#include <cwchar>
#include <cwctype>

// STL

// Block minwindef.h min/max macros to prevent <algorithm> conflict
#define NOMINMAX

#include <algorithm>
#include <atomic>
#include <deque>
#include <list>
#include <memory>
#include <map>
#include <mutex>
#include <shared_mutex>
#include <new>
#include <optional>
#include <queue>
#include <stdexcept>
#include <string>
#include <string_view>
#include <thread>
#include <tuple>
#include <utility>
#include <vector>
#include <unordered_map>
#include <iterator>
#include <math.h>
#include <sstream>
#include <fstream>
#include <iomanip>
#include <filesystem>
#include <functional>
#include <set>
#include <unordered_set>

// WIL
#include <wil/Common.h>
#include <wil/Result.h>
#include <wil/resource.h>
#include <wil/wistd_memory.h>
#include <wil/stl.h>
#include <wil/com.h>
#include <wil/filesystem.h>
#include <wil/win32_helpers.h>

// GSL
// Block GSL Multi Span include because it both has C++17 deprecated iterators
// and uses the C-namespaced "max" which conflicts with Windows definitions.
#ifndef BLOCK_GSL
#define GSL_MULTI_SPAN_H
#include <gsl/gsl>
#endif

// CppCoreCheck
#include <CppCoreCheck/Warnings.h>

// IntSafe
#define ENABLE_INTSAFE_SIGNED_FUNCTIONS
#include <intsafe.h>

// SAL
#include <sal.h>

// WRL
#include <wrl.h>

// TIL - Terminal Implementation Library
#include "til.h"

#pragma warning(pop)

// clang-format on

max_stars_repo_name: Vereis/terminal

max_issues_repo_name: Vereis/terminal

max_forks_repo_name: Vereis/terminal


text:
#pragma once

#include <cuda/runtime_api.hpp>
#include <gsl/gsl-lite.hpp>

namespace thrustshift {

namespace kernel {

template <typename A, class RangeX, class RangeY, class RangeD>
__global__ void axmy(A a, RangeX x, RangeY y, RangeD d) {
	const auto gtid = threadIdx.x + blockIdx.x * blockDim.x;
	gsl_Expects(x.size() == y.size());
	gsl_Expects(y.size() == d.size());
	if (gtid < x.size()) {
		d[gtid] = a * x[gtid] * y[gtid];
	}
}

} // namespace kernel

namespace async {

template <typename A, class RangeX, class RangeY, class RangeD>
void axmy(cuda::stream_t& stream, A a, RangeX&& x, RangeY&& y, RangeD&& d) {

	gsl_Expects(x.size() == y.size());
	gsl_Expects(y.size() == d.size());
	constexpr cuda::grid::block_dimension_t block_dim = 256;
	const cuda::grid::dimension_t grid_dim =
	    (x.size() + block_dim - 1) / block_dim;

	using RangeXD = typename std::remove_reference<RangeX>::type;
	using RangeYD = typename std::remove_reference<RangeY>::type;
	using RangeDD = typename std::remove_reference<RangeD>::type;

	cuda::enqueue_launch(kernel::axmy<A, RangeXD, RangeYD, RangeDD>,
	                     stream,
	                     cuda::make_launch_config(grid_dim, block_dim),
	                     a,
	                     x,
	                     y,
	                     d);
}

} // namespace async

} // namespace thrustshift

max_stars_repo_name: pauleonix/thrustshift

max_issues_repo_name: pauleonix/thrustshift

max_forks_repo_name: pauleonix/thrustshift


text:
#include "allelefreq.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <gsl/gsl_sf_psi.h>
#include <gsl/gsl_errno.h>

void P_update_simple(const uint8_t* G, const double* zetabeta, const double* zetagamma, const double* xi, const double* beta, const double* gamma, double* var_beta, double* var_gamma, long N, long L, long K)
{
    uint8_t genotype;
    long idx, n, l, k;
    double theta_beta_sum, theta_gamma_sum;
    double *var_beta_tmp, *var_gamma_tmp;

    var_beta_tmp = (double*) malloc(K * sizeof(double));
    var_gamma_tmp = (double*) malloc(K * sizeof(double));

    // loop over loci
    for (l=0; l<L; l++) {

        for (k=0; k<K; k++) {

            var_beta_tmp[k] = 0.0;
            var_gamma_tmp[k] = 0.0;
        }

        // loop over samples
        for (n=0; n<N; n++) {

            genotype = G[n*L+l];

            // missing data do not contribute
            if (genotype!=3) {

                // compute xi*zeta_{beta,gamma}
                theta_beta_sum = 0.0;
                theta_gamma_sum = 0.0;
                for (k=0; k<K; k++) {
                    theta_beta_sum += xi[n*K+k] * zetabeta[l*K+k];
                    theta_gamma_sum += xi[n*K+k] * zetagamma[l*K+k];
                }

                // increment var_{beta,gamma}_tmp
                for (k=0; k<K; k++) {
                    var_beta_tmp[k] += (double) genotype * xi[n*K+k] / theta_beta_sum;
                    var_gamma_tmp[k] += (double) (2-genotype) * xi[n*K+k] / theta_gamma_sum;
                }
            }
        }

        // compute var_{beta,gamma}
        for (k=0; k<K; k++) {
            idx = l*K+k;
            var_beta[idx] = beta[idx] + zetabeta[idx] * var_beta_tmp[k];
            var_gamma[idx] = gamma[idx] + zetagamma[idx] * var_gamma_tmp[k];
        }
    }

    free( var_beta_tmp );
    free( var_gamma_tmp );
}

void P_update_logistic(const double* Dvarbeta, const double* Dvargamma, const double* mu, const double* Lambda, double* var_beta, double* var_gamma, double mintol, long L, long K)
{
    /*
    `Dvarbeta` and `Dvargamma` are a function of the values of `var_beta` and `var_gamma`. This
    dependence, however, is explicit on a set of latent populations assignments which in-turn
    depend on the variational parameters estimated at the previous step. So, the variables
    `Dvarbeta` and `Dvargamma` do not have to be updated.
    */

    long l, k, idx, numvar, update;
    long iter = 0;
    double tol = 10.0;
    double tmptol;
    double beta, gamma;
    double A, pbetagamma, pbeta, pgamma, ppbeta, ppgamma;
    double A_1, B_1, C_1, A_2, B_2, C_2;

    /*
    Iterate until succesive estimates are sufficiently similar
    or the number of iterations exceeds 1000.
    */
    while (tol>mintol && iter<1000) {

        numvar = 0;
        tol = 0.0;

        // loop over loci
        for (l=0; l<L; l++) {

            // only update a locus, if all its variational parameters
            // satisfy positivity constraints.
            update = 1;
            for (k=0; k<K; k++) {
                if (var_beta[l*K+k]<=0 || var_gamma[l*K+k]<=0) {
                    update = 0;
                }
            }

            if (update==1) {

                // loop over populations
                for (k=0; k<K; k++) {

                    idx = l*K+k;
                    // compute pseudo-hyperparameters
                    pbetagamma = gsl_sf_psi_1(var_beta[idx]+var_gamma[idx]);
                    pbeta = gsl_sf_psi_1(var_beta[idx]);
                    pgamma = gsl_sf_psi_1(var_gamma[idx]);
                    ppbeta = gsl_sf_psi_n(2, var_beta[idx]);
                    ppgamma = gsl_sf_psi_n(2, var_gamma[idx]);

                    A_1 = pbeta-pbetagamma;
                    B_1 = -1.*pbetagamma;
                    A_2 = -1.*pbetagamma;
                    B_2 = pgamma-pbetagamma;
                    A = (gsl_sf_psi(var_beta[idx]) - gsl_sf_psi(var_gamma[idx]) - mu[l])*Lambda[k];
                    C_1 = -1.*A*pbeta - 0.5*Lambda[k]*ppbeta;
                    C_2 = A*pgamma - 0.5*Lambda[k]*ppgamma;

                    beta = (C_1*B_2-C_2*B_1)/(A_1*B_2-A_2*B_1);
                    gamma = (C_1*A_2-C_2*A_1)/(B_1*A_2-B_2*A_1);

                    // compute var_{beta,gamma}
                    tmptol = var_beta[idx];
                    var_beta[idx] = beta + Dvarbeta[idx];
                    tol += fabs(var_beta[idx]-tmptol);

                    tmptol = var_gamma[idx];
                    var_gamma[idx] = gamma + Dvargamma[idx];
                    tol += fabs(var_gamma[idx]-tmptol);

                    numvar += 1;

                }
            }
        }

        // compute convergence tolerance
        tol = 0.5*tol/numvar;
        iter += 1;
    }

    // printf("tol = %.8f, iter = %lu\n",tol,iter);
}

max_stars_repo_name: gadgrandez/fastStructure

max_issues_repo_name: gadgrandez/fastStructure

max_forks_repo_name: gadgrandez/fastStructure


text:
/*

  simulate NMR peakshape for two-site chemical exchange

  written by Sung-Hun Bae
  Nov. 17, 2008

*/

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <complex>
#include <gsl/gsl_fft_complex.h>

#define TINY  1.0e-12

using namespace std;

enum SITE_LABEL { A, B, N_SITE };
enum QUADRATURE { REAL, IMAG };
typedef struct {
  double Omegaa;// frequency of A site (rad/s)
  double Omegab;// frequency of B site (rad/s)
  double kab;	// reaction rate constant, A -> B (1/s)
  double kba;	// reaction rate constant, B -> A (1/s)
  double R20a;	// transverse relaxation rate of A site (1/s)
  double R20b;	// transverse relaxation rate of B site	(1/s)
  double M0a;	// initial magnetization
  double M0b;	// initial magnetization
  SITE_LABEL site;	// A or B
  QUADRATURE quad;	// REAL or IMAG
  } two_site_exchange;

//
// Magnetization (two-site exchange)
//

void N (double t, double &R, double &I, void * params) {
  two_site_exchange * x = (two_site_exchange *) params;
  complex<double> c1,c2,c3,lp,lm,a11,a12,a21,a22;
  complex<double> Omegaa = sqrt(complex<double>(-1.))*x->Omegaa;
  complex<double> Omegab = sqrt(complex<double>(-1.))*x->Omegab;
  double mag;
  c1 = (-Omegaa - Omegab + x->R20a + x->R20b + x->kab + x->kba);
  c2 = (-Omegaa + Omegab + x->R20a - x->R20b + x->kab - x->kba);
  c3 = sqrt(c2*c2+4.*x->kab*x->kba);
  lp = 0.5*(c1+c3);
  lm = 0.5*(c1-c3);
  if (x->site == A) {
    a11 = 0.5*((1.-c2/c3)*exp(-t*lm)+(1.+c2/c3)*exp(-t*lp));
    a12 = x->kba/c3*(exp(-t*lm)-exp(-t*lp));
    R += real(a11*x->M0a+a12*x->M0b);
    I += imag(a11*x->M0a+a12*x->M0b);
    }
  if (x->site == B) {
    a22 = 0.5*((1.+c2/c3)*exp(-t*lm)+(1.-c2/c3)*exp(-t*lp));
    a21 = x->kab/c3*(exp(-t*lm)-exp(-t*lp));
    R += real(a21*x->M0a+a22*x->M0b);
    I += imag(a21*x->M0a+a22*x->M0b);
    }
}

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

  if (argc < 7 || argc > 17) {
    printf("\t\t\t\t\t\t\tSung-Hun Bae 2008.11\n");
    printf("   ");
    printf("Usage: xpshape [-model AB|ACB]\n");
    printf("   ");
    printf("               [-dw #(Hz)] [-kab #(1/s)] [-kba #(1/s)]\n");
    printf("   ");
    printf("               [-R20a #(1/s)] [-R20b #(1/s)]\n");
    printf("   ");
    printf("               [-a0 conc(M)] [-c0 conc(M)]\n");
    printf("   ");
    printf("               [-swh Hz(600)] [-fid pts(8192)]\n\n");
    printf("   ");
    printf("simulate NMR peakshape for two-site chemical exchange\n\n");
    printf("   ");
    printf("-model,AB  : A = B     (dw,kab,kba,R20a,R20b required)\n");
    printf("   ");
    printf("       ACB : A + C = B (dw,kab,kba,R20a,R20b,a0,c0 required)\n");
    printf("\n");
    printf("   ");
    printf("-dw,  (A.frequency= -dw/2, B.frequency= +dw/2)\n");
    printf("   ");
    printf("-R20a,(exchange-free R2 of A state)\n");
    printf("   ");
    printf("-R20b,(exchange-free R2 of B state)\n");
    printf("   ");
    printf("-kab, (rate constant of A->B)\n");
    printf("   ");
    printf("-kba, (rate constant of B->A)\n");
    printf("   ");
    printf("-a0,  (initial concentration of A)\n");
    printf("   ");
    printf("-c0,  (initial concentration of C)\n");
    printf("   ");
    printf("-swh, (spectrum = -swh/2 ... +swh/2)\n");
    printf("   ");
    printf("-fid, (must be ...,1024,2048,4096,8192,16384,32768,65536,..)\n");
    printf("\n");

    exit (1);
    }

  // default settings

  double SWH = 600;		  // (Hz)
  unsigned int FID_SIZE = 8192;	  // (complex data points)

  // receiving parameters

  unsigned int i = 1;

  string model;
  double dw,R20a,R20b,kab,kba,concA0,concC0;
  bool dw_,R20a_,R20b_,kab_,kba_,concA0_,concC0_;
  double pb,Kd,concA,concB,concC,C,u,v;

  dw_ = R20a_ = R20b_ = kab_ = kba_ = concA0_ = concC0_ = false;

  while ((argc - i) > 0) {
    if ((argc-i) > 0 && !strcmp(argv[i],"-model")) {
      model = string(argv[++i]);
      i++;
      }
    if ((argc-i) > 0 && !strcmp(argv[i],"-dw")) {
      dw = 2.*M_PI*atof(argv[++i]); // (Hz) -> (rad/s)
      dw_ = true;
      i++;
      }
    if ((argc-i) > 0 && !strcmp(argv[i],"-kab")) {
      kab = atof(argv[++i]);
      kab_ = true;
      i++;
      }
    if ((argc-i) > 0 && !strcmp(argv[i],"-kba")) {
      kba = atof(argv[++i]);
      kba_ = true;
      i++;
      }
    if ((argc-i) > 0 && !strcmp(argv[i],"-R20a")) {
      R20a = atof(argv[++i]);
      R20a_ = true;
      i++;
      }
    if ((argc-i) > 0 && !strcmp(argv[i],"-R20b")) {
      R20b = atof(argv[++i]);
      R20b_ = true;
      i++;
      }
    if ((argc-i) > 0 && !strcmp(argv[i],"-a0")) {
      concA0 = atof(argv[++i]);
      concA0_ = true;
      i++;
      }
    if ((argc-i) > 0 && !strcmp(argv[i],"-c0")) {
      concC0 = atof(argv[++i]);
      concC0_ = true;
      i++;
      }
    if ((argc-i) > 0 && !strcmp(argv[i],"-swh")) {
      SWH = atof(argv[++i]);
      i++;
      }
    if ((argc-i) > 0 && !strcmp(argv[i],"-fid")) {
      FID_SIZE = atoi(argv[++i]);
      i++;
      }
    }

  double fid[2*FID_SIZE];	    // complex points
  double dwell = 1./(SWH);	    // (sec)
  double norm = 1./sqrt(FID_SIZE);  // scaling factor

  if (model == "AB" && dw_ && R20a_ && R20b_ && kab_ && kba_) {
    // dw, R20a, R20b, kab, kba
    Kd = kba/kab;
    pb = 1./(1.+Kd);
    C = 1.0;
    }
  else if (model == "ACB" && dw_ && R20a_ && R20b_ && kab_ && kba_
    && concA0_ && concC0_) {
    // dw, R20a, R20b, kab(=kon), kba(=koff), concA0, concC0
    kab = (concC0 > TINY ? kab : 0.0);
    kba = (concC0 > TINY ? kba : 0.0);
    Kd = ((kab > TINY) && (kba > TINY) ? kba/kab : 0.0);
    u = 0.5*(concA0+concC0+Kd);
    v = (concC0 > TINY ? (u-sqrt(u*u-concA0*concC0)) : 0.0);
    concA = (v > TINY ? concA0-v : concA0);
    concB = (v > TINY ? v : 0.0);
    concC = (v > TINY ? concC0 - v : concC0);
    pb = (concC0 > TINY ? kab*concC/(kab*concC+kba) : 0.0);
    C = concC;
    }
  else {
    printf("error: missing parameter(s)\n");
    exit(1);
    }

  two_site_exchange x = {-dw/2.,dw/2.,kab*C,kba,R20a,R20b,(1.-pb),pb,A,REAL};

  // header information
  printf("# peak shape simulation ");
  if (model == "AB")  printf("[A = B]\n");
  if (model == "ACB") printf("[A + C = B]\n");

  printf("# dw     = %g (rad/s)\n", dw);
  printf("# kab    = %g (1/s)\n",kab);
  printf("# kba    = %g (1/s)\n",kba);
  printf("# R20a   = %g (1/s)\n",R20a);
  printf("# R20b   = %g (1/s)\n",R20b);
  if (model == "AB") {
  printf("# pb     = %g\n", pb);
  }
  if (model == "ACB") {
  printf("# kon    = %g (1/s)\n",kab);
  printf("# koff   = %g (1/s)\n",kba);
  printf("# Kd     = %g (M)\n",Kd);
  printf("# pb     = %g\n",pb);
  printf("# [A]ini = %g (M)\n",concA0);
  printf("# [C]ini = %g (M)\n",concC0);
  printf("# [A]    = %g (M)\n",concA);
  printf("# [C]    = %g (M)\n",concC);
  printf("# [B]    = %g (M)\n",concB);
  }
  printf("# FID    = %d (complex points)\n",FID_SIZE);
  printf("# SWH    = %g (Hz)\n",SWH);
  printf("# X-axis unit: (Hz)\n");

  // creating complex data points
  for (i = 0; i < FID_SIZE; i++) {
    fid[2*i] = fid[2*i+1] = 0.0;
    x.site = A;
    N (i*dwell, fid[2*i], fid[2*i+1], &x);  // A
    x.site = B;
    N (i*dwell, fid[2*i], fid[2*i+1], &x);  // B
    }

  // Fast Fourier Transform (FFT)
  gsl_fft_complex_radix2_forward (fid, 1, FID_SIZE);

  // spectrum
  for (i = FID_SIZE/2; i < FID_SIZE; i++)
    printf("%g %g\n",(-0.5*SWH + SWH*(i-FID_SIZE/2)/FID_SIZE),norm*fid[2*i]);
  for (i = 0; i <= FID_SIZE/2; i++)
    printf("%g %g\n",(-0.5*SWH + SWH*(i+FID_SIZE/2)/FID_SIZE),norm*fid[2*i]);
  printf("&\n");

  return 0;
}

/*

  For physical applications it is important to remember that the index
  appearing in the DFT does not correspond directly to a physical frequency.
  If the time-step of the DFT is \Delta then the frequency-domain includes
  both positive and negative frequencies, ranging from -1/(2\Delta) through
  0 to +1/(2\Delta).
  The positive frequencies are stored from the beginning of the array up to
  the middle, and the negative frequencies are stored backwards from the
  end of the array.

     index    z               x = FFT(z)

     0        z(t = 0)        x(f = 0)
     1        z(t = 1)        x(f = 1/(N Delta))
     2        z(t = 2)        x(f = 2/(N Delta))
     .        ........        ..................
     N/2      z(t = N/2)      x(f = +1/(2 Delta),
                                    -1/(2 Delta))
     .        ........        ..................
     N-3      z(t = N-3)      x(f = -3/(N Delta))
     N-2      z(t = N-2)      x(f = -2/(N Delta))
     N-1      z(t = N-1)      x(f = -1/(N Delta))

*/

max_stars_repo_name: sunghunbae/xpshape

max_issues_repo_name: sunghunbae/xpshape

max_forks_repo_name: sunghunbae/xpshape


text:
#pragma once

#include "poll_event.h"
#include "poll_response.h"
#include "poll_target.h"

#include <gsl/span>

#include <optional>
#include <vector>

namespace wrappers::zmq {

[[nodiscard]] std::optional<std::vector<poll_response>>
blocking_poll(gsl::span<poll_target> targets) noexcept;

} // namespace wrappers::zmq

max_stars_repo_name: Longhanks/linkollector-win

max_issues_repo_name: Longhanks/linkollector-win

max_forks_repo_name: Longhanks/linkollector-win


text:
/* specfunc/zeta.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004 Gerard Jungman
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 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.
 */

/* Author:  G. Jungman */

/* This file was taken from the GNU Scientific Library. Some modifications
 * were done in order to make it independent from the rest of GSL
 */

/*
#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_elementary.h>
#include <gsl/gsl_sf_exp.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_pow_int.h>
#include <gsl/gsl_sf_zeta.h>

#include "error.h"

#include "chebyshev.h"
#include "cheb_eval.c"
*/

#include <math.h>
#include <stdio.h>
#include "error.h"

/*-*-*-*-*-*-*-*-*-*- From gsl_machine.h -*-*-*-*-*-*-*-*-*-*-*-*-*/

#define GSL_LOG_DBL_MIN   (-7.0839641853226408e+02)
#define GSL_LOG_DBL_MAX    7.0978271289338397e+02
#define GSL_DBL_EPSILON        2.2204460492503131e-16

/*-*-*-*-*-*-*-*-*-* From gsl_sf_result.h *-*-*-*-*-*-*-*-*-*-*-*/

struct gsl_sf_result_struct {
  double val;
  double err;
};
typedef struct gsl_sf_result_struct gsl_sf_result;

/*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/

/* coefficients for Maclaurin summation in hzeta()
 * B_{2j}/(2j)!
 */
static double hzeta_c[15] = {
  1.00000000000000000000000000000,
  0.083333333333333333333333333333,
 -0.00138888888888888888888888888889,
  0.000033068783068783068783068783069,
 -8.2671957671957671957671957672e-07,
  2.0876756987868098979210090321e-08,
 -5.2841901386874931848476822022e-10,
  1.3382536530684678832826980975e-11,
 -3.3896802963225828668301953912e-13,
  8.5860620562778445641359054504e-15,
 -2.1748686985580618730415164239e-16,
  5.5090028283602295152026526089e-18,
 -1.3954464685812523340707686264e-19,
  3.5347070396294674716932299778e-21,
 -8.9535174270375468504026113181e-23
};

/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/

static int gsl_sf_hzeta_e(const double s, const double q, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(s <= 1.0 || q <= 0.0) {
	PLFIT_ERROR("s must be larger than 1.0 and q must be larger than zero", PLFIT_EINVAL);
  }
  else {
    const double max_bits = 54.0;
    const double ln_term0 = -s * log(q);

    if(ln_term0 < GSL_LOG_DBL_MIN + 1.0) {
	  PLFIT_ERROR("underflow", PLFIT_UNDRFLOW);
    }
    else if(ln_term0 > GSL_LOG_DBL_MAX - 1.0) {
	  PLFIT_ERROR("overflow", PLFIT_OVERFLOW);
    }
    else if((s > max_bits && q < 1.0) || (s > 0.5*max_bits && q < 0.25)) {
      result->val = pow(q, -s);
      result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
      return PLFIT_SUCCESS;
    }
    else if(s > 0.5*max_bits && q < 1.0) {
      const double p1 = pow(q, -s);
      const double p2 = pow(q/(1.0+q), s);
      const double p3 = pow(q/(2.0+q), s);
      result->val = p1 * (1.0 + p2 + p3);
      result->err = GSL_DBL_EPSILON * (0.5*s + 2.0) * fabs(result->val);
      return PLFIT_SUCCESS;
    }
    else {
      /* Euler-Maclaurin summation formula
       * [Moshier, p. 400, with several typo corrections]
       */
      const int jmax = 12;
      const int kmax = 10;
      int j, k;
      const double pmax  = pow(kmax + q, -s);
      double scp = s;
      double pcp = pmax / (kmax + q);
      double ans = pmax*((kmax+q)/(s-1.0) + 0.5);

      for(k=0; k<kmax; k++) {
        ans += pow(k + q, -s);
      }

      for(j=0; j<=jmax; j++) {
        double delta = hzeta_c[j+1] * scp * pcp;
        ans += delta;
        if(fabs(delta/ans) < 0.5*GSL_DBL_EPSILON) break;
        scp *= (s+2*j+1)*(s+2*j+2);
        pcp /= (kmax + q)*(kmax + q);
      }

      result->val = ans;
      result->err = 2.0 * (jmax + 1.0) * GSL_DBL_EPSILON * fabs(ans);
      return PLFIT_SUCCESS;
    }
  }
}

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

double gsl_sf_hzeta(const double s, const double a)
{
  gsl_sf_result result;
  gsl_sf_hzeta_e(s, a, &result);
  return result.val;
}

max_stars_repo_name: jmazon/haskell-igraph

max_issues_repo_name: jmazon/haskell-igraph

max_forks_repo_name: jmazon/haskell-igraph