text:
/**
 * Various binning and weighting routines for aperture pixel tables.
 *
*/

#include <math.h>

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

#include "spc_trace_functions.h"
#include "aXe_utils.h"
#include "aXe_grism.h"
#include "spc_spc.h"
#include "spc_driz.h"
#include "aper_conf.h"
#include "spce_output.h"
#include "spc_optimum.h"

#define SQR(x) ((x)*(x))
#define MIN(x,y) (((x)<(y))?(x):(y))
#define MAX(x,y) (((x)>(y))?(x):(y))

/*
gsl_matrix *
create_weightimage(ap_pixel *ap_p, const beam actbeam,
		   const aperture_conf *conf, const double exptime,
		   const double sky_cps)
{
  drzstamp_dim  dimension;
  drzstamp     *modvar;
  gsl_matrix   *weights;

  dimension = get_resample_dims(ap_p, actbeam);

  if (!dimension.resolution)
    return get_default_weight();

  compute_model_ivar(ap_p, conf->rdnoise, exptime, sky_cps);

  modvar = compute_modvar(ap_p, actbeam, dimension);

  weights = comp_allweight(modvar);

  return weights;
}
*/

/**
 * Function: compute_modvar
 *
 *
 * Parameters:
 * @param ap_p      - the PET table
 * @param actbeam   - beam to compute the modvar structure for
 * @param dimension - dimenasion of the images
 * @param modvar    - the structure to fill
 */
drzstamp *
compute_modvar(ap_pixel *ap_p, const beam actbeam,
	       const drzstamp_dim dimension)
{
  ap_pixel *cur_p;
  ap_pixel *tmp_p;

  quadrangle quad;

  double jacob;

  int jcen, icen;
  int jupp, iupp;
  int jlow, ilow;

  int ii, jj;
  int stpi, stpj;

  double maxarr, arr;
  //int iim, jjm;
  double value, allweig, weig;
  double stpc;
  double totweight;
  drzstamp *modvar;

  // allocate memory for the result structure
  modvar = alloc_drzstamp(dimension);

  // allocate memory
  tmp_p = (ap_pixel *) malloc(sizeof(ap_pixel));

  // go over each pixel
  for (cur_p = ap_p; cur_p->p_x != -1; cur_p++)
    {
      // Skip this pixel if it was not actually used
      if (fabs(cur_p->dist)>actbeam.width+0.5)
        continue;

      // transfer values to the temporary pixel
      tmp_p->lambda  = cur_p->xi;
      tmp_p->dist    = cur_p->dist;
      tmp_p->dxs     = cur_p->dxs;
      tmp_p->dlambda = 1.0;

      // create the quadrangle for the current pixel
      quad = get_quad_from_pixel(tmp_p, actbeam.orient, dimension);

      // get the jacobian (well, easy here)
      // the term "cos(cur_p->dxs)" must be there
      // to correct the enlargement necessary
      // to cover the whole lambda-crossdispersion area!
      // NOT COMPLETELY understood
      jacob = cos(cur_p->dxs);

      // get the central pixel (icen, jcen) of the current PET-pixel
      icen = (int) floor(cur_p->xi   - dimension.xstart+.5);
      jcen = (int) floor(cur_p->dist - dimension.ystart+.5);

      // get the uper and lower extend of the quadrangle in x
      iupp = (int)floor(quad.xmax - (double)icen + 0.5)+1;
      ilow = (int)floor(quad.xmin - (double)icen + 0.5);

      // get the uper and lower extend of the quadrangle in x
      jupp = (int)floor(quad.ymax - (double)jcen + 0.5)+1;
      jlow = (int)floor(quad.ymin - (double)jcen + 0.5);

      maxarr=0.0;
      totweight = 0.0;
      // go over the extend in x
      for (ii=ilow;ii<iupp;ii++)
	{
	  // go over the extend in x
	  for (jj=jlow;jj<jupp;jj++)
	    {

	      // get the coordinates of the current output pixel
	      stpi = icen+ii;
	      stpj = jcen+jj;

	      // check whether the current output pixel is within
	      // the stamp image; continue if not
	      if ((stpi>=dimension.xsize)||(stpi<0)
		  ||(stpj>=dimension.ysize)||(stpj<0))
		continue;

	      // get the area which falls onto the current output pixel
	      arr = boxer(stpi,stpj,quad.x,quad.y);

	      if (arr > 0.0)
		{
		  // get the already existing counts and weights
		  stpc = gsl_matrix_get(modvar->counts,stpi,stpj);
		  weig = gsl_matrix_get(modvar->weight,stpi,stpj);

		  // compute the new, total weight of the current output pixel
		  allweig = weig + arr*cur_p->weight;

		  // do a weighted sum of the count value at the current
		  // output pixel
		  value = (stpc*weig + arr*cur_p->model*cur_p->weight*jacob)
		    / (allweig);

		  // store the new count value and the new weight
		  gsl_matrix_set(modvar->counts,stpi,stpj,value);
		  gsl_matrix_set(modvar->weight,stpi,stpj,allweig);
		}
	    }
	}
    }

  // freep the temporary PET pixel
  free(tmp_p);

  // return the result
  return modvar;
}

void
shift_tracelength(ap_pixel *ap_p, const double xi_shift)
{
  ap_pixel *cur_p;

  // go over each PET pixel
  for (cur_p = ap_p; cur_p->p_x != -1; cur_p++)
    {
      // apply a shift in trace distance
      cur_p->xi += xi_shift;
    }
}

/**
 * Function: prepare_inv_variance
 * The function prepares the inverse model variances
 * for a data set containing an object and a background PET.
 * This is done in subroutines which compute the
 * theoretical inverse pixel variance following
 * different methods for the object and the background PET.
 *
 * Parameters:
 * @param ap_p     - the object PET table
 * @param bg_p     - the background PET table
 * @param dobck    - integer indicating background subtraction
 * @param conf     - the configuration structure
 * @param exptime  - the exposure time
 * @param sky_cps  - the sky background
 */
void
prepare_inv_variance(ap_pixel *ap_p, ap_pixel *bg_p, const int dobck,
		     const aperture_conf *conf, const double exptime,
		     const double sky_cps, const double xi_shift)
{

  if (dobck)
    {
      // shift both, the object and the background PET
      shift_tracelength(ap_p, xi_shift);
      shift_tracelength(bg_p, xi_shift);

      // compute the inverse variance for the object + background PET
      compute_total_ivar(bg_p, bg_p, conf->rdnoise, exptime, sky_cps);
    }
  else
    {
      // shift both, the object PET
      shift_tracelength(ap_p, xi_shift);

      // compute the inverse variance for the object PET
      compute_object_ivar(ap_p, conf->rdnoise, exptime, sky_cps);
    }
}

/**
 * Function: compute_model_ivar
 * The function computes for each PET pixel the
 * associated inverse variance value. The input
 * for variance and inverse variance are the
 * model value, the contamination value, the
 * constant sky background value and the readnoise.
 * The computed inverse variance value is stored
 * in the weight entry of the PET.
 * The function also offers the possibility to
 * shift the trace distance by a fixed amount.
 *
 * Parameters:
 * @param ap_p     - the PET table
 * @param rdnoise  - the readnoise value
 * @param exptime  - the exposure time
 * @param sky_cps  - the sky background
 * @param xi_shift - shift in trace distance
 */
void
compute_object_ivar(ap_pixel *ap_p, const double rdnoise,
		    const double exptime, const double sky_cps)
{
  double variance;
  double sqr_rdnoise;
  double sqr_exptime;
  ap_pixel *cur_p;

  // square the readnoise
  sqr_rdnoise = rdnoise*rdnoise;

  // square the exposure time
  sqr_exptime = exptime*exptime;

  // go over each PET pixel
  for (cur_p = ap_p; cur_p->p_x != -1; cur_p++)
    {
	// compute the variance
	variance = ((cur_p->model + cur_p->contam + sky_cps) * exptime
		    + sqr_rdnoise) / sqr_exptime;

      // store the inverse variance
      cur_p->weight = 1.0/variance;
    }
}

/**
 * Function: compute_value_ivar
 * The function computes for each PET pixel the
 * associated inverse variance value. The input
 * to derive the variance is the count value and
 * the readnoise.
 * The computed inverse variance value is stored
 * in the weight entry of the PET.
 * The function also offers the possibility to
 * shift the trace distance by a fixed amount.
 *
 * Parameters:
 * @param ap_p    - the PET table
 * @param rdnoise  - the readnoise value
 * @param exptime - the exposure time
 * @param sky_cps - the sky background
 * @param xi_shift - shift in trace distance
 */
void
compute_total_ivar(ap_pixel *ap_p, const ap_pixel *bg_p, const double rdnoise,
		   const double exptime, const double sky_cps)
{
  double variance;
  double sqr_rdnoise;
  double sqr_exptime;

  ap_pixel *cur_p;
  const ap_pixel *bac_p;

  // square the readnoise
  sqr_rdnoise = rdnoise*rdnoise;

  // square the exposure time
  sqr_exptime = exptime*exptime;

  // go over each PET pixel
  bac_p = bg_p;
  for (cur_p = ap_p; cur_p->p_x != -1; cur_p++)
    {
      // check that the foreground and background
      // PET element describe the same pixel
      if (bac_p->x != cur_p->x || bac_p->y != cur_p->y)
	aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		     "aXe_PET2SPC: Background PET and Object PET "
		     "have different pixel orders in PET's.\n");
      else
      // compute the variance
      variance = ((cur_p->count + bac_p->count + sky_cps) * exptime + sqr_rdnoise) / sqr_exptime;

      // store the inverse variance
      cur_p->weight = 1.0/variance;

      // increase the background pointer
      bac_p++;
    }
}

/**
 * Function: alloc_drzstamp
 * Allocates memory for a drizzle stamp structure
 * according the the dimension specified in the input.
 *
 * Parameters:
 * @param dimension - the dimesions to allocate
 *
 * Returns:
 * @return ret  - the drizzle stamp image allocated
 */
drzstamp *
alloc_drzstamp(const drzstamp_dim dimension)
{
  drzstamp   *ret;
  gsl_matrix *counts;
  gsl_matrix *weight;

  ret = (drzstamp *) malloc(sizeof(drzstamp));

  /* Allocate the stamp matrix*/
  /* Fill stamp with NaN values */
  counts = gsl_matrix_alloc(dimension.xsize,dimension.ysize);
  //  gsl_matrix_set_all(counts, GSL_NAN);
  gsl_matrix_set_all(counts, 0.0);

  /* Allocate the weight matrix*/
  /* Fill weight with 0.0 values */
  weight = gsl_matrix_alloc(dimension.xsize,dimension.ysize);
  gsl_matrix_set_all(weight, 0.0);

  // fill the output structure
  ret->counts = counts;
  ret->weight = weight;

  // return the allocated structure
  return ret;
}

/**
 * Function: get_all_dims
 * The function derives the dimensional information
 * on the foreground and, if available, on the
 * background PET. The dimensions are compared, and
 * in case of inequality an error is thrown.
 * If the dimensions are equal, one of them is
 * returned.
 *
 * Parameters:
 * @param ap_p    - the table of aperture pixels
 * @param bg_p    - the table of background pixels
 * @param actbeam - the beam to compute the dimensions for
 * @param dobck   - integer indicating background subtraction
 *
 * Returns:
 * @return ret  - the dimension structure
 */
drzstamp_dim
get_all_dims(const ap_pixel *ap_p, const ap_pixel *bg_p, const beam actbeam,
	     const int dobck)
{
  drzstamp_dim ret;
  drzstamp_dim bck;

  // get the aperture PET diemnsion
  ret = get_resample_dims(ap_p, actbeam);

  if (dobck)
    {
      // get the background PET dimension
      bck = get_resample_dims(bg_p, actbeam);

      // check whether the two dimensions
      // are equal, throw an error if not
      if (ret.resolution != bck.resolution
	  || ret.xstart  != bck.xstart
	  || ret.ystart  != bck.ystart
	  || ret.xsize   != bck.xsize
	  || ret.ysize   != bck.ysize)
	aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		     "aXe_PET2SPC: Background PET and Object PET "
		     "do not have the same dimension.\n");
    }

  // return one of the dimensional infos
  return ret;
}

/**
 * Function: get_resample_dims
 * The function parses through a PET table and determines
 * the dimension of the image in the coordinates
 * trace distance - cross-dispersion direction. Those
 * dimensions are stored in a special structure and
 * returned to the calling routine.
 *
 * Parameters:
 * @param ap_p    - the table of aperture pixels
 * @param actbeam - the beam to compute the dimensions for
 *
 * Returns:
 * @return ret  - the dimension structure
 */
drzstamp_dim
get_resample_dims(const ap_pixel *ap_p, const beam actbeam)
{
  drzstamp_dim ret;

  const ap_pixel *cur_p;

  int i_xint, i_dist;
  int imin_xint=0, imin_dist=0;
  int imax_xint=0, imax_dist=0;

  int npixel = 0;

  // immediately return empty PET's
  if (ap_p==NULL)
    return get_default_dim();

  // set the tmp pixel to the table start
  cur_p = ap_p;

  // initialize some relevant numbers
  // based on the first pixel
  i_xint = floor(cur_p->xi   + 0.5);
  i_dist = floor(cur_p->dist + 0.5);
  imin_xint = i_xint;
  imax_xint = i_xint;
  imin_dist = i_dist;
  imax_dist = i_dist;

  // go over all pixels, check
  // for new MIN's ans MAX's
  for (cur_p = ap_p; cur_p->p_x != -1; cur_p++)
    {
      // Skip the pixel if it was not actually used
      if (fabs(cur_p->dist)>actbeam.width+.5)
	continue;

      // compute the rounded value in (x,y)
      i_xint = floor(cur_p->xi   + 0.5);
      i_dist = floor(cur_p->dist + 0.5);

      // check for extrema in x
      imin_xint = MIN(imin_xint, i_xint);
      imax_xint = MAX(imax_xint, i_xint);

      // check for extrema in y
      imin_dist = MIN(imin_dist, i_dist);
      imax_dist = MAX(imax_dist, i_dist);

      // enhance the pixel counter
      npixel++;
    }

  // in case there were no valid pixels,
  // return a dummy structure
  if (!npixel)
    return get_default_dim();

  // fill the structure with the correct values
  ret.resolution = 1.0;
  ret.xstart     = imin_xint - NSPARE_PIX;
  ret.ystart     = imin_dist - NSPARE_PIX;
  ret.xsize      = imax_xint - imin_xint + 2 * NSPARE_PIX;
  ret.ysize      = imax_dist - imin_dist + 2 * NSPARE_PIX;

  // return the structure
  return ret;
}

/**
 * Function: get_default_dim
 * The function creates and returns the default
 * dimension structure. The dewfault dimension
 * structure contains a 0.0 in all fields.
 *
 * Returns:
 * @return res  - the dimension structure
 */
drzstamp_dim
get_default_dim()
{
  drzstamp_dim ret;

  // just set everything to zero
  ret.resolution=0.0;
  ret.xstart = 0;
  ret.ystart = 0;
  ret.xsize  = 0;
  ret.ysize  = 0;

  // return the result
  return ret;
}

/**
 * Function: get_default_modvar
 * The function computes the default drizzle stamp
 * structure, which in this case is used for
 * tha calculation of the resampled variance and
 * profile image.
 *
 * Returns:
 * @return res  - the resulting drzstamp structure
 */
drzstamp *
get_default_modvar()
{
  gsl_matrix *counts;
  gsl_matrix *weight;
  drzstamp *res;

  res = (drzstamp *) malloc(sizeof(drzstamp));

  /* Create a dummy stamp image */
  counts = gsl_matrix_alloc(10,10);
  weight = gsl_matrix_alloc(10,10);
  /* Fill stamp with 0.0 values */
  gsl_matrix_set_all(counts, 0.0);
  gsl_matrix_set_all(weight, 0.0);

  // fill the output
  res->counts = counts;
  res->weight = weight;

  // return the output
  return res;
}

gsl_matrix *
get_default_weight()
{
  gsl_matrix *res;

  res = gsl_matrix_alloc(10,10);
  gsl_matrix_set_all(res, 0.0);

  // return the output
  return res;
}

/*
 * Function: comp_opt_weight
 * The function uses the model map and and the variance
 * map and computes optimale weigths according to the
 * Hoorne method.
 * Pixels outside of the extraction
 * area are set to 1000.0
 *
 * Parameters:
 * @param mod_map - the matrix with the model values
 * @param var_map - the matrix with the variance values
 * @param ob      - the object structure
 *
 * Returns:
 * @return weight - the matrix with the exposure time weights
 */
gsl_matrix *
comp_allweight(drzstamp *modvar)
{
  gsl_matrix *weight;

  double mod_sum, weight_sum;
  double contr, norm, allweight;
  double mod_val;
  double act_weight=0;

  //int beamInt = 0;
  int i, j;

  // allocate the weight matrix and set the default
  weight = gsl_matrix_alloc(modvar->counts->size1, modvar->counts->size2);
  gsl_matrix_set_all(weight, 1000.0);

  //* go over all columns
  for (i=0; i < (int)modvar->counts->size1; i++)
    {
      mod_sum = 0.0;
      contr = 0.0;
      allweight = 0.0;
      weight_sum=0.0;

      // determine for each column the total model counts
      // and the number of pixels with non-zero model_counts
      for (j=0; j < (int)modvar->counts->size2; j++)
	{
	  mod_sum = mod_sum + gsl_matrix_get(modvar->counts, i, j);
	  contr = contr + 1.0;
	}

      // check whether the column has model values.
      // normalize the model values and compute
      // optimal weights if yes.
      if (mod_sum > 0.0)
	{
	  //* determine the mean model value
	  //	  norm = mod_sum / contr;
	  norm = mod_sum;

	  // go over each row
	  for (j=0; j < (int)modvar->counts->size2; j++)
	    {
	      // normalize the model counts
	      mod_val = gsl_matrix_get(modvar->counts, i, j)/norm;

	      // store the normalized model counts
	      gsl_matrix_set(modvar->counts, i, j, mod_val);

	      // add up the normalization value for the optimal weights
	      weight_sum = weight_sum
		+ mod_val*mod_val*gsl_matrix_get(modvar->weight, i, j);
	    }

	  // finally compute and write the weights:
	  // go over each pixel
	  for (j=0; j < (int)modvar->counts->size2; j++)
	    {
	      if (gsl_matrix_get(modvar->counts, i, j) > 0.0 )
		{
		  // compute and set the individual pixel weight
		  act_weight = gsl_matrix_get(modvar->counts,i,j)
		    *gsl_matrix_get(modvar->weight,i,j)/weight_sum;
		  //		  gsl_matrix_set(weight,i,j,gsl_matrix_get(modvar->counts,i,j)
		  //		 *gsl_matrix_get(modvar->weight,i,j)/weight_sum);
		  gsl_matrix_set(weight, i, j, act_weight);
		  allweight = allweight + gsl_matrix_get(weight, i, j);
		}
	      else
		{
		  // set the default value if no model value is zero
		  gsl_matrix_set(weight, i, j,0.0);
		}
	    }
	}
      // if the column does not have model values at all:
      else
	{
	  // go over each row
	  for (j=0; j < (int)modvar->counts->size2; j++)
	    {
	      // set the inside default value
	      gsl_matrix_set(weight, i,j,1.0);
	      allweight = allweight + gsl_matrix_get(modvar->counts, i, j);

	    }
	}
    }

  // return the weight matrix
  return weight;
}

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe


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

#ifndef JROBINSON_H
#define JROBINSON_H

#define JROBINSON_EPS (0.0000001)

int jrobinson( const gsl_matrix *A, double epsilon, gsl_vector *p1, gsl_vector *p2, int &iter );

#endif

ext: h


text:
#include <math.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_eigen.h>
#include "klp_matrix_initializers.h"
#include "klp_matrix_functions.h"
#include "shared/constants.h"

TRANSITION_MATRIX transpose_matrix(TRANSITION_MATRIX matrix) {
  int i, j;
  TRANSITION_MATRIX transposed_matrix = init_transition_matrix(matrix.row_length, matrix.type);

  for (i = 0; i < matrix.row_length; ++i) {
    for (j = 0; j < matrix.row_length; ++j) {
      T_COL_ORDER(transposed_matrix, i, j) = T_ROW_ORDER(matrix, i, j);
    }
  }

  free_transition_matrix(matrix);
  return transposed_matrix;
}

TRANSITION_MATRIX inverse(TRANSITION_MATRIX transition_matrix) {
  int i, j, signum;
  gsl_matrix* matrix_to_invert = gsl_matrix_alloc(transition_matrix.row_length, transition_matrix.row_length);
  gsl_matrix* inversion_matrix = gsl_matrix_alloc(transition_matrix.row_length, transition_matrix.row_length);
  gsl_permutation* permutation = gsl_permutation_alloc(transition_matrix.row_length);

  for (i = 0; i < transition_matrix.row_length; ++i) {
    for (j = 0; j < transition_matrix.row_length; ++j) {
      gsl_matrix_set(matrix_to_invert, i, j, T_ROW_ORDER(transition_matrix, i, j));
    }
  }

  gsl_linalg_LU_decomp(matrix_to_invert, permutation, &signum);
  gsl_linalg_LU_invert(matrix_to_invert, permutation, inversion_matrix);

  for (i = 0; i < transition_matrix.row_length; ++i) {
    for (j = 0; j < transition_matrix.row_length; ++j) {
      T_ROW_ORDER(transition_matrix, i, j) = gsl_matrix_get(inversion_matrix, i, j);
    }
  }

  gsl_matrix_free(matrix_to_invert);
  gsl_matrix_free(inversion_matrix);
  gsl_permutation_free(permutation);

  return transition_matrix;
}

TRANSITION_MATRIX convert_klp_matrix_to_transition_matrix(KLP_MATRIX* klp_matrix, KLP_PARAMS* klp_params) {
  int resolved;
  double* number_of_adjacent_moves;
  transition_probability probability_function = NULL;

  resolved = find_start_and_end_positions_in_klp_matrix(klp_matrix, klp_params);

  if (klp_params->max_dist) {
    if (!klp_params->bp_dist) {
      set_bp_dist_from_start_and_end_positions(*klp_matrix, klp_params, resolved);
    }

    extend_klp_matrix_to_all_possible_positions(klp_matrix, *klp_params);
    populate_remaining_probabilities_in_klp_matrix(klp_matrix, *klp_params);

    if (resolved != 2) {
      find_start_and_end_positions_in_klp_matrix(klp_matrix, klp_params);
    }
  }

  number_of_adjacent_moves = populate_number_of_adjacent_moves(*klp_matrix, *klp_params);

#ifdef DEBUG
  int i;
  printf("\nFull dataset:\n");

  for (i = 0; i < klp_matrix->length; ++i) {
    printf("%d\t%d\t%f\t%d possible move(s)\n", klp_matrix->k[i], klp_matrix->l[i], klp_matrix->p[i], (int)number_of_adjacent_moves[i]);
  }

  printf("\n");
#endif

  switch (10 * klp_params->hastings + klp_params->energy_based) {
    case 0:
      probability_function = &transition_rate_from_probabilities;
#ifdef DEBUG
      printf("probability_function: transition_rate_from_probabilities\n");
#endif
      break;

    case 1:
      probability_function = &transition_rate_from_energies;
#ifdef DEBUG
      printf("probability_function: transition_rate_from_energies\n");
#endif
      break;

    case 10:
      probability_function = &transition_rate_from_probabilities_with_hastings;
#ifdef DEBUG
      printf("probability_function: transition_rate_from_probabilities_with_hastings\n");
#endif
      break;

    case 11:
      probability_function = &transition_rate_from_energies_with_hastings;
#ifdef DEBUG
      printf("probability_function: transition_rate_from_energies_with_hastings\n");
#endif
      break;
  }

  return populate_transition_matrix_from_stationary_matrix(*klp_matrix, *klp_params, number_of_adjacent_moves, probability_function);
}

int find_start_and_end_positions_in_klp_matrix(KLP_MATRIX* klp_matrix, KLP_PARAMS* klp_params) {
  int i, resolved = 0;

  if (klp_params->start_state == -1) {
    for (i = 0; i < klp_matrix->length && klp_params->start_state == -1; ++i) {
      if (klp_matrix->k[i] == 0) {
        klp_params->start_state = i;
        resolved++;
      }
    }
  } else {
    resolved++;
  }

  if (klp_params->end_state == -1) {
    for (i = 0; i < klp_matrix->length && klp_params->end_state == -1; ++i) {
      if (klp_matrix->l[i] == 0) {
        klp_params->end_state = i;
        resolved++;
      }
    }
  } else {
    resolved++;
  }

#ifdef DEBUG
  printf("\nstart_index:\t%d\n", klp_params->start_state);
  printf("end_index:\t%d\n", klp_params->end_state);
  printf("resolved:\t%d\n", resolved);
#endif
  return resolved;
}

void set_bp_dist_from_start_and_end_positions(const KLP_MATRIX klp_matrix, KLP_PARAMS* klp_params, int resolved) {
  int distance_from_start, distance_from_end;

  distance_from_start = distance_from_end = -1;

  if (klp_params->start_state > 0) {
    distance_from_start = klp_matrix.l[klp_params->start_state];
  }

  if (klp_params->end_state > 0) {
    distance_from_end = klp_matrix.k[klp_params->end_state];
  }

  if (distance_from_start == distance_from_end && resolved) {
    klp_params->bp_dist = distance_from_start;
  } else if (distance_from_start >= 0 && distance_from_end == -1) {
    klp_params->bp_dist = distance_from_start;
  } else if (distance_from_end >= 0 && distance_from_start == -1) {
    klp_params->bp_dist = distance_from_end;
  } else {
    fprintf(stderr, "Can't infer the input structure distances for the energy grid. We found (0, %d) and (%d, 0). Consider using the -d flag to manually set the base pair distance between the two structures.\n", distance_from_end, distance_from_start);
    printf("-3\n");
    exit(0);
  }

#ifdef DEBUG
  printf("bp_dist:\t%d\n", klp_params->bp_dist);
#endif
}

void extend_klp_matrix_to_all_possible_positions(KLP_MATRIX* klp_matrix, const KLP_PARAMS klp_params) {
  int i, j, m, position_in_input_data, pointer, valid_positions = 0;

#ifdef DEBUG
  printf("\nAccessible positions (top-left is [0, 0]):\n");
#endif

  for (i = 0; i <= klp_params.max_dist; ++i) {
    for (j = 0; j <= klp_params.max_dist; ++j) {
      if (
        i + j >= klp_params.bp_dist &&
        i + klp_params.bp_dist >= j &&
        j + klp_params.bp_dist >= i &&
        (i + j) % 2 == klp_params.bp_dist % 2
      ) {
#ifdef DEBUG
        position_in_input_data = -1;

        for (m = 0; m < klp_matrix->length && position_in_input_data == -1; ++m) {
          if (klp_matrix->k[m] == i && klp_matrix->l[m] == j) {
            position_in_input_data = m;
          }
        }

        printf(position_in_input_data == -1 ? "X" : "O");
#endif
        valid_positions++;
      } else {
#ifdef DEBUG
        printf(" ");
#endif
      }
    }

#ifdef DEBUG
    printf("\n");
#endif
  }

  klp_matrix->k = realloc(klp_matrix->k, valid_positions * sizeof(int));
  klp_matrix->l = realloc(klp_matrix->l, valid_positions * sizeof(int));
  klp_matrix->p = realloc(klp_matrix->p, valid_positions * sizeof(double));

  pointer = klp_matrix->length;

#ifdef DEBUG
  printf("\nInput dataset:\n");

  for (i = 0; i < klp_matrix->length; ++i) {
    printf("%d\t%d\t%d\t%f\n", i, klp_matrix->k[i], klp_matrix->l[i], klp_matrix->p[i]);
  }

#endif

  for (i = 0; i <= klp_params.max_dist; ++i) {
    for (j = 0; j <= klp_params.max_dist; ++j) {
      if (
        i + j >= klp_params.bp_dist &&
        i + klp_params.bp_dist >= j &&
        j + klp_params.bp_dist >= i &&
        (i + j) % 2 == klp_params.bp_dist % 2
      ) {
        position_in_input_data = -1;

        for (m = 0; m < klp_matrix->length && position_in_input_data == -1; ++m) {
          if (klp_matrix->k[m] == i && klp_matrix->l[m] == j) {
            position_in_input_data = m;
          }
        }

        if (position_in_input_data < 0) {
          klp_matrix->k[pointer] = i;
          klp_matrix->l[pointer] = j;
          klp_matrix->p[pointer] = 0.;
          pointer++;
        }
      }
    }
  }

  klp_matrix->length = valid_positions;
}

void populate_remaining_probabilities_in_klp_matrix(KLP_MATRIX* klp_matrix, const KLP_PARAMS klp_params) {
  int i;
  double epsilon_per_cell;

  if (klp_params.epsilon) {
    // Extend the energy grid by adding an epsilon value to all 0-probability positions.
    epsilon_per_cell = klp_params.epsilon / klp_matrix->length;

    for (i = 0; i < klp_matrix->length; ++i) {
      if (klp_matrix->p[i] > 0) {
        klp_matrix->p[i] = (klp_matrix->p[i] + epsilon_per_cell) / (1. + klp_params.epsilon);
      } else {
        klp_matrix->p[i] = epsilon_per_cell / (1. + klp_params.epsilon);
      }
    }
  }
}

double* populate_number_of_adjacent_moves(const KLP_MATRIX klp_matrix, const KLP_PARAMS klp_params) {
  int i;
  double* number_of_adjacent_moves;

  number_of_adjacent_moves = malloc(klp_matrix.length * sizeof(double));

  for (i = 0; i < klp_matrix.length; ++i) {
    number_of_adjacent_moves[i] = RUN_TYPE(klp_params.run_type, DIAG_MOVES_ONLY_FLAG) ? (double)number_of_permissible_single_bp_moves(klp_matrix, i) : (double)(klp_matrix.length - 1);
  }

  return number_of_adjacent_moves;
}

int number_of_permissible_single_bp_moves(const KLP_MATRIX klp_matrix, int i) {
  int j, x, y, a, b, num_moves = 0;

  x = klp_matrix.k[i];
  y = klp_matrix.l[i];

  for (j = 0; j < klp_matrix.length; ++j) {
    a = klp_matrix.k[j];
    b = klp_matrix.l[j];

    if (
      // Because N(x, y) is restricted to entries in *k and *l, we *assume* the input data satisfies the triangle inequality and bounds.
      (int)abs(x - a) == 1 && (int)abs(y - b) == 1
    ) {
      num_moves++;
    }
  }

  return num_moves;
}

TRANSITION_MATRIX populate_transition_matrix_from_stationary_matrix(const KLP_MATRIX klp_matrix, const KLP_PARAMS klp_params, const double* number_of_adjacent_moves, transition_probability probability_function) {
  int i, j;
  double row_sum;
  TRANSITION_MATRIX transition_matrix;

  transition_matrix = init_transition_matrix(klp_matrix.length, MATRIX_TYPE(klp_params.rate_matrix));

  for (i = 0; i < transition_matrix.row_length; ++i) {
    row_sum = 0.;

    for (j = 0; j < transition_matrix.row_length; ++j) {
      if (i != j) {
        if (RUN_TYPE(klp_params.run_type, FULLY_CONNECTED_FLAG) || (RUN_TYPE(klp_params.run_type, DIAG_MOVES_ONLY_FLAG) && ONE_BP_MOVE(i, j))) {
          if (NONZERO_TO_NONZERO_PROB(i, j)) {
            T_ROW_ORDER(transition_matrix, i, j) = \
                                                   probability_function(klp_matrix, number_of_adjacent_moves, i, j, klp_params.rate_matrix);
          }
        }

        row_sum += T_ROW_ORDER(transition_matrix, i, j);
      }
    }

    T_ROW_ORDER(transition_matrix, i, i) = klp_params.rate_matrix ? -row_sum : 1 - row_sum;
  }

  return transition_matrix;
}

double transition_rate_from_probabilities(const KLP_MATRIX klp_matrix, const double* number_of_adjacent_moves, int i, int j, short rate_matrix) {
  if (rate_matrix) {
    return MIN(1., klp_matrix.p[j] / klp_matrix.p[i]);
  } else {
    return MIN(1., klp_matrix.p[j] / klp_matrix.p[i]) / number_of_adjacent_moves[i];
  }
}

double transition_rate_from_energies(const KLP_MATRIX klp_matrix, const double* number_of_adjacent_moves, int i, int j, short rate_matrix) {
  if (rate_matrix) {
    return MIN(1., exp(-(klp_matrix.p[j] - klp_matrix.p[i]) / RT));
  } else {
    return MIN(1., exp(-(klp_matrix.p[j] - klp_matrix.p[i]) / RT)) / number_of_adjacent_moves[i];
  }
}

double transition_rate_from_probabilities_with_hastings(const KLP_MATRIX klp_matrix, const double* number_of_adjacent_moves, int i, int j, short rate_matrix) {
  if (rate_matrix) {
    return MIN(1., (number_of_adjacent_moves[i] / number_of_adjacent_moves[j]) * (klp_matrix.p[j] / klp_matrix.p[i]));
  } else {
    return MIN(1., (number_of_adjacent_moves[i] / number_of_adjacent_moves[j]) * (klp_matrix.p[j] / klp_matrix.p[i])) / number_of_adjacent_moves[i];
  }
}

double transition_rate_from_energies_with_hastings(const KLP_MATRIX klp_matrix, const double* number_of_adjacent_moves, int i, int j, short rate_matrix) {
  if (rate_matrix) {
    return MIN(1., (number_of_adjacent_moves[i] / number_of_adjacent_moves[j]) * exp(-(klp_matrix.p[j] - klp_matrix.p[i]) / RT));
  } else {
    return MIN(1., (number_of_adjacent_moves[i] / number_of_adjacent_moves[j]) * exp(-(klp_matrix.p[j] - klp_matrix.p[i]) / RT)) / number_of_adjacent_moves[i];
  }
}

ext: c

max_stars_repo_name: evansenter/hermes

max_issues_repo_name: evansenter/hermes

max_forks_repo_name: evansenter/hermes


text:
#pragma once

#include "ui_MaterialVisualizationWidget.h"

#include <mitkDataNode.h>
#include <vtkDataObject.h>

#include <array>
#include <gsl.h>

namespace mitk
{
class DataNode;
class TransferFunction;
}

/*! \brief   A widget for editing transfer function. */
class MaterialVisualizationWidget : public QWidget
{
    Q_OBJECT
public:
    MaterialVisualizationWidget(QWidget* parent = nullptr, Qt::WindowFlags f = 0);
    ~MaterialVisualizationWidget();

    /*!
     * \brief   Sets the data node whose transfer function is to be edited.
     *
     * \param   node                The data node.
     * \param   dataSet             The vtkDataSet stored in the data node (the access may differ for mitk::Surface and mitk::UnstructuredGrid).
     * \param   tfPropertyName      Name of the transfer function property used by the renderer.
     * \param   attributeType       Type of the attribute to be shown (vtkDataObject::POINT or vtkDataObject::CELL).
     * \param   propertyStorageNode If non-null, the node where to store various transfer functions to be re-used later. If null, the 'node' will be used.
     */
    void setDataNode(mitk::DataNode* node, vtkDataSet* dataSet, gsl::cstring_span<> tfPropertyName,
                     vtkDataObject::AttributeTypes attributeType = vtkDataObject::POINT,
                     mitk::DataNode* propertyStorageNode = nullptr);
    /*! \brief   Resets the data node to null. */
    void resetDataNode() { setDataNode(nullptr, nullptr, ""); }

    /*! \brief   Gets data node. */
    mitk::DataNode* getDataNode() const { return _node; }

private slots:
    void _setCurrentArray(const QString& arrayName);
    void _setCurrentComponent(int comboBoxIndex);
    void _resetTF();
    void _setNanColor(QColor c);

    void _setColorSpace(int colorSpace);
    void _loadTF();
    void _saveTF();

private:
    void _dataNodeDeleted();
    std::string _getCurrentTFPropName();
    void _setupDataComboBoxes();
    void _setTransferFunction(bool reset = false);
    static void _setupTFProperties(mitk::TransferFunction* tf, int componentIndex);
    static void _generateTF(mitk::TransferFunction* tf, const std::array<double, 2>& range);
    static void _rescaleTF(mitk::TransferFunction* tf, const std::array<double, 2>& range);
    QColor _getNanColor();
    void _addNodeObserver();
    void _removeNodeObserver();

    mitk::DataNode* _node = nullptr;
    vtkDataSet* _dataSet = nullptr;
    std::string _tfPropertyName;
    vtkDataObject::AttributeTypes _attributeType = vtkDataObject::POINT;
    mitk::DataNode* _propertyStorageNode = nullptr;
    Ui::MaterialVisualizationWidget _UI;
    unsigned long _nodeObserverTag = -1;
};

max_stars_repo_name: carthurs/CRIMSONGUI

max_issues_repo_name: carthurs/CRIMSONGUI

max_forks_repo_name: carthurs/CRIMSONGUI


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

#include <complex.h>
#include <fftw3.h>

#include <gsl/gsl_sf_result.h>
#include <gsl/gsl_sf_gamma.h>
#include "ccl.h"

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

This is the famous FFTLog.

First imlplemented by the living legend Andrew Hamilton:

http://casa.colorado.edu/~ajsh/FFTLog/

This version is a C version that was adapted from the C++ version found
in Copter JWG Carlson, another big loss for the cosmology community.

https://github.com/jwgcarlson/Copter

I've transformed this from C++ to C99 as the lowest common denominator
and provided bindings for C++ and python.

These are the C++ bindings

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

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
#ifndef M_LN2
#define M_LN2 0.69314718056
#endif
/* This code is FFTLog, which is described in arXiv:astro-ph/9905191 */

static double complex lngamma_fftlog(double complex z)
{
  gsl_sf_result lnr, phi;
  gsl_sf_lngamma_complex_e(creal(z), cimag(z), &lnr, &phi);
  return lnr.val + I*phi.val;
}

static double complex gamma_fftlog(double complex z)
{
  return cexp(lngamma_fftlog(z));
}

static double complex polar (double r, double phi)
{
  return (r*cos(phi) +I*(r*sin(phi)));
}

static void lngamma_4(double x, double y, double* lnr, double* arg)
{
  double complex w = lngamma_fftlog(x+y*I);
  if(lnr) *lnr = creal(w);
  if(arg) *arg = cimag(w);
}

static double goodkr(int N, double mu, double q, double L, double kr)
{
  double xp = (mu+1+q)/2;
  double xm = (mu+1-q)/2;
  double y = M_PI*N/(2*L);
  double lnr, argm, argp;
  lngamma_4(xp, y, &lnr, &argp);
  lngamma_4(xm, y, &lnr, &argm);
  double arg = log(2/kr) * N/L + (argp + argm)/M_PI;
  double iarg = round(arg);
  if(arg != iarg)
    kr *= exp((arg - iarg)*L/N);
  return kr;
}

/* Pre-compute the coefficients that appear in the FFTLog implementation of
 * the discrete Hankel transform.  The parameters N, mu, and q here are the
 * same as for the function fht().  The parameter L is defined (for whatever
 * reason) to be N times the logarithmic spacing of the input array, i.e.
 *   L = N * log(r[N-1]/r[0])/(N-1) */
static void compute_u_coefficients(int N, double mu, double q, double L, double kcrc, double complex *u)
{
  double y = M_PI/L;
  double k0r0 = kcrc * exp(-L);
  double t = -2*y*log(k0r0/2);

  if(q == 0) {
    double x = (mu+1)/2;
    double lnr, phi;
    for(int m = 0; m <= N/2; m++) {
      lngamma_4(x, m*y, &lnr, &phi);
      u[m] = polar(1.0,m*t + 2*phi);
    }
  }
  else {
    double xp = (mu+1+q)/2;
    double xm = (mu+1-q)/2;
    double lnrp, phip, lnrm, phim;
    for(int m = 0; m <= N/2; m++) {
      lngamma_4(xp, m*y, &lnrp, &phip);
      lngamma_4(xm,-m*y, &lnrm, &phim);
      u[m] = polar(exp(q*M_LN2 + lnrp - lnrm), m*t + phip - phim);
    }
  }

  for(int m = N/2+1; m < N; m++)
    u[m] = conj(u[N-m]);
  if((N % 2) == 0)
    u[N/2] = (creal(u[N/2]) + I*0.0);
}

/* Compute the discrete Hankel transform of the function a(r).  See the FFTLog
 * documentation (or the Fortran routine of the same name in the FFTLog
 * sources) for a description of exactly what this function computes.
 * If u is NULL, the transform coefficients will be computed anew and discarded
 * afterwards.  If you plan on performing many consecutive transforms, it is
 * more efficient to pre-compute the u coefficients. */
static void fht(int npk, int N,
		double *k, double **pk,
		double *r, double **xi,
		double dim, double mu, double q, double kcrc,
		int noring, double complex* u, int *status)
{
  fftw_plan forward_plan, reverse_plan;
  double L = log(k[N-1]/k[0]) * N/(N-1.);
  double complex* ulocal = NULL;
  if(u == NULL) {
    if(noring)
      kcrc = goodkr(N, mu, q, L, kcrc);

    ulocal = malloc (sizeof(complex double)*N);
    if(ulocal==NULL)
      *status=CCL_ERROR_MEMORY;

    if(*status == 0) {
      compute_u_coefficients(N, mu, q, L, kcrc, ulocal);
      u = ulocal;
    }
  }

  fftw_complex* a_tmp;
  fftw_complex* b_tmp;
  if(*status == 0) {
    a_tmp = fftw_alloc_complex(N);
    if(a_tmp==NULL)
      *status=CCL_ERROR_MEMORY;
  }
  if(*status == 0) {
    b_tmp = fftw_alloc_complex(N);
    if(b_tmp==NULL)
      *status=CCL_ERROR_MEMORY;
  }

  if(*status == 0) {
    /* Compute the convolution b = a*u using FFTs */
    forward_plan = fftw_plan_dft_1d(N,
                                    (fftw_complex*) a_tmp,
                                    (fftw_complex*) b_tmp,
                                    -1, FFTW_ESTIMATE);
    reverse_plan = fftw_plan_dft_1d(N,
                                    (fftw_complex*) b_tmp,
                                    (fftw_complex*) b_tmp,
                                    +1, FFTW_ESTIMATE);
  }

  if(*status == 0) {
    #pragma omp parallel default(none) \
                         shared(npk, N, k, pk, r, xi, \
                                dim, mu, q, kcrc, u, status, \
                                forward_plan, reverse_plan, \
                                L, ulocal)
    {
      int local_status = 0;

      double *prefac_pk=NULL;
      if(local_status == 0) {
        prefac_pk = malloc(N*sizeof(double));
        if(prefac_pk==NULL)
          local_status=CCL_ERROR_MEMORY;
      }

      double *prefac_xi=NULL;
      if(local_status == 0) {
        prefac_xi = malloc(N*sizeof(double));
        if(prefac_xi==NULL)
          local_status=CCL_ERROR_MEMORY;
      }

      fftw_complex* a=NULL;
      fftw_complex* b=NULL;
      if(local_status == 0) {
        a = fftw_alloc_complex(N);
        if(a==NULL)
        local_status=CCL_ERROR_MEMORY;
      }

      if(local_status == 0) {
        b = fftw_alloc_complex(N);
        if(b==NULL)
          local_status=CCL_ERROR_MEMORY;
      }

      if(local_status == 0) {
        for(int i = 0; i < N; i++)
          prefac_pk[i] = pow(k[i], dim/2-q);

        /* Compute k's corresponding to input r's */
        double k0r0 = kcrc * exp(-L);
        r[0] = k0r0/k[0];
        for(int n = 1; n < N; n++)
          r[n] = r[0] * exp(n*L/N);

        double one_over_2pi_dhalf = pow(2*M_PI,-dim/2);
        for(int i = 0; i < N; i++)
          prefac_xi[i] = one_over_2pi_dhalf * pow(r[i], -dim/2-q);

        #pragma omp for
        for(int j = 0; j < npk; j++) {
          for(int i = 0; i < N; i++)
            a[i] = prefac_pk[i] * pk[j][i];

          fftw_execute_dft(forward_plan,a,b);
          for(int m = 0; m < N; m++)
            b[m] *= u[m] / (double)(N);       // divide by N since FFTW doesn't normalize the inverse FFT
          fftw_execute_dft(reverse_plan,b,b);

          /* Reverse b array */
          double complex tmp;
          for(int n = 0; n < N/2; n++) {
            tmp = b[n];
            b[n] = b[N-n-1];
            b[N-n-1] = tmp;
          }

          for(int i = 0; i < N; i++)
            xi[j][i] = prefac_xi[i] * creal(b[i]);
        }
      }

      free(prefac_pk);
      free(prefac_xi);
      fftw_free(a);
      fftw_free(b);

      if (local_status) {
        #pragma omp atomic write
        *status = local_status;
      }
    } //end omp parallel
  }

  if(*status == 0) {
    fftw_destroy_plan(forward_plan);
    fftw_destroy_plan(reverse_plan);
  }

  free(ulocal);
  //TODO: free this up
  fftw_free(a_tmp);
  fftw_free(b_tmp);
}

void ccl_fftlog_ComputeXi2D(double mu, double epsilon,
			    int npk, int N, double *l,double **cl,
			    double *th, double **xi, int *status)
{
  fht(npk, N, l, cl, th, xi, 2., mu, epsilon, 1, 1, NULL, status);
}

void ccl_fftlog_ComputeXi3D(double l, double epsilon,
			    int npk, int N, double *k, double **pk,
			    double *r, double **xi, int *status)
{
  fht(npk, N, k, pk, r, xi, 3., l+0.5, epsilon, 1, 1, NULL, status);
}

max_stars_repo_name: timothydmorton/CCL

max_issues_repo_name: timothydmorton/CCL

max_forks_repo_name: timothydmorton/CCL


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

double epsl = 0.01;

typedef struct {
    double x, y, z;
} position;

typedef struct _vertex {
    position pos;
    gsl_matrix_view *C;
    struct _vertex *p,*l,*r;
} vertex;

void usage(char *exec) {
    printf("%s -s <input_file> <count_vertex>\n", exec);
    printf("%s -c <input_file>\n", exec);
    printf("%s -show <input_file>\n", exec);
}

int* ddd(char *fileName, int len){
    FILE *file = NULL;
    if ((file = fopen(fileName, "r")) == NULL) {
        /* error openning the file */
        perror("fopen: ");
        return NULL;
    }
    int n, i, j = 0, k = 0;
    char buff[6], buffFloat[17];
    memset(buff, 0, sizeof(buff));
    memset(buffFloat, 0, sizeof(buffFloat));
    if (fread(buff, sizeof(char), 5, file) < 5 || strcmp(buff, "E = [")) {
        perror("file out of params: ");
        return NULL;
    }

    double **mat = (double **)calloc(len, sizeof(double *));
    for (i = 0; i < len; i++) {
        mat[i] = (double *)calloc(3, sizeof(double *));
    }
    i = 0;
    while ((n = fread(buff, sizeof(char), 1, file)) == 1) {
        if (buff[0] == ' ') {
            if (k > 0) {
                mat[i][j] = atof(buffFloat);
                memset(buffFloat, 0, sizeof(buffFloat));
                k = 0;
                j++;
            }
            continue;
        }
        else if (buff[0] == ';') {
            if (k) {
                mat[i][j] = atof(buffFloat);
                memset(buffFloat, 0, sizeof(buffFloat));
                k = 0;
            }
            j = 0;
            i++;
        }
        else {
            buffFloat[k++] = buff[0];
        }
    }
    fclose(file);

    int current = 0;
    for (i = 0; i < len; i++) {

        for (j = 0; j < len; j++) {
            printf(" %f ", sqrt(pow(mat[i][0] - mat[j][0], 2) + pow(mat[i][1] - mat[j][1], 2) + pow(mat[i][2] - mat[j][2], 2)));
            fflush(stdout);
        }
        printf("\n");
    }
}

void printMat(double **mat, int len) {
    int i, j;
    for (i = 0; i < len; i++) {
        for (j = 0; j < len; j++) {
            //if (mat[i][j] != 0) {
                printf(" %f ", mat[i][j]);
                fflush(stdout);
            //}
        }
        printf("\n");
    }
}

double **matrixConstruct(char *fileName, int len){
    FILE *file = NULL;
    if ((file = fopen(fileName, "r")) == NULL) {
        /* error openning the file */
        perror("fopen: ");
        return NULL;
    }
    int n, i, j = 0, k = 0;
    char buff[6], buffFloat[17];
    memset(buff, 0, sizeof(buff));
    memset(buffFloat, 0, sizeof(buffFloat));
    if (fread(buff, sizeof(char), 5, file) < 5 || strcmp(buff, "E = [")) {
        perror("file out of params: ");
        return NULL;
    }

    double **mat = (double **)calloc(len, sizeof(double *));
    for (i = 0; i < len; i++) {
        mat[i] = (double *)calloc(len, sizeof(double *));
    }
    i = 0;
    while ((n = fread(buff, sizeof(char), 1, file)) == 1) {
        if (buff[0] == ' ') {
            if (k > 0) {
                mat[i][j] = atof(buffFloat);
                memset(buffFloat, 0, sizeof(buffFloat));
                k = 0;
                j++;
            }
            continue;
        }
        else if (buff[0] == ';') {
            if (k) {
                mat[i][j] = atof(buffFloat);
                memset(buffFloat, 0, sizeof(buffFloat));
                k = 0;
            }
            j = 0;
            i++;
        }
        else {
            buffFloat[k++] = buff[0];
        }
    }
    fclose(file);
    for (i = 0; i < len; i++)
    {
        for (j = 0; j < len; j++)
        {
            if(pow(mat[i][j] - mat[j][i],2) > 0.0000001) return NULL;
        }
    }

    return mat;
}

double dis(double **m, int i){
    if (i-1 < 0) return -1;
    return m[i][i-1];
}

double theta(double **m, int i){
    if (i-1 < 0) return -1000;
    return acos((pow(m[i][i+1],2) + pow(m[i][i-1],2) - pow(m[i+1][i-1],2))/(2*m[i][i+1]*m[i][i-1]));
}

double cosO(double **m, int i){
    //printf("%f %f\n", pow(m[i][i+1], 2) + pow(m[i+1][i+3],2) - 2*m[i][i+1]*m[i+1][i+3]*cos(theta(m, i+1))*cos(theta(m, i+2)) - pow(m[i][i+3],2), 2*m[i][i+1]*m[i+1][i+3]*sin(theta(m, i+1))*sin(theta(m, i+2)) );
    double theta2 = acos((pow(m[i+1][i+2],2) + pow(m[i+1][i+3],2) - pow(m[i+2][i+3],2))/(2*m[i+1][i+3]*m[i+1][i+2]));
    double res = (pow(m[i][i+1], 2) + pow(m[i+1][i+3],2) - 2*m[i][i+1]*m[i+1][i+3]*cos(theta(m, i+1))*cos(theta2) - pow(m[i][i+3],2))/
                    (2*m[i][i+1]*m[i+1][i+3]*sin(theta(m, i+1))*sin(theta2));
    if (res > 1) return 1;
    if (res < -1) return -1;
    return res;
}

double cosOmega(double **m, int i){
    if (i-3 < 0) return 1000;

    double Ai = (2*pow(m[i-2][i-1],2))*(pow(m[i-3][i-2],2) + pow(m[i-2][i] ,2) - pow(m[i-3][i] ,2));
    double Bi = pow(m[i-3][i-2],2) + pow(m[i-2][i-1],2) - pow(m[i-3][i-1],2);
    double Ci = pow(m[i-2][i-1],2) + pow(m[i-2][i],2) - pow(m[i-1][i],2);
    double Di = sqrt(
        4*pow(m[i-3][i-2],2)
        *pow(m[i-2][i-1],2)
        -pow(Bi,2) );
    double Ei = sqrt(4*pow(m[i-2][i-1],2)*pow(m[i-2][i],2)-pow(Ci,2) );
    double res = (Ai - Bi*Ci)/(Di*Ei);
    if (res > 1) return 1;
    if (res < -1) return -1;
    return res;
}

int isFeasible(double **mat, int len, vertex *v, int i){
    if (!v) return 0;
    int j = i-1;
    double pos[i][3];
    vertex *aux = v->p;
    while (aux)
    {
        //printf("\n\tj = %d {%f, %f, %f}\n", j, aux->pos.x, aux->pos.y, aux->pos.z);
        pos[j][0] = aux->pos.x;
        pos[j][1] = aux->pos.y;
        pos[j--][2] = aux->pos.z;
        aux = aux->p;
    }
    for (j = 0; j < i-3; j++)
    {
        if (mat[i][j] > 0)
        {
            double dij = pow(v->pos.x - pos[j][0],2)+pow(v->pos.y - pos[j][1],2)+pow(v->pos.z - pos[j][2],2);
            double diff = pow(mat[i][j],2)- dij;
            if (diff < 0) diff = diff*(-1);
            double dist = sqrt(diff);//pow(diff,2);
            //printf("{%d,%d} {%f, %f, %f}-{%f, %f, %f} \n\t%.20f- %.20f = %.20f \n %.20f \n", i, j, v->pos.x, v->pos.y, v->pos.z, pos[j][0], pos[j][1], pos[j][2], mat[i][j], dij, diff, dist);
            if (dist > epsl)
            {
               return 0;
            }
        }
    }
    return 1;
}

double LDE(double **mat, double **m, int len){
    int i, j, count = 0;
    double dij = 0;
    for (i = 0; i < len; i++)
    {
        for (j = i+1; j < len; j++)
        {
            if (mat[i][j]>0){
                dij = dij + pow((pow(mat[i][j],2) - (pow(m[i][0] - m[j][0], 2) + pow(m[i][1] - m[j][1], 2) + pow(m[i][2] - m[j][2], 2))),2)/mat[i][j];
                count++;
            }
        }
    }
    return dij/count;

}

int qtd = 0;
int BranchAndPrune(double **mat, int len, vertex *v, int i){
    if(i < len){
        double thetai = theta(mat,i-1);
        double cti = cos(thetai);
        double sti = sin(thetai);
        double cwi = cosOmega(mat, i);
        //double cwi2 = cosO(mat, i-3);
        //printf("%f %f \n", cwi, cwi2 );
        //printf("%f %f \n", acos(cwi)*57.295779513, acos(cwi2)*57.295779513 );
        if (cwi*cwi > 1){
            printf("\n %f Error cwi", cwi);
            getchar();
        }
        double swi = sqrt(1-pow(cwi, 2));
        double di = dis(mat,i);
        double bi1[] = {0-cti,0-sti,0,0-di*cti,
                        sti*cwi,0-cti*cwi,0-swi,di*sti*cwi,
                        sti*swi,0-cti*swi,cwi,di*sti*swi,
                        0,0,0,1};
        int j;
        //for (j = 0; j < 16; j++)
        //{
       //     printf("%f ", bi1[j]);
        //}
       // printf("\n");

        gsl_matrix_view Bi1 = gsl_matrix_view_array(bi1, 4, 4);

        swi = 0-swi;
        double bi2[] = {0-cti,0-sti,0,0-di*cti,
                        sti*cwi,0-cti*cwi,0-swi,di*sti*cwi,
                        sti*swi,0-cti*swi,cwi,di*sti*swi,
                        0,0,0,1};

        //for (j = 0; j < 16; j++)
        //{
        //    printf("%f ", bi2[j]);
        //}
        //printf("\n");

        gsl_matrix_view Bi2 = gsl_matrix_view_array(bi2, 4, 4);

        vertex *rv = (vertex*)calloc(1, sizeof(vertex));
        rv->p = v;
        rv->C = (gsl_matrix_view*)calloc(1, sizeof(gsl_matrix_view));
        *(rv->C) = gsl_matrix_view_array((double*)calloc(16,sizeof(double)), 4, 4);

        gsl_blas_dgemm(CblasNoTrans, CblasNoTrans,
                   1.0, &(v->C->matrix), &Bi1.matrix,
                   0.0, &(rv->C->matrix));

        vertex *lv = (vertex*)calloc(1, sizeof(vertex));
        lv->p = v;
        lv->C = (gsl_matrix_view*)calloc(1, sizeof(gsl_matrix_view));
        *(lv->C) = gsl_matrix_view_array((double*)calloc(16,sizeof(double)), 4, 4);

        gsl_blas_dgemm(CblasNoTrans, CblasNoTrans,
                   1.0, &(v->C->matrix), &Bi2.matrix,
                   0.0, &(lv->C->matrix));

        gsl_matrix_view *Xi = (gsl_matrix_view*)calloc(1, sizeof(gsl_matrix_view));
        *Xi = gsl_matrix_view_array((double*)calloc(4,sizeof(double)), 4, 1);

        double *y = (double*)calloc(4,sizeof(double));
        y[3] = 1;
        gsl_matrix_view *Y = (gsl_matrix_view*)calloc(1, sizeof(gsl_matrix_view));
        *Y = gsl_matrix_view_array(y, 4, 1);

        gsl_blas_dgemm(CblasNoTrans, CblasNoTrans,
                    1.0, &(rv->C->matrix), &(Y->matrix),
                    0.0, &(Xi->matrix));
        rv->pos.x = Xi->matrix.data[0];
        rv->pos.y = Xi->matrix.data[1];
        rv->pos.z = Xi->matrix.data[2];

        gsl_blas_dgemm(CblasNoTrans, CblasNoTrans,
                    1.0, &(lv->C->matrix), &(Y->matrix),
                    0.0, &(Xi->matrix));
        lv->pos.x = Xi->matrix.data[0];
        lv->pos.y = Xi->matrix.data[1];
        lv->pos.z = Xi->matrix.data[2];
        free(y);
        free(Y);
        free(Xi->matrix.data);
        free(Xi);

        //printf("[ %f, %f", rv->pos.x, rv->pos.y);
        //printf("  %f]\n", rv->pos.z);
        //printf("[ %f, %f", lv->pos.x, lv->pos.y);
        //printf("  %f ]\n", lv->pos.z);

        //getchar();
            int s;
            //printf("\n");
        if (isFeasible(mat, len, rv, i)){
            v->r = rv;
            //for (s = 0; s < i; s++)
            //{
            //    printf(">");
            //}

            //printf("1\n");
            fflush(stdout);
            BranchAndPrune(mat, len, rv, i+1);
        }else {
            //for (s = 0; s < i; s++)
            //{
            //    printf("<");
            //}
            //
            //printf("2\n");
            free(rv->C->matrix.data);
            free(rv->C);
            free(rv);
            v->r = NULL;
        }

        if (isFeasible(mat, len, lv, i)){
            v->l = lv;

            //for (s = 0; s < i; s++)
            //{
            //    printf(">");
            //}
            //
            //printf("3\n");
            fflush(stdout);
            BranchAndPrune(mat, len, lv, i+1);
        }else {

            //for (s = 0; s < i; s++)
            //{
            //    printf("<");
            //}
            //
            //printf("4\n");
            free(lv->C->matrix.data);
            free(lv->C);
            free(lv);
            v->l = NULL;
        }
    }else{
        //printf(" ---------------- ue\n");
        //fflush(stdout);
        qtd++;
        double **mat2 = calloc(len, sizeof(double*));
        int i;
        for (i = 0; i < len; i++)
        {
            mat2[i] = calloc(3, sizeof(double));
        }
        i = len-1;
        printf("\n[");
        while (v)
        {
            mat2[i][0] = v->pos.x;
            mat2[i][1] = v->pos.y;
            mat2[i--][2] = v->pos.z;
            //printf("%f %f %f; ", v->pos.x, v->pos.y, v->pos.z);
            v = v->p;
        }

        double lde = LDE(mat, mat2, len);

        printf("\n \t LDE = %.20f\n", lde);

        for (i = 0; i < len; i++)
        {
            printf("%f %f %f; ", mat[i][0], mat[i][1], mat[i][2]);
        }

        printf("]\n");
        return 0;
    }
}

int printGrafo(int len, vertex *v, int i){
    int j = -1;
    if (v->l){
        printf("(%d, %d),", i, i+1);
        j = printGrafo(len, v->l, i+1);
    }
    if (v->r){
        if (j != -1){
            printf("(%d, %d),", i, j);
            j = printGrafo(len, v->r, j);
        }else
        {
            printf("(%d, %d),", i, i+1);
            return printGrafo(len, v->r, i+1);
        }
    }
    if (j == -1) return i+1;
    return j;
}

vertex *resul[10];

int solve(char *fileName, int len) {
    double **mat = matrixConstruct(fileName, len);
    if(mat == NULL){
        printf("Error when construct matrix\n");
        fflush(stdout);
        return 0;
    }

    double b1[] = {1,0,0,0,
                    0,1,0,0,
                    0,0,1,0,
                    0,0,0,1};
    gsl_matrix_view B1 = gsl_matrix_view_array(b1, 4, 4);

    double b2[] = {-1,0,0,0-dis(mat,1),
                    0,1,0,0,
                    0,0,-1,0,
                    0,0,0,1};
    gsl_matrix_view B2 = gsl_matrix_view_array(b2, 4, 4);

    double theta2 = theta(mat,1);
    //printf("[ cos %f\n", theta2);

    double ct2 = cos(theta2);
    double st2 = sin(theta2);
    double b3[] = {0-ct2,0-st2,0,0-dis(mat,2)*ct2,
                    st2,0-ct2,0,dis(mat,2)*st2,
                    0,0,1,0,
                    0,0,0,1};
    gsl_matrix_view B3 = gsl_matrix_view_array(b3, 4, 4);

    vertex *v = (vertex*)calloc(1, sizeof(vertex));
    v->C = &B1;
    //printf("[ %f, %f %f\n", v->pos.x, v->pos.y, v->pos.z);

    vertex *v2 = (vertex*)calloc(1, sizeof(vertex));
    v->r = v2;
    v2->C = &B2;
    v2->p = v;
    v2->pos.x = 0-dis(mat,1);

     //   printf("[ %f, %f %f\n", v2->pos.x, v2->pos.y, v2->pos.z);
    vertex *v3 = (vertex*)calloc(1, sizeof(vertex));
    v2->r = v3;
    v3->C = (gsl_matrix_view*)calloc(1, sizeof(gsl_matrix_view));
    *(v3->C) = gsl_matrix_view_array((double*)calloc(16,sizeof(double)), 4, 4);
    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans,
                   1.0, &B2.matrix, &(B3.matrix),
                   0.0, &(v3->C->matrix));
    v3->p = v2;
    v3->pos.x = dis(mat,2)*ct2 -dis(mat,1);
    v3->pos.y = dis(mat,2)*st2;

    //printf("[ %f, %f", v3->pos.x, v3->pos.y);
    //    printf("  %f]\n", v3->pos.z);

    BranchAndPrune(mat, len, v3, 3);
    printf("%d", qtd);
    //printf("[");
    //printGrafo(len, v, 0);
    //printf("]");
    return 1;
}

int countCols(char *fileName) {
    FILE *file = NULL;
    if ((file = fopen(fileName, "r")) == NULL) {
        /* error openning the file */
        perror("fopen: ");
        return 1;
    }

    int n, j = 0, k = 0;
    char buff[6];
    memset(buff, 0, sizeof(buff));
    if (fread(buff, sizeof(char), 5, file) < 5 || strcmp(buff, "E = [")) {
        fclose(file);
        perror("file out of params: ");
        return 1;
    }

    while ((n = fread(buff, sizeof(char), 1, file)) == 1) {
        if (buff[0] == ' ') {
            if (k > 0) {
                k = 0;
                j++;
            }
            continue;
        }
        else if (buff[0] == ';') {
            if (k) j++;
            break;
        }
        else {
            k++;
        }
    }
    fclose(file);
    return j;
}

int main(int argc, char **argv) {
    double time_spent = 0.0;

	clock_t begin = clock();

    if (argc < 3) {
        usage(argv[0]);
    }
    else {
        if (!strcmp(argv[1], "-s")) {
            int len;
            if (argc < 4) {
                len = countCols(argv[2]);
            }else {
                if ((len = atoi(argv[3])) == 0){
                    /* error openning the file */
                    perror("count_vertex param error: ");
                    return 1;
                }
            }

            solve(argv[2], len);
        } else if (!strcmp(argv[1], "-c")) {
            printf("%d\n",countCols(argv[2]));
        } else if (!strcmp(argv[1], "-show")) {
            printMat(matrixConstruct(argv[2], countCols(argv[2])), countCols(argv[2]));
        } else if (!strcmp(argv[1], "-d")) {
            ddd(argv[2],10);
        } else {
            usage(argv[0]);
        }
    }

    //getchar();

    clock_t end = clock();

	// calculate elapsed time by finding difference (end - begin) and
	// divide by CLOCKS_PER_SEC to convert to seconds
	time_spent += (double)(end - begin) / CLOCKS_PER_SEC;

	printf("\nTime elpased is %.20f seconds", time_spent);
    fflush(stdout);
}

max_stars_repo_name: caomem/BP

max_issues_repo_name: caomem/BP

max_forks_repo_name: caomem/BP


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

int main(void) {
    /* coefficients of f(x) = -1 + x^5 */
    double a[6] = {-1, 0, 0, 0, 0, 1};
    printf("x = %f\n", gsl_poly_eval(a, 6, 3.0));
    return 0;
}

// int main(void) {
//     int i;
//     /* coefficients of P(x) = -1 + x^5 */
//     double a[6] = {-1, 0, 0, 0, 0, 1};
//     double z[10];
//     gsl_poly_complex_workspace* w = gsl_poly_complex_workspace_alloc(6);
//     gsl_poly_complex_solve(a, 6, w, z);
//     gsl_poly_complex_workspace_free(w);
//     for (i = 0; i < 5; i++) {
//         printf("z%d = %+.18f %+.18f\n", i, z[2 * i], z[2 * i + 1]);
//     }
//     return 0;
// }

max_stars_repo_name: shakfu/pd-psl

max_issues_repo_name: shakfu/pd-psl

max_forks_repo_name: shakfu/pd-psl


text:
#ifndef TOPICH
#define TOPICH

#include <gsl/gsl_vector.h>
#include <math.h>
#include "utils.h"
#include "typedefs.h"
#include "hyperparameter.h"

#define MH_ETA_STDEV 0.005
#define MH_GAM_STDEV 0.005

/* === topic methods === */

/*
 * update the count of a word in a topic
 *
 */

void topic_update_word(topic* t, int w, double update);

/*
 * update the document count in a topic
 *
 */

void topic_update_doc_cnt(topic* t, double update);

/*
 * write log probability in a file
 *
 */

void topic_write_log_prob(topic* t, FILE* f);

/* === tree methods === */

/*
 * allocate a new tree with a certain depth
 *
 */

tree* tree_new(int depth,
               int nwords,
               gsl_vector* eta,
               gsl_vector* gam,
               double scaling_shape,
               double scaling_scale);

/*
 * add children to a node down to the depth of the tree;
 * return the leaf of that path
 *
 */

topic* tree_fill(topic* t);

/*
 * add a child to a topic
 *
 */

topic* topic_add_child(topic* t);

/*
 * make a new topic
 *
 */

topic* topic_new(int nwords, int level, topic* parent, tree* tr);

/*
 * given a leaf with 0 instances, delete the leaf and all ancestors
 * which also have 0 instances.  (!!! use asserts here)
 *
 */

void tree_prune(topic* t);

/*
 * delete a node from the tree
 *
 */

void delete_node(topic* t);

/*
 * write a tree to file
 *
 */

void tree_write_log_prob(tree* tree, FILE* f);

/*
 * sample a document path from a tree
 *
 */

void populate_prob_dfs(topic* node, doc* d, double* logsum, double* pprob, int root_level);
void tree_sample_doc_path(tree* tr, doc* d, short do_remove, int root_level);

/*
 * sample a new path in the tree for a document
 *
 */

void tree_sample_path_for_doc(tree* t, doc* d);

/*
 * update the tree from an entire document
 *
 */

void tree_update_from_doc(doc* d, double update, int root_level);

/*
 * sample a leaf from the tree with populated probabilties
 *
 */

topic* tree_sample_path(topic* node, double logsum);
topic* tree_sample_dfs(double r, topic* node, double* sum, double logsum);
void tree_add_doc_to_path(topic* node, doc* d, int root_level);
void tree_remove_doc_from_path(tree* tr, doc* d, int root_level);

/*
 * write a tree to a file
 *
 */

void write_tree(tree* tf, FILE* file);
void write_tree_levels(tree* tr, FILE* file);
void write_tree_level_dfs(topic* t, FILE* f);
void write_tree_topics_dfs(topic* t, FILE* f);

/*
 * scores
 *
 */

double gamma_score(topic* t);
double gamma_score_PY(topic* t, double gam_add);
double eta_score(topic* t);
double log_gamma_ratio(doc* d, topic* t, int level);
double log_gamma_ratio_new(doc* d, int level, double eta, int nterms);
void tree_mh_update_eta(tree* tr);
void dfs_sample_scaling(topic* t);

/*
 * copying a tree
 *
 */

void copy_topic(const topic* src, topic* dest);
void copy_tree_dfs(const topic* src, topic* dest);
tree * copy_tree(const tree* tr);

void free_tree(tree * tr);
void free_tree_dfs(topic * t);

int ntopics_in_tree(tree * tr);
int ntopics_in_tree_dfs(topic * t);

#endif

max_stars_repo_name: thatindiandude/PGM-Project

max_issues_repo_name: thatindiandude/PGM-Project

max_forks_repo_name: thatindiandude/PGM-Project


text:
/*
   Copyright [2020] [IBM Corporation]
   Licensed under the Apache License, Version 2.0 (the "License");
   you may not use this file except in compliance with the License.
   You may obtain a copy of the License at
       http://www.apache.org/licenses/LICENSE-2.0
   Unless required by applicable law or agreed to in writing, software
   distributed under the License is distributed on an "AS IS" BASIS,
   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   See the License for the specific language governing permissions and
   limitations under the License.
*/

/*
 * Authors:
 *
 */

#ifndef _NUPM_OPENED_SPACE_H_
#define _NUPM_OPENED_SPACE_H_

#include <common/byte_span.h>
#include <common/fd_locked.h>
#include <common/logging.h> /* log_source */
#include <common/memory_mapped.h>
#include <common/moveable_ptr.h>
#include <common/string_view.h>
#include <common/types.h> /* addr_t */
#include <gsl/pointers>
#include <cassert>
#include <cstddef>
#include <string>
#include <vector>

namespace nupm
{
struct dax_manager;

struct range_use
{
private:
  common::moveable_ptr<dax_manager> _dm;
  /* Note: arena_fs used multiple ranges */
  std::vector<common::memory_mapped> _iovm;

  std::vector<common::memory_mapped> address_coverage_check(std::vector<common::memory_mapped> &&iovm);
  using byte = common::byte;
  using byte_span = common::byte_span;
public:
#if 0
  const auto & operator[](std::size_t i) const { return _iovm.at(i).iov(); }
#else
  byte_span operator[](std::size_t i) const { const auto &iov = _iovm.at(i).iov(); return common::make_byte_span(::base(iov), ::size(iov)); }
#endif
  range_use(dax_manager *dm_, std::vector<common::memory_mapped> &&);
  range_use(const range_use &) = delete;
  range_use &operator=(const range_use &) = delete;
  range_use(range_use &&) noexcept = default;
  void grow(std::vector<common::memory_mapped> &&);
  void shrink(std::size_t size);
  ~range_use();
  gsl::not_null<dax_manager *> dm() const { return _dm; }
  ::off_t size() const;
};

struct space_opened : private common::log_source
{
private:
  using byte_span = common::byte_span;
  common::fd_locked _fd_locked;
  /* Note: arena_fs may someday use multiple ranges */
  range_use _range;

  /* owns the file mapping */
  std::vector<common::memory_mapped> map_dev(int fd, const addr_t base_addr);
  std::vector<common::memory_mapped> map_fs(int fd, const std::vector<byte_span> &mapping, ::off_t offset);
public:
  space_opened(const common::log_source &, dax_manager * dm_, common::fd_locked && fd, const addr_t base_addr);
  space_opened(const common::log_source &, dax_manager * dm_, common::fd_locked && fd, const std::vector<byte_span> &mapping);
  space_opened(space_opened &&) noexcept = default;
  void shrink(std::size_t size);
  void grow(std::vector<byte_span> && iovv);
  int fd() const { return _fd_locked.fd(); }
  range_use &range() { return _range; }
};
}
#endif


text:
#ifndef SDCA_MATH_BLAS_H
#define SDCA_MATH_BLAS_H

#if defined(BLAS_MKL)
  //#include <mkl.h> // produces warnings with -pedantic
  #include <mkl_cblas.h>
#elif defined(BLAS_ACCELERATE)
  #include <Accelerate/Accelerate.h>
#elif defined(BLAS_MATLAB)
  #include <cstddef>
  typedef std::ptrdiff_t ptrdiff_t; // Matlab's blas.h needs ptrdiff_t
  #include <blas.h>
#elif defined(BLAS_DEFAULT)
  #include <cblas.h>
#elif defined(BLAS_DEFAULT_LOCAL_HEADER)
  #include "cblas.h"
#endif

#if defined(BLAS_MATLAB)
  enum CBLAS_TRANSPOSE {CblasNoTrans=111, CblasTrans=112, CblasConjTrans=113};
#endif

namespace sdca {

#if defined(MKL_INT)
  typedef MKL_INT blas_int;
#elif defined(BLAS_MATLAB)
  typedef ptrdiff_t blas_int;
#else
  typedef int blas_int;
#endif

inline void
sdca_blas_scal(
    const blas_int n,
    const float alpha,
    float* X
  ) {
#if defined(sscal)
  blas_int inc(1);
  sscal(const_cast<blas_int*>(&n), const_cast<float*>(&alpha), X, &inc);
#else
  cblas_sscal(n, alpha, X, 1);
#endif
}

inline void
sdca_blas_scal(
    const blas_int n,
    const double alpha,
    double* X
  ) {
#if defined(dscal)
  blas_int inc(1);
  dscal(const_cast<blas_int*>(&n), const_cast<double*>(&alpha), X, &inc);
#else
  cblas_dscal(n, alpha, X, 1);
#endif
}

inline void
sdca_blas_copy(
    const blas_int n,
    const float* X,
    float* Y
  ) {
#if defined(scopy)
  blas_int inc(1);
  scopy(const_cast<blas_int*>(&n), const_cast<float*>(X), &inc, Y, &inc);
#else
  cblas_scopy(n, X, 1, Y, 1);
#endif
}

inline void
sdca_blas_copy(
    const blas_int n,
    const double* X,
    double* Y
  ) {
#if defined(dcopy)
  blas_int inc(1);
  dcopy(const_cast<blas_int*>(&n), const_cast<double*>(X), &inc, Y, &inc);
#else
  cblas_dcopy(n, X, 1, Y, 1);
#endif
}

inline void
sdca_blas_swap(
    const blas_int n,
    float* X,
    float* Y
  ) {
#if defined(scopy)
  blas_int inc(1);
  sswap(const_cast<blas_int*>(&n), X, &inc, Y, &inc);
#else
  cblas_sswap(n, X, 1, Y, 1);
#endif
}

inline void
sdca_blas_swap(
    const blas_int n,
    double* X,
    double* Y
  ) {
#if defined(dcopy)
  blas_int inc(1);
  dswap(const_cast<blas_int*>(&n), X, &inc, Y, &inc);
#else
  cblas_dswap(n, X, 1, Y, 1);
#endif
}

inline void
sdca_blas_axpy(
    const blas_int n,
    const float alpha,
    const float* X,
    float* Y
  ) {
#if defined(saxpy)
  blas_int inc(1);
  saxpy(const_cast<blas_int*>(&n), const_cast<float*>(&alpha),
        const_cast<float*>(X), &inc, Y, &inc);
#else
  cblas_saxpy(n, alpha, X, 1, Y, 1);
#endif
}

inline void
sdca_blas_axpy(
    const blas_int n,
    const double alpha,
    const double* X,
    double* Y
  ) {
#if defined(daxpy)
  blas_int inc(1);
  daxpy(const_cast<blas_int*>(&n), const_cast<double*>(&alpha),
        const_cast<double*>(X), &inc, Y, &inc);
#else
  cblas_daxpy(n, alpha, X, 1, Y, 1);
#endif
}

inline void
sdca_blas_axpby(
    const blas_int n,
    const float alpha,
    const float* X,
    const float beta,
    float* Y
  ) {
#if defined(cblas_saxpby)
  cblas_saxpby(n, alpha, X, 1, beta, Y, 1);
#elif defined(sscal)
  blas_int inc(1);
  sscal(const_cast<blas_int*>(&n), const_cast<float*>(&beta), Y, &inc);
  saxpy(const_cast<blas_int*>(&n), const_cast<float*>(&alpha),
        const_cast<float*>(X), &inc, Y, &inc);
#else
  cblas_sscal(n, beta, Y, 1);
  cblas_saxpy(n, alpha, X, 1, Y, 1);
#endif
}

inline void
sdca_blas_axpby(
    const blas_int n,
    const double alpha,
    const double* X,
    const double beta,
    double* Y
  ) {
#if defined(cblas_daxpby)
  cblas_daxpby(n, alpha, X, 1, beta, Y, 1);
#elif defined(dscal)
  blas_int inc(1);
  dscal(const_cast<blas_int*>(&n), const_cast<double*>(&beta), Y, &inc);
  daxpy(const_cast<blas_int*>(&n), const_cast<double*>(&alpha),
        const_cast<double*>(X), &inc, Y, &inc);
#else
  cblas_dscal(n, beta, Y, 1);
  cblas_daxpy(n, alpha, X, 1, Y, 1);
#endif
}

inline float
sdca_blas_dot(
    const blas_int n,
    const float* X,
    const float* Y
  ) {
#if defined(sdot)
  blas_int inc(1);
  return sdot(const_cast<blas_int*>(&n), const_cast<float*>(X), &inc,
              const_cast<float*>(Y), &inc);
#else
  return cblas_sdot(n, X, 1, Y, 1);
#endif
}

inline double
sdca_blas_dot(
    const blas_int n,
    const double* X,
    const double* Y
  ) {
#if defined(ddot)
  blas_int inc(1);
  return ddot(const_cast<blas_int*>(&n), const_cast<double*>(X), &inc,
              const_cast<double*>(Y), &inc);
#else
  return cblas_ddot(n, X, 1, Y, 1);
#endif
}

inline float
sdca_blas_asum(
    const blas_int n,
    const float* X
  ) {
#if defined(sasum)
  blas_int inc(1);
  return sasum(const_cast<blas_int*>(&n), const_cast<float*>(X), &inc);
#else
  return cblas_sasum(n, X, 1);
#endif
}

inline double
sdca_blas_asum(
    const blas_int n,
    const double* X
  ) {
#if defined(dasum)
  blas_int inc(1);
  return dasum(const_cast<blas_int*>(&n), const_cast<double*>(X), &inc);
#else
  return cblas_dasum(n, X, 1);
#endif
}

inline float
sdca_blas_nrm2(
    const blas_int n,
    const float* X
  ) {
#if defined(snrm2)
  blas_int inc(1);
  return snrm2(const_cast<blas_int*>(&n), const_cast<float*>(X), &inc);
#else
  return cblas_snrm2(n, X, 1);
#endif
}

inline double
sdca_blas_nrm2(
    const blas_int n,
    const double* X
  ) {
#if defined(dnrm2)
  blas_int inc(1);
  return dnrm2(const_cast<blas_int*>(&n), const_cast<double*>(X), &inc);
#else
  return cblas_dnrm2(n, X, 1);
#endif
}

inline void
sdca_blas_gemv(
    const blas_int m,
    const blas_int n,
    const float* A,
    const float* X,
    float* Y,
    const CBLAS_TRANSPOSE transA = CblasNoTrans,
    const float alpha = 1,
    const float beta = 0
    ) {
#if defined(sgemv)
  blas_int inc(1);
  char trans((transA == CblasNoTrans) ? 'N' : 'T');
  sgemv(&trans, const_cast<blas_int*>(&m), const_cast<blas_int*>(&n),
        const_cast<float*>(&alpha),
        const_cast<float*>(A), const_cast<blas_int*>(&m),
        const_cast<float*>(X), &inc, const_cast<float*>(&beta), Y, &inc);
#else
  cblas_sgemv(CblasColMajor, transA, m, n, alpha, A, m,
              X, 1, beta, Y, 1);
#endif
}

inline void
sdca_blas_gemv(
    const blas_int m,
    const blas_int n,
    const double* A,
    const double* X,
    double* Y,
    const CBLAS_TRANSPOSE transA = CblasNoTrans,
    const double alpha = 1,
    const double beta = 0
    ) {
#if defined(dgemv)
  blas_int inc(1);
  char trans((transA == CblasNoTrans) ? 'N' : 'T');
  dgemv(&trans, const_cast<blas_int*>(&m), const_cast<blas_int*>(&n),
        const_cast<double*>(&alpha),
        const_cast<double*>(A), const_cast<blas_int*>(&m),
        const_cast<double*>(X), &inc, const_cast<double*>(&beta), Y, &inc);
#else
  cblas_dgemv(CblasColMajor, transA, m, n, alpha, A, m,
              X, 1, beta, Y, 1);
#endif
}

inline void
sdca_blas_ger(
    const blas_int m,
    const blas_int n,
    const float alpha,
    const float* X,
    const float* Y,
    float* A
    ) {
#if defined(sger)
  blas_int inc(1);
  sger(const_cast<blas_int*>(&m), const_cast<blas_int*>(&n),
       const_cast<float*>(&alpha), const_cast<float*>(X), &inc,
       const_cast<float*>(Y), &inc, A, const_cast<blas_int*>(&m));
#else
  cblas_sger(CblasColMajor, m, n, alpha, X, 1, Y, 1, A, m);
#endif
}

inline void
sdca_blas_ger(
    const blas_int m,
    const blas_int n,
    const double alpha,
    const double* X,
    const double* Y,
    double* A
    ) {
#if defined(dger)
  blas_int inc(1);
  dger(const_cast<blas_int*>(&m), const_cast<blas_int*>(&n),
       const_cast<double*>(&alpha), const_cast<double*>(X), &inc,
       const_cast<double*>(Y), &inc, A, const_cast<blas_int*>(&m));
#else
  cblas_dger(CblasColMajor, m, n, alpha, X, 1, Y, 1, A, m);
#endif
}

inline void
sdca_blas_gemm(
    const blas_int m,
    const blas_int n,
    const blas_int k,
    const float* A,
    const blas_int lda,
    const float* B,
    const blas_int ldb,
    float* C,
    const CBLAS_TRANSPOSE transA = CblasNoTrans,
    const CBLAS_TRANSPOSE transB = CblasNoTrans,
    const float alpha = 1,
    const float beta = 0
    ) {
#if defined(sgemm)
  char transa((transA == CblasNoTrans) ? 'N' : 'T');
  char transb((transB == CblasNoTrans) ? 'N' : 'T');
  sgemm(&transa, &transb, const_cast<blas_int*>(&m), const_cast<blas_int*>(&n),
        const_cast<blas_int*>(&k), const_cast<float*>(&alpha),
        const_cast<float*>(A), const_cast<blas_int*>(&lda),
        const_cast<float*>(B), const_cast<blas_int*>(&ldb),
        const_cast<float*>(&beta), C, const_cast<blas_int*>(&m));
#else
  cblas_sgemm(CblasColMajor, transA, transB, m, n, k,
              alpha, A, lda, B, ldb, beta, C, m);
#endif
}

inline void
sdca_blas_gemm(
    const blas_int m,
    const blas_int n,
    const blas_int k,
    const double* A,
    const blas_int lda,
    const double* B,
    const blas_int ldb,
    double* C,
    const CBLAS_TRANSPOSE transA = CblasNoTrans,
    const CBLAS_TRANSPOSE transB = CblasNoTrans,
    const double alpha = 1,
    const double beta = 0
    ) {
#if defined(dgemm)
  char transa((transA == CblasNoTrans) ? 'N' : 'T');
  char transb((transB == CblasNoTrans) ? 'N' : 'T');
  dgemm(&transa, &transb, const_cast<blas_int*>(&m), const_cast<blas_int*>(&n),
        const_cast<blas_int*>(&k), const_cast<double*>(&alpha),
        const_cast<double*>(A), const_cast<blas_int*>(&lda),
        const_cast<double*>(B), const_cast<blas_int*>(&ldb),
        const_cast<double*>(&beta), C, const_cast<blas_int*>(&m));
#else
  cblas_dgemm(CblasColMajor, transA, transB, m, n, k,
              alpha, A, lda, B, ldb, beta, C, m);
#endif
}

}

#endif

max_stars_repo_name: mlapin/libsdca

max_issues_repo_name: mlapin/libsdca

max_forks_repo_name: mlapin/libsdca


text:
#include <gsl/gsl_math.h>
#include <gsl/gsl_cblas.h>
#include "cblas.h"

void
cblas_sgbmv (const enum CBLAS_ORDER order, const enum CBLAS_TRANSPOSE TransA,
             const int M, const int N, const int KL, const int KU,
             const float alpha, const float *A, const int lda, const float *X,
             const int incX, const float beta, float *Y, const int incY)
{
#define BASE float
#include "source_gbmv_r.h"
#undef BASE
}

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 <stdarg.h>
#include <string.h>
#include <math.h>
#include <gbpLib.h>
#include <gbpRNG.h>
#include <gbpMCMC.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_fit.h>
#include <gsl/gsl_interp.h>

void free_MCMC(MCMC_info *MCMC) {
    int           i_P, i_DS;
    int           i_array;
    MCMC_DS_info *current_DS;
    MCMC_DS_info *next_DS;

    SID_log("Freeing MCMC structure...", SID_LOG_OPEN);

    // Parameter arrays
    for(i_P = 0; i_P < MCMC->n_P; i_P++)
        SID_free(SID_FARG MCMC->P_names[i_P]);
    SID_free(SID_FARG MCMC->P_names);
    SID_free(SID_FARG MCMC->P_init);
    SID_free(SID_FARG MCMC->P_new);
    SID_free(SID_FARG MCMC->P_last);
    SID_free(SID_FARG MCMC->P_chain);
    SID_free(SID_FARG MCMC->P_limit_min);
    SID_free(SID_FARG MCMC->P_limit_max);
    if(MCMC->n_arrays > 0) {
        for(i_array = 0; i_array < MCMC->n_arrays; i_array++) {
            SID_free(SID_FARG MCMC->array[i_array]);
            SID_free(SID_FARG MCMC->array_name[i_array]);
        }
        SID_free(SID_FARG MCMC->array);
        SID_free(SID_FARG MCMC->array_name);
    }

    // Covariance and displacement vector
    free_MCMC_covariance(MCMC);

    // Random number generator
    if(MCMC->RNG != NULL)
        free_RNG(MCMC->RNG);

    // Dataset arrays
    free_MCMC_arrays(MCMC);
    free_MCMC_DS(MCMC);

    // Communicators
    SID_Comm_free(&(MCMC->comm));

    SID_log("Done.", SID_LOG_CLOSE);
}

max_stars_repo_name: gbpoole/gbpCode

max_issues_repo_name: gbpoole/gbpCode

max_forks_repo_name: gbpoole/gbpCode


text:
/* eigen/gensymm.c
 *
 * Copyright (C) 2007 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 <stdlib.h>

#include <config.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>

/*
 * This module computes the eigenvalues of a real generalized
 * symmetric-definite eigensystem A x = \lambda B x, where A and
 * B are symmetric, and B is positive-definite.
 */

/*
gsl_eigen_gensymm_alloc()

Allocate a workspace for solving the generalized symmetric-definite
eigenvalue problem. The size of this workspace is O(2n).

Inputs: n - size of matrices

Return: pointer to workspace
*/

gsl_eigen_gensymm_workspace *
gsl_eigen_gensymm_alloc(const size_t n)
{
  gsl_eigen_gensymm_workspace *w;

  if (n == 0)
    {
      GSL_ERROR_NULL ("matrix dimension must be positive integer",
                      GSL_EINVAL);
    }

  w = (gsl_eigen_gensymm_workspace *) calloc (1, sizeof (gsl_eigen_gensymm_workspace));

  if (w == 0)
    {
      GSL_ERROR_NULL ("failed to allocate space for workspace", GSL_ENOMEM);
    }

  w->size = n;

  w->symm_workspace_p = gsl_eigen_symm_alloc(n);
  if (!w->symm_workspace_p)
    {
      gsl_eigen_gensymm_free(w);
      GSL_ERROR_NULL("failed to allocate space for symm workspace", GSL_ENOMEM);
    }

  return (w);
} /* gsl_eigen_gensymm_alloc() */

/*
gsl_eigen_gensymm_free()
  Free workspace w
*/

void
gsl_eigen_gensymm_free (gsl_eigen_gensymm_workspace * w)
{
  RETURN_IF_NULL (w);

  if (w->symm_workspace_p)
    gsl_eigen_symm_free(w->symm_workspace_p);

  free(w);
} /* gsl_eigen_gensymm_free() */

/*
gsl_eigen_gensymm()

Solve the generalized symmetric-definite eigenvalue problem

A x = \lambda B x

for the eigenvalues \lambda.

Inputs: A    - real symmetric matrix
        B    - real symmetric and positive definite matrix
        eval - where to store eigenvalues
        w    - workspace

Return: success or error
*/

int
gsl_eigen_gensymm (gsl_matrix * A, gsl_matrix * B, gsl_vector * eval,
                   gsl_eigen_gensymm_workspace * w)
{
  const size_t N = A->size1;

  /* check matrix and vector sizes */

  if (N != A->size2)
    {
      GSL_ERROR ("matrix must be square to compute eigenvalues", GSL_ENOTSQR);
    }
  else if ((N != B->size1) || (N != B->size2))
    {
      GSL_ERROR ("B matrix dimensions must match A", GSL_EBADLEN);
    }
  else if (eval->size != N)
    {
      GSL_ERROR ("eigenvalue vector must match matrix size", GSL_EBADLEN);
    }
  else if (w->size != N)
    {
      GSL_ERROR ("matrix size does not match workspace", GSL_EBADLEN);
    }
  else
    {
      int s;

      /* compute Cholesky factorization of B */
      s = gsl_linalg_cholesky_decomp(B);
      if (s != GSL_SUCCESS)
        return s; /* B is not positive definite */

      /* transform to standard symmetric eigenvalue problem */
      gsl_eigen_gensymm_standardize(A, B);

      s = gsl_eigen_symm(A, eval, w->symm_workspace_p);

      return s;
    }
} /* gsl_eigen_gensymm() */

/*
gsl_eigen_gensymm_standardize()
  Reduce the generalized symmetric-definite eigenproblem to
the standard symmetric eigenproblem by computing

C = L^{-1} A L^{-t}

where L L^t is the Cholesky decomposition of B

Inputs: A - (input/output) real symmetric matrix
        B - real symmetric, positive definite matrix in Cholesky form

Return: success

Notes: A is overwritten by L^{-1} A L^{-t}
*/

int
gsl_eigen_gensymm_standardize(gsl_matrix *A, const gsl_matrix *B)
{
  const size_t N = A->size1;
  size_t i;
  double a, b, c;

  for (i = 0; i < N; ++i)
    {
      /* update lower triangle of A(i:n, i:n) */

      a = gsl_matrix_get(A, i, i);
      b = gsl_matrix_get(B, i, i);
      a /= b * b;
      gsl_matrix_set(A, i, i, a);

      if (i < N - 1)
        {
          gsl_vector_view ai = gsl_matrix_subcolumn(A, i, i + 1, N - i - 1);
          gsl_matrix_view ma =
            gsl_matrix_submatrix(A, i + 1, i + 1, N - i - 1, N - i - 1);
          gsl_vector_const_view bi =
            gsl_matrix_const_subcolumn(B, i, i + 1, N - i - 1);
          gsl_matrix_const_view mb =
            gsl_matrix_const_submatrix(B, i + 1, i + 1, N - i - 1, N - i - 1);

          gsl_blas_dscal(1.0 / b, &ai.vector);

          c = -0.5 * a;
          gsl_blas_daxpy(c, &bi.vector, &ai.vector);

          gsl_blas_dsyr2(CblasLower, -1.0, &ai.vector, &bi.vector, &ma.matrix);

          gsl_blas_daxpy(c, &bi.vector, &ai.vector);

          gsl_blas_dtrsv(CblasLower,
                         CblasNoTrans,
                         CblasNonUnit,
                         &mb.matrix,
                         &ai.vector);
        }
    }

  return GSL_SUCCESS;
} /* gsl_eigen_gensymm_standardize() */

max_stars_repo_name: ruslankuzmin/julia


text:
#pragma once

#include <gsl/span>

#include <cstddef>

namespace imageview {

// Implementation of the PixelFormat concept for the 8-bit grayscale pixel
// format.
class PixelFormatGrayscale8 {
 public:
  using color_type = unsigned char;
  static constexpr int kBytesPerPixel = 1;

  constexpr color_type read(gsl::span<const std::byte, kBytesPerPixel> data) const;

  constexpr void write(const color_type& color, gsl::span<std::byte, kBytesPerPixel> data) const;
};

constexpr PixelFormatGrayscale8::color_type PixelFormatGrayscale8::read(
    gsl::span<const std::byte, kBytesPerPixel> data) const {
  return static_cast<color_type>(data[0]);
}

constexpr void PixelFormatGrayscale8::write(const color_type& color, gsl::span<std::byte, kBytesPerPixel> data) const {
  data[0] = static_cast<std::byte>(color);
}

}  // namespace imageview

max_stars_repo_name: alexanderbelous/imageview

max_issues_repo_name: alexanderbelous/imageview

max_forks_repo_name: alexanderbelous/imageview


text:
//          Copyright Jean Pierre Cimalando 2019.
// Distributed under the Boost Software License, Version 1.0.
//    (See accompanying file LICENSE.md or copy at
//          http://www.boost.org/LICENSE_1_0.txt)

#pragma once
#include "keystate.h"
#include "instruments/synth_fx.h"
#include <gsl/gsl>
#include <string>
#include <vector>
#include <bitset>

struct Player_Song_Metadata {
    std::string name;
    std::string author;
    std::vector<std::string> text;
    char format[32] = {};
    unsigned track_count = 0;
};

struct Player_State {
    Keyboard_State kb;
    unsigned repeat_mode = 0;
    double time_position = 0;
    double duration = 0;
    double tempo = 0;
    unsigned speed = 100;
    std::string file_path;
    Player_Song_Metadata song_metadata;
    std::bitset<16> channel_enabled;
    float audio_levels[10] {};
    int fx_parameters[Synth_Fx::Parameter_Count] {};
};

max_stars_repo_name: jpcima/smf-dsp

max_issues_repo_name: jpcima/smf-dsp

max_forks_repo_name: jpcima/smf-dsp


text:
/*
 * MIT License
 *
 * Copyright (c) 2017 Intel Corporation
 *
 * 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.
 *
 *
 * Abstract: misc util functions and macros
 */

#ifndef __UTIL_H__
#define __UTIL_H__

#include <algorithm>
#include <gsl/span>
#include <iterator>
#include <sstream>
#include <string>
#include <vector>

typedef gsl::span<const uint8_t> cbspan;

#ifndef block_round
#define block_round(ODD, BLK)                                                  \
    ((ODD) + (((BLK) - ((ODD) & ((BLK)-1))) & ((BLK)-1)))
#endif

template <class TContainer>
static inline bool starts_with(const TContainer& haystack,
                               const TContainer& needle)
{
    return haystack.size() >= needle.size() &&
           std::equal(needle.begin(), needle.end(), haystack.begin());
}

static inline void split(const std::string& s, char delim,
                         std::vector<std::string>& result)
{
    std::string::size_type prev_pos = 0, pos = 0;
    do
    {
        pos = s.find(delim, pos);
        if (pos == std::string::npos)
        {
            result.emplace_back(std::move(s.substr(prev_pos)));
        }
        else
        {
            result.emplace_back(std::move(s.substr(prev_pos, pos - prev_pos)));
            prev_pos = ++pos;
        }
    } while (pos != std::string::npos);
}

static inline std::vector<std::string> split(const std::string& s, char delim)
{
    std::vector<std::string> elems;
    split(s, delim, elems);
    return elems;
}

static constexpr uint64_t prime = 0x100000001B3ull;
static constexpr uint64_t basis = 0xCBF29CE484222325ull;

static constexpr uint64_t const_hash(const char* str,
                                     uint64_t last_value = basis)
{
    return *str ? const_hash(str + 1, (*str ^ last_value) * prime) : last_value;
}

static inline uint64_t hash(const std::string& str)
{
    uint64_t ret{basis};

    for (auto i : str)
    {
        ret ^= static_cast<uint64_t>(i);
        ret *= prime;
    }

    return ret;
}

#endif /* __UTIL_H__ */

max_stars_repo_name: Intel-BMC/mtd-util

max_issues_repo_name: Intel-BMC/mtd-util

max_forks_repo_name: Intel-BMC/mtd-util


text:
#include "benchmark.hpp"
#include <lapacke.h>

int svd(double *A, int m, int n, double *U, double *s, double *V_t) {
  const int lda = n;
  const int ldu = m;
  const int ldvt = n;
  const char jobu = 'A';
  const char jobvt = 'A';
  const int superb_size = (m < n) ? m : n;
  double *superb = malloc(sizeof(double) * (superb_size - 1));
  lapack_int retval = LAPACKE_dgesvd(LAPACK_ROW_MAJOR,
                                     jobu,
                                     jobvt,
                                     m,
                                     n,
                                     A,
                                     lda,
                                     s,
                                     U,
                                     ldu,
                                     V_t,
                                     ldvt,
                                     superb);
  if (retval > 0) {
    return -1;
  }

  /* Clean up */
  free(superb);

  return 0;
}

int main() {
  for (size_t k = 10; k <= 1000; k += 10) {
    size_t m = k;
    double *A = create_random_sq_matrix(m);
    double *U = malloc(sizeof(double) * m * m);
    double *d = malloc(sizeof(double) * m);
    double *V_t = malloc(sizeof(double) * m * m);

    sleep(0.1);
    struct timespec t = tic();
    svd(A, m, m, U, d, V_t);
    printf("matrix_size: %ld\t: dgesvd: %fs\n", m, toc(&t));
    free(A);
    free(U);
    free(d);
    free(V_t);
  }

  return 0;
}

max_stars_repo_name: daoran/proto

max_issues_repo_name: daoran/proto

max_forks_repo_name: daoran/proto


text:
/* err/strerror.c
 *
 * 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.
 */

#include <config.h>
#include <gsl/gsl_errno.h>

const char *
gsl_strerror (const int gsl_errno)
{
  switch (gsl_errno)
    {
    case GSL_SUCCESS:
      return "success" ;
    case GSL_FAILURE:
      return "failure" ;
    case GSL_CONTINUE:
      return "the iteration has not converged yet";
    case GSL_EDOM:
      return "input domain error" ;
    case GSL_ERANGE:
      return "output range error" ;
    case GSL_EFAULT:
      return "invalid pointer" ;
    case GSL_EINVAL:
      return "invalid argument supplied by user" ;
    case GSL_EFAILED:
      return "generic failure" ;
    case GSL_EFACTOR:
      return "factorization failed" ;
    case GSL_ESANITY:
      return "sanity check failed - shouldn't happen" ;
    case GSL_ENOMEM:
      return "malloc failed" ;
    case GSL_EBADFUNC:
      return "problem with user-supplied function";
    case GSL_ERUNAWAY:
      return "iterative process is out of control";
    case GSL_EMAXITER:
      return "exceeded max number of iterations" ;
    case GSL_EZERODIV:
      return "tried to divide by zero" ;
    case GSL_EBADTOL:
      return "specified tolerance is invalid or theoretically unattainable" ;
    case GSL_ETOL:
      return "failed to reach the specified tolerance" ;
    case GSL_EUNDRFLW:
      return "underflow" ;
    case GSL_EOVRFLW:
      return "overflow" ;
    case GSL_ELOSS:
      return "loss of accuracy" ;
    case GSL_EROUND:
      return "roundoff error" ;
    case GSL_EBADLEN:
      return "matrix/vector sizes are not conformant" ;
    case GSL_ENOTSQR:
      return "matrix not square" ;
    case GSL_ESING:
      return "singularity or extremely bad function behavior detected" ;
    case GSL_EDIVERGE:
      return "integral or series is divergent" ;
    case GSL_EUNSUP:
      return "the required feature is not supported by this hardware platform";
    case GSL_EUNIMPL:
      return "the requested feature is not (yet) implemented";
    case GSL_ECACHE:
      return "cache limit exceeded";
    case GSL_ETABLE:
      return "table limit exceeded";
    case GSL_ENOPROG:
      return "iteration is not making progress towards solution";
    case GSL_ENOPROGJ:
      return "jacobian evaluations are not improving the solution";
    case GSL_ETOLF:
      return "cannot reach the specified tolerance in F";
    case GSL_ETOLX:
      return "cannot reach the specified tolerance in X";
    case GSL_ETOLG:
      return "cannot reach the specified tolerance in gradient";
    case GSL_EOF:
      return "end of file";
    default:
      return "unknown error code" ;
    }
}


text:
#include <lapacke.h>
#include <stdio.h>

/* Auxiliary routine: printing a matrix */
void print_matrix_rowmajor( char* desc, lapack_int m, lapack_int n, double* mat, lapack_int ldm ) {
        lapack_int i, j;
        printf( "\n %s\n", desc );

        for( i = 0; i < m; i++ ) {
                for( j = 0; j < n; j++ ) printf( " %6.2f", mat[i*ldm+j] );
                printf( "\n" );
        }
}

/* Auxiliary routine: printing a matrix */
void print_matrix_colmajor( char* desc, lapack_int m, lapack_int n, double* mat, lapack_int ldm ) {
        lapack_int i, j;
        printf( "\n %s\n", desc );

        for( i = 0; i < m; i++ ) {
                for( j = 0; j < n; j++ ) printf( " %6.2f", mat[i+j*ldm] );
                printf( "\n" );
        }
}

/* Auxiliary routine: printing a vector of integers */
void print_vector( char* desc, lapack_int n, lapack_int* vec ) {
        lapack_int j;
        printf( "\n %s\n", desc );
        for( j = 0; j < n; j++ ) printf( " %6i", vec[j] );
        printf( "\n" );
}

max_stars_repo_name: dnoan/OpenBLAS

max_issues_repo_name: dnoan/OpenBLAS

max_forks_repo_name: dnoan/OpenBLAS


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

#include <config.h>
#include <stdlib.h>
#include <math.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_fft_halfcomplex.h>

#include "hc_pass.h"

int
FUNCTION(gsl_fft_halfcomplex,backward) (BASE data[], const size_t stride,
                                        const size_t n,
                                        const TYPE(gsl_fft_halfcomplex_wavetable) * wavetable,
                                        TYPE(gsl_fft_real_workspace) * work)
{
  int status = FUNCTION(gsl_fft_halfcomplex,transform) (data, stride, n, wavetable, work) ;
  return status ;
}

int
FUNCTION(gsl_fft_halfcomplex,inverse) (BASE data[], const size_t stride,
                                       const size_t n,
                                       const TYPE(gsl_fft_halfcomplex_wavetable) * wavetable,
                                       TYPE(gsl_fft_real_workspace) * work)
{
  int status = FUNCTION(gsl_fft_halfcomplex,transform) (data, stride, n, wavetable, work);

  if (status)
    {
      return status;
    }

  /* normalize inverse fft with 1/n */

  {
    const double norm = 1.0 / n;
    size_t i;
    for (i = 0; i < n; i++)
      {
        data[stride*i] *= norm;
      }
  }
  return status;
}

int
FUNCTION(gsl_fft_halfcomplex,transform) (BASE data[], const size_t stride, const size_t n,
                                         const TYPE(gsl_fft_halfcomplex_wavetable) * wavetable,
                                         TYPE(gsl_fft_real_workspace) * work)
{
  BASE * const scratch = work->scratch;

  BASE * in;
  BASE * out;
  size_t istride, ostride ;

  size_t factor, product, q;
  size_t i;
  size_t nf;
  int state;
  int product_1;
  int tskip;
  TYPE(gsl_complex) *twiddle1, *twiddle2, *twiddle3, *twiddle4;

  if (n == 0)
    {
      GSL_ERROR ("length n must be positive integer", GSL_EDOM);
    }

  if (n == 1)
    {                           /* FFT of one data point is the identity */
      return 0;
    }

  if (n != wavetable->n)
    {
      GSL_ERROR ("wavetable does not match length of data", GSL_EINVAL);
    }

  if (n != work->n)
    {
      GSL_ERROR ("workspace does not match length of data", GSL_EINVAL);
    }

  nf = wavetable->nf;
  product = 1;
  state = 0;

  for (i = 0; i < nf; i++)
    {
      factor = wavetable->factor[i];
      product_1 = product;
      product *= factor;
      q = n / product;

      tskip = (q + 1) / 2 - 1;

      if (state == 0)
        {
          in = data;
          istride = stride;
          out = scratch;
          ostride = 1;
          state = 1;
        }
      else
        {
          in = scratch;
          istride = 1;
          out = data;
          ostride = stride;
          state = 0;
        }

      if (factor == 2)
        {
          twiddle1 = wavetable->twiddle[i];
          FUNCTION(fft_halfcomplex,pass_2) (in, istride, out, ostride,
                                            product, n, twiddle1);
        }
      else if (factor == 3)
        {
          twiddle1 = wavetable->twiddle[i];
          twiddle2 = twiddle1 + tskip;
          FUNCTION(fft_halfcomplex,pass_3) (in, istride, out, ostride,
                                            product, n, twiddle1, twiddle2);
        }
      else if (factor == 4)
        {
          twiddle1 = wavetable->twiddle[i];
          twiddle2 = twiddle1 + tskip;
          twiddle3 = twiddle2 + tskip;
          FUNCTION(fft_halfcomplex,pass_4) (in, istride, out, ostride,
                                            product, n, twiddle1, twiddle2,
                                            twiddle3);
        }
      else if (factor == 5)
        {
          twiddle1 = wavetable->twiddle[i];
          twiddle2 = twiddle1 + tskip;
          twiddle3 = twiddle2 + tskip;
          twiddle4 = twiddle3 + tskip;
          FUNCTION(fft_halfcomplex,pass_5) (in, istride, out, ostride,
                                            product, n, twiddle1, twiddle2,
                                            twiddle3, twiddle4);
        }
      else
        {
          twiddle1 = wavetable->twiddle[i];
          FUNCTION(fft_halfcomplex,pass_n) (in, istride, out, ostride,
                                            factor, product, n, twiddle1);
        }
    }

  if (state == 1)               /* copy results back from scratch to data */
    {
      for (i = 0; i < n; i++)
        {
          data[stride*i] = scratch[i] ;
        }
    }

  return 0;

}