text:
#include <gmock/gmock.h>
#include <gtest/gtest.h>

#include "boost_coroutine/http_client.h"
#include "callback/http_client.h"
#include "coroutine_ts/http_client.h"
#include "sync/http_client.h"
#include "task/http_client.h"
#include "threaded/http_client.h"

#include "executor.h"

#include <boost/algorithm/string.hpp>

#include <gsl/gsl_util>

#if 1
static inline auto log_msg = std::stringstream{};
#else
static inline auto& log_msg = std::cerr;
#endif

static HttpResponse request_get_impl(std::string_view uri) {
  if (uri == "/file1") {
    return Redirect{"/newplaceoffile1"};
  }
  if (uri == "/newplaceoffile1") {
    return Redirect{"/file2"};
  }
  if (uri == "/file2") {
    return Content{"content1"};
  }
  if (uri == "/file3") {
    return Content{"content2"};
  }
  if (boost::algorithm::starts_with(uri, "/extremeredirect")) {
    std::vector<std::string> parts;
    boost::algorithm::split(parts, uri, [](char c) { return c == '/'; });
    if (parts.size() == 2) {
      return Redirect{std::string{uri} + "/1"};
    }
    if (parts.size() == 3) {
      auto num = std::stol(parts[2]);
      if (num < 4) {
        return Redirect{"/extremeredirect/" + std::to_string(num + 1)};
      }
    }
    return Content{"finally here"};
  }
  return Content{""};
}

namespace synchronous {
class ConcreteHttpClient : public HttpClient {
public:
  HttpResponse request_get(std::string_view uri) override {
    log_msg << "request " << uri << "..." << std::endl;
    auto out
      = gsl::finally([&] { log_msg << "...request " << uri << std::endl; });
    return request_get_impl(uri);
  }
};
}

namespace threaded {
class ConcreteHttpClient : public HttpClient {
public:
  HttpResponse request_get(std::string_view uri) override {
    log_msg << "request " << uri << "..." << std::endl;
    auto out
      = gsl::finally([&] { log_msg << "...request " << uri << std::endl; });
    return request_get_impl(uri);
  }
};
}

namespace callback {

class ConcreteHttpClient : public HttpClient {
public:
  void request_get(const std::string&                       uri,
                   const std::function<void(HttpResponse)>& cb) override {
    log_msg << "request " << uri << "..." << std::endl;
    auto handle = [uri]() -> HttpResponse { return request_get_impl(uri); };
    get_current_executor().add_work([=] {
      log_msg << "...request " << uri << std::endl;
      get_current_executor().add_work([=] { cb(handle()); });
    });
  }
};
}

namespace task_based {

class ConcreteHttpClient : public HttpClient {
public:
  task<HttpResponse> request_get(const std::string& uri) override {
    log_msg << "request " << uri << "..." << std::endl;
    const auto handle
      = [uri]() -> HttpResponse { return request_get_impl(uri); };

    auto fut = make_task<HttpResponse>();
    get_current_executor().add_work([=] {
      log_msg << "...request " << uri << std::endl;

      fut->set_result(handle());
    });
    return fut;
  }
};
}

namespace boost_coroutine {

class ConcreteHttpClient : public HttpClient {
public:
  task<HttpResponse> request_get(const std::string& uri) override {
    log_msg << "request " << uri << "..." << std::endl;
    const auto handle = [=]() -> HttpResponse { return request_get_impl(uri); };

    auto fut = make_task<HttpResponse>();
    get_current_executor().add_work([=] {
      log_msg << "...request " << uri << std::endl;
      fut->set_result(handle());
    });
    return fut;
  }
};
}

namespace coroutines_ts {

class ConcreteHttpClient : public HttpClient {
public:
  task<HttpResponse> request_get(const std::string& uri) override {
    log_msg << "request " << uri << "..." << std::endl;
    const auto handle = [=]() -> HttpResponse { return request_get_impl(uri); };

    auto fut = make_task<HttpResponse>();
    get_current_executor().add_work([=] {
      log_msg << "...request " << uri << std::endl;
      fut->set_result(handle());
    });
    return fut;
  }
};
}


text:
/**
 *
 * @file core_ctrtri.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Julien Langou
 * @author Henricus Bouwmeester
 * @author Mathieu Faverge
 * @date 2010-11-15
 * @generated c Tue Jan  7 11:44:46 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 * @ingroup CORE_PLASMA_Complex32_t
 *
 *  CORE_ctrtri - Computes the inverse of a complex upper or lower
 *  triangular matrix A.
 *
 *******************************************************************************
 *
 * @param[in] uplo
 *          = PlasmaUpper: Upper triangle of A is stored;
 *          = PlasmaLower: Lower triangle of A is stored.
 *
 * @param[in] diag
 *          = PlasmaNonUnit: A is non-unit triangular;
 *          = PlasmaUnit:    A is unit triangular.
 *
 * @param[in] N
 *          The order of the matrix A. N >= 0.
 *
 * @param[in,out] A
 *          On entry, the triangular matrix A.  If UPLO = 'U', the
 *          leading N-by-N upper triangular part of the array A
 *          contains the upper triangular matrix, and the strictly
 *          lower triangular part of A is not referenced.  If UPLO =
 *          'L', the leading N-by-N lower triangular part of the array
 *          A contains the lower triangular matrix, and the strictly
 *          upper triangular part of A is not referenced.  If DIAG =
 *          'U', the diagonal elements of A are also not referenced and
 *          are assumed to be 1.  On exit, the (triangular) inverse of
 *          the original matrix.
 *
 * @param[in] LDA
 *          The leading dimension of the array A. LDA >= max(1,N).
 *
 * @param[out] info
 *          - 0 on successful exit
 *          - <0 if -i, the i-th argument had an illegal value
 *          - >0 if i, A(i,i) is exactly zero.  The triangular
 *          matrix is singular and its inverse can not be computed.
 *
 ******************************************************************************/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_ctrtri = PCORE_ctrtri
#define CORE_ctrtri PCORE_ctrtri
#endif
void CORE_ctrtri(PLASMA_enum uplo, PLASMA_enum diag, int N, PLASMA_Complex32_t *A, int LDA, int *info)
{
    *info = LAPACKE_ctrtri_work(
        LAPACK_COL_MAJOR,
        lapack_const(uplo), lapack_const(diag),
        N, A, LDA);
}

ext: c


text:
/*
** hemodynamic modelling for single-subject LISA
**
** G.Lohmann, MPI-KYB,  2018
*/
#include <viaio/Vlib.h>
#include <viaio/VImage.h>
#include <viaio/mu.h>
#include <viaio/option.h>

#include <gsl/gsl_cblas.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_spline.h>

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

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

#define LEN     10000   /* buffer length        */
#define NTRIALS 10000   /* max number of trials */

/* standard parameter values for gamma function, Glover 99 */
double a1 = 6;
double b1 = 0.9;
double a2 = 12;
double b2 = 0.9;
double cc = 0.35;

typedef struct TrialStruct {
  int   id;
  float onset;
  float duration;
  float height;
} Trial;

void printmat(gsl_matrix *R,char *str)
{
  int i,j;
  fprintf(stderr," %s: \n",str);
  for (i=0; i<R->size1; i++) {
    for (j=0; j<R->size2; j++) {
      fprintf(stderr," %9.6f",gsl_matrix_get(R,i,j));
    }
    fprintf(stderr,"\n");
  }
  fprintf(stderr,"\n");
}

void tprintmat(gsl_matrix *R,char *str)
{
  int i,j;
  fprintf(stderr," %s: \n",str);
  for (i=0; i<R->size2; i++) {
    for (j=0; j<R->size1; j++) {
      fprintf(stderr," %9.6f",gsl_matrix_get(R,j,i));
    }
    fprintf(stderr,"\n");
  }
  fprintf(stderr,"\n");
}

void printvec(gsl_vector *x,char *str)
{
  int i;
  fprintf(stderr," %s: \n",str);
  for (i=0; i<x->size; i++) {
    fprintf(stderr," %f\n",x->data[i]);
  }
  fprintf(stderr,"\n");
}

/* output txt file for plotting */
void PlotDesign(gsl_matrix *X,double tr,VString filename)
{
  int i,j;
  double u=0;

  FILE *fp = fopen(filename,"w");
  if (fp == NULL) VError("err opening plot file");
  double t = 0;
  for (i=0; i<X->size1; i++) {
    fprintf(fp," %8.2f",t);
    for (j=0; j<X->size2; j++) {
      u = gsl_matrix_get(X,i,j);
      fprintf(fp," %10.6f",u);
    }
    t += tr;
    fprintf(fp,"\n");
  }
  fclose(fp);
}

/* create design X */
gsl_matrix *VCreateDesign(int ntimesteps,int nevents,int hemomodel,VBoolean firstcol,gsl_matrix *covariates)
{
  int dim = nevents;
  if (hemomodel == 1) dim = nevents*2;
  if (hemomodel == 2) dim = nevents*3;
  if (firstcol == TRUE) dim++;
  if (covariates != NULL) dim += covariates->size2;  /* nuisance covariates without task labels */

  gsl_matrix *X = gsl_matrix_calloc(ntimesteps,dim);
  return X;
}

/* Glover kernel, gamma function */
double xgamma(double xx, double t0)
{
  double x, y, scale = 20;
  double y1, y2;
  double d1, d2;
  x = xx - t0;
  if(x < 0 || x > 50)
    return 0;
  d1 = a1 * b1;
  d2 = a2 * b2;
  y1 = pow(x / d1, a1) * exp(-(x - d1) / b1);
  y2 = pow(x / d2, a2) * exp(-(x - d2) / b2);
  y = y1 - cc * y2;
  y /= scale;
  return y;
}

/* Glover kernel, gamma function, parameters changed for block designs */
double bgamma(double xx, double t0)
{
  double x, y, scale = 120;
  double y1, y2;
  double d1, d2;
  double aa1 = 6;
  double bb1 = 0.9;
  double aa2 = 12;
  double bb2 = 0.9;
  double cx  = 0.1;
  x = xx - t0;
  if(x < 0 || x > 50)
    return 0;
  d1 = aa1 * bb1;
  d2 = aa2 * bb2;
  y1 = pow(x / d1, aa1) * exp(-(x - d1) / bb1);
  y2 = pow(x / d2, aa2) * exp(-(x - d2) / bb2);
  y = y1 - cx * y2;
  y /= scale;
  return y;
}

/* first derivative */
double deriv1_gamma(double x, double t0)
{
  double d1, d2, y1, y2, y, xx;
  double scale = 20.0;
  xx = x - t0;
  if(xx < 0 || xx > 50) return 0;
  d1 = a1 * b1;
  d2 = a2 * b2;
  y1 = pow(d1, -a1) * a1 * pow(xx, (a1 - 1.0)) * exp(-(xx - d1) / b1)
    - (pow((xx / d1), a1) * exp(-(xx - d1) / b1)) / b1;
  y2 = pow(d2, -a2) * a2 * pow(xx, (a2 - 1.0)) * exp(-(xx - d2) / b2)
    - (pow((xx / d2), a2) * exp(-(xx - d2) / b2)) / b2;
  y = y1 - cc * y2;
  y /= scale;
  return y;
}

/* second derivative */
double deriv2_gamma(double x, double t0)
{
  double d1, d2, y1, y2, y3, y4, y, xx;
  double scale = 20.0;
  xx = x - t0;
  if (xx < 0 || xx > 50) return 0;
  d1 = a1 * b1;
  d2 = a2 * b2;
  y1 = pow(d1, -a1) * a1 * (a1 - 1) * pow(xx, a1 - 2) * exp(-(xx - d1) / b1)
    - pow(d1, -a1) * a1 * pow(xx, (a1 - 1)) * exp(-(xx - d1) / b1) / b1;
  y2 = pow(d1, -a1) * a1 * pow(xx, a1 - 1) * exp(-(xx - d1) / b1) / b1
    - pow((xx / d1), a1) * exp(-(xx - d1) / b1) / (b1 * b1);
  y1 = y1 - y2;
  y3 = pow(d2, -a2) * a2 * (a2 - 1) * pow(xx, a2 - 2) * exp(-(xx - d2) / b2)
    - pow(d2, -a2) * a2 * pow(xx, (a2 - 1)) * exp(-(xx - d2) / b2) / b2;
  y4 = pow(d2, -a2) * a2 * pow(xx, a2 - 1) * exp(-(xx - d2) / b2) / b2
    - pow((xx / d2), a2) * exp(-(xx - d2) / b2) / (b2 * b2);
  y2 = y3 - y4;
  y = y1 - cc * y2;
  y /= scale;
  return y;
}

/* Gaussian function */
double xgauss(double xx, double t0)
{
  double sigma=1.0;
  double x, y, z, a = 2.506628273;
  x = (xx - t0);
  z = x / sigma;
  y = exp(-0.5*z*z) / (sigma * a);
  return y;
}

void XConvolve(double *src,double *dst,double *kernel,int nt,int kernelsize)
{
  int i,j,jj,k;
  double sum=0;

  for (i=0; i<nt; i++) {
    sum = 0;
    k=0;
    for (j=i; j<i+kernelsize; j++) {
      jj = i-k;
      if (jj >= 0 && jj < nt) {
	sum += src[jj] * kernel[k];
      }
      k++;
    }
    dst[i] = sum;
  }
}

