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

#define _XOPEN_SOURCE 500

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

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

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

#include "constants.h"
#include "struct.h"
#include "splinecoeffs.h"
#include "fresnel.h"
#include "likelihood.h"

#include "wip.h"

#include <time.h> /* for testing */

/* Number of points to be used in linear integration - hardcoded for now */
#define nbptsintdefault 32768 /* Default number of points to use for linear overlaps */

/********************************* Utilities ****************************************/

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

  /* Checking the length */
  if( freq->size != noisevalues->size) {
    printf("Error: incompatible sizes in EvaluateNoise.\n");
    exit(1);
  }

  /* Checking the boundary frequencies */
  if( gsl_vector_get(freq, 0) - fLow < -1e-15 || gsl_vector_get(freq, nbpts-1) > fHigh + 1e-15 ) {
    printf("Error: incompatible frequency range in EvaluateNoise.\n");
    printf("freq[0]=%g vs fLow=%g, freq[max]=%g vs fHigh=%g\n",gsl_vector_get(freq, 0), fLow, gsl_vector_get(freq, nbpts-1), fHigh);
    printf(" %i, %i\n",gsl_vector_get(freq, 0) < fLow , gsl_vector_get(freq, nbpts-1) > fHigh);
    printf("%g\n",gsl_vector_get(freq, 0) - fLow);
    exit(1);
  }

  for(int i=0; i<nbpts; i++) {
    gsl_vector_set(noisevalues, i, ObjectFunctionCall(Snoise,gsl_vector_get(freq, i)));
  }
}

/* Function building a frequency vector with linear or logarithmic sampling */
void SetLinearFrequencies(
  gsl_vector* freqvector,    /* Output pointer to gsl_vector, already allocated */
  const double fmin,         /* Lower bound of the frequency interval */
  const double fmax,         /* Upper bound of the frequency interval */
  const int nbpts)           /* Number of points */
{
  /* Vector of frequencies with logarithmic spacing */
  double stepf = (fmax-fmin)/(nbpts-1.);
  for(int i=0; i<nbpts; i++) {
    gsl_vector_set(freqvector, i, fmin + i * stepf);
  }
}
void SetLogFrequencies(
  gsl_vector* freqvector,    /* Output pointer to gsl_vector, already allocated */
  const double fmin,         /* Lower bound of the frequency interval */
  const double fmax,         /* Upper bound of the frequency interval */
  const int nbpts)           /* Number of points */
{
  /* Vector of frequencies with logarithmic spacing */
  double lnratio = log(fmax/fmin);
  double lnfmin = log(fmin);
  for(int i=0; i<nbpts; i++) {
    gsl_vector_set(freqvector, i, exp(lnfmin + (double) i/(nbpts-1.) * lnratio));
  }
}

/* Function determining logarithmically spaced frequencies from a waveform given as a list of modes, a minimal frequency and a number of points */
void ListmodesSetFrequencies(
  struct tagListmodesCAmpPhaseFrequencySeries *list,     /* Waveform, list of modes in amplitude/phase form */
  double fLow,                                           /* Additional lower frequency limit */
  double fHigh,                                          /* Additional upper frequency limit */
  int nbpts,                                             /* Number of frequency samples */
  int tagsampling,                                       /* Tag for linear (0) or logarithmic (1) sampling */
  gsl_vector* freqvector)                                /* Output: vector of frequencies, already allocated with nbpts */
{
  /* Checking the length */
  if((int) freqvector->size != nbpts) {
    printf("Error: incompatible sizes in ListmodesSetFrequencies.\n");
    exit(1);
  }

  /* Determining the frequency interval - from the lowest frequency of the 22 mode to the highest frequency covered by at least one mode */
  double minf, maxf;
  ListmodesCAmpPhaseFrequencySeries* listelementmode22 = ListmodesCAmpPhaseFrequencySeries_GetMode(list, 2, 2);
  minf = gsl_vector_get(listelementmode22->freqseries->freq, 0);
  maxf = gsl_vector_get(list->freqseries->freq, (int) list->freqseries->freq->size - 1);
  ListmodesCAmpPhaseFrequencySeries* listelement = list;
  while(listelement) {
    maxf = fmax(maxf, gsl_vector_get(listelement->freqseries->freq, (int) listelement->freqseries->freq->size - 1));
    listelement = listelement->next;
  }

  /* Checking that the waveform covers the fLow */
  if(fLow < minf) printf("Warning: in ListmodesSetFrequencies, fLow not covered by the wave data.\n");

  /* Actual boundaries of the frequency vector - cuts what is below fLow and above fHigh */
  minf = fmax(minf, fLow);
  maxf = fmin(maxf, fHigh);

  /* Setting values of the vector of frequencies, linear or logarithmic sampling */
  if(tagsampling==0) SetLinearFrequencies(freqvector, minf, maxf, nbpts);
  else if(tagsampling==1) SetLogFrequencies(freqvector, minf, maxf, nbpts);
  else {
    printf("Error: incorrect tagsampling in ListmodesSetFrequencies.");
    exit(1);
  }
}

/* Function computing a simple trapeze integration for real data */
static double TrapezeIntegrate(const gsl_vector* x, const gsl_vector* y)
{
  if(x->size!=y->size) {
    printf("Error: trying to apply TrapezeIntegrate on vectors of different lengths.\n");
    exit(1);
  }
  int N = ((int) x->size) - 1;
  double result = 0.;
  for(int i=0; i<N; i++){
    result += (gsl_vector_get(x, i+1) - gsl_vector_get(x, i)) * (gsl_vector_get(y, i) + gsl_vector_get(y, i+1))/2.;
  }
  return result;
}

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

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

  /* Minimal, maximal frequencies */
  double fmin1 = gsl_vector_get(freqseries1->freq, 0);
  double fmin2 = gsl_vector_get(freqseries2->freq, 0);
  double fmax1 = gsl_vector_get(freqseries1->freq, size1 - 1);
  double fmax2 = gsl_vector_get(freqseries2->freq, size2 - 1);
  double fmin0 = fmax(fLow, fmax(fmin1, fmin2));
  double fmax0 = fmin(fHigh, fmin(fmax1, fmax2));

  /* Vector of frequencies with logarithmic spacing */
  gsl_vector* freqvector = gsl_vector_alloc(nbpts);
  SetLogFrequencies(freqvector, fmin0, fmax0, nbpts);

  /* Initializing the splines */
  /* Note: since this must also apply to mode contribution after processing, real and imaginary parts of the amplitude are present - but since they should differ for LLV detectors by a constant factor, they can be interpolated */
  gsl_interp_accel* accel_amp1real = gsl_interp_accel_alloc();
  gsl_interp_accel* accel_amp2real = gsl_interp_accel_alloc();
  gsl_interp_accel* accel_amp1imag = gsl_interp_accel_alloc();
  gsl_interp_accel* accel_amp2imag = gsl_interp_accel_alloc();
  gsl_interp_accel* accel_phase1 = gsl_interp_accel_alloc();
  gsl_interp_accel* accel_phase2 = gsl_interp_accel_alloc();
  gsl_spline* amp1real = gsl_spline_alloc(gsl_interp_cspline, size1);
  gsl_spline* amp2real = gsl_spline_alloc(gsl_interp_cspline, size2);
  gsl_spline* amp1imag = gsl_spline_alloc(gsl_interp_cspline, size1);
  gsl_spline* amp2imag = gsl_spline_alloc(gsl_interp_cspline, size2);
  gsl_spline* phase1 = gsl_spline_alloc(gsl_interp_cspline, size1);
  gsl_spline* phase2 = gsl_spline_alloc(gsl_interp_cspline, size2);
  gsl_vector* valuesvector = gsl_vector_alloc(nbpts);
  gsl_spline_init(amp1real, gsl_vector_const_ptr(freqseries1->freq,0), gsl_vector_const_ptr(freqseries1->amp_real,0), size1);
  gsl_spline_init(amp1imag, gsl_vector_const_ptr(freqseries1->freq,0), gsl_vector_const_ptr(freqseries1->amp_imag,0), size1);
  gsl_spline_init(amp2real, gsl_vector_const_ptr(freqseries2->freq,0), gsl_vector_const_ptr(freqseries2->amp_real,0), size2);
  gsl_spline_init(amp2imag, gsl_vector_const_ptr(freqseries2->freq,0), gsl_vector_const_ptr(freqseries2->amp_imag,0), size2);
  gsl_spline_init(phase1, gsl_vector_const_ptr(freqseries1->freq,0), gsl_vector_const_ptr(freqseries1->phase,0), size1);
  gsl_spline_init(phase2, gsl_vector_const_ptr(freqseries2->freq,0), gsl_vector_const_ptr(freqseries2->phase,0), size2);

  /* Main loop - vector of values to be evaluated */
  double f, phi1, phi2, Sn;
  double complex A1;
  double complex A2;
  double* freqvectordata = freqvector->data;
  double* valuesvectordata = valuesvector->data;
  int i=0;
  for(i=0; i<nbpts; i++){
    if(i==0) {f = fmax(freqvectordata[i], fmin0);}
    else if(i==nbpts-1) {f = fmin(freqvectordata[i], fmax0);}
    else {f = freqvectordata[i];}
    A1 = gsl_spline_eval(amp1real, f, accel_amp1real) + I*gsl_spline_eval(amp1imag, f, accel_amp1imag);
    A2 = gsl_spline_eval(amp2real, f, accel_amp2real) + I*gsl_spline_eval(amp2imag, f, accel_amp2imag);
    phi1 = gsl_spline_eval(phase1, f, accel_phase1);
    phi2 = gsl_spline_eval(phase2, f, accel_phase2);
    Sn = ObjectFunctionCall(Snoise,f);
    valuesvectordata[i] = 4.*creal(A1*conj(A2)*cexp(I*(phi1-phi2))/Sn);
  }

  /* Trapeze integration */
  res = TrapezeIntegrate(freqvector, valuesvector);

  /* Clean up */
  gsl_vector_free(freqvector);
  gsl_vector_free(valuesvector);
  gsl_interp_accel_free(accel_amp1real);
  gsl_interp_accel_free(accel_amp2real);
  gsl_interp_accel_free(accel_amp1imag);
  gsl_interp_accel_free(accel_amp2imag);
  gsl_interp_accel_free(accel_phase1);
  gsl_interp_accel_free(accel_phase2);
  gsl_spline_free(amp1real);
  gsl_spline_free(amp2real);
  gsl_spline_free(amp1imag);
  gsl_spline_free(amp2imag);
  gsl_spline_free(phase1);
  gsl_spline_free(phase2);

  return res;
}

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

  /* Determining the frequency interval - from the lowest frequency of the 22 mode to the highest frequency covered by at least one mode */
  double minf1, maxf1, minf2, maxf2, minf, maxf;
  ListmodesCAmpPhaseFrequencySeries* listelement1mode22 = ListmodesCAmpPhaseFrequencySeries_GetMode(list1, 2, 2);
  ListmodesCAmpPhaseFrequencySeries* listelement2mode22 = ListmodesCAmpPhaseFrequencySeries_GetMode(list2, 2, 2);
  minf1 = gsl_vector_get(listelement1mode22->freqseries->freq, 0);
  minf2 = gsl_vector_get(listelement2mode22->freqseries->freq, 0);
  maxf1 = ListmodesCAmpPhaseFrequencySeries_maxf(list1);
  maxf2 = ListmodesCAmpPhaseFrequencySeries_maxf(list2);
  /* Taking into account fLow, fHigh -- fstartobs, determined from 22 mode, is also used as a fLow */
  if(!(fHigh==0.)) {
    maxf1 = fmin(maxf1, fHigh);
    maxf2 = fmin(maxf2, fHigh);
  }
  if(!(fLow==0.)) {
    minf1 = fmax(minf1, fLow);
    minf2 = fmax(minf2, fLow);
  }
  if(!(fstartobs1==0.)) minf1 = fmax(minf1, fstartobs1);
  if(!(fstartobs2==0.)) minf1 = fmax(minf2, fstartobs2);
  /* Actual boundaries to be used in the overlap - intersection of freqs covered by wf 1 and 2 */
  minf = fmax(minf1, minf2);
  maxf = fmin(maxf1, maxf2);

  /* Vector of frequencies used for the overlap */
  gsl_vector* freqoverlap = gsl_vector_alloc(nbpts);
  SetLogFrequencies(freqoverlap, minf, maxf, nbpts);

  /* Evaluating each frequency series by interpolating and summing the mode contributions */
  ReImFrequencySeries* freqseries1 = NULL;
  ReImFrequencySeries_Init(&freqseries1, nbpts);
  ReImFrequencySeries_SumListmodesCAmpPhaseFrequencySeries(freqseries1, list1, freqoverlap, fLow, fHigh, fstartobs1);
  ReImFrequencySeries* freqseries2 = NULL;
  ReImFrequencySeries_Init(&freqseries2, nbpts);
  ReImFrequencySeries_SumListmodesCAmpPhaseFrequencySeries(freqseries2, list2, freqoverlap, fLow, fHigh, fstartobs2);

  /* Compute the integrand */
  gsl_vector* valuesoverlap = gsl_vector_alloc(nbpts);
  double* hreal1data = freqseries1->h_real->data;
  double* himag1data = freqseries1->h_imag->data;
  double* hreal2data = freqseries2->h_real->data;
  double* himag2data = freqseries2->h_imag->data;
  double* freqdata = freqoverlap->data;
  for(int i=0; i<nbpts; i++) {
    gsl_vector_set(valuesoverlap, i, 4.*creal( (hreal1data[i] + I*himag1data[i]) * (hreal2data[i] - I*himag2data[i]) / ObjectFunctionCall(Snoise,freqdata[i])));
  }

  /* Final trapeze integration */
  double overlap = TrapezeIntegrate(freqoverlap, valuesoverlap);

  /* Clean up */
  ReImFrequencySeries_Cleanup(freqseries1);
  ReImFrequencySeries_Cleanup(freqseries2);
  gsl_vector_free(freqoverlap);
  gsl_vector_free(valuesoverlap);

  return overlap;
}

/* Function computing the overlap (h1|h2) between h1 given in Re/Im form and h2 given as a list of mode contributions (factors sYlm already included), for a given vector of noise values - uses simple trapeze integration - generates the frequency series of h2 in Re/Im form by summing the mode contributions first, then computes the overlap */
double FDOverlapReImvsListmodesCAmpPhase(
  struct tagReImFrequencySeries *freqseries1,            /* First waveform, in Re/Im form */
  struct tagListmodesCAmpPhaseFrequencySeries *list2,    /* Second waveform, list of modes in amplitude/phase form */
  gsl_vector* noisevalues,                               /* Vector for the noise values on the freq of h1 */
  double fLow,                                           /* Minimal frequency - set to 0 to ignore */
  double fHigh,                                          /* Maximal frequency - set to 0 to ignore */
  double fstartobs2)                                     /* Starting frequency for the 22 mode of wf 2 - as determined from a limited duration of the observation - set to 0 to ignore */
{
  /* Check the lengths */
  if(freqseries1->freq->size != noisevalues->size) {
    printf("Error: inconsistent lengths in FDOverlapReImvsListmodesCAmpPhase.\n");
    exit(1);
  }

  /* Frequencies used for the overlap */
  int nbpts = (int) freqseries1->freq->size;
  gsl_vector* freqoverlap = freqseries1->freq;

  /* Evaluating frequency series 2 by interpolating and summing the mode contributions */
  ReImFrequencySeries* freqseries2 = NULL;
  ReImFrequencySeries_Init(&freqseries2, nbpts);
  ReImFrequencySeries_SumListmodesCAmpPhaseFrequencySeries(freqseries2, list2, freqoverlap, fLow, fHigh, fstartobs2);

  /* Compute the integrand */
  gsl_vector* valuesoverlap = gsl_vector_alloc(nbpts);
  double* hreal1data = freqseries1->h_real->data;
  double* himag1data = freqseries1->h_imag->data;
  double* hreal2data = freqseries2->h_real->data;
  double* himag2data = freqseries2->h_imag->data;
  double* noisedata = noisevalues->data;
  for(int i=0; i<nbpts; i++) {
    gsl_vector_set(valuesoverlap, i, 4.*creal( (hreal1data[i] + I*himag1data[i]) * (hreal2data[i] - I*himag2data[i]) / noisedata[i]));
  }

  /* Final trapeze integration */
  double overlap = TrapezeIntegrate(freqoverlap, valuesoverlap);

  /* Clean up */
  ReImFrequencySeries_Cleanup(freqseries2);
  gsl_vector_free(valuesoverlap);

  return overlap;
}

/* Function computing the overlap (h1|h2) between two waveforms given as Re/Im frequency series (common freq values assumed), for a given vector of noise values - uses simple trapeze integration */
double FDOverlapReImvsReIm(
  struct tagReImFrequencySeries *freqseries1,  /* First waveform, frequency series in Re/Im form */
  struct tagReImFrequencySeries *freqseries2,  /* Second waveform, frequency series in Re/Im form */
  gsl_vector* noisevalues)                     /* Vector for the noise values on common freq of the freqseries */
{
  /* Check the lengths */
  if(freqseries1->freq->size != noisevalues->size || freqseries2->freq->size != noisevalues->size) {
    printf("Error: inconsistent lengths in FDOverlapReImvsReIm.\n");
    exit(1);
  }

  /* Frequency vector - assuming they match beyond their mere lengths */
  gsl_vector* freqoverlap = freqseries1->freq;
  int nbpts = (int) freqoverlap->size;

  /* Compute the integrand */
  gsl_vector* valuesoverlap = gsl_vector_alloc((int) freqoverlap->size);
  double* hreal1data = freqseries1->h_real->data;
  double* himag1data = freqseries1->h_imag->data;
  double* hreal2data = freqseries2->h_real->data;
  double* himag2data = freqseries2->h_imag->data;
  double* noisedata = noisevalues->data;
  for(int i=0; i<nbpts; i++) {
    gsl_vector_set(valuesoverlap, i, 4.*creal( (hreal1data[i] + I*himag1data[i]) * (hreal2data[i] - I*himag2data[i]) / noisedata[i]));
  }

  /* Final trapeze integration */
  double overlap = TrapezeIntegrate(freqoverlap, valuesoverlap);

  /* Clean up */
  gsl_vector_free(valuesoverlap);

  return overlap;
}

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

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

  double *f1 = h1->freq->data;
  int n1 = h1->freq->size;
  double *h1Ar = h1->amp_real->data;
  double *h1Ai = h1->amp_imag->data;
  double *h1p = h1->phase->data;

  double *f2 = h2->freq->data;
  int n2 = h2->freq->size;
  double *h2Ar = h2->amp_real->data;
  double *h2Ai = h2->amp_imag->data;
  double *h2p = h2->phase->data;

  /* fLow or fHigh <= 0 means use intersection of signal domains  */
  /* NOTE: factor 4 was previously missing */
  double overlap = 4.*wip_phase(f1, n1, f2, n2, h1Ar, h1Ai, h1p, h2Ar, h2Ai, h2p, Snoise, 1.0, fLow, fHigh);
  return overlap;
}

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

  /* Main loop over the modes - goes through all the modes present */
  ListmodesCAmpPhaseFrequencySeries* listelementh1 = listh1;
  while(listelementh1) {
    ListmodesCAmpPhaseFrequencySeries* listelementh2 = listh2;
    while(listelementh2) {
      /* Scaling fstartobs1/2 with the appropriate factor of m (for the 21 mode we use m=2) - setting fmin in the overlap accordingly */
      int mmax1 = max(2, listelementh1->m);
      int mmax2 = max(2, listelementh2->m);
      double fcutLow = fmax(fLow, fmax(((double) mmax1)/2. * fstartobs1, ((double) mmax2)/2. * fstartobs2));
      overlap += FDSinglemodeWIPOverlap(listelementh1->freqseries, listelementh2->freqseries, Snoise, fcutLow, fHigh);

      listelementh2 = listelementh2->next;
    }
    listelementh1 = listelementh1->next;
  }
  return overlap;
}

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

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

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

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

/* Function computing the log likelihood -1/2(h-s|h-s), with s the signal, h the template (both given as frequency series in Re/Im form) - h, s assumed to be given on the same set of frequencies - same for the vector of noise values - fstartobs for s, h has already been taken into account */
double FDLogLikelihoodReIm(
  struct tagReImFrequencySeries *s,    /* First waveform (injection), frequency series in Re/Im form */
  struct tagReImFrequencySeries *h,    /* Second waveform (template), frequency series in Re/Im form */
  gsl_vector* noisevalues)             /* Vector for the noise values on common freq of the freqseries */
{
  /* Check the lengths */
  if(s->freq->size != noisevalues->size || h->freq->size != noisevalues->size) {
    printf("Error: inconsistent lengths in FDLogLikelihoodReIm.\n");
    exit(1);
  }

  /* Taking the difference between the frequency series: diff = h-s */
  int nbpts = (int) s->freq->size;
  ReImFrequencySeries* diff = NULL;
  ReImFrequencySeries_Init(&diff, nbpts);
  gsl_vector_memcpy(diff->freq, s->freq);
  gsl_vector_memcpy(diff->h_real, h->h_real);
  gsl_vector_memcpy(diff->h_imag, h->h_imag);
  gsl_vector_sub(diff->h_real, s->h_real);
  gsl_vector_sub(diff->h_imag, s->h_imag);

  /* Likelihood lnL = -1/2(h-s|h-s) */
  double lnL;
  lnL = -1./2 * FDOverlapReImvsReIm(diff, diff, noisevalues);

  /* Clean up */
  ReImFrequencySeries_Cleanup(diff);

  return lnL;
}

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

/*  Moved these to splinecoeffs.c

// Note: for the spines in matrix form, the first column contains the x values, so the coeffs start at 1
static double EvalCubic(
  gsl_vector* coeffs,
  double eps,
  double eps2,
  double eps3)
{
  double p0 = gsl_vector_get(coeffs, 1);
  double p1 = gsl_vector_get(coeffs, 2);
  double p2 = gsl_vector_get(coeffs, 3);
  double p3 = gsl_vector_get(coeffs, 4);
  return p0 + p1*eps + p2*eps2 + p3*eps3;
}
static double EvalQuad(
  gsl_vector* coeffs,
  double eps,
  double eps2)
{
  double p0 = gsl_vector_get(coeffs, 1);
  double p1 = gsl_vector_get(coeffs, 2);
  double p2 = gsl_vector_get(coeffs, 3);
  return p0 + p1*eps + p2*eps2;
}
*/

/* Quadratic Legendre approximation to compute values at minf and maxf when they do not fall on the grid of a freqseries, using the two first/last intervals */
static double EstimateBoundaryLegendreQuad(
  gsl_vector* vectx,   /**/
  gsl_vector* vecty,   /**/
  int j,               /**/
  double xvalue)       /**/
{
  double x0 = 0;
  double x1 = gsl_vector_get(vectx, j+1) - gsl_vector_get(vectx, j);
  double x2 = gsl_vector_get(vectx, j+2) - gsl_vector_get(vectx, j);
  double x = xvalue - gsl_vector_get(vectx, j);
  double y0 = gsl_vector_get(vecty, j);
  double y1 = gsl_vector_get(vecty, j+1);
  double y2 = gsl_vector_get(vecty, j+2);
  if(!(x>=x0 && x<=x2)) {
    printf("Error: value out of bounds in EstimateBoundaryLegendreQuad.\n");
    exit(1);
  }
  return y0*(x-x1)*(x-x2)/(x0-x1)/(x0-x2) + y1*(x-x0)*(x-x2)/(x1-x0)/(x1-x2) + y2*(x-x0)*(x-x1)/(x2-x0)/(x2-x1);
}

