text:
/*! \file interp.h
    \brief Functions for interpolation to initialize the root solver
 */

// used by the gsl interpolation functions
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_sort_double.h>

#ifndef REAL_TYPEDEF
#define REAL_TYPEDEF
#ifndef single			// compiler option determines variable size
typedef double REAL;
#else
typedef float REAL;
#endif
#endif

#define N_INTERP_POINTS 3
#define NEAREST_TOL 1.0e-3

#define MIN_IDX(indices,values) ((values[0] < values[1]) ? indices[0]:indices[1])
#define DISTANCE(x1,y1,x2,y2) hypot(x1-x2,y1-y2)
#define MAXIMUM(x,y) ((x) > (y) ? (x) : (y))

// interpolation functions
REAL interpBetw2ptsDist(REAL distx_x0,REAL f0,REAL f1,REAL distx1_x0);
int cartesianToCylindrical(int npts,REAL R[],REAL Z[],REAL zeta[],REAL x[],REAL y[],REAL z[]);
int interp2d(int nx,int nxi,REAL x[],REAL y[],REAL z[],REAL w[],REAL xi[],REAL yi[],REAL zi[],REAL wi[]);
int find_nearestn(int nclosest,int nx,int nxi,size_t idx_closest[nx][nclosest],REAL x[],REAL y[],REAL xi[],REAL yi[]);
int calc_slopes(double mx[],double my[],size_t idxs[],double x[],double y[],double z[]);

// legacy functions used by previous versions of the root finder
int interp(REAL y[],REAL xi[],REAL si[],int ny,int nxi);
int interpxy(REAL x[],REAL y[],REAL xi[],REAL yi[],int nx,int nxi);
int interpgsl(REAL x[],REAL y[],REAL xi[],REAL yi[],int nx,int nxi);
int find_nearest(int nx,int xirows,int xicols,int xi_idx_below[],int yi_idx_below[],REAL x[],REAL y[],REAL xi[][xicols],REAL yi[][xicols]);
int interpReg2d(int nx,int xirows,int xicols,REAL x[],REAL y[],REAL z[],REAL w[],REAL xi[][xicols],REAL yi[][xicols],REAL zi[][xicols],REAL wi[][xicols]);

// ---  local function definitions

//! Linear interpolation using two points to a third point
/*! Uses the distance between the two known points ``distx1_x0'' = (x1-x0),
    the distance between the unknown point to a known point ``distx_x0'' = (x-x0),
    and the values of the function at x1 and x0 (f1 = f(x1),f2 = f(x2)),
    to find the value at x: f= f(x).
 */
REAL interpBetw2ptsDist(REAL distx_x0,REAL f0,REAL f1,REAL distx1_x0){
  REAL f = (f1-f0)/distx1_x0*distx_x0 + f0;
  return f;
}

//! Converts cartesian coordinates (x,y,z) to cylindrical coordinates (R,Z,zeta)
int cartesianToCylindrical(int npts,REAL R[],REAL Z[],REAL zeta[],REAL x[],REAL y[],REAL z[]){
  int i;
  for(i=0;i<npts;i++){
    R[i]=hypot(x[i],z[i]);
    Z[i]=y[i];
    zeta[i]=-atan2(z[i],x[i]);	// negative because in forming right-hand coord system
  }                               //  x: radially out, y: vertically up, z: toroidal
                                  //  esi routines give theta going clockwise,
                                  //  a: out in minor radius, gives z and theta antiparallel

  return 0;
}

//! Finds the nearest n points on a grid of (xi,yi) for each point (x,y)
int find_nearestn(int nclosest,int nx,int nxi,size_t idx_closest[nx][nclosest],REAL x[],REAL y[],REAL xi[],REAL yi[]){
  double dist[nxi],d01,d02,d12,hyp;
  int point,i,n_temp;
  size_t *closest_pts;
  size_t *idx_closest_pt;

#if DEBUG
  fprintf(stderr,"point,   d01,    d02,   d12,    hyp,    collinear\n");
#endif

  for(point=0;point<nx;point++){		// loop for each x,y

    for(i=0;i<nxi;i++){	// find the distance between (x,y) and all the (xi,yi)
      //      dist[i]=sqrt((x[point]-xi[i])*(x[point]-xi[i])+(y[point]-yi[i])*(y[point]-yi[i]));
      dist[i]=DISTANCE(x[point],y[point],xi[i],yi[i]);
    }
    idx_closest_pt=idx_closest[point];
    gsl_sort_smallest_index(idx_closest_pt,nclosest,dist,1,nxi);// find the nearest xi

    // make sure they aren't collinear
    d01=DISTANCE(xi[idx_closest_pt[0]],yi[idx_closest_pt[0]],xi[idx_closest_pt[1]],yi[idx_closest_pt[1]]);
    d02=DISTANCE(xi[idx_closest_pt[0]],yi[idx_closest_pt[0]],xi[idx_closest_pt[2]],yi[idx_closest_pt[2]]);
    d12=DISTANCE(xi[idx_closest_pt[1]],yi[idx_closest_pt[1]],xi[idx_closest_pt[2]],yi[idx_closest_pt[2]]);
    //    hyp = fmax(fmax(d01,d02),d12);
    hyp = MAXIMUM(MAXIMUM(d01,d02),d12);
    n_temp = nclosest;
    while((d01+d02+d12-2.0*hyp < NEAREST_TOL) && (nclosest < nxi)){ // loop until the points aren't collinear
      n_temp++;
      closest_pts = (size_t *)malloc(sizeof(size_t)*n_temp);
      gsl_sort_smallest_index(closest_pts,n_temp,dist,1,nxi);// find the nearest xi
      idx_closest_pt[2]=closest_pts[n_temp-1];
      d02=DISTANCE(xi[idx_closest_pt[0]],yi[idx_closest_pt[0]],xi[idx_closest_pt[2]],yi[idx_closest_pt[2]]);
      d12=DISTANCE(xi[idx_closest_pt[1]],yi[idx_closest_pt[1]],xi[idx_closest_pt[2]],yi[idx_closest_pt[2]]);
      //hyp = fmax(fmax(d01,d02),d12);
      hyp = MAXIMUM(MAXIMUM(d01,d02),d12);

      free(closest_pts);
    }
#if DEBUG
    fprintf(stderr,"%d, %g, %g, %g, %g, %g\n",point,d01,d02,d12,hyp,2.0*hyp-d01-d02-d12);
#endif

  }
    //   for(i=0;i<xirows*xicols;i++){
    //     fprintf(stderr,"xi       yi\n");
    //   fprintf(stderr,"%g, %g\n",xi[0][i],yi[0][i]);
    //   }
    //   for(i=0;i<nx;i++){
    //     fprintf(stderr,"x,      xi_below,      y,    ,yi_below\n");
    //     fprintf(stderr,"%g, %g, %g, %g\n",x[i],xi[xi_idx_below[i]][0],y[i],yi[0][yi_idx_below[i]]);
    //   }

    return 0;
 }

  //! Calculates the slopes (mx,my) for three points (x,y) which each have function value z
  /*! these three points are have the indices idxs */
  int calc_slopes(double mx[],double my[],size_t idxs[],double x[],double y[],double z[]){
    size_t i,ip1,ip2,j,jp1,jp2;
    double dx[N_INTERP_POINTS],dy[N_INTERP_POINTS],dz[N_INTERP_POINTS];

    for(j=0;j<N_INTERP_POINTS;j++){
      ip1 = idxs[(size_t)fmod(j+1,N_INTERP_POINTS)];
      ip2 = idxs[(size_t)fmod(j+2,N_INTERP_POINTS)];
      dx[j] = x[ip1]-x[ip2];
      dy[j] = y[ip1]-y[ip2];
      dz[j] = z[ip1]-z[ip2];
    }
    for(j=0;j<N_INTERP_POINTS;j++){
      jp1 = (size_t)fmod(j+1,N_INTERP_POINTS);
      jp2 = (size_t)fmod(j+2,N_INTERP_POINTS);
      i = idxs[j];
      ip1 = idxs[jp1];
      ip2 = idxs[jp2];
      mx[j] = -1.0*(y[i]*dz[j]+y[ip1]*dz[jp1]+y[ip2]*dz[jp2]);
      mx[j]/=(x[i]*dy[j]+x[ip1]*dy[jp1]+x[ip2]*dy[jp2]);
      my[j] = -1.0*(x[i]*dz[j]+x[ip1]*dz[jp1]+x[ip2]*dz[jp2]);
      my[j]/=(y[i]*dx[j]+y[ip1]*dx[jp1]+y[ip2]*dx[jp2]);
    }

    return 0;
  }

//! High-level 2d interpolation function
/*! Given values zi = f(xi,yi) and wi = g(xi,yi) on an arbitrary, irregular grid (xi,yi)
    finds the values at (x,y).
    The method is to find the nearest 3 grid points for each (x,y), then linearly interpolate (triangulation).
 */
  int interp2d(int nx,int nxi,REAL x[],REAL y[],REAL z[],REAL w[],REAL xi[],REAL yi[],REAL zi[],REAL wi[]){
    int point;
    size_t *idxs;
    double mxi[N_INTERP_POINTS],myi[N_INTERP_POINTS],a,b,c;

    // find the indices of the 3 (xi,yi) closest to x,y
    size_t idx_closest[nx][N_INTERP_POINTS];
    find_nearestn(N_INTERP_POINTS,nx,nxi,idx_closest,x,y,xi,yi);

    // loop around points
    for(point=0;point<nx;point++){
      idxs = idx_closest[point];
      // first get z
      // calculate the slopes (mx,my) for the nearest 3 points
      calc_slopes(mxi,myi,idxs,xi,yi,zi);
      // linear interp at each nearest point
      a = zi[idxs[0]] + mxi[0]*(x[point]-xi[idxs[0]]) + myi[0]*(y[point]-yi[idxs[0]]);
      b = zi[idxs[1]] + mxi[1]*(x[point]-xi[idxs[1]]) + myi[1]*(y[point]-yi[idxs[1]]);
      c = zi[idxs[2]] + mxi[2]*(x[point]-xi[idxs[2]]) + myi[2]*(y[point]-yi[idxs[2]]);
      // average to get the best guess
      z[point] = (a+b+c)/3.0;

      // next get w
      // calculate the slopes (mx,my) for the nearest 3 points
      calc_slopes(mxi,myi,idxs,xi,yi,wi);
      // linear interp at each nearest point
      a = wi[idxs[0]] + mxi[0]*(x[point]-xi[idxs[0]]) + myi[0]*(y[point]-yi[idxs[0]]);
      b = wi[idxs[1]] + mxi[1]*(x[point]-xi[idxs[1]]) + myi[1]*(y[point]-yi[idxs[1]]);
      c = wi[idxs[2]] + mxi[2]*(x[point]-xi[idxs[2]]) + myi[2]*(y[point]-yi[idxs[2]]);
      // average to get the best guess
      w[point] = (a+b+c)/3.0;

    }

    return 0;
  }

// --- legacy functions

int interp(REAL y[],REAL xi[],REAL si[],int ny,int nxi){
  int i,j;
  //  printf("interp:\n%10s,%10s,%10s,%10s,%10s\n","q_orig","q_final","q_i","q_i+1","s_i");
  for(i=0;i<ny;i++){
    for(j=0;j<nxi-1;j++) if(xi[j+1] > y[i]) break;
    //printf("%10g,",y[i]);
    y[i] = xi[j] + si[j]*(y[i]-xi[j]);
    //printf("%10g,%10g,%10g,%10g\n",y[i],xi[j],xi[j+1],si[j]);
  }
  return 0;
}

int interpxy(REAL x[],REAL y[],REAL xi[],REAL yi[],int nx,int nxi){
  int i,j;
  REAL si[nxi-1];
  // use linear interpolation
  //printf("\nstheta:");
  for(i=0;i<nxi-1;i++){
    si[i] = (yi[i+1]-yi[i])/(xi[i+1]-xi[i]);
    //printf("%g,",si[i]);
  }
  si[nxi-2] = si[0];
  //printf("\n");

  //printf("interp:\n%10s,%10s,%10s,%10s,%10s\n","q_orig","q_final","q_i","q_i+1","s_i");
  for(i=0;i<nx;i++){
    for(j=0;j<nxi-1;j++) if(xi[j+1] > x[i]) break;
    //printf("%10g,",x[i]);
    y[i] = yi[j] + si[j]*(x[i]-xi[j]);
    //printf("%10g,%10g,%10g,%10g\n",y[i],xi[j],xi[j+1],si[j]);
  }

  return 0;
}

// interpolate the array x[] of size nx
// on the grid xi[],y[i] of size nxi
int interpgsl(REAL x[],REAL y[],REAL xi[],REAL yi[],int nx,int nxi){
  // make sure periodic
  //  yi[nxi-1]=yi[0];
  // xi[nxi-1]=-xi[0];
  gsl_interp_accel *acc = gsl_interp_accel_alloc ();
  const gsl_interp_type *t = gsl_interp_cspline;
  gsl_spline *spline = gsl_spline_alloc (t, nxi);

  int i;
  gsl_spline_init (spline, xi, yi, nxi);
  //printf("\ninterpgsl:\n");
  for (i = 0; i < nx; i++)
    {
      y[i] = gsl_spline_eval (spline, x[i], acc);
      //printf ("%g, %g\n", x[i], y[i]);
    }

  gsl_spline_free (spline);
  gsl_interp_accel_free (acc);

  return 0;
}

int find_nearest1d(int nx,int nxi,int nnearest,int stride,size_t xi_idx[nx][nnearest],REAL x[],REAL xi[]){
  int i,j;
  REAL xerr[nxi];
  for(i=0;i<nx;i++){
    for(j=0;j<nxi;j++){
      xerr[j]=fabs(x[i]-xi[j*stride]);
    }
    gsl_sort_smallest_index(xi_idx[i],nnearest,xerr,stride,nxi);// find the n nearest xi
  }

  return 0;
}

// xi varies slowly (each row), yi varies quickly (each col)
int find_nearest(int nx,int xirows,int xicols,int xi_idx_below[],int yi_idx_below[],REAL x[],REAL y[],REAL xi[][xicols],REAL yi[][xicols]){
  size_t nearest2_idxs[2];
  double nearest2_vals[2],xerr[xirows],yerr[xirows];
  int i,j;

  for(i=0;i<nx;i++){		// loop for each x,y

    for(j=0;j<xirows;j++){	// find the error between x and all the xi
      xerr[j]=fabs(x[i]-xi[j][0]);
    }
    gsl_sort_smallest_index(nearest2_idxs,2,xerr,1,xirows);// find the two nearest xi
    nearest2_vals[0] = xi[nearest2_idxs[0]][0];
    nearest2_vals[1] = xi[nearest2_idxs[1]][0];
    xi_idx_below[i] = MIN_IDX(nearest2_idxs,nearest2_vals); // this is the nearest xi less than x

    for(j=0;j<xicols;j++){	// find the error between y and all the yi
      yerr[j]=fabs(y[i]-yi[0][j]);
    }
    gsl_sort_smallest_index(nearest2_idxs,2,yerr,1,xicols);// find the two nearest yi
    nearest2_vals[0] = yi[0][nearest2_idxs[0]];
    nearest2_vals[1] = yi[0][nearest2_idxs[1]];
    yi_idx_below[i] = MIN_IDX(nearest2_idxs,nearest2_vals); // this is the nearest yi less than y
  }

  return 0;
}

// limitation: assumes regularly-spaced grids xi and yi
int interpReg2d(int nx,int xirows,int xicols,REAL x[],REAL y[],REAL z[],REAL w[],REAL xi[][xicols],REAL yi[][xicols],REAL zi[][xicols],REAL wi[][xicols]){
  int i,j;

  // form the matrix of of slopes for y interpolation
  double slopeziyi[xirows][xicols-1],slopewiyi[xirows][xicols-1];
  double invdi;
  for(i=0;i<xirows;i++){
    for(j=0;j<xicols-1;j++){
      invdi = 1.0/(yi[i][j+1]-yi[i][j]);
      slopeziyi[i][j] = (zi[i][j+1]-zi[i][j])*invdi;
      slopewiyi[i][j] = (wi[i][j+1]-wi[i][j])*invdi;
    }
  }

  // find the indices of the nearest points (xi[i],yi[j]); note always exists (xi[i+1],yi[j+1])
  int xi_idx_below[nx],yi_idx_below[nx];
  find_nearest(nx,xirows,xicols,xi_idx_below,yi_idx_below,x,y,xi,yi);

  // for fixed x, interpolate in y
  double a,b,dist_y_yi,x_frac;
  for(i=0;i<nx;i++){
    dist_y_yi=y[i]-yi[0][yi_idx_below[i]];
    x_frac = (x[i]-xi[xi_idx_below[i]][0])/(xi[xi_idx_below[i]+1][0]-xi[xi_idx_below[i]][0]);

    // compute z[i]
    // for fixed xi = i, interpolate in y between yi[j] and yi[j+1]
    a = slopeziyi[xi_idx_below[i]][yi_idx_below[i]]*dist_y_yi + zi[xi_idx_below[i]][yi_idx_below[i]];
    // for fixed xi = i+1, interpolate in y between yi[j] and yi[j+1]
    b = slopeziyi[xi_idx_below[i]+1][yi_idx_below[i]]*dist_y_yi + zi[xi_idx_below[i]+1][yi_idx_below[i]];
    // now interpolate in x between xi[i] and xi[i+1]
    z[i] = (b-a)*x_frac + a;

    // compute w[i]
    // for fixed xi = i, interpolate in y between yi[j] and yi[j+1]
    a = slopewiyi[xi_idx_below[i]][yi_idx_below[i]]*dist_y_yi + wi[xi_idx_below[i]][yi_idx_below[i]];
    // for fixed xi = i+1, interpolate in y between yi[j] and yi[j+1]
    b = slopewiyi[xi_idx_below[i]+1][yi_idx_below[i]]*dist_y_yi + wi[xi_idx_below[i]+1][yi_idx_below[i]];
    // now interpolate in x between xi[i] and xi[i+1]
    w[i] = (b-a)*x_frac + a;

  }

  return 0;
}

max_stars_repo_name: LeiShi/Synthetic-Diagnostics-Platform


text:
#ifndef UTIL_H_
#define UTIL_H_

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

/*
 * constants
 */

#define TYPE_REQ 1
#define TYPE_REQ_FOLLOW 2
#define TYPE_RES 0

#define NB_MBUF 8191
#define MBUF_SIZE (2048 + sizeof(struct rte_mbuf) + RTE_PKTMBUF_HEADROOM)
#define MBUF_CACHE_SIZE 32
#define MAX_BURST_SIZE 32
#define MAX_LCORES 32
#define NB_RXD 128  // RX descriptors
#define NB_TXD 512  // TX descriptors
#define RX_QUEUE_PER_LCORE 4

#define IP_SRC "10.1.0.12"
#define IP_DST "10.1.0.1"
#define CLIENT_PORT 11234
#define SERVICE_PORT 1234

static const struct rte_eth_conf port_conf = {
    .rxmode =
        {
            .max_rx_pkt_len = ETHER_MAX_LEN,
            .split_hdr_size = 0,
            .header_split = 0,    // Header Split disabled
            .hw_ip_checksum = 0,  // IP checksum offload disabled
            .hw_vlan_filter = 0,  // VLAN filtering disabled
            .jumbo_frame = 0,     // Jumbo Frame Support disabled
            .hw_strip_crc = 0,    // CRC stripped by hardware
        },
    .txmode =
        {
            .mq_mode = ETH_MQ_TX_NONE,
        },
};

/*
 * custom types
 */

/********* Misc *********/

// Get current time in ns
uint64_t get_cur_ns() {
  struct timespec ts;
  clock_gettime(CLOCK_REALTIME, &ts);
  uint64_t t = ts.tv_sec * 1000 * 1000 * 1000 + ts.tv_nsec;
  return t;
}

// For storing results
typedef struct LatencyResults_ {
  uint64_t *sjrn_times;
  uint64_t *sjrn_times_short;
  uint64_t *sjrn_times_long;
  uint64_t *work_ratios;
  uint64_t *work_ratios_short;
  uint64_t *work_ratios_long;
  uint32_t (*queue_lengths)[3];
  size_t count;
  size_t count_short;
  size_t count_long;
} LatencyResults;

/********* Distributions *********/

// Exponential Distribution
typedef struct ExpDist_ {
  gsl_rng *r;
  double mu;
  uint64_t cur_ns;
} ExpDist;

void init_exp_dist(ExpDist *exp_dist, double mu) {
  gsl_rng_env_setup();
  const gsl_rng_type *T = gsl_rng_default;
  gsl_rng *r = gsl_rng_alloc(T);
  uint64_t cur_ns = get_cur_ns();
  ExpDist temp_exp_dist = {r, mu, cur_ns};
  memcpy(exp_dist, &temp_exp_dist, sizeof(ExpDist));
}

uint64_t exp_dist_next_arrival_ns(ExpDist *exp_dist) {
  exp_dist->cur_ns += gsl_ran_exponential(exp_dist->r, exp_dist->mu);
  return exp_dist->cur_ns;
}

uint64_t exp_dist_work_ns(ExpDist *exp_dist) {
  return gsl_ran_exponential(exp_dist->r, exp_dist->mu);
}

void free_exp_dist(const ExpDist *exp_dist) {
  gsl_rng_free(exp_dist->r);
  free((void *)exp_dist);
}

// Lognormal Distribution
typedef struct LognormalDist_ {
  gsl_rng *r;
  double mu;
  double sigma;
} LognormalDist;

void init_lognormal_dist(LognormalDist *lognormal_dist, double mu,
                         double sigma) {
  gsl_rng_env_setup();
  const gsl_rng_type *T = gsl_rng_default;
  gsl_rng *r = gsl_rng_alloc(T);
  LognormalDist temp_lognormal_dist = {r, mu, sigma};
  memcpy(lognormal_dist, &temp_lognormal_dist, sizeof(LognormalDist));
}

uint64_t lognormal_dist_work_ns(LognormalDist *lognormal_dist) {
  return gsl_ran_lognormal(lognormal_dist->r, lognormal_dist->mu,
                           lognormal_dist->sigma);
}

void free_lognormal_dist(LognormalDist *lognormal_dist) {
  gsl_rng_free(lognormal_dist->r);
  free((void *)lognormal_dist);
}

// Bimodal Distribution
typedef struct BimodalDist_ {
  gsl_rng *r;
  uint64_t work_1;
  uint64_t work_2;
  double ratio;
} BimodalDist;

void init_bimodal_dist(BimodalDist *bimodal_dist, uint64_t work_1_ns,
                       uint64_t work_2_ns, double ratio) {
  gsl_rng_env_setup();
  const gsl_rng_type *T = gsl_rng_default;
  gsl_rng *r = gsl_rng_alloc(T);
  BimodalDist temp_bimodal_dist = {r, work_1_ns, work_2_ns, ratio};
  memcpy(bimodal_dist, &temp_bimodal_dist, sizeof(BimodalDist));
}

uint64_t bimodal_dist_work_ns(BimodalDist *bimodal_dist) {
  double num = gsl_ran_flat(bimodal_dist->r, 0.0, 1.0);
  if (num < bimodal_dist->ratio) {
    return bimodal_dist->work_1;
  } else {
    return bimodal_dist->work_2;
  }
}

void free_bimodal_dist(BimodalDist *bimodal_dist) {
  gsl_rng_free(bimodal_dist->r);
  free((void *)bimodal_dist);
}

// Trimodal Distribution
typedef struct TrimodalDist_ {
  gsl_rng *r;
  uint64_t work_1;
  uint64_t work_2;
  uint64_t work_3;
  double ratio1;
  double ratio2;
} TrimodalDist;

void init_trimodal_dist(TrimodalDist *trimodal_dist, uint64_t work_1_ns,
                        uint64_t work_2_ns, uint64_t work_3_ns, double ratio1,
                        double ratio2) {
  gsl_rng_env_setup();
  const gsl_rng_type *T = gsl_rng_default;
  gsl_rng *r = gsl_rng_alloc(T);
  TrimodalDist temp_trimodal_dist = {
      r, work_1_ns, work_2_ns, work_3_ns, ratio1, ratio2,
  };
  memcpy(trimodal_dist, &temp_trimodal_dist, sizeof(TrimodalDist));
}

uint64_t trimodal_dist_work_ns(TrimodalDist *trimodal_dist) {
  double num = gsl_ran_flat(trimodal_dist->r, 0.0, 1.0);
  if (num < trimodal_dist->ratio1) {
    return trimodal_dist->work_1;
  } else if (num < trimodal_dist->ratio1 + trimodal_dist->ratio2) {
    return trimodal_dist->work_2;
  } else {
    return trimodal_dist->work_3;
  }
}

void free_trimodal_dist(TrimodalDist *trimodal_dist) {
  gsl_rng_free(trimodal_dist->r);
  free((void *)trimodal_dist);
}

typedef struct Message_ {
  uint8_t recir_flag;
  uint8_t core_idx;
  uint8_t type;
  uint16_t seq_num;
  uint32_t queue_length[3];
  uint16_t client_id;
  uint32_t req_id;
  uint32_t pkts_length;
  uint64_t run_ns;
  uint64_t gen_ns;
} __attribute__((__packed__)) Message;

struct mbuf_table {
  uint32_t len;
  struct rte_mbuf *m_table[MAX_BURST_SIZE];
};

struct lcore_configuration {
  uint32_t vid;                                // virtual core id
  uint32_t port;                               // one port
  uint32_t tx_queue_id;                        // one TX queue
  uint32_t n_rx_queue;                         // number of RX queues
  uint32_t rx_queue_list[RX_QUEUE_PER_LCORE];  // list of RX queues
  struct mbuf_table tx_mbufs;                  // mbufs to hold TX queue
} __rte_cache_aligned;

/*
 * global variables
 */

uint32_t enabled_port_mask = 1;
uint32_t enabled_ports[RTE_MAX_ETHPORTS];
uint32_t n_enabled_ports = 0;
uint32_t n_rx_queues = 0;
uint32_t n_lcores = 0;

struct rte_mempool *pktmbuf_pool = NULL;
struct ether_addr port_eth_addrs[RTE_MAX_ETHPORTS];
struct lcore_configuration lcore_conf[MAX_LCORES];

uint8_t header_template[sizeof(struct ether_hdr) + sizeof(struct ipv4_hdr) +
                        sizeof(struct udp_hdr)];

/*
 * functions for generation
 */

// send packets, drain TX queue
static void send_pkt_burst(uint32_t lcore_id) {
  struct lcore_configuration *lconf = &lcore_conf[lcore_id];
  struct rte_mbuf **m_table = (struct rte_mbuf **)lconf->tx_mbufs.m_table;

  uint32_t n = lconf->tx_mbufs.len;
  uint32_t ret = rte_eth_tx_burst(lconf->port, lconf->tx_queue_id, m_table,
                                  lconf->tx_mbufs.len);
  if (unlikely(ret < n)) {
    do {
      rte_pktmbuf_free(m_table[ret]);
    } while (++ret < n);
  }
  lconf->tx_mbufs.len = 0;
}

// put packet into TX queue
static void enqueue_pkt(uint32_t lcore_id, struct rte_mbuf *mbuf) {
  struct lcore_configuration *lconf = &lcore_conf[lcore_id];
  lconf->tx_mbufs.m_table[lconf->tx_mbufs.len++] = mbuf;

  // enough packets in TX queue
  if (unlikely(lconf->tx_mbufs.len == MAX_BURST_SIZE)) {
    send_pkt_burst(lcore_id);
  }
}

/*
 * functions for initialization
 */

// init header template
static void init_header_template(void) {
  memset(header_template, 0, sizeof(header_template));
  struct ether_hdr *eth = (struct ether_hdr *)header_template;
  struct ipv4_hdr *ip =
      (struct ipv4_hdr *)((uint8_t *)eth + sizeof(struct ether_hdr));
  struct udp_hdr *udp =
      (struct udp_hdr *)((uint8_t *)ip + sizeof(struct ipv4_hdr));
  uint32_t pkt_len = sizeof(header_template) + sizeof(Message);

  // eth header
  eth->ether_type = rte_cpu_to_be_16(ETHER_TYPE_IPv4);
  struct ether_addr src_addr = {
      .addr_bytes = {0xaa, 0xbb, 0xcc, 0xdd, 0xee, 0xff}};
  struct ether_addr dst_addr = {
      .addr_bytes = {0x00, 0x11, 0x22, 0x33, 0x44, 0x55}};
  ether_addr_copy(&src_addr, &eth->s_addr);
  ether_addr_copy(&dst_addr, &eth->d_addr);

  // ip header
  char src_ip[] = IP_SRC;
  char dst_ip[] = IP_DST;
  int st1 = inet_pton(AF_INET, src_ip, &(ip->src_addr));
  int st2 = inet_pton(AF_INET, dst_ip, &(ip->dst_addr));
  if (st1 != 1 || st2 != 1) {
    fprintf(stderr, "inet_pton() failed.Error message: %s %s", strerror(st1),
            strerror(st2));
    exit(EXIT_FAILURE);
  }
  ip->total_length = rte_cpu_to_be_16(pkt_len - sizeof(struct ether_hdr));
  ip->version_ihl = 0x45;
  ip->type_of_service = 0;
  ip->packet_id = 0;
  ip->fragment_offset = 0;
  ip->time_to_live = 64;
  ip->next_proto_id = IPPROTO_UDP;
  uint32_t ip_cksum;
  uint16_t *ptr16 = (uint16_t *)ip;
  ip_cksum = 0;
  ip_cksum += ptr16[0];
  ip_cksum += ptr16[1];
  ip_cksum += ptr16[2];
  ip_cksum += ptr16[3];
  ip_cksum += ptr16[4];
  ip_cksum += ptr16[6];
  ip_cksum += ptr16[7];
  ip_cksum += ptr16[8];
  ip_cksum += ptr16[9];
  ip_cksum = ((ip_cksum & 0xffff0000) >> 16) + (ip_cksum & 0x0000ffff);
  if (ip_cksum > 65535) {
    ip_cksum -= 65535;
  }
  ip_cksum = (~ip_cksum) & 0x0000ffff;
  if (ip_cksum == 0) {
    ip_cksum = 0xffff;
  }
  ip->hdr_checksum = (uint16_t)ip_cksum;

  // udp header
  udp->src_port = htons(CLIENT_PORT);
  udp->dst_port = htons(SERVICE_PORT);
  udp->dgram_len = rte_cpu_to_be_16(pkt_len - sizeof(struct ether_hdr) -
                                    sizeof(struct ipv4_hdr));
  udp->dgram_cksum = 0;
}

// check link status
static void check_link_status(void) {
  const uint32_t check_interval_ms = 100;
  const uint32_t check_iterations = 90;
  uint32_t i, j;
  struct rte_eth_link link;
  for (i = 0; i < check_iterations; i++) {
    uint8_t all_ports_up = 1;
    for (j = 0; j < n_enabled_ports; j++) {
      uint32_t portid = enabled_ports[j];
      memset(&link, 0, sizeof(link));
      rte_eth_link_get_nowait(portid, &link);
      if (link.link_status) {
        printf("\tport %u link up - speed %u Mbps - %s\n", portid,
               link.link_speed,
               (link.link_duplex == ETH_LINK_FULL_DUPLEX) ? "full-duplex"
                                                          : "half-duplex");
      } else {
        all_ports_up = 0;
      }
    }

    if (all_ports_up == 1) {
      printf("check link status finish: all ports are up\n");
      break;
    } else if (i == check_iterations - 1) {
      printf("check link status finish: not all ports are up\n");
    } else {
      rte_delay_ms(check_interval_ms);
    }
  }
}

// initialize all status
static void init(void) {
  uint32_t i, j;

  // create mbuf pool
  printf("create mbuf pool\n");
  pktmbuf_pool = rte_mempool_create(
      "mbuf_pool", NB_MBUF, MBUF_SIZE, MBUF_CACHE_SIZE,
      sizeof(struct rte_pktmbuf_pool_private), rte_pktmbuf_pool_init, NULL,
      rte_pktmbuf_init, NULL, rte_socket_id(), 0);
  if (pktmbuf_pool == NULL) {
    rte_exit(EXIT_FAILURE, "cannot init mbuf pool\n");
  }

  // determine available ports
  printf("create enabled ports\n");
  uint32_t n_total_ports = 0;
  n_total_ports = rte_eth_dev_count();

  if (n_total_ports == 0) {
    rte_exit(EXIT_FAILURE, "cannot detect ethernet ports\n");
  }
  if (n_total_ports > RTE_MAX_ETHPORTS) {
    n_total_ports = RTE_MAX_ETHPORTS;
  }

  // get info for each enabled port
  struct rte_eth_dev_info dev_info;
  n_enabled_ports = 0;
  printf("\tports: ");
  for (i = 0; i < n_total_ports; i++) {
    if ((enabled_port_mask & (1 << i)) == 0) {
      continue;
    }
    enabled_ports[n_enabled_ports++] = i;
    rte_eth_dev_info_get(i, &dev_info);
    printf("%u ", i);
  }
  printf("\n");

  // find number of active lcores
  printf("create enabled cores\n\tcores: ");
  n_lcores = 0;
  for (i = 0; i < MAX_LCORES; i++) {
    if (rte_lcore_is_enabled(i)) {
      n_lcores++;
      printf("%u ", i);
    }
  }

  // ensure numbers are correct
  if (n_lcores % n_enabled_ports != 0) {
    rte_exit(EXIT_FAILURE,
             "number of cores (%u) must be multiple of ports (%u)\n", n_lcores,
             n_enabled_ports);
  }

  uint32_t rx_queues_per_lcore = RX_QUEUE_PER_LCORE;
  uint32_t rx_queues_per_port =
      rx_queues_per_lcore * n_lcores / n_enabled_ports;
  uint32_t tx_queues_per_port = n_lcores / n_enabled_ports;

  if (rx_queues_per_port < rx_queues_per_lcore) {
    rte_exit(EXIT_FAILURE,
             "rx_queues_per_port (%u) must be >= rx_queues_per_lcore (%u)\n",
             rx_queues_per_port, rx_queues_per_lcore);
  }

  // assign each lcore some RX queues and a port
  printf("set up %d RX queues per port and %d TX queues per port\n",
         rx_queues_per_port, tx_queues_per_port);
  uint32_t portid_offset = 0;
  uint32_t rx_queue_id = 0;
  uint32_t tx_queue_id = 0;
  uint32_t vid = 0;
  for (i = 0; i < MAX_LCORES; i++) {
    if (rte_lcore_is_enabled(i)) {
      lcore_conf[i].vid = vid++;
      lcore_conf[i].n_rx_queue = rx_queues_per_lcore;
      for (j = 0; j < rx_queues_per_lcore; j++) {
        lcore_conf[i].rx_queue_list[j] = rx_queue_id++;
      }
      lcore_conf[i].port = enabled_ports[portid_offset];
      lcore_conf[i].tx_queue_id = tx_queue_id++;
      if (rx_queue_id % rx_queues_per_port == 0) {
        portid_offset++;
        rx_queue_id = 0;
        tx_queue_id = 0;
      }
    }
  }

  // initialize each port
  for (portid_offset = 0; portid_offset < n_enabled_ports; portid_offset++) {
    uint32_t portid = enabled_ports[portid_offset];

    int32_t ret = rte_eth_dev_configure(portid, rx_queues_per_port,
                                        tx_queues_per_port, &port_conf);
    if (ret < 0) {
      rte_exit(EXIT_FAILURE, "cannot configure device: err=%d, port=%u\n", ret,
               portid);
    }
    rte_eth_macaddr_get(portid, &port_eth_addrs[portid]);

    // initialize RX queues
    for (i = 0; i < rx_queues_per_port; i++) {
      ret = rte_eth_rx_queue_setup(
          portid, i, NB_RXD, rte_eth_dev_socket_id(portid), NULL, pktmbuf_pool);
      if (ret < 0) {
        rte_exit(EXIT_FAILURE, "rte_eth_rx_queue_setup: err=%d, port=%u\n", ret,
                 portid);
      }
    }

    // initialize TX queues
    for (i = 0; i < tx_queues_per_port; i++) {
      ret = rte_eth_tx_queue_setup(portid, i, NB_TXD,
                                   rte_eth_dev_socket_id(portid), NULL);
      if (ret < 0) {
        rte_exit(EXIT_FAILURE, "rte_eth_tx_queue_setup: err=%d, port=%u\n", ret,
                 portid);
      }
    }

    // start device
    ret = rte_eth_dev_start(portid);
    if (ret < 0) {
      rte_exit(EXIT_FAILURE, "rte_eth_dev_start: err=%d, port=%u\n", ret,
               portid);
    }

    rte_eth_promiscuous_enable(portid);

    char mac_buf[ETHER_ADDR_FMT_SIZE];
    ether_format_addr(mac_buf, ETHER_ADDR_FMT_SIZE, &port_eth_addrs[portid]);
    printf("initiaze queues and start port %u, MAC address:%s\n", portid,
           mac_buf);
  }

  if (!n_enabled_ports) {
    rte_exit(EXIT_FAILURE,
             "all available ports are disabled. Please set portmask.\n");
  }
  check_link_status();
  init_header_template();
}

/*
 * misc
 */

static uint64_t timediff_in_us(uint64_t new_t, uint64_t old_t) {
  return (new_t - old_t) * 1000000UL / rte_get_tsc_hz();
}

#endif  // UTIL_H_

ext: h

max_stars_repo_name: netx-repo/RackSched

max_issues_repo_name: netx-repo/RackSched

max_forks_repo_name: netx-repo/RackSched


text:
//
//  mcmc_solar.c
//
//
//  Created by Juan Pablo Molano on 11/11/15.
//
//

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <math.h>
#define USAGE "./mcmc_solar.x n_burn n_steps"

//------------------------------------------------------------------------------------------------------------

float *load_matrix(char *filename, int *n, int *m);
void get_data( float *array , float *x, int n_row);
float *reserva(int n);
void print_data(float *array, float *array2, float *array3,  float *array4 , float *array5 , float *array6 , float *array7 , float *array8,  int n_puntos);
void get_patch( float *array,  float *x , int n_row);
void my_model( float *x_obs, float *array, float a, float b, float c , float d,  int n_row);
float likelihood(float *y_obs, float *y_model, int n_row);
int min_likelihood( float *array , int n_row);

//------------------------------------------------------------------------------------------------------------

int main(int argc, char **argv){
    float *matrix , *tim, *man, *count, *a_walk, *b_walk, *c_walk, *d_walk, *l_walk, *y_init, *y_prime, *best_y;
    srand48(time(NULL));
    float a_prime, b_prime, c_prime, d_prime, l_prime, l_init, gamma, alpha, beta , best_a , best_b, best_c, best_d;
    int n_row, n_cols , i, rlin = 0;

    a_walk = reserva(atof(argv[1])+atof(argv[2]));
    b_walk = reserva(atof(argv[1])+atof(argv[2]));
    c_walk = reserva(atof(argv[1])+atof(argv[2]));
    d_walk = reserva(atof(argv[1])+atof(argv[2]));
    l_walk = reserva(atof(argv[1])+atof(argv[2]));
    y_init = reserva(n_row);
    y_prime = reserva(n_row);

    a_walk[0] = drand48();
    b_walk[0] = drand48();
    c_walk[0] = drand48();
    d_walk[0] = drand48();

    matrix = load_matrix( "sol.dat" , &n_row, &n_cols);

    man = reserva( n_row);
    tim = reserva(n_row);
    man = reserva(n_row);
    best_y = reserva( n_row);
    get_data( matrix , tim , n_row);
    get_patch( matrix ,  man  ,  n_row);

    for (i=0;i<n_row;i++) {
        if (man[i] == -99) {
            man[i] = 0.0;
        }
    }

    my_model( tim , y_init , a_walk[0] , b_walk[0] , c_walk[0] , d_walk[0] , n_row );

    l_walk[0] = likelihood( man , y_init, n_row);

    const gsl_rng_type * T;
    gsl_rng * r;
    gsl_rng_env_setup();
    T = gsl_rng_default;
    r = gsl_rng_alloc (T);

    for (i=0; i < atof(argv[1]) + atof(argv[2]); i++) {

        a_prime = gsl_ran_gaussian(r, 0.1) + a_walk[i];
        b_prime = gsl_ran_gaussian(r, 0.1) + b_walk[i];
        c_prime = gsl_ran_gaussian(r, 0.1) + c_walk[i];
        d_prime = gsl_ran_gaussian(r, 0.1) + d_walk[i];

        my_model( tim, y_init, a_walk[i], b_walk[i], c_walk[i], d_walk[i], n_row);
        my_model(tim, y_prime ,a_prime, b_prime, c_prime, d_prime, n_row);

        l_prime = likelihood( man, y_prime , n_row);
        l_init = likelihood( man, y_init ,  n_row);

        gamma = l_init - l_prime;
        if (gamma>=0.0) {
            a_walk[i+1] = a_prime;
            b_walk[i+1] = b_prime;
            c_walk[i+1] = c_prime;
            d_walk[i+1] = d_prime;
            l_walk[i+1] = l_prime;
        }
        else{

            beta = drand48();
            alpha = exp(gamma);
            if (beta <= alpha) {

                a_walk[i+1]  = a_prime;
                b_walk[i+1]  = b_prime;
                c_walk[i+1]  = c_prime;
                d_walk[i+1]  = d_prime;
                l_walk[i+1]  = l_prime;
            }
            else{

                a_walk[i+1] = a_walk[i];
                b_walk[i+1] = b_walk[i];
                c_walk[i+1] = c_walk[i];
                d_walk[i+1] = d_walk[i];
                l_walk[i+1] = l_init;

            }

        }
    }

    i = min_likelihood( l_walk , n_row);
    best_a = a_walk[i];
    best_b = b_walk[i];
    best_c = c_walk[i];
    best_d = d_walk[i];

    my_model( tim , best_y , best_a, best_b , best_c , best_d , n_row);

    print_data( tim , man,  l_walk , a_walk , b_walk , c_walk , d_walk , best_y, n_row);

    gsl_rng_free (r);

    return 0;
}

//------------------------------------------------------------------------------------------------------------

float *load_matrix(char *filename, int *n, int *m){
    float *matrix;
    FILE *in;
    int n_row = 4632, n_cols = 5;
    int i;
    int j;

    if(!(in=fopen(filename, "r"))){
        printf("Problem opening file %s\n", filename);
        exit(1);
    }

    matrix = malloc(n_row * n_cols * sizeof(float));

    for(i=0;i<n_row;i++){
        for(j=0;j<n_cols;j++){
            fscanf(in, "%f", &matrix[i*n_cols + j]);
        }
    }
    *n = n_row;
    *m = n_cols;
    return matrix;
}

//------------------------------------------------------------------------------------------------------------

void get_data( float *array , float *x , int n_row){

    int i,j;

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

            x[i] = array[i*5] + array[i*5 + 1]/12 + array[i*5 + 2]/365;

    }

}

//------------------------------------------------------------------------------------------------------------

void get_patch( float *array,  float *x , int n_row){

    int i,j;

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

        x[i] = array[i*5 + 3];
    }

}

//------------------------------------------------------------------------------------------------------------

float *reserva(int n){
    float *array;
    int i;
    if(!(array = malloc(n * sizeof(float)))){
        printf("Problema en reserva\n");
        exit(1);
    }
    for(i=0;i<n ;i++){
        array[i] = 0.0;
    }
    return array;
}

//------------------------------------------------------------------------------------------------------------

void print_data(float *array, float *array2 , float *array3 , float *array4 , float *array5 , float *array6, float *array7, float *array8 , int n_puntos){
    int i;
    for(i=0;i<n_puntos;i++){
        printf("%f \t %f \t %f \t %f \t %f \t %f \t %f \t %f \n", array[i] , array2[i] , array3[i] , array4[i] , array5[i] , array6[i], array7[i], array8[i]);
    }
}

//------------------------------------------------------------------------------------------------------------

float likelihood(float *y_obs, float *y_model, int n_row){

    int i;
    float c=0;

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

        c = c + ( y_obs[i] - y_model[i]);
    }
    c = (1.0/2.0)*pow(c,2);

    return c;

}

//------------------------------------------------------------------------------------------------------------

void my_model( float *x_obs, float *array , float a, float b, float c , float d,  int n_row){

    int i;
    for (i=0; i<n_row; i++) {
        array[i] = a*cos((2*3.141592653589793/d)*x_obs[i] + b) + c;
    }
}

//------------------------------------------------------------------------------------------------------------

int min_likelihood( float *array , int n_row){

    int i;
    int ind = 0;

    for(i = 1; i < n_row; i ++){
        if ( array[i] < ind) {
            ind = i;
        }

    }

    return ind;
}

ext: c


text:
//====---- Sudoku/Board_Section.h                                     ----====//
//
// Internal in-between object for provide access to a section of Board.
//====--------------------------------------------------------------------====//
#pragma once

#include "Board_Section_iterator.h"
#include "Board_Section_traits.h" // Convenience include
#include "Location.h"
#include "Location_Utilities.h"
#include "exceptions.h"
#include "traits.h"

#include <gsl/gsl> // index

#include <type_traits>

#include "Board.fwd.h" // Forward declarations

#include <cassert>

namespace Sudoku::Board_Section
{
template<typename T, int N, Section S, bool is_const>
class [[nodiscard]] Board_Section_
{
	using index          = gsl::index;
	using Location       = ::Sudoku::Location<N>;
	using Location_Block = ::Sudoku::Location_Block<N>;
	using OwnerT =
		std::conditional_t<is_const, Board<T, N> const&, Board<T, N>&>;

	// explicit friends for conversion
	friend class Board_Section_<T, N, Section::row, true>;
	friend class Board_Section_<T, N, Section::row, false>;
	friend class Board_Section_<T, N, Section::col, true>;
	friend class Board_Section_<T, N, Section::col, false>;
	friend class Board_Section_<T, N, Section::block, true>;
	friend class Board_Section_<T, N, Section::block, false>;

public:
	using value_type             = T;
	using pointer                = std::conditional_t<is_const, T const*, T*>;
	using reference              = std::conditional_t<is_const, T const&, T&>;
	using iterator               = Section_iterator<T, N, S, is_const, false>;
	using const_iterator         = Section_iterator<T, N, S, true, false>;
	using reverse_iterator       = Section_iterator<T, N, S, is_const, true>;
	using const_reverse_iterator = Section_iterator<T, N, S, true, true>;

	template<typename idT, typename = std::enable_if_t<is_int_v<idT>>>
	constexpr Board_Section_(OwnerT board, idT id) noexcept
		: board_(board), id_(id)
	{
		assert(is_valid_size<N>(id));
	}
	constexpr Board_Section_(OwnerT board, Location loc) noexcept
		: board_(board), id_(to_id(loc))
	{
		assert(is_valid<N>(loc));
	}
	// Conversion to other Section type (maintaining or gaining const)
	template<
		Section Sx,
		bool B,
		typename = std::enable_if_t<!B || B == is_const>>
	constexpr Board_Section_(
		Board_Section_<T, N, Sx, B> other, index pivot_elem) noexcept
		: board_(other.board_), id_(to_id(other.location(pivot_elem)))
	{
	}
	// [[implicit]] conversion to const
	// NOLINTNEXTLINE(google-explicit-constructor, hicpp-explicit-conversions)
	constexpr operator Board_Section_<T, N, S, true>() const noexcept
	{
		return Board_Section_<T, N, S, true>(board_, id_);
	}

	[[nodiscard]] static constexpr int size() noexcept { return elem_size<N>; }
	[[nodiscard]] constexpr index id() const noexcept { return id_; }
	[[nodiscard]] constexpr Location location(index elem) const noexcept;

	[[nodiscard]] constexpr reference front() noexcept { return (*this)[0]; }
	[[nodiscard]] constexpr reference back() noexcept
	{
		return (*this)[static_cast<index>(size()) - 1];
	}
	[[nodiscard]] constexpr reference operator[](const index elem) noexcept
	{
		return board_[location(elem)];
	}
	[[nodiscard]] constexpr T const& operator[](const index elem) const noexcept
	{
		return board_[location(elem)];
	}
	// Checked access
	[[nodiscard]] constexpr reference at(const index elem);

	constexpr iterator begin() noexcept { return iterator(&board_, id_, 0); }
	constexpr iterator end() noexcept { return iterator(&board_, id_, size()); }
	[[nodiscard]] constexpr const_iterator cbegin() const noexcept
	{
		return const_iterator(&board_, id_, 0);
	}
	[[nodiscard]] constexpr const_iterator cend() const noexcept
	{
		return const_iterator(&board_, id_, size());
	}
	[[nodiscard]] constexpr const_iterator begin() const noexcept
	{
		return cbegin();
	}
	[[nodiscard]] constexpr const_iterator end() const noexcept
	{
		return cend();
	}
	constexpr reverse_iterator rbegin() noexcept
	{
		return reverse_iterator(&board_, id_, size() - 1);
	}
	constexpr reverse_iterator rend() noexcept
	{
		return reverse_iterator(&board_, id_, -1);
	}
	[[nodiscard]] constexpr const_reverse_iterator crbegin() const noexcept
	{
		return const_reverse_iterator(&board_, id_, size() - 1);
	}
	[[nodiscard]] constexpr const_reverse_iterator crend() const noexcept
	{
		return const_reverse_iterator(&board_, id_, -1);
	}
	[[nodiscard]] constexpr const_reverse_iterator rbegin() const noexcept
	{
		return crbegin();
	}
	[[nodiscard]] constexpr const_reverse_iterator rend() const noexcept
	{
		return crend();
	}

private:
	OwnerT board_;
	const gsl::index id_{};

	// Internal helper
	static constexpr gsl::index to_id(Location loc) noexcept;
};

//====--------------------------------------------------------------------====//
// Member Functions

// Checked access
template<typename T, int N, Section S, bool is_const>
[[nodiscard]] inline constexpr
	typename Board_Section_<T, N, S, is_const>::reference
	Board_Section_<T, N, S, is_const>::at(const index elem)
{
	if (!is_valid_size<N>(elem))
	{
		throw error::invalid_Location{"Board_Section::at(elem)"};
	}
	return board_.at(location(elem));
}

template<typename T, int N, Section S, bool is_const>
[[nodiscard]] constexpr Location<N>
	Board_Section_<T, N, S, is_const>::location(index elem) const noexcept
{
	assert(is_valid_size<N>(elem));
	switch (S)
	{
	case Section::row: return Location(id_, elem);
	case Section::col: return Location(elem, id_);
	case Section::block: return Location_Block(id_, elem);
	}
}

// Internal helper
template<typename T, int N, Section S, bool is_const>
inline constexpr gsl::index
	Board_Section_<T, N, S, is_const>::to_id(Location loc) noexcept
{
	switch (S)
	{
	case Section::row: return loc.row();
	case Section::col: return loc.col();
	case Section::block: return loc.block();
	}
}

} // namespace Sudoku::Board_Section

max_stars_repo_name: FeodorFitsner/fwkSudoku

max_issues_repo_name: FeodorFitsner/fwkSudoku

max_forks_repo_name: FeodorFitsner/fwkSudoku


text:
/* Copyright (c) Facebook, Inc. and its affiliates.
 * This source code is licensed under the MIT license found in the
 * LICENSE file in the root directory of this source tree.
 */

#include <cblas.h>
#include <complex.h>

#defmacro DEFINE_SVV(TNAME, TYPE, OPNAME, OPTEXT)

  #define OPER(a,b) OPTEXT

  TYPE TNAME##Vector_##OPNAME(TYPE *a, TYPE *b, int n)
  {
    /* try cblas */
  #if #TYPE == "float" && #OPNAME == "dot"
    return cblas_sdot(n, a, 1, b, 1);
  #elif #TYPE == "double" && #OPNAME == "dot"
    return cblas_ddot(n, a, 1, b, 1);
  #else
    int i;
    TYPE s = 0;
  # pragma unroll(i)
    for(i=0;i<n;i++)
      s += OPER(a[i],b[i]);
    return s;
  #endif
  }

#endmacro

#defmacro FORALLTYPES(macro,...)
  macro(Int, int, __VA_ARGS__)
  macro(Long, long, __VA_ARGS__)
  macro(Float, float, __VA_ARGS__)
  macro(Double, double, __VA_ARGS__)
  macro(Complex, complex double, __VA_ARGS__)
#endmacro

FORALLTYPES(DEFINE_SVV,dot,(a)*(b));
FORALLTYPES(DEFINE_SVV,sqrdist,((a)-(b))*((a)-(b)))

#define DEFINE_TYPEINFO(TNAME,TYPE) { #TNAME, #TYPE },
struct { const char *tname, *type; } typeInfo[] = {
  FORALLTYPES(DEFINE_TYPEINFO)
};

max_stars_repo_name: zydeco/CParser

max_issues_repo_name: zydeco/CParser

max_forks_repo_name: zydeco/CParser


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

#pragma once

#ifndef applicationwin32_3971e4cb_1be6_4066_b7ea_da679d40c987_h
#define applicationwin32_3971e4cb_1be6_4066_b7ea_da679d40c987_h

#include <gslib/string.h>
#include <gslib/std.h>
#include <ariel/application.h>
#include <windows.h>

#undef min
#undef max

__ariel_begin__

typedef LRESULT(__stdcall *fnwndproc)(HWND, UINT, WPARAM, LPARAM);
typedef list<string> arg_list;

extern void set_execute_path_as_directory();
extern void init_application_environment(app_env& env);
extern void init_application_environment(app_env& env, HINSTANCE hinst, HINSTANCE hprevinst, const gchar* argv[], int argc);

struct app_data
{
    HINSTANCE           hinst = nullptr;
    HWND                hwnd = nullptr;
    fnwndproc           wndproc = nullptr;
    arg_list            arglist;

public:
    bool install(const app_config& cfg, const app_env& env);
    int run();
};

struct app_env
{
    HINSTANCE           hinst = nullptr;
    HINSTANCE           hprevinst = nullptr;
    arg_list            arglist;

public:
    void init(HINSTANCE hinst, HINSTANCE hprevinst, const gchar* argv[], int argc) { init_application_environment(*this, hinst, hprevinst, argv, argc); }
};

__ariel_end__

#endif

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


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

#include <config.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_blas.h>

int
gsl_multimin_test_gradient (const gsl_vector *g, double epsabs)
{
  double norm;

  if (epsabs < 0.0)
    {
      GSL_ERROR ("absolute tolerance is negative", GSL_EBADTOL);
    }

  norm = gsl_blas_dnrm2(g);

  if (norm < epsabs)
    {
      return GSL_SUCCESS;
    }

  return GSL_CONTINUE;
}

int
gsl_multimin_test_size (const double size, double epsabs)
{
  if (epsabs < 0.0)
    {
      GSL_ERROR ("absolute tolerance is negative", GSL_EBADTOL);
    }

  if (size < epsabs)
    {
      return GSL_SUCCESS;
    }

  return GSL_CONTINUE;
}

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 <math.h>
#include <gsl/gsl_randist.h>

int main (void)
{
  double x;
  const gsl_rng_type * T;
  gsl_rng * r;

  gsl_rng_env_setup();

  T = gsl_rng_default;
  r = gsl_rng_alloc(T);

  for (x = 0.1 ; x < 2; x+= 0.1)
    {
      double y0 = exp(x);
      double err = 0.1*y0;

      printf("%g %g %g\n", x, y0 + gsl_ran_gaussian(r, err), err);
    }
}

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
#include <math.h>
#include <stdlib.h>
#include <gsl/gsl_siman.h>

/* set up parameters for this simulated annealing run */

/* how many points do we try before stepping */
#define N_TRIES 200

/* how many iterations for each T? */
#define ITERS_FIXED_T 10

/* max step size in random walk */
#define STEP_SIZE 10

/* Boltzmann constant */
#define K 1.0

/* initial temperature */
#define T_INITIAL 0.002

/* damping factor for temperature */
#define MU_T 1.005
#define T_MIN 2.0e-6

gsl_siman_params_t params
  = {N_TRIES, ITERS_FIXED_T, STEP_SIZE,
     K, T_INITIAL, MU_T, T_MIN};

/* now some functions to test in one dimension */
double E1(void *xp)
{
  double x = * ((double *) xp);

  return exp(-pow((x-1.0),2.0))*sin(8*x);
}

double M1(void *xp, void *yp)
{
  double x = *((double *) xp);
  double y = *((double *) yp);

  return fabs(x - y);
}

void S1(const gsl_rng * r, void *xp, double step_size)
{
  double old_x = *((double *) xp);
  double new_x;

  double u = gsl_rng_uniform(r);
  new_x = u*2*step_size - step_size + old_x;

  memcpy(xp, &new_x, sizeof(new_x));
}

void P1(void *xp)
{
  printf("%12g", *((double *) xp));
}

int
main(int argc, char *argv[])
{
  gsl_rng_type * T;
  gsl_rng * r;

  double x_initial = 15.5;

  gsl_rng_env_setup();

  T = gsl_rng_default;
  r = gsl_rng_alloc(T);

  gsl_siman_solve(r, &x_initial, E1, S1, M1, P1,
                  NULL, NULL, NULL,
                  sizeof(double), params);
  return 0;
}


text:
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <complex.h>

#include <float.h>

#include <fftw3.h>

#include "gauss_conv.h"
#include "seis_data.h"

#include "calculate_source.h"

void calculate_source(seis_data *obs, seis_data *syn, double *syn_src, double *obs_src, int src_no)
{
    double lambda1, lambda2; /* <---Regularization parameters for Tikhonov regularization */

    int npts = obs->num_samples;
    int num_rec = obs->num_rec;

    int N = npts;       /* fft size */
    int Nfreq = N/2 +1; /* size of freq vector */

    int i, j; /* counters for loops */

    /*************************************************************************/

    /* allocate space for vector of sizes */
    int *nn = malloc(sizeof(*nn) * num_rec);

    /* allocate space for syn/obs and Fourier transforms */
    double *sd = fftw_malloc(sizeof(*sd) * N * num_rec);
    fftw_complex *SD = fftw_malloc(sizeof(*SD) * Nfreq * num_rec);

    /* allocate space for syn/obs source and Fourier transforms  */
    double *k = fftw_malloc(sizeof(*k) * N);
    fftw_complex *K = fftw_malloc(sizeof(*K) * Nfreq);

    /* Allocate space for Fourier transform of Green's function (impulse response) */
    fftw_complex *G = fftw_malloc(sizeof(*G) * Nfreq * num_rec);

    /* allocate some temporary working space */
    fftw_complex *temp = fftw_malloc(sizeof(*temp) * 3 * N);
    fftw_complex *num = temp + N;
    fftw_complex *den = temp + 2 * N;

    /* check memory allocations */
    if (!nn || !sd || !SD || !k || !K || !G || !temp)
    {
        fprintf(stderr,"Error: calculate_source: memory allocation failure\n");
        exit(1);
    }

    /*************************************************************************/

    /* fill array of sizes - all sizes are the same */
    for (i=0;i<num_rec;i++) nn[i] = N;

    /* make FFTW plans */
    fftw_plan plan_forward_sd, plan_forward_k, plan_backward_K;

    /* synth/obs data - forward only */
    plan_forward_sd = fftw_plan_many_dft_r2c(1,nn,num_rec,sd,NULL,1,N,SD,NULL,1,Nfreq,FFTW_MEASURE);

    /* synth/obs source - forward and backward */
    plan_forward_k = fftw_plan_dft_r2c_1d(N,k,K,FFTW_MEASURE);
    plan_backward_K = fftw_plan_dft_c2r_1d(N,K,k,FFTW_MEASURE);

    /*************************************************************************/
    /* copy synth source and data over */
    for (j=0;j<npts;j++)
    {
        k[j] = syn_src[j];
    }
    for (j=npts;j<N;j++)
    {
        /* zero padding */
        k[j] = 0.0;
    }

    for (i=0;i<num_rec;i++)
    {
        for (j=0;j<npts;j++)
        {
            sd[i*N + j] = syn->traces[src_no*npts*num_rec + i*npts + j];
        }
        for (j=npts;j<N;j++)
        {
            /* zero padding */
            sd[i*N + j] = 0.0;
        }
    }

    /* Execute FFTW plans */
    fftw_execute(plan_forward_sd);
    fftw_execute(plan_forward_k);

    /*************************************************************************/

    /* calculate denominator of first fraction */
    /* zero out num and den of second fraction */
    for (j=0;j<N;j++)
    {
        /* pre computing this reduces the number of divisions by a factor of num_rec */
        temp[j] = 1.0 / (conj(K[j]) * K[j] + lambda1);
        num[j] = 0.0;
        den[j] = 0.0;
    }

    /* solve for Green's function (freq domain) */
    for (i=0;i<num_rec;i++)
    {
        /* alternatively, call spectral_division_wl or spectral_division_tik */
        /* in these cases temp does not need to be the inverse as above */
        for (j=0;j<N;j++)
        {
            G[i*N + j] = conj(K[j]) * SD[i*N + j] * temp[j];
        }
    }

    /*************************************************************************/

    /* copy over observed data and take Fourier transform */
    for (i=0;i<num_rec;i++)
    {
        for (j=0;j<npts;j++)
        {
            sd[i*N + j] = obs->traces[src_no*npts*num_rec + i*npts + j];
        }
        for (j=npts;j<N;j++)
        {
            /* zero padding */
            sd[i*N + j] = 0.0;
        }
    }

    fftw_execute(plan_forward_sd);

    /*************************************************************************/
    /* solve for obs source in least squares sense */

    for (i=0;i<num_rec;i++)
    {
        for (j=0;j<N;j++)
        {
            num[j] += conj(G[i*N + j]) * SD[j];
            den[j] += conj(G[i*N + j]) * G[i*N + j];
        }
    }

    /* alternatively, call spectral_division_wl or spectral_division_tik */
    for (j=0;j<N;j++)
    {
        K[j] = num[j] / (den[j] +lambda2);
    }

    /* Execute FFTW plan for backwards transform of k */
    fftw_execute(plan_backward_K);

    /* scale k */
    double sf = 1.0 / (double)N;
    for (j=0;i<N;j++)
    {
        k[j] = sf * k[j];
    }

    /* get minimum phase reconstruction. write result in output vector */
    min_phase_reconstruction(k, obs_src, N);

}

void min_phase_reconstruction(const double *signal_in, double *signal_out, const int npts)
{
    /* construct a minimum phase wavelet from the data in signal_in                 */
    /* npts is the actual length of the data without any padding                    */
    /* signal_out should be allocated already and able to hold at least npts values */

    /* this restricts us to even N */
    int N = next_power2(npts*8);

    /* allocate temporary arrays */
    fftw_complex *x1 = malloc(sizeof(*x1) * N);
    fftw_complex *x2 = malloc(sizeof(*x2) * N);
    if (!x1 || !x2)
    {
        fprintf(stderr, "Error: min_phase_reconstruction: Memory allocation failure\n");
        exit(1);
    }

    /* make fftw plans */
    fftw_plan ft_plan_for, ft_plan_bac;
    ft_plan_for = fftw_plan_dft_1d(N, x1, x2, FFTW_FORWARD, FFTW_MEASURE);
    ft_plan_bac = fftw_plan_dft_1d(N, x2, x1, FFTW_BACKWARD, FFTW_MEASURE);

    int i;
    /* copy data over */
    for (i=0;i<npts;i++)
    {
        x1[i] = signal_in[i];
    }
    /* put zero padding at the end */
    for (i=npts;i<N;i++)
    {
        x1[i] = 0.0;
    }

    /* forward fft */
    fftw_execute(ft_plan_for);

    /* take log of abs */ /* <---TODO: need to handle values with 0 magnitude */
    for (i=0;i<N;i++)
    {
        x2[i] = log(cabs(x2[i]));
    }

    /* inverse fft */
    fftw_execute(ft_plan_bac);

    /* for inverse fft, need to divide values by fft length */
    double scale_fact = 1.0 / (double)N;
    int half_ind = N/2;

    /* window signal */
    /* first element (zero frequency) is unchanged - take real part and scale */
    x1[0] = scale_fact * creal(x1[0]);

    /* take real part, scale, and window - positive freqs by 2, negative freqs by 0 */
    for (i=1;i<N/2;i++)
    {
        x1[i] = 2.0 * scale_fact * creal(x1[i]);
        x1[i + half_ind] = 0.0;
    }
    /* element at N/2 (nyquist frequency) is unchanged - take real part and scale */
    x1[half_ind] = scale_fact * creal(x1[half_ind]);

    /* execute forward fft of windowed signal */
    fftw_execute(ft_plan_for);

    /* run signal through (complex) exponential */
    for (i=0;i<N;i++)
    {
        x2[i] = cexp(x2[i]);
    }

    /* take inverse fft */
    fftw_execute(ft_plan_bac);

    /* take real part and scale data. place result in output array*/
    for (i=0;i<npts;i++)
    {
        signal_out[i] = scale_fact * creal(x1[i]);
    }

    /* clean up */
    fftw_free(x1);
    fftw_free(x2);
    fftw_destroy_plan(ft_plan_for);
    fftw_destroy_plan(ft_plan_bac);

}

void spectral_division_wl(fftw_complex *num, fftw_complex *den, int N, double wtr, fftw_complex *out)
{
    int i;
    double tol = 2 * DBL_EPSILON;

    /* temporary to hold magnitude of den */
    double *d_mag = malloc(sizeof(*d_mag) * N);
    if (!d_mag)
    {
        fprintf(stderr, "Error: spectral_division_wl: Memory allocation failure\n");
        exit(1);
    }

    double d_amp_max = cabs(den[0]);
    /* get magnitude of den and pick max */
    for (i=0;i<N;i++)
    {
        d_mag[i] = cabs(den[i]);
        d_amp_max = d_mag[i] > d_amp_max ? d_mag[i] : d_amp_max;
    }

    double wtr_use = wtr * d_amp_max;

    fftw_complex den_use;
    for (i=0;i<N;i++)
    {
        if (d_mag[i] < tol)
        {
            /* den magnitude equal to zero */
            /* set value to wtr_use */
            den_use = wtr_use;
        }
        else if (d_mag[i] < wtr_use)
        {
            /* den magnitude less than water level */
            /* keep phase the same, but change magnitude to wtr_use */
            den_use = wtr_use * (den[i] / d_mag[i]);
        }
        else
        {
            /* den magnitude greater than water level */
            /* no change */
            den_use = den[i];
        }

        out[i] = num[i] / den_use;
    }

    free(d_mag);
}

void spectral_division_tik(fftw_complex *num, fftw_complex *den, int N, double lambda, fftw_complex *out)
{
    int i;
    for (i=0;i<N;i++)
    {
        out[i] = (conj(den[i]) * num[i]) / (conj(den[i]) * den[i] + lambda);
    }
}

void fft_convolution(double *x, double *h, double *y, int npts)
{
    /* convolution of two equal length signals, x and h.    */
    /* x and h are both npts long                           */
    /* the result (y) is length 2*npts-1                    */
    /* y should already be allocated                        */

    /* fft size needs to be at least 2*npts -1 */
    int N = next_power2(2 * npts);

    int fft_freq_sz = N/2 + 1;

    double *x_in = fftw_malloc(sizeof(*x_in) * N * 2);
    fftw_complex *X = fftw_malloc(sizeof(*X) * fft_freq_sz * 2);
    if (!x_in || !X)
    {
        fprintf(stderr,"Error: fft_convolution: memory allocation failure\n");
        exit(1);
    }

    double *h_in = x_in + N;
    fftw_complex *H = X + fft_freq_sz;

    int how_many = 2;
    int n[2] = {N,N};
    int stride = 1;
    int rank = 1;
    fftw_plan fplan, bplan;
    fplan = fftw_plan_many_dft_r2c(rank,n,how_many,x_in,NULL,stride,N,X,NULL,stride,fft_freq_sz,FFTW_MEASURE);
    bplan = fftw_plan_dft_c2r_1d(N,X,x_in,FFTW_MEASURE);

    int i;
    /* copy data */
    for (i=0;i<npts;i++)
    {
        x_in[i] = x[i];
        h_in[i] = h[i];
    }
    /* add zero padding */
    for (i=npts;i<N;i++)
    {
        x_in[i] = 0.0;
        h_in[i] = 0.0;
    }

    /* execute forward plan */
    fftw_execute(fplan);

    /* perform multiplication - keep result in X */
    for (i=0;i<fft_freq_sz;i++)
    {
        X[i] *= H[i];
    }

    /* execute backward plan */
    fftw_execute(bplan);

    /* scale and copy data to output array */
    double sf = 1.0 / (double)N;
    for (i=0;i<2*npts-1;i++)
    {
        y[i] = sf * x_in[i];
    }

    /* clean up */
    fftw_destroy_plan(fplan);
    fftw_destroy_plan(bplan);

    fftw_free(x_in);
    fftw_free(X);

}

#include <lapacke.h>

void deconvolution_system(const double *x, double *h, const double *y, const int npts, const double lambda)
{
    /* NOTE: This is really only practical when npts is fairly small        */
    /*                                                                      */
    /* the convolution problem can be formulated as a matrix multiplication */
    /* X*h = y, where the columns of X contain shifted copies of x          */
    /* Therefore, the deconvolution problem is solving this overdetermined  */
    /* system for h. To do this, we solve (X'*X + lambda*I)*h = X'*y        */
    /* where ' denotes matrix transpose, I is the identity matrix, and      */
    /* lambda is a Tikhonov regularization parameter                        */
    /* x and h should be length npts and y is length (2*npts-1)             */

    int N = next_power2(2*npts);
    int freq_sz = N/2 + 1;
    int conv_sz = 2*npts -1;

    double *X = malloc(sizeof(*X) * npts * conv_sz);
    double *xx = malloc(sizeof(*xx) * 3 * npts);
    double *yy = malloc(sizeof(*yy) * conv_sz);
    double *s = malloc(sizeof(*s) * npts);  /* s holds the singular values */

    /* copy data and zero pad */
    int i, j;
    for (i=0;i<npts;i++)
    {
        xx[i] = 0.0;
        xx[i+npts] = x[i];
        xx[i+2*npts] = 0.0;
        yy[i] = y[i];
    }
    for (i=npts;i<conv_sz;i++)
    {
        yy[i] = y[i];
    }

    for (i=0;i<npts;i++)
    {
        for (j=0;j<conv_sz;j++)
        {
            X[i*conv_sz + j] = xx[(npts - i) + j];
        }
    }

    char trans = 'N';

    int m = conv_sz; /* number of rows in X */
    int n = npts; /* number of cols in X */

    int nrhs = 1; /* number of right hand side (1) */

    int lda = m; /* leading dim of A, m for col major layout */
    int ldb = m; /* leading dim of B */

    double rcond = lambda;
    int rank;

    int info;

    info = LAPACKE_dgelss(LAPACK_COL_MAJOR, m, n, nrhs, X, lda, yy, ldb, s, rcond, &rank);

    if (info != 0)
    {
        fprintf(stderr, "Error: deconvolution_system: Lapack failure\n");
        exit(1);
    }

    /* copy solution over to h */
    for (i=0;i<npts;i++)
    {
        h[i] = yy[i];
    }

}

max_stars_repo_name: dt-jackson/FWI

max_issues_repo_name: dt-jackson/FWI

max_forks_repo_name: dt-jackson/FWI


text:
//
//  Extrap1d.h

#ifndef __NeutronDetectorSim__Extrap1d__
#define __NeutronDetectorSim__Extrap1d__

#include <iostream>
#include <vector>
#include <eigen3/Eigen/Dense>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_spline.h>
#include "ConsolePrint.hpp"

class Extrap1d {
public:
	Extrap1d(Eigen::ArrayXd x, Eigen::ArrayXd y);
	Extrap1d(std::vector<double> x, std::vector<double> y);
	~Extrap1d();
	double operator()(const double &x);
	Eigen::ArrayXd operator()(const Eigen::ArrayXd &x);
private:
	double min, max;
	Eigen::ArrayXd x;
	Eigen::ArrayXd y;
	gsl_spline *interpolator;
	gsl_interp_accel *accelerator;
};

#endif /* defined(__NeutronDetectorSim__Extrap1d__) */

max_stars_repo_name: SkyToGround/DetectorSensitivity

max_issues_repo_name: SkyToGround/DetectorSensitivity

max_forks_repo_name: SkyToGround/DetectorSensitivity


text:
#ifndef GLOBALS

#define GLOBALS

#include <gsl/gsl_matrix.h>

extern int NumEl;
extern int NumEdges;
extern int NumNodes;
extern int P;
extern int Np;
extern gsl_matrix *LIFT;
extern gsl_matrix *VolMat;
extern gsl_matrix *MassMatrix;

extern double *x;
extern double *y;
extern double *z;
extern double *b;
extern double *m1;
extern double *m2;
extern double *db;
extern double *dm1;
extern double *dm2;
extern double *dz;
extern double *dh;
extern double *nFriction;

extern double *NodalX;
extern double *NodalY;
extern double *NodalB;
extern double *Nodalz;
extern double *Nodalm1;
extern double *Nodalm2;
extern double *NodalnFriction;

extern double *A;
extern double *Q;

//extern double *Fhat1L;
//extern double *Fhat1R;
//extern double *Fhat2L;
//extern double *Fhat2R;
extern double *RHSA;
extern double *RHSQ;

extern double dt;
extern double max_lambda;
extern double mindh;

extern const double g;
extern const double H0;
extern const double VELZERO;

extern int *WD;		// keeps track of the wet/dry status of an element

#endif

max_stars_repo_name: evalseth/DG-RAIN

max_issues_repo_name: evalseth/DG-RAIN

max_forks_repo_name: evalseth/DG-RAIN


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

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

int main() {
  int n = 4;
  double** a = alloc_dmatrix(n, n);
  mat_elem(a, 0, 0) = 0.959291425205444;
  mat_elem(a, 0, 1) = 0.257508254123736;
  mat_elem(a, 0, 2) = 0.243524968724989;
  mat_elem(a, 0, 3) = 0.251083857976031;
  mat_elem(a, 1, 0) = 0.547215529963803;
  mat_elem(a, 1, 1) = 0.840717255983663;
  mat_elem(a, 1, 2) = 0.929263623187228;
  mat_elem(a, 1, 3) = 0.616044676146639;
  mat_elem(a, 2, 0) = 0.138624442828679;
  mat_elem(a, 2, 1) = 0.254282178971531;
  mat_elem(a, 2, 2) = 0.349983765984809;
  mat_elem(a, 2, 3) = 0.473288848902729;
  mat_elem(a, 3, 0) = 0.149294005559057;
  mat_elem(a, 3, 1) = 0.814284826068816;
  mat_elem(a, 3, 2) = 0.196595250431208;
  mat_elem(a, 3, 3) = 0.351659507062997;
  double* x = alloc_dvector(n);
  x[0] = 0.830828627896291;
  x[1] = 0.585264091152724;
  x[2] = 0.549723608291140;
  x[3] = 0.917193663829810;
  double* y = alloc_dvector(n);
  y[0] = 0.961898080855054;
  y[1] = 0.00463422413406744;
  y[2] = 0.774910464711502;
  y[3] = 0.817303220653433;
  double alpha = 2.3;
  double beta = 0.5;

  printf("a: "); fprint_dmatrix(stdout, n, n, a);
  printf("x: "); fprint_dvector(stdout, n, x);
  printf("y: "); fprint_dvector(stdout, n, y);

  cblas_dgemv(CblasColMajor, CblasNoTrans, n, n, alpha, mat_ptr(a), n, vec_ptr(x), 1,
              beta, vec_ptr(y), 1);

  printf("%10.5f * a * x + %10.5f * y: ", alpha, beta); fprint_dvector(stdout, n, y);

  free_dmatrix(a);
  free_dvector(x);
  free_dvector(y);
}

/* gemv.m

A = [0.959291425205444,0.257508254123736,0.243524968724989,0.251083857976031;0.547215529963803,0.840717255983663,0.929263623187228,0.616044676146639;0.138624442828679,0.254282178971531,0.349983765984809,0.473288848902729;0.149294005559057,0.814284826068816,0.196595250431208,0.351659507062997]
x = [0.830828627896291;0.585264091152724;0.549723608291140;0.917193663829810]
y = [0.961898080855054;0.00463422413406744;0.774910464711502;0.817303220653433]
2.3 * A * x + 0.5 * y

*/

max_stars_repo_name: t-sakashita/rokko

max_issues_repo_name: t-sakashita/rokko

max_forks_repo_name: t-sakashita/rokko


text:
/* stable/stable.h
 *
 * Main header file of Libstable. Contains all declarations of the
 * usable functions in the library.
 *
 * Copyright (C) 2013. Javier Royuela del Val
 *                     Federico Simmross Wattenberg
 *
 * 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; version 3 of the License.
 *
 * 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, see <http://www.gnu.org/licenses/>.
 *
 *
 *  Javier Royuela del Val.
 *  E.T.S.I. Telecomunicación
 *  Universidad de Valladolid
 *  Paseo de Belén 15, 47002 Valladolid, Spain.
 *  jroyval@lpi.tel.uva.es
 */

#ifndef _stable_H_
#define _stable_H_

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

#include <gsl/gsl_integration.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_vector.h>

#define TINY 1e-50
#define EPS 2.2204460492503131E-16

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

/******************************************************************************/
/*          Library parameters                                                */
/******************************************************************************/

extern FILE * FLOG;            // Log file (optional)
extern FILE * FINTEG;          // Integrand evaluations output file (debug purposes)

extern unsigned short THREADS;    // threads of execution (0 => total available)
extern unsigned short IT_MAX;     // Maximum # of iterations in quadrature methods
extern unsigned short SUBS;       // # of integration subintervals
extern unsigned short METHOD1;    // Integration method on main subinterval
extern unsigned short METHOD2;    // Integration method on second subinterval
extern unsigned short METHOD3;    // Integration method on third subinterval

extern unsigned short INV_MAXITER; // Maximum # of iterations inversion method

extern double relTOL;     // Relative error tolerance
extern double absTOL;     // Absolut error tolerance
//extern double FACTOR;   //
extern double ALPHA_TH;    // Alpha threshold
extern double BETA_TH;    // Beta threshold
extern double EXP_MAX;    // Exponent maximum value
extern double XXI_TH;     // Zeta threshold
extern double THETA_TH;   // Theta threshold

extern double AUX1; // Auxiliary values
extern double AUX2;

#ifdef DEBUG
extern unsigned int integ_eval; // # of integrand evaluations
#endif

/******************************************************************************/
/******************************************************************************/
// Particular cases
enum
  {
    NOVALID =-1,
    STABLE,
    ALPHA_1,
    GAUSS ,
    CAUCHY,
    LEVY,
    STABLE_B1,
    ALPHA_1_B1
  };

// Function to evaluate
enum
  {
    CDF,
    PDF
  };

// Quadrature methods
enum
  {
    STABLE_QAG2 = 0,
    STABLE_QUADSTEP,
    STABLE_QROMBPOL,
    STABLE_QROMBRAT,
    STABLE_QNG,
    STABLE_QAG1,
    STABLE_QAG5,
    STABLE_VECT
  };

/************************************************************************
 ************************************************************************
 * Scalar methods                                                       *
 ************************************************************************
 ************************************************************************/
/* Scalar methods are those that, for each thread of execution, obtain a
   single evaluation of the desired function, at a single point.
*/

/******************************************************************************/
/*    Stable distribution structure.                                          */
/******************************************************************************/
struct StableDistStruct
  {
    /* Parameters:
    0-parametrization describen in Nolan, 1997 is employed by default
        alpha : stability index
        beta : skewness parameter
        scale: scale parameter
        mu_0 : 0-parametrization location parameter
        mu_1 : correspondig 1-parametrization location parameter    */
    double alpha;
    double beta;
    double sigma;
    double mu_0;
    double mu_1;

    /* Particular cases indicator (Gauss, Cauchy, Levy distribution, alpha==1, etc.) */
    int ZONE;

    /* Pointers to pdf and cdf evaluation functions */
    double(*stable_pdf_point)(struct StableDistStruct *, const double, double *);
    double(*stable_cdf_point)(struct StableDistStruct *, const double, double *);

    /* Precalculated values. */
    double alphainvalpha1;  /* alpha/(alpha-1)*/
    double xi;            /* -beta*tan(alpha*pi/2)*/
    double theta0;        /* 1/alpha*atan(beta*(tan(alpha*pi/2))=atan(-xi)/alpha;*/
    double c1, c2_part, c3;  /* additive and multiplicative constants*/
    double k1;     /* cos(alpha*theta0)^(1/(alpha-1)) = (1+xi^2)^(-0.5/(alpha-1));*/
    double S;     /* (1+xi^2)^(1/(2*alpha));*/
    double Vbeta1; /*pow(1/dist->alpha,dist->alphainvalpha1) *
                     (dist->alpha-1)*pow(-cos(dist->alpha*PI_2),1/(dist->alpha-1))*/

    /* These ones change from point to point of evaluation */
    double theta0_; /* theta0_ = +-theta0 */
    double beta_;
    double xxipow;  /* (x-xi)^(alpha/(alpha-1))*/

    /* gsl integration workspace */
    gsl_integration_workspace * gslworkspace;

    /* gsl random numbers generator */
    gsl_rng * gslrand;
  };

typedef struct StableDistStruct StableDist;
/******************************************************************************/
/*        Auxiliary functions                                                 */
/******************************************************************************/
unsigned int stable_get_THREADS();
void   stable_set_THREADS(unsigned int threads);

unsigned int stable_get_IT_MAX();
void   stable_set_IT_MAX(unsigned int itmax);

unsigned int stable_get_INV_MAXITER();
void   stable_set_INV_MAXITER(unsigned int invmaxiter);

int stable_get_METHOD1();
void stable_set_METHOD1(int method);

int stable_get_METHOD2();
void stable_set_METHOD2(int method);

int stable_get_METHOD3();
void stable_set_METHOD3(int method);

double stable_get_relTOL();
void   stable_set_relTOL(double reltol);

double stable_get_absTOL();
void   stable_set_absTOL(double abstol);

double stable_get_ALPHA_TH();
void   stable_set_ALPHA_TH(double alphath);

double stable_get_BETA_TH();
void   stable_set_BETA_TH(double betath);

double stable_get_XXI_TH();
void   stable_set_XXI_TH(double xxith);

double stable_get_THETA_TH();
void   stable_set_THETA_TH(double thetath);

FILE * stable_get_FINTEG();
FILE * stable_set_FINTEG(char * filename);

FILE * stable_get_FLOG();
FILE * stable_set_FLOG(char * filename);

StableDist *stable_create(double alpha, double beta, double sigma, double mu,
                          int parametrization);

StableDist *stable_copy(StableDist *src_dist);

void stable_free(StableDist *dist);

int stable_setparams(StableDist *dist,
                     double alpha, double beta, double sigma, double mu,
                     int parametrization);

int stable_checkparams(double alpha, double beta, double sigma, double mu,
                       int parametrization);

void error_handler (const char * reason, const char * file,
               int line, int gsl_errno);

/******************************************************************************/
/*   PDF in particular cases                                                  */
/******************************************************************************/

double stable_pdf_point_GAUSS(StableDist *dist, const double x, double *err);

double stable_pdf_point_CAUCHY(StableDist *dist, const double x, double *err);

double stable_pdf_point_LEVY(StableDist *dist, const double x, double *err);

/******************************************************************************/
/*   PDF otherwise                                                            */
/******************************************************************************/

double stable_pdf_point_STABLE(StableDist *dist, const double x, double *err);

double stable_pdf_point_ALPHA_1(StableDist *dist, const double x, double *err);

double stable_pdf_point(StableDist *dist, const double x, double *err);

void stable_pdf(StableDist *dist, const double* x, const unsigned int Nx,
                double *pdf, double *err);

/******************************************************************************/
/*   PDF integrand functions                                                  */
/******************************************************************************/

double stable_pdf_g(double theta, void *dist);
double stable_pdf_g_aux1(double theta, void *args);
double stable_pdf_g_aux2(double theta, void *args);

/******************************************************************************/
/*   CDF in particular cases                                                  */
/******************************************************************************/

double stable_cdf_point_GAUSS(StableDist *dist, const double x, double *err);

double stable_cdf_point_CAUCHY(StableDist *dist, const double x, double *err);

double stable_cdf_point_LEVY(StableDist *dist, const double x, double *err);

/******************************************************************************/
/*   CDF otherwise                                                            */
/******************************************************************************/

double stable_cdf_point_STABLE(StableDist *dist, const double x, double *err);

double stable_cdf_point_ALPHA_1(StableDist *dist, const double x, double *err);

double stable_cdf_point(StableDist *dist, const double x, double *err);

void stable_cdf(StableDist *dist, const double* x, const unsigned int Nx,
                double *cdf, double *err);

/******************************************************************************/
/*   CDF integrad functions                                                   */
/******************************************************************************/

double stable_cdf_g(double theta, void *dist);

/******************************************************************************/
/*   CDF^{-1} (quantiles)                                                     */
/******************************************************************************/

double stable_q_point(StableDist * dist, const double q, double * err);
void   stable_q(StableDist *dist, const double* q, const unsigned int Nq,
                double * inv, double * err);

/************************************************************************
 ************************************************************************
 * Vectorial methods                                                    *
 ************************************************************************
 ************************************************************************/
/*
  Alternative non-parallelized methods of evaluation have been implemented.
  These methods exploit the fact that some calculations are shared between
  different points of evaluation when evaluating the PDF or CDF, so these
  calculations can be realized just once.

  The performance achieved is high, sometimes comparable with parallel
  methods when little precision is required. However, achievable precision
  with these methods is low and non desired behavior of the PDF and CDF
  evaluation is observed.
*/

/* Stable distribution structure for vectorial methods*/

typedef struct
  {
/* Parameters:
    0-parametrization describen in Nolan, 1997 is employed by default
        alpha : stability index
        beta : skewness parameter
        scale: scale parameter
        mu_0 : 0-parametrization location parameter
        mu_1 : correspondig 1-parametrization location parameter    */
    double alpha;
    double beta;
    double sigma;
    double mu_0;
    double mu_1;

    /* Particular cases indicator (Gauss, Cauchy, Levy distribution, alpha==1, etc.) */
    int ZONE;

    /* Precalculated values */
    double alphainvalpha1;     /* alpha/(alpha-1)*/
    double xi;               /* -beta*tan(alpha*pi/2)*/
    double theta0;           /* 1/alpha*atan(beta*(tan(alpha*pi/2))=atan(-xi)/alpha;*/
    double c1, c2_part, c3;  /* additive and multiplicative constants*/
    double k1;               /* cos(alpha*theta0)^(1/(alpha-1)) = (1+xi^2)^(-0.5/(alpha-1));*/
    double S;                /* (1+xi^2)^(1/(2*alpha));*/
    double Vbeta1;           /*pow(1/dist->alpha,dist->alphainvalpha1) *
                                   (dist->alpha-1)*pow(-cos(dist->alpha*PI_2),1/(dist->alpha-1))*/

     /* These ones change from point to point of evaluation */
    double theta0_; /* theta0_ = +-theta0 */
    double beta_;
    double *xxipow;  /* (x-xi)^(alpha/(alpha-1))*/

    gsl_integration_workspace * gslworkspace;
    gsl_rng * gslrand;
  }
StableDistV;

typedef struct
  {
    double (*ptr_funcion)(StableDist *dist, const double x, double *err);
    StableDist *dist;
    const double *x;
    int Nx;
    double *pdf;
    double *err;
  }
StableArgsPdf;

typedef struct
  {
    double (*ptr_funcion)(StableDist *dist, const double x, double *err);
    StableDist *dist;
    const double *x;
    int Nx;
    double *cdf;
    double *err;
  }
StableArgsCdf;

/************************************************************************
 ************************************************************************
 * Parameter estimation                                                 *
 ************************************************************************
 ************************************************************************/

/******************************************************************************/
/*        Parameter estimation structure                                      */
/******************************************************************************/
typedef struct
  {
    StableDist *dist;
    double *data;
    unsigned int length;
    double nu_c;
    double nu_z;
  }
stable_like_params;

/* Estimation functions */

void stable_fit_init(StableDist *dist, const double *data,
             const unsigned int length,  double *nu_c,double *nu_z);

int stable_fit_koutrouvelis(StableDist *dist, const double *data, const unsigned int length);

int stable_fit(StableDist *dist, const double *data, const unsigned int length);

int stable_fit_mle(StableDist *dist, const double *data, const unsigned int length);

int stable_fit_mle2d(StableDist *dist, const double *data, const unsigned int length);

int stable_fit_whole(StableDist *dist, const double *data, const unsigned int length);

/* Auxiliary functions */

gsl_complex stable_samplecharfunc_point(const double* x,
             const unsigned int N, double t);

void stable_samplecharfunc(const double* x, const unsigned int Nx,
             const double* t, const unsigned int Nt, gsl_complex * z);

void stable_fft(double *data, const unsigned int length, double * y);

double stable_loglikelihood(StableDist *dist, double *data, const unsigned int length);

//stable_like_params

int stable_fit_iter_whole(StableDist *dist, const double * data, const unsigned int length);

int stable_fit_iter(StableDist *dist, const double * data,
             const unsigned int length, const double nu_c, const double nu_z);

double stable_loglike_p(stable_like_params *params);

double stable_minusloglikelihood(const gsl_vector * theta, void * p);

/************************************************************************
 ************************************************************************
 * Random numbers generation                                            *
 ************************************************************************
 ************************************************************************/
void stable_rnd(StableDist *dist, double* rnd, unsigned int n);

double stable_rnd_point(StableDist *dist);

void stable_rnd_seed(StableDist * dist, unsigned long int s);
#endif //stable_H

max_stars_repo_name: TantraLabs/gonum

max_issues_repo_name: TantraLabs/gonum

max_forks_repo_name: TantraLabs/gonum