/* hemodynamic modelling at high temporal resolution */
void VHemoModel(Trial *trial,int ntrials,int nevents,int nt,double xtr,int hemomodel,VBoolean firstcol,
		gsl_matrix *X,gsl_matrix *covariates)
{
  int i,j,jj,k;
  double t,t0,t1,h;

  /* set temporal resolution to 0.5 sec, reduce discretization artefacts */
  double hfactor = floor(xtr/0.5);
  if (hfactor < 1.0) hfactor = 1.0;

  double tr = xtr/hfactor;
  int jfactor = (int)hfactor;
  int ntimesteps = jfactor*nt;

  for(i = 0; i < nevents; i++) {
    double xmin = VRepnMaxValue(VFloatRepn);
    for(j = 0; j < ntrials; j++) {
      if(trial[j].id != i) continue;
      if(trial[j].duration < xmin) xmin = trial[j].duration;
    }
  }

  /* get kernels */
  t1 = 30.0;    /* kernel duration = 30 secs */
  int kernelsize = t1 / tr;
  if (tr < 0.01) VWarning(" implausible TR (%f seconds)",tr);

  /* fprintf(stderr," hf= %f,  tr= %f  %f, kernelsize: %d\n",hfactor,xtr,tr,kernelsize); */

  double *kernel1=NULL,*kernel2=NULL;
  double *bkernel  = (double *)VCalloc(kernelsize,sizeof(double));
  double *kernel0  = (double *)VCalloc(kernelsize,sizeof(double));
  if (hemomodel == 1 || hemomodel == 2) {
    kernel1 = (double *)VCalloc(kernelsize,sizeof(double));
  }
  if (hemomodel == 2) {
    kernel2 = (double *)VCalloc(kernelsize,sizeof(double));
  }

  i = 0;
  for (t = 0; t < t1; t += tr) {
    if (i >= kernelsize) break;

    bkernel[i] = xgauss(t, 5.0*tr);
    kernel0[i] = xgamma(t, 0.0);

    if(hemomodel == 1 || hemomodel == 2)
      kernel1[i] = deriv1_gamma(t, 0);
    if(hemomodel == 2)
      kernel2[i] = deriv2_gamma(t, 0);
    i++;
  }

  /* tmp storage */
  double *x = (double *) VCalloc(ntimesteps,sizeof(double));
  double *y = (double *) VCalloc(ntimesteps,sizeof(double));

  /* constant in column 0 */
  int col=0;
  if (firstcol == TRUE) {
    gsl_matrix_set_zero(X);
    for (j=0; j<nt; j++) gsl_matrix_set(X,j,0,1.0);
    col=1;
  }

  /* for each trial,event, do... */
  for (i = 0; i < nevents; i++) {
    for (k=0; k<ntimesteps; k++) x[k] = y[k] = 0;

    /* read design info */
    int trialcount = 0;
    for (j = 0; j < ntrials; j++) {
      if(trial[j].id != i+1) continue;
      trialcount++;
      t0 = trial[j].onset;
      t1 = trial[j].onset + trial[j].duration;
      h  = trial[j].height;
      int k0 = (int) (t0/tr + 0.5);
      int k1 = (int) (t1/tr + 0.5);

      int klen=0;
      for (k=k0; k<=k1; k++) {
	if (k < 0 || k >= ntimesteps) continue;
	x[k] = h;
	klen++;
      }
      if(trialcount < 1)
	VError(" no trials in event %d, please re-number event-ids,  ntrials= %d", i + 1,ntrials);
    }

    /* convolve */
    if (hemomodel == 3)  {       /* block design, gaussian kernel */
      XConvolve(x,y,bkernel,ntimesteps,kernelsize);
    }
    else {                       /* gamma function kernel */
      XConvolve(x,y,kernel0,ntimesteps,kernelsize);
    }

    for (j=0; j<X->size1; j++) {
      jj = j*jfactor;
      if (jj >= ntimesteps) continue;
      gsl_matrix_set(X,j,col,y[jj]/hfactor);
    }
    col++;

    if (hemomodel == 1 || hemomodel == 2) {   /* gamma function kernel, first derivative */
      XConvolve(x,y,kernel1,ntimesteps,kernelsize);
      for (j=0; j<X->size1; j++) {
	jj = j*jfactor;
	if (jj >= ntimesteps) continue;
	gsl_matrix_set(X,j,col,y[jj]/hfactor);
      }
      col++;
    }
    if (hemomodel == 2) {   /* gamma function kernel, second derivative */
      XConvolve(x,y,kernel2,ntimesteps,kernelsize);
      for (j=0; j<X->size1; j++) {
	jj = j*jfactor;
	if (jj >= ntimesteps) continue;
	gsl_matrix_set(X,j,col,y[jj]/hfactor);
      }
      col++;
    }
  }

  /* add further covariates that have no task labels, will not be permuted */
  if (covariates != NULL) {
    for (i=0; i<covariates->size2; i++) {
      for (j=0; j<covariates->size1; j++) {
	if (col >= X->size2) VError(" VHemoModel, column= %d  %d\n",col,(int)X->size2);
	gsl_matrix_set(X,j,col,gsl_matrix_get(covariates,j,i));
      }
      col++;
    }
  }

  /* cleanup */
  VFree(x);
  VFree(y);
  VFree(kernel0);
  VFree(bkernel);
  if (kernel1 != NULL) VFree(kernel1);
  if (kernel2 != NULL) VFree(kernel2);
}

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_roots.h>
#include <gsl/gsl_math.h>
#include "zero_crossing.h"

/* \fn double zero_crossing(double (*model)(double, void *), void *params, double *x_lo, double *x_hi)
 * \brief Find a zero crossing of the supplied model. */
double zero_crossing(double (*model)(double, void *), void *params, double *x_lo, double *x_hi)
{
	int status;
	int iter=0;
	int max_iter = 100;

	double x;	//root

	const gsl_root_fsolver_type *T;
	gsl_root_fsolver *s;

	gsl_function F;

	F.function = model;
	F.params   = params;

	T = gsl_root_fsolver_brent;
	s = gsl_root_fsolver_alloc(T);

	gsl_root_fsolver_set(s,&F, *x_lo, *x_hi);

	do{
		iter++;
		status = gsl_root_fsolver_iterate(s);
		x      = gsl_root_fsolver_root(s);

		*x_lo   = gsl_root_fsolver_x_lower(s);
		*x_hi   = gsl_root_fsolver_x_upper(s);

		status = gsl_root_test_interval(*x_lo, *x_hi,0,0.001);

		//printf("iter %d x_lo %e x %e x_hi %e\n",iter,*x_lo,x,*x_hi);

	}while((status==GSL_CONTINUE)&&(iter<max_iter));

	gsl_root_fsolver_free(s);

	return x;
}

max_stars_repo_name: brantr/zero-crossing

max_issues_repo_name: brantr/zero-crossing

max_forks_repo_name: brantr/zero-crossing


text:
#include <cblas.h>
#include <clapack.h>
#ifdef NoChange
   #define dgeqrf_ dgeqrf
#elif defined(UpCase)
   #define dgeqrf_ DGEQRF
#endif
main(int nargs, char **args)
{
   extern void dgeqrf_(F77_INTEGER*,F77_INTEGER*,double*,F77_INTEGER*,
                       double*,double*,F77_INTEGER*,F77_INTEGER*);
   double A[1]={1.0}, b[1]={1.0}, c[1];
   F77_INTEGER N=1, info;
   dgeqrf_(&N, &N, A, &N, b, c, &N, &info);
}

max_stars_repo_name: kevleyski/math-atlas

max_issues_repo_name: kevleyski/math-atlas

max_forks_repo_name: kevleyski/math-atlas


text:
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <string.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_histogram.h>

#include "allvars.h"
#include "proto.h"

void get_cov_matrix(double *theta, int nstep, int ntheta)
{
  int i, j;
  double cov;

  for(i=0; i<ntheta; i++)
  {
    cov = gsl_stats_covariance(&theta[i*nstep], 1, &theta[i*nstep], 1, nstep);
    cov = cov>1.0e-6 ? cov : 1.0e-6;
    cov_matrix[i*ntheta + i] = cov;
    for(j=0; j<i; j++)
    {
      cov = gsl_stats_covariance(&theta[i*nstep], 1, &theta[j*nstep], 1, nstep);
      cov = cov>0.0001?cov:0.0;
      cov_matrix[i*ntheta + j] = cov_matrix[j*ntheta + i] = cov;
    }
  }
}
void get_cov_matrix_diag(double *theta, int nstep, int ntheta)
{
  int i, j;
  double cov;

  for(i=0; i<ntheta; i++)
  {
    cov = gsl_stats_covariance(&theta[i*nstep], 1, &theta[i*nstep], 1, nstep);
    cov = cov>1.0e-6 ? cov : 1.0e-6;
    cov_matrix[i*ntheta + i] = cov;
    for(j=0; j<i; j++)
    {
      cov_matrix[i*ntheta + j] = cov_matrix[j*ntheta + i] = 0.0;
    }
  }
}

void mcmc_stats(char * fname)
{
  FILE *fp;
  double **theta, tmp, err_max1, err_max2, theta_up, theta_low;
  int i, nstep, istep;

  theta = malloc(ntheta*sizeof(double));
  for(i=0; i<ntheta; i++)
  {
    theta[i] = malloc((nmcmc)*sizeof(double));
  }

  fp = fopen(fname, "r");
  if(fp==NULL)
  {
    strcpy(str_error_exit, fname);
    error_exit(2);
  }

  istep = 0;
  nstep = 0;
  while(!feof(fp) && nstep < nmcmc)
  {
    fscanf(fp, "%d", &istep);
    for(i=0; i<ntheta; i++)
    {
      fscanf(fp, "%lf", &theta[i][nstep]);
    }
    fscanf(fp,"%lf\n", &tmp);
    nstep++;
  }
  printf("nstep: %d\n", nstep);
  if(nstep <= nbuilt)
  {
    strcpy(str_error_exit, "mcmc.txt");
    error_exit(10);
  }

  printf("Sts: ID    par   err1   err2\n");
  for(i=0; i<ntheta; i++)
  {
    gsl_sort(&theta[i][nbuilt-1], 1, nstep-nbuilt);
    //printf("ss:%f %f %f\n", theta[i*n_mcmc + 40000], theta[i*n_mcmc+ 40100], theta[i*n_mcmc+ 40200]);
    theta_best[i] = gsl_stats_mean(&theta[i][nbuilt-1], 1, nstep-nbuilt);
    //theta_best[i] = gsl_stats_quantile_from_sorted_data(&theta[i][nbuilt-1], 1, nstep-nbuilt, 0.500);
    theta_low = gsl_stats_quantile_from_sorted_data(&theta[i][nbuilt-1], 1, nstep-nbuilt, 0.1585);
    theta_up = gsl_stats_quantile_from_sorted_data(&theta[i][nbuilt-1], 1, nstep-nbuilt, 1.0-0.1585);

    err_max1 = (theta_best[i] - theta_range[i][0])/2.35;
    err_max2 = (theta_range[i][1] - theta_best[i])/2.35;

    theta_best_var[i*2] = max(min( theta_best[i] - theta_low, err_max1 ), 0.0);
    theta_best_var[i*2+1] = max( min( theta_up - theta_best[i], err_max2 ), 0.0);

    printf("Sts: %d %f %f %f\n", i, theta_best[i], theta_best_var[i*2], theta_best_var[i*2+1]);

  }

  fclose(fp);
  return;
}

max_stars_repo_name: LiyrAstroph/MICA

max_issues_repo_name: LiyrAstroph/MICA

max_forks_repo_name: LiyrAstroph/MICA


text:
/**
 * @file beamformer.h
 * @brief Beamforming in the subband domain.
 * @author John McDonough and Kenichi Kumatani
 */
#ifndef BEAMFORMER_H
#define BEAMFORMER_H

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

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

#include "stream/stream.h"
#include "beamformer/spectralinfoarray.h"
#include "modulated/modulated.h"

#define SSPEED 343740.0

class BeamformerWeights {
public:
  BeamformerWeights( unsigned fftLen, unsigned chanN, bool halfBandShift, unsigned NC = 1 );
  ~BeamformerWeights();

  void calcMainlobe(  float samplerate, const gsl_vector* delays,  bool isGSC );
  void calcMainlobe2( float samplerate, const gsl_vector* delaysT, const gsl_vector* delaysJ, bool isGSC );
  void calcMainlobeN( float samplerate, const gsl_vector* delaysT, const gsl_matrix* delaysIs, unsigned NC, bool isGSC );
  void calcSidelobeCancellerP_f( unsigned fbinX, const gsl_vector* packedWeight );
  void calcSidelobeCancellerU_f( unsigned fbinX, const gsl_vector_complex* wa );
  void calcBlockingMatrix( unsigned fbinX );

  bool write_fir_coeff(const String& fn, unsigned winType);

#ifdef ENABLE_LEGACY_BTK_API
  bool writeFIRCoeff(const String& fn, unsigned winType);
#endif

  void setSidelobeCanceller_f( unsigned fbinX, gsl_vector_complex* wl_f ){
    gsl_vector_complex_memcpy( wl_[fbinX], wl_f );
  }
  void setQuiescentVector( unsigned fbinX, gsl_vector_complex *wq_f, bool isGSC=false );
  void setQuiescentVectorAll( gsl_complex z, bool isGSC=false );
  void setTimeAlignment();

  bool     isHalfBandShift() const {return(halfBandShift_);}
  unsigned NC() const {return(NC_);}
  unsigned fftLen() const {return(fftLen_);}
  unsigned chanN() const {return(chanN_);}

  gsl_vector_complex** arrayManifold() const { return (ta_); }
  gsl_vector_complex* wq_f( unsigned fbinX ) const { return wq_[fbinX]; }
  gsl_vector_complex* wl_f( unsigned fbinX ) const { return wl_[fbinX]; }
  gsl_vector_complex** wq() const { return (wq_); }
  gsl_matrix_complex** B() const { return (B_); }
  gsl_vector_complex** wa() const { return (wa_); }
  gsl_vector_complex** CSDs() const { return CSDs_; }
  gsl_vector_complex* wp1() const { return wp1_; }

private:
  void alloc_weights_();
  void free_weights_();

  unsigned fftLen_;
  unsigned chanN_;
  bool halfBandShift_;
  unsigned NC_; // the numbef of constraints
  gsl_vector_complex** wq_; // a quiescent weight vector for each frequency bin, wq_[fbinX][chanN]
  gsl_matrix_complex** B_;  // a blocking matrix for each frequency bin,         B_[fbinX][chanN][chanN-NC]
  gsl_vector_complex** wa_; // an active weight vector for each frequency bin,   wa_[fbinX][chanN-NC]
  gsl_vector_complex** wl_; // wl_[fbinX] = B_[fbinX] * wa_[fbinX]
  gsl_vector_complex** ta_; // do time alignment for multi-channel waves. It is also called an array manifold. _ta[fbinX][chanN].
  gsl_vector_complex*  wp1_;  // a weight vector of postfiltering,   _wp[fbinX]
  gsl_vector_complex** CSDs_; // cross spectral density for the post-filtering
};

typedef refcount_ptr<BeamformerWeights>     BeamformerWeightsPtr;

// ----- definition for class `SubbandBeamformer' -----
//
class SubbandBeamformer : public VectorComplexFeatureStream {
 public:
  SubbandBeamformer(unsigned fftLen = 512, bool halfBandShift = false, const String& nm = "SubbandBeamformer");
  ~SubbandBeamformer();

  virtual const gsl_vector_complex* next(int frame_no = -5);
  virtual void reset();

  unsigned fftLen() const { return fftLen_; }
  unsigned fftLen2() const { return fftLen2_; }
  unsigned chanN() const { return channelList_.size(); }
  virtual unsigned dim() const { return chanN();}

  bool is_end() const {return is_end_;}
  const gsl_vector_complex* snapshot_array_f(unsigned fbinX) const { return (snapshot_array_->snapshot(fbinX)); }
  virtual SnapShotArrayPtr  snapshot_array() const { return(snapshot_array_); }
  void         set_channel(VectorComplexFeatureStreamPtr& chan);
  virtual void clear_channel();

#ifdef ENABLE_LEGACY_BTK_API
  bool isEnd() { return is_end(); }
  const gsl_vector_complex* snapShotArray_f(unsigned fbinX){ return snapshot_array_f(fbinX); }
  virtual SnapShotArrayPtr  getSnapShotArray(){ return(snapshot_array()); }
  void         setChannel(VectorComplexFeatureStreamPtr& chan){ set_channel(chan); }
  virtual void clearChannel(){ clear_channel(); }
#endif

protected:
  typedef list<VectorComplexFeatureStreamPtr>	ChannelList_;
  typedef ChannelList_::iterator		ChannelIterator_;

  SnapShotArrayPtr				snapshot_array_;
  unsigned					fftLen_;
  unsigned					fftLen2_;
  bool						halfBandShift_;
  ChannelList_					channelList_;
};

// ----- definition for class `SubbandDS' -----
//

class SubbandDS : public SubbandBeamformer {
 public:
  SubbandDS(unsigned fftLen = 512, bool halfBandShift = false, const String& nm = "SubbandDS");
  ~SubbandDS();

  virtual const gsl_vector_complex* next(int frame_no = -5);
  virtual void reset();
  virtual void clear_channel();
  virtual const gsl_vector_complex *get_weights(unsigned fbinX) const { return bfweight_vec_[0]->wq_f(fbinX); }
  virtual BeamformerWeights* beamformer_weight_object(unsigned srcX=0) const { return bfweight_vec_[srcX]; }

  virtual void calc_array_manifold_vectors(float samplerate, const gsl_vector* delays);
  virtual void calc_array_manifold_vectors_2(float samplerate, const gsl_vector* delaysT, const gsl_vector* delaysJ);
  virtual void calc_array_manifold_vectors_n(float samplerate, const gsl_vector* delaysT, const gsl_matrix* delaysJ, unsigned NC=2);

#ifdef ENABLE_LEGACY_BTK_API
  virtual void clearChannel(){ clear_channel(); }
  virtual const gsl_vector_complex *getWeights(unsigned fbinX) const { return get_weights(fbinX); }
  virtual BeamformerWeights* getBeamformerWeightObject(unsigned srcX=0) const { return beamformer_weight_object(srcX); }
  virtual void calcArrayManifoldVectors(float sampleRate, const gsl_vector* delays){
    calc_array_manifold_vectors(sampleRate, delays);
  }
  virtual void calcArrayManifoldVectors2(float sampleRate, const gsl_vector* delaysT, const gsl_vector* delaysJ){
    calc_array_manifold_vectors_2(sampleRate, delaysT, delaysJ);
  }
  virtual void calcArrayManifoldVectorsN(float sampleRate, const gsl_vector* delaysT, const gsl_matrix* delaysJ, unsigned NC=2){
    calc_array_manifold_vectors_n(sampleRate, delaysT, delaysJ, NC);
  }
#endif /* #ifdef ENABLE_LEGACY_BTK_API */

protected:
  void alloc_image_();
  void alloc_bfweight_(int nSrc, int NC);

  vector<BeamformerWeights *>                   bfweight_vec_; // weights of a beamformer per source.
};

#define NO_PROCESSING 0x00
#define SCALING_MDP   0x01
class SubbandGSC : public SubbandDS {
public:
  SubbandGSC(unsigned fftLen = 512, bool halfBandShift = false, const String& nm = "SubbandGSC")
    : SubbandDS( fftLen, halfBandShift, nm ),normalize_weight_(false){}
  ~SubbandGSC();

  virtual const gsl_vector_complex* next(int frame_no = -5);

  void normalize_weight(bool flag){ normalize_weight_ = flag; }
  void set_quiescent_weights_f(unsigned fbinX, const gsl_vector_complex* srcWq);
  void set_active_weights_f(unsigned fbinX, const gsl_vector* packedWeight);
  void zero_active_weights();
  void calc_gsc_weights(float samplerate, const gsl_vector* delaysT);
  void calc_gsc_weights_2(float samplerate, const gsl_vector* delaysT, const gsl_vector* delaysJ);
  void calc_gsc_weights_n(float samplerate, const gsl_vector* delaysT, const gsl_matrix* delaysJ, unsigned NC=2);

  bool write_fir_coeff(const String& fn, unsigned winType=1);
  gsl_matrix_complex* blocking_matrix(unsigned srcX, unsigned fbinX){
    return (bfweight_vec_[srcX]->B())[fbinX];
  }

#ifdef ENABLE_LEGACY_BTK_API
  void normalizeWeight(bool flag){ normalize_weight(flag); }
  void setQuiescentWeights_f(unsigned fbinX, const gsl_vector_complex * srcWq){ set_quiescent_weights_f(fbinX, srcWq); }
  void setActiveWeights_f(unsigned fbinX, const gsl_vector* packedWeight){ set_active_weights_f(fbinX, packedWeight); }
  void zeroActiveWeights(){ zero_active_weights(); }
  void calcGSCWeights(float sampleRate, const gsl_vector* delaysT){ calc_gsc_weights(sampleRate, delaysT); }
  void calcGSCWeights2(float sampleRate, const gsl_vector* delaysT, const gsl_vector* delaysJ){ calc_gsc_weights_2(sampleRate, delaysT, delaysJ); }
  void calcGSCWeightsN(float sampleRate, const gsl_vector* delaysT, const gsl_matrix* delaysJ, unsigned NC=2){ calc_gsc_weights_n(sampleRate, delaysT, delaysJ, NC); }
  bool writeFIRCoeff(const String& fn, unsigned winType=1){ return write_fir_coeff(fn, winType); }
  gsl_matrix_complex* getBlockingMatrix(unsigned srcX, unsigned fbinX){ return blocking_matrix(srcX, fbinX); }
#endif /* #ifdef ENABLE_LEGACY_BTK_API */

protected:
  bool normalize_weight_;
};

/**
   @class SubbandGSCRLS
   @brief implementation of recursive least squares of a GSC
   @usage
   1. calcGSCWeights()
   2. initPrecisionMatrix() or setPrecisionMatrix()
   3. update_sctive_weight_vecotrs( false ) if you want to stop adapting the active weight vectors.
   @note notations are  based on Van Trees, "Optimum Array Processing", pp. 766-767.
 */
typedef enum {
  CONSTANT_NORM           = 0x01,
  THRESHOLD_LIMITATION    = 0x02,
  NO_QUADRATIC_CONSTRAINT = 0x00
} QuadraticConstraintType;

// ----- definition for class `SubbandGSCRLS' -----
//

class SubbandGSCRLS : public SubbandGSC {
 public:
  SubbandGSCRLS(unsigned fftLen = 512, bool halfBandShift = false, float mu = 0.9, float sigma2=0.0, const String& nm = "SubbandGSCRLS");
  ~SubbandGSCRLS();

  virtual const gsl_vector_complex* next(int frame_no = -5);
  virtual void reset();

  void init_precision_matrix(float sigma2 = 0.01);
  void set_precision_matrix(unsigned fbinX, gsl_matrix_complex *Pz);
  void update_active_weight_vecotrs(bool flag){ is_wa_updated_ = flag; }
  void set_quadratic_constraint(float alpha, int qctype=1){ alpha_=alpha; qctype_=(QuadraticConstraintType)qctype; }

#ifdef ENABLE_LEGACY_BTK_API
  void initPrecisionMatrix(float sigma2 = 0.01){ init_precision_matrix(sigma2); }
  void setPrecisionMatrix(unsigned fbinX, gsl_matrix_complex *Pz){ set_precision_matrix(fbinX, Pz); }
  void updateActiveWeightVecotrs(bool flag){ update_active_weight_vecotrs(flag); }
  void setQuadraticConstraint(float alpha, int qctype=1){ set_quadratic_constraint(alpha, qctype); }
#endif /* #ifdef ENABLE_LEGACY_BTK_API */

private:
  void update_active_weight_vector2_(int frame_no); /* the case of the half band shift = False */
  bool alloc_subbandGSCRLS_image_();
  void free_subbandGSCRLS_image_();

  gsl_vector_complex** gz_; /* Gain vectors */
  gsl_matrix_complex** Pz_; /* Precision matrices */
  gsl_vector_complex* Zf_;  /* output of the blocking matrix at each frequency */
  gsl_vector_complex* wa_;
  float  mu_;              /* Exponential factor for the covariance matrix */
  float* diagonal_weights_;
  float  alpha_;            /* Weight for the quadratic constraint*/
  QuadraticConstraintType qctype_;
  bool is_wa_updated_;

  /* work space for updating active weight vectors */
  gsl_vector_complex* PzH_Z_;
  gsl_matrix_complex* _I;
  gsl_matrix_complex* mat1_;
};

// ----- definition for class `SubbandMMI' -----
//

class SubbandMMI : public SubbandDS {
public:
  SubbandMMI(unsigned fftLen = 512, bool halfBandShift = false, unsigned targetSourceX=0, unsigned nSource=2, int pfType=0, float alpha=0.9, const String& nm = "SubbandMMI")
    : SubbandDS( fftLen, halfBandShift, nm ),
    targetSourceX_(targetSourceX),
    nSource_(nSource),
    pftype_(pfType),
    alpha_(alpha),
    use_binary_mask_(false),
    binary_mask_type_(0),
    interference_outputs_(NULL),
    avg_output_(NULL)
  {}

  ~SubbandMMI();

  virtual const gsl_vector_complex* next(int frame_no = -5);

  void use_binary_mask(float avgFactor=-1.0, unsigned fwidth=1, unsigned type=0);
  void calc_weights(  float samplerate, const gsl_matrix* delays);
  void calc_weights_n( float samplerate, const gsl_matrix* delays, unsigned NC=2);
  void set_hi_active_weights_f(unsigned fbinX, const gsl_vector* pkdWa, const gsl_vector* pkdwb, int option=0);
  void set_active_weights_f(unsigned fbinX, const gsl_matrix* packedWeights, int option=0);

#ifdef ENABLE_LEGACY_BTK_API
  void useBinaryMask(float avgFactor=-1.0, unsigned fwidth=1, unsigned type=0){ use_binary_mask(avgFactor, fwidth, type); }
  void calcWeights(  float sampleRate, const gsl_matrix* delays){ calc_weights(sampleRate, delays); }
  void calcWeightsN( float sampleRate, const gsl_matrix* delays, unsigned NC=2){ calc_weights_n(sampleRate, delays, NC); }
  void setHiActiveWeights_f(unsigned fbinX, const gsl_vector* pkdWa, const gsl_vector* pkdwb, int option=0){
    set_hi_active_weights_f(fbinX, pkdWa, pkdwb, option);
  }
  void setActiveWeights_f(unsigned fbinX, const gsl_matrix* packedWeights, int option=0){
    set_active_weights_f(fbinX, packedWeights, option);
  }
#endif /* #ifdef ENABLE_LEGACY_BTK_API */

private:
  void calc_interference_outputs_();
  void binary_masking_( gsl_vector_complex** interferenceOutputs, unsigned targetSourceX, gsl_vector_complex* output );

  unsigned                                      targetSourceX_; // the n-th source will be emphasized
  unsigned                                      nSource_;       // the number of sound sources
  int                                           pftype_;
  float                                        alpha_;
  bool                                          use_binary_mask_; // true if you use a binary mask
  unsigned                                      binary_mask_type_;// 0:use GSC's outputs, 1:use outputs of the upper branch.
  gsl_vector_complex**                          interference_outputs_;
  gsl_vector_complex*                           avg_output_;
  float                                        avg_factor_;
  unsigned                                      fwidth_;
};

// ----- definition for class `SubbandMVDR' -----
//

/**
   @class SubbandMVDR

   @usage
   1. setChannel()
   2. calc_array_manifold_vectors(), calc_array_manifold_vectors2() or calc_array_manifold_vectorsN().
   3. set_noise_spatial_spectral_matrix() or set_diffuse_noise_model()
   4. calc_mvdr_weights()
 */
class SubbandMVDR : public SubbandDS {
 public:
  /**
     @brief Basic MVDR beamformer implementation
     @param int fftLen[in]
     @param bool halfBandShift[in]
   */
  SubbandMVDR(unsigned fftLen = 512, bool halfBandShift = false, const String& nm = "SubbandMVDR");
  ~SubbandMVDR();

  virtual const gsl_vector_complex* next(int frame_no = -5);
  virtual void clear_channel();
  bool calc_mvdr_weights(float samplerate, float dThreshold = 1.0E-8, bool calcInverseMatrix = true);
  const gsl_vector_complex* mvdr_weights(unsigned fbinX) const { return wmvdr_[fbinX]; }

  const gsl_matrix_complex *noise_spatial_spectral_matrix(unsigned fbinX) const { return R_[fbinX]; }
  bool set_noise_spatial_spectral_matrix(unsigned fbinX, gsl_matrix_complex* Rnn);
  bool set_diffuse_noise_model(const gsl_matrix* micPositions, float samplerate, float sspeed = 343740.0); /* micPositions[][x,y,z] */
  void set_all_diagonal_loading(float diagonalWeight);
  void set_diagonal_looading(unsigned fbinX, float diagonalWeight);
  /**
     @brief Divide each non-diagonal elemnt by 1 + mu instead of diagonal loading. mu can be interpreted as the ratio of the sensor noise to the ambient noise power.
     @param float mu[in]
   */
  void divide_all_nondiagonal_elements(float mu){
    for(unsigned fbinX=0;fbinX<=fftLen_/2;fbinX++)
      divide_nondiagonal_elements( fbinX, mu );
  }
  void divide_nondiagonal_elements(unsigned fbinX, float mu);
  gsl_matrix_complex**  noise_spatial_spectral_matrix() const { return R_; }

#ifdef ENABLE_LEGACY_BTK_API
  void clearChannel(){ clear_channel(); }
  bool calcMVDRWeights( float sampleRate, float dThreshold = 1.0E-8, bool calcInverseMatrix = true ){ return calc_mvdr_weights(sampleRate, dThreshold, calcInverseMatrix); }
  const gsl_vector_complex* getMVDRWeights(unsigned fbinX){ return mvdr_weights(fbinX); }
  const gsl_matrix_complex *getNoiseSpatialSpectralMatrix(unsigned fbinX){ return noise_spatial_spectral_matrix(fbinX); }
  bool setNoiseSpatialSpectralMatrix(unsigned fbinX, gsl_matrix_complex* Rnn){ return set_noise_spatial_spectral_matrix(fbinX, Rnn); }
  bool setDiffuseNoiseModel(const gsl_matrix* micPositions, float sampleRate, float sspeed = 343740.0){ return set_diffuse_noise_model(micPositions, sampleRate, sspeed); }
  void setAllLevelsOfDiagonalLoading(float diagonalWeight){ set_all_diagonal_loading(diagonalWeight); }
  void setLevelOfDiagonalLoading(unsigned fbinX, float diagonalWeight){ set_diagonal_looading(fbinX, diagonalWeight); }
  void divideAllNonDiagonalElements( float mu ){ divide_all_nondiagonal_elements(mu); }
  void divideNonDiagonalElements( unsigned fbinX, float mu ){ divide_nondiagonal_elements(fbinX, mu); }
  gsl_matrix_complex**  getNoiseSpatialSpectralMatrix(){ return noise_spatial_spectral_matrix(); }
#endif /* #ifdef ENABLE_LEGACY_BTK_API */

protected:
  gsl_matrix_complex**                           R_; /* Noise spatial spectral matrices */
  gsl_matrix_complex**                           invR_;
  gsl_vector_complex**                           wmvdr_;
  float*                                         diagonal_weights_;
};

// ----- definition for class `SubbandMVDRGSC' -----
//

/**
   @class SubbandMVDRGSC

   @usage
   1. setChannel()
   2. calc_array_manifold_vectors(), calc_array_manifold_vectors2() or calc_array_manifold_vectorsN().
   3. set_noise_spatial_spectral_matrix() or set_diffuse_noise_model()
   4. calc_mvdr_weights()
   5. calc_blocking_matrix1() or calc_blocking_matrix2()
   6. set_active_weights_f()
 */
class SubbandMVDRGSC : public SubbandMVDR {
 public:
  /**
     @brief MVDR beamforming implementation
     @param int fftLen[in]
     @param bool halfBandShift[in]
   */
  SubbandMVDRGSC(unsigned fftLen = 512, bool halfBandShift = false, const String& nm = "SubbandMVDR");
  ~SubbandMVDRGSC();

  virtual const gsl_vector_complex* next(int frame_no = -5);

  void set_active_weights_f(unsigned fbinX, const gsl_vector* packedWeight);
  void zero_active_weights();
  bool calc_blocking_matrix1(float samplerate, const gsl_vector* delaysT);
  bool calc_blocking_matrix2();
  void upgrade_blocking_matrix();
  const gsl_vector_complex* blocking_matrix_output(int outChanX=0);

#ifdef ENABLE_LEGACY_BTK_API
  void setActiveWeights_f(unsigned fbinX, const gsl_vector* packedWeight){ set_active_weights_f(fbinX, packedWeight); }
  void zeroActiveWeights(){ zero_active_weights(); }
  bool calcBlockingMatrix1(float sampleRate, const gsl_vector* delaysT){ return calc_blocking_matrix1(sampleRate, delaysT); }
  bool calcBlockingMatrix2(){ return calc_blocking_matrix2(); }
  void upgradeBlockingMatrix(){ upgrade_blocking_matrix(); }
  const gsl_vector_complex* blockingMatrixOutput(int outChanX=0){ return blocking_matrix_output(outChanX); }
#endif

protected:
  bool normalize_weight_;
};

typedef Inherit<SubbandBeamformer, VectorComplexFeatureStreamPtr> SubbandBeamformerPtr;
typedef Inherit<SubbandDS, SubbandBeamformerPtr> SubbandDSPtr;
typedef Inherit<SubbandGSC, SubbandDSPtr> SubbandGSCPtr;
typedef Inherit<SubbandGSCRLS, SubbandGSCPtr> SubbandGSCRLSPtr;
typedef Inherit<SubbandMMI, SubbandDSPtr> SubbandMMIPtr;
typedef Inherit<SubbandMVDR, SubbandDSPtr> SubbandMVDRPtr;
typedef Inherit<SubbandMVDRGSC, SubbandMVDRPtr> SubbandMVDRGSCPtr;

// ----- members for class `SubbandOrthogonalizer' -----
//

class SubbandOrthogonalizer : public VectorComplexFeatureStream {
public:
  SubbandOrthogonalizer(SubbandMVDRGSCPtr &beamformer, int outChanX=0, const String& nm = "SubbandOrthogonalizer");
  ~SubbandOrthogonalizer();
  virtual const gsl_vector_complex* next(int frame_no = -5);

private:
  SubbandMVDRGSCPtr beamformer_;
  int outChanX_;
};

typedef Inherit<SubbandOrthogonalizer, VectorComplexFeatureStreamPtr> SubbandOrthogonalizerPtr;

class SubbandBlockingMatrix : public SubbandGSC {
public:
  SubbandBlockingMatrix(unsigned fftLen=512, bool halfBandShift=false, const String& nm = "SubbandBlockingMatrix")
    :SubbandGSC(fftLen, halfBandShift, nm ){;}

  ~SubbandBlockingMatrix();
  virtual const gsl_vector_complex* next(int frame_no = -5);
};

// ----- definition for class DOAEstimatorSRPBase' -----
//
class DOAEstimatorSRPBase {
public:
  DOAEstimatorSRPBase( unsigned nBest, unsigned fbinMax );
  virtual ~DOAEstimatorSRPBase();

  const gsl_vector *nbest_rps() const { return nBestRPs_; }
  const gsl_matrix *nbest_doas() const { return argMaxDOAs_;}
  const gsl_matrix *response_power_matrix() const { return rpMat_;}
  float energy() const {return energy_;}
  void  final_nbest_hypotheses(){get_nbest_hypotheses_from_accrp_();}
  void  set_energy_threshold(float engeryThreshold){ engery_threshold_ = engeryThreshold; }
  void  set_frequency_range(unsigned fbinMin, unsigned fbinMax){ fbinMin_ = fbinMin; fbinMax_ = fbinMax;}
  void  init_accs(){ init_accs_(); }
  void  set_search_param(float minTheta=-M_PI/2, float maxTheta=M_PI/2,
                         float minPhi=-M_PI/2,   float maxPhi=M_PI/2,
                         float widthTheta=0.1,   float widthPhi=0.1);

#ifdef ENABLE_LEGACY_BTK_API
  const gsl_vector *getNBestRPs(){ return nbest_rps(); }
  const gsl_matrix *getNBestDOAs(){ return nbest_doas(); }
  const gsl_matrix *getResponsePowerMatrix(){ return response_power_matrix(); }
  float getEnergy(){return energy();}
  void getFinalNBestHypotheses(){ final_nbest_hypotheses(); }
  void setEnergyThreshold(float engeryThreshold){ set_energy_threshold(engeryThreshold); }
  void setFrequencyRange(unsigned fbinMin, unsigned fbinMax){ set_frequency_range(fbinMin, fbinMax); }
  void initAccs(){ init_accs(); }
  void setSearchParam(float minTheta=-M_PI/2, float maxTheta=M_PI/2,
                      float minPhi=-M_PI/2,   float maxPhi=M_PI/2,
                      float widthTheta=0.1, float widthPhi=0.1)
  {
    set_search_param(minTheta, maxTheta, minPhi, maxPhi, widthTheta, widthPhi);
  }
#endif

protected:
  void clear_table_();
  virtual void get_nbest_hypotheses_from_accrp_();
  virtual void init_accs_();

  float widthTheta_;
  float widthPhi_;
  float minTheta_;
  float maxTheta_;
  float minPhi_;
  float maxPhi_;
  unsigned nTheta_;
  unsigned nPhi_;
  unsigned fbinMin_;
  unsigned fbinMax_;
  unsigned nBest_;
  bool   table_initialized_;

  gsl_vector *accRPs_;
  gsl_vector *nBestRPs_;
  gsl_matrix *argMaxDOAs_;
  vector<gsl_vector_complex **> svTbl_; // [][fftL2+1][_dim]
  gsl_matrix         *rpMat_;

  float engery_threshold_;
  float energy_;

#ifdef  __MBDEBUG__
  void allocDebugWorkSapce();
#endif /* #ifdef __MBDEBUG__ */
};

// ----- definition for class DOAEstimatorSRPDSBLA' -----
//
/**
   @brief estimate the direction of arrival based on the maximum steered response power

   @usage
   1. construct an object
   2. set the geometry of the linear array
   3. call next()
 */
class DOAEstimatorSRPDSBLA :
  public DOAEstimatorSRPBase, public SubbandDS {
public:
  DOAEstimatorSRPDSBLA( unsigned nBest, unsigned samplerate, unsigned fftLen, const String& nm="DOAEstimatorSRPDSBLA" );
  ~DOAEstimatorSRPDSBLA();

  const gsl_vector_complex* next(int frame_no = -5);
  void reset();

  void set_array_geometry(gsl_vector *positions);

#ifdef ENABLE_LEGACY_BTK_API
  void setArrayGeometry(gsl_vector *positions){ set_array_geometry(positions); }
#endif

protected:
  virtual void   calc_steering_unit_table_();
  virtual float calc_response_power_( unsigned uttX );

private:
  virtual void set_look_direction_( int nChan, float theta );

  unsigned    samplerate_;
  gsl_matrix *arraygeometry_; // [micX][x,y,z]
};

typedef refcount_ptr<DOAEstimatorSRPBase> DOAEstimatorSRPBasePtr;
typedef Inherit<DOAEstimatorSRPDSBLA, SubbandDSPtr> DOAEstimatorSRPDSBLAPtr;

// ----- definition for functions' -----
//

float calc_energy(SnapShotArrayPtr snapShotArray, unsigned fbinMin, unsigned fbinMax, unsigned fftLen2, bool  halfBandShift=false);

void calc_gsc_output(const gsl_vector_complex* snapShot,
                     gsl_vector_complex* wl_f, gsl_vector_complex* wq_f,
                     gsl_complex *pYf, bool normalizeWeight=false );

bool pseudoinverse( gsl_matrix_complex *A, gsl_matrix_complex *invA, float dThreshold =  1.0E-8 );

void calc_all_delays(float x, float y, float z, const gsl_matrix* mpos, gsl_vector* delays);

void calc_product(gsl_vector_complex* synthesisSamples, gsl_matrix_complex* gs_W, gsl_vector_complex* product);

#ifdef ENABLE_LEGACY_BTK_API
inline void calcAllDelays(float x, float y, float z, const gsl_matrix* mpos, gsl_vector* delays)
{
  calc_all_delays(x, y, z, mpos, delays);
}

inline void calcProduct(gsl_vector_complex* synthesisSamples, gsl_matrix_complex* gs_W, gsl_vector_complex* product)
{
  calc_product(synthesisSamples, gs_W, product);
}
#endif

#endif

ext: h


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

int main( int argc, const char *argv[]  ) {
        double x = 5.0 ;
            double y = gsl_sf_bessel_J0( x  ) ;
                printf("J0(%g) = % .18e\n", x, y ) ;
                    return 0 ;

}

max_stars_repo_name: subramon/qlu

max_issues_repo_name: subramon/qlu

max_forks_repo_name: subramon/qlu


text:
/* rng/transputer.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 James Theiler, Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 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 <gsl/gsl_rng.h>

/* This is the INMOS Transputer Development System generator. The sequence is,

   x_{n+1} = (a x_n) mod m

   with a = 1664525 and m = 2^32. The seed specifies the initial
   value, x_1.

   The theoretical value of x_{10001} is 1244127297.

   The period of this generator is 2^30. */

