text:
/*
 *  A rng module implementation which can return a sample as an array.
 *  Author : Pierre Schnizer <schnizer@users.sourceforge.net>
 *  Date   : 1. July. 2003
 *
 */

#include <pygsl/error_helpers.h>
#include <pygsl/block_helpers.h>
#include <pygsl/utils.h>

#ifdef PyGSL_NO_IMPORT_API
#undef PyGSL_NO_IMPORT_API
#endif
#include <pygsl/rng_helpers.h>
#include <pygsl/rng.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

/*
 * All doc strings
 */
static PyObject *module = NULL;

#include "rng_helpers.c"
#include "rngmodule_docs.h"

static void rng_delete(PyGSL_rng *self);
static PyObject * rng_call(PyGSL_rng *self, PyObject *args);

#undef PyGSL_RNG_Check

#define PyGSL_RNG_Check(op) ((op)->ob_type == &PyGSL_rng_pytype)

static PyObject *		/* on "instance.attr" */
rng_getattr (PyGSL_rng *self, char *name);

PyTypeObject PyGSL_rng_pytype = {
  PyObject_HEAD_INIT(NULL)	/* fix up the type slot in initcrng */
  0,				/* ob_size */
  "PyGSL_rng",			/* tp_name */
  sizeof(PyGSL_rng),	        /* tp_basicsize */
  0,				/* tp_itemsize */

  /* standard methods */
  (destructor)  rng_delete,       /* tp_dealloc  ref-count==0  */
  (printfunc)   0,		   /* tp_print    "print x"     */
  (getattrfunc) rng_getattr,       /* tp_getattr  "x.attr"      */
  (setattrfunc) 0,		   /* tp_setattr  "x.attr=v"    */
  (cmpfunc)     0,		   /* tp_compare  "x > y"       */
  (reprfunc)    0,                 /* tp_repr     `x`, print x  */

  /* type categories */
  0,				/* tp_as_number   +,-,*,/,%,&,>>,pow...*/
  0,				/* tp_as_sequence +,[i],[i:j],len, ...*/
  0,				/* tp_as_mapping  [key], len, ...*/

  /* more methods */
  (hashfunc)     0,		/* tp_hash    "dict[x]" */
  (ternaryfunc)  rng_call,      /* tp_call    "x()"     */
  (reprfunc)     0,             /* tp_str     "str(x)"  */
  (getattrofunc) 0,		/* tp_getattro */
  (setattrofunc) 0,		/* tp_setattro */
  0,				/* tp_as_buffer */
  0L,				/* tp_flags */
  rng_type_doc		/* tp_doc */
};

static PyObject *
PyGSL_rng_init(PyObject *self, PyObject *args, const gsl_rng_type * rng_type)
{

     PyGSL_rng *rng = NULL;

     FUNC_MESS_BEGIN();
     rng =  (PyGSL_rng *) PyObject_NEW(PyGSL_rng, &PyGSL_rng_pytype);
     if(rng == NULL){
	  return NULL;
     }
     if(rng_type == NULL){
	  rng->rng = gsl_rng_alloc(gsl_rng_default);
	  gsl_rng_set(rng->rng, gsl_rng_default_seed);
     }else{
	  rng->rng = gsl_rng_alloc(rng_type);
     }
     FUNC_MESS_END();
     return (PyObject *) rng;
}

#define RNG_ARNG(name)                                                       \
static PyObject* PyGSL_rng_init_ ## name (PyObject *self, PyObject *args)    \
{                                                                            \
     PyObject *tmp = NULL;                                                   \
     FUNC_MESS_BEGIN();                                                      \
     tmp = PyGSL_rng_init(self, args,  gsl_rng_ ## name);                    \
     if (tmp == NULL){                                                       \
	  PyGSL_add_traceback(module, __FILE__, __FUNCTION__, __LINE__);     \
     }                                                                       \
     FUNC_MESS_END();                                                        \
     return tmp;                                                             \
}

#include "rng_list.h"

static PyObject *
PyGSL_rng_init_default(PyObject *self, PyObject *args)
{
     PyObject *tmp = NULL;
     FUNC_MESS_BEGIN();
     tmp = PyGSL_rng_init(self, args,  NULL);
     if (tmp == NULL){
	  PyGSL_add_traceback(module, __FILE__, __FUNCTION__, __LINE__ - 1);
     }
     FUNC_MESS_END();
     return tmp;
}

static void
rng_delete(PyGSL_rng *self)
{
     FUNC_MESS_BEGIN();
     assert(PyGSL_RNG_Check(self));
     if(self->rng != NULL){
	  DEBUG_MESS(5, "Freeing gsl_rng @ %p", self->rng);
	  gsl_rng_free(self->rng);
	  self->rng = NULL;
     }
     DEBUG_MESS(1, " self %p\n",(void *) self);
     PyObject_Del(self);
     self = NULL;
     FUNC_MESS_END();
}

static PyObject *		/* on "instance.uniform()" and "instance()" */
rng_call (PyGSL_rng *self, PyObject *args)
{
     PyObject *tmp;

     FUNC_MESS_BEGIN();
     assert(PyGSL_RNG_Check(self));
     tmp = PyGSL_rng_to_double(self, args, gsl_rng_uniform);
     if(!tmp)
	  PyGSL_add_traceback(module, __FILE__, "rng.__call__", __LINE__ - 2);
     FUNC_MESS_END();
     return tmp;
}

static PyObject *
rng_uniform_pos (PyGSL_rng *self, PyObject *args)
{
     PyObject *tmp;

     FUNC_MESS_BEGIN();
     assert(PyGSL_RNG_Check(self));
     tmp = PyGSL_rng_to_double(self, args, gsl_rng_uniform_pos);
     if(!tmp)
	  PyGSL_add_traceback(module, __FILE__, "rng.uniform_pos", __LINE__-2);
     FUNC_MESS_END();
     return tmp;
}

static PyObject *
rng_uniform_int (PyGSL_rng *self, PyObject *args)
{
     PyObject *tmp = NULL;

     FUNC_MESS_BEGIN();
     assert(PyGSL_RNG_Check(self));
     tmp = PyGSL_rng_ul_to_ulong(self, args, gsl_rng_uniform_int);
     if(!tmp)
	  PyGSL_add_traceback(module, __FILE__, "rng.uniform_int", __LINE__-2);
     FUNC_MESS_END();
     return tmp;
}

static PyObject *
rng_get(PyGSL_rng *self, PyObject *args)
{
     PyObject  *tmp = NULL;

     FUNC_MESS_BEGIN();
     assert(PyGSL_RNG_Check(self));
     tmp = PyGSL_rng_to_ulong(self, args, gsl_rng_get);
     if(!tmp)
	  PyGSL_add_traceback(module, __FILE__, "rng.get", __LINE__ - 2);
     FUNC_MESS_END();
     return tmp;
}

static PyObject *
rng_set(PyGSL_rng *self, PyObject *args)
{
     PyObject *tmp = NULL, *seed = NULL;
     unsigned long int useed;
     int lineno;

     FUNC_MESS_BEGIN();
     assert(PyGSL_RNG_Check(self));
     if(0 == PyArg_ParseTuple(args, "O", &tmp)){
	  lineno = __LINE__; goto fail;
     }
     assert(tmp != NULL);
     seed = PyNumber_Long(tmp);
     if(!seed){lineno = __LINE__ - 1; goto fail;}
     useed =  PyLong_AsUnsignedLong(seed);
     gsl_rng_set(self->rng, useed);

     Py_INCREF(Py_None);
     FUNC_MESS_END();
     return Py_None;

 fail:
     FUNC_MESS("FAIL");
     PyGSL_add_traceback(module, __FILE__, "rng.set", lineno);
     return NULL;
}

static PyObject *
rng_name(PyGSL_rng *self, PyObject *args)
{
     PyObject *tmp = NULL;
     FUNC_MESS_BEGIN();
     assert(PyGSL_RNG_Check(self));
     if(0 == PyArg_ParseTuple(args, ":name"))
	  return NULL;
     tmp = PyString_FromString(gsl_rng_name(self->rng));
     FUNC_MESS_END();
     return tmp;
}

static PyObject *
rng_max(PyGSL_rng *self, PyObject *args)
{
     PyObject *tmp = NULL;
     FUNC_MESS_BEGIN();
     assert(PyGSL_RNG_Check(self));
     if(0 == PyArg_ParseTuple(args, ":max"))
	  return NULL;
     tmp = PyLong_FromUnsignedLong(gsl_rng_max(self->rng));
     FUNC_MESS_END();
     return tmp;
}

static PyObject *
rng_min(PyGSL_rng *self, PyObject *args)
{
     PyObject *tmp = NULL;
     FUNC_MESS_BEGIN();
     assert(PyGSL_RNG_Check(self));
     if(0 == PyArg_ParseTuple(args, ":min"))
	  return NULL;
     tmp = PyLong_FromUnsignedLong(gsl_rng_min(self->rng));
     FUNC_MESS_END();
     return tmp;
}

static PyObject *
rng_clone(PyGSL_rng *self, PyObject *args)
{
     PyGSL_rng * rng = NULL;

     FUNC_MESS_BEGIN();
     assert(PyGSL_RNG_Check(self));
     if(0 == PyArg_ParseTuple(args, ":clone"))
	  return NULL;
     rng =  (PyGSL_rng *) PyObject_NEW(PyGSL_rng, &PyGSL_rng_pytype);
     rng->rng = gsl_rng_clone(self->rng);
     FUNC_MESS_END();
     return (PyObject *) rng;
}

/*
 * Is #name a standard macro definition?
 */
#define RNG_DISTRIBUTION(name, function)                                     \
static PyObject* rng_ ## name (PyGSL_rng *self, PyObject *args)              \
{                                                                            \
     PyObject *tmp = NULL;                                                   \
     FUNC_MESS_BEGIN();                                                      \
     tmp = PyGSL_rng_ ## function (self, args,  gsl_ran_ ## name);           \
     if (tmp == NULL){                                                       \
       /* PyGSL_add_traceback(module, __FILE__, "rng." #name, __LINE__); */  \
	  PyGSL_add_traceback(module, __FILE__, __FUNCTION__, __LINE__);     \
     }                                                                       \
     FUNC_MESS_END();                                                        \
     return tmp;                                                             \
}

#include "rng_distributions.h"

/* Redefine to trigger emacs into correct coloring */
/*
 * This list is not optimized. I guess... Do you know how to optimize it?
 */
static struct PyMethodDef rng_methods[] = {
  {"get", (PyCFunction) rng_get, METH_VARARGS, rng_get_doc},
  {"set", (PyCFunction) rng_set, METH_VARARGS, rng_set_doc},
  {"uniform", (PyCFunction) rng_call, METH_VARARGS, rng_uniform_doc},
  {"uniform_pos", (PyCFunction) rng_uniform_pos, METH_VARARGS, rng_uniform_pos_doc},
  {"uniform_int", (PyCFunction) rng_uniform_int, METH_VARARGS, rng_uniform_int_doc},
  {"name", (PyCFunction) rng_name, METH_VARARGS, NULL},
  {"max", (PyCFunction) rng_max, METH_VARARGS, rng_max_doc},
  {"min", (PyCFunction) rng_min, METH_VARARGS, rng_min_doc},
  {"clone", (PyCFunction) rng_clone, METH_VARARGS, rng_clone_doc},
#if  (PY_MAJOR_VERSION == 2) && (PY_MINOR_VERSION == 3)
  /* RNG clone can not be used to copy a rng type in python 2.3 No idea how to do that correctly */
#else
  {"__copy__", (PyCFunction) rng_clone, METH_VARARGS, rng_clone_doc},
#endif
  /* distributions */
  {"gaussian", (PyCFunction) rng_gaussian,METH_VARARGS, rng_gaussian_doc},
  {"gaussian_ratio_method", (PyCFunction) rng_gaussian_ratio_method,METH_VARARGS, rng_gaussian_ratio_doc},
  {"ugaussian", (PyCFunction) rng_ugaussian,METH_VARARGS, rng_ugaussian_doc},
  {"ugaussian_ratio_method",(PyCFunction) rng_ugaussian_ratio_method,METH_VARARGS, rng_ugaussian_ratio_doc},
  {"gaussian_tail", (PyCFunction) rng_gaussian_tail,METH_VARARGS, rng_gaussian_tail_doc},
  {"ugaussian_tail",(PyCFunction) rng_ugaussian_tail,METH_VARARGS, rng_ugaussian_tail_doc},
  {"bivariate_gaussian", (PyCFunction) rng_bivariate_gaussian,METH_VARARGS, rng_bivariate_gaussian_doc},
  {"exponential",(PyCFunction)rng_exponential,METH_VARARGS, rng_exponential_doc},
  {"laplace",(PyCFunction)rng_laplace,METH_VARARGS, rng_laplace_doc},
  {"exppow",(PyCFunction)rng_exppow,METH_VARARGS, rng_exppow_doc},
  {"cauchy",(PyCFunction)rng_cauchy,METH_VARARGS, rng_cauchy_doc},
  {"rayleigh",(PyCFunction)rng_rayleigh,METH_VARARGS, rng_rayleigh_doc},
  {"rayleigh_tail",(PyCFunction)rng_rayleigh_tail,METH_VARARGS, rng_rayleigh_tail_doc},
  {"levy",(PyCFunction)rng_levy,METH_VARARGS, rng_levy_doc},
  {"levy_skew",(PyCFunction)rng_levy_skew,METH_VARARGS, rng_levy_skew_doc},
  {"gamma",(PyCFunction)rng_gamma,METH_VARARGS, rng_gamma_doc},
  {"gamma_int",(PyCFunction)rng_gamma_int,METH_VARARGS, NULL},
  {"flat",(PyCFunction)rng_flat,METH_VARARGS, rng_flat_doc},
  {"lognormal",(PyCFunction)rng_lognormal,METH_VARARGS, rng_lognormal_doc},
  {"chisq",(PyCFunction)rng_chisq,METH_VARARGS, rng_chisq_doc},
  {"fdist",(PyCFunction)rng_fdist,METH_VARARGS, rng_fdist_doc},
  {"tdist",(PyCFunction)rng_tdist,METH_VARARGS, rng_tdist_doc},
  {"beta",(PyCFunction)rng_beta,METH_VARARGS, rng_beta_doc},
  {"logistic",(PyCFunction)rng_logistic,METH_VARARGS, rng_logistic_doc},
  {"pareto",(PyCFunction)rng_pareto,METH_VARARGS, rng_pareto_doc},
  {"dir_2d",(PyCFunction)rng_dir_2d,METH_VARARGS, rng_dir_2d_doc},
  {"dir_2d_trig_method",(PyCFunction)rng_dir_2d_trig_method,METH_VARARGS, rng_dir_2d_trig_method_doc},
  {"dir_3d",(PyCFunction)rng_dir_3d,METH_VARARGS, rng_dir_3d_doc},
  {"dir_nd",(PyCFunction)rng_dir_nd,METH_VARARGS, rng_dir_nd_doc},
  {"weibull",(PyCFunction)rng_weibull,METH_VARARGS, rng_weibull_doc},
  {"gumbel1",(PyCFunction)rng_gumbel1,METH_VARARGS, rng_gumbel1_doc},
  {"gumbel2",(PyCFunction)rng_gumbel2,METH_VARARGS, rng_gumbel2_doc},
  {"poisson",(PyCFunction)rng_poisson,METH_VARARGS, rng_poisson_doc},
  {"bernoulli",(PyCFunction)rng_bernoulli,METH_VARARGS, rng_bernoulli_doc},
  {"binomial",(PyCFunction)rng_binomial,METH_VARARGS, rng_binomial_doc},
  {"negative_binomial",(PyCFunction)rng_negative_binomial,METH_VARARGS, rng_negative_binomial_doc},
  {"pascal",(PyCFunction)rng_pascal,METH_VARARGS, rng_pascal_doc},
  {"geometric",(PyCFunction)rng_geometric,METH_VARARGS, rng_geometric_doc},
  {"hypergeometric",(PyCFunction)rng_hypergeometric,METH_VARARGS, rng_hypergeometric_doc},
  {"logarithmic",(PyCFunction)rng_logarithmic,METH_VARARGS, rng_logarithmic_doc},
  {"landau",(PyCFunction)rng_landau,METH_VARARGS, rng_landau_doc},
  {"erlang",(PyCFunction)rng_erlang,METH_VARARGS, NULL},
  {"multinomial",(PyCFunction)rng_multinomial,METH_VARARGS, multinomial_doc},
  {"dirichlet",(PyCFunction)rng_dirichlet,METH_VARARGS, rng_dirichlet_doc},
  {NULL, NULL,}
};

static PyObject *
rng_getattr(PyGSL_rng *self, char *name)
{
     PyObject *tmp = NULL;

     FUNC_MESS_BEGIN();
     assert(PyGSL_RNG_Check(self));

     tmp = Py_FindMethod(rng_methods, (PyObject *) self, name);
     if(NULL == tmp){
	  PyGSL_add_traceback(module, __FILE__, "rng.__attr__", __LINE__ - 1);
	  return NULL;
     }
     return tmp;
}

static PyObject *
rng_create_list(PyObject *self, PyObject *args)
{
     const gsl_rng_type** thisRNGType=gsl_rng_types_setup();
     PyObject* list = NULL, *item=NULL;

     FUNC_MESS_BEGIN();
     /* provide other rng types as subclasses of rng */

     list = PyList_New(0);
     while ((*thisRNGType)!=NULL) {
	  item = PyString_FromString((*thisRNGType)->name);
	  Py_INCREF(item);
	  assert(item);
	  PyGSL_clear_name(PyString_AsString(item), PyString_Size(item));
	  if(PyList_Append(list, item) != 0)
	       goto fail;
	  thisRNGType++;
     }
     FUNC_MESS_END();
     return list;
 fail:
     Py_XDECREF(list);
     Py_XDECREF(item);
     return NULL;
}

/*---------------------------------------------------------------------------*/
/* Module set up */
#define RNG_GENERATE_PDF
#undef RNG_DISTRIBUTION
#define RNG_DISTRIBUTION(name, function)                                   \
static PyObject* rng_ ## name ## _pdf (PyObject *self, PyObject *args)     \
{                                                                          \
     PyObject * tmp;                                                       \
     FUNC_MESS_BEGIN();                                                    \
     tmp =  PyGSL_pdf_ ## function (self, args, gsl_ran_ ## name ## _pdf); \
     if (tmp == NULL){                                                     \
	  PyGSL_add_traceback(module, __FILE__, #name  "_pdf", __LINE__);  \
     }                                                                     \
     FUNC_MESS_END();                                                      \
     return tmp;                                                           \
}

#include "rng_distributions.h"

static PyObject* rng_dirichlet_lnpdf (PyObject *self, PyObject *args)
{
     PyObject *tmp;
     FUNC_MESS_BEGIN();
     tmp = PyGSL_pdf_dA_to_dA(self, args, gsl_ran_dirichlet_lnpdf);
     FUNC_MESS_END();
     return tmp;
}

static PyObject* rng_multinomial_lnpdf (PyObject *self, PyObject *args)
{
     PyObject *tmp;
     FUNC_MESS_BEGIN();
     tmp = PyGSL_pdf_uidA_to_uiA(self, args, gsl_ran_multinomial_lnpdf);
     FUNC_MESS_END();
     return tmp;
}

static const char rng_env_setup_doc[] =
"This function reads the environment variables `GSL_RNG_TYPE' and\n\
`GSL_RNG_SEED'.\n\
The environment variable `GSL_RNG_TYPE' should be the name of a\n\
generator, such as `taus' or `mt19937'. The environment variable\n\
`GSL_RNG_SEED' should contain the desired seed value. It is\n\
converted to an `unsigned long int' using the C library function\n\
`strtoul'.\n";

static PyObject *
PyGSL_rng_env_setup(PyObject *self, PyObject *args)
{
     gsl_rng_env_setup();
     Py_INCREF(Py_None);
     return Py_None;
}

static PyMethodDef PyGSL_rng_module_functions[] = {
     {"borosh13"        , PyGSL_rng_init_borosh13        , METH_NOARGS, NULL},
     {"cmrg"            , PyGSL_rng_init_cmrg            , METH_NOARGS, NULL},
     {"coveyou"         , PyGSL_rng_init_coveyou         , METH_NOARGS, NULL},
     {"fishman18"       , PyGSL_rng_init_fishman18       , METH_NOARGS, NULL},
     {"fishman20"       , PyGSL_rng_init_fishman20       , METH_NOARGS, NULL},
     {"fishman2x"       , PyGSL_rng_init_fishman2x       , METH_NOARGS, NULL},
     {"gfsr4"           , PyGSL_rng_init_gfsr4           , METH_NOARGS, NULL},
     {"knuthran"        , PyGSL_rng_init_knuthran        , METH_NOARGS, NULL},
     {"knuthran2"       , PyGSL_rng_init_knuthran2       , METH_NOARGS, NULL},
#ifdef  _PYGSL_GSL_HAS_RNG_KNUTHRAN2002
     {"knuthran2002"    , PyGSL_rng_init_knuthran2002    , METH_NOARGS, NULL},
#endif
     {"lecuyer21"       , PyGSL_rng_init_lecuyer21       , METH_NOARGS, NULL},
     {"minstd"          , PyGSL_rng_init_minstd          , METH_NOARGS, NULL},
     {"mrg"             , PyGSL_rng_init_mrg             , METH_NOARGS, NULL},
     {"mt19937"         , PyGSL_rng_init_mt19937         , METH_NOARGS, NULL},
     {"mt19937_1999"    , PyGSL_rng_init_mt19937_1999    , METH_NOARGS, NULL},
     {"mt19937_1998"    , PyGSL_rng_init_mt19937_1998    , METH_NOARGS, NULL},
     {"r250"	        , PyGSL_rng_init_r250            , METH_NOARGS, NULL},
     {"ran0"	        , PyGSL_rng_init_ran0            , METH_NOARGS, NULL},
     {"ran1"	        , PyGSL_rng_init_ran1            , METH_NOARGS, NULL},
     {"ran2"	        , PyGSL_rng_init_ran2            , METH_NOARGS, NULL},
     {"ran3"	        , PyGSL_rng_init_ran3            , METH_NOARGS, NULL},
     {"rand"	        , PyGSL_rng_init_rand            , METH_NOARGS, NULL},
     {"rand48"          , PyGSL_rng_init_rand48          , METH_NOARGS, NULL},
     {"random128_bsd"   , PyGSL_rng_init_random128_bsd   , METH_NOARGS, NULL},
     {"random128_glibc2", PyGSL_rng_init_random128_glibc2, METH_NOARGS, NULL},
     {"random128_libc5" , PyGSL_rng_init_random128_libc5 , METH_NOARGS, NULL},
     {"random256_bsd"   , PyGSL_rng_init_random256_bsd   , METH_NOARGS, NULL},
     {"random256_glibc2", PyGSL_rng_init_random256_glibc2, METH_NOARGS, NULL},
     {"random256_libc5" , PyGSL_rng_init_random256_libc5 , METH_NOARGS, NULL},
     {"random32_bsd"    , PyGSL_rng_init_random32_bsd    , METH_NOARGS, NULL},
     {"random32_glibc2" , PyGSL_rng_init_random32_glibc2 , METH_NOARGS, NULL},
     {"random32_libc5"  , PyGSL_rng_init_random32_libc5  , METH_NOARGS, NULL},
     {"random64_bsd"    , PyGSL_rng_init_random64_bsd    , METH_NOARGS, NULL},
     {"random64_glibc2" , PyGSL_rng_init_random64_glibc2 , METH_NOARGS, NULL},
     {"random64_libc5"  , PyGSL_rng_init_random64_libc5  , METH_NOARGS, NULL},
     {"random8_bsd"     , PyGSL_rng_init_random8_bsd     , METH_NOARGS, NULL},
     {"random8_glibc2"  , PyGSL_rng_init_random8_glibc2  , METH_NOARGS, NULL},
     {"random8_libc5"   , PyGSL_rng_init_random8_libc5   , METH_NOARGS, NULL},
     {"random_bsd"      , PyGSL_rng_init_random_bsd      , METH_NOARGS, NULL},
     {"random_glibc2"   , PyGSL_rng_init_random_glibc2   , METH_NOARGS, NULL},
     {"random_libc5"    , PyGSL_rng_init_random_libc5    , METH_NOARGS, NULL},
     {"randu"           , PyGSL_rng_init_randu           , METH_NOARGS, NULL},
     {"ranf"            , PyGSL_rng_init_ranf            , METH_NOARGS, NULL},
     {"ranlux"          , PyGSL_rng_init_ranlux          , METH_NOARGS, NULL},
     {"ranlux389"       , PyGSL_rng_init_ranlux389       , METH_NOARGS, NULL},
     {"ranlxd1"         , PyGSL_rng_init_ranlxd1         , METH_NOARGS, NULL},
     {"ranlxd2"         , PyGSL_rng_init_ranlxd2         , METH_NOARGS, NULL},
     {"ranlxs0"         , PyGSL_rng_init_ranlxs0         , METH_NOARGS, NULL},
     {"ranlxs1"         , PyGSL_rng_init_ranlxs1         , METH_NOARGS, NULL},
     {"ranlxs2"         , PyGSL_rng_init_ranlxs2         , METH_NOARGS, NULL},
     {"ranmar"          , PyGSL_rng_init_ranmar          , METH_NOARGS, NULL},
     {"slatec"          , PyGSL_rng_init_slatec          , METH_NOARGS, NULL},
     {"taus"            , PyGSL_rng_init_taus            , METH_NOARGS, NULL},
     {"taus2"           , PyGSL_rng_init_taus2           , METH_NOARGS, NULL},
     {"taus113"         , PyGSL_rng_init_taus113         , METH_NOARGS, NULL},
     {"transputer"      , PyGSL_rng_init_transputer      , METH_NOARGS, NULL},
     {"tt800"           , PyGSL_rng_init_tt800           , METH_NOARGS, NULL},
     {"uni"             , PyGSL_rng_init_uni             , METH_NOARGS, NULL},
     {"uni32"           , PyGSL_rng_init_uni32           , METH_NOARGS, NULL},
     {"vax"             , PyGSL_rng_init_vax             , METH_NOARGS, NULL},
     {"waterman14"      , PyGSL_rng_init_waterman14      , METH_NOARGS, NULL},
     {"zuf"             , PyGSL_rng_init_zuf             , METH_NOARGS, NULL},
     {"rng", PyGSL_rng_init_default, METH_NOARGS, rng_doc},
     {"list_available_rngs", rng_create_list, METH_VARARGS},
     /*densities*/
     {"gaussian_pdf",rng_gaussian_pdf,METH_VARARGS, rng_gaussian_pdf_doc},
     {"ugaussian_pdf",rng_ugaussian_pdf,METH_VARARGS, rng_ugaussian_pdf_doc},
     {"gaussian_tail_pdf",rng_gaussian_tail_pdf,METH_VARARGS, rng_gaussian_tail_pdf_doc},
     {"ugaussian_tail_pdf",rng_ugaussian_tail_pdf,METH_VARARGS, rng_ugaussian_tail_pdf_doc},
     {"bivariate_gaussian_pdf",rng_bivariate_gaussian_pdf,METH_VARARGS, rng_bivariate_gaussian_pdf_doc},
     {"exponential_pdf",rng_exponential_pdf,METH_VARARGS, rng_exponential_pdf_doc},
     {"laplace_pdf",rng_laplace_pdf,METH_VARARGS, rng_laplace_pdf_doc},
     {"exppow_pdf",rng_exppow_pdf,METH_VARARGS, rng_exppow_pdf_doc},
     {"cauchy_pdf",rng_cauchy_pdf,METH_VARARGS, rng_cauchy_pdf_doc},
     {"rayleigh_pdf",rng_rayleigh_pdf,METH_VARARGS, rng_rayleigh_pdf_doc},
     {"rayleigh_tail_pdf",rng_rayleigh_tail_pdf,METH_VARARGS, rng_rayleigh_tail_pdf_doc},
     {"gamma_pdf",rng_gamma_pdf,METH_VARARGS, rng_gamma_pdf_doc},
     {"flat_pdf",rng_flat_pdf,METH_VARARGS, rng_flat_pdf_doc},
     {"lognormal_pdf",rng_lognormal_pdf,METH_VARARGS, rng_lognormal_pdf_doc},
     {"chisq_pdf",rng_chisq_pdf,METH_VARARGS, rng_chisq_pdf_doc},
     {"fdist_pdf",rng_fdist_pdf,METH_VARARGS, rng_fdist_pdf_doc},
     {"tdist_pdf",rng_tdist_pdf,METH_VARARGS, rng_tdist_pdf_doc},
     {"beta_pdf",rng_beta_pdf,METH_VARARGS, rng_beta_pdf_doc},
     {"logistic_pdf",rng_logistic_pdf,METH_VARARGS, rng_logistic_pdf_doc},
     {"pareto_pdf",rng_pareto_pdf,METH_VARARGS, rng_pareto_pdf_doc},
     {"weibull_pdf",rng_weibull_pdf,METH_VARARGS, rng_weibull_pdf_doc},
     {"gumbel1_pdf",rng_gumbel1_pdf,METH_VARARGS, rng_gumbel1_pdf_doc},
     {"gumbel2_pdf",rng_gumbel2_pdf,METH_VARARGS, rng_gumbel2_pdf_doc},
     {"poisson_pdf",rng_poisson_pdf,METH_VARARGS, rng_poisson_pdf_doc},
     {"bernoulli_pdf",rng_bernoulli_pdf,METH_VARARGS, rng_bernoulli_pdf_doc},
     {"binomial_pdf",rng_binomial_pdf,METH_VARARGS,  rng_binomial_pdf_doc},
     {"negative_binomial_pdf",rng_negative_binomial_pdf,METH_VARARGS, rng_negative_binomial_pdf_doc},
     {"pascal_pdf",rng_pascal_pdf,METH_VARARGS, rng_pascal_pdf_doc},
     {"geometric_pdf",rng_geometric_pdf,METH_VARARGS, rng_geometric_pdf_doc},
     {"hypergeometric_pdf",rng_hypergeometric_pdf,METH_VARARGS, rng_hypergeometric_pdf_doc},
     {"logarithmic_pdf",rng_logarithmic_pdf,METH_VARARGS, rng_logarithmic_pdf_doc},
     {"landau_pdf",rng_landau_pdf,METH_VARARGS, rng_landau_pdf_doc},
     {"erlang_pdf",rng_erlang_pdf,METH_VARARGS, NULL},
     {"multinomial_pdf",rng_multinomial_pdf,METH_VARARGS,multinomial_pdf_doc},
     {"dirichlet_pdf",rng_dirichlet_pdf,METH_VARARGS, rng_dirichlet_pdf_doc},
     {"multinomial_lnpdf",rng_multinomial_lnpdf,METH_VARARGS, NULL},
     {"dirichlet_lnpdf",rng_dirichlet_lnpdf,METH_VARARGS, rng_dirichlet_lnpdf_doc},
     {"env_setup",PyGSL_rng_env_setup,METH_NOARGS, (char*) rng_env_setup_doc},
     {NULL, NULL, 0}        /* Sentinel */
};

static void
set_api_pointer(void)
{

     FUNC_MESS_BEGIN();
     PyGSL_API[PyGSL_RNG_ObjectType_NUM] = (void *) &PyGSL_rng_pytype;
     DEBUG_MESS(2, "__PyGSL_RNG_API   @ %p,  ", (void *) PyGSL_API);
     DEBUG_MESS(2, "PyGSL_rng_pytype  @ %p,  ", (void *) &PyGSL_rng_pytype);
     /* fprintf(stderr, "__PyGSL_RNG_API @ %p\n", (void *) __PyGSL_RNG_API); */
     FUNC_MESS_END();
}

void
initrng(void)
{
     PyObject *m=NULL, *item=NULL, *dict=NULL;
     PyObject *api=NULL;

     m = Py_InitModule("rng", PyGSL_rng_module_functions);
     assert(m);
     /* import_array(); */
     init_pygsl();

     /* create_rng_types(m); */
     module = m;

     dict = PyModule_GetDict(m);
     if(!dict)
	  goto fail;

     if (!(item = PyString_FromString(rng_module_doc))){
	  PyErr_SetString(PyExc_ImportError,
			  "I could not generate module doc string!");
	  goto fail;
     }
     if (PyDict_SetItemString(dict, "__doc__", item) != 0){
	  PyErr_SetString(PyExc_ImportError,
			  "I could not init doc string!");
	  goto fail;
     }

     PyGSL_rng_pytype.ob_type = &PyType_Type;

     set_api_pointer();
     api = PyCObject_FromVoidPtr((void *) PyGSL_API, NULL);
     assert(api);
     if (PyDict_SetItemString(dict, "_PYGSL_RNG_API", api) != 0){
	  PyErr_SetString(PyExc_ImportError,
			  "I could not add  _PYGSL_RNG_API!");
	  goto fail;
     }

     return;

 fail:
     if(!PyErr_Occurred()){
	  PyErr_SetString(PyExc_ImportError, "I could not init rng module!");
     }

}

max_stars_repo_name: juhnowski/FishingRod

max_issues_repo_name: juhnowski/FishingRod

max_forks_repo_name: juhnowski/FishingRod


text:
#include <viaio/VImage.h>
#include <string.h>
#include <gsl/gsl_matrix.h>
#include <viaio/Vlib.h>
#include "gsl_utils.h"

void whitecov(VImage rho_vol, gsl_matrix_float *Y,
	      gsl_matrix_float* invM, gsl_matrix_float* pinvX,
	      gsl_matrix_float* X, int numlags, int slice)
{

  /* set up variables */
  int n = Y->size1;
  /* these values represent the matlab interval 2:n */
  int k1 = 1;
  int k2 = n-1;
  /* loop counter */
  int i,j;

  /* matrix buffers */
  gsl_matrix_float *result,*buffer;

  /* : X^T */
  gsl_matrix_float* transX = gsl_matrix_float_alloc(X->size2,X->size1);
  gsl_matrix_float_transpose_memcpy(transX,X);

  /* : betahat_ls = pinvX * Y; */
  gsl_matrix_float* betahat_ls = fmat_x_mat(pinvX, Y, NULL);

  /* gsl_blas_sgemm(CblasNoTrans,CblasNoTrans,1.0,pinvX,Y,0.0,betahat_ls); */

  /* : resid = Y - X' * betahat_ls; */
  result = fmat_x_mat(transX, betahat_ls, NULL);
  gsl_matrix_float_sub(Y,result);
  gsl_matrix_float_free(result);

  /* rename */
  gsl_matrix_float* resid = Y;

  if (numlags == 1) {

    /* : Cov0=sum(resid.*resid,1); */
    result = gsl_matrix_float_alloc(resid->size1, resid->size2);
    gsl_matrix_float_memcpy (result, resid);
    gsl_matrix_float_mul_elements (result, result);
    gsl_vector_float* Cov0 = fsum(result, 1,NULL);
    gsl_matrix_float_free(result);

    /* : Cov1=sum(resid(k1,:).*resid(k1-1,:),1); */
    gsl_matrix_float_view sub1 =
      gsl_matrix_float_submatrix(resid, k1,0,k2,resid->size2);
    gsl_matrix_float_view sub2 =
      gsl_matrix_float_submatrix(resid, k1-1, 0,k2,resid->size2);
    result = gsl_matrix_float_alloc(k2, resid->size2);
    gsl_matrix_float_memcpy (result, &sub1.matrix);
    gsl_matrix_float_mul_elements (result, &sub2.matrix);
    gsl_vector_float* Cov1 = fsum(result,1,NULL);
    gsl_matrix_float_free(result);

    /* : Covadj=invM*[Cov0; Cov1]; */
    buffer = gsl_matrix_float_alloc(2,Cov0->size);
    float* pBuffer = buffer->data;
    float* pCov0 = Cov0->data;
    float* pCov1 = Cov1->data;
    for(i=0;i<Cov0->size;i++) {
      *pBuffer++ = *pCov0++;
    }
    for(i=0;i<Cov1->size;i++) {
      *pBuffer++ = *pCov1++;
    }
    gsl_matrix_float* Covadj = fmat_x_mat(invM, buffer, NULL);
    gsl_matrix_float_free(buffer);

    /* :  rho_vol(:,slice,1)=(Covadj(2,:)./ ...
     * (Covadj(1,:)+(Covadj(1,:)<=0)).*(Covadj(1,:)>0))'; */
    int cols = Covadj->size2;
    float* p = Covadj->data;
    float val;
    for(i=0;i<cols;i++) {
      val = *(p+cols) / ((*p <= 0) ? *p +1 : *p) * (*p > 0);
      VPixel(rho_vol, i, slice, 0, VFloat) = val;
      p++;
    }

    gsl_matrix_float_free(Covadj);
    gsl_vector_float_free(Cov0);
    gsl_vector_float_free(Cov1);
  }

  else {
    /* :
     * for lag=0:numlags
     *    Cov(lag+1,:)=sum(resid(1:(n-lag),:).*resid((lag+1):n,:));
     * end
     */
    int lag;
    gsl_matrix_float* Cov = gsl_matrix_float_alloc(numlags+1,resid->size2);
    gsl_vector_float* v;
    for(lag=0;lag<=numlags;lag++) {
      gsl_matrix_float_view sub1 =
	gsl_matrix_float_submatrix(resid,0,0,n-lag,resid->size2);
      gsl_matrix_float_view sub2 =
	gsl_matrix_float_submatrix(resid,lag,0,n-lag,resid->size2 );
      result = gsl_matrix_float_alloc(n-lag, resid->size2);
      gsl_matrix_float_memcpy (result, &sub1.matrix);
      gsl_matrix_float_mul_elements (result, &sub2.matrix);
      v = fsum(result, 1,NULL);
      gsl_matrix_float_free(result);

      memcpy((float*)(Cov->data+(lag*Cov->size2)),
	     (float*)v->data,
	     sizeof(float)*v->size);

      gsl_vector_float_free(v);
    }

    /* :  Covadj=invM*Cov; */
    gsl_matrix_float* Covadj = fmat_x_mat(invM, Cov, NULL);

    /* : rho_vol(:,slice,:)= ( Covadj(0:(numlags+1),:) ...
     *  .*( ones(numlags,1)*((Covadj(1,:)>0)./ ...
     * (Covadj(1,:)+(Covadj(1,:)<=0)))) )'; */
    gsl_matrix_float* line = gsl_matrix_float_alloc(1,Covadj->size2);
    float* pC = Covadj->data;
    float* pB = line->data;
    for(i=0;i<line->size2;i++) {
      *pB = (float)(*pC > 0) / (float)(*pC+(*pC<=0));
      pB++;
      pC++;
    }

    buffer = gsl_matrix_float_alloc(numlags,1);
    gsl_matrix_float_set_all(buffer,1);
    result = fmat_x_mat(buffer, line,NULL);

    gsl_matrix_float_view sub =
      gsl_matrix_float_submatrix(Covadj,1,0,numlags,Covadj->size2);
    gsl_matrix_float_free(buffer);
    buffer = gsl_matrix_float_alloc(result->size1, result->size2);
    gsl_matrix_float_memcpy (buffer, &sub.matrix);
    gsl_matrix_float_mul_elements (buffer, result);

    float* p = buffer->data;
    for(i=0;i<buffer->size1;i++) {
      for(j=0;j<buffer->size2;j++) {
	VPixel(rho_vol,j,slice,i, VFloat) = (float)*p++;
      }
    }

    gsl_matrix_float_free(buffer);
    gsl_matrix_float_free(result);
    gsl_matrix_float_free(line);
    gsl_matrix_float_free(Covadj);
    gsl_matrix_float_free(Cov);
  }

  gsl_matrix_float_free(betahat_ls);
  gsl_matrix_float_free(transX);
}

ext: c

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
/**
 */

#ifndef _DISP_CONF_H
#define _DISP_CONF_H
#include <gsl/gsl_vector.h>
#include "aXe_grism.h"

/**
 * A structure to contain the coefficient of a 2D
 * polynomial that can be used to compute the given
 * coefficient of the polynomial dispersion relation
 * as at a particular i.j location in the image
 */
typedef struct
{
  gsl_vector *pol; /* A vector containing the 2D polynomial coefficients */
  int for_grism;   /* If true, pol is of form sum(a_i x^i, i=0..order)
                      otherwise, it is of form sum(a_i 1/x^i, i=0..order) */
  d_point cpoint;  /* The detector location at which this structure
                      was computed */
  int ID;          /* The ID of the beam for which this coefficient
                      is defined */
  char file[MAXCHAR];
}
dispstruct;

typedef struct
{
  char file[MAXCHAR];

  int ID;
  int for_grism;

  int n_order;
  gsl_vector **all_coeffs;
}
global_disp;

extern float
get_beam_mmag_extract (char *filename, int beamID);

extern float
get_beam_mmag_mark (char *filename, int beamID);

extern int
get_beam_disp_norder (char *filename, int beamID);

extern gsl_vector *
get_beam_disp_order (char *filename, const int for_grism,
                                 int beamID, int order);

extern float
get_disp_coeff_at_pos (char *filename, const int for_grism, int beamID,
                             int order, d_point p);

extern gsl_vector *
get_disp_coeffs_at_pos (char *filename, const int for_grism,
                                    int beamID, d_point p);
extern void
free_dispstruct(dispstruct *p);

extern dispstruct *
get_dispersion_at_pos (char *filename, int beamID, d_point p);

extern void
dispstruct_fprintf (FILE * file, dispstruct * disp);

extern dispstruct *
get_dispstruct_at_pos (char *filename, const int for_grism,
                                   int beamID, d_point p);

extern int
check_for_grism (char *filename, int beamID);

extern int
check_disp_coeff (char *filename,const int for_grism,int beamID,int order);

extern gsl_vector *
get_prange (char *filename, int beamID);

extern global_disp *
get_global_disp(char *filename, const int for_grism, int beamID);

extern gsl_vector *
get_calvector_from_gdisp(const global_disp *gdisp, const d_point p);

extern double
get_2Ddep_value(const gsl_vector *coeffs, const d_point p);

extern void
print_global_disp(const global_disp *gdisp);

extern int
check_disp_order(const gsl_vector *tmp, const int beamID, const int order);

extern void
free_global_disp(global_disp *gdisp);

#endif

ext: h

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe


text:
/****************************************************************************
 * Constants
 ***************************************************************************/

#ifndef CONST_H
#define CONST_H

#include <string>
#include <random>
#include <chrono>
#include <iostream>
#include <gsl/gsl_vector_float.h>
#include <gsl/gsl_blas.h>

using namespace std;

namespace MD
{
  //GLOBAL ALIASES
  typedef float myCoor;
  typedef gsl_vector_float myVect;
  typedef gsl_matrix_float myMatr;

  static const auto& myVectMalloc = gsl_vector_float_calloc;
  static const auto& myVectMemdel = gsl_vector_float_free;
  static const auto& myVectSet = gsl_vector_float_set;
  static const auto& myVectGet = gsl_vector_float_get;
  static const auto& myVectMin = gsl_vector_float_min;
  static const auto& myVectMax = gsl_vector_float_max;
  static const auto& myVectCpy = gsl_vector_float_memcpy;
  static const auto& myVectZro = gsl_vector_float_set_zero;

  static const auto& myMatrMalloc = gsl_matrix_float_calloc;
  static const auto& myMatrMemdel = gsl_matrix_float_free;
  static const auto& myMatrSet = gsl_matrix_float_set;
  static const auto& myMatrSetrow = gsl_matrix_float_set_row;
  static const auto& myMatrGet = gsl_matrix_float_get;
  static const auto& myMatrZro = gsl_matrix_float_set_zero;

  static const auto& myVectAdd = gsl_vector_float_add;
  static const auto& myVectSub = gsl_vector_float_sub;
  static const auto& myVectScl = gsl_vector_float_scale;
  static const auto& myVectMul = gsl_vector_float_mul;
  static const auto& myVectDiv = gsl_vector_float_div;
  static const auto& myMatrAdd = gsl_matrix_float_add;
  static const auto& myMatrScl = gsl_matrix_float_scale;
  static const auto& myVectDot = gsl_blas_sdot;
  static const auto& myVectNrm = gsl_blas_snrm2;
  static const auto& myMatrVecmul = gsl_blas_sgemv;
  //GLOBAL VARIABLES
  static ostream& COUTPT = cout;
  static istream& CINPUT = cin;
  extern const myCoor ONEPIE, TWOPIE, ROOTWO, IVRTWO, COULMB, BOLTZM, MINZRO;
  extern const myCoor STPSZE, STPHLF, STPINV, STPIHF, STEPSQ, STPIRT;
  extern const int DIM;
  extern const int MAXNPART, MAXNBOND, MAXNANGL, MAXNDHED;
  extern const int BONDFORC;
  extern const bool ANGLPART, LJONPART, DPDPARTK, LONGMEMO;
  extern const char *paraFl, *topoFl, *strcFl, *coorFl, *trajFl, *memoFl,
                    *convFl, *defNme, *datNme, *cgStFl;
  extern const char *crdFmt, *typFmt, *nmeFmt, *sNmFmt;
  extern const string ctrlAtom, ctrlMass, ctrlBond, ctrlAngl, ctrlDhed,
                      ctrlRsid, ctrlTime, ctrlStop, ctrlIgnr, ctrlEnd,
                      emptyNme;
  //GLOBAL FUNCTIONS
  static default_random_engine ENGINE;
  static normal_distribution<myCoor> NRMDIS{0.0, 1.0};
  static uniform_real_distribution<myCoor> UNIDIS(0.0,1.0);
  void myVectInv(myVect* vec1, const myCoor leng);
  void myVectRev(myVect* vec1);
  void myVectCrs(myVect* vec1, const myVect* vec2);
  void myVectVecCrs(const myVect* vec1, const myVect* vec2, myVect* vec3);
  void myVectUnt(myVect* uVec, myCoor& dist);
  void myVectRnd(myVect* uVec);
  void myVectOrd(myVect* sVec, const myVect* uVec);
  void myVectInp(myVect* vect, istream& iput);
  void rndVec(myVect* uVec);
  void myVectOut(const myVect* vect);
  void dbuggr();
  //GLOBAL DUMMY VECTOR
  static myVect *FOR1 = myVectMalloc(DIM), *FOR2 = myVectMalloc(DIM),
                *FOR3 = myVectMalloc(DIM), *FTOT = myVectMalloc(DIM),
                *FORC = myVectMalloc(DIM),
                *VELO = myVectMalloc(DIM), *POST = myVectMalloc(DIM),
                *PBC1 = myVectMalloc(DIM), *PBC2 = myVectMalloc(DIM);
  //NOTE: be careful when using these vectors, as they may be
  //      used elsewhere
  static myVect *NRM1 = myVectMalloc(DIM), *NRM2 = myVectMalloc(DIM),
                *LEN1 = myVectMalloc(DIM), *LEN2 = myVectMalloc(DIM),
                *LEN3 = myVectMalloc(DIM), *UNIT = myVectMalloc(DIM);
}
#endif