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

#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 "LLVgeometry.h"
#include "struct.h"
#include "LLVFDresponse.h"
#include "timeconversion.h"

/************************************************************/
/********* Functions setting detector geometry **************/

/* Function setting the response matrix of a given detector, in cartesian coordinates */
void SetMatrixD(
  gsl_matrix* D,                       /* Output: matrix of the detector response Dij */
  const Detectortag tag)               /* Tag identifying the detector */
{
  /* Allocating and defining the cartesian unit vectors along the two arms */
  gsl_vector* nx = gsl_vector_alloc(3);
  gsl_vector* ny = gsl_vector_alloc(3);
  if(tag==LHO) {
    gsl_vector_set(nx, 0, LAL_LHO_4K_ARM_X_DIRECTION_X);
    gsl_vector_set(nx, 1, LAL_LHO_4K_ARM_X_DIRECTION_Y);
    gsl_vector_set(nx, 2, LAL_LHO_4K_ARM_X_DIRECTION_Z);
    gsl_vector_set(ny, 0, LAL_LHO_4K_ARM_Y_DIRECTION_X);
    gsl_vector_set(ny, 1, LAL_LHO_4K_ARM_Y_DIRECTION_Y);
    gsl_vector_set(ny, 2, LAL_LHO_4K_ARM_Y_DIRECTION_Z);
  }
  else if(tag==LLO) {
    gsl_vector_set(nx, 0, LAL_LLO_4K_ARM_X_DIRECTION_X);
    gsl_vector_set(nx, 1, LAL_LLO_4K_ARM_X_DIRECTION_Y);
    gsl_vector_set(nx, 2, LAL_LLO_4K_ARM_X_DIRECTION_Z);
    gsl_vector_set(ny, 0, LAL_LLO_4K_ARM_Y_DIRECTION_X);
    gsl_vector_set(ny, 1, LAL_LLO_4K_ARM_Y_DIRECTION_Y);
    gsl_vector_set(ny, 2, LAL_LLO_4K_ARM_Y_DIRECTION_Z);
  }
  else if(tag==VIRGO) {
    gsl_vector_set(nx, 0, LAL_VIRGO_ARM_X_DIRECTION_X);
    gsl_vector_set(nx, 1, LAL_VIRGO_ARM_X_DIRECTION_Y);
    gsl_vector_set(nx, 2, LAL_VIRGO_ARM_X_DIRECTION_Z);
    gsl_vector_set(ny, 0, LAL_VIRGO_ARM_Y_DIRECTION_X);
    gsl_vector_set(ny, 1, LAL_VIRGO_ARM_Y_DIRECTION_Y);
    gsl_vector_set(ny, 2, LAL_VIRGO_ARM_Y_DIRECTION_Z);
  }
  else {
    printf("Error: invalid detector tag\n");
    exit(1);
  }

  /* Setting the components of the matrix D = 1/2 (nx nx - ny ny) */
  for(int i=0; i<3; i++) {
    for(int j=0; j<3; j++) {
      gsl_matrix_set(D, i, j, 1./2*(gsl_vector_get(nx, i)*gsl_vector_get(nx, j) - gsl_vector_get(ny, i)*gsl_vector_get(ny, j)));
    }
  }
  gsl_vector_free(nx);
  gsl_vector_free(ny);
}

/* Function setting the position of a detector, in cartesian coordinates */
void SetVectorXd(
  gsl_vector* Xd,                      /* Output: position vector of the detector */
  const Detectortag tag)               /* Tag identifying the detector */
{
  if(tag==LHO) {
    gsl_vector_set(Xd, 0, LAL_LHO_4K_VERTEX_LOCATION_X_SI);
    gsl_vector_set(Xd, 1, LAL_LHO_4K_VERTEX_LOCATION_Y_SI);
    gsl_vector_set(Xd, 2, LAL_LHO_4K_VERTEX_LOCATION_Z_SI);
  }
  else if(tag==LLO) {
    gsl_vector_set(Xd, 0, LAL_LLO_4K_VERTEX_LOCATION_X_SI);
    gsl_vector_set(Xd, 1, LAL_LLO_4K_VERTEX_LOCATION_Y_SI);
    gsl_vector_set(Xd, 2, LAL_LLO_4K_VERTEX_LOCATION_Z_SI);
  }
  else if(tag==VIRGO) {
    gsl_vector_set(Xd, 0, LAL_VIRGO_VERTEX_LOCATION_X_SI);
    gsl_vector_set(Xd, 1, LAL_VIRGO_VERTEX_LOCATION_Y_SI);
    gsl_vector_set(Xd, 2, LAL_VIRGO_VERTEX_LOCATION_Z_SI);
  }
  else {
    printf("Error: invalid detector tag\n");
    exit(1);
  }
}

/* Function setting the cartesian coordinates of the wave frame vectors (X,Y,Z), given the position in the sky and polarization */
void SetVectorsXYZ(
  gsl_vector* X,                       /* Output: cartesian vector of the wave frame unit vector X */
  gsl_vector* Y,                       /* Output: cartesian vector of the wave frame unit vector Y */
  gsl_vector* Z,                       /* Output: cartesian vector of the wave frame unit vector Z */
  const double theta,                  /* First angle for the position in the sky (Earth-based spherical angle) */
  const double phi,                    /* Second angle for the position in the sky (Earth-based spherical angle) */
  const double psi)                    /* Polarization angle */
{
  double cospsi = cos(psi);
  double sinpsi = sin(psi);
  double costheta = cos(theta);
  double sintheta = sin(theta);
  double cosphi = cos(phi);
  double sinphi = sin(phi);
  /* Unit vector X */
  gsl_vector_set(X, 0, -cospsi*sinphi+sinpsi*costheta*cosphi);
  gsl_vector_set(X, 1, cospsi*cosphi+sinpsi*costheta*sinphi);
  gsl_vector_set(X, 2, -sinpsi*sintheta);
  /* Unit vector Y */
  gsl_vector_set(Y, 0, sinpsi*sinphi+cospsi*costheta*cosphi);
  gsl_vector_set(Y, 1, -sinpsi*cosphi+cospsi*costheta*sinphi);
  gsl_vector_set(Y, 2, -cospsi*sintheta);
  /* Unit vector Z */
  gsl_vector_set(Z, 0, -sintheta*cosphi);
  gsl_vector_set(Z, 1, -sintheta*sinphi);
  gsl_vector_set(Z, 2, -costheta);
}

/***************************************/
/********* Core functions **************/

/* Function to convert string input network string to Networktag */
Networktag ParseNetworktag(char* string) {
  Networktag tag;
  if(strcmp(string, "L")==0) tag = L;
  else if(strcmp(string, "H")==0) tag = H;
  else if(strcmp(string, "V")==0) tag = V;
  else if(strcmp(string, "LH")==0) tag = LH;
  else if(strcmp(string, "LV")==0) tag = LV;
  else if(strcmp(string, "HV")==0) tag = HV;
  else if(strcmp(string, "LHV")==0) tag = LHV;
  else {
    printf("Error in ParseNetworktag: string not recognized.\n");
    exit(1);
  }
  return tag;
}

/* Function evaluating the Fourier-domain factors for LHV detectors */
/* Note: in case only one or two detectors are considered, amplitudes for the missing ones are set to 0 */
/* Contrarily to the LISA case, here returns trivially 0 or 1 */
/* (allows minimal changes from the old structure that assumed LHV - but not optimal) */
static int EvaluateDetectorFactor3Det(
  double* factor1,                       /* Output for factor for detector 1 */
  double* factor2,                       /* Output for factor for detector 2 */
  double* factor3,                       /* Output for factor for detector 3 */
  const Networktag networktag)                  /* Selector for the detector network */
{
  switch(networktag) {
    case L:
    *factor1 = 1.;
    *factor2 = 0;
    *factor3 = 0;
    break;
    case H:
    *factor1 = 0;
    *factor2 = 1.;
    *factor3 = 0;
    break;
    case V:
    *factor1 = 0;
    *factor2 = 0;
    *factor3 = 1.;
    break;
    case LH:
    *factor1 = 1.;
    *factor2 = 1.;
    *factor3 = 0;
    break;
    case LV:
    *factor1 = 1.;
    *factor2 = 0;
    *factor3 = 1.;
    break;
    case HV:
    *factor1 = 0;
    *factor2 = 1.;
    *factor3 = 1.;
    break;
    case LHV:
    *factor1 = 1.;
    *factor2 = 1.;
    *factor3 = 1.;
    break;
    default:
    printf("Error in EvaluateDetectorFactor3Det: networktag not recognized.\n");
    exit(1);
  }
  return SUCCESS;
}

/* Core function processing a signal (in the form of a list of modes) through the Fourier-domain LLV response (for a given detector), for given values of the inclination, position in the sky and polarization angle */
int LLVSimFDResponse(
  struct tagListmodesCAmpPhaseFrequencySeries **listhlm,  /* Input: list of modes in Frequency-domain amplitude and phase form as produced by the ROM */
  struct tagListmodesCAmpPhaseFrequencySeries **lists,    /* Output: list of contribution of each mode in the detector signal, in Frequency-domain amplitude and phase form, for the given detector and sky position */
  const double gpstime,                                   /* GPS time (s) when the signal at coalescence reaches geocenter */
  const double ra,                                        /* Position in the sky: J2000.0 right ascension (rad) */
  const double dec,                                       /* Position in the sky: J2000.0 declination (rad) */
  const double inclination,                               /* Inclination of the source (rad) */
  const double psi,                                       /* Polarization angle (rad) */
  const Detectortag tag)                                  /* Tag identifying the detector */
{
  /* Define matrix D and position vector Xd of the detector */
  gsl_matrix* D = gsl_matrix_alloc(3,3);
  gsl_vector* Xd = gsl_vector_alloc(3);
  SetVectorXd(Xd, tag);
  SetMatrixD(D, tag);

  /* Conversion from (ra, dec) to the Earth-based spherical angles (theta, phi) - neglecting nutation and precession, and identifying UT1 and UTC, so accurate roughly to a second of time */
  double gmst_angle = gmst_angle_from_gpstime(gpstime);
  double theta = PI/2 - dec;
  double phi = ra - gmst_angle;

  /* Define waveframe unit vectors (X,Y,Z) */
  gsl_vector* X = gsl_vector_alloc(3);
  gsl_vector* Y = gsl_vector_alloc(3);
  gsl_vector* Z = gsl_vector_alloc(3);
  SetVectorsXYZ(X, Y, Z, theta, phi, psi);

  /* Compute the delay from geocenter to the detector */
  double delaylength;
  gsl_blas_ddot(Xd, Z, &delaylength);
  double twopidelay = 2*PI*delaylength/C_SI;

  /* Compute the value of pattern functions Fplus, Fcross */
  gsl_vector* DX = gsl_vector_calloc(3); /* Temporary vector D.X, initialized to 0 */
  gsl_vector* DY = gsl_vector_calloc(3); /* Temporary vector D.Y, initialized to 0 */
  gsl_blas_dgemv( CblasNoTrans, 1., D, X, 0, DX);
  gsl_blas_dgemv( CblasNoTrans, 1., D, Y, 0, DY);
  double XDX;
  double XDY;
  double YDY;
  gsl_blas_ddot(X, DX, &XDX);
  gsl_blas_ddot(X, DY, &XDY);
  gsl_blas_ddot(Y, DY, &YDY);
  double Fplus = XDX - YDY;
  double Fcross = 2*XDY;

  /* Main loop over the modes - goes through all the modes present, stopping when encountering NULL */
  ListmodesCAmpPhaseFrequencySeries* listelement = *listhlm;
  while(listelement) {

    /* Definitions: l,m, frequency series and length */
    int l = listelement->l;
    int m = listelement->m;
    CAmpPhaseFrequencySeries* freqseries = listelement->freqseries;
    gsl_vector* freq = freqseries->freq;
    gsl_vector* amp_real = freqseries->amp_real;
    gsl_vector* amp_imag = freqseries->amp_imag;
    gsl_vector* phase = freqseries->phase;
    int len = (int) freq->size;
    double f;
    double complex camp;
    double complex camps;

    /* Computing the Ylm combined factors for plus and cross for this mode */
    /* Capital Phi is set to 0 by convention */
    double complex Yfactorplus;
    double complex Yfactorcross;
    if (!(l%2)) {
      Yfactorplus = 1./2 * (SpinWeightedSphericalHarmonic(inclination, 0., -2, l, m) + conj(SpinWeightedSphericalHarmonic(inclination, 0., -2, l, -m)));
      Yfactorcross = I/2 * (SpinWeightedSphericalHarmonic(inclination, 0., -2, l, m) - conj(SpinWeightedSphericalHarmonic(inclination, 0., -2, l, -m)));
    }
    else {
      Yfactorplus = 1./2 * (SpinWeightedSphericalHarmonic(inclination, 0., -2, l, m) - conj(SpinWeightedSphericalHarmonic(inclination, 0., -2, l, -m)));
      Yfactorcross = I/2 * (SpinWeightedSphericalHarmonic(inclination, 0., -2, l, m) + conj(SpinWeightedSphericalHarmonic(inclination, 0., -2, l, -m)));
    }

    /* Initializing frequency series structure for this mode, for the signal s = F+ h+ + Fx hx */
    CAmpPhaseFrequencySeries *modefreqseriess = NULL;
    CAmpPhaseFrequencySeries_Init(&modefreqseriess, len);
    gsl_vector* freqs = modefreqseriess->freq;
    gsl_vector* amp_reals = modefreqseriess->amp_real;
    gsl_vector* amp_imags = modefreqseriess->amp_imag;
    gsl_vector* phases = modefreqseriess->phase;

    /* Loop over the frequencies - multiplying by F+, Fx, and add phase due to the delay from geocenter to the detector */
    double complex factorcamp = Fplus*Yfactorplus + Fcross*Yfactorcross;
    for(int j=0; j<len; j++) {
      f = gsl_vector_get(freq, j);
      camp = gsl_vector_get(amp_real, j) + I*gsl_vector_get(amp_imag, j);
      camps = camp * factorcamp;
      gsl_vector_set(amp_reals, j, creal(camps));
      gsl_vector_set(amp_imags, j, cimag(camps));
      gsl_vector_set(phases, j, gsl_vector_get(phase, j) + twopidelay*f);
    }
    /* Copying the vectors of frequencies */
    gsl_vector_memcpy(freqs, freq);

    /* Append the modes to the ouput list-of-modes structures */
    *lists = ListmodesCAmpPhaseFrequencySeries_AddModeNoCopy(*lists, modefreqseriess, l, m);

    /* Going to the next mode in the list */
    listelement = listelement->next;
  }

  /* Cleaning */
  gsl_matrix_free(D);
  gsl_vector_free(Xd);
  gsl_vector_free(X);
  gsl_vector_free(Y);
  gsl_vector_free(Z);
  gsl_vector_free(DX);
  gsl_vector_free(DY);

  return SUCCESS;
}

/* Core function processing a signal (in the form of a list of modes) through the Fourier-domain LLV response for a given detector network, for given values of the inclination, position in the sky and polarization angle */
/* Note: as for now, asssumes the three detectors are L,H,V - amplitudes simply set to 0 in those that are not selected by networktag */
int LLVSimFDResponse3Det(
  struct tagListmodesCAmpPhaseFrequencySeries **list1,    /* Output: list of contributions of each mode in the signal of detector 1, in Frequency-domain amplitude and phase form */
  struct tagListmodesCAmpPhaseFrequencySeries **list2,    /* Output: list of contributions of each mode in the signal of detector 1, in Frequency-domain amplitude and phase form */
  struct tagListmodesCAmpPhaseFrequencySeries **list3,    /* Output: list of contributions of each mode in the signal of detector 1, in Frequency-domain amplitude and phase form */
  struct tagListmodesCAmpPhaseFrequencySeries **listhlm,  /* Input: list of modes in Frequency-domain amplitude and phase form as produced by the ROM */
  const double gpstime,                                   /* GPS time (s) when the signal at coalescence reaches geocenter */
  const double ra,                                        /* Position in the sky: J2000.0 right ascension (rad) */
  const double dec,                                       /* Position in the sky: J2000.0 declination (rad) */
  const double inclination,                               /* Inclination of the source (rad) */
  const double psi,                                       /* Polarization angle (rad) */
  const Networktag tag)                                   /* Tag identifying the network to use */
{
  /* Read which detectors are to be included in the network */
  double factor1 = 0;
  double factor2 = 0;
  double factor3 = 0;
  EvaluateDetectorFactor3Det(&factor1, &factor2, &factor3, tag);

  /* Define matrix D and position vector Xd for the 3 detectors */
  gsl_matrix* D1 = gsl_matrix_alloc(3,3);
  gsl_matrix* D2 = gsl_matrix_alloc(3,3);
  gsl_matrix* D3 = gsl_matrix_alloc(3,3);
  gsl_vector* Xd1 = gsl_vector_alloc(3);
  gsl_vector* Xd2 = gsl_vector_alloc(3);
  gsl_vector* Xd3 = gsl_vector_alloc(3);
  SetVectorXd(Xd1, LHO);
  SetVectorXd(Xd2, LLO);
  SetVectorXd(Xd3, VIRGO);
  SetMatrixD(D1, LHO);
  SetMatrixD(D2, LLO);
  SetMatrixD(D3, VIRGO);

  /* Conversion from (ra, dec) to the Earth-based spherical angles (theta, phi) - neglecting nutation and precession, and identifying UT1 and UTC, so accurate roughly to a second of time */
  double gmst_angle = gmst_angle_from_gpstime(gpstime);
  double theta = PI/2 - dec;
  double phi = ra - gmst_angle;

  /* Define waveframe unit vectors (X,Y,Z) */
  gsl_vector* X = gsl_vector_alloc(3);
  gsl_vector* Y = gsl_vector_alloc(3);
  gsl_vector* Z = gsl_vector_alloc(3);
  SetVectorsXYZ(X, Y, Z, theta, phi, psi);

  /* Compute the delays from geocenter to each detector */
  double delaylength1 = 0;
  double delaylength2 = 0;
  double delaylength3 = 0;
  gsl_blas_ddot(Xd1, Z, &delaylength1);
  gsl_blas_ddot(Xd2, Z, &delaylength2);
  gsl_blas_ddot(Xd3, Z, &delaylength3);
  double twopidelay1 = 2*PI*delaylength1/C_SI;
  double twopidelay2 = 2*PI*delaylength2/C_SI;
  double twopidelay3 = 2*PI*delaylength3/C_SI;

  /* Compute the value of pattern functions Fplus, Fcross for each detector */
  gsl_vector* DX1 = gsl_vector_calloc(3); /* Temporary vector D.X, initialized to 0 */
  gsl_vector* DY1 = gsl_vector_calloc(3); /* Temporary vector D.Y, initialized to 0 */
  gsl_vector* DX2 = gsl_vector_calloc(3); /* Temporary vector D.X, initialized to 0 */
  gsl_vector* DY2 = gsl_vector_calloc(3); /* Temporary vector D.Y, initialized to 0 */
  gsl_vector* DX3 = gsl_vector_calloc(3); /* Temporary vector D.X, initialized to 0 */
  gsl_vector* DY3 = gsl_vector_calloc(3); /* Temporary vector D.Y, initialized to 0 */
  gsl_blas_dgemv( CblasNoTrans, 1., D1, X, 0, DX1);
  gsl_blas_dgemv( CblasNoTrans, 1., D1, Y, 0, DY1);
  gsl_blas_dgemv( CblasNoTrans, 1., D2, X, 0, DX2);
  gsl_blas_dgemv( CblasNoTrans, 1., D2, Y, 0, DY2);
  gsl_blas_dgemv( CblasNoTrans, 1., D3, X, 0, DX3);
  gsl_blas_dgemv( CblasNoTrans, 1., D3, Y, 0, DY3);
  double XDX1 = 0.;
  double XDY1 = 0.;
  double YDY1 = 0.;
  double XDX2 = 0.;
  double XDY2 = 0.;
  double YDY2 = 0.;
  double XDX3 = 0.;
  double XDY3 = 0.;
  double YDY3 = 0.;
  gsl_blas_ddot(X, DX1, &XDX1);
  gsl_blas_ddot(X, DY1, &XDY1);
  gsl_blas_ddot(Y, DY1, &YDY1);
  gsl_blas_ddot(X, DX2, &XDX2);
  gsl_blas_ddot(X, DY2, &XDY2);
  gsl_blas_ddot(Y, DY2, &YDY2);
  gsl_blas_ddot(X, DX3, &XDX3);
  gsl_blas_ddot(X, DY3, &XDY3);
  gsl_blas_ddot(Y, DY3, &YDY3);
  double Fplus1 = XDX1 - YDY1;
  double Fcross1 = 2*XDY1;
  double Fplus2 = XDX2 - YDY2;
  double Fcross2 = 2*XDY2;
  double Fplus3 = XDX3 - YDY3;
  double Fcross3 = 2*XDY3;

  /* Main loop over the modes - goes through all the modes present, stopping when encountering NULL */
  ListmodesCAmpPhaseFrequencySeries* listelement = *listhlm;
  while(listelement) {

    /* Definitions: l,m, frequency series and length */
    int l = listelement->l;
    int m = listelement->m;
    CAmpPhaseFrequencySeries* freqseries = listelement->freqseries;
    gsl_vector* freq = freqseries->freq;
    gsl_vector* amp_real = freqseries->amp_real;
    gsl_vector* amp_imag = freqseries->amp_imag;
    gsl_vector* phase = freqseries->phase;
    int len = (int) freq->size;
    double f;
    double complex camp;
    double complex camp1;
    double complex camp2;
    double complex camp3;

    /* Computing the Ylm combined factors for plus and cross for this mode */
    /* Capital Phi is set to 0 by convention */
    double complex Yfactorplus;
    double complex Yfactorcross;
    if (!(l%2)) {
      Yfactorplus = 1./2 * (SpinWeightedSphericalHarmonic(inclination, 0., -2, l, m) + conj(SpinWeightedSphericalHarmonic(inclination, 0., -2, l, -m)));
      Yfactorcross = I/2 * (SpinWeightedSphericalHarmonic(inclination, 0., -2, l, m) - conj(SpinWeightedSphericalHarmonic(inclination, 0., -2, l, -m)));
    }
    else {
      Yfactorplus = 1./2 * (SpinWeightedSphericalHarmonic(inclination, 0., -2, l, m) - conj(SpinWeightedSphericalHarmonic(inclination, 0., -2, l, -m)));
      Yfactorcross = I/2 * (SpinWeightedSphericalHarmonic(inclination, 0., -2, l, m) + conj(SpinWeightedSphericalHarmonic(inclination, 0., -2, l, -m)));
    }

    /* Initializing frequency series structure for this mode, for the signal s = F+ h+ + Fx hx in each detector */
    CAmpPhaseFrequencySeries *modefreqseries1 = NULL;
    CAmpPhaseFrequencySeries *modefreqseries2 = NULL;
    CAmpPhaseFrequencySeries *modefreqseries3 = NULL;
    CAmpPhaseFrequencySeries_Init(&modefreqseries1, len);
    CAmpPhaseFrequencySeries_Init(&modefreqseries2, len);
    CAmpPhaseFrequencySeries_Init(&modefreqseries3, len);
    gsl_vector* freq1 = modefreqseries1->freq;
    gsl_vector* freq2 = modefreqseries2->freq;
    gsl_vector* freq3 = modefreqseries3->freq;
    gsl_vector* amp_real1 = modefreqseries1->amp_real;
    gsl_vector* amp_real2 = modefreqseries2->amp_real;
    gsl_vector* amp_real3 = modefreqseries3->amp_real;
    gsl_vector* amp_imag1 = modefreqseries1->amp_imag;
    gsl_vector* amp_imag2 = modefreqseries2->amp_imag;
    gsl_vector* amp_imag3 = modefreqseries3->amp_imag;
    gsl_vector* phase1 = modefreqseries1->phase;
    gsl_vector* phase2 = modefreqseries2->phase;
    gsl_vector* phase3 = modefreqseries3->phase;

    /* Loop over the frequencies - multiplying by F+, Fx, and add phase due to the delay from geocenter to the detector */
    double complex factorcamp1 = Fplus1*Yfactorplus + Fcross1*Yfactorcross;
    double complex factorcamp2 = Fplus2*Yfactorplus + Fcross2*Yfactorcross;
    double complex factorcamp3 = Fplus3*Yfactorplus + Fcross3*Yfactorcross;
    for(int j=0; j<len; j++) {
      f = gsl_vector_get(freq, j);
      camp = gsl_vector_get(amp_real, j) + I*gsl_vector_get(amp_imag, j);
      /* Note: include factor123 as detector selectors, which are simpy 0 or 1 whether the detector is in the array or not */
      camp1 = factor1 * camp * factorcamp1;
      camp2 = factor2 * camp * factorcamp2;
      camp3 = factor3 * camp * factorcamp3;
      gsl_vector_set(amp_real1, j, creal(camp1));
      gsl_vector_set(amp_real2, j, creal(camp2));
      gsl_vector_set(amp_real3, j, creal(camp3));
      gsl_vector_set(amp_imag1, j, cimag(camp1));
      gsl_vector_set(amp_imag2, j, cimag(camp2));
      gsl_vector_set(amp_imag3, j, cimag(camp3));
      gsl_vector_set(phase1, j, gsl_vector_get(phase, j) + twopidelay1*f);
      gsl_vector_set(phase2, j, gsl_vector_get(phase, j) + twopidelay2*f);
      gsl_vector_set(phase3, j, gsl_vector_get(phase, j) + twopidelay3*f);
    }
    /* Copying the vectors of frequencies */
    gsl_vector_memcpy(freq1, freq);
    gsl_vector_memcpy(freq2, freq);
    gsl_vector_memcpy(freq3, freq);

    /* Append the modes to the ouput list-of-modes structures */
    *list1 = ListmodesCAmpPhaseFrequencySeries_AddModeNoCopy(*list1, modefreqseries1, l, m);
    *list2 = ListmodesCAmpPhaseFrequencySeries_AddModeNoCopy(*list2, modefreqseries2, l, m);
    *list3 = ListmodesCAmpPhaseFrequencySeries_AddModeNoCopy(*list3, modefreqseries3, l, m);

    /* Going to the next mode in the list */
    listelement = listelement->next;
  }

  /* Cleaning */
  gsl_matrix_free(D1);
  gsl_matrix_free(D2);
  gsl_matrix_free(D3);
  gsl_vector_free(Xd1);
  gsl_vector_free(Xd2);
  gsl_vector_free(Xd3);
  gsl_vector_free(X);
  gsl_vector_free(Y);
  gsl_vector_free(Z);
  gsl_vector_free(DX1);
  gsl_vector_free(DX2);
  gsl_vector_free(DX3);
  gsl_vector_free(DY1);
  gsl_vector_free(DY2);
  gsl_vector_free(DY3);

  return SUCCESS;
}

max_stars_repo_name: titodalcanton/flare

max_issues_repo_name: titodalcanton/flare

max_forks_repo_name: titodalcanton/flare


text:
/* specfunc/gamma_inc.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include "gsl_sf_erf.h"
#include "gsl_sf_exp.h"
#include "gsl_sf_log.h"
#include "gsl_sf_gamma.h"

#include "error.h"

/* The dominant part,
 * D(a,x) := x^a e^(-x) / Gamma(a+1)
 */
static
int
gamma_inc_D(const double a, const double x, gsl_sf_result * result)
{
  if(a < 10.0) {
    double lnr;
    gsl_sf_result lg;
    gsl_sf_lngamma_e(a+1.0, &lg);
    lnr = a * log(x) - x - lg.val;
    result->val = exp(lnr);
    result->err = 2.0 * GSL_DBL_EPSILON * (fabs(lnr) + 1.0) * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    double mu = (x-a)/a;
    double term1;
    gsl_sf_result gstar;
    gsl_sf_result ln_term;
    gsl_sf_log_1plusx_mx_e(mu, &ln_term);  /* log(1+mu) - mu */
    gsl_sf_gammastar_e(a, &gstar);
    term1 = exp(a*ln_term.val)/sqrt(2.0*M_PI*a);
    result->val  = term1/gstar.val;
    result->err  = 2.0 * GSL_DBL_EPSILON * (fabs(a*ln_term.val) + 1.0) * fabs(result->val);
    result->err += gstar.err/fabs(gstar.val) * fabs(result->val);
    return GSL_SUCCESS;
  }
}

/* P series representation.
 */
static
int
gamma_inc_P_series(const double a, const double x, gsl_sf_result * result)
{
  const int nmax = 5000;

  gsl_sf_result D;
  int stat_D = gamma_inc_D(a, x, &D);

  double sum  = 1.0;
  double term = 1.0;
  int n;
  for(n=1; n<nmax; n++) {
    term *= x/(a+n);
    sum  += term;
    if(fabs(term/sum) < GSL_DBL_EPSILON) break;
  }

  result->val  = D.val * sum;
  result->err  = D.err * fabs(sum);
  result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);

  if(n == nmax)
    GSL_ERROR ("error", GSL_EMAXITER);
  else
    return stat_D;
}

/* Q large x asymptotic
 */
static
int
gamma_inc_Q_large_x(const double a, const double x, gsl_sf_result * result)
{
  const int nmax = 5000;

  gsl_sf_result D;
  const int stat_D = gamma_inc_D(a, x, &D);

  double sum  = 1.0;
  double term = 1.0;
  double last = 1.0;
  int n;
  for(n=1; n<nmax; n++) {
    term *= (a-n)/x;
    if(fabs(term/last) > 1.0) break;
    if(fabs(term/sum)  < GSL_DBL_EPSILON) break;
    sum  += term;
    last  = term;
  }

  result->val  = D.val * (a/x) * sum;
  result->err  = D.err * fabs((a/x) * sum);
  result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);

  if(n == nmax)
    GSL_ERROR ("error", GSL_EMAXITER);
  else
    return stat_D;
}

/* Uniform asymptotic for x near a, a and x large.
 * See [Temme, p. 285]
 * FIXME: need c1 coefficient
 */
static
int
gamma_inc_Q_asymp_unif(const double a, const double x, gsl_sf_result * result)
{
  const double rta = sqrt(a);
  const double eps = (x-a)/a;

  gsl_sf_result ln_term;
  const int stat_ln = gsl_sf_log_1plusx_mx_e(eps, &ln_term);  /* log(1+eps) - eps */
  const double eta  = eps * sqrt(-2.0*ln_term.val/(eps*eps));

  gsl_sf_result erfc;

  double R;
  double c0, c1;

  gsl_sf_erfc_e(eta*M_SQRT2*rta, &erfc);

  if(fabs(eps) < GSL_ROOT5_DBL_EPSILON) {
    c0 = -1.0/3.0 + eps*(1.0/12.0 - eps*(23.0/540.0 - eps*(353.0/12960.0 - eps*589.0/30240.0)));
    c1 = 0.0;
  }
  else {
    double rt_term;
    rt_term = sqrt(-2.0 * ln_term.val/(eps*eps));
    c0 = (1.0 - 1.0/rt_term)/eps;
    c1 = 0.0;
  }

  R = exp(-0.5*a*eta*eta)/(M_SQRT2*M_SQRTPI*rta) * (c0 + c1/a);

  result->val  = 0.5 * erfc.val + R;
  result->err  = GSL_DBL_EPSILON * fabs(R * 0.5 * a*eta*eta) + 0.5 * erfc.err;
  result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);

  return stat_ln;
}

/* Continued fraction for Q.
 *
 * Q(a,x) = D(a,x) a/x F(a,x)
 *             1   (1-a)/x  1/x  (2-a)/x   2/x  (3-a)/x
 *  F(a.x) =  ---- ------- ----- -------- ----- -------- ...
 *            1 +   1 +     1 +   1 +      1 +   1 +
 *
 * Uses Gautschi equivalent series method for the CF evaluation.
 *
 * Assumes a != x + 1, so that the first term of the
 * CF recursion is not undefined. This is why we need
 * gamma_inc_Q_CF_protected() below. Based on a problem
 * report by Teemu Ikonen [Tue Oct 10 12:17:19 MDT 2000].
 */
static
int
gamma_inc_Q_CF(const double a, const double x, gsl_sf_result * result)
{
  const int kmax = 5000;

  gsl_sf_result D;
  const int stat_D = gamma_inc_D(a, x, &D);

  double sum  = 1.0;
  double tk   = 1.0;
  double rhok = 0.0;
  int k;

  for(k=1; k<kmax; k++) {
    double ak;
    if(GSL_IS_ODD(k))
      ak = (0.5*(k+1.0)-a)/x;
    else
      ak = 0.5*k/x;
    rhok  = -ak*(1.0 + rhok)/(1.0 + ak*(1.0 + rhok));
    tk   *= rhok;
    sum  += tk;
    if(fabs(tk/sum) < GSL_DBL_EPSILON) break;
  }

  result->val  = D.val * (a/x) * sum;
  result->err  = D.err * fabs((a/x) * sum);
  result->err += GSL_DBL_EPSILON * (2.0 + 0.5*k) * fabs(result->val);

  if(k == kmax)
    GSL_ERROR ("error", GSL_EMAXITER);
  else
    return stat_D;
}

/* See note above for  gamma_inc_Q_CF(). */
static
int
gamma_inc_Q_CF_protected(const double a, const double x, gsl_sf_result * result)
{
  if(fabs(1.0 - a + x) < 2.0*GSL_DBL_EPSILON) {
    /*
     * This is a problem region because when
     * 1.0 - a + x = 0 the first term of the
     * CF recursion is undefined.
     *
     * I missed this condition when I first
     * implemented gamma_inc_Q_CF() function,
     * so now I have to fix it by side-stepping
     * this point, using the recursion relation
     *   Q(a,x) = Q(a-1,x) + x^(a-1) e^(-z) / Gamma(a)
     *          = Q(a-1,x) + D(a-1,x)
     * to lower 'a' by one, giving an a=x point,
     * which is fine.
     */

    gsl_sf_result D;
    gsl_sf_result tmp_CF;

    const int stat_tmp_CF = gamma_inc_Q_CF(a-1.0, x, &tmp_CF);
    const int stat_D = gamma_inc_D(a-1.0, x, &D);
    result->val  = tmp_CF.val + D.val;
    result->err  = tmp_CF.err + D.err;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_ERROR_SELECT_2(stat_tmp_CF, stat_D);
  }
  else {
    return gamma_inc_Q_CF(a, x, result);
  }
}

/* Useful for small a and x. Handles the subtraction analytically.
 */
static
int
gamma_inc_Q_series(const double a, const double x, gsl_sf_result * result)
{
  double term1;  /* 1 - x^a/Gamma(a+1) */
  double sum;    /* 1 + (a+1)/(a+2)(-x)/2! + (a+1)/(a+3)(-x)^2/3! + ... */
  int stat_sum;
  double term2;  /* a temporary variable used at the end */

  {
    /* Evaluate series for 1 - x^a/Gamma(a+1), small a
     */
    const double pg21 = -2.404113806319188570799476;  /* PolyGamma[2,1] */
    const double lnx  = log(x);
    const double el   = M_EULER+lnx;
    const double c1 = -el;
    const double c2 = M_PI*M_PI/12.0 - 0.5*el*el;
    const double c3 = el*(M_PI*M_PI/12.0 - el*el/6.0) + pg21/6.0;
    const double c4 = -0.04166666666666666667
                       * (-1.758243446661483480 + lnx)
                       * (-0.764428657272716373 + lnx)
		       * ( 0.723980571623507657 + lnx)
		       * ( 4.107554191916823640 + lnx);
    const double c5 = -0.0083333333333333333
                       * (-2.06563396085715900 + lnx)
		       * (-1.28459889470864700 + lnx)
		       * (-0.27583535756454143 + lnx)
		       * ( 1.33677371336239618 + lnx)
		       * ( 5.17537282427561550 + lnx);
    const double c6 = -0.0013888888888888889
                       * (-2.30814336454783200 + lnx)
                       * (-1.65846557706987300 + lnx)
                       * (-0.88768082560020400 + lnx)
                       * ( 0.17043847751371778 + lnx)
                       * ( 1.92135970115863890 + lnx)
                       * ( 6.22578557795474900 + lnx);
    const double c7 = -0.00019841269841269841
                       * (-2.5078657901291800 + lnx)
                       * (-1.9478900888958200 + lnx)
                       * (-1.3194837322612730 + lnx)
                       * (-0.5281322700249279 + lnx)
                       * ( 0.5913834939078759 + lnx)
                       * ( 2.4876819633378140 + lnx)
                       * ( 7.2648160783762400 + lnx);
    const double c8 = -0.00002480158730158730
                       * (-2.677341544966400 + lnx)
                       * (-2.182810448271700 + lnx)
                       * (-1.649350342277400 + lnx)
                       * (-1.014099048290790 + lnx)
                       * (-0.191366955370652 + lnx)
                       * ( 0.995403817918724 + lnx)
                       * ( 3.041323283529310 + lnx)
                       * ( 8.295966556941250 + lnx);
    const double c9 = -2.75573192239859e-6
                       * (-2.8243487670469080 + lnx)
                       * (-2.3798494322701120 + lnx)
                       * (-1.9143674728689960 + lnx)
                       * (-1.3814529102920370 + lnx)
                       * (-0.7294312810261694 + lnx)
                       * ( 0.1299079285269565 + lnx)
                       * ( 1.3873333251885240 + lnx)
                       * ( 3.5857258865210760 + lnx)
                       * ( 9.3214237073814600 + lnx);
    const double c10 = -2.75573192239859e-7
                       * (-2.9540329644556910 + lnx)
                       * (-2.5491366926991850 + lnx)
                       * (-2.1348279229279880 + lnx)
                       * (-1.6741881076349450 + lnx)
                       * (-1.1325949616098420 + lnx)
                       * (-0.4590034650618494 + lnx)
                       * ( 0.4399352987435699 + lnx)
                       * ( 1.7702236517651670 + lnx)
                       * ( 4.1231539047474080 + lnx)
                       * ( 10.342627908148680 + lnx);

    term1 = a*(c1+a*(c2+a*(c3+a*(c4+a*(c5+a*(c6+a*(c7+a*(c8+a*(c9+a*c10)))))))));
  }

  {
    /* Evaluate the sum.
     */
    const int nmax = 5000;
    double t = 1.0;
    int n;
    sum = 1.0;

    for(n=1; n<nmax; n++) {
      t *= -x/(n+1.0);
      sum += (a+1.0)/(a+n+1.0)*t;
      if(fabs(t/sum) < GSL_DBL_EPSILON) break;
    }

    if(n == nmax)
      stat_sum = GSL_EMAXITER;
    else
      stat_sum = GSL_SUCCESS;
  }

  term2 = (1.0 - term1) * a/(a+1.0) * x * sum;
  result->val  = term1 + term2;
  result->err  = GSL_DBL_EPSILON * (fabs(term1) + 2.0*fabs(term2));
  result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  return stat_sum;
}

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

int
gsl_sf_gamma_inc_Q_e(const double a, const double x, gsl_sf_result * result)
{
  if(a <= 0.0 || x < 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(x == 0.0) {
    result->val = 1.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(x <= 0.5*a) {
    /* If the series is quick, do that. It is
     * robust and simple.
     */
    gsl_sf_result P;
    int stat_P = gamma_inc_P_series(a, x, &P);
    result->val  = 1.0 - P.val;
    result->err  = P.err;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return stat_P;
  }
  else if(a >= 1.0e+06 && (x-a)*(x-a) < a) {
    /* Then try the difficult asymptotic regime.
     * This is the only way to do this region.
     */
    return gamma_inc_Q_asymp_unif(a, x, result);
  }
  else if(a < 0.2 && x < 5.0) {
    /* Cancellations at small a must be handled
     * analytically; x should not be too big
     * either since the series terms grow
     * with x and log(x).
     */
    return gamma_inc_Q_series(a, x, result);
  }
  else if(a <= x) {
    if(x <= 1.0e+06) {
      /* Continued fraction is excellent for x >~ a.
       * We do not let x be too large when x > a since
       * it is somewhat pointless to try this there;
       * the function is rapidly decreasing for
       * x large and x > a, and it will just
       * underflow in that region anyway. We
       * catch that case in the standard
       * large-x method.
       */
      return gamma_inc_Q_CF(a, x, result);
    }
    else {
      return gamma_inc_Q_large_x(a, x, result);
    }
  }
  else {
    if(a < 0.8*x) {
      /* Continued fraction again. The convergence
       * is a little slower here, but that is fine.
       * We have to trade that off against the slow
       * convergence of the series, which is the
       * only other option.
       */
      return gamma_inc_Q_CF(a, x, result);
    }
    else {
      gsl_sf_result P;
      int stat_P = gamma_inc_P_series(a, x, &P);
      result->val  = 1.0 - P.val;
      result->err  = P.err;
      result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
      return stat_P;
    }
  }
}

int
gsl_sf_gamma_inc_P_e(const double a, const double x, gsl_sf_result * result)
{
  if(a <= 0.0 || x < 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(x == 0.0) {
    result->val = 0.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(x < 20.0 || x < 0.5*a) {
    /* Do the easy series cases. Robust and quick.
     */
    return gamma_inc_P_series(a, x, result);
  }
  else if(a > 1.0e+06 && (x-a)*(x-a) < a) {
    /* Crossover region. Note that Q and P are
     * roughly the same order of magnitude here,
     * so the subtraction is stable.
     */
    gsl_sf_result Q;
    int stat_Q = gamma_inc_Q_asymp_unif(a, x, &Q);
    result->val  = 1.0 - Q.val;
    result->err  = Q.err;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return stat_Q;
  }
  else if(a <= x) {
    /* Q <~ P in this area, so the
     * subtractions are stable.
     */
    gsl_sf_result Q;
    int stat_Q;
    if(a > 0.2*x) {
      stat_Q = gamma_inc_Q_CF(a, x, &Q);
    }
    else {
      stat_Q = gamma_inc_Q_large_x(a, x, &Q);
    }
    result->val  = 1.0 - Q.val;
    result->err  = Q.err;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return stat_Q;
  }
  else {
    if((x-a)*(x-a) < a) {
      /* This condition is meant to insure
       * that Q is not very close to 1,
       * so the subtraction is stable.
       */
      gsl_sf_result Q;
      int stat_Q = gamma_inc_Q_CF_protected(a, x, &Q);
      result->val  = 1.0 - Q.val;
      result->err  = Q.err;
      result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
      return stat_Q;
    }
    else {
      return gamma_inc_P_series(a, x, result);
    }
  }
}

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

#include "eval.h"

double gsl_sf_gamma_inc_P(const double a, const double x)
{
  EVAL_RESULT(gsl_sf_gamma_inc_P_e(a, x, &result));
}

double gsl_sf_gamma_inc_Q(const double a, const double x)
{
  EVAL_RESULT(gsl_sf_gamma_inc_Q_e(a, x, &result));
}

ext: c


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

#ifndef __GSL_VECTOR_COMPLEX_FLOAT_H__
#define __GSL_VECTOR_COMPLEX_FLOAT_H__

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector_float.h>
#include <gsl/gsl_vector_complex.h>
#include <gsl/gsl_block_complex_float.h>

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

__BEGIN_DECLS

typedef struct
{
  size_t size;
  size_t stride;
  float *data;
  gsl_block_complex_float *block;
  int owner;
} gsl_vector_complex_float;

typedef struct
{
  gsl_vector_complex_float vector;
} _gsl_vector_complex_float_view;

typedef _gsl_vector_complex_float_view gsl_vector_complex_float_view;

typedef struct
{
  gsl_vector_complex_float vector;
} _gsl_vector_complex_float_const_view;

typedef const _gsl_vector_complex_float_const_view gsl_vector_complex_float_const_view;

/* Allocation */

gsl_vector_complex_float *gsl_vector_complex_float_alloc (const size_t n);
gsl_vector_complex_float *gsl_vector_complex_float_calloc (const size_t n);

gsl_vector_complex_float *
gsl_vector_complex_float_alloc_from_block (gsl_block_complex_float * b,
                                           const size_t offset,
                                           const size_t n,
                                           const size_t stride);

gsl_vector_complex_float *
gsl_vector_complex_float_alloc_from_vector (gsl_vector_complex_float * v,
                                             const size_t offset,
                                             const size_t n,
                                             const size_t stride);

void gsl_vector_complex_float_free (gsl_vector_complex_float * v);

/* Views */

_gsl_vector_complex_float_view
gsl_vector_complex_float_view_array (float *base,
                                     size_t n);

_gsl_vector_complex_float_view
gsl_vector_complex_float_view_array_with_stride (float *base,
                                                 size_t stride,
                                                 size_t n);

_gsl_vector_complex_float_const_view
gsl_vector_complex_float_const_view_array (const float *base,
                                           size_t n);

_gsl_vector_complex_float_const_view
gsl_vector_complex_float_const_view_array_with_stride (const float *base,
                                                       size_t stride,
                                                       size_t n);

_gsl_vector_complex_float_view
gsl_vector_complex_float_subvector (gsl_vector_complex_float *base,
                                         size_t i,
                                         size_t n);

_gsl_vector_complex_float_view
gsl_vector_complex_float_subvector_with_stride (gsl_vector_complex_float *v,
                                                size_t i,
                                                size_t stride,
                                                size_t n);

_gsl_vector_complex_float_const_view
gsl_vector_complex_float_const_subvector (const gsl_vector_complex_float *base,
                                               size_t i,
                                               size_t n);

_gsl_vector_complex_float_const_view
gsl_vector_complex_float_const_subvector_with_stride (const gsl_vector_complex_float *v,
                                                      size_t i,
                                                      size_t stride,
                                                      size_t n);

_gsl_vector_float_view
gsl_vector_complex_float_real (gsl_vector_complex_float *v);

_gsl_vector_float_view
gsl_vector_complex_float_imag (gsl_vector_complex_float *v);

_gsl_vector_float_const_view
gsl_vector_complex_float_const_real (const gsl_vector_complex_float *v);

_gsl_vector_float_const_view
gsl_vector_complex_float_const_imag (const gsl_vector_complex_float *v);

/* Operations */

gsl_complex_float
gsl_vector_complex_float_get (const gsl_vector_complex_float * v, const size_t i);

void gsl_vector_complex_float_set (gsl_vector_complex_float * v, const size_t i,
                                   gsl_complex_float z);

gsl_complex_float
*gsl_vector_complex_float_ptr (gsl_vector_complex_float * v, const size_t i);

const gsl_complex_float
*gsl_vector_complex_float_const_ptr (const gsl_vector_complex_float * v, const size_t i);

void gsl_vector_complex_float_set_zero (gsl_vector_complex_float * v);
void gsl_vector_complex_float_set_all (gsl_vector_complex_float * v,
                                       gsl_complex_float z);
int gsl_vector_complex_float_set_basis (gsl_vector_complex_float * v, size_t i);

int gsl_vector_complex_float_fread (FILE * stream,
                                    gsl_vector_complex_float * v);
int gsl_vector_complex_float_fwrite (FILE * stream,
                                     const gsl_vector_complex_float * v);
int gsl_vector_complex_float_fscanf (FILE * stream,
                                     gsl_vector_complex_float * v);
int gsl_vector_complex_float_fprintf (FILE * stream,
                                      const gsl_vector_complex_float * v,
                                      const char *format);

int gsl_vector_complex_float_memcpy (gsl_vector_complex_float * dest, const gsl_vector_complex_float * src);

int gsl_vector_complex_float_reverse (gsl_vector_complex_float * v);

int gsl_vector_complex_float_swap (gsl_vector_complex_float * v, gsl_vector_complex_float * w);
int gsl_vector_complex_float_swap_elements (gsl_vector_complex_float * v, const size_t i, const size_t j);

int gsl_vector_complex_float_isnull (const gsl_vector_complex_float * v);
int gsl_vector_complex_float_ispos (const gsl_vector_complex_float * v);
int gsl_vector_complex_float_isneg (const gsl_vector_complex_float * v);

#ifdef HAVE_INLINE

extern inline
gsl_complex_float
gsl_vector_complex_float_get (const gsl_vector_complex_float * v,
                              const size_t i)
{
#if GSL_RANGE_CHECK
  if (i >= v->size)
    {
      gsl_complex_float zero = {{0, 0}};
      GSL_ERROR_VAL ("index out of range", GSL_EINVAL, zero);
    }
#endif
  return *GSL_COMPLEX_FLOAT_AT (v, i);
}

extern inline
void
gsl_vector_complex_float_set (gsl_vector_complex_float * v,
                              const size_t i, gsl_complex_float z)
{
#if GSL_RANGE_CHECK
  if (i >= v->size)
    {
      GSL_ERROR_VOID ("index out of range", GSL_EINVAL);
    }
#endif
  *GSL_COMPLEX_FLOAT_AT (v, i) = z;
}

extern inline
gsl_complex_float *
gsl_vector_complex_float_ptr (gsl_vector_complex_float * v,
                              const size_t i)
{
#if GSL_RANGE_CHECK
  if (i >= v->size)
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return GSL_COMPLEX_FLOAT_AT (v, i);
}

extern inline
const gsl_complex_float *
gsl_vector_complex_float_const_ptr (const gsl_vector_complex_float * v,
                                    const size_t i)
{
#if GSL_RANGE_CHECK
  if (i >= v->size)
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return GSL_COMPLEX_FLOAT_AT (v, i);
}

#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_VECTOR_COMPLEX_FLOAT_H__ */

ext: h

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_statistics_double.h>
#include "tz_error.h"
#include "tz_constant.h"
#include "tz_iarray.h"
#include "tz_stack_attribute.h"
#include "tz_image_lib.h"
#include "tz_objdetect.h"
#include "tz_imatrix.h"
#include "tz_stack_math.h"
#include "tz_stack_bwdist.h"
#include "tz_voxel_linked_list.h"
#include "tz_stack_sampling.h"
#include "tz_matlabio.h"
#include "tz_darray.h"
#include "tz_stack_draw.h"
#include "tz_stack_code.h"
#include "tz_stack_bwmorph.h"
#include "tz_stack.h"
#include "tz_voxel_graphics.h"

INIT_EXCEPTION_MAIN(e)

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

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

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

  char file_path[100];
  sprintf(file_path, "../data/%s/mask.tif", neuron_name);

  Stack *stack = Read_Stack(file_path);
  //Stack_Not(stack, stack);

  //Set_Matlab_Path(matlab_path);
  //Stack *dist = Stack_Bwdist(stack);
  /*
  Objlabel_Workspace *ow = New_Objlabel_Workspace();
  ow->conn = 18;
  Stack_Not(stack, stack);
  Stack *dist = Stack_Boundary_Code(stack, NULL, ow);
  Stack_Not(stack, stack);
  */

  Stack *dist = Stack_Bwdist_L_U16(stack, NULL, 1);

  sprintf(file_path, "../data/%s/dist.tif", neuron_name);
  Write_Stack(file_path, dist);
  printf("Distance map %s created.\n", file_path);

  /*
  Write_Stack(file_path,
	      Scale_Double_Stack((double *)dist->array, stack->width,
				 stack->height, stack->depth, GREY));
  */

  Kill_Stack(stack);

  //Stack *seeds = Stack_Local_Max(dist, NULL, STACK_LOCMAX_FLAT);
  //Stack_Clean_Locmax(dist, seeds);

  Stack *seeds = Stack_Locmax_Region(dist, 26);
  //Write_Stack("../data/test.tif", seeds);

  Object_3d_List *objs = Stack_Find_Object_N(seeds, NULL, 1, 0, 26);
  Zero_Stack(seeds);
  int objnum = 0;
  while (objs != NULL) {
    Object_3d *obj = objs->data;
    Voxel_t center;
    Object_3d_Central_Voxel(obj, center);
    Set_Stack_Pixel(seeds, center[0], center[1], center[2], 0, 1);
    objs = objs->next;
    objnum++;
  }
  printf("objnum: %d\n", objnum);
  sprintf(file_path, "../data/%s/seeds.tif", neuron_name);
  Write_Stack(file_path, seeds);
  printf("%s (seed mask) created.\n", file_path);

  /*
  int max_idx = 0;
  double max_dist = darray_max((double *) (dist->array),
			       Stack_Voxel_Number(dist),
			       &max_idx);
  */

  Voxel_List *list = Stack_To_Voxel_List(seeds);
  printf("%d\n", Voxel_List_Length(list));
  int *seed_offset = (int *) malloc(sizeof(int) * Voxel_List_Length(list));
  int i;
  int idx = Voxel_List_Length(list) - 1;
  int stack_length = Stack_Voxel_Number(seeds);
  for (i = 0; i < stack_length; i++) {
    if (seeds->array[i] > 0) {
      seed_offset[idx--] = i;
    }
  }

  if (idx != -1) {
    TZ_WARN(ERROR_DATA_VALUE);
  }

  sprintf(file_path, "../data/%s/seed_offset.ar", neuron_name);
  iarray_write(file_path, seed_offset, Voxel_List_Length(list));
  printf("%s (seed indices) created.\n", file_path);

  Pixel_Array *pa = Voxel_List_Sampling(dist, list);

  sprintf(file_path, "../data/%s/seeds.pa", neuron_name);
  Pixel_Array_Write(file_path, pa);
  printf("%s (seed values d^2) created.\n", file_path);

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

  sprintf(file_path, "../data/%s/seeds.ar", neuron_name);
  darray_write(file_path, pa_array, pa->size);
  printf("%s (seed values d) created.\n", file_path);

  /*
  printf("mean: %g, std: %g\n", gsl_stats_mean(pa_array, 1, pa->size),
	 sqrt(gsl_stats_variance(pa_array, 1, pa->size)));
  printf("%g\n", gsl_stats_max(pa_array, 1, pa->size));
  */
  //int max_idx2;
  //double max_dist2 = darray_max(pa_array, pa->size, &max_idx2);

  int n = Stack_Voxel_Number(seeds);
  //double *dist_array = (double *) dist->array;
  uint16 *dist_array = (uint16 *) dist->array;
  for (i = 0; i < n; i++) {
    if (dist_array[i] == 1) {
      seeds->array[i] = 0;
    }
  }
  sprintf(file_path, "../data/%s/large_seeds.tif", neuron_name);
  Write_Stack(file_path, seeds);
  printf("%s (seeds for d > 1) created.\n", file_path);

  Kill_Stack(dist);
  Kill_Stack(seeds);
  free(seed_offset);
  Kill_Pixel_Array(pa);

  return 0;
}


text:
/*
 * global.h
 *
 *  Created on: Apr 5, 2013
 *      Authors: vgomez, Sep Thijssen
 */

#ifndef GLOBAL_H_
#define GLOBAL_H_

#include <vector>
#include <math.h>
#include <iostream>

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

#include <ctime>

#include <hal_quadrotor/State.h>

#define SIZE_MSG	5

using namespace std;

// Message structure
//	IIIII.XXX.xxx.YYY.yyy.ZZZ.zzz.UUU.uuu.VVV.vvv.WWW.vvv
//	IIIII		identifier of the UVA
//	XXX.xxx	position x
//	YYY.yyy	position y
//	ZZZ.zzz	position z
//	UUU.uuu	velocity x
//	VVV.vvv	velocity y
//	WWW.www	velocity z

typedef vector<double> vec;
typedef vector<vec> vvec;

template<typename G>
ostream& operator<<(ostream& os, const vector<G>& v)
{
	typename vector<G>::const_iterator it;
	for (it=v.begin(); it!=v.end(); it++)
		os << *it << " ";
	os << endl;
	return os;
}

string getState(const hal_quadrotor::State::ConstPtr& msg);

#endif /* GLOBAL_H_ */


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

int
main (void)
{
  gsl_permutation * p = gsl_permutation_alloc (3);

  gsl_permutation_init (p);

  do
   {
      gsl_permutation_fprintf (stdout, p, " %u");
      printf ("\n");
   }
  while (gsl_permutation_next(p) == GSL_SUCCESS);

  gsl_permutation_free (p);

  return 0;
}

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#ifndef OPENMC_CELL_H
#define OPENMC_CELL_H

#include <cstdint>
#include <functional> // for hash
#include <limits>
#include <memory> // for unique_ptr
#include <string>
#include <unordered_map>
#include <vector>

#include <gsl/gsl>
#include "hdf5.h"
#include "pugixml.hpp"
#include "dagmc.h"

#include "openmc/constants.h"
#include "openmc/neighbor_list.h"
#include "openmc/position.h"
#include "openmc/surface.h"

namespace openmc {

//==============================================================================
// Constants
//==============================================================================

enum class Fill {
  MATERIAL,
  UNIVERSE,
  LATTICE
};

// TODO: Convert to enum
constexpr int32_t OP_LEFT_PAREN   {std::numeric_limits<int32_t>::max()};
constexpr int32_t OP_RIGHT_PAREN  {std::numeric_limits<int32_t>::max() - 1};
constexpr int32_t OP_COMPLEMENT   {std::numeric_limits<int32_t>::max() - 2};
constexpr int32_t OP_INTERSECTION {std::numeric_limits<int32_t>::max() - 3};
constexpr int32_t OP_UNION        {std::numeric_limits<int32_t>::max() - 4};

//==============================================================================
// Global variables
//==============================================================================

class Cell;
class ParentCell;
class CellInstance;
class Universe;
class UniversePartitioner;

namespace model {
  extern std::unordered_map<int32_t, int32_t> cell_map;
  extern std::vector<std::unique_ptr<Cell>> cells;

  extern std::unordered_map<int32_t, int32_t> universe_map;
  extern std::vector<std::unique_ptr<Universe>> universes;
} // namespace model

//==============================================================================
//! A geometry primitive that fills all space and contains cells.
//==============================================================================

class Universe
{
public:
  int32_t id_;                  //!< Unique ID
  std::vector<int32_t> cells_;  //!< Cells within this universe

  //! \brief Write universe information to an HDF5 group.
  //! \param group_id An HDF5 group id.
  void to_hdf5(hid_t group_id) const;

  BoundingBox bounding_box() const;

  std::unique_ptr<UniversePartitioner> partitioner_;
};

//==============================================================================
//==============================================================================

class Cell {
public:
  //----------------------------------------------------------------------------
  // Constructors, destructors, factory functions

  explicit Cell(pugi::xml_node cell_node);
  Cell() {};
  virtual ~Cell() = default;

  //----------------------------------------------------------------------------
  // Methods

  //! \brief Determine if a cell contains the particle at a given location.
  //!
  //! The bounds of the cell are detemined by a logical expression involving
  //! surface half-spaces. At initialization, the expression was converted
  //! to RPN notation.
  //!
  //! The function is split into two cases, one for simple cells (those
  //! involving only the intersection of half-spaces) and one for complex cells.
  //! Simple cells can be evaluated with short circuit evaluation, i.e., as soon
  //! as we know that one half-space is not satisfied, we can exit. This
  //! provides a performance benefit for the common case. In
  //! contains_complex, we evaluate the RPN expression using a stack, similar to
  //! how a RPN calculator would work.
  //! \param r The 3D Cartesian coordinate to check.
  //! \param u A direction used to "break ties" the coordinates are very
  //!   close to a surface.
  //! \param on_surface The signed index of a surface that the coordinate is
  //!   known to be on.  This index takes precedence over surface sense
  //!   calculations.
  virtual bool
  contains(Position r, Direction u, int32_t on_surface) const = 0;

  //! Find the oncoming boundary of this cell.
  virtual std::pair<double, int32_t>
  distance(Position r, Direction u, int32_t on_surface, Particle* p) const = 0;

  //! Write all information needed to reconstruct the cell to an HDF5 group.
  //! \param group_id An HDF5 group id.
  virtual void to_hdf5(hid_t group_id) const = 0;

  //! Get the BoundingBox for this cell.
  virtual BoundingBox bounding_box() const = 0;

  //----------------------------------------------------------------------------
  // Accessors

  //! Get the temperature of a cell instance
  //! \param[in] instance Instance index. If -1 is given, the temperature for
  //!   the first instance is returned.
  //! \return Temperature in [K]
  double temperature(int32_t instance = -1) const;

  //! Set the temperature of a cell instance
  //! \param[in] T Temperature in [K]
  //! \param[in] instance Instance index. If -1 is given, the temperature for
  //!   all instances is set.
  //! \param[in] set_contained If this cell is not filled with a material,
  //!   collect all contained cells with material fills and set their
  //!   temperatures.
  void set_temperature(double T, int32_t instance = -1, bool set_contained = false);

  //! Get the name of a cell
  //! \return Cell name
  const std::string& name() const { return name_; };

  //! Set the temperature of a cell instance
  //! \param[in] name Cell name
  void set_name(const std::string& name) { name_ = name; };

  //! Get all cell instances contained by this cell
  //! \return Map with cell indexes as keys and instances as values
  std::unordered_map<int32_t, std::vector<int32_t>>
  get_contained_cells() const;

protected:
  void
  get_contained_cells_inner(std::unordered_map<int32_t, std::vector<int32_t>>& contained_cells,
                            std::vector<ParentCell>& parent_cells) const;

public:
  //----------------------------------------------------------------------------
  // Data members

  int32_t id_;                //!< Unique ID
  std::string name_;          //!< User-defined name
  Fill type_;                  //!< Material, universe, or lattice
  int32_t universe_;          //!< Universe # this cell is in
  int32_t fill_;              //!< Universe # filling this cell
  int32_t n_instances_{0};    //!< Number of instances of this cell

  //! \brief Index corresponding to this cell in distribcell arrays
  int distribcell_index_{C_NONE};

  //! \brief Material(s) within this cell.
  //!
  //! May be multiple materials for distribcell.
  std::vector<int32_t> material_;

  //! \brief Temperature(s) within this cell.
  //!
  //! The stored values are actually sqrt(k_Boltzmann * T) for each temperature
  //! T. The units are sqrt(eV).
  std::vector<double> sqrtkT_;

  //! Definition of spatial region as Boolean expression of half-spaces
  std::vector<std::int32_t> region_;
  //! Reverse Polish notation for region expression
  std::vector<std::int32_t> rpn_;
  bool simple_;  //!< Does the region contain only intersections?

  //! \brief Neighboring cells in the same universe.
  NeighborList neighbors_;

  Position translation_ {0, 0, 0}; //!< Translation vector for filled universe

  //! \brief Rotational tranfsormation of the filled universe.
  //
  //! The vector is empty if there is no rotation. Otherwise, the first 9 values
  //! give the rotation matrix in row-major order. When the user specifies
  //! rotation angles about the x-, y- and z- axes in degrees, these values are
  //! also present at the end of the vector, making it of length 12.
  std::vector<double> rotation_;

  std::vector<int32_t> offset_;  //!< Distribcell offset table
};

struct CellInstanceItem {
  int32_t index {-1};        //! Index into global cells array
  int     lattice_indx{-1};  //! Flat index value of the lattice cell
};

//==============================================================================

class CSGCell : public Cell
{
public:
  CSGCell();

  explicit CSGCell(pugi::xml_node cell_node);

  bool
  contains(Position r, Direction u, int32_t on_surface) const;

  std::pair<double, int32_t>
  distance(Position r, Direction u, int32_t on_surface, Particle* p) const;

  void to_hdf5(hid_t group_id) const;

  BoundingBox bounding_box() const;

protected:
  bool contains_simple(Position r, Direction u, int32_t on_surface) const;
  bool contains_complex(Position r, Direction u, int32_t on_surface) const;
  BoundingBox bounding_box_simple() const;
  static BoundingBox bounding_box_complex(std::vector<int32_t> rpn);

  //! Applies DeMorgan's laws to a section of the RPN
  //! \param start Starting point for token modification
  //! \param stop Stopping point for token modification
  static void apply_demorgan(std::vector<int32_t>::iterator start,
                             std::vector<int32_t>::iterator stop);

  //! Removes complement operators from the RPN
  //! \param rpn The rpn to remove complement operators from.
  static void remove_complement_ops(std::vector<int32_t>& rpn);

  //! Returns the beginning position of a parenthesis block (immediately before
  //! two surface tokens) in the RPN given a starting position at the end of
  //! that block (immediately after two surface tokens)
  //! \param start Starting position of the search
  //! \param rpn The rpn being searched
  static std::vector<int32_t>::iterator
  find_left_parenthesis(std::vector<int32_t>::iterator start,
                        const std::vector<int32_t>& rpn);

};

//==============================================================================

#ifdef DAGMC
class DAGCell : public Cell
{
public:
  DAGCell();

  bool contains(Position r, Direction u, int32_t on_surface) const;

  std::pair<double, int32_t>
  distance(Position r, Direction u, int32_t on_surface, Particle* p) const;

  BoundingBox bounding_box() const;

  void to_hdf5(hid_t group_id) const;

  moab::DagMC* dagmc_ptr_; //!< Pointer to DagMC instance
  int32_t dag_index_;      //!< DagMC index of cell
};
#endif

//==============================================================================
//! Speeds up geometry searches by grouping cells in a search tree.
//
//! Currently this object only works with universes that are divided up by a
//! bunch of z-planes.  It could be generalized to other planes, cylinders,
//! and spheres.
//==============================================================================

class UniversePartitioner
{
public:
  explicit UniversePartitioner(const Universe& univ);

  //! Return the list of cells that could contain the given coordinates.
  const std::vector<int32_t>& get_cells(Position r, Direction u) const;

private:
  //! A sorted vector of indices to surfaces that partition the universe
  std::vector<int32_t> surfs_;

  //! Vectors listing the indices of the cells that lie within each partition
  //
  //! There are n+1 partitions with n surfaces.  `partitions_.front()` gives the
  //! cells that lie on the negative side of `surfs_.front()`.
  //! `partitions_.back()` gives the cells that lie on the positive side of
  //! `surfs_.back()`.  Otherwise, `partitions_[i]` gives cells sandwiched
  //! between `surfs_[i-1]` and `surfs_[i]`.
  std::vector<std::vector<int32_t>> partitions_;
};

//==============================================================================
//! Define a containing (parent) cell
//==============================================================================

struct ParentCell {
  gsl::index cell_index;
  gsl::index lattice_index;
};

//==============================================================================
//! Define an instance of a particular cell
//==============================================================================

struct CellInstance {
  //! Check for equality
  bool operator==(const CellInstance& other) const
  { return index_cell == other.index_cell && instance == other.instance; }

  gsl::index index_cell;
  gsl::index instance;
};

struct CellInstanceHash {
  std::size_t operator()(const CellInstance& k) const
  {
    return 4096*k.index_cell + k.instance;
  }
};

//==============================================================================
// Non-member functions
//==============================================================================

void read_cells(pugi::xml_node node);

#ifdef DAGMC
int32_t next_cell(DAGCell* cur_cell, DAGSurface* surf_xed);
#endif

} // namespace openmc
#endif // OPENMC_CELL_H

max_stars_repo_name: Hit-Weixg/openmc

max_issues_repo_name: sinkaa-ife/openmc

max_forks_repo_name: sinkaa-ife/openmc


text:
#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION

#include <Python.h>
#include <numpy/arrayobject.h>

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

#include <math.h>

// Forward function declaration.
static PyObject *zsampler_sample_region(PyObject *self, PyObject *args);
static PyObject *zsampler_sample_categorical(PyObject *self, PyObject *args);

// Boilerplate: method list.
static PyMethodDef methods[] = {
  { "sample_region", zsampler_sample_region, METH_VARARGS, "Doc string."},
  { "sample_bulk_categorical", zsampler_sample_categorical, METH_VARARGS, "Doc string."},
  { NULL, NULL, 0, NULL } /* Sentinel */
};

static struct PyModuleDef sampler_module =
{
    PyModuleDef_HEAD_INIT,
    "zsampler", /* name of module */
    "",          /* module documentation, may be NULL */
    -1,          /* size of per-interpreter state of the module, or -1 if the module keeps state in global variables. */
    methods
};

PyMODINIT_FUNC PyInit_zsampler(void)
{
    import_array(); // crucial command if NumPy API is used
    return PyModule_Create(&sampler_module);
}

/*****************************************************************************
 * Helper methods                                                            *
 *****************************************************************************/
void print_int_array(npy_int64* array, int length) {
  for(int i = 0; i < length; i++) {
    printf("%ld ", array[i]);
  }
  printf("\n");
}

void print_f_array(npy_float64* array, int length) {
  for(int i = 0; i < length; i++) {
    printf("%lf ", array[i]);
  }
  printf("\n");
}

static PyObject* zsampler_sample_categorical(PyObject *self, PyObject *args) {

    PyArrayObject *py_Z, *py_prob;
    if (!PyArg_ParseTuple(args, "O!O!",
                            &PyArray_Type, &py_Z,
                            &PyArray_Type, &py_prob)) {
        return NULL;
    }

    npy_int64 C_N, C_K;
    C_N = PyArray_SHAPE(py_prob)[0];
    C_K = PyArray_SHAPE(py_prob)[1];

    const gsl_rng_type * C_T;
    gsl_rng * C_r;
    gsl_rng_env_setup();
    C_T = gsl_rng_taus;
    C_r = gsl_rng_alloc (C_T);

    double *_prob; // probabilities for sampling from multinomial
    _prob = (double*)malloc(C_K * sizeof(double));

    unsigned int *_new_sample;
    _new_sample = (unsigned int*)malloc(C_K * sizeof(unsigned int));

    PyArrayObject* out_Z = (PyArrayObject*)PyArray_NewCopy(py_Z, NPY_ANYORDER);
    for (int i=0; i < C_N; i++) {

        for (int k=0; k < C_K; k++) {
            _prob[k] = (*(npy_float64*)PyArray_GETPTR2(py_prob, i, k));
         }

        gsl_ran_multinomial(C_r, (size_t)C_K, 1, _prob, _new_sample);

        for (int k = 0; k < C_K; k++) {
            npy_int64* Zk = (npy_int64*)PyArray_GETPTR2(out_Z, i, k);
            Zk[0] = (npy_int64)_new_sample[k];
         }
    }

    free(_prob);
    free(_new_sample);
    gsl_rng_free(C_r);
    return (PyObject*)out_Z;  // no need to IncRef out_Z
}

//      Mixture weight priors is block-specific
static PyObject* zsampler_sample_region(PyObject *self, PyObject *args) {

  // Input variables
  PyArrayObject *py_Z, *py_counts, *py_covariates, *py_beta, *py_alpha;
  PyArrayObject *py_region_alloc, *py_region_label_counts;

  /*  parse single numpy array argument */
    if (!PyArg_ParseTuple(args, "O!O!O!O!O!O!O!",
                            &PyArray_Type, &py_Z,
                            &PyArray_Type, &py_counts,
                            &PyArray_Type, &py_covariates,
                            &PyArray_Type, &py_beta,
                            &PyArray_Type, &py_alpha,
                            &PyArray_Type, &py_region_alloc,
                            &PyArray_Type, &py_region_label_counts)) {
        return NULL;
    }

    Py_INCREF(py_region_alloc);
    Py_INCREF(py_region_label_counts);

    const gsl_rng_type * C_T;
    gsl_rng * C_r;
    gsl_rng_env_setup();
    C_T = gsl_rng_taus;
    C_r = gsl_rng_alloc (C_T);

    // Temporary variables
    npy_int64 C_N, C_K, C_J;
    C_N = PyArray_SHAPE(py_Z)[0];
    C_K = PyArray_SHAPE(py_Z)[1];
    C_J = PyArray_SHAPE(py_beta)[0];

    PyArrayObject* out_Z = (PyArrayObject*)PyArray_NewCopy(py_Z, NPY_ANYORDER);

    npy_int64 *C_counts = (npy_int64 *) PyArray_DATA(py_counts);
    npy_float64 *C_alpha = (npy_float64 *) PyArray_DATA(py_alpha);
    npy_int64 *C_region_alloc = (npy_int64 *) PyArray_DATA(py_region_alloc);

    npy_int64 C_labelsum;
    npy_int64 *C_label_counts;
    C_label_counts = (npy_int64*)malloc(C_K * sizeof(npy_int64));
    double *_mu;
    _mu = (double*)malloc(C_K * sizeof(double));
    double *_logpoi;
    _logpoi = (double*)malloc(C_K * sizeof(double));
    double *_logcat;
    _logcat =(double*)malloc(C_K * sizeof(double));
    double *_lik;
    _lik =(double*)malloc(C_K * sizeof(double));
    double *_prob; // probabilities for sampling from multinomial
    _prob = (double*)malloc(C_K * sizeof(double));
    unsigned int *_new_sample;
    _new_sample = (unsigned int*)malloc(C_K * sizeof(unsigned int));

    npy_float64 C_alpha_sum = 0.0;
    for(int k=0; k < C_K; k++) {
      C_alpha_sum = C_alpha_sum + C_alpha[k];
    }

    for (int i=0; i < C_N; i++) {

      C_labelsum = 0;
      for (int k = 0; k < C_K; k++) {
        // subtract current cell counts from the region counts
        npy_int64* accessor = (npy_int64*)PyArray_GETPTR2(py_region_label_counts, C_region_alloc[i], k);
        accessor[0] = accessor[0] - (*(npy_int64*)PyArray_GETPTR2(out_Z, i, k));  // TODO: make this simpler

        C_label_counts[k] = accessor[0];
        C_labelsum = C_labelsum + C_label_counts[k];
      }

      // computation of likelihood and interim calculations
      for (int k=0; k < C_K; k++) {
        _mu[k] = 0;
        for (int j=0; j < C_J; j++) {
            _mu[k] = _mu[k] + ((*(npy_float64*)PyArray_GETPTR2(py_covariates, i, j)) * (*(npy_float64*)PyArray_GETPTR2(py_beta, j, k)));
        }

        _logpoi[k] = ((double)C_counts[i])*_mu[k] - exp(_mu[k]);
        _logcat[k] = log(((double)C_label_counts[k] + C_alpha[k])  / ((double)C_labelsum + C_alpha_sum));

        _lik[k] = (_logpoi[k] + _logcat[k]);
      }

      double maxval = _lik[0];
      for (int k=0; k < C_K; k++) {
        if(_lik[k] > maxval) {
          maxval = _lik[k];
        }
      }

      double sumexp = 0;
      for (int k=0; k < C_K; k++) {
        sumexp = sumexp + exp((double)_lik[k] - maxval);
      }

      double logsumexp = maxval + log(sumexp);

      for (int k=0; k < C_K; k++) {
        _prob[k] = exp((double)_lik[k] - logsumexp);
      }

      gsl_ran_multinomial(C_r, (size_t)C_K, 1, _prob, _new_sample);
      for (int k = 0; k < C_K; k++) {
        npy_int64* Zk = (npy_int64*)PyArray_GETPTR2(out_Z, i, k);
        Zk[0] = (npy_int64)_new_sample[k];

        // update the region label counts
        npy_int64* accessor = (npy_int64*)PyArray_GETPTR2(py_region_label_counts, C_region_alloc[i], k);
        accessor[0] = accessor[0] + Zk[0];
      }
    }

    free(C_label_counts);
    free(_mu);
    free(_logpoi);
    free(_logcat);
    free(_lik);
    free(_prob);
    free(_new_sample);

    gsl_rng_free(C_r);
    Py_DECREF(py_region_alloc);
    Py_DECREF(py_region_label_counts);
    return (PyObject*)out_Z;  // no need to IncRef out_Z
}


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

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

/*
gsl_movstat_apply_accum()
  Apply moving window statistic to input vector. This is a generalized
routine to handle window endpoints and apply a given accumulator to
the input vector.

Inputs: endtype      - end point handling criteria
        x            - input vector, size n
        accum        - accumulator to apply moving window statistic
        accum_params - parameters to pass to accumulator
        y            - output vector, size n
        z            - second output vector (i.e. minmax), size n; can be NULL
        w            - workspace

Notes:
1) It is allowed to have x = y for in-place moving statistics
*/

int
gsl_movstat_apply_accum(const gsl_movstat_end_t endtype,
                        const gsl_vector * x,
                        const gsl_movstat_accum * accum,
                        void * accum_params,
                        gsl_vector * y,
                        gsl_vector * z,
                        gsl_movstat_workspace * w)
{
  if (x->size != y->size)
    {
      GSL_ERROR("input and output vectors must have same length", GSL_EBADLEN);
    }
  else if (z != NULL && x->size != z->size)
    {
      GSL_ERROR("input and output vectors must have same length", GSL_EBADLEN);
    }
  else
    {
      const int n = (int) x->size;
      const int H = w->H; /* number of samples to left of current sample */
      const int J = w->J; /* number of samples to right of current sample */
      int i;
      double x1 = 0.0;    /* pad values for data edges */
      double xN = 0.0;
      double result[2];

      /* initialize accumulator */
      (accum->init)(w->K, w->state);

      /* pad initial window if necessary */
      if (endtype != GSL_MOVSTAT_END_TRUNCATE)
        {
          if (endtype == GSL_MOVSTAT_END_PADZERO)
            {
              x1 = 0.0;
              xN = 0.0;
            }
          else if (endtype == GSL_MOVSTAT_END_PADVALUE)
            {
              x1 = gsl_vector_get(x, 0);
              xN = gsl_vector_get(x, n - 1);
            }

          /* pad initial windows with H values */
          for (i = 0; i < H; ++i)
            (accum->insert)(x1, w->state);
        }
      else if (accum->delete_oldest == NULL) /* FIXME XXX */
        {
          /* save last K - 1 samples of x for later (needed for in-place input/output) */
          int idx1 = GSL_MAX(n - J - H, 0);
          int idx2 = n - 1;

          for (i = idx1; i <= idx2; ++i)
            w->work[i - idx1] = gsl_vector_get(x, i);
        }

      /* process input vector and fill y(0:n - J - 1) */
      for (i = 0; i < n; ++i)
        {
          double xi = gsl_vector_get(x, i);
          int idx = i - J;

          (accum->insert)(xi, w->state);

          if (idx >= 0)
            {
              (accum->get)(accum_params, result, w->state);
              gsl_vector_set(y, idx, result[0]);

              if (z != NULL)
                gsl_vector_set(z, idx, result[1]);
            }
        }

      if (endtype == GSL_MOVSTAT_END_TRUNCATE)
        {
          /* need to fill y(n-J:n-1) using shrinking windows */
          int idx1 = GSL_MAX(n - J, 0);
          int idx2 = n - 1;

          if (accum->delete_oldest == NULL)
            {
              int wsize = n - GSL_MAX(n - J - H, 0); /* size of work array */

              for (i = idx1; i <= idx2; ++i)
                {
                  int nsamp = n - GSL_MAX(i - H, 0); /* number of samples in this window */
                  int j;

                  (accum->init)(w->K, w->state);

                  for (j = wsize - nsamp; j < wsize; ++j)
                    (accum->insert)(w->work[j], w->state);

                  /* yi = acc_get [ work(i:K-2) ] */
                  (accum->get)(accum_params, result, w->state);
                  gsl_vector_set(y, i, result[0]);

                  if (z != NULL)
                    gsl_vector_set(z, i, result[1]);
                }
            }
          else
            {
              for (i = idx1; i <= idx2; ++i)
                {
                  if (i - H > 0)
                    {
                      /* delete oldest window sample as we move closer to edge */
                      (accum->delete_oldest)(w->state);
                    }

                  /* yi = acc_get [ work(i:K-2) ] */
                  (accum->get)(accum_params, result, w->state);
                  gsl_vector_set(y, i, result[0]);

                  if (z != NULL)
                    gsl_vector_set(z, i, result[1]);
                }
            }
        }
      else
        {
          /* pad final windows and fill y(n-J:n-1) */
          for (i = 0; i < J; ++i)
            {
              int idx = n - J + i;

              (accum->insert)(xN, w->state);

              if (idx >= 0)
                {
                  (accum->get)(accum_params, result, w->state);
                  gsl_vector_set(y, idx, result[0]);

                  if (z != NULL)
                    gsl_vector_set(z, idx, result[1]);
                }
            }
        }

      return GSL_SUCCESS;
    }
}

/*
gsl_movstat_apply()
  Apply user-defined moving window function to input vector

Inputs: endtype - end point handling criteria
        F       - user-defined function
        x       - input vector, size n
        y       - output vector, size n
        w       - workspace
*/

int
gsl_movstat_apply(const gsl_movstat_end_t endtype, const gsl_movstat_function * F,
                  const gsl_vector * x, gsl_vector * y, gsl_movstat_workspace * w)
{
  int status = gsl_movstat_apply_accum(endtype, x, gsl_movstat_accum_userfunc, (void *) F, y, NULL, w);
  return status;
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
#ifndef EIGEN_BENCH_UTIL_H
#define EIGEN_BENCH_UTIL_H

#include <Eigen/Core>
#include "BenchTimer.h"

using namespace std;
using namespace Eigen;

#include <boost/preprocessor/repetition/enum_params.hpp>
#include <boost/preprocessor/repetition.hpp>
#include <boost/preprocessor/seq.hpp>
#include <boost/preprocessor/array.hpp>
#include <boost/preprocessor/arithmetic.hpp>
#include <boost/preprocessor/comparison.hpp>
#include <boost/preprocessor/punctuation.hpp>
#include <boost/preprocessor/punctuation/comma.hpp>
#include <boost/preprocessor/stringize.hpp>

template<typename MatrixType> void initMatrix_random(MatrixType& mat) __attribute__((noinline));
template<typename MatrixType> void initMatrix_random(MatrixType& mat)
{
  mat.setRandom();// = MatrixType::random(mat.rows(), mat.cols());
}

template<typename MatrixType> void initMatrix_identity(MatrixType& mat) __attribute__((noinline));
template<typename MatrixType> void initMatrix_identity(MatrixType& mat)
{
  mat.setIdentity();
}

#ifndef __INTEL_COMPILER
#define DISABLE_SSE_EXCEPTIONS()  { \
  int aux; \
  asm( \
  "stmxcsr   %[aux]           \n\t" \
  "orl       $32832, %[aux]   \n\t" \
  "ldmxcsr   %[aux]           \n\t" \
  : : [aux] "m" (aux)); \
}
#else
#define DISABLE_SSE_EXCEPTIONS()
#endif

#ifdef BENCH_GMM
#include <gmm/gmm.h>
template <typename EigenMatrixType, typename GmmMatrixType>
void eiToGmm(const EigenMatrixType& src, GmmMatrixType& dst)
{
  dst.resize(src.rows(),src.cols());
  for (int j=0; j<src.cols(); ++j)
    for (int i=0; i<src.rows(); ++i)
      dst(i,j) = src.coeff(i,j);
}
#endif

#ifdef BENCH_GSL
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_eigen.h>
template <typename EigenMatrixType>
void eiToGsl(const EigenMatrixType& src, gsl_matrix** dst)
{
  for (int j=0; j<src.cols(); ++j)
    for (int i=0; i<src.rows(); ++i)
      gsl_matrix_set(*dst, i, j, src.coeff(i,j));
}
#endif

#ifdef BENCH_UBLAS
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/vector.hpp>
template <typename EigenMatrixType, typename UblasMatrixType>
void eiToUblas(const EigenMatrixType& src, UblasMatrixType& dst)
{
  dst.resize(src.rows(),src.cols());
  for (int j=0; j<src.cols(); ++j)
    for (int i=0; i<src.rows(); ++i)
      dst(i,j) = src.coeff(i,j);
}
template <typename EigenType, typename UblasType>
void eiToUblasVec(const EigenType& src, UblasType& dst)
{
  dst.resize(src.size());
  for (int j=0; j<src.size(); ++j)
      dst[j] = src.coeff(j);
}
#endif

#endif // EIGEN_BENCH_UTIL_H


text:
#ifndef ARRAY_H
#define ARRAY_H

#include <cstdarg>
#include <cstdio>
#include <cstdlib>
#include <cassert>
#include <iostream>
#include <algorithm>
#include <charm++.h>
#if CMK_HAS_CBLAS
#include <cblas.h>
#endif

namespace CharjArray {
  class Range {
  public:
    int size, start, stop;
    Range() {}
    Range(int size_) : size(size_), start(0), stop(size) {}
    Range(int start_, int stop_) :
    size(stop_ - start_), start(start_), stop(stop_) {
      assert(stop >= start);
    }
    void pup(PUP::er& p) {
        p | size;
        p | start;
        p | stop;
    }
  };

  template<int dims>
  class Domain {
  public:
    Range ranges[dims];

    Domain() {}

    Domain(Range ranges_[]) {
      for (int i = 0; i < dims; i++)
	ranges[i] = ranges_[i];
    }

    Domain(Range range) {
      ranges[0] = range;
    }

    Domain(Range range1, Range range2) {
      // TODO: fix Charj generator so it uses the array
      ranges[0] = range1;
      ranges[1] = range2;
    }

    int size() const {
      int total = 0;

      for (int i = 0; i < dims; i++)
	if (total == 0)
	  total = ranges[i].size;
	else
	  total *= ranges[i].size;

      return total;
    }

    void pup(PUP::er& p) {
        for (int i=0; i<dims; ++i) p | ranges[i];
    }
  };

  template<int dims>
  class RowMajor {
  public:
    static int access(const int i, const Domain<dims> &d) {
      return i - d.ranges[0].start;
    }
    static int access(const int i, const int j, const Domain<dims> &d) {
      return (i - d.ranges[0].start) * d.ranges[1].size + j -
        d.ranges[1].start;
    }
    // Generic access method, not used right now.
    // static int access(const int *i, const Domain<dims> &d) {
    //   int off = i[0];
    //   int dimoff = 1;
    //   for (int j = ndims-1; j > 0; --j) {
    //     dimoff *= d.ranges[j].size;
    //     off += dimoff * (i[j] - d.ranges[j].start);
    //   }
    //   return off;
    // }
  };

  template<int dims>
  class ColMajor {
  public:
    static int access(const int i, const Domain<dims> &d) {
      return i - d.ranges[0].start;
    }
    static int access(const int i, const int j, const Domain<dims> &d) {
      return (j - d.ranges[1].start) * d.ranges[1].size + i -
        d.ranges[0].start;
    }
  };

  template<class type, int dims = 1, class atype = RowMajor<dims> >
  class Array {
  private:
    Domain<dims> domain;
    type *block;
    int ref_cout;
    bool did_init;
    Array* ref_parent;

  public:
    Array(Domain<dims> domain_) : ref_parent(0), did_init(false) {
      init(domain_);
    }

    Array(type **block_) : did_init(false) {
      block = *block_;
    }

    Array(type& block_) : did_init(false) {
      block = &block_;
    }

    Array() : ref_parent(0), did_init(false) {

    }

    Array(Array* parent, Domain<dims> domain_)
        : ref_parent(parent), did_init(false) {
      domain = domain_;
      block = parent->block;
    }

    void init(Domain<dims> &domain_) {
      domain = domain_;
      //if (atype == ROW_MAJOR)
      block = new type[domain.size()];
      //printf("Array: allocating memory, size=%d, base pointer=%p\n",
      //       domain.size(), block);
      did_init = true;
    }

    type* raw() { return block; }

    ~Array() {
      if (did_init) delete[] block;
    }

    /*type* operator[] (const Domain<dims> &domain) {
      return block[domain.ranges[0].size];
      }*/

    type& operator[] (const int i) {
      return block[atype::access(i, domain)];
    }

    const type& operator[] (const int i) const {
      return block[atype::access(i, domain)];
    }

    type& access(const int i) {
      return this->operator[](i);
    }

    type& access(const int i, const int j) {
      //printf("Array: accessing, index (%d,%d), offset=%d, base pointer=%p\n",
      //i, j, atype::access(i, j, domain), block);
      return block[atype::access(i, j, domain)];
    }

    type& access(const int i, const Range r) {
      Domain<1> d(r);
      //printf("Array: accessing subrange, size = %d, range (%d,%d), base pointer=%p\n",
      //d.size(), r.start, r.stop, block);
      type* buf = new type[d.size()];
      for (int j = 0; j < d.size(); j++) {
        //printf("Array: copying element (%d,%d), base pointer=%p\n", i, j, block);
        buf[j] = block[atype::access(i, j, domain)];
      }
      return *buf;
    }

    const type& access(const int i, const int j) const {
      return block[atype::access(i, j, domain)];
    }

    Array<type, dims, atype>* operator[] (const Domain<dims> &domain) {
      return new Array<type, dims, atype>(this, domain);
    }

    int size() const {
      return domain.size();
    }

    int size(int dim) const {
      return domain.ranges[dim].size;
    }

    void pup(PUP::er& p) {
        p | domain;
        if (p.isUnpacking()) {
            block = new type[domain.size()];
        }
        PUParray(p, block, domain.size());
    }

    void fill(const type &t) {
      for (int i = 0; i < domain.size(); ++i)
	block[i] = t;
    }

    /// Do these arrays have the same shape and contents?
    bool operator==(const Array &rhs) const
    {
      for (int i = 0; i < dims; ++i)
	if (this->size(i) != rhs.size(i))
	  return false;

      for (int i = 0; i < this->size(); ++i)
	if (this->block[i] != rhs.block[i])
	  return false;

      return true;
    }
    bool operator!=(const Array &rhs) const
    {
      return !(*this == rhs);
    }
  };

  /**
     A local Matrix class for various sorts of linear-algebra work.

     Indexed from 0, to reflect the C-heritage of Charj.
   */
  template <typename V, class atype = RowMajor<2> >
  class Matrix : public Array<V, 2, atype>
  {
  public:
    Matrix() { }
    /// A square matrix
    Matrix(unsigned int n) : Array<V,2,atype>(Domain<2>(n,n)) { }

    /// A identity matrix
    static Matrix* ident(int n)
    {
      Matrix *ret = new Matrix(n);
      ret->fill(0);

      for (int i = 0; i < n; ++i)
	ret->access(i,i) = 1;

      return ret;
    }
  };

  template <typename T, class atype = RowMajor<1> >
  class Vector : public Array<T, 1, atype>
  {
  public:
    Vector() { }
    Vector(unsigned int n) : Array<T, 1, atype>(Range(n)) { }
  };

  /// Compute the inner (dot) product v1^T * v2
  // To compute v1^H * v2, call as dot(v1.C(), v2)
  template<typename T, class atype1, class atype2>
  T dot(const Vector<T, atype1> *pv1, const Vector<T, atype2> *pv2)
  {
    const Vector<T, atype1> &v1 = *pv1, &v2 = *pv2;
    assert(v1.size() == v2.size());
    // XXX: This default initialization worries me some, since it
    // won't necessarily be an additive identity for all T. - Phil
    T ret = T();
    int size = v1.size();
    for (int i = 0; i < size; ++i)
      ret += v1[i] * v2[i];
    return ret;
  }
#if CMK_HAS_CBLAS
  template <>
  float dot<float, RowMajor<1>, RowMajor<1> >(const Vector<float, RowMajor<1> > *pv1,
					      const Vector<float, RowMajor<1> > *pv2)
  {
    const Vector<float, RowMajor<1> > &v1 = *pv1, &v2 = *pv2;
    assert(v1.size() == v2.size());
    return cblas_sdot(v1.size(), &(v1[0]), 1, &(v2[0]), 1);
  }
  template <>
  double dot<double, RowMajor<1>, RowMajor<1> >(const Vector<double, RowMajor<1> > *pv1,
						const Vector<double, RowMajor<1> > *pv2)
  {
    const Vector<double, RowMajor<1> > &v1 = *pv1, &v2 = *pv2;
    assert(v1.size() == v2.size());
    return cblas_ddot(v1.size(), &(v1[0]), 1, &(v2[0]), 1);
  }
#endif

  /// Computer the 1-norm of the given vector
  template<typename T, class atype>
    T norm1(const Vector<T, atype> *pv)
  {
    const Vector<T, atype> &v = *pv;
    // XXX: See comment about additive identity in dot(), above
    T ret = T();
    int size = v.size();
    for (int i = 0; i < size; ++i)
      ret += v[i];
    return ret;
  }

  /// Compute the Euclidean (2) norm of the given vector
  template<typename T, class atype>
    T norm2(const Vector<T, atype> *pv)
  {
    const Vector<T, atype> &v = *pv;
    // XXX: See comment about additive identity in dot(), above
    T ret = T();
    int size = v.size();
    for (int i = 0; i < size; ++i)
      ret += v[i] * v[i];
    return sqrt(ret);
  }
#if CMK_HAS_CBLAS
  template<>
    float norm2<float, RowMajor<1> >(const Vector<float, RowMajor<1> > *pv)
  {
    const Vector<float, RowMajor<1> > &v = *pv;
    return cblas_snrm2(v.size(), &(v[0]), 1);
  }
  template<>
    double norm2<double, RowMajor<1> >(const Vector<double, RowMajor<1> > *pv)
  {
    const Vector<double, RowMajor<1> > &v = *pv;
    return cblas_dnrm2(v.size(), &(v[0]), 1);
  }
#endif

  /// Compute the infinity (max) norm of the given vector
  // Will fail on zero-length vectors
  template<typename T, class atype>
    T normI(const Vector<T, atype> *pv)
  {
    const Vector<T, atype> &v = *pv;
    T ret = v[0];
    int size = v.size();
    for (int i = 1; i < size; ++i)
      ret = max(ret, v[i]);
    return ret;
  }

  /// Scale a vector by some constant
  template<typename T, typename U, class atype>
    void scale(const T &t, Vector<U, atype> *pv)
  {
    const Vector<T, atype> &v = *pv;
    int size = v.size();
    for (int i = 0; i < size; ++i)
      v[i] = t * v[i];
  }
#if CMK_HAS_CBLAS
  template<>
    void scale<float, float, RowMajor<1> >(const float &t,
					   Vector<float, RowMajor<1> > *pv)
  {
    Vector<float, RowMajor<1> > &v = *pv;
    cblas_sscal(v.size(), t, &(v[0]), 1);
  }
  template<>
    void scale<double, double, RowMajor<1> >(const double &t,
					     Vector<double, RowMajor<1> > *pv)
  {
    Vector<double, RowMajor<1> > &v = *pv;
    cblas_dscal(v.size(), t, &(v[0]), 1);
  }
#endif

  /// Add one vector to a scaled version of another
  template<typename T, typename U, class atype>
    void axpy(const T &a, const Vector<U, atype> *px, Vector<U, atype> *py)
  {
    Vector<T, atype> &x = *px;
    const Vector<T, atype> &y = *py;
    int size = x.size();
    assert(size == y.size());
    for (int i = 0; i < size; ++i)
      x[i] = a * x[i] + y[i];
  }
#if CMK_HAS_CBLAS
  template<>
    void axpy<float, float, RowMajor<1> >(const float &a,
					  const Vector<float, RowMajor<1> > *px,
					  Vector<float, RowMajor<1> > *py)
  {
    const Vector<float, RowMajor<1> > &x = *px;
    Vector<float, RowMajor<1> > &y = *py;
    int size = x.size();
    assert(size == y.size());
    cblas_saxpy(size, a, &(x[0]), 1, &(y[0]), 1);
  }
  template<>
    void axpy<double, double, RowMajor<1> >(const double &a,
					  const Vector<double, RowMajor<1> > *px,
					  Vector<double, RowMajor<1> > *py)
  {
    const Vector<double, RowMajor<1> > &x = *px;
    Vector<double, RowMajor<1> > &y = *py;
    int size = x.size();
    assert(size == y.size());
    cblas_daxpy(size, a, &(x[0]), 1, &(y[0]), 1);
  }
#endif
}

#endif

max_stars_repo_name: scottkwarren/config-db

max_issues_repo_name: scottkwarren/config-db

max_forks_repo_name: scottkwarren/config-db


text:
#pragma once
#include "iparam.h"
#include "configaccess.h"
#include "optioninfo.h"
#include "format.h"
#include "streamreader.h"
#include <gsl/gsl>
#include <cmdlime/errors.h>
#include <cmdlime/customnames.h>
#include <sstream>
#include <optional>
#include <memory>
#include <functional>

namespace cmdlime::detail{

template<typename T>
class Param : public IParam{
public:
    Param(std::string name,
          std::string shortName,
          std::string type,
          std::function<T&()> paramGetter)
        : info_(std::move(name), std::move(shortName), std::move(type))
        , paramGetter_(std::move(paramGetter))
    {
    }

    void setDefaultValue(const T& value)
    {
        hasValue_ = true;
        defaultValue_ = value;
    }

    OptionInfo& info() override
    {
        return info_;
    }

    const OptionInfo& info() const override
    {
        return info_;
    }

private:
    bool read(const std::string& data) override
    {
        auto stream = std::stringstream{data};
        if (!readFromStream(stream, paramGetter_()))
            return false;
        hasValue_ = true;
        return true;
    }

    bool hasValue() const override
    {
        return hasValue_;
    }

    bool isOptional() const override
    {
        return defaultValue_.has_value();
    }

    std::string defaultValue() const override
    {
        if (!defaultValue_.has_value())
            return {};
        auto stream = std::stringstream{};
        stream << defaultValue_.value();
        return stream.str();
    }

private:
    OptionInfo info_;
    std::function<T&()> paramGetter_;
    std::optional<T> defaultValue_;
    bool hasValue_ = false;
};

template <>
inline bool Param<std::string>::read(const std::string& data)
{
    paramGetter_() = data;
    hasValue_ = true;
    return true;
}

template<typename T, typename TConfig>
class ParamCreator{
    using NameProvider = typename Format<ConfigAccess<TConfig>::format()>::nameProvider;
public:
    ParamCreator(TConfig& cfg,
                 const std::string& varName,
                 const std::string& type,
                 std::function<T&()> paramGetter)
        : cfg_(cfg)
    {
        Expects(!varName.empty());
        Expects(!type.empty());
        param_ = std::make_unique<Param<T>>(NameProvider::name(varName),
                                            NameProvider::shortName(varName),
                                            NameProvider::valueName(type),
                                            std::move(paramGetter));
    }

    ParamCreator<T, TConfig>& operator<<(const std::string& info)
    {
        param_->info().addDescription(info);
        return *this;
    }

    ParamCreator<T, TConfig>& operator<<(const Name& customName)
    {
        param_->info().resetName(customName.value());
        return *this;
    }

    ParamCreator<T, TConfig>& operator<<(const ShortName& customName)
    {
        static_assert(Format<ConfigAccess<TConfig>::format()>::shortNamesEnabled,
                      "Current command line format doesn't support short names");
        param_->info().resetShortName(customName.value());
        return *this;
    }

    ParamCreator<T, TConfig>& operator<<(const WithoutShortName&)
    {
        static_assert(Format<ConfigAccess<TConfig>::format()>::shortNamesEnabled,
                      "Current command line format doesn't support short names");
        param_->info().resetShortName({});
        return *this;
    }

    ParamCreator<T, TConfig>& operator<<(const ValueName& valueName)
    {
        param_->info().resetValueName(valueName.value());
        return *this;
    }

    ParamCreator<T, TConfig>& operator()(T defaultValue = {})
    {
        defaultValue_ = std::move(defaultValue);
        param_->setDefaultValue(defaultValue_);
        return *this;
    }

    operator T()
    {
        ConfigAccess<TConfig>{cfg_}.addParam(std::move(param_));
        return defaultValue_;
    }

private:
    std::unique_ptr<Param<T>> param_;
    T defaultValue_;
    TConfig& cfg_;
};

template <typename T, typename TConfig>
ParamCreator<T, TConfig> makeParamCreator(TConfig& cfg,
                                          const std::string& varName,
                                          const std::string& type,
                                          std::function<T&()> paramGetter)
{
    return ParamCreator<T, TConfig>{cfg, varName, type, std::move(paramGetter)};
}

}

max_stars_repo_name: GerHobbelt/hypertextcpp

max_issues_repo_name: GerHobbelt/hypertextcpp

max_forks_repo_name: GerHobbelt/hypertextcpp