static inline unsigned long int transputer_get (void *vstate);
static double transputer_get_double (void *vstate);
static void transputer_set (void *state, unsigned long int s);

typedef struct
  {
    unsigned long int x;
  }
transputer_state_t;

static unsigned long int
transputer_get (void *vstate)
{
  transputer_state_t *state = (transputer_state_t *) vstate;

  state->x = (1664525 * state->x) & 0xffffffffUL;

  return state->x;
}

static double
transputer_get_double (void *vstate)
{
  return transputer_get (vstate) / 4294967296.0 ;
}

static void
transputer_set (void *vstate, unsigned long int s)
{
  transputer_state_t *state = (transputer_state_t *) vstate;

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

  state->x = s;

  return;
}

static const gsl_rng_type transputer_type =
{"transputer",                          /* name */
 0xffffffffUL,                  /* RAND_MAX */
 1,                             /* RAND_MIN */
 sizeof (transputer_state_t),
 &transputer_set,
 &transputer_get,
 &transputer_get_double};

const gsl_rng_type *gsl_rng_transputer = &transputer_type;

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/* Authors:
   G. Jungman
   J. Scott (james.scott@lexifi.com)
 */
/* Implementation for Sobol generator.
 * See
 *   [Bratley+Fox, TOMS 14, 88 (1988)]
 *   [Antonov+Saleev, USSR Comput. Maths. Math. Phys. 19, 252 (1980)]
 */
#include <config.h>
#include <gsl/gsl_qrng.h>
#include <gsl/gsl_rng.h>

#define SOBOL_BIT_COUNT 30

/* prototypes for generator type functions */
static size_t sobol_state_size(unsigned int dimension);
static int sobol_init(void * state, unsigned int dimension);
static int sobol_get(void * state, unsigned int dimension, double * v);

/* global Sobol generator type object */
static const gsl_qrng_type sobol_type =
{
  "sobol",
  0,
  sobol_state_size,
  sobol_init,
  sobol_get
};
const gsl_qrng_type * gsl_qrng_sobol = &sobol_type;

/* Sobol generator state.
 *   sequence_count       = number of calls with this generator
 *   last_numerator_vec   = last generated numerator vector
 *   last_denominator_inv = 1/denominator for last numerator vector
 *   v_direction          = direction number table
 */
typedef struct
{
  unsigned int  sequence_count;
  double        last_denominator_inv;
  int           *last_numerator_vec;
  int           *v_direction[SOBOL_BIT_COUNT];
} sobol_state_t;

static size_t sobol_state_size(unsigned int dimension)
{
  return
    sizeof(sobol_state_t) +      /* The struct */
    sizeof(int) * dimension +    /* for last_numerator_vec */
    sizeof(int) * dimension * SOBOL_BIT_COUNT; /* for the direction no.s */
}

/* l is a degree number, between 1 and the degree of the polynomial
   associated with the dimension being initialized */
static int sobol_init_direction(gsl_rng *r, int l)
{
  /* See Peter Jaeckel, Monte Carlo Methods in Finance, Wiley 2002, p86.
   */
  int wkl = gsl_rng_uniform_int(r, 1<<l) | 1;
  return wkl;
}

void get_primitive_polynomials(int dimension, int *degree_table, int *primitive_polynomials);

static int sobol_init(void * state, unsigned int dimension)
{
  sobol_state_t * s_state = (sobol_state_t *) state;
  unsigned int i_dim;
  int j, k;
  int ell;
  int *degree_table;
  int *primitive_polynomials;
  int *includ;
  int max_degree;
  gsl_rng *r = gsl_rng_alloc(gsl_rng_default);

  if(dimension < 1) {
    return GSL_EINVAL;
  }

  { /* initialize degree_table, primitive_polynomials, max_degree and includ */
    degree_table = (int *) malloc(sizeof(int) * dimension);
    if (degree_table==0)
      GSL_ERROR_NULL ("allocation of degree table failed for sobol init", GSL_ENOMEM);
    primitive_polynomials = (int *) malloc(sizeof(int) * dimension);
    if (primitive_polynomials==0) {
      free(degree_table);
      GSL_ERROR_NULL ("allocation of primitives failed for sobol init", GSL_ENOMEM);
    }

    /* Generate the primitive polynomials */
    get_primitive_polynomials(dimension, degree_table, primitive_polynomials);

    max_degree = degree_table[dimension-1];
    includ = (int *) malloc(sizeof(int) * max_degree);
    if (includ==0) {
      free(degree_table);
      free(primitive_polynomials);
      GSL_ERROR_NULL ("allocation of 'includ' failed for sobol init", GSL_ENOMEM);
    }
  }

  s_state->last_numerator_vec = (int *) ((char *) s_state + sizeof(sobol_state_t));

  /* Initialize direction table in dimension 0. */
  for(k=0; k<SOBOL_BIT_COUNT; k++) {
    s_state->v_direction[k] =
      (int *) (((char *) s_state) +
	       sizeof(sobol_state_t) +
	       sizeof(int) * dimension +
	       sizeof(int) * dimension * k);
    s_state->v_direction[k][0] = 1;
  }

  /* Initialize in remaining dimensions. */
  for(i_dim=1; i_dim<dimension; i_dim++) {

    const int poly_index = i_dim;
    const int degree_i = degree_table[poly_index];

    /* Expand the polynomial bit pattern to separate
     * components of the logical array includ[].
     */
    int p_i = primitive_polynomials[poly_index];
    for(k = degree_i-1; k >= 0; k--) {
      includ[k] = ((p_i % 2) == 1);
      p_i /= 2;
    }

    /* Leading elements for dimension i are randomly initialized */
    for(j=0; j<degree_i; j++) s_state->v_direction[j][i_dim] = sobol_init_direction(r, j+1);

    /* Calculate remaining elements for this dimension,
     * as explained in Bratley+Fox, section 2.
     */
    for(j=degree_i; j<SOBOL_BIT_COUNT; j++) {
      int newv = s_state->v_direction[j-degree_i][i_dim];
      ell = 1;
      for(k=0; k<degree_i; k++) {
        ell *= 2;
        if(includ[k]) newv ^= (ell * s_state->v_direction[j-k-1][i_dim]);
      }
      s_state->v_direction[j][i_dim] = newv;
    }
  }

  /* Multiply columns of v by appropriate power of 2. */
  ell = 1;
  for(j=SOBOL_BIT_COUNT-1-1; j>=0; j--) {
    ell *= 2;
    for(i_dim=0; i_dim<dimension; i_dim++) {
      s_state->v_direction[j][i_dim] *= ell;
    }
  }

  /* 1/(common denominator of the elements in v_direction) */
  s_state->last_denominator_inv = 1.0 /(2.0 * ell);

  /* final setup */
  s_state->sequence_count = 0;
  for(i_dim=0; i_dim<dimension; i_dim++) s_state->last_numerator_vec[i_dim] = 0;

  free(degree_table);
  free(primitive_polynomials);
  free(includ);
  gsl_rng_free(r);

  return GSL_SUCCESS;
}

static int sobol_get(void * state, unsigned int dimension, double * v)
{
  sobol_state_t * s_state = (sobol_state_t *) state;

  unsigned int i_dimension;

  /* Find the position of the least-significant zero in count. */
  int ell = 0;
  int c = s_state->sequence_count;
  while(1) {
    ++ell;
    if((c % 2) == 1) c /= 2;
    else break;
  }

  /* Check for exhaustion. */
  if(ell > SOBOL_BIT_COUNT) return GSL_EFAILED; /* FIXME: good return code here */

  for(i_dimension=0; i_dimension<dimension; i_dimension++) {
    const int direction_i     = s_state->v_direction[ell-1][i_dimension];
    const int old_numerator_i = s_state->last_numerator_vec[i_dimension];
    const int new_numerator_i = old_numerator_i ^ direction_i;
    s_state->last_numerator_vec[i_dimension] = new_numerator_i;
    v[i_dimension] = new_numerator_i * s_state->last_denominator_inv;
  }

  s_state->sequence_count++;

  return GSL_SUCCESS;
}

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
/* matrixio.c */

// 8.4.13 Example programs for matrices of gsl-ref.pdf
// GNU GSL GNU Scientific Library reference book
// pp. 99,100

/*
  Compiling advice:
  gcc matrixio.c -lgsl -lgslcblas

 */

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

int main (void) {
  int i, j, k = 0;
  gsl_matrix * m = gsl_matrix_alloc (100, 100);
  gsl_matrix * a = gsl_matrix_alloc (100, 100);

  for (i = 0; i < 100; i++)
    for (j = 0; j < 100; j++)
      gsl_matrix_set (m, i, j, 0.23 + 100*i + j);

  {
    FILE * f = fopen("test.dat", "wb");
    gsl_matrix_fwrite (f, m);
    fclose (f);
  }

  {
    FILE * f = fopen ("test.dat", "rb") ;
    gsl_matrix_fread(f, a);
    fclose (f);
  }

  for (i = 0; i < 100; i++)
    for (j = 0; j < 100; j++)
      {
	double mij = gsl_matrix_get (m, i, j);
	double aij = gsl_matrix_get (a, i, j);
	if (mij != aij) k++;
      }

  gsl_matrix_free(m);
  gsl_matrix_free(a);

  printf("differences = %d (should be zero)\n", k);

  return (k > 0);
}

max_stars_repo_name: ernestyalumni/CompPhys

max_issues_repo_name: ernestyalumni/CompPhys

max_forks_repo_name: ernestyalumni/CompPhys


text:
#include "lah.h"

#ifdef HAVE_LAPACK
#include <lapacke.h>
lah_Return lah_solveLU(lah_mat *B, lah_mat const *L, lapack_int *ipiv)
{
    if ( L == NULL ||  B == NULL || L->nR != B->nR
        || B->nR != L->nC)
    {
        return lahReturnParameterError;
    }
    /*
    if (matrix_layout == LAPACK_COL_MAJOR)
    {
        lda = L->nR;
        ldb = B->nR;
    }
    else
    {
        lda = L->nC;
        ldb = B->nC;
    }
    */

    return (!GETRS(LAH_LAPACK_LAYOUT, 'N', B->nR, B->nC, L->data , LAH_LEADING_DIM(L), ipiv, B->data, LAH_LEADING_DIM(B))) ? lahReturnOk : lahReturnExternError;
}
#endif


text:
/*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*
**                                                                                **
** This file forms part of the Underworld geophysics modelling application.         **
**                                                                                **
** For full license and copyright information, please refer to the LICENSE.md file  **
** located at the project root, or contact the authors.                             **
**                                                                                **
**~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*/
#include <petsc.h>
#include <petscvec.h>
#include <petscmat.h>
#include <petscksp.h>

#include "common-driver-utils.h"

/*

Stokes output:
  ---------------------------------
  Operator summary:
    K
    G
    f,
    h
    u
    p
  ---------------------------------
  Solution summary:
    max_u
    min_u
    average_u

    |r_1|
    |r_2|
---------------------------------
  Solver summary:
    name

---------------------------------
  Petsc build summary:

*/
PetscErrorCode BSSCR_stokes_output( PetscViewer v, Mat stokes_A, Vec stokes_b, Vec stokes_x, KSP ksp, PetscInt monitor_index )
{
	Mat K,G,D,C;
	Vec f,h, u,p;

	K = G = D = C = PETSC_NULL;
	f = h = PETSC_NULL;
	u = p = PETSC_NULL;

	MatNestGetSubMat( stokes_A, 0,0, &K );
	MatNestGetSubMat( stokes_A, 0,1, &G );
	MatNestGetSubMat( stokes_A, 1,0, &D );
	MatNestGetSubMat( stokes_A, 1,1, &C );

	VecNestGetSubVec( stokes_x, 0, &u );
	VecNestGetSubVec( stokes_x, 1, &p );

	VecNestGetSubVec( stokes_b, 0, &f );
	VecNestGetSubVec( stokes_b, 1, &h );

	PetscViewerASCIIPrintf( v, "Stokes Output:\n");
	PetscViewerASCIIPushTab( v );
	/*--------------------------------------------------------------------------------------------*/
	PetscViewerASCIIPrintf( v, "--------------------------------------------------\n");
	PetscViewerASCIIPrintf( v, "Operator summary:\n");
	PetscViewerASCIIPushTab( v );

	if (K) {	BSSCR_MatInfoLog(v,K, "stokes_A11");		PetscViewerASCIIPrintf( v, "\n");	}
	if (G) {	BSSCR_MatInfoLog(v,G, "stokes_A12");		PetscViewerASCIIPrintf( v, "\n");	}
	if (D){		BSSCR_MatInfoLog(v,D, "stokes_A21"); 		PetscViewerASCIIPrintf( v, "\n");	}
	if (C) {	BSSCR_MatInfoLog(v,C, "stokes_A22"); 		PetscViewerASCIIPrintf( v, "\n");	}

	if (f) {	BSSCR_VecInfoLog(v,f,"stokes_b1");		PetscViewerASCIIPrintf( v, "\n");	}
	if (h) {	BSSCR_VecInfoLog(v,h,"stokes_b2");		PetscViewerASCIIPrintf( v, "\n");	}

	PetscViewerASCIIPopTab( v );

	/*--------------------------------------------------------------------------------------------*/
	PetscViewerASCIIPrintf( v, "--------------------------------------------------\n");
	PetscViewerASCIIPrintf( v, "Solution summary:\n");
	PetscViewerASCIIPushTab( v );

	if (u) {	BSSCR_VecInfoLog(v,u,"x1");		PetscViewerASCIIPrintf( v, "\n");	}
	if (p) {	BSSCR_VecInfoLog(v,p,"x2");		PetscViewerASCIIPrintf( v, "\n");	}

	{
		PetscScalar s,sum;
		PetscReal r,max,min;
		PetscInt N, loc;
		double r1,r2;
		Vec K_d;
		PetscInt loc_max, loc_min;

		VecGetSize( u, &N );
		VecMax( u, &loc, &r );		PetscViewerASCIIPrintf( v, "u_max: %1.12e [%d] \n", r, loc );
		VecMin( u, &loc, &r );		PetscViewerASCIIPrintf( v, "u_min: %1.12e [%d] \n", r, loc );

		VecDot( u,u, &s );			PetscViewerASCIIPrintf( v, "u_rms: %1.12e \n", sqrt( PetscRealPart(s) )/N );

		VecDuplicate( u, &K_d );
		MatGetDiagonal( K, K_d );
		VecMax( K_d, &loc_max, &max );
		VecMin( K_d, &loc_min, &min );
		PetscViewerASCIIPrintf( v,"Algebraic contrast: max(K_d)=%.3e [%d] , min(K_d)=%.3e [%d] , max(K_d)/min(K_d) = %.8e \n", max,loc_max, min,loc_min, max/min );

		MatGetRowMax( K, K_d, PETSC_NULL );
		VecMax( K_d, &loc_max, &max );
		MatGetRowMin( K, K_d, PETSC_NULL );
		VecAbs( K_d );
		VecMin( K_d, &loc_min, &min );
		PetscViewerASCIIPrintf( v,"Algebraic contrast:   max(K)=%.3e [%d] , |min(K)|=%.3e [%d]  ,   max(K)/|min(K)| = %.8e \n", max,loc_max, min,loc_min, max/min );
		Stg_VecDestroy(&K_d );

		PetscViewerASCIIPrintf( v, "\n");

		VecGetSize( p, &N );
		VecMax( p, &loc, &r );		PetscViewerASCIIPrintf( v, "p_max:  %1.12e [%d] \n", r, loc );
		VecMin( p, &loc, &r );		PetscViewerASCIIPrintf( v, "p_min:  %1.12e [%d] \n", r, loc );

		VecDot( p,p, &s );			PetscViewerASCIIPrintf( v, "p_rms:  %1.12e \n", sqrt( PetscRealPart(s) )/N );

		VecSum( p, &sum );			PetscViewerASCIIPrintf( v, "sum(p): %1.12e \n", sum );

		PetscViewerASCIIPrintf( v, "\n");

		r1 = BSSCR_StokesMomentumResidual( K,G,f, u,p );
		PetscViewerASCIIPrintf( v, "|r1| = %1.12e <momentum> \n", r1 );
		r2 = BSSCR_StokesContinuityResidual( G,C,h, u,p );
		PetscViewerASCIIPrintf( v, "|r2| = %1.12e <continuity> \n", r2 );
		PetscViewerASCIIPrintf( v, "\n");
	}
	PetscViewerASCIIPopTab( v );

	/*--------------------------------------------------------------------------------------------*/
	if (ksp) {
		PetscViewerASCIIPrintf( v, "--------------------------------------------------\n");
		PetscViewerASCIIPrintf( v, "Solver summary:\n");
		PetscViewerASCIIPushTab( v );

		BSSCR_KSPLogSolve( v, monitor_index, ksp );
		BSSCR_BSSCR_KSPLogSolveSummary( v, monitor_index, ksp );
		PetscViewerASCIIPrintf( v, "\n");
		PetscViewerASCIIPopTab( v );
	}

	/*--------------------------------------------------------------------------------------------*/
	PetscViewerASCIIPrintf( v, "--------------------------------------------------\n");
	PetscViewerASCIIPrintf( v, "Petsc build summary:\n");
	PetscViewerASCIIPushTab( v );

	BSSCR_GeneratePetscHeader_for_viewer( v );
	PetscViewerASCIIPrintf( v, "\n");

	PetscViewerASCIIPopTab( v );
	/*--------------------------------------------------------------------------------------------*/

	PetscViewerASCIIPopTab(v);

	PetscFunctionReturn(0);
}


text:
#ifndef CWANNIER_DOS_VALUES_H
#define CWANNIER_DOS_VALUES_H

#include <stdbool.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_sort_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_eigen.h>
#include "HTightBinding.h"
#include "ctetra/dos.h"

double* DosValues(HTightBinding *Hrs, gsl_matrix *R, int na, int nb, int nc, double *Es, double num_dos, bool all_Es);

double* DosEnergyDerivValues(HTightBinding *Hrs, gsl_matrix *R, int na, int nb, int nc, double *Es, int num_dos, double num_electrons, double *fermi, double *dos_fermi, double *dos_deriv_fermi);

#endif // CWANNIER_DOS_VALUES_H

max_stars_repo_name: tflovorn/cwannier

max_issues_repo_name: tflovorn/cwannier

max_forks_repo_name: tflovorn/cwannier


text:
/* $Id$      */
/*--------------------------------------------------------------------*/
/*;  Copyright (C) 2012-2016                                          */
/*;  Associated Universities, Inc. Washington DC, USA.                */
/*;                                                                   */
/*;  This program is free software; you can redistribute it and/or    */
/*;  modify it under the terms of the GNU General Public License as   */
/*;  published by the Free Software Foundation; either version 2 of   */
/*;  the License, or (at your option) any later version.              */
/*;                                                                   */
/*;  This program is distributed in the hope that it will be useful,  */
/*;  but WITHOUT ANY WARRANTY; without even the implied warranty of   */
/*;  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the    */
/*;  GNU General Public License for more details.                     */
/*;                                                                   */
/*;  You should have received a copy of the GNU General Public        */
/*;  License along with this program; if not, write to the Free       */
/*;  Software Foundation, Inc., 675 Massachusetts Ave, Cambridge,     */
/*;  MA 02139, USA.                                                   */
/*;                                                                   */
/*;Correspondence about this software should be addressed as follows: */
/*;         Internet email: bcotton@nrao.edu.                         */
/*;         Postal address: William Cotton                            */
/*;                         National Radio Astronomy Observatory      */
/*;                         520 Edgemont Road                         */
/*;                         Charlottesville, VA 22903-2475 USA        */
/*--------------------------------------------------------------------*/
#ifndef OBITPOLNCALFIT_H
#define OBITPOLNCALFIT_H

#include "Obit.h"
#include "ObitErr.h"
#include "ObitUV.h"
#include "ObitThread.h"
#include "ObitInfoList.h"
#include "ObitSource.h"
#include "ObitSourceList.h"
#include "ObitAntennaList.h"
#include "ObitComplex.h"
#include "ObitTableCP.h"
#include "ObitTablePD.h"
#include "ObitTableBP.h"
#ifdef HAVE_GSL
#include <gsl/gsl_multifit_nlin.h>
#endif /* HAVE_GSL */

/*-------- Obit:  Merx mollis mortibus nuper ------------------*/
/**
 * \file ObitPolnCalFit.h
 *
 * ObitPolnCalFit Class for fitting polarization calibration to UV data.
 *
 * This class does least squares fitting of a full nonlinear model
 * To instrumental polarization and calibration terms and optionally source
 * polarization parameters.
 * Results are stored in a combination of AIPS PD (instrumental poln),
 * AIPS CP (source poln) and AIPS BP (gain corrections including R-L phase).
 *
 * \section ObitPolnCalFitaccess Creators and Destructors
 * An ObitPolnCalFit will usually be created using ObitPolnCalFitCreate which allows
 * specifying a name for the object as well as other information.
 *
 * A copy of a pointer to an ObitPolnCalFit should always be made using the
 * #ObitPolnCalFitRef function which updates the reference count in the object.
 * Then whenever freeing an ObitPolnCalFit or changing a pointer, the function
 * #ObitPolnCalFitUnref will decrement the reference count and destroy the object
 * when the reference count hits 0.
 * There is no explicit destructor.
 */
/*---------------Private structures----------------*/
/**
 * \enum polnParmType
 * enum for I/O file type.
 * This specifies the type of underlying data file.
 */
enum polnParmType {
  /** Unspecified */
  polnParmUnspec=0,
  /** Source parameter */
  polnParmSou,
  /** Antenna parameter */
  polnParmAnt,
  /** Antenna gain */
  polnParmGain,
  /** Phase difference  */
  polnParmPD
}; /* end enum polnParmType */

/** typedef for enum for polnParmType . */
typedef enum polnParmType PolnParmType;

/* Threaded function argument */
typedef struct {
  /* Obit error stack object */
  ObitErr      *err;
  /** thread    */
  ObitThread    *thread;
  /** thread number, >0 -> no threading   */
  olong        ithread;
  /** Do error analysis? */
  gboolean doError;
  /** First 0-rel datum in Data/Wt arrays */
  olong lo;
  /** Last 0-rel datum  in Data/Wt arrays */
  olong hi;
  /** selected antenna, 0-> all */
  olong selAnt;
  /** selected source, 0-> all */
  olong selSou;
 /** Number of visibilities being fitted */
  olong nvis;
 /** Number of data points being fitted (nvis*4) */
  olong ndata;
  /** Array of input data arrays, [10 x nvis]
      each row: 2*parallactic angle (rad), blank,
      RR_r, RR_i, LL_r, LL_i, RL_r, RL_i,  LR_r, LR_i   */
  ofloat *inData;
  /** Array of input data weights (~ 1/var), [4 x nvis]
      each row: RR, LL, RL,  LR   */
  ofloat *inWt;
  /** Source no of visibility, in range [0,nsou-1] */
  olong *souNo;
  /** Antenna numbers  of visibility, in range [0,nant-1] */
  olong *antNo;
  /** Number of antennas */
  olong nant;
  /** Reference antenna */
  olong refAnt;
  /** Fit only reference antenna R-L phase */
  gboolean doFitRL;
  /** Fit global X & Y feed gains?: */
  gboolean doFitGain;
  /** Are the feeds circularly polarized? */
  gboolean isCircFeed;
  /** Reference frequency (Hz) */
  odouble refFreq;
  /** Current frequency (Hz) */
  odouble curFreq;
  /** R-L (or X-Y) phase difference */
  odouble PD;
  /** Error estimate R-L (or X-Y) phase difference */
  odouble PDerr;
  /** Fit Stokes I poln, per cal */
  gboolean *doFitI;
  /** Fit fit source linear poln, per cal */
  gboolean *doFitPol;
  /** Fit fit Stokes V poln, per cal */
  gboolean *doFitV;
  /** Antenna parameters 4 x nant,
      each row: D1_r, D1_i, D2_r, D2_i */
  odouble *antParm;
  /** Antenna error estimates 4 x nant */
  odouble *antErr;
  /** Antenna parameters fit flags, 4 x nant */
  gboolean **antFit;
  /** Antenna parameters number, 4 x nant */
  olong **antPNumb;
  /** Antenna gains 2 x nant, each row: Gain_X, gain_Y */
  odouble *antGain;
  /** Antenna gains error estimates 2 x nant */
  odouble *antGainErr;
  /* Antenna parameters fit flags, 2 x nant */
  gboolean **antGainFit;
  /** Antenna gain parameters number, 2 x nant */
  olong **antGainPNumb;
  /** Number of calibrator sources */
  olong nsou;
  /** Source parameters 4 x nsou,
     each row: I, Poln_amp, Poln_RL (rad), V */
  odouble *souParm;
  /** Source error estimates 4 x nsou, */
  odouble *souErr;
  /** Source parameters fit flags, 4 x nsou */
  gboolean **souFit;
  /** Source parameters number, 4 x nsou */
  olong **souPNumb;
  /** Source Ids in SU table, Data */
  olong *souIDs;
  /** Inverse Source Ids lookup, index of data source ID-1
      gives calibrator number */
  olong *isouIDs;
  /** Number of parameters being fitted */
  olong nparam;
  /** Number of valid parallel observations */
  olong nPobs;
  /** Number of valid cross pol observations */
  olong nXobs;
  /** Chi squared of fit */
  ofloat ChiSq;
  /** Sum of parallel residuals */
  odouble sumParResid;
  /** Sum of cross residuals */
  odouble sumXResid;
  /** Parameter type */
  PolnParmType paramType;
  /** Parameter number */
  olong paramNumber;
  /** Sum of first derivatives */
  odouble sumDeriv;
  /** Sum of second derivatives */
  odouble sumDeriv2;
  /** Sum of weights */
  odouble sumWt;
  /** Frequency of data */
  odouble freq;
  /** R-L Phase of calibrators (rad) at Freq */
  ofloat *RLPhase;
  /** Fractional polarization of calibrators */
  ofloat *PPol;
  /** Derivative of PPol (GHz) of calibrators */
  ofloat *dPPol;
  /** Central channel (1-rel) of fit */
  olong Chan;
  /** IF number (1-rel) */
  olong IFno;
  /** complex work arrays for antenna based parameters, max 100 ant */
  dcomplex RS[100], RD[100], LS[100], LD[100],
    RSc[100], RDc[100], LSc[100], LDc[100],
    PR[100], PRc[100], PL[100], PLc[100];
  odouble SR[100], DR[100], SL[100], DL[100];
  /** Max antenna number */
  olong maxAnt;
  /** Solution method */
  gchar solnType[5];
} PolnFitArg;

/*--------------Class definitions-------------------------------------*/
/** ObitPolnCalFit Class structure. */
typedef struct {
#include "ObitPolnCalFitDef.h"   /* this class definition */
} ObitPolnCalFit;

/*----------------- Macroes ---------------------------*/
/**
 * Macro to unreference (and possibly destroy) an ObitPolnCalFit
 * returns a ObitPolnCalFit*.
 * in = object to unreference
 */
#define ObitPolnCalFitUnref(in) ObitUnref (in)

/**
 * Macro to reference (update reference count) an ObitPolnCalFit.
 * returns a ObitPolnCalFit*.
 * in = object to reference
 */
#define ObitPolnCalFitRef(in) ObitRef (in)

/**
 * Macro to determine if an object is the member of this or a
 * derived class.
 * Returns TRUE if a member, else FALSE
 * in = object to reference
 */
#define ObitPolnCalFitIsA(in) ObitIsA (in, ObitPolnCalFitGetClass())

/*---------------Public functions---------------------------*/
/** Public: Class initializer. */
void ObitPolnCalFitClassInit (void);

/** Public: Default Constructor. */
ObitPolnCalFit* newObitPolnCalFit (gchar* name);

/** Public: Create/initialize ObitPolnCalFit structures */
ObitPolnCalFit* ObitPolnCalFitCreate (gchar* name);
/** Typedef for definition of class pointer structure */
typedef ObitPolnCalFit* (*ObitPolnCalFitCreateFP) (gchar* name);

/** Public: ClassInfo pointer */
gconstpointer ObitPolnCalFitGetClass (void);

/** Public: Copy (deep) constructor. */
ObitPolnCalFit* ObitPolnCalFitCopy  (ObitPolnCalFit *in,
				     ObitPolnCalFit *out, ObitErr *err);

/** Public: Copy structure. */
void ObitPolnCalFitClone (ObitPolnCalFit *in, ObitPolnCalFit *out,
			  ObitErr *err);

/** Public: Fit a UV data */
void ObitPolnCalFitFit (ObitPolnCalFit* in, ObitUV *inUV,
			ObitUV *outUV, ObitErr *err);
/** Typedef for definition of class pointer structure */
typedef void(*ObitPolnCalFitFitFP) (ObitPolnCalFit* in, ObitUV *inUV,
				    ObitUV *outUV, ObitErr *err);

/*----------- ClassInfo Structure -----------------------------------*/
/**
 * ClassInfo Structure.
 * Contains class name, a pointer to any parent class
 * (NULL if none) and function pointers.
 */
typedef struct  {
#include "ObitPolnCalFitClassDef.h"
} ObitPolnCalFitClassInfo;

#endif /* OBITPOLNCALFIT_H */

max_stars_repo_name: sarrvesh/Obit

max_issues_repo_name: sarrvesh/Obit

max_forks_repo_name: sarrvesh/Obit


text:
/*
 * mutation.c
 *
 *  Created on: 5.5.2017
 *      Author: heine
 */

#define _XOPEN_SOURCE 600
#include <math.h>
#include <stdlib.h>
#include <stdint.h>
#include <stdio.h>
#include <time.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <mpi.h>

//static const gsl_rng_type * T; // Generator type
//static const gsl_rng * rgen;   // Generator

/*void set_mutationseed( int id ) {

  struct timespec spec;

  clock_gettime( CLOCK_REALTIME, &spec );

  gsl_rng_env_setup();
  T = gsl_rng_default;
  rgen = gsl_rng_alloc( T );

  gsl_rng_set( rgen, ( unsigned long int ) ( round( spec.tv_nsec / 1.0e6 ) + 34 * id ) );
  }*/

void
mutation(
	 long N,
	 double *x,
	 int state_dim,
	 gsl_rng * rgen )
{

  for ( long i = 0; i < N; i++ )
    for ( int d = 0; d  < state_dim; d++ )
      x[ state_dim *  i + d ] = x[ state_dim * i + d ] +
	gsl_ran_gaussian_ziggurat( rgen, 1.0 );

}

max_stars_repo_name: heinekmp/AIRPF

max_issues_repo_name: heinekmp/AIRPF

max_forks_repo_name: heinekmp/AIRPF


text:
#include <math.h>
#define FORCE_OPENBLAS_COMPLEX_STRUCT
#include <cblas.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

typedef int MKL_INT;
typedef openblas_complex_double MKL_Complex16;

#define AVX512_ALIGN 64

#include "gene_bruit_rayleigh_scalaire.c"

int main(void)
{

  int     val;
  int 	  M, N, K;
  int 	  lda, ldb, ldc;
  int     kboucle;
  float   charge;
  char    *p;

  FILE *fichier1;
  fichier1 = fopen("../results/results_openblas_gemm_fp64.dat", "w");

  int has_param_m = 0, has_param_n = 0, has_param_k = 0;
  char transa = CblasNoTrans, transb = CblasNoTrans;

  if ((p = getenv("MATMUL_M"))) {
    M = atoi(p);
    has_param_m=1;
  }
  if ((p = getenv("MATMUL_N"))) {
    N = atoi(p);
    has_param_n=1;
  }
  if ((p = getenv("MATMUL_K"))) {
    K = atoi(p);
    has_param_k=1;
  }

  if ((p = getenv("TRANSA"))) {
    if(*p == 'T')
		transa = CblasTrans;
  }
  if ((p = getenv("TRANSB"))) {
    if(*p == 'T')
		transb = CblasTrans;
  }

  /* Begin huge loop */
  for (kboucle = 1; kboucle < 31; kboucle++) {

    val = 100 * kboucle;
    M = has_param_m?M:val;
    N = has_param_n?N:val;
	K = has_param_k?K:val;

	lda = (transa == CblasNoTrans)?M:K;
    ldb = (transb == CblasNoTrans)?K:N;
    ldc = M;

    printf(">>>>> Matrix size %dx%d (K=%d, TransA=%d TransB=%d) <<<<<<\n", M, N, K,\
		transa - CblasNoTrans, transb - CblasNoTrans);

    int    iboucle, jboucle;
    int    param = 20;
    float *mat_real, *mat_imag;

    // Timers
    struct timespec tpdeb, tpfin, tpcour;
    clockid_t clock_id = CLOCK_REALTIME;
    int status2;

    double dureeloc, dureetot;
    dureetot = 0.0;

    // BLAS
    MKL_Complex16 alpha, beta;
    MKL_Complex16 *A, *B, *C;

    alpha.real = 1.0;
    alpha.imag = 0.0;
    beta.real  = 0.0;
    beta.imag  = 0.0;

    mat_real = (float *)aligned_alloc(AVX512_ALIGN, val * val * sizeof(float));
    mat_imag = (float *)aligned_alloc(AVX512_ALIGN, val * val * sizeof(float));

    A = (MKL_Complex16 *)aligned_alloc(AVX512_ALIGN, M * K * sizeof(MKL_Complex16));
    B = (MKL_Complex16 *)aligned_alloc(AVX512_ALIGN, K * N * sizeof(MKL_Complex16));
    C = (MKL_Complex16 *)aligned_alloc(AVX512_ALIGN, M * N * sizeof(MKL_Complex16));

    for (iboucle = 0; iboucle < val * val; iboucle++) {
      gene_bruit_rayleigh_scalaire(param, mat_real + iboucle,
                                   mat_imag + iboucle);
    }

    // Force hermitian matrix
    for (iboucle = 0; iboucle < K; iboucle++) {
      for (jboucle = 0; jboucle < N; jboucle++) {
        B[(iboucle * N) + jboucle].real = (double)(*(mat_real + (jboucle * K) + iboucle));

        B[(iboucle * N) + jboucle].imag =
            -1.0 * (double)(*(mat_imag + (jboucle * K) + iboucle));
      }
    }

    for (iboucle = 0; iboucle < M * K; iboucle++) {
      A[iboucle].real = (double)mat_real[iboucle];
      A[iboucle].imag = (double)mat_imag[iboucle];
    }

    status2 = clock_gettime(clock_id, &tpdeb); // start timer;

    cblas_zgemm(CblasColMajor, transa, transb, \
        M, N, K, \
        &alpha, A, lda, \
        B, ldb, &beta, C, ldc);

    status2 = clock_gettime(clock_id, &tpfin); // stop timer

    dureeloc = (float)(tpfin.tv_sec - tpdeb.tv_sec) +
               (float)(tpfin.tv_nsec - tpdeb.tv_nsec) * 1.e-9;
    dureetot = dureetot + dureeloc;

    printf("execution time = %f ms\n", dureeloc * 1000);
    charge = (float)M;
    charge = 8.0f * charge * charge * charge;
    charge = charge / dureeloc;
    printf("%f GFLOPS\n", charge / 1e9);
    fprintf(fichier1, "%5.2f\n", charge / 1e9);

    free(mat_real);
    free(mat_imag);

    free(A);
    free(B);
    free(C);
  }
  fclose(fichier1);

  return 0;
}


text:
/* specfunc/bessel_J0.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 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 */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_mode.h>
#include "bessel.h"
#include "bessel_amp_phase.h"
#include <gsl/gsl_sf_trig.h>
#include <gsl/gsl_sf_bessel.h>

#include "cheb_eval.c"

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

/* based on SLATEC besj0, 1977 version, w. fullerton */

/* chebyshev expansions for Bessel functions

 series for bj0        on the interval  0.          to  1.60000d+01
                                        with weighted error   7.47e-18
                                         log weighted error  17.13
                               significant figures required  16.98
                                    decimal places required  17.68

*/

static double bj0_data[13] = {
   0.100254161968939137,
  -0.665223007764405132,
   0.248983703498281314,
  -0.0332527231700357697,
   0.0023114179304694015,
  -0.0000991127741995080,
   0.0000028916708643998,
  -0.0000000612108586630,
   0.0000000009838650793,
  -0.0000000000124235515,
   0.0000000000001265433,
  -0.0000000000000010619,
   0.0000000000000000074,
};
static cheb_series bj0_cs = {
  bj0_data,
  12,
  -1, 1,
  9
};

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

int gsl_sf_bessel_J0_e(const double x, gsl_sf_result * result)
{
  double y = fabs(x);

  /* CHECK_POINTER(result) */

  if(y < 2.0*GSL_SQRT_DBL_EPSILON) {
    result->val = 1.0;
    result->err = y*y;
    return GSL_SUCCESS;
  }
  else if(y <= 4.0) {
    return cheb_eval_e(&bj0_cs, 0.125*y*y - 1.0, result);
  }
  else {
    const double z = 32.0/(y*y) - 1.0;
    gsl_sf_result ca;
    gsl_sf_result ct;
    gsl_sf_result cp;
    const int stat_ca = cheb_eval_e(&_gsl_sf_bessel_amp_phase_bm0_cs,  z, &ca);
    const int stat_ct = cheb_eval_e(&_gsl_sf_bessel_amp_phase_bth0_cs, z, &ct);
    const int stat_cp = gsl_sf_bessel_cos_pi4_e(y, ct.val/y, &cp);
    const double sqrty = sqrt(y);
    const double ampl  = (0.75 + ca.val) / sqrty;
    result->val  = ampl * cp.val;
    result->err  = fabs(cp.val) * ca.err/sqrty + fabs(ampl) * cp.err;
    result->err += GSL_DBL_EPSILON * fabs(result->val);
    return GSL_ERROR_SELECT_3(stat_ca, stat_ct, stat_cp);
  }
}

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

#include "eval.h"

double gsl_sf_bessel_J0(const double x)
{
  EVAL_RESULT(gsl_sf_bessel_J0_e(x, &result));
}


text:
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 * Copyright 2018 - Matteo Ragni, Matteo Cocetti - University of Trento
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

#include <cblas.h>
#include "libeuler.h"

/**
 * @brief Internal: Strut that will be passed to euler functions callbacks
 */
typedef struct euler_passthrough
{
  double ts;              /**< Integration steps. Taken from options struct */
  double alpha;           /**< Tustin coefficient. Taken from options struct */
  lapack_int u_offset;    /**< Input offset for \f$t+h\f$ callbacks. Taken from options struct */
  euler_ode_function f;   /**< Vector field to integrate. Taken from input struct */
  euler_ode_jacobian df;  /**< Jacobian of the vector field. Taken from options struct */
  const double *xk;       /**< Current state. Taken from parameters */
  lapack_int x_size;      /**< Ode dimension, taken from the input struct */
  double *work_f;         /**< Working space. Allocated in integration step */
  double *work_df;        /**< Working space. Allocated in integration step */
  void *data;             /**< User supplied data. Taken from parameters */
} euler_passtrough;

/**
 * @brief Euler implicit step wrapper
 *
 * The vector field wrapper implements the following:
 * \f{
 *   g(\cdot) = -x(t+h) + x(t) + (1-\alpha) h f(x(t), u_{1..u_{off}, p} +
 *            \alpha h f(x(t+h), u_{u_{off}..dim(u)}, p)
 * \f}
 * by using user supplied vector field callback. The implicit step search the root
 * of \f$g(\cdot)\f$.
 */
void euler_function_wrapper(double *f, const double t, const double *x, const double *u, const double **p, void *data);

/**
 * @brief Euler implicit step wrapper jacobian
 *
 * The vector field wrapper implements the following jacobian matrix:
 * \f{
 *   \nabla_{x(t+h)} g(\cdot) = -I + \alpha h \nabla_{x(t+h)} f(x(t+h), u_{u_{off}..dim(u)}, p)
 * \f}
 * by using user supplied vector field callback
 */
void euler_jacobian_wrapper(double *f, const double t, const double *x, const double *u, const double **p, void *data);

    euler_ret euler(const euler_options *opt, double *xp, const double t, const double *x, const double *u, const double **p, void *data)
{

  /* EXPLICIT IMPLEMENTATION */
/* Performin a very simple step if alpha == 0 */
if (opt->alpha == 0)
{
  opt->f(xp, t, x, u, p, opt->data);
  cblas_dscal(opt->x_size, opt->ts, xp, 1);
  cblas_daxpy(opt->x_size, 1, x, 1, xp, 1);
  return EULER_SUCCESS;
  }

  /* IMPLICIT IMPLEMENTTION */
  /* Allocating working memory */
  double *work_f = (double *)calloc(2 * opt->x_size, sizeof(double));
  if (!work_f)
    return EULER_EMALLOC;
  double *work_df = (double *)calloc(2 * opt->x_size * opt->x_size, sizeof(double));
  if(!work_df) {
    free(work_f);
    return EULER_EMALLOC;
  }
  /* Setting up the identity matrix in work_df second space, once forever */
  for (lapack_int i = 0; i < opt->x_size; i++)
    work_df[(opt->x_size * opt->x_size) + i + i * opt->x_size] = -1.0;

  /* Setting up options for Euler step */
  newton_options newton_opts = {
    opt->ordering, opt->x_size, opt->x_size,
    opt->s_tol, opt->x_tol, opt->max_iter,
    euler_function_wrapper,
    euler_jacobian_wrapper
  };

  euler_passtrough pt = {
    opt->ts, opt->alpha, opt->u_offset,
    opt->f, opt->df, x, opt->x_size,
    work_f, work_df,
    opt->data
  };

  cblas_dcopy(opt->x_size, x, 1, xp, 1);
  newton_ret nwt = newton_solve(&newton_opts, t, xp, u, p, ((void *)&pt));

  /* Freeing space */
  free(work_f);
  free(work_df);
  if (nwt > NEWTON_MAX_ITER)
    return EULER_GENERIC;
  return EULER_SUCCESS;
}

void euler_function_wrapper(double *f, const double t, const double *x, const double *u, const double **p, void *data) {
  euler_passtrough *_data = ((euler_passtrough *)data);

  /* Evaluating f(x(k), u(k)) and f(x(k+1), u(k+1)), storing result in work_f */
  _data->f(_data->work_f, t, _data->xk, u, p, _data->data);
  _data->f(_data->work_f + _data->x_size, t, x, u + _data->u_offset, p, _data->data);

  /* Computing: x(k) - x(k+1) + (1-alpha) ts f(x(k), u(k)) + alpha ts f(x(k+1), u(k+1)) */
  cblas_dscal(_data->x_size, (1 - _data->alpha) * _data->ts, _data->work_f, 1);
  cblas_dscal(_data->x_size, _data->alpha * _data->ts, _data->work_f + _data->x_size, 1);
  cblas_daxpy(_data->x_size, 1, _data->work_f + _data->x_size, 1, _data->work_f, 1);
  cblas_daxpy(_data->x_size, -1, x, 1, _data->work_f, 1);
  cblas_daxpy(_data->x_size, 1, _data->xk, 1, _data->work_f, 1);

  /* Copying result in output */
  cblas_dcopy(_data->x_size, _data->work_f, 1, f, 1);
}

void euler_jacobian_wrapper(double *df, const double t, const double *x, const double *u, const double **p, void *data) {
  euler_passtrough *_data = ((euler_passtrough *)data);

  /* Evaluating JAC(f)(x(k+1), u(k+1)) */
  _data->df(_data->work_df, t, x, u + _data->u_offset, p, _data->data);

  /* Computing: -I + alpha ts JAC(f)(x(k+1), u(k+1)) (still using level 1 Blas) */
  cblas_dscal(_data->x_size * _data->x_size, _data->alpha * _data->ts, _data->work_df, 1);
  cblas_daxpy(_data->x_size * _data->x_size, 1, _data->work_df + _data->x_size * _data->x_size, 1, _data->work_df, 1);

  /* Copying result in output */
  cblas_dcopy(_data->x_size * _data->x_size, _data->work_df, 1, df, 1);
}

max_stars_repo_name: MatteoRagni/libeuler

max_issues_repo_name: MatteoRagni/libeuler

max_forks_repo_name: MatteoRagni/libeuler


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

int main (void)
{
  double x, y;
  x = 5.0;
  y = gsl_sf_bessel_J0 (x);
  printf ("J0(%g) = %.18e\n", x, y);
  return 0;
}

max_stars_repo_name: kevintee/Predicting-Gene-Networks

max_issues_repo_name: kevintee/Predicting-Gene-Networks

max_forks_repo_name: kevintee/Predicting-Gene-Networks


text:
/* Copyright (c) 2014, Giuseppe Argentieri <giuseppe.argentieri@ts.infn.it>

 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice,
 * this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright notice,
 * this list of conditions and the following disclaimer in the documentation
 * and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY 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.
 *
 */
/*
 *
 *
 *       Filename:  polar.c
 *
 *    Description:  Cartesian <-> Polar coordinates transformation
 *
 *        Version:  1.0
 *        Created:  06/05/2014 14:21:40
 *       Revision:  none
 *        License:  BSD
 *
 *         Author:  Giuseppe Argentieri (ga), giuseppe.argentieri@ts.infn.it
 *   Organization:  Università degli Studi di Trieste
 *
 *
 */

#include <math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_math.h>
#include "funcs.h"

/*
 *      FUNCTION
 *         Name:  bloch_vector
 *  Description:  Return the Bloch vector corresponding to (r,theta,phi) polar coords
 *
 */
int bloch_vector ( gsl_vector* v, double r, double theta, double phi )
{
	/* Check the validity of coordinates */
	if ( r > 1 )                   /* Non physical state */
	{
		v = NULL ;
		return -1;
	}
	theta = fmod (theta, M_PI ) ;          /* Taking the modulo */
	phi   = fmod ( phi, 2*M_PI )  ;

	gsl_vector_set ( v, 1, r*sin(theta)*cos(phi) ) ;
	gsl_vector_set ( v, 2, r*sin(theta)*sin(phi) ) ;
	gsl_vector_set ( v, 3, r*cos(theta) ) ;

	return 0;
}		/* -----  end of function bloch_vector  ----- */

/*
 *      FUNCTION
 *         Name:  polars
 *  Description:  From the Bloch vector representation to the polar coordinates
 *
 */
int polars ( double* r, double* theta, double* phi, const gsl_vector* v )
{
	*r = gsl_hypot3( VECTOR(v,1), VECTOR(v,2), VECTOR(v,3) ) ;
	if ( *r > 1 )
		return -1 ;

	*theta = acos(VECTOR(v,3)) ;

	if ( gsl_fcmp(fabs(*theta),1,1e-9) )
		*phi = 0 ;
	else
		*phi = atan2(VECTOR(v,2),VECTOR(v,1)) ;

	return 0;
}		/* -----  end of function polars  ----- */


text:
//
// cngsl.h
//
// Copy Number Utilities that use gsl
//
// copyright 2017 Peter Andrews
//

#ifndef PAA_CNGSL_H_
#define PAA_CNGSL_H_

#include <gsl/gsl_cdf.h>

#include <algorithm>
#include <limits>
#include <string>
#include <vector>

#include "cn.h"

namespace paa {

class CNStateCall {
 public:
  CNStateCall(const std::vector<double> & state_probs_,
              const unsigned int expected_state_,
              const double segment_count_) :
      state_probs{state_probs_},
    not_state_probs(state_probs.size()),
    expected_state{expected_state_},
    segment_count{segment_count_} {
    double loss_prob{0.0};
    double not_loss_prob{0.0};
    double expected_prob{0.0};
    double not_expected_prob{0.0};
    double gain_prob{0.0};
    double not_gain_prob{0.0};
    for (unsigned int s{0}; s != state_probs.size(); ++s) {
      if (state_probs[s] > state_probs[best_call]) {
        best_call = s;
      }
      if (s < expected_state) {
        loss_prob += state_probs[s];
        not_gain_prob += state_probs[s];
        not_expected_prob += state_probs[s];
      } else if (s > expected_state) {
        gain_prob += state_probs[s];
        not_loss_prob += state_probs[s];
        not_expected_prob += state_probs[s];
      } else {
        expected_prob += state_probs[s];
        not_loss_prob += state_probs[s];
        not_gain_prob += state_probs[s];
      }
    }
    for (unsigned int n{0}; n != state_probs.size(); ++n) {
      for (unsigned int s{0}; s != state_probs.size(); ++s) {
        if (n == s) continue;
        not_state_probs[n] += state_probs[s];
      }
    }

    type = best_call < expected_state ? "loss" :
        (best_call > expected_state ? "gain" : "norm");
    const double min_prob{std::numeric_limits<double>::min()};
    score = -log10(max(best_call < expected_state ? not_loss_prob :
                       (best_call > expected_state ? not_gain_prob :
                        not_expected_prob),
                       min_prob));
    prob_not_expected = max(not_expected_prob, min_prob);
    prob_not_loss = max(not_loss_prob, min_prob);
    prob_not_gain = max(not_gain_prob, min_prob);
    }

  std::vector<double> state_probs;
  std::vector<double> not_state_probs;
  unsigned int expected_state;
  unsigned int best_call{0};
  double segment_count{0.0};
  std::string type{""};
  double score{0.0};
  double prob_not_expected{0.0};
  double prob_not_loss{0.0};
  double prob_not_gain{0.0};
};

class CNStateCaller {
 public:
  CNStateCaller(const Reference & ref_,
                const std::vector<Bin> & bins_,
                const CN_Bins & profile_) :
      ref{ref_}, bins{bins_}, profile{profile_},
    x_chr{ref.find_x_chromosome()},
    y_chr{ref.find_y_chromosome()},
    is_male{is_profile_male()},
    expected_states{get_expected_states()},
    average_autosome_count{get_average_autosome_count()},
    zero_state_count{average_autosome_count / 10} { }

  double get_average_autosome_count() const {
    double count{0};
    unsigned int n_count{0};
    for (unsigned int b{0}; b != bins.size(); ++b) {
      const unsigned int chr{bins[b].chromosome()};
      if (chr == x_chr || chr == y_chr) continue;
      count += profile[b].norm_count();
      ++n_count;
    }
    return count / n_count;
  }

  std::vector<unsigned char> get_expected_states() const {
    std::vector<unsigned char> result;
    for (unsigned int c{0}; c != ref.n_chromosomes(); ++c) {
      const bool is_x{c == x_chr};
      const bool is_y{c == y_chr};
      if (is_male && (is_x || is_y)) {
        result.push_back(1);
      } else if (!is_male && is_y) {
        result.push_back(0);
      } else {
        result.push_back(2);
      }
    }
    return result;
  }

  bool is_profile_male() const {
    double ratio{0};
    unsigned int n_ratio{0};
    for (unsigned int b{0}; b != bins.size(); ++b) {
      const unsigned int chr{bins[b].chromosome()};
      if (chr == y_chr) {
        ratio += profile[b].ratio();
        ++n_ratio;
      }
    }
    return ratio / n_ratio > 0.3333;
  }

  CNStateCall call(const unsigned int start,
                   const unsigned int stop) const {
    const double par_fraction{fraction_par(start, stop)};
    const unsigned int expected_par_copy{2};
    const unsigned int expected_state{static_cast<unsigned int>(
        0.5 + par_fraction * expected_par_copy +
        (1 - par_fraction) * expected_states[bins[start].chromosome()])};
    const double segment_count{[start, stop, this]() {
        double result{0};
        for (unsigned int b{start}; b != stop; ++b) {
          result += profile[b].norm_count();
        }
        return result;
      }()};

    std::vector<double> state_chi_squareds(n_states);
    std::vector<double> state_probs(n_states);
    double total_prob{0.0};
    double diff_power{2.0};
    while (total_prob <= 0) {
      for (unsigned int s{0}; s != n_states; ++s) {
        const double expected_count{(stop - start) *
              (s ? s * average_autosome_count / 2 : zero_state_count)};
        const double diff{expected_count > segment_count ?
              expected_count - segment_count : segment_count - expected_count};
        const double chi_squared{pow(diff, diff_power) / expected_count};
        state_chi_squareds[s] = chi_squared;
        const double prob{gsl_cdf_chisq_Q(chi_squared, 1)};
        total_prob += state_probs[s] = prob;
      }
      diff_power -= 0.1;
    }
    diff_power += 0.1;

    for (unsigned int s{0}; s != n_states; ++s) {
      state_probs[s] /= total_prob;
    }

    return CNStateCall{move(state_probs), expected_state, segment_count};
  }

  double fraction_par(const unsigned int start,
                      const unsigned int stop) const {
    if (bins[start].chromosome() != x_chr) {
      return 0.0;
    }
    const unsigned int par_coords[2][2]{{60000, 2699520},
      {154931043, 155260560}};
    if (bins[stop - 1].stop_position() <= par_coords[0][0] ||
        bins[start].start_position() >= par_coords[1][1] ||
        (bins[start].start_position() >= par_coords[0][1] &&
         bins[stop - 1].stop_position() <= par_coords[1][0])) {
      return 0.0;
    }
    double result{0.0};
    for (unsigned int b{start}; b != stop; ++b) {
      double frac_pos{0.0};
      for (const unsigned int p : {0, 1}) {
        const unsigned int par_len{par_coords[p][1] - par_coords[p][0]};
        if (bins[b].start_position() <= par_coords[p][0]) {
          if (bins[b].stop_position() >= par_coords[p][1]) {
            frac_pos += 1.0 * par_len / bins[b].length();
          } else if (bins[b].stop_position() >= par_coords[p][0]) {
            frac_pos += 1.0 * (bins[b].stop_position() - par_coords[p][0]) /
                bins[b].length();
          }
        } else if (bins[b].start_position() <= par_coords[p][1]) {
          if (bins[b].stop_position() <= par_coords[p][1]) {
            frac_pos += 1.0;
          } else {
            frac_pos += 1.0 * (par_coords[p][1] - bins[b].start_position()) /
                bins[b].length();
          }
        }
      }
      result += frac_pos;
    }
    return result / (stop - start);
  }

  const Reference & ref;
  const std::vector<Bin> & bins;
  const CN_Bins & profile;
  const unsigned int x_chr;
  const unsigned int y_chr;
  const bool is_male;
  const std::vector<unsigned char> expected_states;
  const double average_autosome_count;
  const double zero_state_count;
  const unsigned int n_states{20};
};

}  // namespace paa

#endif  // PAA_CNGSL_H_

max_stars_repo_name: docpaa/mumdex

max_stars_repo_licenses: Zlib MIT

max_issues_repo_name: docpaa/mumdex

max_issues_repo_licenses: Zlib MIT

max_forks_repo_name: docpaa/mumdex

max_forks_repo_licenses: Zlib MIT