/* Function computing the integrand values */
void ComputeIntegrandValues(
  CAmpPhaseFrequencySeries** integrand,     /* Output: values of the integrand on common frequencies (initialized in the function) */
  CAmpPhaseFrequencySeries* freqseries1,    /* Input: frequency series for wf 1 */
  CAmpPhaseSpline* splines2,                /* Input: splines in matrix form for wf 2 */
  ObjectFunction * Snoise,           /* Noise function */
  double fLow,                              /* Lower bound of the frequency - 0 to ignore */
  double fHigh)                             /* Upper bound of the frequency - 0 to ignore */
{
  gsl_set_error_handler(&Err_Handler);

  /* Determining the boundaries of indices */
  gsl_vector* freq1 = freqseries1->freq;
  int imin1 = 0;
  int imax1 = freq1->size - 1;
  double* f1 = freq1->data;
  double f2min = gsl_matrix_get(splines2->quadspline_phase, 0, 0);
  double f2max = gsl_matrix_get(splines2->quadspline_phase, splines2->quadspline_phase->size1 - 1, 0);
  if((fLow>0 && (f1[imax1]<=fLow || f2max<=fLow)) || (fHigh>0 && (f1[imin1]>=fHigh || f2min>=fHigh))) {
    printf("Error: range of frequencies incompatible with fLow, fHigh in IntegrandValues.\n");
    printf("need one of {%g, %g} <= %g and one of {%g, %g}>=%g\n",f1[imax1],f2max,fLow,f1[imin1],f2min,fHigh);
    exit(1);
  }
  /* If starting outside, move the ends of the frequency series to be just outside the final minf and maxf */
  double minf = fmax(f1[imin1], f2min);
  double maxf = fmin(f1[imax1], f2max);
  if(fLow>0) {minf = fmax(fLow, minf);}
  if(fHigh>0) {maxf = fmin(fHigh, maxf);}
  while(f1[imin1+1]<=minf) imin1++;
  while(f1[imax1-1]>=maxf) imax1--;
  /* Estimate locally values for freqseries1 at the boundaries */
  double areal1minf = EstimateBoundaryLegendreQuad(freq1, freqseries1->amp_real, imin1, minf);
  double aimag1minf = EstimateBoundaryLegendreQuad(freq1, freqseries1->amp_imag, imin1, minf);
  double phi1minf = EstimateBoundaryLegendreQuad(freq1, freqseries1->phase, imin1, minf);
  double areal1maxf = EstimateBoundaryLegendreQuad(freq1, freqseries1->amp_real, imax1-2, maxf); /* Note the imax1-2 */
  double aimag1maxf = EstimateBoundaryLegendreQuad(freq1, freqseries1->amp_imag, imax1-2, maxf); /* Note the imax1-2 */
  double phi1maxf = EstimateBoundaryLegendreQuad(freq1, freqseries1->phase, imax1-2, maxf); /* Note the imax1-2 */

  /* Initializing output structure */
  int nbpts = imax1 + 1 - imin1;
  CAmpPhaseFrequencySeries_Init(integrand, nbpts);

  /* Loop computing integrand values */
  gsl_vector* freq = (*integrand)->freq;
  gsl_vector* ampreal = (*integrand)->amp_real;
  gsl_vector* ampimag = (*integrand)->amp_imag;
  gsl_vector* phase = (*integrand)->phase;
  double f, eps, eps2, eps3, ampreal1, ampimag1, phase1, ampreal2, ampimag2, phase2, invSn;
  double complex camp;
  double* areal1 = freqseries1->amp_real->data;
  double* aimag1 = freqseries1->amp_imag->data;
  double* phi1 = freqseries1->phase->data;
  gsl_matrix* splineAreal2 = splines2->spline_amp_real;
  gsl_matrix* splineAimag2 = splines2->spline_amp_imag;
  gsl_matrix* quadsplinephase2 = splines2->quadspline_phase;
  int i2 = 0; int j = 0;
  for(int i=imin1; i<=imax1; i++) {
    /* Distinguish the case where we are at minf or maxf */
    if(i==imin1) {
      f = minf;
      ampreal1 = areal1minf;
      ampimag1 = aimag1minf;
      phase1 = phi1minf;
    }
    else if(i==imax1) {
      f = maxf;
      ampreal1 = areal1maxf;
      ampimag1 = aimag1maxf;
      phase1 = phi1maxf;
    }
    else {
      f = gsl_vector_get(freq1, i);
      ampreal1 = areal1[i];
      ampimag1 = aimag1[i];
      phase1 = phi1[i];
    }
    /* Adjust the index in the spline if necessary and compute */
    while(gsl_matrix_get(splines2->quadspline_phase, i2+1, 0)<f) i2++;
    eps = f - gsl_matrix_get(splines2->quadspline_phase, i2, 0);
    eps2 = eps*eps;
    eps3 = eps2*eps;
    gsl_vector_view coeffsampreal2 = gsl_matrix_row(splineAreal2, i2);
    gsl_vector_view coeffsampimag2 = gsl_matrix_row(splineAimag2, i2);
    gsl_vector_view coeffsphase2 = gsl_matrix_row(quadsplinephase2, i2);
    ampreal2 = EvalCubic(&coeffsampreal2.vector, eps, eps2, eps3);
    ampimag2 = EvalCubic(&coeffsampimag2.vector, eps, eps2, eps3);
    phase2 = EvalQuad(&coeffsphase2.vector, eps, eps2);
    invSn = 1./ObjectFunctionCall(Snoise,f);
    camp = invSn * (ampreal1 + I*ampimag1) * (ampreal2 - I*ampimag2);
    gsl_vector_set(freq, j, f);
    gsl_vector_set(ampreal, j, creal(camp));
    gsl_vector_set(ampimag, j, cimag(camp));
    gsl_vector_set(phase, j, phase1 - phase2);
    j++;
  }
}

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

  /* Determining the boundaries of indices - frequency vectors assumed to be the same for channels 1,2,3 */
  gsl_vector* freq1 = freqseries1chan1->freq;
  int imin1 = 0;
  int imax1 = freq1->size - 1;
  double* f1 = freq1->data;
  double f2min = gsl_matrix_get(splines2chan1->quadspline_phase, 0, 0);
  double f2max = gsl_matrix_get(splines2chan1->quadspline_phase, splines2chan1->quadspline_phase->size1 - 1, 0);
  if((fLow>0 && (f1[imax1]<=fLow || f2max<=fLow)) || (fHigh>0 && (f1[imin1]>=fHigh || f2min>=fHigh))) {
    //printf("Error: range of frequencies incompatible with fLow, fHigh in IntegrandValues.\n");
    //printf("need both {%g, %g} > %g and both {%g, %g} < %g\n",f1[imax1],f2max,fLow,f1[imin1],f2min,fHigh);
    return -1;
  }
  /* If starting outside, move the ends of the frequency series to be just outside the final minf and maxf */
  double minf = fmax(f1[imin1], f2min);
  double maxf = fmin(f1[imax1], f2max);
  if(fLow>0) {minf = fmax(fLow, minf);}
  if(fHigh>0) {maxf = fmin(fHigh, maxf);}
  while(f1[imin1+1]<=minf) imin1++;
  while(f1[imax1-1]>=maxf) imax1--;
  //printf("imin=%i, imax=%i\n",imin1,imax1);
  int nbpts = imax1 + 1 - imin1;
  //printf("nbpts=%i\n",nbpts);
  if(nbpts<4)return -1;
  /* Estimate locally values for freqseries1 at the boundaries - phase vectors assumed to be the same for channels 1,2,3 - this is still true now that the response-processed phase includes the signal phase + R-delay phase, which is the same for all channels */
  double areal1chan1minf = EstimateBoundaryLegendreQuad(freq1, freqseries1chan1->amp_real, imin1, minf);
  double aimag1chan1minf = EstimateBoundaryLegendreQuad(freq1, freqseries1chan1->amp_imag, imin1, minf);
  double areal1chan2minf = EstimateBoundaryLegendreQuad(freq1, freqseries1chan2->amp_real, imin1, minf);
  double aimag1chan2minf = EstimateBoundaryLegendreQuad(freq1, freqseries1chan2->amp_imag, imin1, minf);
  double areal1chan3minf = EstimateBoundaryLegendreQuad(freq1, freqseries1chan3->amp_real, imin1, minf);
  double aimag1chan3minf = EstimateBoundaryLegendreQuad(freq1, freqseries1chan3->amp_imag, imin1, minf);
  double phi1minf = EstimateBoundaryLegendreQuad(freq1, freqseries1chan1->phase, imin1, minf);
  double areal1chan1maxf = EstimateBoundaryLegendreQuad(freq1, freqseries1chan1->amp_real, imax1-2, maxf); /* Note the imax1-2 */
  double aimag1chan1maxf = EstimateBoundaryLegendreQuad(freq1, freqseries1chan1->amp_imag, imax1-2, maxf); /* Note the imax1-2 */
  double areal1chan2maxf = EstimateBoundaryLegendreQuad(freq1, freqseries1chan2->amp_real, imax1-2, maxf); /* Note the imax1-2 */
  double aimag1chan2maxf = EstimateBoundaryLegendreQuad(freq1, freqseries1chan2->amp_imag, imax1-2, maxf); /* Note the imax1-2 */
  double areal1chan3maxf = EstimateBoundaryLegendreQuad(freq1, freqseries1chan3->amp_real, imax1-2, maxf); /* Note the imax1-2 */
  double aimag1chan3maxf = EstimateBoundaryLegendreQuad(freq1, freqseries1chan3->amp_imag, imax1-2, maxf); /* Note the imax1-2 */
  double phi1maxf = EstimateBoundaryLegendreQuad(freq1, freqseries1chan1->phase, imax1-2, maxf); /* Note the imax1-2 */

  /* Initializing output structure */
  CAmpPhaseFrequencySeries_Init(integrand, nbpts);

  /* Loop computing integrand values - phases are the same for chan1, chan2 and chan3 */
  gsl_vector* freq = (*integrand)->freq;
  gsl_vector* ampreal = (*integrand)->amp_real;
  gsl_vector* ampimag = (*integrand)->amp_imag;
  gsl_vector* phase = (*integrand)->phase;
  double f, eps, eps2, eps3, ampreal1chan1, ampimag1chan1, ampreal1chan2, ampimag1chan2, ampreal1chan3, ampimag1chan3, phase1, ampreal2chan1, ampimag2chan1, ampreal2chan2, ampimag2chan2, ampreal2chan3, ampimag2chan3, phase2, invSnchan1, invSnchan2, invSnchan3;
  double complex camp;
  double* areal1chan1 = freqseries1chan1->amp_real->data;
  double* aimag1chan1 = freqseries1chan1->amp_imag->data;
  double* areal1chan2 = freqseries1chan2->amp_real->data;
  double* aimag1chan2 = freqseries1chan2->amp_imag->data;
  double* areal1chan3 = freqseries1chan3->amp_real->data;
  double* aimag1chan3 = freqseries1chan3->amp_imag->data;
  double* phi1 = freqseries1chan1->phase->data;
  gsl_matrix* splinechan1real2chan1 = splines2chan1->spline_amp_real;
  gsl_matrix* splinechan1imag2chan1 = splines2chan1->spline_amp_imag;
  gsl_matrix* splinechan1real2chan2 = splines2chan2->spline_amp_real;
  gsl_matrix* splinechan1imag2chan2 = splines2chan2->spline_amp_imag;
  gsl_matrix* splinechan1real2chan3 = splines2chan3->spline_amp_real;
  gsl_matrix* splinechan1imag2chan3 = splines2chan3->spline_amp_imag;
  gsl_matrix* quadsplinephase2 = splines2chan1->quadspline_phase;
  int i2 = 0; int j = 0;
  for(int i=imin1; i<=imax1; i++) {
    /* Distinguish the case where we are at minf or maxf */
    if(i==imin1) {
      f = minf;
      ampreal1chan1 = areal1chan1minf;
      ampimag1chan1 = aimag1chan1minf;
      ampreal1chan2 = areal1chan2minf;
      ampimag1chan2 = aimag1chan2minf;
      ampreal1chan3 = areal1chan3minf;
      ampimag1chan3 = aimag1chan3minf;
      phase1 = phi1minf;
    }
    else if(i==imax1) {
      f = maxf;
      ampreal1chan1 = areal1chan1maxf;
      ampimag1chan1 = aimag1chan1maxf;
      ampreal1chan2 = areal1chan2maxf;
      ampimag1chan2 = aimag1chan2maxf;
      ampreal1chan3 = areal1chan3maxf;
      ampimag1chan3 = aimag1chan3maxf;
      phase1 = phi1maxf;
    }
    else {
      f = gsl_vector_get(freq1, i);
      ampreal1chan1 = areal1chan1[i];
      ampimag1chan1 = aimag1chan1[i];
      ampreal1chan2 = areal1chan2[i];
      ampimag1chan2 = aimag1chan2[i];
      ampreal1chan3 = areal1chan3[i];
      ampimag1chan3 = aimag1chan3[i];
      phase1 = phi1[i];
    }
    /* Adjust the index in the spline if necessary and compute */
    while(gsl_matrix_get(splines2chan1->quadspline_phase, i2+1, 0)<f) i2++;
    eps = f - gsl_matrix_get(splines2chan1->quadspline_phase, i2, 0);
    eps2 = eps*eps;
    eps3 = eps2*eps;
    gsl_vector_view coeffsampreal2chan1 = gsl_matrix_row(splinechan1real2chan1, i2);
    gsl_vector_view coeffsampimag2chan1 = gsl_matrix_row(splinechan1imag2chan1, i2);
    gsl_vector_view coeffsampreal2chan2 = gsl_matrix_row(splinechan1real2chan2, i2);
    gsl_vector_view coeffsampimag2chan2 = gsl_matrix_row(splinechan1imag2chan2, i2);
    gsl_vector_view coeffsampreal2chan3 = gsl_matrix_row(splinechan1real2chan3, i2);
    gsl_vector_view coeffsampimag2chan3 = gsl_matrix_row(splinechan1imag2chan3, i2);
    gsl_vector_view coeffsphase2 = gsl_matrix_row(quadsplinephase2, i2);
    ampreal2chan1 = EvalCubic(&coeffsampreal2chan1.vector, eps, eps2, eps3);
    ampimag2chan1 = EvalCubic(&coeffsampimag2chan1.vector, eps, eps2, eps3);
    ampreal2chan2 = EvalCubic(&coeffsampreal2chan2.vector, eps, eps2, eps3);
    ampimag2chan2 = EvalCubic(&coeffsampimag2chan2.vector, eps, eps2, eps3);
    ampreal2chan3 = EvalCubic(&coeffsampreal2chan3.vector, eps, eps2, eps3);
    ampimag2chan3 = EvalCubic(&coeffsampimag2chan3.vector, eps, eps2, eps3);
    phase2 = EvalQuad(&coeffsphase2.vector, eps, eps2);
    invSnchan1 = 1./ObjectFunctionCall(Snoise1,f);
    invSnchan2 = 1./ObjectFunctionCall(Snoise2,f);
    invSnchan3 = 1./ObjectFunctionCall(Snoise3,f);

    camp = invSnchan1 * (ampreal1chan1 + I*ampimag1chan1) * (ampreal2chan1 - I*ampimag2chan1) + invSnchan2 * (ampreal1chan2 + I*ampimag1chan2) * (ampreal2chan2 - I*ampimag2chan2) + invSnchan3 * (ampreal1chan3 + I*ampimag1chan3) * (ampreal2chan3 - I*ampimag2chan3);
    //dump
    /*
    printf("j=%i, im=%g\n a11=(%g,%g),  a12=(%g,%g),  a13=(%g,%g)\n a11=(%g,%g),  a12=(%g,%g),  a13=(%g,%g)\n",j,cimag(camp),
	   ampreal1chan1,ampimag1chan1,ampreal1chan2,ampimag1chan2,ampreal1chan3,ampimag1chan3,
	   ampreal2chan1,ampimag2chan1,ampreal2chan2,ampimag2chan2,ampreal2chan3,ampimag2chan3);
    printf("in1=%g, in2=%g, in3=%g\n",invSnchan1,invSnchan2,invSnchan3);
    */

    gsl_vector_set(freq, j, f);
    gsl_vector_set(ampreal, j, creal(camp));
    gsl_vector_set(ampimag, j, cimag(camp));
    gsl_vector_set(phase, j, phase1 - phase2);
    j++;
  }
  return 0;
}

/* Function computing the overlap (h1|h2) between two given modes in amplitude/phase form, one being already interpolated, for a given noise function - uses the amplitude/phase representation (Fresnel) */
double FDSinglemodeFresnelOverlap(
  struct tagCAmpPhaseFrequencySeries *freqseries1, /* First mode h1, in amplitude/phase form */
  struct tagCAmpPhaseSpline *splines2,             /* Second mode h2, already interpolated in matrix form */
  ObjectFunction * Snoise,                  /* Noise function */
  double fLow,                                     /* Lower bound of the frequency window for the detector */
  double fHigh)                                    /* Upper bound of the frequency window for the detector */
{
  /* Computing the integrand values, on the frequency grid of h1 */
  CAmpPhaseFrequencySeries* integrand = NULL;
  ComputeIntegrandValues(&integrand, freqseries1, splines2, Snoise, fLow, fHigh);

  /* Rescaling the integrand */
  double scaling = 10./gsl_vector_get(integrand->freq, integrand->freq->size-1);
  gsl_vector_scale(integrand->freq, scaling);
  gsl_vector_scale(integrand->amp_real, 1./scaling);
  gsl_vector_scale(integrand->amp_imag, 1./scaling);

  /* Interpolating the integrand */
  CAmpPhaseSpline* integrandspline = NULL;
  BuildSplineCoeffs(&integrandspline, integrand);

  /* Computing the integral - including here the factor 4 and the real part */
  double overlap = 4.*creal(ComputeInt(integrandspline->spline_amp_real, integrandspline->spline_amp_imag, integrandspline->quadspline_phase));

  /* Clean up */
  CAmpPhaseSpline_Cleanup(integrandspline);
  CAmpPhaseFrequencySeries_Cleanup(integrand);

  return overlap;
}

/* Function computing the overlap (h1|h2) between two given modes in amplitude/phase form for each non-correlated channel 1,2,3, one being already interpolated, for a given noise function - uses the amplitude/phase representation (Fresnel) */
double FDSinglemodeFresnelOverlap3Chan(
  struct tagCAmpPhaseFrequencySeries *freqseries1chan1, /* First mode h1 for channel 1, in amplitude/phase form */
  struct tagCAmpPhaseFrequencySeries *freqseries1chan2, /* First mode h1 for channel 2, in amplitude/phase form */
  struct tagCAmpPhaseFrequencySeries *freqseries1chan3, /* First mode h1 for channel 3, in amplitude/phase form */
  struct tagCAmpPhaseSpline *splines2chan1,             /* Second mode h2 for channel 1, already interpolated in matrix form */
  struct tagCAmpPhaseSpline *splines2chan2,             /* Second mode h2 for channel 2, already interpolated in matrix form */
  struct tagCAmpPhaseSpline *splines2chan3,             /* Second mode h2 for channel 3, already interpolated in matrix form */
  ObjectFunction * Snoisechan1,                  /* Noise function */
  ObjectFunction * Snoisechan2,                  /* Noise function */
  ObjectFunction * Snoisechan3,                  /* Noise function */
  double fLow,                                      /* Lower bound of the frequency window for the detector */
  double fHigh)                                     /* Upper bound of the frequency window for the detector */
{
  /* Computing the integrand values, on the frequency grid of h1 */
  CAmpPhaseFrequencySeries* integrand = NULL;
  if(0>ComputeIntegrandValues3Chan(&integrand, freqseries1chan1, freqseries1chan2, freqseries1chan3, splines2chan1, splines2chan2, splines2chan3, Snoisechan1, Snoisechan2, Snoisechan3, fLow, fHigh))return 0;//if allowed freq range does not exist, return 0 for overlap

  /* Rescaling the integrand */
  double scaling = 10./gsl_vector_get(integrand->freq, integrand->freq->size-1);

//TEST
//scaling = 1.;

  gsl_vector_scale(integrand->freq, scaling);
  gsl_vector_scale(integrand->amp_real, 1./scaling);
  gsl_vector_scale(integrand->amp_imag, 1./scaling);

  //dump
  /*
  for(int ii=0;ii<integrand->freq->size;ii++)
    printf("ii=%i, f=%g, integrand = ( %g, %g, %g )\n",ii,gsl_vector_get(integrand->freq,ii),gsl_vector_get(integrand->amp_real,ii),gsl_vector_get(integrand->amp_imag,ii),gsl_vector_get(integrand->phase,ii));
  */

  /* Interpolating the integrand */
  CAmpPhaseSpline* integrandspline = NULL;
  BuildSplineCoeffs(&integrandspline, integrand);

  /* Computing the integral - including here the factor 4 and the real part */
  double overlap = 4.*creal(ComputeInt(integrandspline->spline_amp_real, integrandspline->spline_amp_imag, integrandspline->quadspline_phase));

  /* Clean up */
  CAmpPhaseSpline_Cleanup(integrandspline);
  CAmpPhaseFrequencySeries_Cleanup(integrand);

  return overlap;
}

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

  /* Main loop over the modes - goes through all the modes present */
  ListmodesCAmpPhaseFrequencySeries* listelementh1 = listh1;
  while(listelementh1) {
    ListmodesCAmpPhaseSpline* listelementsplines2 = listsplines2;
    while(listelementsplines2) {
      /* Scaling fstartobs1/2 with the appropriate factor of m (for the 21 mode we use m=2) - setting fmin in the overlap accordingly */
      int mmax1 = max(2, listelementh1->m);
      int mmax2 = max(2, listelementsplines2->m);
      double fcutLow = fmax(fLow, fmax(((double) mmax1)/2. * fstartobs1, ((double) mmax2)/2. * fstartobs2));
      overlap += FDSinglemodeFresnelOverlap(listelementh1->freqseries, listelementsplines2->splines, Snoise, fcutLow, fHigh);

      listelementsplines2 = listelementsplines2->next;
    }
    listelementh1 = listelementh1->next;
  }
  return overlap;
}

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

  /* Main loop over the modes - goes through all the modes present, the same for all three channels 1,2,3 */
  ListmodesCAmpPhaseFrequencySeries* listelementh1chan1 = listh1chan1;
  while(listelementh1chan1) { /* We use the structure for channel 1 to loop through modes */
    ListmodesCAmpPhaseFrequencySeries* listelementh1chan2 = ListmodesCAmpPhaseFrequencySeries_GetMode(listh1chan2, listelementh1chan1->l, listelementh1chan1->m);
    ListmodesCAmpPhaseFrequencySeries* listelementh1chan3 = ListmodesCAmpPhaseFrequencySeries_GetMode(listh1chan3, listelementh1chan1->l, listelementh1chan1->m);
    ListmodesCAmpPhaseSpline* listelementsplines2chan1 = listsplines2chan1;
    while(listelementsplines2chan1) { /* We use the structure for channel 1 to loop through modes */
      ListmodesCAmpPhaseSpline* listelementsplines2chan2 = ListmodesCAmpPhaseSpline_GetMode(listsplines2chan2, listelementsplines2chan1->l, listelementsplines2chan1->m);
      ListmodesCAmpPhaseSpline* listelementsplines2chan3 = ListmodesCAmpPhaseSpline_GetMode(listsplines2chan3, listelementsplines2chan1->l, listelementsplines2chan1->m);
      /* Scaling fstartobs1/2 with the appropriate factor of m (for the 21 mode we use m=2) - setting fmin in the overlap accordingly */
      int mmax1 = max(2, listelementh1chan1->m);
      int mmax2 = max(2, listelementsplines2chan1->m);
      double fcutLow = fmax(fLow, fmax(((double) mmax1)/2. * fstartobs1, ((double) mmax2)/2. * fstartobs2));
      double overlapmode = FDSinglemodeFresnelOverlap3Chan(listelementh1chan1->freqseries, listelementh1chan2->freqseries, listelementh1chan3->freqseries, listelementsplines2chan1->splines, listelementsplines2chan2->splines, listelementsplines2chan3->splines, Snoise1, Snoise2, Snoise3, fcutLow, fHigh);
      overlap += overlapmode;
      listelementsplines2chan1 = listelementsplines2chan1->next;
    }
    listelementh1chan1 = listelementh1chan1->next;
  }

  return overlap;
}

