text:
/*
	STOCHASTIC LOTKA-VOLTERRA MODEL WITH MULTIDIMENSIONAL NICHE SPACES
	by Tommaso Biancalani <tommasob@mit.edu>

	Supporting Material to the paper:
	"Framework for analyzing ecological trait-based models in multidimensional niche spaces" - Phys. Rev. E 91, 052107

	ABSTRACT
	We develop a theoretical framework for analyzing ecological models with a multidimensional niche space. Our approach relies on the fact that ecological niches are described by sequences of symbols, which allows us to include multiple phenotypic traits. Ecological drivers, such as competitive exclusion, are modeled by introducing the Hamming distance between two sequences. We show that a suitable transform diagonalizes the community interaction matrix of these models, making it possible to predict the conditions for niche differentiation and, close to the instability onset, the asymptotically long time population distributions of niches. We exemplify our method using the Lotka-Volterra equations with an exponential competition kernel.

	INSTALL
 	The code requires the GSL (GNU Scientific Library). On linux, I simply compile it with:
 	$ gcc lotka-volterra.c -lgsl -lgslcblas -lm -Wall -o "lotka-volterra"
 */

#include <stdio.h>
#include <string.h>
#include <unistd.h>
#include <stdlib.h>
#include <sys/types.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h>
#include <sys/time.h>
#include <time.h>
#include <sys/select.h>
#include <math.h>

// Headers for special functions
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf_exp.h>

// Headers for Mersenne-Twister random generator
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

// Global variables
#define R_CONST		1.			// Competition kernel length scale
#define	SIGMA		2			// Competition kernel exponent
#define	V_CONST		100.		// System size of single compartment
#define	T_SAMP		1. 			// Sampling time (outputs every T_SAMP)
#define	TOT_TIME	100. 		// Evolve the system until t/V = TOT_TIME

// Sequence space definitions (sequences are called 'genomes')
typedef char GENOME_S;
typedef unsigned int GENOME_N;
#define	C_CONST	5.52486			// Competition kernel normalization constant
#define L 	3					// Phenotype length
#define NUM_GENOMES	54			// Total number of genome in the systems
const int DELTA[L] = {3,9,2};  	// It specifies the sequence space
GENOME_S *genome_0 = "000";
GENOME_S genome_db[NUM_GENOMES][L+1];
#define A(k1,k2) (k1)+2*(k2)  	// reaction k1 in genome k2
FILE *out;
unsigned int *org;				// org[I]: number of organisms with genome I

// GSL rn definitions (Mersenne-Twister)
const gsl_rng_type * rgT;
gsl_rng * rg;
double rn;
unsigned long int seed;
/**********************************************************/

GENOME_S *genome_num2str(GENOME_N I)
{
/*
	Helper for converting genomes from numbers (GENOME_N) to string (GENOME_S).
	Return converted genome in string format.
*/
	int l;
	GENOME_S *tbl = "0123456789abcdefghijklmnopqrstuvwxyz";
	GENOME_S *out;
	out = (GENOME_S *) malloc(sizeof(GENOME_S)*L);
	memcpy(out, genome_0, sizeof(GENOME_S)*L);
	int len = L-1;
	for(l=0; 1; l++)
	{
		 out[len--] = tbl[I % DELTA[l]];
		 if((I /= DELTA[l]) == 0)
		 	break;
	}
	return out;
}

unsigned int Hamming_dist(GENOME_S *genome_I, GENOME_S *genome_J)
{
/*
	Return Hamming distance between genomes I and J
*/
	int l;
	unsigned int ham_dist=L;
	for(l=0; l<L; l++)
		if(genome_I[l] - genome_J[l] == 0)
			ham_dist--;
	return ham_dist;
}

double G_IJ(GENOME_S *genome_I, GENOME_S *genome_J) {
/*
	Return competition rate G on genome I due to genome J.
*/
	double arg, ham_dist;
	ham_dist = ((double)Hamming_dist(genome_I, genome_J)) / R_CONST;
	arg = gsl_pow_int(ham_dist, SIGMA);
	return gsl_sf_exp(-arg);
}

int gillespie_evolve()
{
/*
	Print model dynamic on standard output. Dynamic is computed using Gillespie's algorithm. Return zero when time reaches TOT_TIME.
*/
	GENOME_N I, J;		// Genome indexes
	int k;				// Reaction index
	long int step; 		// Time iteration index
	int sample_step;	// Sampled time iteration index
	double tau, t; 		// Stochastic times
	double *T_rates; 	// Transition rates
	double T_sum, Tr2, T_tmp;
	double rn; 			// rand number

	// Allocate memory for transition rates
	T_rates = (double *)malloc(sizeof(double) * NUM_GENOMES * 2);

	// Main Loop
	t=0.; step=0; sample_step=0;
	while(t < TOT_TIME)
	{
		// Output as time reaches next sampled time
		if(t > (double) (T_SAMP * sample_step))
		{
			for(I=0; I<NUM_GENOMES; I++)
				fprintf(out, "%.3f ", org[I]/V_CONST);
			fprintf(out, "\n");
			sample_step++;
		}

		// Compute transiton rates
		for (I=0; I<NUM_GENOMES; I++)
		{
			T_rates[A(0, I)] = org[I]; 	// reaction 0 = birth
			T_rates[A(1, I)] = 0.; 		// reaction 1 = death
			for(J=0; J<NUM_GENOMES; J++)
				T_rates[A(1,I)] += G_IJ(genome_db[I], genome_db[J]) * org[J];
			T_rates[A(1,I)] *= org[I] / (V_CONST * C_CONST);
		}

		// Compute total transition rate
		T_sum = 0.;
		for(I=0; I<NUM_GENOMES; I++)
			T_sum += T_rates[A(0, I)] + T_rates[A(1, I)];

		// Generate next reaction time (tau)
		rn = gsl_rng_uniform_pos(rg);
		tau = (1./T_sum) * log(1./rn);

		// Update time
		t += tau;

		// Generate what reaction occurs and on which species
		Tr2 = gsl_rng_uniform_pos(rg) * T_sum;
		T_tmp = 0.;
		for (I=0; I<NUM_GENOMES; I++)
		{
			for(k=0; k<2; k++)
			{
				T_tmp += T_rates[A(k,I)];
				if(T_tmp >= Tr2)  goto gotten; // this is a fine use of goto :-)
			}
		}
		gotten:

		// Execute reaction i on species k
		switch(k) // Wich reaction?
		{
			// BIRTH
			case 0:
				org[I]++;
			break;

			// DEATH
			case 1:
				org[I]--;
			break;
		}
		// update step
		step = step+1;
	} // main loop
	return 0;
} // gillespe_evolve()

int main(int argc, char * argv[])
{
/*
	Main function. Arguments are not processed. The model parameters are set by the constants and globals defined at the beginning of the code.
*/
	GENOME_N I; // Genome index

	// Initialize random generator
	srand(time(NULL));
	seed = (unsigned long int)rand();
	rgT = gsl_rng_default; // default = Mersenne Twister
	rg = gsl_rng_alloc (rgT);
	gsl_rng_set (rg, seed);

	// Initialize organisms, e.g. org[23] = number of organisms with genome 23
	org = (unsigned int *)malloc(sizeof(unsigned int) * NUM_GENOMES);

	// Initialize system around homogeneous state
	for (I=0; I<NUM_GENOMES; I++)
		org[I] = V_CONST;

	// Set file descriptor for output
	out = stdout;

	// Initialize genome space
	for(I=0; I<NUM_GENOMES; I++)
		strcpy(genome_db[I], genome_num2str(I));

	// Compute model dynamic
	gillespie_evolve();

	// Free descriptors and quit
	fclose(out);
	gsl_rng_free(rg);
	return 0;
} // main()

max_stars_repo_name: Llewlyn/Stochastic-multidimensional-LV

max_issues_repo_name: Llewlyn/Stochastic-multidimensional-LV

max_forks_repo_name: Llewlyn/Stochastic-multidimensional-LV


text:
//
// ePDF
//

#pragma once

#include "ePDF/ndistributions.h"

#include <yaml-cpp/yaml.h>
#include <complex>
#include <vector>
#include <memory>
#include <gsl/gsl_integration.h>
#include <functional>

namespace ePDF
{
  /**
   * @brief The "evol_params" structure contains the evolution
   * parameters.
   */
  struct evol_params
  {
    double                          x;
    double                          Q;
    int                             id;
    std::shared_ptr<NDistributions> Ndist;
    evol_params operator = (evol_params const& p)
    {
      x     = p.x;
      Q     = p.Q;
      id    = p.id;
      Ndist = p.Ndist;
      return *this;
    }
  };

  /**
   * @brief The "xDistribution" class
   */
  class xDistributions
  {
  public:
    /**
     * @brief The "xDistribution" constructor.
     * @param config: the YAML:Node with the parameters
     */
    xDistributions(YAML::Node const& config);

    /**
     * @brief The "xDistribution" destructor.
     */
    ~xDistributions();

    /**
     * @brief Function that returns the PDFs in x space.
     * @param x: Bjorken x
     * @param Q: the final scale
     */
    std::vector<double> Evolve(double const& x, double const& Q);

    /**
     * @brief Interfaces to the integrand functions according to the
     * path chosen to perform the inverse Mellin transform.
     */
    std::vector<double> integrand(double const& y) const;
    std::vector<double> talbot(double const& y) const;
    std::vector<double> straight(double const& y) const;

    /**
     * @brief Integrators functions
     */
    std::vector<double> trapezoid(double const& a, double const& b) const;
    std::vector<double> gauss(double const& a, double const& b) const;
    std::vector<double> gaussGSL(double const& a, double const& b);

    /**
     * @brief Function that returns the N-th (complex) moment of PDFs.
     * @param N: moment
     * @param Q: the final scale
     */
    std::vector<std::complex<double>> Moments(std::complex<double> const& N, double const& Q) const;

    /**
     * @brief Function that sets the parameter structure externally
     * @param p: the input structure
     */
    void SetParameters(evol_params const& p) { _p = p; };

  private:
    std::string               const _contour;
    std::string               const _integrator;
    double                    const _eps;
    evol_params                     _p;
    gsl_integration_workspace      *_w;
  };
}

ext: h

max_stars_repo_name: vbertone/ePDF

max_issues_repo_name: gstagnit/ePDF

max_forks_repo_name: gstagnit/ePDF


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

/* zsolve.c - finds the complex roots of  = 0 */

#include <config.h>
#include <math.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_poly.h>

/* C-style matrix elements */
#define MAT(m,i,j,n) ((m)[(i)*(n) + (j)])

/* Fortran-style matrix elements */
#define FMAT(m,i,j,n) ((m)[((i)-1)*(n) + ((j)-1)])

#include "companion.c"
#include "balance.c"
#include "qr.c"

int
gsl_poly_complex_solve (const double *a, size_t n,
                        gsl_poly_complex_workspace * w,
                        gsl_complex_packed_ptr z)
{
  int status;
  double *m;

  if (n == 0)
    {
      GSL_ERROR ("number of terms must be a positive integer", GSL_EINVAL);
    }

  if (n == 1)
    {
      GSL_ERROR ("cannot solve for only one term", GSL_EINVAL);
    }

  if (a[n - 1] == 0)
    {
      GSL_ERROR ("leading term of polynomial must be non-zero", GSL_EINVAL) ;
    }

  if (w->nc != n - 1)
    {
      GSL_ERROR ("size of workspace does not match polynomial", GSL_EINVAL);
    }

  m = w->matrix;

  set_companion_matrix (a, n - 1, m);

  balance_companion_matrix (m, n - 1);

  status = qr_companion (m, n - 1, z);

  if (status)
    {
      GSL_ERROR("root solving qr method failed to converge", GSL_EFAILED);
    }

  return GSL_SUCCESS;
}

ext: c

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#define _GNU_SOURCE
#include <stdio.h>
#include <stdarg.h>
#include <math.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv2.h>

/* don't check ranges in Matrix operations and make Matrix call inline */
#define GSL_RANGE_CHECK_OFF
//#define HAVE_INLINE

#define FALSE 0
#define TRUE 1

#ifndef M_PI
#define M_PI            3.141592653589793238462643383280        /* pi */
#endif

#define LN2PI 1.837877066409345339

#define ODEdriver gsl_odeiv2_driver_alloc_y_new
#define rkf45 gsl_odeiv2_step_rkf45
#define rk8pd gsl_odeiv2_step_rk8pd
#define rkck  gsl_odeiv2_step_rkck
#define rk2   gsl_odeiv2_step_rk2

#define EOL -1

#define max(A, B) ((A) > (B) ? (A) : (B))
#define min(A, B) ((A) < (B) ? (A) : (B))
#define MIN(A, B) ((A) < (B) ? (A) : (B))
#define MAX(A, B) ((A) > (B) ? (A) : (B))
#define V(x) (x)*(x)

#define FIT_TIME_POINT(x) (x->mode == FIT)
#define RECORD_TIME_POINT(x, u) (x->mode == OUTPUT || (x->mode == FIT && u == x->t[x->i]))
#define NEXT_TIME_POINT(x) if (x->mode == FIT) x->i++
#define RESET_TIME_POINT(x) x->i = 0

enum {FIT, OUTPUT};

enum {SUCCESS=0, FAIL=1, ERROR1=1, ERROR2, ERROR3, ERROR4, ERROR5, ERROR6,
      ERROR7, ERROR8, ERROR9, ERROR10};

typedef uint boolean;

// for internal use
// type of parameter
enum {NOTUSED=0, DERIVED=21, FITTED=22, CONSTANT=23, AUXILLARY=24, INDIVIDUAL_CONSTANT=28};
// whether a fitted or auxillary parameter has a hyperparameter associated with it
enum {HYPER=26};
// prior of a fitted parameter that does not have a hyperparameters
enum {NOTAPPLICABLE=0, NORMAL=1, GAMMA=2, SICHI2=3, UNIFORM=4, PARETO=5,
      TRUNCNORMAL=6, EXPONENTIAL=7, BINOMIAL=8, POISSON=9, INVGAMMA=10,
      UNIFORM0=11, BETA=12, NORMALMEAN=13, CATEGORICAL=14, GEOMETRIC=15,
      USERDEFINED=16, AUXILLARY_HYPER=27, CIRCULAR=29};
// number type of the parameter
enum {REAL=17, INTEGER=18, CATEGORY=19, INT = INTEGER};
// type of data
enum {NOT_AVAILABLE=0, ABOVE_DETECTION=1, BELOW_DETECTION=2, MEASURED=3};

// for use in model files
// prior or type of the declared parameter
enum {Normal        =NORMAL,
      Gamma         =GAMMA,
      Sichi2        =SICHI2,
      Uniform       =UNIFORM,
      Circular      =CIRCULAR,
      Pareto        =PARETO,
      TruncNormal   =TRUNCNORMAL,
      Exponential   =EXPONENTIAL,
      Binomial      =BINOMIAL,
      Poisson       =POISSON,
      InvGamma      =INVGAMMA,
      Uniform0      =UNIFORM0,
      Beta          =BETA,
      NormalMean    =NORMALMEAN,
      Categorical   =CATEGORICAL,
      Geometric     =GEOMETRIC,
      UserDefined   =USERDEFINED,
      HyperP         =HYPER,
      Constant      =CONSTANT,
      Auxillary     =AUXILLARY,
      Auxillary_Hyper =AUXILLARY_HYPER,
      Derived       =DERIVED,
      Individual_Constant=INDIVIDUAL_CONSTANT};
enum {Nor           =NORMAL,
      Gam 	    =GAMMA,
      Sic 	    =SICHI2,
      Unif 	    =UNIFORM,
      Circ          =CIRCULAR,
      Pare 	    =PARETO,
      TNorm 	    =TRUNCNORMAL,
      Exp 	    =EXPONENTIAL,
      Bin 	    =BINOMIAL,
      Poi 	    =POISSON,
      IGam 	    =INVGAMMA,
      Unif0 	    =UNIFORM0,
      Be 	    =BETA,
      NorM 	    =NORMALMEAN,
      Cat 	    =CATEGORICAL,
      Geom 	    =GEOMETRIC,
      UD 	    =USERDEFINED,
      Hyp 	    =HYPER,
      Cst 	    =CONSTANT,
      Aux 	    =AUXILLARY,
      Aux_Hyp       =AUXILLARY_HYPER,
      Der	    =DERIVED,
      IndCst        =INDIVIDUAL_CONSTANT};

// number type of the parameter
enum {Real=REAL, Integer=INTEGER, Category=CATEGORY, Int=INTEGER};
// type of data
enum {NA=0, Above_Detection=1, Below_Detection=2, Measured=3};

typedef struct {
  uint mode, n, i;
  double *t;
} TimePoints;

struct line {
  int n;
  double *x;
  double *y;
  double *m;
  double *c;
};

struct data_individual {
  int n;
  int np;
  uint **value;
  double **Y;
  double saturated;
};

struct data_treatment {
  int nind;
  int *usevar;
  char levels[1000];
  struct data_individual *Individual;
};

struct data {
  int ntrt;
  int nvar;
  int maxnind;
  int maxnp;
  char filename[100];
  struct data_treatment *Treatment;
};

typedef struct treatmentlist {
  int n;
  int *treatment;
} Treatments;

typedef struct hyperblock {
  /* set in CreateBlock */
  int use;

  /* set in AddHyperToTreatment */
  int nhypers_used;
  int *hyper_used;
  int **nhypers_with_treatment;
  int ***hyper_with_treatment;
  int nassociated_pars;
  int *associated_par;
  int *ntreatments_with_parameter;
  int **treatment_with_parameter;

  /* set in Bayes */
  int hdistribution;
  int pdistribution;
  int *dependent_block;
} Hyperblock;

struct hyper {
  int index;
  int use;
  int parno;
  int hdist;
  int number;
  char name[100];
  double location;
  double scale;
  double shape;
};

typedef struct allhypers {
  int init;
  int ntrt;
  int *nind;
  int *trtused;
  int nhyper_parameters;
  struct hyper *H;
  int nblocks;
  Hyperblock *block;
  double *hyper;         //used for initialisation of thread hyper values
  double **hyper_store;
} Hyperparameters;

struct variable {
  int index;
  char name[100];
};

struct variables {
  int    variables;
  double scaletime;
  struct variable *W;
};

struct parameter {
  int index;
  int number;
  int distribution;
  int circular;
  int type;
  int hashyper;
  int hyper_location_block_index;
  int hyper_shape_block_index;
  int hyper_scale_block_index;
  int ncategories;
  double *p; // used for categorical parameters
  double location;
  double scale;
  double shape;
  char name[100];
  double (*(*link))(double);
};

struct blocks {
  int nblocks;
  int **pblock;
};

struct parameterset {
  int    pass;
  int    parameters;   // maximum parameter index
  int    npar;         // number of fitted parameters
  int    ninteger;
  int    *integer;
  int    nauxillary;
  int    *auxillary;
  struct parameter *P;
  struct blocks Blocks;
};

struct parameter_block {
  int categorical;
  int done;
  int npar;
  int *pblock;
  int covfinal;
  int covstart;
  int inrange;
  int storeno;
  int error;
  double accept;      /*mcmc acceptances*/
  double proposed;    /*mcmc iterations*/
  double covscale;    /*covariance matrix scaling factor*/
  //  double **L;         /*covariance matrix*/
  gsl_matrix *L;
  double **theta_store;
};

struct chain {
  int done;       /*if current individual is ready for next phase*/
  int loop;
  /* parameters for threads */
  double logmaxL;
  double logL;
  double *logP;
  double logpost;
  double *theta;       /*current parameter values*/
  double minlogL;
  double prop_logL;
  double prop_logP;
  double prop_logpost;
  int storeno;
  int laststoreno;
  double accept;      /*mcmc acceptances*/
  double proposed;    /*mcmc iterations*/
  double **theta_store;
  double *logLstore;
  double *logpoststore;
  double exchangeaccept;      /*mcmc acceptances*/
  double exchangeproposed;    /*mcmc iterations*/
  struct parameter_block *PB;
};

struct sim_individual {
  int done;
  int chains;
  int *chainno;
  double *temp;
  struct chain *Chain;
  struct parameterset Q;
};

typedef struct {
  int n;
  int m;
  int **a;
} Matrix;

struct sim_treatment {
  struct sim_individual *Individual;
};

struct thread {
  double *hyper;
  struct sim_treatment *Treatment;
  gsl_rng *stream;
};

typedef struct simulation {
  Hyperparameters HP;
  int *trtused;
  int **indused;
  int ntp;            // individual with most parameters
  int nrp;            // individual with most parameters
  int rdp;
  int sim;
  int maxc;
  int seed;
  int thin;
  int mcmc;
  int popd;
  int marg;
  int maxnp;
  int error;
  int range;
  int chains;
  int notime;
  int copyind;
  int threads;
  int verbose; /* 0: no output, 1: output chain 0+lowest exchangerate,
		  2: all chains*/
  int derived;
  int discard;
  int samples;
  int storeno;
  int gsl_error;
  int population;
  int covsamples;
  int growthcurves;
  int parsperblock;
  int phase1adjust;
  int phase2adjust;
  int mincovsamples;
  int ignorewarning;
  int usehypermeans;
  int sample_parameters;
  char expno[50];
  char outexpno[50];
  char *parsfile;
  char *tempfile;
  char logfile[1000];
  char factors[1000];
  double xrate;
  double timeend;
  double lowrate;
  double highrate;
  double minprior;
  double initscale;
  double covthresh;
  double timestart;
  double scalefactor;
  double saturatedlogL;
  double minexchangerate;
  struct sim_individual Sample;
  struct sim_treatment *Treatment;
  struct variables V;
} Simulation;

extern char expno[];

extern int* integervector( int );
extern uint* uintegervector( int );
extern uint** uintegermatrix( int, int );
extern double* doublevector( int );
extern double** doublematrix( int, int );
extern double*** doublematrix3( int, int, int );
extern double**** doublematrix4( int, int, int, int );
extern char* charvector( int );
extern char** charmatrix( int, int );
extern double factln( int );
extern double gammln( double );
extern void binorm( double*, double, double, double, double, double );
extern long int lrint(double);
extern double pow10(double);
extern double round(double);
extern double rtnewt(void(*)(double, double*, double*, double*),
		     double, double, double, double*);
extern double rtbis( double (*)(double, double*), double, double,
		     double, int*, double*);

extern void *malloc(size_t);
extern void exit(int);
extern void free(void*);
extern void Error(char*);
extern void ReadParam( char*, void*, char* ) ;

extern uint   binomialRNG(uint, double, const gsl_rng*);
extern void   multinomialRNG(uint, const double*, size_t, unsigned int*, const gsl_rng*);
extern int    poissonRNG(double, const gsl_rng*);
extern int    uniform_intRNG(int, int, const gsl_rng*);
extern double normalRNG(double, double, const gsl_rng*);
extern double betaRNG(double, double, const gsl_rng*);
extern double chisqRNG(double, const gsl_rng*);
extern double scaleinvchisqRNG(double, double, const gsl_rng*);
extern double invgammaRNG(double, double, const gsl_rng*);
extern double gammaRNG(double, double, const gsl_rng*);
extern double uniformRNG(double, double, const gsl_rng*);
extern double paretoRNG(double, double, const gsl_rng*);
extern void   dirichlet_multinomialRNG(uint, const double*, size_t, uint*, const gsl_rng*);
extern void   dirichletRNG(const double*, size_t, double*, const gsl_rng*);
extern uint   beta_binomialRNG(uint, double, double, const gsl_rng*);
extern uint   negative_binomialRNG(double, double, const gsl_rng*);
extern uint   beta_negative_binomialRNG(uint, double, double, const gsl_rng*);
/* Matrix *wishartRNG(int, int, Matrix*, Matrix*, const gsl_rng*); */
/* void   MVNRNG(int, Matrix*, Matrix*, Matrix*, const gsl_rng*); */

extern double poissonCMF(const double, const double, int);
extern double gamma_scalePDF(const double, const double, const double);
extern double betaPDF(const double, const double, const double);
extern double exponentialPDF(const double, const double);
extern double poissonPMF(const double, const double);
extern double geometricPDF(const double, const double);
extern double trnormalPDF(const double, const double, const double);
extern double normalCDF(const double, const double, const double, int);
extern double normalPDF(const double, const double, const double);
extern double binomialPMF(const double, const double, const double);
extern double multinomialPMF(const unsigned int, const double*, const double*);
extern double uniformPDF(const double, const double, const double);
extern double logisticPDF(const double, const double, const double);
extern double betabinomialPMF(const double, const double, const double, const double);
//extern double MVNPDF(Matrix*, Matrix*, Matrix*, const double, int);
extern double von_mises_cdf ( double x, double a, double b );
extern double von_mises_cdf_inv ( double cdf, double a, double b );

extern Matrix *MatrixInit(Matrix*, int, int, int*);
extern Matrix *MatrixDestroy(Matrix*);
extern double identity(double);
extern double logit(double);
extern double invlogit(double);

int    function(int, int, double*, double*, double*, TimePoints*);
int    Model(const uint, const uint, const uint, double*, double**,
	     TimePoints*);
//int    Model(int, int, int, double*, double**);
void   Block(struct parameterset*, int, ...);
void   Par(struct parameterset*, int, int, char*, int, ...);
void Create(Hyperparameters *H, int *h, int block, int par,
	    Matrix M, char *name, int type, double p1, double p2);
void   Hyper(Hyperparameters *HP, int index, int number, char *name, int distribution,
	     double shape_or_loc, double scale);
void   CreateHyperBlock(Hyperparameters *HP, int index);
void   AddHyperToATreatment(Hyperparameters *HP, int index, int h, int pindex,
			    int treatment);
void   AddHyperToTreatments(Hyperparameters *HP, int index, int h, int pindex, ...);
void   HyperParameters(Treatments, Hyperparameters*);
void   Var(struct variables*, int, char*);
void   GlobalParameters(void);
void   Variables(struct variables*);
void   OutputModel(int, int, double*, double*);
void   OutputData(int, int, double*, double*, double*, uint*);
void   DerivedHyperParameters(int trt, double *p);
void   DerivedParameters(int, int, int, double*, double**, TimePoints*);
void   Parameters(int, int, double*, struct parameterset*);
void   PredictData(int, int, double*, double*, double*, gsl_rng*);
void   SimulateData(int, int, double*, double*, double*, double*, uint*, gsl_rng*);
void   SaturatedModel(int, int, double*, double*, double*, uint*);
void   Residual(int, int, double*, double*, double*, double*, uint*);
void   Equations(double*, double*, int, double*, double*, int, double);
double timestep(void);
double logLikelihood(const uint, const uint, const double*,
		     const double*, const double*, const uint*);
double logLikelihoodI(const uint trt, const uint ind, const double *p,
		      double **V, double *Data, const uint *value,
		      const uint N);
void WAIC(int trt, int ind, double *lnL, double *V,
	  double *p, double *Data, uint *value);
//double logLikelihood(int, int, double*, double*, double*, int*);
double UserPDF(double);
void   ScaleTime(struct variables*, double);


text:
static char help[] =
"1D Stokes problem with DMDA and SNES.  Option prefix -wel_.\n\n";

// show solution (staggered):
// ./well -snes_monitor -snes_converged_reason -da_refine 7 -snes_monitor_solution draw -draw_pause 1

// try:
//   -ksp_type preonly -pc_type svd                         [default]
//   -ksp_type minres -pc_type none -ksp_converged_reason
//   -ksp_type gmres -pc_type none -ksp_converged_reason
// and fieldsplit below

// Jacobian is correct and symmetric:
// ./well -snes_monitor -da_refine 3 -snes_type test
// ./well -snes_monitor -da_refine 3 -mat_is_symmetric

// generate matrix in matlab
// ./well -snes_converged_reason -da_refine 1 -mat_view ascii:foo.m:ascii_matlab
// then:
// >> M = [whos to get name]
// >> A = M(1:2:end,1:2:end);  BT = M(1:2:end,2:2:end);  B = M(2:2:end,1:2:end);  C = M(2:2:end,2:2:end);  T = [A BT; B C]

/* Schur-complement preconditioning:

./well -snes_converged_reason -snes_monitor -da_refine 7 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -fieldsplit_pressure_pc_type none -ksp_converged_reason

   * see snes example ex70.c
   * note these -ksp_type also work:  gmres, cgs, richardson
   * above is same as  -fieldsplit_velocity_pc_type ilu

with minres need positive definite preconditioner so use -pc_fieldsplit_schur_fact_type diag (see http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/PC/PCFieldSplitSetSchurFactType.html):

./well -snes_converged_reason -snes_monitor -da_refine 7 -ksp_type minres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type diag -fieldsplit_pressure_pc_type none -ksp_converged_reason

quote from http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/PC/PCFieldSplitSetSchurFactType.html:

    If applied exactly, FULL factorization is a direct solver. The preconditioned
    operator with LOWER or UPPER has all eigenvalues equal to 1 and minimal polynomial
    of degree 2, so KSPGMRES converges in 2 iterations. If the iteration count is very
    low, consider using KSPFGMRES or KSPGCR which can use one less preconditioner
    application in this case. Note that the preconditioned operator may be highly
    non-normal, so such fast convergence may not be observed in practice. With DIAG,
    the preconditioned operator has three distinct nonzero eigenvalues and minimal
    polynomial of degree at most 4, so KSPGMRES converges in at most 4 iterations.

   * indeed we see FULL converge in 1 iterations ... it IS a direct solver like SVD
   * and LOWER/UPPER + GMRES give convergence in 2 iterations as advertised
   * Murphy, Golub, Wathen (2000) explains final comment re DIAG and gmres
   * in practice DIAG + GMRES

schur as direct solver:

./well -snes_converged_reason -snes_monitor -da_refine 7 -ksp_type preonly -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -fieldsplit_pressure_pc_type none -ksp_converged_reason

view all matrices (both interlaced and blocks) as dense with fieldsplit:

./well -snes_converged_reason -snes_monitor -da_refine 1 -ksp_type preonly -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -fieldsplit_pressure_pc_type none -mat_view ::ascii_dense

put diagonal blocks only in matlab files foov.m, foop.m (off-diagonal: how?)
./well -snes_converged_reason -snes_monitor -da_refine 1 -ksp_type preonly -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -fieldsplit_pressure_pc_type none -fieldsplit_velocity_mat_view :foov.m:ascii_matlab -fieldsplit_pressure_mat_view :foop.m:ascii_matlab

40 second run (linux-c-opt) on 8 million grid points:

timer ./well -snes_converged_reason -snes_monitor -da_refine 22 -ksp_type gmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type diag -fieldsplit_pressure_pc_type none -ksp_converged_reason -fieldsplit_velocity_ksp_type cg -fieldsplit_velocity_pc_type icc

   * also works with -ksp_type minres
   * also in parallel (no speedup) if -fieldsplit_velocity_pc_type bjacobi -fieldsplit_velocity_sub_pc_type icc
*/

#include <petsc.h>
#include <time.h>

typedef enum {STAGGERED, REGULAR} SchemeType;
static const char *SchemeTypes[] = {"staggered","regular",
                                    "SchemeType", "", NULL};

typedef struct {
    double     L,       // length (height) of well (m)
               rho,     // density of water (kg m-3)
               g,       // acceleration of gravity (m s-2)
               mu;      // dynamic viscosity of water (Pa s)
    SchemeType scheme;
} AppCtx;

typedef struct {
  double u, p;
} Field;

PetscErrorCode ExactSolution(DMDALocalInfo *info, Vec X, AppCtx *user) {
    PetscErrorCode ierr;
    double h, x;
    Field  *aX;
    int    i;
    h  = user->L / (info->mx-1);
    ierr = DMDAVecGetArray(info->da,X,&aX); CHKERRQ(ierr);
    for (i=info->xs; i<info->xs+info->xm; i++) {
        aX[i].u = 0.0;
        if (user->scheme == STAGGERED)
            x = h * (i + 0.5);
        else
            x = h * i;
        if (i < info->mx - 1)
            aX[i].p = user->rho * user->g * (user->L - x);
        else
            aX[i].p = 0.0;
    }
    ierr = DMDAVecRestoreArray(info->da,X,&aX); CHKERRQ(ierr);
    return 0;
}

// the staggered version of residual evaluation F(X)
// grid has p at staggered locations "O" where incompressibility u_x=0 is enforced
//    x_i-1        O        x_i        O        x_i+1        O
//    u_i-1                 u_i                 u_i+1
//               p_i-1                p_i                  p_i+1
PetscErrorCode FormFunctionStaggeredLocal(DMDALocalInfo *info, Field *X,
                                          Field *F, AppCtx *user) {
    const double h  = user->L / (info->mx-1),
                 h2 = h * h;
    int          i;
    for (i=info->xs; i<info->xs+info->xm; i++) {
        if (i == 0) { // bottom of well
            F[i].u = X[i].u;                                              // u(0) = 0
            F[i].p = - (X[i+1].u - 0.0) / h;                              // -u_x(0+1/2) = 0
        } else if (i == 1) {
            F[i].u = - user->mu * (X[i+1].u - 2 * X[i].u + 0.0) / h2      // -mu u_xx(x1) + p_x(x1) = - rho g
                       + (X[i].p - X[i-1].p) / h
                       + user->rho * user->g;
            F[i].p = - (X[i+1].u - X[i].u) / h;                           // - u_x(x1+1/2) = 0
        // generic cases:  - mu u_xx(xi) + p_x(xi) = - rho g  AND  - u_x(xi+1/2) = 0
        } else if (i > 1 && i < info->mx - 1) {
            F[i].u = - user->mu * (X[i+1].u - 2 * X[i].u + X[i-1].u) / h2
                       + (X[i].p - X[i-1].p) / h
                       + user->rho * user->g;
            F[i].p = - (X[i+1].u - X[i].u) / h;
        } else if (i == info->mx - 1) { // top of well
            F[i].u = - user->mu * (- 2 * X[i].u + 2 * X[i-1].u) / h2      // -mu u_xx(xm-1) + p_x(xm-1) = - rho g
                       + (- 2 * X[i-1].p) / h + user->rho * user->g;      // and  u_x(xm-1) = 0  and  p(xm-1) = 0
            F[i].u /= 2;                                                  // for symmetry
            F[i].p = X[i].p;                                              // no actual d.o.f. here
        } else {
            SETERRQ(PETSC_COMM_WORLD,1,"no way to get here");
        }
    }
    return 0;
}

PetscErrorCode FormFunctionRegularLocal(DMDALocalInfo *info, Field *X,
                                        Field *F, AppCtx *user) {
    const double h  = user->L / (info->mx-1),
                 h2 = h * h;
    int          i;
    for (i=info->xs; i<info->xs+info->xm; i++) {
        if (i == 0) { // bottom of well
            F[i].u = X[i].u;                                              // u(0) = 0
            F[i].p = - (X[i+1].u - 0.0) / (2 * h);                        // -u_x(0+1/2) = 0  (and / 2 for symmetry)
        } else if (i == 1) {
            F[i].u = - user->mu * (X[i+1].u - 2 * X[i].u + 0.0) / h2      // -mu u_xx(x1) + p_x(x1) = - rho g
                       + (X[i+1].p - X[i-1].p) / (2 * h)
                       + user->rho * user->g;
            F[i].p = - (X[i+1].u - 0.0) / (2 * h);                        // - u_x(x1) = 0
        // generic cases:  - mu u_xx(xi) + p_x(xi) = - rho g  AND  - u_x(xi) = 0
//STARTREGFUNC
        } else if (i > 1 && i < info->mx - 2) {
            F[i].u = - user->mu * (X[i+1].u - 2 * X[i].u + X[i-1].u) / h2
                       + (X[i+1].p - X[i-1].p) / (2 * h)
                       + user->rho * user->g;
            F[i].p = - (X[i+1].u - X[i-1].u) / (2 * h);
//ENDREGFUNC
        } else if (i == info->mx - 2) {
            F[i].u = - user->mu * (X[i+1].u - 2 * X[i].u + X[i-1].u) / h2 // -mu u_xx(xm-2) + p_x(xm-2) = - rho g
                       + (0.0 - X[i-1].p) / (2 * h)
                       + user->rho * user->g;
            F[i].p = - (X[i+1].u - X[i-1].u) / (2 * h);                   // - u_x(xm-2) = 0
        } else if (i == info->mx - 1) { // top of well
            F[i].u = - user->mu * (- 2 * X[i].u + 2 * X[i-1].u) / h2      // -mu u_xx(xm-1) + p_x(xm-1) = - rho g
                       + (- X[i-1].p) / h + user->rho * user->g;          // and  u_x(xm-1) = 0  and  p(xm-1) = 0
            F[i].u /= 2;                                                  // for symmetry
            F[i].p = X[i].p;                                              // p(xm-1) = 0
        } else {
            SETERRQ(PETSC_COMM_WORLD,1,"no way to get here");
        }
    }
    return 0;
}

PetscErrorCode FormJacobianStaggeredLocal(DMDALocalInfo *info, double *X,
                                          Mat J, Mat P, AppCtx *user) {
    PetscErrorCode ierr;
    MatStencil   col[5],row;
    double       v[5];
    const double h  = user->L / (info->mx-1),
                 h2 = h * h;
    int          i;
    for (i=info->xs; i<info->xs+info->xm; i++) {
        row.i = i;
        if (i == 0) {
            row.c = 0;  col[0].i = i;    col[0].c = 0;  v[0] = 1.0;
            ierr = MatSetValuesStencil(P,1,&row,1,col,v,INSERT_VALUES); CHKERRQ(ierr);
            row.c = 1;  col[0].i = i+1;  col[0].c = 0;  v[0] = - 1.0 / h;
            ierr = MatSetValuesStencil(P,1,&row,1,col,v,INSERT_VALUES); CHKERRQ(ierr);
        } else if (i == 1) {
            row.c = 0;
            col[0].c = 0;  col[0].i = i;    v[0] = 2.0 * user->mu / h2;
            col[1].c = 0;  col[1].i = i+1;  v[1] = - user->mu / h2;
            col[2].c = 1;  col[2].i = i;    v[2] = 1.0 / h;
            col[3].c = 1;  col[3].i = i-1;  v[3] = - 1.0 / h;
            ierr = MatSetValuesStencil(P,1,&row,4,col,v,INSERT_VALUES); CHKERRQ(ierr);
            row.c = 1;
            col[0].c = 0;  col[0].i = i;    v[0] = 1.0 / h;
            col[1].c = 0;  col[1].i = i+1;  v[1] = - 1.0 / h;
            ierr = MatSetValuesStencil(P,1,&row,2,col,v,INSERT_VALUES); CHKERRQ(ierr);
        } else if (i > 1 && i < info->mx - 1) {
            row.c = 0;
            col[0].c = 0;  col[0].i = i;    v[0] = 2.0 * user->mu / h2;
            col[1].c = 0;  col[1].i = i-1;  v[1] = - user->mu / h2;
            col[2].c = 0;  col[2].i = i+1;  v[2] = - user->mu / h2;
            col[3].c = 1;  col[3].i = i;    v[3] = 1.0 / h;
            col[4].c = 1;  col[4].i = i-1;  v[4] = - 1.0 / h;
            ierr = MatSetValuesStencil(P,1,&row,5,col,v,INSERT_VALUES); CHKERRQ(ierr);
            row.c = 1;
            col[0].c = 0;  col[0].i = i;    v[0] = 1.0 / h;
            col[1].c = 0;  col[1].i = i+1;  v[1] = - 1.0 / h;
            ierr = MatSetValuesStencil(P,1,&row,2,col,v,INSERT_VALUES); CHKERRQ(ierr);
        } else if (i == info->mx - 1) {
            row.c = 0;
            col[0].c = 0;  col[0].i = i;    v[0] = user->mu / h2;
            col[1].c = 0;  col[1].i = i-1;  v[1] = - user->mu / h2;
            col[2].c = 1;  col[2].i = i-1;  v[2] = - 1.0 / h;
            ierr = MatSetValuesStencil(P,1,&row,3,col,v,INSERT_VALUES); CHKERRQ(ierr);
            row.c = 1;
            col[0].c = 1;  col[0].i = i;    v[0] = 1.0;
            ierr = MatSetValuesStencil(P,1,&row,1,col,v,INSERT_VALUES); CHKERRQ(ierr);
        } else {
            SETERRQ(PETSC_COMM_WORLD,1,"no way to get here");
        }
    }

    ierr = MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
    ierr = MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
    if (J != P) {
        ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
        ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    }
    return 0;
}

PetscErrorCode FormJacobianRegularLocal(DMDALocalInfo *info, double *X,
                                        Mat J, Mat P, AppCtx *user) {
    PetscErrorCode ierr;
    MatStencil   col[5],row;
    double       v[5];
    const double h  = user->L / (info->mx-1),
                 h2 = h * h;
    int          i;
    for (i=info->xs; i<info->xs+info->xm; i++) {
        row.i = i;
        if (i == 0) {
            row.c = 0;  col[0].i = i;    col[0].c = 0;  v[0] = 1.0;
            ierr = MatSetValuesStencil(P,1,&row,1,col,v,INSERT_VALUES); CHKERRQ(ierr);
            row.c = 1;  col[0].i = i+1;  col[0].c = 0;  v[0] = - 1.0 / (2.0 * h);
            ierr = MatSetValuesStencil(P,1,&row,1,col,v,INSERT_VALUES); CHKERRQ(ierr);
        } else if (i == 1) {
            row.c = 0;
            col[0].c = 0;  col[0].i = i;    v[0] = 2.0 * user->mu / h2;
            col[1].c = 0;  col[1].i = i+1;  v[1] = - user->mu / h2;
            col[2].c = 1;  col[2].i = i+1;  v[2] = 1.0 / (2.0 * h);
            col[3].c = 1;  col[3].i = i-1;  v[3] = - 1.0 / (2.0 * h);
            ierr = MatSetValuesStencil(P,1,&row,4,col,v,INSERT_VALUES); CHKERRQ(ierr);
            row.c = 1;
            col[0].c = 0;  col[0].i = i+1;    v[0] = - 1.0 / (2.0 * h);
            ierr = MatSetValuesStencil(P,1,&row,1,col,v,INSERT_VALUES); CHKERRQ(ierr);
//STARTREGMAT
        } else if (i > 1 && i < info->mx - 2) {
            row.c = 0;
            col[0].c = 0;  col[0].i = i;    v[0] = 2.0 * user->mu / h2;
            col[1].c = 0;  col[1].i = i-1;  v[1] = - user->mu / h2;
            col[2].c = 0;  col[2].i = i+1;  v[2] = - user->mu / h2;
            col[3].c = 1;  col[3].i = i+1;  v[3] = 1.0 / (2.0 * h);
            col[4].c = 1;  col[4].i = i-1;  v[4] = - 1.0 / (2.0 * h);
            ierr = MatSetValuesStencil(P,1,&row,5,col,v,INSERT_VALUES); CHKERRQ(ierr);
            row.c = 1;
            col[0].c = 0;  col[0].i = i-1;  v[0] = 1.0 / (2.0 * h);
            col[1].c = 0;  col[1].i = i+1;  v[1] = - 1.0 / (2.0 * h);
            ierr = MatSetValuesStencil(P,1,&row,2,col,v,INSERT_VALUES); CHKERRQ(ierr);
//ENDREGMAT
        } else if (i == info->mx - 2) {
            row.c = 0;
            col[0].c = 0;  col[0].i = i;    v[0] = 2.0 * user->mu / h2;
            col[1].c = 0;  col[1].i = i-1;  v[1] = - user->mu / h2;
            col[2].c = 0;  col[2].i = i+1;  v[2] = - user->mu / h2;
            col[3].c = 1;  col[3].i = i-1;  v[3] = - 1.0 / (2.0 * h);
            ierr = MatSetValuesStencil(P,1,&row,4,col,v,INSERT_VALUES); CHKERRQ(ierr);
            row.c = 1;
            col[0].c = 0;  col[0].i = i-1;  v[0] = 1.0 / (2.0 * h);
            col[1].c = 0;  col[1].i = i+1;  v[1] = - 1.0 / (2.0 * h);
            ierr = MatSetValuesStencil(P,1,&row,2,col,v,INSERT_VALUES); CHKERRQ(ierr);
        } else if (i == info->mx - 1) {
            row.c = 0;
            col[0].c = 0;  col[0].i = i;    v[0] = user->mu / h2;
            col[1].c = 0;  col[1].i = i-1;  v[1] = - user->mu / h2;
            col[2].c = 1;  col[2].i = i-1;  v[2] = - 1.0 / (2.0 * h);
            ierr = MatSetValuesStencil(P,1,&row,3,col,v,INSERT_VALUES); CHKERRQ(ierr);
            row.c = 1;
            col[0].c = 1;  col[0].i = i;    v[0] = 1.0;
            ierr = MatSetValuesStencil(P,1,&row,1,col,v,INSERT_VALUES); CHKERRQ(ierr);
        } else {
            SETERRQ(PETSC_COMM_WORLD,1,"no way to get here");
        }
    }

    ierr = MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
    ierr = MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
    if (J != P) {
        ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
        ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    }
    return 0;
}

int main(int argc,char **args) {
    PetscErrorCode ierr;
    DM            da;
    SNES          snes;
    KSP           ksp;
    PC            pc;
    AppCtx        user;
    Vec           X, Xexact;
    PetscBool     randominit = PETSC_FALSE, shorterrors = PETSC_FALSE;
    double        uerrnorm, perrnorm, pnorm, setol = 1.0e-8;
    DMDALocalInfo info;

    PetscInitialize(&argc,&args,NULL,help);
    user.rho = 1000.0;
    user.g   = 9.81;
    user.L   = 10.0;
    user.mu  = 1.0;   // Pa s; = 1.0 for corn syrup; = 10^-3 for liquid water
    user.scheme = STAGGERED;

    ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"well_","options for well",""); CHKERRQ(ierr);
    ierr = PetscOptionsEnum("-scheme","finite difference scheme type",
           "well.c",SchemeTypes,
           (PetscEnum)user.scheme,(PetscEnum*)&user.scheme,NULL); CHKERRQ(ierr);
    ierr = PetscOptionsBool("-randominit","initialize u,p with random",
           "well.c",randominit,&randominit,NULL);CHKERRQ(ierr);
    ierr = PetscOptionsBool("-shorterrors","abbreviated error output (e.g. for regression testing)",
           "well.c",shorterrors,&shorterrors,NULL);CHKERRQ(ierr);
    ierr = PetscOptionsEnd(); CHKERRQ(ierr);

    ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,3,2,1,NULL,&da); CHKERRQ(ierr);
    ierr = DMSetFromOptions(da); CHKERRQ(ierr);
    ierr = DMSetUp(da); CHKERRQ(ierr);
    ierr = DMSetApplicationContext(da,&user); CHKERRQ(ierr);
    ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);
    ierr = DMDASetUniformCoordinates(da,0.0,user.L,0.0,1.0,0.0,1.0); CHKERRQ(ierr);
    ierr = DMDASetCoordinateName(da,0,"x"); CHKERRQ(ierr);
    ierr = DMDASetFieldName(da,0,"velocity"); CHKERRQ(ierr);
    ierr = DMDASetFieldName(da,1,"pressure"); CHKERRQ(ierr);

    ierr = DMCreateGlobalVector(da,&X); CHKERRQ(ierr);

    ierr = SNESCreate(PETSC_COMM_WORLD,&snes); CHKERRQ(ierr);
    ierr = SNESSetDM(snes,da); CHKERRQ(ierr);
    if (user.scheme == STAGGERED) {
        ierr = DMDASNESSetFunctionLocal(da,INSERT_VALUES,
                 (DMDASNESFunction)FormFunctionStaggeredLocal,&user); CHKERRQ(ierr);
        ierr = DMDASNESSetJacobianLocal(da,
                 (DMDASNESJacobian)FormJacobianStaggeredLocal,&user); CHKERRQ(ierr);
    } else if (user.scheme == REGULAR) {
        ierr = DMDASNESSetFunctionLocal(da,INSERT_VALUES,
                 (DMDASNESFunction)FormFunctionRegularLocal,&user); CHKERRQ(ierr);
        ierr = DMDASNESSetJacobianLocal(da,
                 (DMDASNESJacobian)FormJacobianRegularLocal,&user); CHKERRQ(ierr);
    } else {
        SETERRQ(PETSC_COMM_WORLD,1,"no way to get here");
    }
    // set defaults to -ksp_type preonly -pc_type svd ... which does not scale or parallelize but is robust
    ierr = SNESGetKSP(snes,&ksp); CHKERRQ(ierr);
    ierr = KSPSetType(ksp,KSPPREONLY); CHKERRQ(ierr);
    ierr = KSPGetPC(ksp,&pc); CHKERRQ(ierr);
    ierr = PCSetType(pc,PCSVD); CHKERRQ(ierr);
    ierr = SNESSetFromOptions(snes); CHKERRQ(ierr);

    if (randominit) {
        PetscRandom   rctx;
        ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rctx); CHKERRQ(ierr);
        ierr = PetscRandomSetSeed(rctx,time(NULL)); CHKERRQ(ierr);
        ierr = PetscRandomSeed(rctx); CHKERRQ(ierr);
        ierr = VecSetRandom(X,rctx); CHKERRQ(ierr);
        ierr = PetscRandomDestroy(&rctx); CHKERRQ(ierr);
    } else {
        ierr = VecSet(X,0.0); CHKERRQ(ierr);
    }
//ierr = VecView(X,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);

    ierr = SNESSolve(snes,NULL,X); CHKERRQ(ierr);

    ierr = VecDuplicate(X,&Xexact); CHKERRQ(ierr);
    ierr = ExactSolution(&info,Xexact,&user); CHKERRQ(ierr);

//ierr = VecView(X,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
//ierr = VecView(Xexact,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);

    ierr = VecAXPY(X,-1.0,Xexact); CHKERRQ(ierr);    // X <- X + (-1.0) Xexact
    ierr = VecStrideNorm(X,0,NORM_INFINITY,&uerrnorm); CHKERRQ(ierr);
    ierr = VecStrideNorm(X,1,NORM_INFINITY,&perrnorm); CHKERRQ(ierr);
    ierr = VecStrideNorm(Xexact,1,NORM_INFINITY,&pnorm); CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,
           "on %d point grid with h=%g and scheme = '%s':\n",
           info.mx,user.L/(info.mx-1),SchemeTypes[user.scheme]); CHKERRQ(ierr);
    if (shorterrors && uerrnorm < setol && perrnorm/pnorm < setol) {
        ierr = PetscPrintf(PETSC_COMM_WORLD,
               "  |u-uexact|_inf < %.1e,  |p-pexact|_inf / |pexact|_inf < %.1e\n",
               setol,setol); CHKERRQ(ierr);
    } else {
        ierr = PetscPrintf(PETSC_COMM_WORLD,
               "  |u-uexact|_inf = %.3e,  |p-pexact|_inf / |pexact|_inf = %.3e\n",
               uerrnorm,perrnorm/pnorm); CHKERRQ(ierr);
    }

    VecDestroy(&X);  VecDestroy(&Xexact);
    SNESDestroy(&snes);  DMDestroy(&da);
    return PetscFinalize();
}


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

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_movstat.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_statistics.h>

typedef double qqracc_type_t;
typedef qqracc_type_t ringbuf_type_t;

#include "ringbuf.c"

typedef struct
{
  qqracc_type_t *window; /* linear array for current window */
  ringbuf *rbuf;         /* ring buffer storing current window */
} qqracc_state_t;

static size_t
qqracc_size(const size_t n)
{
  size_t size = 0;

  size += sizeof(qqracc_state_t);
  size += n * sizeof(qqracc_type_t);
  size += ringbuf_size(n);

  return size;
}

static int
qqracc_init(const size_t n, void * vstate)
{
  qqracc_state_t * state = (qqracc_state_t *) vstate;

  state->window = (qqracc_type_t *) ((unsigned char *) vstate + sizeof(qqracc_state_t));
  state->rbuf = (ringbuf *) ((unsigned char *) state->window + n * sizeof(qqracc_type_t));

  ringbuf_init(n, state->rbuf);

  return GSL_SUCCESS;
}

static int
qqracc_insert(const qqracc_type_t x, void * vstate)
{
  qqracc_state_t * state = (qqracc_state_t *) vstate;

  /* add new element to ring buffer */
  ringbuf_insert(x, state->rbuf);

  return GSL_SUCCESS;
}

static int
qqracc_delete(void * vstate)
{
  qqracc_state_t * state = (qqracc_state_t *) vstate;

  if (!ringbuf_is_empty(state->rbuf))
    ringbuf_pop_back(state->rbuf);

  return GSL_SUCCESS;
}

/* FIXME XXX: this is inefficient - could be improved by maintaining a sorted ring buffer */
static int
qqracc_get(void * params, qqracc_type_t * result, const void * vstate)
{
  const qqracc_state_t * state = (const qqracc_state_t *) vstate;
  double q = *(double *) params;
  size_t n = ringbuf_copy(state->window, state->rbuf);
  double quant1, quant2;

  gsl_sort(state->window, 1, n);

  /* compute q-quantile and (1-q)-quantile */
  quant1 = gsl_stats_quantile_from_sorted_data(state->window, 1, n, q);
  quant2 = gsl_stats_quantile_from_sorted_data(state->window, 1, n, 1.0 - q);

  /* compute q-quantile range */
  *result = quant2 - quant1;

  return GSL_SUCCESS;
}

static const gsl_movstat_accum qqr_accum_type =
{
  qqracc_size,
  qqracc_init,
  qqracc_insert,
  qqracc_delete,
  qqracc_get
};

const gsl_movstat_accum *gsl_movstat_accum_qqr = &qqr_accum_type;

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
/* Copyright 2016 The TensorFlow Authors. All Rights Reserved.

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.
==============================================================================*/

// This is a set of different implementations for the basic matrix by matrix
// multiply function, commonly known as GEMM after the BLAS library's naming.
// Having a standard interface enables us to swap out implementations on
// different platforms, to make sure we're using the optimal version. They are
// implemented as C++ template functors, so they're easy to swap into all of the
// different kernels that use them.

#if !defined(EIGEN_USE_THREADS)
#error "EIGEN_USE_THREADS must be enabled by all .cc files including this."
#endif  // EIGEN_USE_THREADS

#include <string.h>
#include <map>
#include <vector>

#include "tensorflow/core/common_runtime/threadpool_device.h"
#include "tensorflow/core/framework/op_kernel.h"
#include "tensorflow/core/framework/tensor.h"
#include "tensorflow/core/framework/tensor_types.h"

// Apple provides an optimized BLAS library that is better than Eigen for their
// devices, so use that if possible.
#if defined(__APPLE__) && defined(USE_GEMM_FOR_CONV)
#include <Accelerate/Accelerate.h>
#define USE_CBLAS_GEMM
#endif  // __APPLE__

// Older Raspberry Pi systems don't have NEON SIMD acceleration, so Eigen falls
// back to scalar code, but OpenBLAS has much faster support so prefer that.
#if defined(RASPBERRY_PI) && defined(USE_GEMM_FOR_CONV) && defined(USE_OPENBLAS)
#include <cblas.h>
#define USE_CBLAS_GEMM
#endif

// A readable but slow implementation of matrix multiplication, useful for
// debugging and understanding the algorithm. Use instead of FastGemmFunctor in
// the Im2ColConvFunctor template definition inside the op registration to
// enable. Assumes row-major ordering of the values in memory.
template <class T1, class T2, class T3>
class ReferenceGemmFunctor {
 public:
  void operator()(tensorflow::OpKernelContext* ctx, size_t m, size_t n,
                  size_t k, const T1* a, size_t lda, const T2* b, size_t ldb,
                  T3* c, size_t ldc) {
    const size_t a_i_stride = lda;
    const size_t a_l_stride = 1;
    const size_t b_j_stride = 1;
    const size_t b_l_stride = ldb;
    const size_t c_i_stride = ldc;
    const size_t c_j_stride = 1;
    size_t i, j, l;
    for (j = 0; j < n; j++) {
      for (i = 0; i < m; i++) {
        T3 total(0);
        for (l = 0; l < k; l++) {
          const size_t a_index = ((i * a_i_stride) + (l * a_l_stride));
          const T1 a_value = a[a_index];
          const size_t b_index = ((j * b_j_stride) + (l * b_l_stride));
          const T2 b_value = b[b_index];
          total += (a_value * b_value);
        }
        const size_t c_index = ((i * c_i_stride) + (j * c_j_stride));
        c[c_index] = total;
      }
    }
  }
};

// Uses the optimized EigenTensor library to implement the matrix multiplication
// required by the Im2ColConvFunctor class. We supply the two input and one
// output types so that the accumulator can potentially be higher-precision than
// the inputs, even though we don't currently take advantage of this.
template <class T1, class T2, class T3>
class FastGemmFunctor {
 public:
  void operator()(tensorflow::OpKernelContext* ctx, size_t m, size_t n,
                  size_t k, const T1* a, size_t lda, const T2* b, size_t ldb,
                  T3* c, size_t ldc) {
    typename tensorflow::TTypes<const T1>::Matrix a_matrix(a, m, k);
    typename tensorflow::TTypes<const T2>::Matrix b_matrix(b, k, n);
    typename tensorflow::TTypes<T3>::Matrix c_matrix(c, m, n);

    Eigen::array<Eigen::IndexPair<Eigen::DenseIndex>, 1> dim_pair;
    dim_pair[0].first = 1;
    dim_pair[0].second = 0;
    c_matrix.device(ctx->eigen_device<Eigen::ThreadPoolDevice>()) =
        a_matrix.contract(b_matrix, dim_pair);
  }
};

// If we have a fast CBLAS library, use its implementation through a wrapper.
#if defined(USE_CBLAS_GEMM)
template <>
class FastGemmFunctor<float, float, float> {
 public:
  void operator()(tensorflow::OpKernelContext* ctx, size_t m, size_t n,
                  size_t k, const float* a, size_t lda, const float* b,
                  size_t ldb, float* c, size_t ldc) {
    cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, n, k, 1.0f, a,
                lda, b, ldb, 0.0f, c, ldc);
  }
};
#endif  // USE_CBLAS_GEMM

max_stars_repo_name: tianyapiaozi/tensorflow

max_issues_repo_name: shrikunjsarda/tensorflow

max_forks_repo_name: shrikunjsarda/tensorflow


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

#include <tree_2k.h>
#include <tree_2k_utils.h>

#include "cl_params.h"

void tree_spatial_dims_alloc(Params params, double **center,
                             double **extent);
void tree_spatial_dims_free(double *center, double *extent);
void insert_points(tree_2k_t *tree, Params params);

int main(int argc, char *argv[]) {
    tree_2k_err_t status;
    double *center, *extent;
    Params params;
    initCL(&params);
    parseCL(&params, &argc, &argv);
    if (params.verbose)
        dumpCL(stderr, "# ", &params);
    tree_spatial_dims_alloc(params, &center, &extent);
    tree_2k_t * tree;
    status = tree_2k_alloc(&tree, params.rank, center, extent,
                           params.nr_points, params.bucket_size);
    if (status != TREE_2K_SUCCESS)
        errx(EXIT_FAILURE, "can not allocate tree");
    insert_points(tree, params);
    status = tree_2k_fwrite(tree, stdout);
    if (status != TREE_2K_SUCCESS)
        errx(EXIT_FAILURE, "tree write failed");
    tree_2k_free(&tree);
    tree_spatial_dims_free(center, extent);
    finalizeCL(&params);
    return EXIT_SUCCESS;
}

void tree_spatial_dims_alloc(Params params, double **center,
                             double **extent) {
    *center = (double *) malloc(params.rank*sizeof(double));
    *extent = (double *) malloc(params.rank*sizeof(double));
    if (*center == NULL || *extent == NULL)
        errx(EXIT_FAILURE, "can not allocate center or extent");
    for (int i = 0; i < params.rank; i++) {
        (*center)[i] = 0.0;
        (*extent)[i] = 1.0;
    }
}

void tree_spatial_dims_free(double *center, double *extent) {
    free(center);
    free(extent);
}

void insert_points(tree_2k_t *tree, Params params) {
    double *coords;
    gsl_rng *rng;
    gsl_rng_env_setup();
    rng = gsl_rng_alloc(gsl_rng_default);
    if ((coords = (double *) malloc(params.rank*sizeof(double))) == NULL)
        errx(EXIT_FAILURE, "can not allocate coords");
    for (int point_nr = 0; point_nr < params.nr_points; point_nr++) {
        for (int i = 0; i < params.rank; i++)
            coords[i] = 2.0*gsl_rng_uniform(rng) - 1.0;
        tree_2k_err_t status = tree_2k_insert(tree, coords, NULL);
        if (status != TREE_2K_SUCCESS)
            errx(EXIT_FAILURE, "failed to insert point %d", point_nr);
    }
    gsl_rng_free(rng);
}

max_stars_repo_name: Gjacquenot/training-material

max_issues_repo_name: Gjacquenot/training-material

max_forks_repo_name: Gjacquenot/training-material


text:
#pragma once

#include <algorithm>
#include <cstddef>
#include <cstdint>
#include <tuple>
#include <type_traits>

#include <gsl/gsl>

#include "chainerx/array.h"
#include "chainerx/axes.h"
#include "chainerx/backend_util.h"
#include "chainerx/constant.h"
#include "chainerx/dtype.h"
#include "chainerx/indexer.h"
#include "chainerx/macro.h"
#include "chainerx/shape.h"
#include "chainerx/strides.h"

namespace chainerx {
namespace indexable_array_detail {

template <typename To, typename From>
using WithConstnessOf = dtype_detail::WithConstnessOf<To, From>;

static inline std::tuple<const uint8_t*, const uint8_t*> GetDataRange(const Array& a) {
    std::tuple<int64_t, int64_t> range = chainerx::GetDataRange(a.shape(), a.strides(), a.GetItemSize());
    int64_t lower = std::get<0>(range);
    int64_t upper = std::get<1>(range);
    const uint8_t* base = static_cast<const uint8_t*>(internal::GetRawOffsetData(a));
    return std::tuple<const uint8_t*, const uint8_t*>{base + lower, base + upper};
}

}  // namespace indexable_array_detail

template <typename T, int8_t kNdim = kDynamicNdim>
class IndexableArray {
public:
    using ElementType = T;

private:
    template <typename U>
    using WithConstnessOfT = indexable_array_detail::WithConstnessOf<U, T>;
    using VoidType = WithConstnessOfT<void>;
    using DeviceStorageType = TypeToDeviceStorageType<T>;

public:
    IndexableArray(VoidType* data, const Strides& strides) : data_{data} { std::copy(strides.begin(), strides.end(), strides_); }

    IndexableArray(const Array& array, const Strides& strides) : IndexableArray{internal::GetRawOffsetData(array), strides} {
        CHAINERX_ASSERT(TypeToDtype<T> == array.dtype());

#if CHAINERX_DEBUG
        std::tie(first_, last_) = indexable_array_detail::GetDataRange(array);
#endif  // CHAINERX_DEBUG
    }

    explicit IndexableArray(const Array& array) : IndexableArray{array, array.strides()} {}

    CHAINERX_HOST_DEVICE int8_t ndim() const { return kNdim; }

    CHAINERX_HOST_DEVICE const int64_t* strides() const { return strides_; }

    CHAINERX_HOST_DEVICE VoidType* data() const { return data_; }

    CHAINERX_HOST_DEVICE DeviceStorageType& operator[](const int64_t* index) const {
        auto data_ptr = static_cast<WithConstnessOfT<uint8_t>*>(data_);
        for (int8_t dim = 0; dim < kNdim; ++dim) {
            data_ptr += strides_[dim] * index[dim];
        }
#if CHAINERX_DEBUG
        CHAINERX_ASSERT(first_ == nullptr || first_ <= data_ptr);
        CHAINERX_ASSERT(last_ == nullptr || data_ptr <= last_ - sizeof(T));
#endif  // CHAINERX_DEBUG
        return *static_cast<DeviceStorageType*>(static_cast<VoidType*>(data_ptr));
    }

    CHAINERX_HOST_DEVICE WithConstnessOfT<DeviceStorageType>& operator[](const IndexIterator<kNdim>& it) const {
        return operator[](it.index());
    }

    // Permutes the axes.
    //
    // It is the caller's responsibility to ensure validity of permutation.
    // If the permutation is invalid, the behavior is undefined.
    IndexableArray<T, kNdim>& Permute(const Axes& axes) {
        CHAINERX_ASSERT(axes.size() == static_cast<size_t>(kNdim));
        int64_t c[kNdim]{};
        std::copy(std::begin(strides_), std::end(strides_), c);
        for (size_t i = 0; i < kNdim; ++i) {
            strides_[i] = c[axes[i]];
        }
        return *this;
    }

private:
    WithConstnessOfT<void>* data_;
#if CHAINERX_DEBUG
    const uint8_t* first_{nullptr};
    const uint8_t* last_{nullptr};
#endif  // CHAINERX_DEBUG
    int64_t strides_[kNdim];
};

// Static 0-dimensional specialization.
template <typename T>
class IndexableArray<T, 0> {
public:
    using ElementType = T;

private:
    template <typename U>
    using WithConstnessOfT = indexable_array_detail::WithConstnessOf<U, T>;
    using VoidType = WithConstnessOfT<void>;
    using DeviceStorageType = TypeToDeviceStorageType<T>;

public:
    IndexableArray(VoidType* data, const Strides& strides) : data_{data} { CHAINERX_ASSERT(0 == strides.ndim()); }

    IndexableArray(const Array& array, const Strides& strides) : IndexableArray{internal::GetRawOffsetData(array), strides} {
        CHAINERX_ASSERT(TypeToDtype<T> == array.dtype());
    }

    explicit IndexableArray(const Array& array) : IndexableArray{array, array.strides()} {}

    CHAINERX_HOST_DEVICE constexpr int8_t ndim() const { return 0; }

    CHAINERX_HOST_DEVICE constexpr const int64_t* strides() const { return nullptr; }

    CHAINERX_HOST_DEVICE VoidType* data() const { return data_; }

    CHAINERX_HOST_DEVICE DeviceStorageType& operator[](const int64_t* index) const {
        CHAINERX_ASSERT(index == nullptr || index[0] == 0);
        return *static_cast<WithConstnessOfT<DeviceStorageType>*>(data_);
    }

    CHAINERX_HOST_DEVICE DeviceStorageType& operator[](const IndexIterator<0>& it) const { return operator[](it.index()); }

    IndexableArray<T, 0>& Permute(const Axes& /*axes*/) {
        // NOOP for 1-dimensional array.
        return *this;
    }

private:
    WithConstnessOfT<void>* data_;
};

// Static 1-dimensional specialization.
template <typename T>
class IndexableArray<T, 1> {
public:
    using ElementType = T;

private:
    template <typename U>
    using WithConstnessOfT = indexable_array_detail::WithConstnessOf<U, T>;
    using VoidType = WithConstnessOfT<void>;
    using DeviceStorageType = TypeToDeviceStorageType<T>;

public:
    IndexableArray(VoidType* data, const Strides& strides) : data_{data}, stride_{strides[0]} { CHAINERX_ASSERT(1 == strides.ndim()); }

    IndexableArray(const Array& array, const Strides& strides) : IndexableArray{internal::GetRawOffsetData(array), strides} {
        CHAINERX_ASSERT(TypeToDtype<T> == array.dtype());

#if CHAINERX_DEBUG
        std::tie(first_, last_) = indexable_array_detail::GetDataRange(array);
#endif  // CHAINERX_DEBUG
    }

    explicit IndexableArray(const Array& array) : IndexableArray{array, array.strides()} {}

    CHAINERX_HOST_DEVICE constexpr int8_t ndim() const { return 1; }

    CHAINERX_HOST_DEVICE const int64_t* strides() const { return &stride_; }

    CHAINERX_HOST_DEVICE VoidType* data() const { return data_; }

    CHAINERX_HOST_DEVICE DeviceStorageType& operator[](const int64_t* index) const {
        auto data_ptr = reinterpret_cast<WithConstnessOfT<uint8_t>*>(data_) +
                        stride_ * index[0];  // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast)
#if CHAINERX_DEBUG
        CHAINERX_ASSERT(first_ == nullptr || first_ <= data_ptr);
        CHAINERX_ASSERT(last_ == nullptr || data_ptr <= last_ - sizeof(T));
#endif  // CHAINERX_DEBUG
        return *static_cast<DeviceStorageType*>(static_cast<VoidType*>(data_ptr));
    }

    CHAINERX_HOST_DEVICE DeviceStorageType& operator[](const IndexIterator<1>& it) const { return operator[](it.index()); }

    IndexableArray<T, 1>& Permute(const Axes& /*axes*/) {
        // NOOP for 1-dimensional array.
        return *this;
    }

private:
    WithConstnessOfT<void>* data_;
#if CHAINERX_DEBUG
    const uint8_t* first_{nullptr};
    const uint8_t* last_{nullptr};
#endif  // CHAINERX_DEBUG
    int64_t stride_{};
};

// Runtime determined dynamic dimension specialization.
template <typename T>
class IndexableArray<T, kDynamicNdim> {
public:
    using ElementType = T;

private:
    template <typename U>
    using WithConstnessOfT = indexable_array_detail::WithConstnessOf<U, T>;
    using VoidType = WithConstnessOfT<void>;
    using DeviceStorageType = TypeToDeviceStorageType<T>;

public:
    IndexableArray(WithConstnessOfT<void>* data, const Strides& strides) : data_{data}, ndim_{strides.ndim()} {
        std::copy(strides.begin(), strides.end(), strides_);
    }

    IndexableArray(const Array& array, const Strides& strides) : IndexableArray{internal::GetRawOffsetData(array), strides} {
        CHAINERX_ASSERT(TypeToDtype<T> == array.dtype());

#if CHAINERX_DEBUG
        std::tie(first_, last_) = indexable_array_detail::GetDataRange(array);
#endif  // CHAINERX_DEBUG
    }

    explicit IndexableArray(const Array& array) : IndexableArray{array, array.strides()} {}

    CHAINERX_HOST_DEVICE int8_t ndim() const { return ndim_; }

    CHAINERX_HOST_DEVICE const int64_t* strides() const { return strides_; }

    CHAINERX_HOST_DEVICE VoidType* data() const { return data_; }

    CHAINERX_HOST_DEVICE DeviceStorageType& operator[](const int64_t* index) const {
        auto data_ptr = static_cast<WithConstnessOfT<uint8_t>*>(data_);
        for (int8_t dim = 0; dim < ndim_; ++dim) {
            data_ptr += strides_[dim] * index[dim];
        }
#if CHAINERX_DEBUG
        CHAINERX_ASSERT(first_ == nullptr || first_ <= data_ptr);
        CHAINERX_ASSERT(last_ == nullptr || data_ptr <= last_ - sizeof(T));
#endif  // CHAINERX_DEBUG
        return *static_cast<DeviceStorageType*>(static_cast<VoidType*>(data_ptr));
    }

    CHAINERX_HOST_DEVICE DeviceStorageType& operator[](const IndexIterator<kDynamicNdim>& it) const { return operator[](it.index()); }

    // Permutes the axes.
    //
    // Given axes may be fewer than that held by the array.
    // In that case, the axes in the array will be reduced.
    //
    // It is the caller's responsibility to ensure validity of permutation.
    // If the permutation is invalid, the behavior is undefined.
    IndexableArray<T, kDynamicNdim>& Permute(const Axes& axes) {
        CHAINERX_ASSERT(axes.size() <= static_cast<size_t>(ndim_));
        int64_t c[kMaxNdim]{};
        std::copy(std::begin(strides_), std::end(strides_), c);
        for (size_t i = 0; i < axes.size(); ++i) {
            strides_[i] = c[axes[i]];
        }
        ndim_ = static_cast<int8_t>(axes.size());
        return *this;
    }

private:
    WithConstnessOfT<void>* data_;
#if CHAINERX_DEBUG
    const uint8_t* first_{nullptr};
    const uint8_t* last_{nullptr};
#endif  // CHAINERX_DEBUG
    int64_t strides_[kMaxNdim];
    int8_t ndim_;
};

}  // namespace chainerx

max_stars_repo_name: hikjik/chainer

max_issues_repo_name: nolfwin/chainer

max_forks_repo_name: nolfwin/chainer


text:
// Copyright 2021 atframework
// Created by owent on 2020-08-18

#pragma once

#include <gsl/select-gsl.h>

#include <config/atframe_utils_build_feature.h>
#include <config/compiler_features.h>

#include <design_pattern/nomovable.h>
#include <design_pattern/noncopyable.h>

#include <random/random_generator.h>

#include "atframe/atapp_conf.h"

namespace atapp {
struct LIBATAPP_MACRO_API_HEAD_ONLY etcd_discovery_action_t {
  enum type {
    EN_NAT_UNKNOWN = 0,
    EN_NAT_PUT,
    EN_NAT_DELETE,
  };
};

class etcd_discovery_node {
 public:
  using on_destroy_fn_t = std::function<void(etcd_discovery_node &)>;
  using ptr_t = std::shared_ptr<etcd_discovery_node>;

  UTIL_DESIGN_PATTERN_NOCOPYABLE(etcd_discovery_node)
  UTIL_DESIGN_PATTERN_NOMOVABLE(etcd_discovery_node)

 public:
  LIBATAPP_MACRO_API etcd_discovery_node();
  LIBATAPP_MACRO_API ~etcd_discovery_node();

  UTIL_FORCEINLINE const atapp::protocol::atapp_discovery &get_discovery_info() const { return node_info_; }
  LIBATAPP_MACRO_API void copy_from(const atapp::protocol::atapp_discovery &input);

  UTIL_FORCEINLINE const std::pair<uint64_t, uint64_t> &get_name_hash() const { return name_hash_; }

  UTIL_FORCEINLINE void set_private_data_ptr(void *input) { private_data_ptr_ = input; }
  UTIL_FORCEINLINE void *get_private_data_ptr() const { return private_data_ptr_; }
  UTIL_FORCEINLINE void set_private_data_u64(uint64_t input) { private_data_u64_ = input; }
  UTIL_FORCEINLINE uint64_t get_private_data_u64() const { return private_data_u64_; }
  UTIL_FORCEINLINE void set_private_data_i64(int64_t input) { private_data_i64_ = input; }
  UTIL_FORCEINLINE int64_t get_private_data_i64() const { return private_data_i64_; }
  UTIL_FORCEINLINE void set_private_data_uptr(uintptr_t input) { private_data_uptr_ = input; }
  UTIL_FORCEINLINE uintptr_t get_private_data_uptr() const { return private_data_uptr_; }
  UTIL_FORCEINLINE void set_private_data_iptr(intptr_t input) { private_data_iptr_ = input; }
  UTIL_FORCEINLINE intptr_t get_private_data_iptr() const { return private_data_iptr_; }

  LIBATAPP_MACRO_API void set_on_destroy(on_destroy_fn_t fn);
  LIBATAPP_MACRO_API const on_destroy_fn_t &get_on_destroy() const;
  LIBATAPP_MACRO_API void reset_on_destroy();

  LIBATAPP_MACRO_API const atapp::protocol::atapp_gateway &next_ingress_gateway() const;
  LIBATAPP_MACRO_API int32_t get_ingress_size() const;

 private:
  atapp::protocol::atapp_discovery node_info_;
  std::pair<uint64_t, uint64_t> name_hash_;
  union {
    void *private_data_ptr_;
    uint64_t private_data_u64_;
    int64_t private_data_i64_;
    uintptr_t private_data_uptr_;
    intptr_t private_data_iptr_;
  };
  on_destroy_fn_t on_destroy_fn_;
  mutable int32_t ingress_index_;
  mutable atapp::protocol::atapp_gateway ingress_for_listen_;
};

class etcd_discovery_set {
 public:
  using node_by_name_t = std::unordered_map<std::string, etcd_discovery_node::ptr_t>;
  using node_by_id_t = std::unordered_map<uint64_t, etcd_discovery_node::ptr_t>;
  using ptr_t = std::shared_ptr<etcd_discovery_set>;

  struct node_hash_t {
    enum { HASH_POINT_PER_INS = 80 };

    etcd_discovery_node::ptr_t node;
    std::pair<uint64_t, uint64_t> hash_code;
  };

  UTIL_DESIGN_PATTERN_NOCOPYABLE(etcd_discovery_set)
  UTIL_DESIGN_PATTERN_NOMOVABLE(etcd_discovery_set)
 public:
  LIBATAPP_MACRO_API etcd_discovery_set();
  LIBATAPP_MACRO_API ~etcd_discovery_set();

  LIBATAPP_MACRO_API bool empty() const;

  LIBATAPP_MACRO_API etcd_discovery_node::ptr_t get_node_by_id(uint64_t id) const;
  LIBATAPP_MACRO_API etcd_discovery_node::ptr_t get_node_by_name(gsl::string_view name) const;

  LIBATAPP_MACRO_API etcd_discovery_node::ptr_t get_node_by_consistent_hash(const void *buf, size_t bufsz) const;
  LIBATAPP_MACRO_API etcd_discovery_node::ptr_t get_node_by_consistent_hash(uint64_t key) const;
  LIBATAPP_MACRO_API etcd_discovery_node::ptr_t get_node_by_consistent_hash(int64_t key) const;
  LIBATAPP_MACRO_API etcd_discovery_node::ptr_t get_node_by_consistent_hash(gsl::string_view key) const;

  LIBATAPP_MACRO_API etcd_discovery_node::ptr_t get_node_by_random() const;
  LIBATAPP_MACRO_API etcd_discovery_node::ptr_t get_node_by_round_robin() const;

  LIBATAPP_MACRO_API const std::vector<etcd_discovery_node::ptr_t> &get_sorted_nodes() const;
  LIBATAPP_MACRO_API std::vector<etcd_discovery_node::ptr_t>::const_iterator lower_bound_sorted_nodes(
      uint64_t id, gsl::string_view name) const;
  LIBATAPP_MACRO_API std::vector<etcd_discovery_node::ptr_t>::const_iterator upper_bound_sorted_nodes(
      uint64_t id, gsl::string_view name) const;

  LIBATAPP_MACRO_API void add_node(const etcd_discovery_node::ptr_t &node);
  LIBATAPP_MACRO_API void remove_node(const etcd_discovery_node::ptr_t &node);
  LIBATAPP_MACRO_API void remove_node(uint64_t id);
  LIBATAPP_MACRO_API void remove_node(gsl::string_view name);

 private:
  void rebuild_cache() const;
  void clear_cache() const;

 private:
  node_by_name_t node_by_name_;
  node_by_id_t node_by_id_;

  mutable std::vector<node_hash_t> hashing_cache_;
  mutable std::vector<etcd_discovery_node::ptr_t> round_robin_cache_;
  mutable util::random::xoshiro256_starstar random_generator_;
  mutable size_t round_robin_index_;
};
}  // namespace atapp

max_stars_repo_name: atframework/libatapp

max_issues_repo_name: atframework/libatapp

max_forks_repo_name: atframework/libatapp


text:
// include all the libraries that are needed in our simulations

#include <iostream>
#include <vector>
#include <cmath>
#include <new>
#include <math.h>
#include <time.h>
#include <algorithm>
#include <fstream>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <string>
#include <sstream>
#include <string.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_statistics_double.h>

#include "theConstants.h"

using namespace std;

max_stars_repo_name: MattHartfield/FacSexBGSSims

max_issues_repo_name: MattHartfield/FacSexBGSSims

max_forks_repo_name: MattHartfield/FacSexBGSSims


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

/* Author:  G. Jungman */

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

#define REAL double

/* [Engeln-Mullges + Uhlig, Alg. 4.42]
 */

int
gsl_linalg_HH_solve (gsl_matrix * A, const gsl_vector * b, gsl_vector * x)
{
  if (A->size1 > A->size2)
    {
      /* System is underdetermined. */

      GSL_ERROR ("System is underdetermined", GSL_EINVAL);
    }
  else if (A->size2 != x->size)
    {
      GSL_ERROR ("matrix and vector sizes must be equal", GSL_EBADLEN);
    }
  else
    {
      int status ;

      gsl_vector_memcpy (x, b);

      status = gsl_linalg_HH_svx (A, x);

      return status ;
    }
}

int
gsl_linalg_HH_svx (gsl_matrix * A, gsl_vector * x)
{
  if (A->size1 > A->size2)
    {
      /* System is underdetermined. */

      GSL_ERROR ("System is underdetermined", GSL_EINVAL);
    }
  else if (A->size2 != x->size)
    {
      GSL_ERROR ("matrix and vector sizes must be equal", GSL_EBADLEN);
    }
  else
    {
      const size_t N = A->size1;
      const size_t M = A->size2;
      size_t i, j, k;
      REAL *d = (REAL *) malloc (N * sizeof (REAL));

      if (d == 0)
        {
          GSL_ERROR ("could not allocate memory for workspace", GSL_ENOMEM);
        }

      /* Perform Householder transformation. */

      for (i = 0; i < N; i++)
        {
          const REAL aii = gsl_matrix_get (A, i, i);
          REAL alpha;
          REAL f;
          REAL ak;
          REAL max_norm = 0.0;
          REAL r = 0.0;

          for (k = i; k < M; k++)
            {
              REAL aki = gsl_matrix_get (A, k, i);
              r += aki * aki;
            }

          if (r == 0.0)
            {
              /* Rank of matrix is less than size1. */
              free (d);
              GSL_ERROR ("matrix is rank deficient", GSL_ESING);
            }

          alpha = sqrt (r) * GSL_SIGN (aii);

          ak = 1.0 / (r + alpha * aii);
          gsl_matrix_set (A, i, i, aii + alpha);

          d[i] = -alpha;

          for (k = i + 1; k < N; k++)
            {
              REAL norm = 0.0;
              f = 0.0;
              for (j = i; j < M; j++)
                {
                  REAL ajk = gsl_matrix_get (A, j, k);
                  REAL aji = gsl_matrix_get (A, j, i);
                  norm += ajk * ajk;
                  f += ajk * aji;
                }
              max_norm = GSL_MAX (max_norm, norm);

              f *= ak;

              for (j = i; j < M; j++)
                {
                  REAL ajk = gsl_matrix_get (A, j, k);
                  REAL aji = gsl_matrix_get (A, j, i);
                  gsl_matrix_set (A, j, k, ajk - f * aji);
                }
            }

          if (fabs (alpha) < 2.0 * GSL_DBL_EPSILON * sqrt (max_norm))
            {
              /* Apparent singularity. */
              free (d);
              GSL_ERROR("apparent singularity detected", GSL_ESING);
            }

          /* Perform update of RHS. */

          f = 0.0;
          for (j = i; j < M; j++)
            {
              f += gsl_vector_get (x, j) * gsl_matrix_get (A, j, i);
            }
          f *= ak;
          for (j = i; j < M; j++)
            {
              REAL xj = gsl_vector_get (x, j);
              REAL aji = gsl_matrix_get (A, j, i);
              gsl_vector_set (x, j, xj - f * aji);
            }
        }

      /* Perform back-substitution. */

      for (i = N; i > 0 && i--;)
        {
          REAL xi = gsl_vector_get (x, i);
          REAL sum = 0.0;
          for (k = i + 1; k < N; k++)
            {
              sum += gsl_matrix_get (A, i, k) * gsl_vector_get (x, k);
            }

          gsl_vector_set (x, i, (xi - sum) / d[i]);
        }

      free (d);
      return GSL_SUCCESS;
    }
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified