text:
/* matrix/gsl_matrix_complex_float.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 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 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.
 */

#ifndef __GSL_MATRIX_COMPLEX_FLOAT_H__
#define __GSL_MATRIX_COMPLEX_FLOAT_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector_complex_float.h>

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

__BEGIN_DECLS

typedef struct
{
  size_t size1;
  size_t size2;
  size_t tda;
  float * data;
  gsl_block_complex_float * block;
  int owner;
} gsl_matrix_complex_float ;

typedef struct
{
  gsl_matrix_complex_float matrix;
} _gsl_matrix_complex_float_view;

typedef _gsl_matrix_complex_float_view gsl_matrix_complex_float_view;

typedef struct
{
  gsl_matrix_complex_float matrix;
} _gsl_matrix_complex_float_const_view;

typedef const _gsl_matrix_complex_float_const_view gsl_matrix_complex_float_const_view;

/* Allocation */

GSL_FUN gsl_matrix_complex_float *
gsl_matrix_complex_float_alloc (const size_t n1, const size_t n2);

GSL_FUN gsl_matrix_complex_float *
gsl_matrix_complex_float_calloc (const size_t n1, const size_t n2);

GSL_FUN gsl_matrix_complex_float *
gsl_matrix_complex_float_alloc_from_block (gsl_block_complex_float * b,
                                           const size_t offset,
                                           const size_t n1, const size_t n2, const size_t d2);

GSL_FUN gsl_matrix_complex_float *
gsl_matrix_complex_float_alloc_from_matrix (gsl_matrix_complex_float * b,
                                            const size_t k1, const size_t k2,
                                            const size_t n1, const size_t n2);

GSL_FUN gsl_vector_complex_float *
gsl_vector_complex_float_alloc_row_from_matrix (gsl_matrix_complex_float * m,
                                                const size_t i);

GSL_FUN gsl_vector_complex_float *
gsl_vector_complex_float_alloc_col_from_matrix (gsl_matrix_complex_float * m,
                                                const size_t j);

GSL_FUN void gsl_matrix_complex_float_free (gsl_matrix_complex_float * m);

/* Views */

GSL_FUN _gsl_matrix_complex_float_view
gsl_matrix_complex_float_submatrix (gsl_matrix_complex_float * m,
                            const size_t i, const size_t j,
                            const size_t n1, const size_t n2);

GSL_FUN _gsl_vector_complex_float_view
gsl_matrix_complex_float_row (gsl_matrix_complex_float * m, const size_t i);

GSL_FUN _gsl_vector_complex_float_view
gsl_matrix_complex_float_column (gsl_matrix_complex_float * m, const size_t j);

GSL_FUN _gsl_vector_complex_float_view
gsl_matrix_complex_float_diagonal (gsl_matrix_complex_float * m);

GSL_FUN _gsl_vector_complex_float_view
gsl_matrix_complex_float_subdiagonal (gsl_matrix_complex_float * m, const size_t k);

GSL_FUN _gsl_vector_complex_float_view
gsl_matrix_complex_float_superdiagonal (gsl_matrix_complex_float * m, const size_t k);

GSL_FUN _gsl_vector_complex_float_view
gsl_matrix_complex_float_subrow (gsl_matrix_complex_float * m,
                                 const size_t i, const size_t offset,
                                 const size_t n);

GSL_FUN _gsl_vector_complex_float_view
gsl_matrix_complex_float_subcolumn (gsl_matrix_complex_float * m,
                                    const size_t j, const size_t offset,
                                    const size_t n);

GSL_FUN _gsl_matrix_complex_float_view
gsl_matrix_complex_float_view_array (float * base,
                             const size_t n1,
                             const size_t n2);

GSL_FUN _gsl_matrix_complex_float_view
gsl_matrix_complex_float_view_array_with_tda (float * base,
                                      const size_t n1,
                                      const size_t n2,
                                      const size_t tda);

GSL_FUN _gsl_matrix_complex_float_view
gsl_matrix_complex_float_view_vector (gsl_vector_complex_float * v,
                              const size_t n1,
                              const size_t n2);

GSL_FUN _gsl_matrix_complex_float_view
gsl_matrix_complex_float_view_vector_with_tda (gsl_vector_complex_float * v,
                                       const size_t n1,
                                       const size_t n2,
                                       const size_t tda);

GSL_FUN _gsl_matrix_complex_float_const_view
gsl_matrix_complex_float_const_submatrix (const gsl_matrix_complex_float * m,
                                  const size_t i, const size_t j,
                                  const size_t n1, const size_t n2);

GSL_FUN _gsl_vector_complex_float_const_view
gsl_matrix_complex_float_const_row (const gsl_matrix_complex_float * m,
                            const size_t i);

GSL_FUN _gsl_vector_complex_float_const_view
gsl_matrix_complex_float_const_column (const gsl_matrix_complex_float * m,
                               const size_t j);

GSL_FUN _gsl_vector_complex_float_const_view
gsl_matrix_complex_float_const_diagonal (const gsl_matrix_complex_float * m);

GSL_FUN _gsl_vector_complex_float_const_view
gsl_matrix_complex_float_const_subdiagonal (const gsl_matrix_complex_float * m,
                                    const size_t k);

GSL_FUN _gsl_vector_complex_float_const_view
gsl_matrix_complex_float_const_superdiagonal (const gsl_matrix_complex_float * m,
                                      const size_t k);

GSL_FUN _gsl_vector_complex_float_const_view
gsl_matrix_complex_float_const_subrow (const gsl_matrix_complex_float * m,
                                       const size_t i, const size_t offset,
                                       const size_t n);

GSL_FUN _gsl_vector_complex_float_const_view
gsl_matrix_complex_float_const_subcolumn (const gsl_matrix_complex_float * m,
                                          const size_t j, const size_t offset,
                                          const size_t n);

GSL_FUN _gsl_matrix_complex_float_const_view
gsl_matrix_complex_float_const_view_array (const float * base,
                                   const size_t n1,
                                   const size_t n2);

GSL_FUN _gsl_matrix_complex_float_const_view
gsl_matrix_complex_float_const_view_array_with_tda (const float * base,
                                            const size_t n1,
                                            const size_t n2,
                                            const size_t tda);

GSL_FUN _gsl_matrix_complex_float_const_view
gsl_matrix_complex_float_const_view_vector (const gsl_vector_complex_float * v,
                                    const size_t n1,
                                    const size_t n2);

GSL_FUN _gsl_matrix_complex_float_const_view
gsl_matrix_complex_float_const_view_vector_with_tda (const gsl_vector_complex_float * v,
                                             const size_t n1,
                                             const size_t n2,
                                             const size_t tda);

/* Operations */

GSL_FUN void gsl_matrix_complex_float_set_zero (gsl_matrix_complex_float * m);
GSL_FUN void gsl_matrix_complex_float_set_identity (gsl_matrix_complex_float * m);
GSL_FUN void gsl_matrix_complex_float_set_all (gsl_matrix_complex_float * m, gsl_complex_float x);

GSL_FUN int gsl_matrix_complex_float_fread (FILE * stream, gsl_matrix_complex_float * m) ;
GSL_FUN int gsl_matrix_complex_float_fwrite (FILE * stream, const gsl_matrix_complex_float * m) ;
GSL_FUN int gsl_matrix_complex_float_fscanf (FILE * stream, gsl_matrix_complex_float * m);
GSL_FUN int gsl_matrix_complex_float_fprintf (FILE * stream, const gsl_matrix_complex_float * m, const char * format);

GSL_FUN int gsl_matrix_complex_float_memcpy(gsl_matrix_complex_float * dest, const gsl_matrix_complex_float * src);
GSL_FUN int gsl_matrix_complex_float_swap(gsl_matrix_complex_float * m1, gsl_matrix_complex_float * m2);

GSL_FUN int gsl_matrix_complex_float_swap_rows(gsl_matrix_complex_float * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_complex_float_swap_columns(gsl_matrix_complex_float * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_complex_float_swap_rowcol(gsl_matrix_complex_float * m, const size_t i, const size_t j);

GSL_FUN int gsl_matrix_complex_float_transpose (gsl_matrix_complex_float * m);
GSL_FUN int gsl_matrix_complex_float_transpose_memcpy (gsl_matrix_complex_float * dest, const gsl_matrix_complex_float * src);

GSL_FUN int gsl_matrix_complex_float_isnull (const gsl_matrix_complex_float * m);
GSL_FUN int gsl_matrix_complex_float_ispos (const gsl_matrix_complex_float * m);
GSL_FUN int gsl_matrix_complex_float_isneg (const gsl_matrix_complex_float * m);
GSL_FUN int gsl_matrix_complex_float_isnonneg (const gsl_matrix_complex_float * m);

GSL_FUN int gsl_matrix_complex_float_add (gsl_matrix_complex_float * a, const gsl_matrix_complex_float * b);
GSL_FUN int gsl_matrix_complex_float_sub (gsl_matrix_complex_float * a, const gsl_matrix_complex_float * b);
GSL_FUN int gsl_matrix_complex_float_mul_elements (gsl_matrix_complex_float * a, const gsl_matrix_complex_float * b);
GSL_FUN int gsl_matrix_complex_float_div_elements (gsl_matrix_complex_float * a, const gsl_matrix_complex_float * b);
GSL_FUN int gsl_matrix_complex_float_scale (gsl_matrix_complex_float * a, const gsl_complex_float x);
GSL_FUN int gsl_matrix_complex_float_add_constant (gsl_matrix_complex_float * a, const gsl_complex_float x);
GSL_FUN int gsl_matrix_complex_float_add_diagonal (gsl_matrix_complex_float * a, const gsl_complex_float x);

/***********************************************************************/
/* The functions below are obsolete                                    */
/***********************************************************************/
GSL_FUN int gsl_matrix_complex_float_get_row(gsl_vector_complex_float * v, const gsl_matrix_complex_float * m, const size_t i);
GSL_FUN int gsl_matrix_complex_float_get_col(gsl_vector_complex_float * v, const gsl_matrix_complex_float * m, const size_t j);
GSL_FUN int gsl_matrix_complex_float_set_row(gsl_matrix_complex_float * m, const size_t i, const gsl_vector_complex_float * v);
GSL_FUN int gsl_matrix_complex_float_set_col(gsl_matrix_complex_float * m, const size_t j, const gsl_vector_complex_float * v);
/***********************************************************************/

/* inline functions if you are using GCC */

GSL_FUN INLINE_DECL gsl_complex_float gsl_matrix_complex_float_get(const gsl_matrix_complex_float * m, const size_t i, const size_t j);
GSL_FUN INLINE_DECL void gsl_matrix_complex_float_set(gsl_matrix_complex_float * m, const size_t i, const size_t j, const gsl_complex_float x);

GSL_FUN INLINE_DECL gsl_complex_float * gsl_matrix_complex_float_ptr(gsl_matrix_complex_float * m, const size_t i, const size_t j);
GSL_FUN INLINE_DECL const gsl_complex_float * gsl_matrix_complex_float_const_ptr(const gsl_matrix_complex_float * m, const size_t i, const size_t j);

#ifdef HAVE_INLINE

INLINE_FUN
gsl_complex_float
gsl_matrix_complex_float_get(const gsl_matrix_complex_float * m,
                     const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      gsl_complex_float zero = {{0,0}};

      if (i >= m->size1)
        {
          GSL_ERROR_VAL("first index out of range", GSL_EINVAL, zero) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_VAL("second index out of range", GSL_EINVAL, zero) ;
        }
    }
#endif
  return *(gsl_complex_float *)(m->data + 2*(i * m->tda + j)) ;
}

INLINE_FUN
void
gsl_matrix_complex_float_set(gsl_matrix_complex_float * m,
                     const size_t i, const size_t j, const gsl_complex_float x)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_VOID("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_VOID("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  *(gsl_complex_float *)(m->data + 2*(i * m->tda + j)) = x ;
}

INLINE_FUN
gsl_complex_float *
gsl_matrix_complex_float_ptr(gsl_matrix_complex_float * m,
                             const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  return (gsl_complex_float *)(m->data + 2*(i * m->tda + j)) ;
}

INLINE_FUN
const gsl_complex_float *
gsl_matrix_complex_float_const_ptr(const gsl_matrix_complex_float * m,
                                   const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  return (const gsl_complex_float *)(m->data + 2*(i * m->tda + j)) ;
}

#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_MATRIX_COMPLEX_FLOAT_H__ */

ext: h

max_stars_repo_name: berkus/music-cs

max_issues_repo_name: berkus/music-cs

max_forks_repo_name: berkus/music-cs


text:
#ifndef __DUMMY_H__
#define __DUMMY_H__
#include<stdio.h>
//#define NDEBUG
#include<stdint.h>
#include<omp.h>
#include<float.h>
#include<assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <time.h>
#include <math.h>
#include<string.h>
#include<dirent.h>
#include<sys/stat.h>
#include<sys/param.h>
#include<sys/types.h>
#include<signal.h>
//GSL library
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_blas.h>
#include<gsl/gsl_matrix.h>
#include <utility>
#include <string>
#include <iostream>
#include<iterator>
#include<map>
#include<vector>
#include <cctype>
#include<algorithm>
#include<climits>
using namespace std ;

const uint16_t loop_max=500;
uint16_t loop=0;
const double  epsilon= 0.05; //error associated with user variable !
const double  epsilon_m=    0.0000001; //error associated with user variable !
const double epsilon_phi=   0.0000001; //error associated with user variable !
const double epsilon_theta= 0.0000001; //error associated with user variable !
const double lower_limit_theta= 0.0000000001;
const double lower_limit_m =    0.0000000001;
const double lower_limit_phi =    0.000000001;
clock_t  start;
clock_t  end;
gsl_rng *rand_obj ;

double maxima_phi_i_possible=0.00000000000000000000000000000000000001;
double log_Fx_maxima_phi_i_possible=-DBL_MAX;

double * diversity;
uint32_t *freq_single_den; //  freqency of theta_w in denominator

typedef struct {
  bool* word_caps; // keeps tag if word at this pos was
  bool* comma_present_before;
  bool* period_present_before;
  vector<string> word_list ;
  uint32_t *int_word_list ; // Its an array to words
  uint32_t  *topic_usage;
  int real_word_count ; //Number of words in document
  int true_label ;
  int sampled_label;
  uint8_t *topic_list ; //Topic associated with each word
  char* document_name ;
  bool* priority_flag ;
  bool  hasPriorityWord;
  int mygamma;
}document ;

#define topic_count 5   // 1 to 30 !
#define document_count 790
#define total_words_count 10000

document *d_ptr ;
//documet
//Number of topics defined in begining
int updated_doc=0 ;

//GIBBS  samples count
//
int gibbs_counter=0;
const int gibbs_sample_count =25 ;
double *M_gibbs  ;
double *THETA_gibbs  ;
double *PHI_gibbs ;
void init_gibbs_params();
void write_model_gibbs();
void do_gibbs_expectation();
double alpha_comma_Y_gibbs=0;
double alpha_comma_N_gibbs=0;
double alpha_period_Y_gibbs=0;
double alpha_period_N_gibbs=0;
double alpha_upper_Y_gibbs=0;
double alpha_upper_N_gibbs=0;

double alpha_scale=1;
double alpha_shape=1;

const double  log_ratio=log(1000);

double *M ;
double *THETA ;
double *PHI;

//File names
const char* M_file="M.bin";
const char* THETA_file="THETA.bin";
const char* PHI_file="PHI.bin";
const char* TOPIC_file="topics_list.bin";

//For keeping last maxima , this will  speed up  search
double *last_maxima_theta_iw;
double *last_maxima_m_ij;
double *last_maxima_phi ;

//Thread info
const  int m_thread_cnt=8;
const  int theta_thread_cnt=8;

//Priors
const double m_diag=500;
const double m_nondiag=1;
const double m_beta=1;

const double theta_k=1;
const double theta_beta=1;
/*
double phi_sigma_square=10000;
double phi_sigma_square=1;
double phi_mu=0;
*/
//FUNCTIONS
//INIT  functions
void init(); //wrapper   for  all the inits
void init_model_parameters() ; //Inits model params
void init_document_params(); //INIT  document structure
void init_variables();
void init_topic_assignment(document* current_doc);
void test_model();
void free_things();

//Topic update for each document
inline void  update_topic_usage(document* current_doc_ptr, const uint8_t  new_topic,const int old_topic  );
void topic_update(document* current_doc_ptr);
void update_TOPIC();
double get_factor12(const int word_curr,const int word_next,const int t_prime,const int t_last,const int t_next,document* doc,const uint32_t mypos);
//double get_factor12(const uint32_t word_curr,const int word_next,const uint8_t t_prime,const int t_last,const int t_next,document* doc,const uint32_t mypos);
inline double get_factor2_t_i(document* current_doc_ptr, double sum);
double get_current_sum(document* doc,const uint32_t N_d, int sz);
inline void scale_pdf(double* pdf ,const int size,double sum );
inline void make_cdf_from_pdf(double* pdf ,double* cdf , const int size );

//update m_i_j
void  update_M();
void  update_m_ij(const int row,const int col,double * sum_den,uint32_t* freq_den,uint32_t*  freq_num,uint32_t* ref);
void get_max_f_m_ij(const uint8_t i, const uint8_t j,double*log_F_max,double* maxima,double*sum_den,uint32_t*freq_den,uint32_t* freq_num,uint32_t* ref,double * maxima_possible,double * val_maxima);
double evaluate_m_ij_ll(double x0,const uint8_t i,const uint8_t j,double* sum_den,uint32_t* freq_den,uint32_t* freq_num,uint32_t* ref, double *maxima_m_ij_possible,double *log_Fx_maxima_m_ij_possible);
void m_compute_double_term_word_stat(const uint8_t i,double* sum_den ,uint32_t* freq_den,uint32_t* freq_num_i,uint32_t* w_topic_count_map);

//UPDATE THETA matrix
void update_THETA();
void  update_theta_iw(const uint8_t i,const uint32_t w, double* den_single_term_header,double * sum_den,uint32_t* freq_den,const uint32_t freq_single_num,uint32_t *  freq_num,const uint32_t freq_single_den_w,uint32_t * ref);
void get_max_f_theta_iw(const uint8_t i, const uint32_t w,double*log_F_max,double* maxima, double* den_single_term_header,double * sum_den,uint32_t* freq_den,const uint32_t freq_single_num,uint32_t *  freq_num,const uint32_t freq_single_den_w, uint32_t* ref, double * maxima_possible,double *value_possible);
double evaluate_theta_iw_ll(double x0,const uint8_t i,const uint32_t w,const double den_single_term_header,double * sum_den,uint32_t * freq_den,const uint32_t freq_single_num,uint32_t* freq_num,const  uint32_t freq_single_den_w,uint32_t* ref , double *maxima_theta_iw_possible,double *log_Fx_maxima_theta_iw_possible);
double get_den_theta_jw(const uint8_t j , const uint32_t w);
double compute_theta_single_term_denominator(const uint32_t  w );
void theta_compute_double_term_word_stat(const uint32_t w,double* sum_den,uint32_t* freq_den,uint32_t* freq_single_num_w,uint32_t* freq_double_num,uint32_t* ref );

//Common codes
double inline  get_number_between_two_numbers(double low, double high);
inline double rand_0_1();
void ctrl_c_handler(int single_num);
void  write_model();
double timetaken(clock_t start,clock_t  end , const char* mesaage);
double gamma_approx_sampler(double r,double log_F_max,double log_f_r,double maxima);
void   guess_lebel(document* current_doc_ptr);
double m_diagonistics(const int i,const int j,double* sum_den,uint32_t* freq_den,uint32_t* freq_num , uint32_t* ref ,double *maxima_possible,double*value_possible,const char*msg , double maxima);

double theta_diagonistics(const int i, const int w,const double den_single_term_header,double * sum_den,uint32_t* freq_den,const uint32_t freq_single_num,uint32_t *  freq_num,const  uint32_t freq_single_den_w,uint32_t* ref, double *possible_maxima,double * val_at_max,const char* msg,double maxima);
int linear_scan(double*cdf , const int size,double data);
double inline get_d_2D_m(double *ptr, const int row,const int col,const int row_length);
void inline set_d_2D_m(double *ptr, const int row,const int col,double value,const int row_length);
uint32_t  get_i_3D_m(const uint32_t *ptr ,int x,int y,int z,int y_size,int z_size);
void  set_i_3D_m(uint32_t *ptr,int x,int y,int z,int value,int x_size,int y_size);
uint32_t inline get_i_2D_m(uint32_t *ptr, int row,int col,int row_length);
void inline set_i_2D_m(uint32_t *ptr, int row,int col,int value,int row_length);
void load_model_params();
void load_topics();

//update PHI
double  evaluate_phi_i(const double x0,const uint8_t i );
void update_PHI();
void   get_max_f_phi_i(const uint8_t i, double* log_Fx,double* maxima ) ; // calculate maximum of function of log !

void print_topic_corpus();
void print_model_params();
void print_system_config();

void reset_topic_usage();

uint32_t  get_key_3D_m(int x,int y,int z,int x_size,int y_size);
void print_time();

//Makes word  of document
//void make_int_words(document* doc);

void make_int_words(document* doc,vector<string>& word_list );

//These things cam because of new data interpretation
//Now theyi are in strings so,I will use these data structures for
//string interpretation.
void make_word_id(string word);
std::map<string , int > word_id_map;
std::map<int , string > id_word_map;
string  get_word_from_id(int id);

int  get_id_from_word(string word);
//while we read document we will silently populate  this map for furtehr uses !
void  create_word_id_id_word_map(string word);

//helper  function for day to word
bool isWordStartsWithCaps(string word);
bool isCommaBeforeMe(string word) ;
bool isPeriodBeforeMe(string word) ;

//New variables  for  \alpha_c_Y and so on !
double alpha_comma_Y=0;
double alpha_comma_N=0;
double alpha_period_Y=0;
double alpha_period_N=0;
double alpha_upper_Y=0;
double alpha_upper_N=0;
// shape  and scale of \alpha

//make string lower case
string make_lowercase_string(string word);
bool doesWordStartsUpper(string data);

const string highlight_tag_start="<:t:>";
const string highlight_tag_end="<:/t:>";
const string comma=",";
const string period=".";

double get_punctuation_multipler_num(document* doc, int pos,int t_last,int t_next );
double get_punctuation_multipler_den(document* doc, int pos,int t_last, int t_next );

#define  total_configs_punctuation 8
double multiplier_values[total_configs_punctuation];
void fill_multiplier_values();

bool  didReturnDBL_MAX=false;

//alpha_comma_Y
void sample_alpha_comma_Y();
void  make_alpha_comma_Y_statistics(int &num_count,double* sum_den,map<int,vector<int> >& ref);
double evaluate_alpha_comma_Y_log_posterior(double x0,const int &num_count, double* sum_den,const map<int , vector<int> > & ref,double *maxima_alpha_comma_Y_possible,double *log_Fx_maxima_alpha_comma_Y_possible);
void get_max_f_alpha_comma_Y(double* log_F_max,double *maxima ,int &num_count,double*sum_den,map<int,vector<int> >& ref,double * maxima_possible,double * val_maxima);
double alpha_comma_Y_diagonistics(int &num_count,double*sum_den,map<int,vector<int> >& ref,double * possible_maxima,double * val_at_max ,const char* msg,double maxima);

//alpha_comma_N
void sample_alpha_comma_N();
void  make_alpha_comma_N_statistics(int &num_count,double* sum_den,map<int,vector<int> >& ref);
double evaluate_alpha_comma_N_log_posterior(double x0,const int &num_count, double* sum_den,const map<int , vector<int> > & ref,double *maxima_alpha_comma_N_possible,double *log_Fx_maxima_alpha_comma_N_possible);
void get_max_f_alpha_comma_N(double* log_F_max,double *maxima ,int &num_count,double*sum_den,map<int,vector<int> >& ref,double * maxima_possible,double * val_maxima);
double alpha_comma_N_diagonistics(int &num_count,double*sum_den,map<int,vector<int> >& ref,double * possible_maxima,double * val_at_max ,const char* msg,double maxima);

//alpha_upper_Y
void sample_alpha_upper_Y();
void  make_alpha_upper_Y_statistics(int &num_count,double* sum_den,map<int,vector<int> >& ref);
double evaluate_alpha_upper_Y_log_posterior(double x0,const int &num_count, double* sum_den,const map<int , vector<int> > & ref,double *maxima_alpha_upper_Y_possible,double *log_Fx_maxima_alpha_upper_Y_possible);
void get_max_f_alpha_upper_Y(double* log_F_max,double *maxima ,int &num_count,double*sum_den,map<int,vector<int> >& ref,double * maxima_possible,double * val_maxima);
double alpha_upper_Y_diagonistics(int &num_count,double*sum_den,map<int,vector<int> >& ref,double * possible_maxima,double * val_at_max ,const char* msg,double maxima);

//alpha_upper_N
void sample_alpha_upper_N();
void  make_alpha_upper_N_statistics(int &num_count,double* sum_den,map<int,vector<int> >& ref);
double evaluate_alpha_upper_N_log_posterior(double x0,const int &num_count, double* sum_den,const map<int , vector<int> > & ref,double *maxima_alpha_upper_N_possible,double *log_Fx_maxima_alpha_upper_N_possible);
void get_max_f_alpha_upper_N(double* log_F_max,double *maxima ,int &num_count,double*sum_den,map<int,vector<int> >& ref,double * maxima_possible,double * val_maxima);
double alpha_upper_N_diagonistics(int &num_count,double*sum_den,map<int,vector<int> >& ref,double * possible_maxima,double * val_at_max ,const char* msg,double maxima);

//alpha_upper_Y
void sample_alpha_period_Y();
void  make_alpha_period_Y_statistics(int &num_count,double* sum_den,map<int,vector<int> >& ref);
double evaluate_alpha_period_Y_log_posterior(double x0,const int &num_count, double* sum_den,const map<int , vector<int> > & ref,double *maxima_alpha_period_Y_possible,double *log_Fx_maxima_alpha_period_Y_possible);
void get_max_f_alpha_period_Y(double* log_F_max,double *maxima ,int &num_count,double*sum_den,map<int,vector<int> >& ref,double * maxima_possible,double * val_maxima);
double alpha_period_Y_diagonistics(int &num_count,double*sum_den,map<int,vector<int> >& ref,double * possible_maxima,double * val_at_max ,const char* msg,double maxima);

//alpha_period_N
void sample_alpha_period_N();
void  make_alpha_period_N_statistics(int &num_count,double* sum_den,map<int,vector<int> >& ref);
double evaluate_alpha_period_N_log_posterior(double x0,const int &num_count, double* sum_den,const map<int , vector<int> > & ref,double *maxima_alpha_period_N_possible,double *log_Fx_maxima_alpha_period_N_possible);
void get_max_f_alpha_period_N(double* log_F_max,double *maxima ,int &num_count,double*sum_den,map<int,vector<int> >& ref,double * maxima_possible,double * val_maxima);
double alpha_period_N_diagonistics(int &num_count,double*sum_den,map<int,vector<int> >& ref,double * possible_maxima,double * val_at_max ,const char* msg,double maxima);

///home/rkumar2/final_test/data
///data1/rahul/projects/text_mining/branch/punctuation/chris_data
const char* data="/data1/rahul/projects/text_mining/branch/punctuation/chris_data/AB" ;
const char* dict_file="/data1/rahul/projects/text_mining/branch/punctuation/chris_data/word_list.txt" ;
#endif


text:
//
//  Revision.hpp
//  PretendToWork
//
//  Created by tqtifnypmb on 14/12/2017.
//  Copyright © 2017 tqtifnypmb. All rights reserved.
//

#pragma once

#include "../rope/Range.h"
#include "../types.h"

#include <gsl/gsl>
#include <vector>
#include <string>

namespace brick
{

class Rope;
class Revision {
public:
    enum class Operation {
        insert,
        erase,
    };

    Revision(size_t authorId, size_t revId, Operation op, const Range& range);
    Revision(size_t authorId, size_t revId, Operation op, const Range& range, const detail::CodePointList& text);
    Revision(const Revision&) = default;
    Revision() = default;

    void apply(gsl::not_null<Rope*> rope) const;
    bool canApply(gsl::not_null<const Rope*> rope) const;

    size_t authorId() const {
        return authorId_;
    }

    bool prior(const Revision& rev) const {
        return authorId() <= rev.authorId();
    }

    size_t revId() const {
        return revId_;
    }

    const Range& range() const {
        return range_;
    }

    Range& range() {
        return range_;
    }

    int affectLength() const {
        if (op() == Operation::insert) {
            return static_cast<int>(cplist().size());
        } else {
            return range().length;
        }
    }

    Operation op() const {
        return op_;
    }

    bool valid() const {
        return !range_.empty() && range_.location >= 0;
    }

    void setInvalid() {
        range_.length = 0;
    }

    const detail::CodePointList& cplist() const {
        return cplist_;
    }

private:

    size_t authorId_;

    // revId_ only identify individual revision
    // newer revision doesn't have to have greater revId_
    size_t revId_;
    detail::CodePointList cplist_;
    Operation op_;
    Range range_;
};

}   // namespace brick

max_stars_repo_name: tqtifnypmb/brick

max_issues_repo_name: tqtifnypmb/brick

max_forks_repo_name: tqtifnypmb/brick


text:
/* histogram/init2d.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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

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

gsl_histogram2d *
gsl_histogram2d_alloc (const size_t nx, const size_t ny)
{
  gsl_histogram2d *h;

  if (nx == 0)
    {
      GSL_ERROR_VAL ("histogram2d length nx must be positive integer",
			GSL_EDOM, 0);
    }

  if (ny == 0)
    {
      GSL_ERROR_VAL ("histogram2d length ny must be positive integer",
			GSL_EDOM, 0);
    }

  h = (gsl_histogram2d *) malloc (sizeof (gsl_histogram2d));

  if (h == 0)
    {
      GSL_ERROR_VAL ("failed to allocate space for histogram2d struct",
			GSL_ENOMEM, 0);
    }

  h->xrange = (double *) malloc ((nx + 1) * sizeof (double));

  if (h->xrange == 0)
    {
      free (h);		/* exception in constructor, avoid memory leak */

      GSL_ERROR_VAL ("failed to allocate space for histogram2d x ranges",
			GSL_ENOMEM, 0);
    }

  h->yrange = (double *) malloc ((ny + 1) * sizeof (double));

  if (h->yrange == 0)
    {
      free (h->xrange);
      free (h);		/* exception in constructor, avoid memory leak */

      GSL_ERROR_VAL ("failed to allocate space for histogram2d y ranges",
			GSL_ENOMEM, 0);
    }

  h->bin = (double *) malloc (nx * ny * sizeof (double));

  if (h->bin == 0)
    {
      free (h->xrange);
      free (h->yrange);
      free (h);		/* exception in constructor, avoid memory leak */

      GSL_ERROR_VAL ("failed to allocate space for histogram bins",
			GSL_ENOMEM, 0);
    }

  h->nx = nx;
  h->ny = ny;

  return h;
}

gsl_histogram2d *
gsl_histogram2d_calloc_uniform (const size_t nx, const size_t ny,
				const double xmin, const double xmax,
				const double ymin, const double ymax)
{
  gsl_histogram2d *h;

  if (xmin >= xmax)
    {
      GSL_ERROR_VAL ("xmin must be less than xmax", GSL_EINVAL, 0);
    }

  if (ymin >= ymax)
    {
      GSL_ERROR_VAL ("ymin must be less than ymax", GSL_EINVAL, 0);
    }

  h = gsl_histogram2d_calloc (nx, ny);

  if (h == 0)
    {
      return h;
    }

  {
    size_t i;

    for (i = 0; i < nx + 1; i++)
      {
	h->xrange[i] = xmin + ((double) i / (double) nx) * (xmax - xmin);
      }

    for (i = 0; i < ny + 1; i++)
      {
	h->yrange[i] = ymin + ((double) i / (double) ny) * (ymax - ymin);
      }
  }

  return h;
}

gsl_histogram2d *
gsl_histogram2d_calloc (const size_t nx, const size_t ny)
{
  gsl_histogram2d *h;

  if (nx == 0)
    {
      GSL_ERROR_VAL ("histogram2d length nx must be positive integer",
			GSL_EDOM, 0);
    }

  if (ny == 0)
    {
      GSL_ERROR_VAL ("histogram2d length ny must be positive integer",
			GSL_EDOM, 0);
    }

  h = (gsl_histogram2d *) malloc (sizeof (gsl_histogram2d));

  if (h == 0)
    {
      GSL_ERROR_VAL ("failed to allocate space for histogram2d struct",
			GSL_ENOMEM, 0);
    }

  h->xrange = (double *) malloc ((nx + 1) * sizeof (double));

  if (h->xrange == 0)
    {
      free (h);		/* exception in constructor, avoid memory leak */

      GSL_ERROR_VAL ("failed to allocate space for histogram2d x ranges",
			GSL_ENOMEM, 0);
    }

  h->yrange = (double *) malloc ((ny + 1) * sizeof (double));

  if (h->yrange == 0)
    {
      free (h->xrange);
      free (h);		/* exception in constructor, avoid memory leak */

      GSL_ERROR_VAL ("failed to allocate space for histogram2d y ranges",
			GSL_ENOMEM, 0);
    }

  h->bin = (double *) malloc (nx * ny * sizeof (double));

  if (h->bin == 0)
    {
      free (h->xrange);
      free (h->yrange);
      free (h);		/* exception in constructor, avoid memory leak */

      GSL_ERROR_VAL ("failed to allocate space for histogram bins",
			GSL_ENOMEM, 0);
    }

  {
    size_t i;

    for (i = 0; i < nx + 1; i++)
      {
	h->xrange[i] = i;
      }

    for (i = 0; i < ny + 1; i++)
      {
	h->yrange[i] = i;
      }

    for (i = 0; i < nx * ny; i++)
      {
	h->bin[i] = 0;
      }
  }

  h->nx = nx;
  h->ny = ny;

  return h;
}

void
gsl_histogram2d_free (gsl_histogram2d * h)
{
  free (h->xrange);
  free (h->yrange);
  free (h->bin);
  free (h);
}

int
gsl_histogram2d_set_ranges_uniform (gsl_histogram2d * h,
                                    double xmin, double xmax,
                                    double ymin, double ymax)
{
  size_t i;
  const size_t nx = h->nx, ny = h->ny;

  if (xmin >= xmax)
    {
      GSL_ERROR_VAL ("xmin must be less than xmax", GSL_EINVAL, 0);
    }

  if (ymin >= ymax)
    {
      GSL_ERROR_VAL ("ymin must be less than ymax", GSL_EINVAL, 0);
    }

  /* initialize ranges */

  for (i = 0; i <= nx; i++)
    {
      h->xrange[i] = xmin + ((double) i / (double) nx) * (xmax - xmin);
    }

  for (i = 0; i <= ny; i++)
    {
      h->yrange[i] = ymin + ((double) i / (double) ny) * (ymax - ymin);
    }

  /* clear contents */

  for (i = 0; i < nx * ny; i++)
    {
      h->bin[i] = 0;
    }

  return GSL_SUCCESS;
}