max_stars_repo_name: titodalcanton/flare

max_issues_repo_name: titodalcanton/flare

max_forks_repo_name: titodalcanton/flare


text:
#ifndef neo_hookean_h
#define neo_hookean_h

#include <petsc.h>
#include "../include/structs.h"

#ifndef PHYSICS_STRUCT_NH
#define PHYSICS_STRUCT_NH
typedef struct Physics_NH_ *Physics_NH;
struct Physics_NH_ {
  CeedScalar   nu;      // Poisson's ratio
  CeedScalar   E;       // Young's Modulus
};
#endif // PHYSICS_STRUCT_NH

// Create context object
PetscErrorCode PhysicsContext_NH(MPI_Comm comm, Ceed ceed, Units *units,
                                 CeedQFunctionContext *ctx);
PetscErrorCode PhysicsSmootherContext_NH(MPI_Comm comm, Ceed ceed,
    CeedQFunctionContext ctx, CeedQFunctionContext *ctx_smoother);

// Process physics options
PetscErrorCode ProcessPhysics_NH(MPI_Comm comm, Physics_NH phys, Units units);

#endif // neo_hookean_h

ext: h

max_stars_repo_name: AdelekeBankole/libCEED

max_issues_repo_name: AdelekeBankole/libCEED

max_forks_repo_name: AdelekeBankole/libCEED


text:
#ifndef Dark_Arts_Exception_h
#define Dark_Arts_Exception_h

#include "Exception_Types.h"
#include <stdlib.h>
#include <iostream>
#include <sstream>
#include <petsc.h>

class Dark_Arts_Exception
{
public:
  Dark_Arts_Exception( EXCEPTION_TYPE ex_type , const std::string& ex_message);
  Dark_Arts_Exception( EXCEPTION_TYPE ex_type , const std::stringstream& ex_message);
  Dark_Arts_Exception();
  virtual ~Dark_Arts_Exception(){}

  void message(void) const;
  void testing_message(void) const;

private:
  std::string error_message;
};

#endif


text:
#ifndef ROUNDING_H_ZL4GHC5J
#define ROUNDING_H_ZL4GHC5J

#include <cmath>
#include <gsl/gsl>
#include <type_traits>

namespace sens_loc::math {

/// Round the value \c val to the nth \c digit.
///
/// \tparam Number arithmetic type that represents a number
/// \param val number to be rounded
/// \param digit The decimal digit the number shall be rounded to. Positive
/// values indicate the position behind the comma for rational numbers, a
/// negative number indicates a position before the comma.
/// \pre \c std::is_arithmetic_v<Number>
/// \post the result is rounded to the nth digit
/// \post \f$digit == 0 => result == val\f$
/// \post \f$!std::isfinite(val) => result == val\f$
/// \code
/// roundn(10.055, 2) == 10.06;
/// roundn(10.055, 0) == 10.055;
/// roundn(10.055, -1) == 10.00;
/// roundn(1052.4, -2) == 1100.0;
/// roundn(1052, -2) == 1100;
/// roundn(1052, 2) == 1052;
/// \endcode
template <typename Number>
Number roundn(Number val, int digit) noexcept {
    static_assert(std::is_arithmetic_v<Number>);

    if (digit == 0)
        return val;

    if constexpr (std::is_floating_point_v<Number>) {
        if (!std::isfinite(val))
            return val;
        Number potence_10 = std::pow(Number(10), digit);
        return std::round(val * potence_10) / potence_10;
    } else {
        double potence_10 = std::pow(10., digit);
        double rounded    = std::round(val * potence_10) / potence_10;
        return Number(std::round(rounded));
    }
}

}  // namespace sens_loc::math

#endif /* end of include guard: ROUNDING_H_ZL4GHC5J */

max_stars_repo_name: JonasToth/depth-conversions

max_issues_repo_name: JonasToth/depth-conversions

max_forks_repo_name: JonasToth/depth-conversions


text:
/**
 * @file aec.h
 * @brief Cancelation of a voice prompt based on either NLMS or Kalman filtering algorithms
 * @author John McDonough, Wei Chu and Kenichi Kumatani
 */

#ifndef AEC_H
#define AEC_H

#include <stdio.h>
#include <assert.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_blas.h>
#include "common/jexception.h"
#include <gsl/gsl_eigen.h>

#include "stream/stream.h"
#include "btk.h"
#include "beamformer/tracker.h"

/**
* \defgroup NLMSAcousticEchoCancellationFeature NLMS Echo Cancellation Feature
*/
/*@{*/

// ----- definition for class `NLMSAcousticEchoCancellationFeature' -----
//
class NLMSAcousticEchoCancellationFeature : public VectorComplexFeatureStream {
 public:
   NLMSAcousticEchoCancellationFeature(const VectorComplexFeatureStreamPtr& played, const VectorComplexFeatureStreamPtr& recorded, double delta = 100.0, double epsilon = 1.0E-04, double threshold = 100.0, const String& nm = "AEC");
  virtual ~NLMSAcousticEchoCancellationFeature();

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

  virtual void reset() { played_->reset(); recorded_->reset(); gsl_vector_complex_set_zero(filterCoefficient_); }

private:
  bool update_(const gsl_complex Vk);

  VectorComplexFeatureStreamPtr         played_;                    // v(n)
  VectorComplexFeatureStreamPtr         recorded_;                  // a(n)

  unsigned                              fftLen_;
  unsigned                              fftLen2_;
  gsl_vector_complex*                   filterCoefficient_;

  const double                          delta_;
  const double                          epsilon_;
  const double                          threshold_;
};

typedef Inherit<NLMSAcousticEchoCancellationFeature, VectorComplexFeatureStreamPtr> NLMSAcousticEchoCancellationFeaturePtr;
/*@}*/

/**
* \defgroup KalmanFilterEchoCancellationFeature Kalman Filer Echo Cancellation Feature
*/
/*@{*/

// ----- definition for class `KalmanFilterEchoCancellationFeature' -----
//
class KalmanFilterEchoCancellationFeature : public VectorComplexFeatureStream {
 public:
  KalmanFilterEchoCancellationFeature(const VectorComplexFeatureStreamPtr& played, const VectorComplexFeatureStreamPtr& recorded, double beta = 0.95, double sigma2 = 100.0, double threshold = 100.0, const String& nm = "KFEchoCanceller");
  virtual ~KalmanFilterEchoCancellationFeature();

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

  virtual void reset() { played_->reset(); recorded_->reset(); gsl_vector_complex_set_zero(filterCoefficient_); }

private:
  bool update_(const gsl_complex Vk);

  VectorComplexFeatureStreamPtr         played_;                    // v(n)
  VectorComplexFeatureStreamPtr         recorded_;                  // a(n)

  unsigned                              fftLen_;
  unsigned                              fftLen2_;

  gsl_vector_complex*                   filterCoefficient_;
  gsl_vector*                           sigma2_v_;
  gsl_vector*                           K_k_;

  const double                          beta_;
  const double                          threshold_;
  const double                          sigma2_u_;
};

typedef Inherit<KalmanFilterEchoCancellationFeature, VectorComplexFeatureStreamPtr> KalmanFilterEchoCancellationFeaturePtr;
/*@}*/

// ----- definition for class `BlockKalmanFilterEchoCancellationFeature' -----
//
class BlockKalmanFilterEchoCancellationFeature : public VectorComplexFeatureStream {
  public:
    BlockKalmanFilterEchoCancellationFeature(const VectorComplexFeatureStreamPtr& played, const VectorComplexFeatureStreamPtr& recorded, unsigned sampleN = 1, double beta = 0.95, double sigmau2 = 10e-4, double sigmauk2 = 5.0, double threshold = 100.0, double amp4play = 1.0, const String& nm = "KFEchoCanceller");
    virtual ~BlockKalmanFilterEchoCancellationFeature();

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

    virtual void reset()
    {
      played_->reset(); recorded_->reset();
    }

  protected:
    class ComplexBuffer_ {
      public:
    /*
        @brief Construct a circular buffer to hold past and current subband samples
        It keeps nsamp arrays which is completely updated with the period 'nsamp'.
        Each array holds actual values of the samples.
        @param unsigned len[in] is the size of each vector of samples
        @param unsigned nsamp[in] is the period of the circular buffer
    */
      ComplexBuffer_(unsigned len, unsigned sampleN)
        : len_(len), sampleN_(sampleN), zero_(sampleN_ - 1),
          samples_(new gsl_vector_complex*[sampleN_]), subbandSamples_(gsl_vector_complex_calloc(sampleN_))
      {
        for (unsigned i = 0; i < sampleN_; i++)
          samples_[i] = gsl_vector_complex_calloc(len_);
      }

      ~ComplexBuffer_()
      {
        for (unsigned i = 0; i < sampleN_; i++)
          gsl_vector_complex_free(samples_[i]);
        delete[] samples_;
        gsl_vector_complex_free(subbandSamples_);
      }

      gsl_complex sample(unsigned timeX, unsigned binX) const {
        unsigned idx = index_(timeX);
        const gsl_vector_complex* vec = samples_[idx];
        return gsl_vector_complex_get(vec, binX);
      }

      const gsl_vector_complex* get_samples(unsigned m)
      {
        for (unsigned timeX = 0; timeX < sampleN_; timeX++)
          gsl_vector_complex_set(subbandSamples_, timeX, sample(timeX, m));

        return subbandSamples_;
      }

      void next_sample(const gsl_vector_complex* s = NULL, double amp4play = 1.0 ) {
        zero_ = (zero_ + 1) % sampleN_;

        gsl_vector_complex* nextBlock = samples_[zero_];

        if (s == NULL) {
          gsl_vector_complex_set_zero(nextBlock);
        } else {
          if (s->size != len_)
            throw jdimension_error("'ComplexBuffer_': Sizes do not match (%d vs. %d)", s->size, len_);
          assert( s->size == len_ );
          gsl_vector_complex_memcpy(nextBlock, s);
          if( amp4play != 1.0 )
            gsl_blas_zdscal( amp4play, nextBlock );
        }
      }

      void zero() {
        for (unsigned i = 0; i < sampleN_; i++)
          gsl_vector_complex_set_zero(samples_[i]);
        zero_ = sampleN_ - 1;
      }

    private:
      unsigned index_(unsigned idx) const {
        assert(idx < sampleN_);
        unsigned ret = (zero_ + sampleN_ - idx) % sampleN_;
        return ret;
      }

      const unsigned                             len_;
      const unsigned                             sampleN_;
      unsigned                                   zero_; // index of most recent sample
      gsl_vector_complex**                       samples_;
      gsl_vector_complex*                        subbandSamples_;
    };

  static gsl_complex                             ComplexOne_;
  static gsl_complex                             ComplexZero_;

  bool update_(const gsl_vector_complex* Vk);
  void conjugate_(gsl_vector_complex* dest, const gsl_vector_complex* src) const;

  VectorComplexFeatureStreamPtr                 played_;   // v(n)
  VectorComplexFeatureStreamPtr                 recorded_; // a(n)

  unsigned                                      fftLen_;
  unsigned                                      fftLen2_;
  unsigned                                      sampleN_;

  ComplexBuffer_                                buffer_;

  gsl_vector_complex**                          filterCoefficient_;
  gsl_vector*                                   sigma2_v_;
  gsl_matrix_complex**                          K_k_;
  gsl_matrix_complex*                           K_k_k1_;

  const double                                  beta_;
  const double                                  threshold_;
  gsl_matrix_complex**                          Sigma2_u_;
  gsl_vector_complex*                           Gk_;
  gsl_vector_complex*                           scratch_;
  gsl_vector_complex*                           scratch2_;
  gsl_matrix_complex*                           scratchMatrix_;
  gsl_matrix_complex*                           scratchMatrix2_;
  double                                        amp4play_;
  double                                        floorVal_;
  int                                           skippedN_;
  int                                           maxSkippedN_;
};

typedef Inherit<BlockKalmanFilterEchoCancellationFeature, VectorComplexFeatureStreamPtr> BlockKalmanFilterEchoCancellationFeaturePtr;
/*@}*/

/**
* \defgroup DTDBlockKalmanFilterEchoCancellationFeature Block Kalman Filer Echo Cancellation Feature
*/
/*@{*/

// ----- definition for class `InformationFilterEchoCancellationFeature' -----
//
class InformationFilterEchoCancellationFeature : public BlockKalmanFilterEchoCancellationFeature {
 public:
  InformationFilterEchoCancellationFeature(const VectorComplexFeatureStreamPtr& played, const VectorComplexFeatureStreamPtr& recorded, unsigned sampleN = 1, double beta = 0.95, double sigmau2 = 10e-4, double sigmauk2 = 5.0, double snrTh = 2.0, double engTh = 100.0, double smooth = 0.9, double loading = 1.0e-02,  double amp4play = 1.0, const String& nm = "Information Echo Canceller");

  virtual ~InformationFilterEchoCancellationFeature();

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

protected:
  static double _EigenValueThreshold;

  double update_band_(const gsl_complex Ak, const gsl_complex Ek, int frame_no, unsigned m);
  void invert_(gsl_matrix_complex* matrix);
  static void print_matrix_(const gsl_matrix_complex* mat);
  static void print_vector_(const gsl_vector_complex* vec);

  const double                smoothEk_; // for smoothing the error signal, ->1, less smooth
  const double                smoothSk_; // for smoothing the estimated signal
  const double                engTh_;    // threshold of energy

  double*                     snr_;
  double*                     EkEnergy_;
  double*                     SkEnergy_;
  const double                loading_;

  gsl_matrix_complex*         inverse_;
  gsl_eigen_hermv_workspace*  eigenWorkSpace_;
  gsl_vector*                 evalues_;

  gsl_vector_complex*         scratchInverse_;
  gsl_matrix_complex*         scratchMatrixInverse_;
  gsl_matrix_complex*         scratchMatrixInverse2_;
  gsl_matrix_complex*         matrixCopy_;
};

typedef Inherit<InformationFilterEchoCancellationFeature, BlockKalmanFilterEchoCancellationFeaturePtr> InformationFilterEchoCancellationFeaturePtr;

// ----- definition for class `SquareRootInformationFilterEchoCancellationFeature' -----
//
class SquareRootInformationFilterEchoCancellationFeature : public InformationFilterEchoCancellationFeature {
 public:
  SquareRootInformationFilterEchoCancellationFeature(const VectorComplexFeatureStreamPtr& played, const VectorComplexFeatureStreamPtr& recorded, unsigned sampleN = 1, double beta = 0.95, double sigmau2 = 10e-4, double sigmauk2 = 5.0, double snrTh = 2.0, double engTh = 100.0, double smooth = 0.9, double loading = 1.0e-02, double amp4play = 1.0, const String& nm = "Square Root Information Echo Canceller");

  virtual ~SquareRootInformationFilterEchoCancellationFeature();

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

private:
  const gsl_complex load_;
  static gsl_complex calc_givens_rotation_(const gsl_complex& v1, const gsl_complex& v2, gsl_complex& c, gsl_complex& s);
  static void apply_givens_rotation_(const gsl_complex& v1, const gsl_complex& v2, const gsl_complex& c, const gsl_complex& s,
                                   gsl_complex& v1p, gsl_complex& v2p);
  void negative_(gsl_matrix_complex* dest, const gsl_matrix_complex* src);
  void extract_covariance_state_(const gsl_matrix_complex* K_k, const gsl_vector_complex* sk, gsl_vector_complex* xk);

  void temporal_update_(unsigned m);
  void observational_update_(unsigned m, const gsl_complex& Ak, double sigma2_v);
  void diagonal_loading_(unsigned m);

  gsl_vector_complex**                                informationState_;
};

typedef Inherit<SquareRootInformationFilterEchoCancellationFeature, InformationFilterEchoCancellationFeaturePtr> SquareRootInformationFilterEchoCancellationFeaturePtr;

// ----- definition for class `DTDBlockKalmanFilterEchoCancellationFeature' -----
//
class DTDBlockKalmanFilterEchoCancellationFeature : public BlockKalmanFilterEchoCancellationFeature {
 public:
   DTDBlockKalmanFilterEchoCancellationFeature(const VectorComplexFeatureStreamPtr& played, const VectorComplexFeatureStreamPtr& recorded, unsigned sampleN = 1, double beta = 0.95, double sigmau2 = 10e-4, double sigmauk2 = 5.0, double snrTh = 2.0, double engTh = 100.0, double smooth = 0.9, double amp4play = 1.0, const String& nm = "DTDKFEchoCanceller");
  virtual ~DTDBlockKalmanFilterEchoCancellationFeature();

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

private:
  double update_band_(const gsl_complex Ak, const gsl_complex Ek, int frame_no);
  void conjugate_(gsl_vector_complex* dest, const gsl_vector_complex* src) const;

  const double        smoothSk_; // for smoothing the estimated signal
  const double        smoothEk_; // for smoothing the error signal, ->1, less smooth
  const double        engTh_;    // threshold of energy
  double              snr_;

  double              EkEnergy_;
  double              SkEnergy_;
  FILE*               fdb_;
};

typedef Inherit<DTDBlockKalmanFilterEchoCancellationFeature, BlockKalmanFilterEchoCancellationFeaturePtr> DTDBlockKalmanFilterEchoCancellationFeaturePtr;
/*@}*/

#endif // AEC_H


text:
#pragma once

#include <imageview/IsPixelFormat.h>
#include <imageview/internal/ImageViewStorage.h>

#include <gsl/span>

namespace imageview {
namespace detail {

template <class PixelFormat>
class PixelRef {
 public:
  static_assert(IsPixelFormat<PixelFormat>::value, "Not a PixelFormat.");

  using color_type = typename PixelFormat::color_type;

  constexpr PixelRef(gsl::span<std::byte, PixelFormat::kBytesPerPixel> pixel_data,
                     std::reference_wrapper<const PixelFormat> pixel_format);

  constexpr PixelRef(const PixelRef& other) = default;
  constexpr PixelRef(PixelRef&&) = default;
  ~PixelRef() = default;

  // Implicit conversion to color_type.
  // \return the color of the referenced pixel.
  constexpr operator color_type() const;

  // Assigns the specified color to the referenced pixel.
  // \param color - color to assign.
  // \return *this.
  constexpr PixelRef& operator=(const color_type& color);

  // Assigns the specified color to the referenced pixel.
  // Note: PixelRef has reference semantics: copy and move assignment operators change the value of the
  // referenced pixel rather than PixelRef object itself.
  constexpr PixelRef& operator=(const PixelRef& other);
  constexpr PixelRef& operator=(PixelRef&& other);

 private:
  // TODO: store PixelFormat by reference if it's not eligible for empty base optimization,
  // otherwise by value (as an empty base).
  detail::ImageViewStorage<PixelFormat, true> storage_;
};

template <class PixelFormat>
constexpr PixelRef<PixelFormat>::PixelRef(gsl::span<std::byte, PixelFormat::kBytesPerPixel> pixel_data,
                                          std::reference_wrapper<const PixelFormat> pixel_format)
    : storage_(pixel_data.data(), pixel_format.get()) {}

template <class PixelFormat>
constexpr PixelRef<PixelFormat>::operator color_type() const {
  constexpr std::size_t kBytesPerPixel = PixelFormat::kBytesPerPixel;
  const gsl::span<const std::byte, kBytesPerPixel> pixel_data(storage_.data_, kBytesPerPixel);
  return storage_.pixelFormat().read(pixel_data);
}

template <class PixelFormat>
constexpr PixelRef<PixelFormat>& PixelRef<PixelFormat>::operator=(const color_type& color) {
  constexpr std::size_t kBytesPerPixel = PixelFormat::kBytesPerPixel;
  const gsl::span<std::byte, kBytesPerPixel> pixel_data(storage_.data_, kBytesPerPixel);
  storage_.pixelFormat().write(color, pixel_data);
  return *this;
}

template <class PixelFormat>
constexpr PixelRef<PixelFormat>& PixelRef<PixelFormat>::operator=(const PixelRef& other) {
  // TODO: In theory, we could just memcpy() the binary data. However, if PixelFormat is stateful,
  // and the state of our PixelFormat differs from the state of @other, then simply copying the binary
  // data might lead to the wrong result.
  // Possible workaround: add a constexpr if for the case when PixelFormat is an empty class.
  return *this = static_cast<color_type>(other);
}

template <class PixelFormat>
constexpr PixelRef<PixelFormat>& PixelRef<PixelFormat>::operator=(PixelRef&& other) {
  return *this = static_cast<color_type>(other);
}

}  // namespace detail
}  // namespace imageview

max_stars_repo_name: alexanderbelous/imageview

max_issues_repo_name: alexanderbelous/imageview

max_forks_repo_name: alexanderbelous/imageview


text:
#pragma once

#if WIN32
    #define LAPACK_COMPLEX_CUSTOM
    #define lapack_complex_float float
    #define lapack_complex_double double
#endif

#include <cblas.h>
#include <lapacke.h>
#undef I // undefine this def due to complex.h that causes issues later

#define BLAS(x,X) cblas_ ## x
#define LAPACK(x,X) LAPACKE_ ## x

#define CLAPACK_INTERFACE
#define UNUSED(expr) do{(void)(expr);}while(0)

#define DLANGE(X1,X2,X3,X4,X5,X6)       LAPACK(dlange,DLANGE)(LAPACK_COL_MAJOR,X1,X2,X3,X4,X5);UNUSED(X6)

#define DSPTRF(X1,X2,X3,X4,X5)          LAPACK(dsptrf,DSPTRF)(LAPACK_COL_MAJOR,X1,X2,X3,X4)
#define DSPTRS(X1,X2,X3,X4,X5,X6,X7,X8) LAPACK(dsptrs,DSPTRS)(LAPACK_COL_MAJOR,X1,X2,X3,X4,X5,X6,X7)
#define DSPTRI(X1,X2,X3,X4,X5,X6)       LAPACK(dsptri,DSPTRI)(LAPACK_COL_MAJOR,X1,X2,X3,X4)
#define DPPTRF(X1,X2,X3,X4)             LAPACK(dpptrf,DPPTRF)(LAPACK_COL_MAJOR,X1,X2,X3)
#define DPPTRI(X1,X2,X3,X4)             LAPACK(dpptri,DPPTRI)(LAPACK_COL_MAJOR,X1,X2,X3)
#define DGETRF(X1,X2,X3,X4,X5)          LAPACK(dgetrf,DGETRF)(LAPACK_COL_MAJOR,X1,X2,X3,X4,X5)
#define DGETRI(X1,X2,X3,X4)             LAPACK(dgetri,DGETRI)(LAPACK_COL_MAJOR,X1,X2,X3,X4)

#define DGESDD(X1,X2,X3,X4,X5,X6,X7,X8,X9,X10,X11,X12,X13,X14) LAPACK(dgesdd,DGESDD)(LAPACK_COL_MAJOR,X1,X2,X3,X4,X5,X6,X7,X8,X9,X10)

#include <BlasLapackImplementations/OpenMEEGMathsCBlasLapack.h>

max_stars_repo_name: Tobitta/openmeeg

max_stars_repo_licenses: Intel CECILL-B

max_issues_repo_name: Tobitta/openmeeg

max_issues_repo_licenses: Intel CECILL-B

max_forks_repo_name: Tobitta/openmeeg

max_forks_repo_licenses: Intel CECILL-B


text:
/**
 *
 * @file example_sgelqf.c
 *
 *  PLASMA testing routines
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @brief Example using LQ factorization
 *
 * @version 2.6.0
 * @author Bilel Hadri
 * @date 2010-11-15
 * @generated s Tue Jan  7 11:45:20 2014
 *
 **/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

#include <plasma.h>
#include <cblas.h>
#include <lapacke.h>
#include <core_blas.h>

#ifndef max
#define max(a, b) ((a) > (b) ? (a) : (b))
#endif
#ifndef min
#define min(a, b) ((a) < (b) ? (a) : (b))
#endif

int check_orthogonality(int, int, int, float*);
int check_factorization(int, int, float*, float*, int, float*);

int IONE=1;
int ISEED[4] = {0,0,0,1};   /* initial seed for slarnv() */

int main ()
{

    int cores = 2;
    int M     = 10;
    int N     = 15;
    int LDA   = 10;
    int K = min(M, N);
    int info;
    int info_ortho, info_factorization;
    int i,j;
    int LDAxN = LDA*N;

    float *A1 = (float *)malloc(LDA*N*sizeof(float));
    float *A2 = (float *)malloc(LDA*N*sizeof(float));
    float *Q  = (float *)malloc(LDA*N*sizeof(float));
    PLASMA_desc *T;

    /* Check if unable to allocate memory */
    if ((!A1)||(!A2)||(!Q)){
        printf("Out of Memory \n ");
        return EXIT_SUCCESS;
    }

    /* Plasma Initialization */
    PLASMA_Init(cores);
    printf("-- PLASMA is initialized to run on %d cores. \n",cores);

    /* Allocate T */
    PLASMA_Alloc_Workspace_sgelqf(M, N, &T);

    /* Initialize A1 and A2 */
    LAPACKE_slarnv_work(IONE, ISEED, LDAxN, A1);
    for (i = 0; i < M; i++)
        for (j = 0; j < N; j++)
            A2[LDA*j+i] = A1[LDA*j+i] ;

    /* Factorization QR of the matrix A2 */
    info = PLASMA_sgelqf(M, N, A2, LDA, T);

    /* Building the economy-size Q */
    memset((void*)Q, 0, LDA*N*sizeof(float));
    for (i = 0; i < K; i++)
        Q[LDA*i+i] = 1.0;

    PLASMA_sorglq(M, N, K, A2, LDA, T, Q, LDA);

    /* Check the orthogonality, factorization and the solution */
    info_ortho = check_orthogonality(M, N, LDA, Q);
    info_factorization = check_factorization(M, N, A1, A2, LDA, Q);

    if ((info_ortho != 0)|(info_factorization != 0)|(info != 0))
       printf("-- Error in SGELQF example ! \n");
    else
       printf("-- Run of SGELQF example successful ! \n");

    free(A1); free(A2); free(Q); free(T);

    PLASMA_Finalize();

    return EXIT_SUCCESS;
}

/*-------------------------------------------------------------------
 * Check the orthogonality of Q
 */

int check_orthogonality(int M, int N, int LDQ, float *Q)
{
    float alpha, beta;
    float normQ;
    int info_ortho;
    int i;
    int minMN = min(M, N);
    float eps;
    float *work = (float *)malloc(minMN*sizeof(float));

    eps = LAPACKE_slamch_work('e');
    alpha = 1.0;
    beta  = -1.0;

    /* Build the idendity matrix USE DLASET?*/
    float *Id = (float *) malloc(minMN*minMN*sizeof(float));
    memset((void*)Id, 0, minMN*minMN*sizeof(float));
    for (i = 0; i < minMN; i++)
        Id[i*minMN+i] = (float)1.0;

    /* Perform Id - Q'Q */
    if (M >= N)
        cblas_ssyrk(CblasColMajor, CblasUpper, CblasTrans, N, M, alpha, Q, LDQ, beta, Id, N);
    else
        cblas_ssyrk(CblasColMajor, CblasUpper, CblasNoTrans, M, N, alpha, Q, LDQ, beta, Id, M);

    normQ = LAPACKE_slansy_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), 'u', minMN, Id, minMN, work);

    printf("============\n");
    printf("Checking the orthogonality of Q \n");
    printf("||Id-Q'*Q||_oo / (N*eps) = %e \n",normQ/(minMN*eps));

    if ( isnan(normQ / (minMN * eps)) || (normQ / (minMN * eps) > 10.0) ) {
        printf("-- Orthogonality is suspicious ! \n");
        info_ortho=1;
    }
    else {
        printf("-- Orthogonality is CORRECT ! \n");
        info_ortho=0;
    }

    free(work); free(Id);

    return info_ortho;
}

/*------------------------------------------------------------
 *  Check the factorization QR
 */

int check_factorization(int M, int N, float *A1, float *A2, int LDA, float *Q)
{
    float Anorm, Rnorm;
    float alpha, beta;
    int info_factorization;
    int i,j;
    float eps;

    eps = LAPACKE_slamch_work('e');

    float *Ql       = (float *)malloc(M*N*sizeof(float));
    float *Residual = (float *)malloc(M*N*sizeof(float));
    float *work              = (float *)malloc(max(M,N)*sizeof(float));

    alpha=1.0;
    beta=0.0;

    if (M >= N) {
        /* Extract the R */
        float *R = (float *)malloc(N*N*sizeof(float));
        memset((void*)R, 0, N*N*sizeof(float));
        LAPACKE_slacpy_work(LAPACK_COL_MAJOR,'u', M, N, A2, LDA, R, N);

        /* Perform Ql=Q*R */
        memset((void*)Ql, 0, M*N*sizeof(float));
        cblas_sgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, M, N, N, (alpha), Q, LDA, R, N, (beta), Ql, M);
        free(R);
    }
    else {
        /* Extract the L */
        float *L = (float *)malloc(M*M*sizeof(float));
        memset((void*)L, 0, M*M*sizeof(float));
        LAPACKE_slacpy_work(LAPACK_COL_MAJOR,'l', M, N, A2, LDA, L, M);

    /* Perform Ql=LQ */
        memset((void*)Ql, 0, M*N*sizeof(float));
        cblas_sgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, M, N, M, (alpha), L, M, Q, LDA, (beta), Ql, M);
        free(L);
    }

    /* Compute the Residual */
    for (i = 0; i < M; i++)
        for (j = 0 ; j < N; j++)
            Residual[j*M+i] = A1[j*LDA+i]-Ql[j*M+i];

    Rnorm = LAPACKE_slange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), M, N, Residual, M, work);
    Anorm = LAPACKE_slange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), M, N, A2, LDA, work);

    if (M >= N) {
        printf("============\n");
        printf("Checking the QR Factorization \n");
        printf("-- ||A-QR||_oo/(||A||_oo.N.eps) = %e \n",Rnorm/(Anorm*N*eps));
    }
    else {
        printf("============\n");
        printf("Checking the LQ Factorization \n");
        printf("-- ||A-LQ||_oo/(||A||_oo.N.eps) = %e \n",Rnorm/(Anorm*N*eps));
    }

    if (isnan(Rnorm / (Anorm * N *eps)) || (Rnorm / (Anorm * N * eps) > 10.0) ) {
        printf("-- Factorization is suspicious ! \n");
        info_factorization = 1;
    }
    else {
        printf("-- Factorization is CORRECT ! \n");
        info_factorization = 0;
    }

    free(work); free(Ql); free(Residual);

    return info_factorization;
}

ext: c


text:
/*****************************************************************************
*
* Rokko: Integrated Interface for libraries of eigenvalue decomposition
*
* Copyright (C) 2012-2019 by Synge Todo <wistaria@phys.s.u-tokyo.ac.jp>
*
* Distributed under the Boost Software License, Version 1.0. (See accompanying
* file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
*
*****************************************************************************/

#include <stdio.h>
#include <cblas.h>
#include <rokko/cmatrix.h>

int main() {
  int n = 4;
  double** a = alloc_dmatrix(n, n);
  mat_elem(a, 0, 0) = 0.959291425205444;
  mat_elem(a, 0, 1) = 0.257508254123736;
  mat_elem(a, 0, 2) = 0.243524968724989;
  mat_elem(a, 0, 3) = 0.251083857976031;
  mat_elem(a, 1, 0) = 0.547215529963803;
  mat_elem(a, 1, 1) = 0.840717255983663;
  mat_elem(a, 1, 2) = 0.929263623187228;
  mat_elem(a, 1, 3) = 0.616044676146639;
  mat_elem(a, 2, 0) = 0.138624442828679;
  mat_elem(a, 2, 1) = 0.254282178971531;
  mat_elem(a, 2, 2) = 0.349983765984809;
  mat_elem(a, 2, 3) = 0.473288848902729;
  mat_elem(a, 3, 0) = 0.149294005559057;
  mat_elem(a, 3, 1) = 0.814284826068816;
  mat_elem(a, 3, 2) = 0.196595250431208;
  mat_elem(a, 3, 3) = 0.351659507062997;
  double** b = alloc_dmatrix(n, n);
  mat_elem(b, 0, 0) = 0.567821640725221;
  mat_elem(b, 0, 1) = 0.779167230102011;
  mat_elem(b, 0, 2) = 0.469390641058206;
  mat_elem(b, 0, 3) = 0.794284540683907;
  mat_elem(b, 1, 0) = 0.0758542895630636;
  mat_elem(b, 1, 1) = 0.934010684229183;
  mat_elem(b, 1, 2) = 0.0119020695012414;
  mat_elem(b, 1, 3) = 0.311215042044805;
  mat_elem(b, 2, 0) = 0.0539501186666072;
  mat_elem(b, 2, 1) = 0.129906208473730;
  mat_elem(b, 2, 2) = 0.337122644398882;
  mat_elem(b, 2, 3) = 0.528533135506213;
  mat_elem(b, 3, 0) = 0.530797553008973;
  mat_elem(b, 3, 1) = 0.568823660872193;
  mat_elem(b, 3, 2) = 0.162182308193243;
  mat_elem(b, 3, 3) = 0.165648729499781;
  double** c = alloc_dmatrix(n, n);
  mat_elem(c, 0, 0) = 0.601981941401637;
  mat_elem(c, 0, 1) = 0.748151592823710;
  mat_elem(c, 0, 2) = 0.913337361501670;
  mat_elem(c, 0, 3) = 0.996134716626886;
  mat_elem(c, 1, 0) = 0.262971284540144;
  mat_elem(c, 1, 1) = 0.450541598502498;
  mat_elem(c, 1, 2) = 0.152378018969223;
  mat_elem(c, 1, 3) = 0.0781755287531837;
  mat_elem(c, 2, 0) = 0.654079098476782;
  mat_elem(c, 2, 1) = 0.0838213779969326;
  mat_elem(c, 2, 2) = 0.825816977489547;
  mat_elem(c, 2, 3) = 0.442678269775446;
  mat_elem(c, 3, 0) = 0.689214503140008;
  mat_elem(c, 3, 1) = 0.228976968716819;
  mat_elem(c, 3, 2) = 0.538342435260057;
  mat_elem(c, 3, 3) = 0.106652770180584;
  double alpha = 2.1;
  double beta = 0.8;

  printf("a: "); fprint_dmatrix(stdout, n, n, a);
  printf("b: "); fprint_dmatrix(stdout, n, n, b);
  printf("c: "); fprint_dmatrix(stdout, n, n, c);

  cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, n, n, n, alpha, mat_ptr(a),
              n, mat_ptr(b), n, beta, mat_ptr(c), n);

  printf("%10.5f * a * b + %10.5f * c: ", alpha, beta); fprint_dmatrix(stdout, n, n, c);

  free_dmatrix(a);
  free_dmatrix(b);
  free_dmatrix(c);
}

/* gemm.m

A = [0.959291425205444,0.257508254123736,0.243524968724989,0.251083857976031;0.547215529963803,0.840717255983663,0.929263623187228,0.616044676146639;0.138624442828679,0.254282178971531,0.349983765984809,0.473288848902729;0.149294005559057,0.814284826068816,0.196595250431208,0.351659507062997]
B = [0.567821640725221,0.779167230102011,0.469390641058206,0.794284540683907;0.0758542895630636,0.934010684229183,0.0119020695012414,0.311215042044805;0.0539501186666072,0.129906208473730,0.337122644398882,0.528533135506213;0.530797553008973,0.568823660872193,0.162182308193243,0.165648729499781]
C = [0.601981941401637,0.748151592823710,0.913337361501670,0.996134716626886;0.262971284540144,0.450541598502498,0.152378018969223,0.0781755287531837;0.654079098476782,0.0838213779969326,0.825816977489547,0.442678269775446;0.689214503140008,0.228976968716819,0.538342435260057,0.106652770180584]
2.1 * A * B + 0.8 * C

*/

max_stars_repo_name: t-sakashita/rokko

max_issues_repo_name: t-sakashita/rokko

max_forks_repo_name: t-sakashita/rokko


text:
/**
 * \file BaseFilter.h
 */

#ifndef ATK_CORE_BASEFILTER_H
#define ATK_CORE_BASEFILTER_H

#include <cstdint>
#include <string>
#include <vector>

#include <boost/dynamic_bitset/dynamic_bitset.hpp>
#include <gsl/gsl>

#include <ATK/config.h>
#include <ATK/Core/config.h>

#if ATK_PROFILING == 1
#include <chrono>
#endif

#if ATK_USE_THREADPOOL == 1
#include <tbb/queuing_mutex.h>
#endif

namespace ATK
{
  /// Base class for all filters
  class BaseFilter
  {
  public:
    BaseFilter(const BaseFilter&) = delete;
    BaseFilter& operator=(const BaseFilter&) = delete;
    /*!
     * @brief Constructor for the core filter
     * @param nb_input_ports is the total number of input ports of this filter
     * @param nb_output_ports is the total number of output ports of this filter
     */
    ATK_CORE_EXPORT BaseFilter(gsl::index nb_input_ports, gsl::index nb_output_ports);
    /// Move constructor
    ATK_CORE_EXPORT BaseFilter(BaseFilter&& other) noexcept;
    /// Destructor
    virtual ATK_CORE_EXPORT ~BaseFilter();

    /*!
     * @brief Links this filter port to another filter's output port
     * @param input_port is the port number where this filter will receive something
     * @param filter is a pointer to the previous filter
     * @param output_port is the port number where this filter will be connected
     */
    ATK_CORE_EXPORT virtual void set_input_port(gsl::index input_port, gsl::not_null<BaseFilter*> filter, gsl::index output_port) = 0;
    ATK_CORE_EXPORT virtual void set_input_port(gsl::index input_port, BaseFilter& filter, gsl::index output_port);

    /// Starts processing after calling reset
    ATK_CORE_EXPORT void process(gsl::index size);
    /// As process, but doesn't call process_impl
    ATK_CORE_EXPORT void dryrun(gsl::index size);

#if ATK_USE_THREADPOOL == 1
    /// Allows threaded processing
    ATK_CORE_EXPORT void process_parallel(gsl::index size);
#endif

    /*!
     * @brief Setup the input sampling rate, must be in sync with the rest of the pipeline
     * @param rate is the input sampling rate of the plugin
     */
    ATK_CORE_EXPORT void set_input_sampling_rate(gsl::index rate);
    /// Returns this filter internal input sampling rate
    ATK_CORE_EXPORT gsl::index get_input_sampling_rate() const;
    /*!
    * @brief Setup the output sampling rate, must be in sync with the rest of the pipeline
    * @param rate is the output sampling rate of the plugin
    */
    ATK_CORE_EXPORT void set_output_sampling_rate(gsl::index rate);
    /// Returns this filter internal output sampling rate
    ATK_CORE_EXPORT gsl::index get_output_sampling_rate() const;

    /// Returns this filter number of input ports
    ATK_CORE_EXPORT gsl::index get_nb_input_ports() const;
    /*!
    * @brief Changes the number of input ports
    * Will trigger a full setup
    * @param nb_ports is the new number of input ports of the plugin
    */
    ATK_CORE_EXPORT virtual void set_nb_input_ports(gsl::index nb_ports);
    /// Returns this filter number of input ports
    ATK_CORE_EXPORT gsl::index get_nb_output_ports() const;
    /*!
    * @brief Changes the number of output ports
    * Will trigger a full setup
    * @param nb_ports is the new number of output ports of the plugin
    */
    ATK_CORE_EXPORT virtual void set_nb_output_ports(gsl::index nb_ports);
    /// Returns this filter input delay (additional pre-0 samples)
    ATK_CORE_EXPORT void set_input_delay(gsl::index delay);
    /// Returns this filter input delay (additional pre-0 samples)
    ATK_CORE_EXPORT gsl::index get_input_delay() const;
    /// Returns this filter output delay (additional pre-0 samples)
    ATK_CORE_EXPORT void set_output_delay(gsl::index delay);
    /// Returns this filter output delay (additional pre-0 samples)
    ATK_CORE_EXPORT gsl::index get_output_delay() const;
    /*!
    * @brief Changes the filter's latency
    * @param latency is the new latency
    */
    ATK_CORE_EXPORT virtual void set_latency(gsl::index latency);
    /// Returns this filter latency
    ATK_CORE_EXPORT gsl::index get_latency() const;
    /// Returns the pipeline global latency from this plugin
    ATK_CORE_EXPORT gsl::index get_global_latency() const;

    /// Resets the internal state of the filter (mandatory before processing a new clip in a DAW for instance)
    ATK_CORE_EXPORT virtual void full_setup();

    /// Resets the filter so that it will process something if needed
    void reset();
    /// Returns the type that the filter processes
    virtual int get_type() const = 0;
    /// Starts processing without calling reset
    template<bool must_process>
    void process_conditionnally(gsl::index size);
#if ATK_USE_THREADPOOL == 1
    /// Starts parallel processing without calling reset
    void process_conditionnally_parallel(gsl::index size);
#endif

  protected:
    /// The actual filter processing part
    virtual void process_impl(gsl::index size) const = 0;

    /// Prepares the filter by retrieving the inputs arrays
    virtual void prepare_process(gsl::index size) = 0;
    /// Prepares the filter by resizing the outputs arrays
    virtual void prepare_outputs(gsl::index size) = 0;
    /// Changes the internal check to allow a disconnected input port
    void allow_inactive_connection(unsigned int port);

    /// Use this call to recompute internal parameters
    ATK_CORE_EXPORT virtual void setup();

    /// Number of input ports
    gsl::index nb_input_ports = 0;
    /// Number of output ports
    gsl::index nb_output_ports = 0;
    /// Input sampling rate of the plugin
    gsl::index input_sampling_rate = 0;
    /// Output sampling rate of the plugin
    gsl::index output_sampling_rate = 0;
    /// The connections to the output pins of some filters
    std::vector<std::pair<gsl::index, BaseFilter*> > connections;

    /// Input delay of the input port
    gsl::index input_delay = 0;
    /// Output delay of the input port
    gsl::index output_delay = 0;

    /// Latency of the plugin
    gsl::index latency = 0;
    /// Last processed size
    gsl::index last_size = 0;

  private:
    boost::dynamic_bitset<> input_mandatory_connection;
    bool is_reset = false;
#if ATK_PROFILING == 1
    std::string class_name;
    std::chrono::steady_clock::duration input_conversion_time;
    std::chrono::steady_clock::duration output_conversion_time;
    std::chrono::steady_clock::duration process_time;
#endif
#if ATK_USE_THREADPOOL == 1
    tbb::queuing_mutex mutex;
#endif
  };
}

#endif

max_stars_repo_name: D-J-Roberts/AudioTK

max_issues_repo_name: D-J-Roberts/AudioTK

max_forks_repo_name: D-J-Roberts/AudioTK


text:
#include <gsl/gsl_randist.h>
#include <time.h>

gsl_rng *r;

void initialiseStatistics()
{
  r = gsl_rng_alloc(gsl_rng_ranlxs0);
  gsl_rng_set(r, time(NULL));
}

void cleanupStatistics() { free(r); }

double denormalise(double standardNormalValue, double mean, double stddev)
{
  return (standardNormalValue * stddev) + mean;
}

double getRandomGaussian(double mean, double stddev, double *min, double *max)
{
  double v = gsl_ran_gaussian(r, 1.0);
  double d = denormalise(v, mean, stddev);

  if (min != NULL)
  {
    if (d <= *min)
    {
      return getRandomGaussian(mean, stddev, min, max);
    }
  }
  if (max != NULL)
  {
    if (d >= *max)
    {
      return getRandomGaussian(mean, stddev, min, max);
    }
  }
  return d;
}

max_stars_repo_name: kunal-mandalia/economic-model

max_issues_repo_name: kunal-mandalia/economic-model

max_forks_repo_name: kunal-mandalia/economic-model


text:
#ifndef SCANMATCHERPROCESSOR_H
#define SCANMATCHERPROCESSOR_H

// #include <log/sensorlog.h>
#include "sensor/sensor_range/rangereading.h"
#include "sensor/sensor_range/rangesensor.h"
//#include <gsl/gsl_eigen.h>
#include "scanmatcher.h"

namespace GMapping {

class ScanMatcherProcessor {
 public:
  ScanMatcherProcessor(const ScanMatcherMap& m);
  ScanMatcherProcessor(double xmin, double ymin, double xmax, double ymax,
                       double delta, double patchdelta, bool use_icp);
  virtual ~ScanMatcherProcessor();
  virtual void processScan(const RangeReading& reading);
  void setSensorMap(const SensorMap& smap, std::string sensorName = "FLASER");
  void init();
  void setMatchingParameters(double urange, double range, double sigma,
                             int kernsize, double lopt, double aopt,
                             int iterations, bool computeCovariance = false);
  void setRegistrationParameters(double regScore, double critScore);
  OrientedPoint getPose() const;
  inline const ScanMatcherMap& getMap() const { return m_map; }
  inline ScanMatcher& matcher() { return m_matcher; }
  inline void setmaxMove(double mmove) { m_maxMove = mmove; }
  bool useICP;

 protected:
  ScanMatcher m_matcher;
  bool m_computeCovariance;
  bool m_first;
  SensorMap m_sensorMap;
  double m_regScore, m_critScore;
  unsigned int m_beams;
  double m_maxMove;
  // state
  ScanMatcherMap m_map;
  OrientedPoint m_pose;
  OrientedPoint m_odoPose;
  int m_count;
  // gsl_eigen_symmv_workspace * m_eigenspace;
};

};  // namespace GMapping

#endif