int
gsl_histogram2d_set_ranges (gsl_histogram2d * h,
                            const double xrange[], size_t xsize,
                            const double yrange[], size_t ysize)
{
  size_t i;
  const size_t nx = h->nx, ny = h->ny;

  if (xsize != (nx + 1))
    {
      GSL_ERROR_VAL ("size of xrange must match size of histogram",
                     GSL_EINVAL, 0);
    }

  if (ysize != (ny + 1))
    {
      GSL_ERROR_VAL ("size of yrange must match size of histogram",
                     GSL_EINVAL, 0);
    }

  /* initialize ranges */

  for (i = 0; i <= nx; i++)
    {
      h->xrange[i] = xrange[i];
    }

  for (i = 0; i <= ny; i++)
    {
      h->yrange[i] = yrange[i];
    }

  /* clear contents */

  for (i = 0; i < nx * ny; i++)
    {
      h->bin[i] = 0;
    }

  return GSL_SUCCESS;
}

ext: c


text:
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <mpfr.h>
#include <petsc.h>

#include "../ellipsoid/ellipsoid.h"
#include "../ellipsoid/ellSolv.h"
#include "../sphere/sphere.h"
#include "../sphere/sphSolv.h"
#include "../constants.h"

#undef __FUNCT__
#define __FUNCT__ "CoulombExact"
PetscErrorCode CoulombExact(PetscReal eps1, Vec srcXYZ, Vec srcMag, Vec tarXYZ, Vec sol)
{
  PetscErrorCode ierr;
  PetscReal srcX, srcY, srcZ;
  PetscReal tarX, tarY, tarZ;
  PetscReal dx2, dy2, dz2;
  PetscReal val, mag;
  const PetscScalar *srcXYZArray;
  const PetscScalar *srcMagArray;
  const PetscScalar *tarXYZArray;
  PetscScalar *solArray;
  PetscInt nSrc, nTar;
  PetscInt i, k;
  PetscFunctionBeginUser;

  ierr = VecGetSize(srcMag, &nSrc);CHKERRQ(ierr);
  ierr = VecGetSize(sol, &nTar);CHKERRQ(ierr);

  ierr = VecGetArrayRead(srcXYZ, &srcXYZArray);CHKERRQ(ierr);
  ierr = VecGetArrayRead(tarXYZ, &tarXYZArray);CHKERRQ(ierr);
  ierr = VecGetArrayRead(srcMag, &srcMagArray);CHKERRQ(ierr);
  ierr = VecGetArray    (sol   , &solArray   );CHKERRQ(ierr);
  for(i=0; i<nTar; ++i) {
    /* get target X,Y,Z */
    tarX = tarXYZArray[3*i+0];
    tarY = tarXYZArray[3*i+1];
    tarZ = tarXYZArray[3*i+2];
    for(k=0; k<nSrc; ++k) {
      /* get source X,Y,Z */
      srcX = srcXYZArray[3*k+0];
      srcY = srcXYZArray[3*k+1];
      srcZ = srcXYZArray[3*k+2];
      mag  = srcMagArray[k];

      dx2 = (tarX-srcX)*(tarX-srcX);
      dy2 = (tarY-srcY)*(tarY-srcY);
      dz2 = (tarZ-srcZ)*(tarZ-srcZ);

      val = mag/(eps1*PetscSqrtReal(dx2 + dy2 + dz2));
      solArray[i] += val;
    }
  }
  ierr = VecRestoreArrayRead(srcXYZ, &srcXYZArray);CHKERRQ(ierr);
  ierr = VecRestoreArrayRead(tarXYZ, &tarXYZArray);CHKERRQ(ierr);
  ierr = VecRestoreArrayRead(srcMag, &srcMagArray);CHKERRQ(ierr);
  ierr = VecRestoreArray    (sol   , &solArray   );CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

max_stars_repo_name: tom-klotz/ellipsoid-solvation

max_issues_repo_name: tom-klotz/ellipsoid-solvation

max_forks_repo_name: tom-klotz/ellipsoid-solvation


text:
#include <gsl/gsl_odeiv2.h>
#include "func.h"
#include "make_grid.h"

int
make_grid (struct param params, double r, double r1, double y[])
{
  int status;
  int *fake_jac;
  /* The integrator RK8PD (Runge-Kutta Prince-Dormand) doesn't need the
   * jacobian (some more sophisticated integrators do) so make this a
   * null pointer.  The linker will complain that the pointer type is
   * incompatable since it doesn't point to a function with a bunch of
   * arguments, but the pointer won't be used anyway so it doesn't
   * matter. */
  fake_jac = 0;

  gsl_odeiv2_system sys = { func, fake_jac, 2, &params };
  gsl_odeiv2_driver *d = gsl_odeiv2_driver_alloc_y_new
    (&sys, gsl_odeiv2_step_rk8pd,
     1.0e+0, 1.0e-8, 1.0e-8);

  status = gsl_odeiv2_driver_apply (d, &r, r1, y);

  gsl_odeiv2_driver_free (d);
  return 0;
}


text:
/**
 *
 * @precisions normal z -> c d s
 *
 **/
#define _TYPE  PLASMA_Complex64_t
#define _PREC  double
#define _LAMCH LAPACKE_dlamch_work

#define _NAME  "PLASMA_zgeqrf"
/* See Lawn 41 page 120 */
#define _FMULS FMULS_GEQRF(M, N)
#define _FADDS FADDS_GEQRF(M, N)

#include "./timing.c"
#include <cblas.h>

#define COMPLEX

/* computes the residual ||A*P - Q*R||
 * similar to zqpt01.f in testing.
 * A    is m x n
 * AF   is m x n, partial QR factorization of k columns of A.
 * tau  for Householder transformations.
 * jpvt column pivots.
 * work array, dimension (lwork)
 * lwork >= m*n + n.
 */
double qp_residual(
    int m, int n, int k,
    const PLASMA_Complex64_t *A,  int lda,
    const PLASMA_Complex64_t *AF, int ldaf,
    const PLASMA_Complex64_t *tau,
    const int *jpvt,
    PLASMA_Complex64_t *work, int lwork )
{
    PLASMA_Complex64_t mzone = (PLASMA_Complex64_t)-1.0;
    PLASMA_Complex64_t zzero = (PLASMA_Complex64_t) 0.0;
    int j, ldw;
    double rwork[1];

    ldw = m;
    if ( lwork < ldw*n + n ) {
        fprintf( stderr, "Error: lwork %d too small, requires %d\n", lwork, m*n + n );
        return 0.;
    }

    /* copy R to work, with zeros below diagonal */
    LAPACKE_zlaset( LAPACK_COL_MAJOR, 'L', m-1, n, zzero, zzero, &work[1], ldw );
    LAPACKE_zlacpy( LAPACK_COL_MAJOR, 'U', m, n, AF, lda, work, ldw );

    /* form W = QR */
    LAPACKE_zunmqr_work( LAPACK_COL_MAJOR, 'L', 'N', m, n, k,
                         AF, ldaf, tau, work, ldw,
                         &work[n*ldw], lwork-n*ldw );

    //printf( "\n" );
    //printf( "m %d, n %d, k %d\n", m, n, k );
    //CORE_zprint( m, n, A,    lda,  "A"   );
    //CORE_zprint( m, n, AF,   ldaf, "QR" );
    //CORE_zprint( 1, min(m,n), tau, 1, "tau" );
    //CORE_zprint( m, n, work, ldw,  "Q*R" );

    /* compare (subtract) j-th column of W to jpvt(j)-th column of A (jpvt is 1-based) */
    for( j = 0; j < n; ++j ) {
        cblas_zaxpy( m, CBLAS_SADDR(mzone), &A[(jpvt[j]-1)*lda], 1, &work[j*ldw], 1 );
    }

    //CORE_zprint( m, n, work, ldw,  "Q*R - AP" );

    return LAPACKE_zlange_work( LAPACK_COL_MAJOR, 'O', m, n, work, ldw, rwork );
}

static int
RunTest(int *iparam, double *dparam, real_Double_t *t_)
{
    int info, lwork, lwork2, i;
    int *jpvt;
    PLASMA_Complex64_t *tau, *work, *work2;
    double *rwork;

    PASTE_CODE_IPARAM_LOCALS( iparam );

    /* Allocate Data */
    PASTE_CODE_ALLOCATE_MATRIX( A, 1, PLASMA_Complex64_t, LDA, N );

    /* LAPACK needs larger size in real;
     * PLASMA currently uses work and rwork for both real and complex. */
    #ifdef COMPLEX
    lwork = (N+1)*NB;
    #else
    lwork = (N+1)*NB + 2*N;
    #endif
    work  = (PLASMA_Complex64_t*) malloc( lwork * sizeof(PLASMA_Complex64_t) );
    rwork = (double*)             malloc( 2*2*N * sizeof(double)             );
    jpvt  = (int*)                malloc( N     * sizeof(int)                );
    tau   = (PLASMA_Complex64_t*) malloc( 2*N   * sizeof(PLASMA_Complex64_t) );

    if ( jpvt == NULL || tau == NULL || work == NULL || rwork == NULL ) {
        fprintf( stderr, "malloc failed\n" );
        return -1;
    }

    /* zero out pivots (required by LAPACK) */
    for( i = 0; i < N; ++i ) {
        jpvt[i] = 0;
    }

    /* Initialize Data */
    PLASMA_zplrnt(M, N, A, LDA, 123456);

    /* Save A in lapack layout for check */
    PASTE_CODE_ALLOCATE_COPY( Acpy, check, PLASMA_Complex64_t, A, LDA, N );

    START_TIMING();
    info = PLASMA_zgeqp3( M, N, A, LDA, jpvt, tau, work, rwork );
    STOP_TIMING();

    /* Check the solution */
    if ( info != 0 ) {
        printf( "\nPLASMA_zgeqp3 returned error %d.\n", info );
    }
    else if ( check ) {
        lwork2 = (M*N + N);
        work2  = (PLASMA_Complex64_t*) malloc( lwork2 * sizeof(PLASMA_Complex64_t) );
        if ( work2 == NULL ) {
            fprintf( stderr, "test malloc failed\n" );
            return -1;
        }

        dparam[IPARAM_ANORM] = LAPACKE_zlange_work( LAPACK_COL_MAJOR, 'F', N, N, Acpy, LDA, rwork );
        dparam[IPARAM_XNORM] = 1.;
        dparam[IPARAM_BNORM] = 0.;
        dparam[IPARAM_RES]   = qp_residual( M, N, min(M,N), Acpy, LDA,
                                            A, LDA, tau, jpvt, work2, lwork2 );

        // /* compute result with LAPACK */
        // int *jpvt2;
        // PLASMA_Complex64_t *tau2;
        // jpvt2  = (int*)                malloc( N * sizeof(int)    );
        // tau2   = (PLASMA_Complex64_t*) malloc( N * sizeof(PLASMA_Complex64_t) );
        // if ( jpvt2 == NULL || tau2 == NULL ) {
        //     fprintf( stderr, "test malloc failed\n" );
        //     return -1;
        // }
        //
        // /* zero out pivots (required by LAPACK) */
        // for( i = 0; i < N; ++i ) {
        //     jpvt2[i] = 0;
        // }
        //
        // double time = cWtime();
        // #ifdef COMPLEX
        // info = LAPACKE_zgeqp3_work( LAPACK_COL_MAJOR, N, N, Acpy, LDA, jpvt2, tau2, work, lwork, rwork );
        // #else
        // info = LAPACKE_zgeqp3_work( LAPACK_COL_MAJOR, N, N, Acpy, LDA, jpvt2, tau2, work, lwork );
        // #endif
        // time = cWtime() - time;
        // if ( info != 0 ) {
        //     printf( "qp3 returned error %d\n", info );
        // }
        // /* printf( "   %7.3f", time ); */
        //
        // CORE_zprint( M, N, Acpy, LDA, "QR_L" );
        // CORE_zprint( 1, min(M,N), tau2, 1, "tau_L" );

        // cblas_zaxpy( LDA*N, CBLAS_SADDR(mzone), A, 1, Acpy, 1 );
        // dparam[IPARAM_RES] = LAPACKE_zlange_work( LAPACK_COL_MAJOR, 'F', N, N, Acpy, LDA, rwork );
        //
        // cblas_zaxpy( N, CBLAS_SADDR(mzone), tau, 1, tau2, 1 );
        // tnorm = cblas_dznrm2( N, tau2, 1 );
        /* printf( "|t|=%.2e\n", tnorm ); */

        // for( i = 0; i < N; ++i ) {
        //     if ( jpvt[i]+1 != jpvt2[i] ) {
        //         printf( "pivot mis-match jpvt[%2d]+1=%2d, jpvt2[%2d]=%2d\n", i, jpvt[i]+1, i, jpvt2[i] );
        //     }
        // }

        // free( jpvt2 );
        // free( tau2  );

        free( Acpy  );
        free( work2 );
    }

    /* Free data */
    free( A     );
    free( jpvt  );
    free( tau   );
    free( work  );
    free( rwork );

    return 0;
}


text:
/*
 * File:   rng.h
 * Author: fedro
 *
 * Created on November 3, 2014, 5:33 PM
 */

#ifndef RNG_H
#define	RNG_H

#include <gsl/gsl_rng.h>

class RNG
{
public:
    RNG(const gsl_rng_type* rng_type, const int seed=0)
    {
        rng_ = gsl_rng_alloc(rng_type);
        gsl_rng_set(rng_, seed);
    }
    RNG(const int seed=0)
    {
        rng_ = gsl_rng_alloc(gsl_rng_default);
        gsl_rng_set(rng_, seed);
    }
    virtual ~RNG()
    {
        gsl_rng_free(rng_);
    }
    gsl_rng* get() { return rng_; }
private:
    gsl_rng* rng_;
};

#endif	/* RNG_H */


text:
/* -----------------------------------------------------------------------------
 * Copyright 2021 Jonathan Haigh
 * SPDX-License-Identifier: MIT
 * ---------------------------------------------------------------------------*/

#ifndef SQ_INCLUDE_GUARD_system_linux_SqTypeSchemaImpl_h_
#define SQ_INCLUDE_GUARD_system_linux_SqTypeSchemaImpl_h_

#include "core/typeutil.h"
#include "system/SqTypeSchema.gen.h"

#include "system/schema.h"

#include <gsl/gsl>

namespace sq::system::linux {

class SqTypeSchemaImpl : public SqTypeSchema<SqTypeSchemaImpl> {
public:
  explicit SqTypeSchemaImpl(const TypeSchema &type_schema);

  SQ_ND Result get_name() const;
  SQ_ND Result get_doc() const;
  SQ_ND Result get_fields() const;

  SQ_ND Primitive to_primitive() const override;

private:
  gsl::not_null<const TypeSchema *> type_schema_;
};

} // namespace sq::system::linux

#endif // SQ_INCLUDE_GUARD_system_linux_SqTypeSchemaImpl_h_

max_stars_repo_name: jonathanhaigh/sq

max_issues_repo_name: jonathanhaigh/sq

max_forks_repo_name: jonathanhaigh/sq


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

#pragma once

#ifndef res_92f78d87_1896_40fa_bc63_c52e2793c970_h
#define res_92f78d87_1896_40fa_bc63_c52e2793c970_h

#include <gslib/config.h>
#include <gslib/std.h>

__gslib_begin__

class resvdir
{
public:
    resvdir(const gchar* path);
    resvdir(const gchar* dir, const gchar* file);
    const gchar* get_dir() const { return _dir.c_str(); }
    const gchar* get_file() const { return _file.c_str(); }
    const gchar* get_postfix() const;
    const gchar* get_path(string& path) const;

public:
    static bool is_valid_path(const gchar* path);
    static bool is_valid_dir(const gchar* dir);
    static const gchar* get_root(const gchar* path);
    static const gchar* get_next_node(const gchar* path);

protected:
    string      _dir;
    string      _file;
};

typedef unsigned char res;

class resnode
{
public:
    struct hash
    {
        size_t operator()(const resnode& that) const
        { return string_hash(that._key); }
    };
    struct equal_to
    {
        bool operator()(const resnode& r1, const resnode& r2) const
        { return string_hash(r1._key) == string_hash(r2._key); }
    };
    typedef unordered_set<resnode, hash, equal_to> resset;
    typedef resset::iterator iterator;
    typedef resset::const_iterator const_iterator;
    typedef void (*fndestroy)(res*, int);

public:
    resnode();
    resnode(const gchar* key);
    resnode(const resnode& that) { assign(const_cast<resnode&>(that)); }
    ~resnode() { destroy(); }
    void set_destroy(fndestroy del) { _del = del; }
    void set_resource(res* ptr, int len) { _res = ptr; _len = len; }
    res* get_resource() const { return _res; }
    const gchar* get_key() const { return _key.size() ? _key.c_str() : 0; }
    void destroy();
    void assign(resnode& that);
    bool operator<(const resnode& that) const { return _key < that._key; }
    iterator begin() { return _subs.begin(); }
    const_iterator begin() const { return _subs.end(); }
    iterator end() { return _subs.end(); }
    const_iterator end() const { return _subs.end(); }
    resnode* add(const gchar* key, res* ptr = 0, int len = 0);
    resnode* add(resnode* node);
    resnode* find(const gchar* key);
    const resnode* find(const gchar* key) const;

protected:
    string      _key;
    res*        _res;
    int         _len;
    fndestroy   _del;
    resset      _subs;
};

class respack
{
public:
    typedef resnode::resset resset;
    typedef resnode::iterator iterator;
    typedef resnode::const_iterator const_iterator;
    typedef resnode::fndestroy fndestroy;

public:
    respack(): _root(_t(":")) {}
    resnode* get_root() { return &_root; }
    const resnode* get_root() const { return &_root; }
    resnode* add(const gchar* path) { return add(resvdir(path)); }
    resnode* add(const resvdir& vdir) { return add(get_root(), resvdir::get_root(vdir.get_dir()), vdir.get_file()); }
    resnode* add(resnode* node, const gchar* dir, const gchar* file);
    resnode* find(const gchar* path) { return find(get_root(), resvdir::get_root(path)); }
    resnode* find(resnode* node, const gchar* path);
    const resnode* find(const gchar* path) const { return find(get_root(), resvdir::get_root(path)); }
    const resnode* find(const resnode* node, const gchar* path) const;
    resnode* reg(const gchar* path, res* ptr, int len, fndestroy del = 0)
    {
        if(resnode* node = add(path)) {
            node->set_resource(ptr, len);
            node->set_destroy(del);
            return node;
        }
        return 0;
    }

protected:
    resnode     _root;
};

extern respack* get_default_respack();

__gslib_end__

#endif

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


text:
/*
Copyright 2010-2011, D. E. Shaw Research.
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:

* Redistributions of source code must retain the above copyright
  notice, this list of conditions, and the following disclaimer.

* Redistributions in binary form must reproduce the above copyright
  notice, this list of conditions, and the following disclaimer in the
  documentation and/or other materials provided with the distribution.

* Neither the name of D. E. Shaw Research nor the names of its
  contributors may be used to endorse or promote products derived from
  this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
#ifndef __r123_gslmicrorng_dot_h__
#define __r123_gslmicrorng_dot_h__

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

/**   The macro: GSL_MICRORNG(NAME, CBRNGNAME) is the GSL
   analog analog of the C++ r123::MicroURNG template.  It declares a gsl_rng
   type named gsl_rng_NAME which uses the underlying CBRNGNAME
   and can be invoked a limited number of times between calls to NAME_reset.

   When the underlying CBRNG's \c ctr_t is an \ref arrayNxW "r123arrayNxW",
   and the gsl_rng_NAME may called up to \c N*2^32 times
   between calls to \c NAME_reset.

   \c NAME_reset takes a gsl_rng_NAME type, a counter and a key as arguments.
   It restarts the micro-rng with a new base counter and key.

   Note that you must call NAME_reset before the first use
   of a gsl_rng.  NAME_reset is not called automatically by
   gsl_rng_alloc().

   @code
   #include <Random123/threefry.h>
   #include <Random123/gsl_microrng.h> // this file
   GSL_MICRORNG(microcbrng, threefry4x64, 20)	// creates gsl_rng_microcbrng

   int main(int argc, char** argv) {
	gsl_rng *r = gsl_rng_alloc(gsl_rng_microcbrng);
	threefry4x64_ctr_t c = {{}};
	threefry4x64_key_t k = {{}};

	for (...) {
	    c.v[0] = ??; //  some application variable
	    microcbrng_reset(r, c, k);
	    for (...) {
		// gaussian calls r several times.  It is safe for
		// r to be used upto 2^20 times in this loop
		something[i] = gsl_ran_gaussian(r, 1.5);
	    }
	}
   }
   @endcode

*/

#define GSL_MICRORNG(NAME, CBRNGNAME)                                   \
const gsl_rng_type *gsl_rng_##NAME;                                     \
                                                                        \
typedef struct{                                                         \
    CBRNGNAME##_ctr_t ctr;                                              \
    CBRNGNAME##_ctr_t r;                                                \
    CBRNGNAME##_key_t key;                                              \
    R123_ULONG_LONG n;                                                  \
    int elem;                                                           \
} NAME##_state;                                                         \
                                                                        \
static unsigned long int NAME##_get(void *vstate){                      \
    NAME##_state *st = (NAME##_state *)vstate;                          \
    const int N=sizeof(st->ctr.v)/sizeof(st->ctr.v[0]);                 \
    if( st->elem == 0 ){                                                \
        CBRNGNAME##_ctr_t c = st->ctr;                                  \
        c.v[N-1] |= st->n<<(R123_W(CBRNGNAME##_ctr_t)-32);              \
        st->n++;                                                        \
        st->r = CBRNGNAME(c, st->key);                                  \
        st->elem = N;                                                   \
    }                                                                   \
    return 0xffffffff & st->r.v[--st->elem];                            \
}                                                                       \
                                                                        \
static double                                                           \
NAME##_get_double (void * vstate)                                       \
{                                                                       \
    return NAME##_get (vstate)/4294967296.;                             \
}                                                                       \
                                                                        \
static void NAME##_set(void *vstate, unsigned long int s){              \
    NAME##_state *st = (NAME##_state *)vstate;                          \
    st->elem = 0;                                                       \
    st->n = ~0; /* will abort if _reset is not called */                \
}                                                                       \
                                                                        \
static const gsl_rng_type NAME##_type = {                               \
    #NAME,                                                              \
    0xffffffffUL,                                                       \
    0,                                                                  \
    sizeof(NAME##_state),                                               \
    &NAME##_set,                                                        \
    &NAME##_get,                                                        \
    &NAME##_get_double                                                  \
};                                                                      \
                                                                        \
R123_STATIC_INLINE void NAME##_reset(const gsl_rng* gr, CBRNGNAME##_ctr_t c, CBRNGNAME##_key_t k) { \
    NAME##_state* state = (NAME##_state *)gr->state;                    \
    state->ctr = c;                                                     \
    state->key = k;                                                     \
    state->n = 0;                                                       \
    state->elem = 0;                                                    \
}                                                                       \
                                                                        \
const gsl_rng_type *gsl_rng_##NAME = &NAME##_type

#endif

max_stars_repo_name: havogt/philox-wrapper

max_issues_repo_name: havogt/philox-wrapper

max_forks_repo_name: havogt/philox-wrapper


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>

int generate_MCMC_chain(MCMC_info *MCMC) {
    int              flag_success;
    int              i_P, i_DS;
    static char      format_string[256];
    FILE *           fp_report_props;
    static char      filename_report_props[SID_MAX_FILENAME_LENGTH];
    static int       n_P;
    static double *  P_new;
    static double ** M_new;
    static double *  P_last;
    static double *  P_chain;
    static double ** M_last;
    static double *  P_best;
    static double *  P_limit_min;
    static double *  P_limit_max;
    static RNG_info *RNG;
    static int       flag_report_props;
    static int       n_DS;
    static int *     n_M;

    // Initialize a few things on the first call
    switch(MCMC->first_chain_call) {
        case GBP_TRUE:
            n_P                    = MCMC->n_P;
            P_new                  = MCMC->P_new;
            M_new                  = MCMC->M_new;
            P_last                 = MCMC->P_last;
            M_last                 = MCMC->M_last;
            P_chain                = MCMC->P_chain;
            P_best                 = MCMC->P_best;
            P_limit_min            = MCMC->P_limit_min;
            P_limit_max            = MCMC->P_limit_max;
            RNG                    = MCMC->RNG;
            n_DS                   = MCMC->n_DS;
            n_M                    = MCMC->n_M;
            flag_report_props      = SID_CHECK_BITFIELD_SWITCH(MCMC->mode, MCMC_MODE_REPORT_PROPS);
            MCMC->first_chain_call = GBP_FALSE;
            sprintf(filename_report_props, "%s/chains/report_props_%06d.dat", MCMC->filename_output_dir, MCMC->my_chain);
            break;
    }

    // We need to compute the likelihood of the initial conditions
    //   if we are just starting a chain
    if(MCMC->flag_init_chain)
        generate_MCMC_proposition(MCMC, GBP_TRUE);

    // Set the last new state to the new last state
    MCMC->ln_likelihood_last = MCMC->ln_likelihood_new;
    MCMC->ln_Pr_last         = MCMC->ln_Pr_new;
    memcpy(P_last, P_new, n_P * sizeof(double));
    for(i_DS = 0; i_DS < n_DS; i_DS++)
        memcpy(M_last[i_DS], M_new[i_DS], n_M[i_DS] * sizeof(double));

    // Keep generating parameter sets until the mapping function is satisfied
    generate_MCMC_proposition(MCMC, GBP_FALSE);

    // Decide if this is a successful proposition or not...
    if(MCMC->ln_likelihood_new > MCMC->ln_likelihood_chain) {
        MCMC->ln_Pr_new = 0.;
        flag_success    = GBP_TRUE;
    } else if(MCMC->my_chain == SID.My_rank) {
        MCMC->ln_Pr_new = MCMC->ln_likelihood_new - MCMC->ln_likelihood_chain;
        if((double)random_number(RNG) <= exp(MCMC->ln_Pr_new))
            flag_success = GBP_TRUE;
        else
            flag_success = GBP_FALSE;
    }
    // fprintf(stderr,"%d %le %d\n",flag_success,MCMC->ln_Pr_new,SID_CHECK_BITFIELD_SWITCH(MCMC->mode,MCMC_MODE_PARALLEL));
    // if(!SID_CHECK_BITFIELD_SWITCH(MCMC->mode,MCMC_MODE_PARALLEL))
    SID_Bcast(&flag_success, 1, SID_INT, SID_MASTER_RANK, MCMC->comm);

    // ... if it is, then update the chain ...
    if(flag_success) {
        memcpy(P_chain, P_new, (size_t)n_P * sizeof(double));
        MCMC->ln_likelihood_chain = MCMC->ln_likelihood_new;
        MCMC->ln_Pr_chain         = MCMC->ln_Pr_new;
        MCMC->n_success++;
    } else
        MCMC->n_fail++;
    MCMC->n_propositions++;

    // Report the proposition if asked to
    if(flag_report_props && MCMC->my_chain == SID.My_rank) {
        fp_report_props = fopen(filename_report_props, "w");
        if(flag_success)
            fprintf(fp_report_props, "Proposal #%09d: SUCCEEDED\n\n", MCMC->n_propositions);
        else
            fprintf(fp_report_props, "Proposal #%09d: REJECTED\n\n", MCMC->n_propositions);
        sprintf(format_string, "   %s = %%13.6le (was %%13.6le) (best is %%13.6le) (flat prior =%%13.6le -> %%13.6le)\n", MCMC->P_name_format);
        for(i_P = 0; i_P < n_P; i_P++)
            fprintf(fp_report_props, format_string, MCMC->P_names[i_P], P_new[i_P], P_last[i_P], P_best[i_P], P_limit_min[i_P], P_limit_max[i_P]);
        fprintf(fp_report_props, "\n");
        fprintf(fp_report_props,
                "ln(likelihood) =%13.6le+constant (was %13.6le+constant) (best is %13.6le+constant)\n",
                MCMC->ln_likelihood_new,
                MCMC->ln_likelihood_last,
                MCMC->ln_likelihood_best);
        fprintf(fp_report_props, "ln(probability)=%13.6le\n", MCMC->ln_Pr_new);
        fprintf(fp_report_props, "n_success      =%09d\n", MCMC->n_success);
        fprintf(fp_report_props, "n_fail         =%09d\n", MCMC->n_fail);
        fprintf(fp_report_props, "success rate   =%.2f %%\n", (double)MCMC->n_success / (double)(MCMC->n_fail + MCMC->n_success) * 100.0);
        fclose(fp_report_props);
    }

    return (flag_success);
}

max_stars_repo_name: gbpoole/gbpCode

max_issues_repo_name: gbpoole/gbpCode

max_forks_repo_name: gbpoole/gbpCode


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

int main(int argc, char *argv[]) {
    double mu = 1.0;
    if (argc > 1)
        mu = atof(argv[1]);
    int n = 1;
    if (argc > 2)
        n = atoi(argv[2]);
    unsigned long seed;
    if (argc > 3)
        seed = atol(argv[3]);
    else
        seed = time(NULL);
    gsl_rng_env_setup();
    gsl_rng *rng = gsl_rng_alloc(gsl_rng_default);
    gsl_rng_set(rng, seed);
    for (int i = 0; i < n; i++) {
        double x = gsl_ran_exponential(rng, mu);
        printf("%.12f\n", x);
    }
    gsl_rng_free(rng);
    return EXIT_SUCCESS;
}

max_stars_repo_name: Gjacquenot/training-material

max_issues_repo_name: Gjacquenot/training-material

max_forks_repo_name: Gjacquenot/training-material


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

/* #include <cblas.h> */
#include <lapacke_utils.h>

#include <jni.h>
#include "org_bitbucket_ioplus_openblas_Lapacke.h"

/*
 * Class:     org_bitbucket_ioplus_openblas_Lapacke
 * Method:    sgetrf
 * Signature: (II[FI[I)V
 */
    JNIEXPORT void JNICALL Java_org_bitbucket_ioplus_openblas_Lapacke_sgetrf
    (JNIEnv * env, jclass class,
     jint m, jint n, jfloatArray a, jint lda, jintArray ipiv){

        float *A;
        A = (*env)->GetPrimitiveArrayCritical(env, a, 0);
        int *IPIV;
        IPIV = (*env)->GetPrimitiveArrayCritical(env, ipiv, 0);

        LAPACKE_sgetrf(LAPACK_ROW_MAJOR,
                       m, n, A, lda, IPIV);

        (*env)->ReleasePrimitiveArrayCritical(env, a, A, 0);
        (*env)->ReleasePrimitiveArrayCritical(env, ipiv, IPIV, 0);

    }

/*
 * Class:     org_bitbucket_ioplus_openblas_Lapacke
 * Method:    sgetri
 * Signature: (I[FI[I)V
 */
    JNIEXPORT void JNICALL
    Java_org_bitbucket_ioplus_openblas_Lapacke_sgetri
    (JNIEnv * env, jclass class,
     jint n, jfloatArray a, jint lda, jintArray ipiv){

        float *A;
        A = (*env)->GetPrimitiveArrayCritical(env, a, 0);
        int *IPIV;
        IPIV = (*env)->GetPrimitiveArrayCritical(env, ipiv, 0);

        LAPACKE_sgetri(LAPACK_ROW_MAJOR,
                       n, A, lda, IPIV);

        (*env)->ReleasePrimitiveArrayCritical(env, a, A, 0);
        (*env)->ReleasePrimitiveArrayCritical(env, ipiv, IPIV, 0);

    }

/*
 * Class:     org_bitbucket_ioplus_openblas_Lapacke
 * Method:    sdeev
 * Signature: (CCI[FI[F[F[FI[FI)V
 */
JNIEXPORT void JNICALL Java_org_bitbucket_ioplus_openblas_Lapacke_sdeev
(JNIEnv * env, jclass class,
 jchar jobvl,
 jchar jobvr,
 jint n,
 jfloatArray a,
 jint lda,
 jfloatArray wr,
 jfloatArray wi,
 jfloatArray vl,
 jint ldvl,
 jfloatArray vr,
 jint ldvr){
    /* jni args -> c mapping  */

    float *A;
    float *WR;
    float *WI;
    float *VL;
    float *VR;

    A = (*env)->GetPrimitiveArrayCritical(env, a, 0);
    WR = (*env)->GetPrimitiveArrayCritical(env, wr, 0);
    WI = (*env)->GetPrimitiveArrayCritical(env, wi, 0);
    VL = (*env)->GetPrimitiveArrayCritical(env, vl, 0);
    VR = (*env)->GetPrimitiveArrayCritical(env, vr, 0);

    LAPACKE_sgeev(LAPACK_ROW_MAJOR,
                  jobvl,
                  jobvr,
                  n,
                  A,
                  lda,
                  WR, // values
                  WI,
                  VL,
                  ldvl,
                  VR, // vectors
                  ldvr);

    (*env)->ReleasePrimitiveArrayCritical(env, a, A, 0);
    (*env)->ReleasePrimitiveArrayCritical(env, wr, WR, 0);
    (*env)->ReleasePrimitiveArrayCritical(env, wi, WI, 0);
    (*env)->ReleasePrimitiveArrayCritical(env, vr, VR, 0);
    (*env)->ReleasePrimitiveArrayCritical(env, vl, VL, 0);

}

/*
 * Class:     org_bitbucket_ioplus_openblas_Lapacke
 * Method:    sgehrd
 * Signature: (III[FI[F)V
 */
JNIEXPORT void JNICALL Java_org_bitbucket_ioplus_openblas_Lapacke_sgehrd
(JNIEnv * env, jclass class,
 jint n,
 jint ilo,
 jint ihi,
 jfloatArray a,
 jint lda,
 jfloatArray tau
    ){
    float *A  = (*env)->GetPrimitiveArrayCritical(env, a, NULL);
    float *TAU = (*env)->GetPrimitiveArrayCritical(env, tau, NULL);
    LAPACKE_sgehrd(LAPACK_ROW_MAJOR,n, ilo, ihi, A, lda, TAU);
    (*env)->ReleasePrimitiveArrayCritical(env, a, A, 0);
    (*env)->ReleasePrimitiveArrayCritical(env, tau, TAU, 0);
}

/*
 * Class:     org_bitbucket_ioplus_openblas_Lapacke
 * Method:    shseqr
 * Signature: (CCIII[FI[F[F[FI)V
 */
JNIEXPORT void JNICALL Java_org_bitbucket_ioplus_openblas_Lapacke_shseqr
(JNIEnv * env, jclass class,
 jchar job,
 jchar compz,
 jint n,
 jint ilo,
 jint ihi,
 jfloatArray h,
 jint ldh,
 jfloatArray wr,
 jfloatArray wi,
 jfloatArray z,
 jint ldz){

    float *H = (*env)->GetPrimitiveArrayCritical(env, h, 0);
    float *WR = (*env)->GetPrimitiveArrayCritical(env, wr, 0);
    float *WI = (*env)->GetPrimitiveArrayCritical(env, wi, 0);
    float *Z = (*env)->GetPrimitiveArrayCritical(env, z, 0);

    LAPACKE_shseqr(LAPACK_ROW_MAJOR,
                   job,
                   compz,
                   n,
                   ilo,
                   ihi,
                   H,
                   ldh,
                   WR,
                   WI,
                   Z,
                   ldz);

    (*env)->ReleasePrimitiveArrayCritical(env, h, H, 0);
    (*env)->ReleasePrimitiveArrayCritical(env, wr, WR, 0);
    (*env)->ReleasePrimitiveArrayCritical(env, wi, WI, 0);
    (*env)->ReleasePrimitiveArrayCritical(env, z, Z, 0);

}

max_stars_repo_name: sidec/openblas-java

max_issues_repo_name: sidec/openblas-java

max_forks_repo_name: sidec/openblas-java


text:
/*
 * Copyright (C) 2018-2020
 * HERE Europe B.V.
 * Utrecht University (The Netherlands).
 *
 * 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.
 *
 * SPDX-License-Identifier: Apache-2.0
 * License-Filename: LICENSE
 */

//
// Created by Wouter Jongeling (wouter.jongeling@gmail.com)
// Modified by Mitra, Aniket on 2019-06-18.
//

/*!
 * @file BrownianBridge.h
 * @brief
 * @authors Aniket Mitra (aniket.mitra@here.com)
 *  Wouter Jongeling (wouter.jongeling@gmail.com)
 */

#ifndef MOVETK_BROWNIANBRIDGE__H
#define MOVETK_BROWNIANBRIDGE__H

#include <iterator>
#include <cmath>
#include "movetk/utils/Requirements.h"
#include "movetk/utils/Iterators.h"
#include "movetk/geom/trajectory_to_interface.h"
#include "movetk/utils/Iterators.h"
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include "gsl/gsl_min.h"

/*!
 * @file AlgorithmTraits.h
 * @brief The Brownian Bridge Movement Model (BBMM)
 * @authors Aniket Mitra (aniket.mitra@here.com)
 */

/*!
 *  @brief a collection of algorithms in movetk
 */
namespace movetk_algorithms {

    // based on
    // https://www.semanticscholar.org/paper/Computational-Movement-Analysis-Using-Brownian-Sijben/c6104979e609e12e1514b8b605a7d14fa6b5d81a

    /*!
    *  @brief a collection of classes for BBMM
    */
    namespace brownian_bridge {

        /*!
         * @struct ParameterTraits
         * @brief This traits class serves as a collection of types
         * for parameterization of BBMM
         * @tparam _GeometryTraits -  This class is a collection of movetk
         *  geometry types. For example @refitem movetk_core::MovetkGeometryKernel
         * @tparam TrajectoryIterator - An iterator type.
         * For example @refitem TabularTrajectory::TrajectoryIterator
         */
        template<class _GeometryTraits, class TrajectoryIterator>
        struct ParameterTraits {
            typedef _GeometryTraits GeometryTraits;
            /*!*
             * @typedef ::NT
             * @brief typedef  of the number type defined in GeometryTraits.
             * For example @refitem  movetk_core::MovetkGeometryKernel::NT
             * */
            typedef typename GeometryTraits::NT NT;
            /*!*
             * @typedef ::Point
             * @brief typedef of the point type defined in GeometryTraits.
             * For example @refitem movetk_core::MovetkGeometryKernel::MovetkPoint
             *  */
            typedef typename GeometryTraits::MovetkPoint Point;
            /*!* @typedef ::Vector
             * @brief typedef of the vector type defined in GeometryTraits.
             * For example @refitem  movetk_core::MovetkGeometryKernel::MovetkVector
             * */
            typedef typename GeometryTraits::MovetkVector Vector;
            /*!* @typedef ::Parameters
             */
            typedef std::tuple<Point, Point, NT, TrajectoryIterator, TrajectoryIterator> Parameters;
            /*!
             * @enum ParameterColumns
             *
             */
            enum ParameterColumns {
                POINT,  /*!< global enum value 1 */
                MU,     /*!< global enum value 1 */
                SIGMA_SQUARED,  /*!< global enum value 1 */
                BEGIN,  /*!< global enum value 1 */
                END    /*!< global enum value 1 */
            };

        };

        /*!
         * @class template<class GeometryTraits, class ParameterTraits,
         *  class Norm, class InputIterator, std::size_t MaxIter> MLE
         * @brief The Maximum Likelihood Estimator
         * @tparam GeometryTraits -   This class is a collection of movetk
         *  geometry types. For example @refitem movetk_core::MovetkGeometryKernel
         * @tparam ParameterTraits -  This traits class serves as a collection of types
         * for parameterization of BBMM
         * @tparam Norm - The type that models Euclidean distance
         * For example @refitem movetk_support::FiniteNorm
         * @tparam InputIterator - A random access iterator type over a collection of
         * tuples where each tuple is a model of type ParameterTraits::Parameters
         * @tparam MaxIter - A positive integer that bounds the number of iterations in the MLE
         * algorithm
         */
        template<class GeometryTraits, class ParameterTraits,
                class Norm, class InputIterator, std::size_t MaxIter,
                typename = movetk_core::requires_random_access_iterator<InputIterator>,
                typename = movetk_core::requires_tuple<
                        typename InputIterator::value_type>,
                typename = movetk_core::requires_tuple_element_as_movetk_point<GeometryTraits,
                        ParameterTraits::ParameterColumns::POINT,
                        typename InputIterator::value_type>,
                typename = movetk_core::requires_tuple_element_as_movetk_point<GeometryTraits,
                        ParameterTraits::ParameterColumns::MU,
                        typename InputIterator::value_type>,
                typename = movetk_core::requires_tuple_element_as_NT<GeometryTraits,
                        ParameterTraits::ParameterColumns::SIGMA_SQUARED,
                        typename InputIterator::value_type>>
        class MLE {
        private:
            typedef typename GeometryTraits::NT NT;
            const gsl_min_fminimizer_type *T = gsl_min_fminimizer_goldensection;
            gsl_min_fminimizer *s = gsl_min_fminimizer_alloc(T);
            gsl_function F;
            NT estimated_parameter;
            std::size_t iter = 0, status;

            struct fn_params {
                InputIterator first;
                InputIterator beyond;
            };

            static double fn(double sigma_squared, void *p) {
                Norm norm;
                struct fn_params *params = (fn_params *) p;
                InputIterator it = params->first;
                NT log_likelihood = 0;
                while (it != params->beyond) {
                    typename GeometryTraits::MovetkVector v = std::get<ParameterTraits::ParameterColumns::POINT>(*it)
                                                              - std::get<ParameterTraits::ParameterColumns::MU>(*it);
                    NT operand1 = -LOG_TWO_PI - log(static_cast<NT>(sigma_squared));
                    NT operand2 = -norm(v) / (2 * sigma_squared);
                    log_likelihood += operand1 + operand2;
                    it++;
                }
                return -1 * log_likelihood;
            }

            void operator()(InputIterator first, InputIterator beyond,
                            NT result, NT x_lower, NT x_upper, NT eps) {
                struct fn_params params;
                params.first = first;
                params.beyond = beyond;
                F.function = &MLE::fn;
                F.params = &params;
                status = gsl_min_fminimizer_set(s, &F, result, x_lower, x_upper);
                do {
                    iter++;
                    status = gsl_min_fminimizer_iterate(s);
                    result = gsl_min_fminimizer_x_minimum(s);
                    x_lower = gsl_min_fminimizer_x_lower(s);
                    x_upper = gsl_min_fminimizer_x_upper(s);

                    status = gsl_min_test_interval(x_lower, x_upper, eps, 0.0);

                    if (status == GSL_SUCCESS) {
                        estimated_parameter = result;
                        break;
                    }
                } while (status == GSL_CONTINUE && iter < MaxIter);

                if (iter >= MaxIter)
                    estimated_parameter = result;
            }

        public:
            /*!
             *
             * @param first
             * @param beyond
             * @param result
             * @param x_lower
             * @param x_upper
             * @param eps
             */
            MLE(InputIterator first, InputIterator beyond,
                NT result, NT x_lower, NT x_upper, NT eps) {
                Norm norm;
                std::size_t num_elements = std::distance(first, beyond);
                if (num_elements == 1) {
                    typename GeometryTraits::MovetkVector v = std::get<ParameterTraits::ParameterColumns::POINT>(*first)
                                                              - std::get<ParameterTraits::ParameterColumns::MU>(*first);
                    NT l = norm(v);
                    if (l < MOVETK_EPS) {
                        estimated_parameter = MOVETK_EPS;
                        return;
                    }
                }
                (*this)(first, beyond, result, x_lower, x_upper, eps);
            }

            /*!
             *
             * @param first
             * @param beyond
             */
            MLE(InputIterator first, InputIterator beyond) {
                Norm norm;
                InputIterator it = first;
                NT upper_bound = 0, squared_length = 0;
                std::size_t num_elements = std::distance(first, beyond);
                if (num_elements == 1) {
                    typename GeometryTraits::MovetkVector v = std::get<ParameterTraits::ParameterColumns::POINT>(*first)
                                                              - std::get<ParameterTraits::ParameterColumns::MU>(*first);
                    NT l = norm(v);
                    if (l < MOVETK_EPS) {
                        estimated_parameter = MOVETK_EPS;
                        return;
                    }
                }
                while (it != beyond) {
                    typename GeometryTraits::MovetkVector v = std::get<ParameterTraits::ParameterColumns::POINT>(*it)
                                                              - std::get<ParameterTraits::ParameterColumns::MU>(*it);
                    NT l = norm(v);
                    squared_length += l;
                    if (l > upper_bound)
                        upper_bound = l;
                    it++;
                }
                NT initial_estimate = squared_length / (2 * num_elements);
                (*this)(first, beyond, initial_estimate, MOVETK_EPS, upper_bound, MOVETK_EPS);
            }

            /*!
             *
             * @return
             */
            NT operator()() {
                return estimated_parameter;
            }

            /*!
             *
             */
            ~MLE() {
                gsl_min_fminimizer_free(s);
            }
        };

        /*!
         *
         * @tparam GeometryTraits
         * @tparam ProbeTraits
         * @tparam ParameterTraits
         * @tparam Norm
         * @tparam GeoProjection
         */
        template<class GeometryTraits, class ProbeTraits, class ParameterTraits, class Norm,
                class GeoProjection>
        class Model {
        private:

            typedef typename GeometryTraits::NT NT;
            typedef typename GeometryTraits::MovetkPoint Point;
            typedef typename GeometryTraits::MovetkVector Vector;
            movetk_core::MakePoint<GeometryTraits> make_point;

            NT joint_log_likelihood = 0;

            Point mean(Point &p1, Point &p2, NT alpha) {
                Vector v = p2 - p1;
                v *= alpha;
                return (p1 + v);
            }

        public:
            /*!
             *
             * @tparam TrajectoryIterator
             * @tparam OutputIterator
             * @param first
             * @param beyond
             * @param result
             */
            template<class TrajectoryIterator, class OutputIterator,
                    typename = movetk_core::requires_random_access_iterator<TrajectoryIterator>,
                    typename = movetk_core::requires_tuple<
                            typename TrajectoryIterator::value_type>,
                    typename = movetk_core::requires_tuple_element_as_arithmetic<
                            ProbeTraits::ProbeColumns::LAT, typename TrajectoryIterator::value_type>,
                    typename = movetk_core::requires_tuple_element_as_arithmetic<
                            ProbeTraits::ProbeColumns::LON, typename TrajectoryIterator::value_type>,
                    typename = movetk_core::requires_tuple_element_as_size_t<
                            ProbeTraits::ProbeColumns::SAMPLE_DATE, typename TrajectoryIterator::value_type>,
                    typename = movetk_core::requires_output_iterator<OutputIterator>,
                    typename = movetk_core::requires_tuple<
                            typename OutputIterator::value_type>,
                    typename = movetk_core::requires_tuple_element_as_movetk_point<GeometryTraits,
                            ParameterTraits::ParameterColumns::POINT,
                            typename OutputIterator::value_type>,
                    typename = movetk_core::requires_tuple_element_as_movetk_point<GeometryTraits,
                            ParameterTraits::ParameterColumns::MU,
                            typename OutputIterator::value_type>,
                    typename = movetk_core::requires_tuple_element_as_NT<GeometryTraits,
                            ParameterTraits::ParameterColumns::SIGMA_SQUARED,
                            typename OutputIterator::value_type>>
            Model(TrajectoryIterator first, TrajectoryIterator beyond, OutputIterator result) {

                auto reflat = std::get<ProbeTraits::ProbeColumns::LAT>(*first);
                auto reflon = std::get<ProbeTraits::ProbeColumns::LON>(*first);
                GeoProjection ref(reflat, reflon);

                TrajectoryIterator tit = first;
                std::size_t NumPoints = std::distance(first, beyond);

                if (NumPoints == 1) {
                    auto lat = std::get<ProbeTraits::ProbeColumns::LAT>(*first);
                    auto lon = std::get<ProbeTraits::ProbeColumns::LON>(*first);
                    auto projected_point = ref.project(lat, lon);
                    Point p1 = make_point(std::cbegin(projected_point), std::cend(projected_point));
                    *result = std::make_tuple(p1, p1, 0, first, first);
                    return;
                }

                std::size_t multiple = (NumPoints - 1) / 2;
                std::size_t remainder = (NumPoints - 1) % 2;
                TrajectoryIterator last = first + 2 * multiple;

                while (tit != last) {
                    auto lat = std::get<ProbeTraits::ProbeColumns::LAT>(*tit);
                    auto lon = std::get<ProbeTraits::ProbeColumns::LON>(*tit);
                    auto projected_point = ref.project(lat, lon);
                    Point p1 = make_point(std::cbegin(projected_point), std::cend(projected_point));

                    lat = std::get<ProbeTraits::ProbeColumns::LAT>(*(tit + 1));
                    lon = std::get<ProbeTraits::ProbeColumns::LON>(*(tit + 1));
                    projected_point = ref.project(lat, lon);
                    Point p2 = make_point(std::cbegin(projected_point), std::cend(projected_point));

                    lat = std::get<ProbeTraits::ProbeColumns::LAT>(*(tit + 2));
                    lon = std::get<ProbeTraits::ProbeColumns::LON>(*(tit + 2));
                    projected_point = ref.project(lat, lon);
                    Point p3 = make_point(std::cbegin(projected_point), std::cend(projected_point));

                    auto ts1 = std::get<ProbeTraits::ProbeColumns::SAMPLE_DATE>(*tit);
                    auto ts2 = std::get<ProbeTraits::ProbeColumns::SAMPLE_DATE>(*(tit + 1));
                    auto ts3 = std::get<ProbeTraits::ProbeColumns::SAMPLE_DATE>(*(tit + 2));
                    NT alpha = static_cast<NT>(  ts3 - ts2 ) / static_cast<NT>( ts3 - ts1 );

                    Point mu = mean(p1, p3, alpha);
                    *result = std::make_tuple(p2, mu, 0, tit, tit + 2);
                    tit += 2;
                }

                if (remainder == 1) {
                    auto lat = std::get<ProbeTraits::ProbeColumns::LAT>(*last);
                    auto lon = std::get<ProbeTraits::ProbeColumns::LON>(*last);
                    auto projected_point = ref.project(lat, lon);
                    Point p1 = make_point(std::cbegin(projected_point), std::cend(projected_point));

                    lat = std::get<ProbeTraits::ProbeColumns::LAT>(*(last + 1));
                    lon = std::get<ProbeTraits::ProbeColumns::LON>(*(last + 1));
                    projected_point = ref.project(lat, lon);
                    Point p2 = make_point(std::cbegin(projected_point), std::cend(projected_point));
                    Point mu = mean(p1, p2, 0.5);
                    *result = std::make_tuple(mu, mu, 0, last, last + 1);
                }

            }

        };

        /*!
         *
         * @tparam GeometryTraits
         * @tparam ParameterTraits
         */
        template<class GeometryTraits, class ParameterTraits>
        class ParameterSelector {
            std::size_t SIZE;
        public:
            /*!
             *
             * @param size
             */
            ParameterSelector(std::size_t size) : SIZE(size) {}

            /*!
             *
             * @tparam InputIterator
             * @tparam OutputIterator
             * @param first
             * @param beyond
             * @param result
             */
            template<class InputIterator, class OutputIterator,
                    typename = movetk_core::requires_random_access_iterator<InputIterator>,
                    typename = movetk_core::requires_output_iterator<OutputIterator>,
                    typename = movetk_core::requires_tuple<typename InputIterator::value_type>,
                    typename = movetk_core::requires_tuple_element_as_NT<GeometryTraits,
                            ParameterTraits::ParameterColumns::SIGMA_SQUARED,
                            typename InputIterator::value_type>,
                    typename = movetk_core::requires_NT<GeometryTraits,
                            typename OutputIterator::value_type >>
            void operator()(InputIterator first, InputIterator beyond, OutputIterator result) {
                assert(std::distance(first, beyond) >= SIZE);
                std::vector<typename ParameterTraits::NT> coeffs;
                std::transform(first, beyond, std::back_inserter(coeffs), [](auto a) {
                    return std::get<ParameterTraits::ParameterColumns::SIGMA_SQUARED>(a);
                });
                std::sort(std::begin(coeffs), std::end(coeffs));
                std::size_t coeff_interval = std::distance(first, beyond) / SIZE;
                for (std::size_t i = 0; i < SIZE; i++) {
                    std::size_t index = coeff_interval / 2 + coeff_interval * i;
                    *result = coeffs[index];
                }
            }
        };

        template<class GeometryTraits, class ParameterTraits, class Norm>
        class LogLikelihood {
            typedef typename GeometryTraits::NT NT;
            typedef typename ParameterTraits::Parameters Parameters;
        public:
            LogLikelihood() = default;

            template<class InputIterator, class OutputIterator,
                    typename = movetk_core::requires_random_access_iterator<InputIterator>,
                    typename = movetk_core::requires_output_iterator<OutputIterator>,
                    typename = movetk_core::requires_NT<GeometryTraits,
                            typename InputIterator::value_type>,
                    typename = movetk_core::requires_NT<GeometryTraits,
                            typename OutputIterator::value_type >>
            void operator()(const Parameters &params, InputIterator first,
                            InputIterator beyond, OutputIterator result) {
                Norm norm;
                typename GeometryTraits::MovetkVector v = std::get<ParameterTraits::ParameterColumns::POINT>(params)
                                                          - std::get<ParameterTraits::ParameterColumns::MU>(params);
                InputIterator pit = first;
                NT squared_length = norm(v);
                while (pit != beyond) {
                    NT operand1 = -LOG_TWO_PI - log(*pit);
                    NT operand2 = -squared_length / (2 * (*pit));
                    NT log_likelihood = operand1 + operand2;
                    *result = log_likelihood;
                    pit++;
                }
            }

            NT operator()(const Parameters &params, NT sigma_squared) {
                Norm norm;
                typename GeometryTraits::MovetkVector v = std::get<ParameterTraits::ParameterColumns::POINT>(params)
                                                          - std::get<ParameterTraits::ParameterColumns::MU>(params);
                NT squared_length = norm(v);
                NT operand1 = -LOG_TWO_PI - log(sigma_squared);
                NT operand2 = -squared_length / (2 * sigma_squared);
                return (operand1 + operand2);
            }
        };

    }

}

#endif //MOVETK_BROWNIANBRIDGE__H

max_stars_repo_name: hiteshbedre/movetk

max_issues_repo_name: bacusters/movetk


text:
/* rng/knuthran2002.c
 *
 * Copyright (C) 2007 Brian Gough
 * Copyright (C) 2001 Brian Gough, Carlo Perassi
 *
 * 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.
 */

/*
 * This generator is taken from
 *
 * Donald E. Knuth, The Art of Computer Programming, Volume 2, Section 3.6
 * Third Edition, Addison-Wesley,
 *
 * The modifications introduced in the 9th printing (2002) are
 * included here; there's no backwards compatibility with the
 * original.  [ see http://www-cs-faculty.stanford.edu/~knuth/taocp.html ]
 *
 */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_rng.h>

#define BUFLEN 1009             /* length of the buffer aa[] */
#define KK 100                  /* the long lag */
#define LL 37                   /* the short lag */
#define MM (1L << 30)           /* the modulus */
#define TT 70                   /* guaranteed separation between streams */

#define is_odd(x) ((x) & 1)                       /* the units bit of x */
#define mod_diff(x, y) (((x) - (y)) & (MM - 1))   /* (x - y) mod MM */

static inline void ran_array (long int aa[], unsigned int n,
                              long int ran_x[]);
static inline unsigned long int ran_get (void *vstate);
static double ran_get_double (void *vstate);
static void ran_set (void *state, unsigned long int s);

typedef struct
{
  unsigned int i;
  long int aa[BUFLEN];
  long int ran_x[KK];  /* the generator state */
}
ran_state_t;

static inline void
ran_array (long int aa[], unsigned int n, long int ran_x[])
{
  unsigned int i;
  unsigned int j;

  for (j = 0; j < KK; j++)
    aa[j] = ran_x[j];

  for (; j < n; j++)
    aa[j] = mod_diff (aa[j - KK], aa[j - LL]);

  for (i = 0; i < LL; i++, j++)
    ran_x[i] = mod_diff (aa[j - KK], aa[j - LL]);

  for (; i < KK; i++, j++)
    ran_x[i] = mod_diff (aa[j - KK], ran_x[i - LL]);
}

static inline unsigned long int
ran_get (void *vstate)
{
  ran_state_t *state = (ran_state_t *) vstate;

  unsigned int i = state->i;
  unsigned long int v;

  if (i  == 0)
    {
      /* fill buffer with new random numbers */
      ran_array (state->aa, BUFLEN, state->ran_x);
    }

  v = state->aa[i];

  state->i = (i + 1) % KK;

  return v;
}

static double
ran_get_double (void *vstate)
{
  ran_state_t *state = (ran_state_t *) vstate;

  return ran_get (state) / 1073741824.0;        /* RAND_MAX + 1 */
}

static void
ran_set (void *vstate, unsigned long int s)
{
  ran_state_t *state = (ran_state_t *) vstate;

  long x[KK + KK - 1];          /* the preparation buffer */

  register int j;
  register int t;
  register long ss;

  if (s == 0 )
    s = 314159;                 /* default seed used by Knuth */

  ss = (s + 2)&(MM-2);

  for (j = 0; j < KK; j++)
    {
      x[j] = ss;                /* bootstrap the buffer */
      ss <<= 1;
      if (ss >= MM)             /* cyclic shift 29 bits */
        ss -= MM - 2;
    }
  x[1]++;                       /* make x[1] (and only x[1]) odd */

  ss = s & (MM - 1);
  t = TT - 1;
  while (t)
    {
      for (j = KK - 1; j > 0; j--)      /* square */
        {
          x[j + j] = x[j];
          x[j + j - 1] = 0;
        }

      for (j = KK + KK - 2; j >= KK; j--)
        {
          x[j - (KK - LL)] = mod_diff (x[j - (KK - LL)], x[j]);
          x[j - KK] = mod_diff (x[j - KK], x[j]);
        }

      if (is_odd (ss))
        {                       /* multiply by "z" */
          for (j = KK; j > 0; j--)
            {
              x[j] = x[j - 1];
            }
          x[0] = x[KK];         /* shift the buffer cyclically */
          x[LL] = mod_diff (x[LL], x[KK]);
        }

      if (ss)
        ss >>= 1;
      else
        t--;
    }

  for (j = 0; j < LL; j++)
    state->ran_x[j + KK - LL] = x[j];
  for (; j < KK; j++)
    state->ran_x[j - LL] = x[j];

  for (j = 0; j< 10; j++)
    ran_array(x, KK+KK-1, state->ran_x);  /* warm things up */

  state->i = 0;

  return;
}

static const gsl_rng_type ran_type = {
  "knuthran2002",                /* name */
  0x3fffffffUL,                 /* RAND_MAX = (2 ^ 30) - 1 */
  0,                            /* RAND_MIN */
  sizeof (ran_state_t),
  &ran_set,
  &ran_get,
  &ran_get_double
};

const gsl_rng_type *gsl_rng_knuthran2002 = &ran_type;

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified